## 3 Computational Aspects  
First-principles calculations of FeSi\({}_{2}\) structures were performed using the quantum ESPRESSO code [38] imposing periodic boundary conditions. Vanderbilt ultra-soft pseudopotentials were used [39], and the generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE) for the exchange-correlation energy and potential [40] was selected. The criterion of convergence for the total energy was 1.36 \(\times\) 10\({}^{-6}\) eV/formula unit. To speed up the convergence, each eigenvalue was convoluted by a Gaussian with a width \(\sigma\) = 0.272 eV. The cut-off energy was 489.6 eV. Structural optimization was carried out  
Figure 3: Approximant structures for \(\alpha\)-FeSi\({}_{2}\), \(\beta\)-FeSi\({}_{2}\), and \(\Theta\)-C16 (top to bottom order) described on a bcc lattice (left), on an fcc lattice (middle), and with bcc \(z\)-axis projection (right). Black (red) circles refer to TM (Si), whereas open circles correspond to vacancies, \(\mathcal{V}\). For the structural representations on the right, half-open half black circles correspond to alternate occupation of TM (or \(\mathcal{V}\) ) and \(\mathcal{V}\) (or TM) along the [001] direction. (For references to colors in this figure, the reader is referred to the web version of this article).  
using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [41] by simultaneously relaxing both the cell basis vectors and the atomic coordinates. Structural optimization was considered to be complete when the atomic forces were less than 7 meV/A, the stresses smaller than 0.05 GPa, and the total energy during the structural optimization iterative process was changed by less than 1.36 meV. Brillouin zone integration was carried out using the Monkhorst-Pack special sets of \(\mathbf{k}\) points with (4 4 4), (8 8 6), (8 8 8), and (4 4 8) grids for the 48- and 45-atom supercells, tetragonal 3-atom cell, and 12-atom cubic and monoclinic cells, respectively.  
In the first-principles molecular dynamics (FPMD) simulations, we considered the initial 48-atom tetragonal supercells constructed of the tetragonal 6-atom \(\gamma\)- and \(\alpha\)-cells via the (2 \(\times\) 2 \(\times\) 2) translation. One-step FPMD simulations at different temperatures were carried out using the NPT ensemble (constant number of atoms-pressure-temperature). Pressure was applied via the Parrinello and Rahman method [42]. The time step was about 10\({}^{-15}\) s. The system temperature was kept constant by rescaling the velocity with a tolerance of \(\pm\)50 and \(\pm\)100 K for low- and high-temperature simulations. All structures reached their time-averaged equilibrium, and afterwards, the total energy and cell volume of the structures fluctuated only slightly around that equilibrium value. The \(\Gamma\) point was used to save computing time without compromising on accuracy. After FPMD simulations, all the structures were relaxed. The \(\beta\)-FeSi\({}_{2}\) structure was generated after FPMD simulations of the 48-atom \(\gamma\)-supercell at 200 K with subsequent relaxation.  
The quantum ESPRESSO code was employed to study the phonon spectra of the \(\alpha\)- and \(\gamma\)-phases in the framework of the density-functional perturbation theory (DFPT) [43]. The phonon densities of states (PHDOS) of \(\beta\)-FeSi\({}_{2}\) and \(\alpha\)-FeSi\({}_{2}\) were determined using the Phonopy code [44]. The elastic moduli were calculated using the ElaStic code [45]. Given the calculated elastic moduli and cell volumes, the fracture toughness was estimated using the procedure described in Ref. [46].  
Intermediate phases were identified at certain stages of FPMD simulations. The identification of these structures was carried out via group-theoretical analysis using the \"ISOTROPY\" code [47] and following the rule discussed in Ref. [48]. The XcrydSden code [49] was used to visualize the structures. A detailed description of the applied procedures for the identification of different structures was given in our previous papers [48,50].