Computational Atomistic Modeling in Carbon Flatland and Other 2D Nanomaterials

: As in many countries, the rise of nanosciences in Belgium has been triggered in the eighties in the one hand, by the development of scanning tunneling and atomic force microscopes offering an unprecedented possibility to visualize and manipulate the atoms, and in the other hand, by the synthesis of nano-objects in particular carbon nanostructures such as fullerene and nanotubes. Concomitantly, the increasing calculating power and the emergence of computing facilities together with the development of DFT-based ab initio softwares have brought to nanosciences ﬁeld powerful simulation tools to analyse and predict properties of nano-objects. Starting with 0D and 1D nanostructures, the ﬂoor is now occupied by the 2D materials with graphene being the bow of this 2D ship. In this review article, some speciﬁc examples of 2D systems has been chosen to illustrate how not only density functional theory (DFT) but also tight-binding (TB) techniques can be daily used to investigate theoretically the electronic, phononic, magnetic, and transport properties of these atomically thin layered materials.


Introduction
The invention of the scanning tunneling microscope (STM) in 1981 [1,2], winning of the Nobel Prize in Physics by the inventors in 1986, and the invention of the atomic force microscope (AFM) in the same year [3], opened the doors to the Small World described by Richard Feynman in his visionary lecture There's Plenty of Room at the Bottom [4]. This new ability to see and eventually manipulate individual atoms was soon complemented by other prominent experimental achievements-the discovery of fullerenes in 1985 [5], the development of nanoscale transistors, the structural assignement of carbon nanotubes in 1991 [6]-that definitely established nanosciences as the major research field that we know today. Around the same time, the first quantum-simulation codes were developed based on different flavors of theory such as Hartree-Fock (HF), Density Functional Theory (DFT), Møller-Plesset perturbation (MP2), Configuration Interaction (CI), Coupled Cluster (CC), or multi-reference methods [7][8][9][10][11], facilitating the theoretical investigation of materials and nanostructures at the atomic scale. A major advantage of these quantum-mechanical modelling methods compared for example to the semi-emirical tight-binding (TB) techniques [12] is their ab initio, i.e. parameter-free, character. Nevertheless, the computational cost of ab initio methods remains much higher than the one of TB techniques. In UCLouvain, all these new exciting tools for nanosciences arrived in the nineties. The first AFM and STM were acquired in 1994 [13,14]. The first DFT and TB simulators were developed and used in the early nineties [15][16][17][18] and the ABINIT project was launched officially in 1997 [19].
In this review article, published in the special issue entitled 'State of the art of nanosciences in Belgium', a small country being sometimes called the flatland [22], the predictive power of DFT and TB techniques for various properties of 2D materials is illustrated. Through a selection of monolayered systems, namely graphene, borophene, MXene, and h-BN, the potential use of 2D materials in many fields of applications including electronics, electron optics, valleytronics, phononics, and spintronics, will be exposed. The article includes four sections, each investigating a property of a given 2D system using a specific modeling technique. The first section describes the electronic transport properties in polycrystalline graphene using the Green's function approach within the Landauer-Büttiker formalism. The second section presents the optical properties of borophene and graphene monolayers including optical Hall effect. The third section introduces the recent and large family of MXenes and in particular their versatile electronic and phononic properties. Last but not least, the genuine use of h-BN in innovative magnetic tunnel junctions is discussed in the fourth section. With the present review paper, we wish you an exciting modeling experience in a quite localized region of the still hugely unexplored Flatlands!

Electronic Transport in Polycrystalline Graphene
Although graphene was first isolated by mechanical exfoliation [20], it is currently synthesized for the most part by chemical vapor deposition (CVD) techniques. This approach is presently by far the most viable and industry-compatible production method [23,24]. Importantly, CVD yields polycrystalline graphene (PG) where single-crystal domains are separated by graphene grain boundaries (GBs) of irregular shapes. Grains have a roughness comparable to the substrate and a lateral size ranging from a few hundred nanometers to several centimeters [25]. At the grains interface, topological defects accommodate for the lattice misorientation between domains. Those GBs consist of non-hexagonal rings (mostly pentagon-heptagon pairs) randomly distributed along the interface [26,27]. It was revealed by scanning probe microscopy that such extended structural defects induce major perturbations in the electronic structure of graphene [28,29]. The transport of charge carriers is degraded on the one hand by structural defects acting as scattering sites [30][31][32] and on the other hand by electrical potential barriers formed at GBs which tend to be negativey-doped (n-doped) [33,34]. A strong correlation was demonstrated between the quality and the resistivity of GBs [30,32,35]. Macroscopically, the carrier mobility of wafer-scale graphene is drastically limited by GBs [36][37][38]. Nonetheless, Banszerus et al. measured mobilities for CVD graphene rivaling those reported for similar devices based on exfoliated graphene [39]. This suggests that CVD samples are comparable in quality to large-area exfoliated single-crystals that contain intragranular line defects and other structural imperfections [40].
In this context, an atomic-scale analysis of the GBs detrimental effects on the transport properties of PG is required, with so far only a few first-principles studies characterizing charge transport through GBs. Importantly, Yazyev et al. [41] showed that GBs degrading role could be alleviated through atomically precise manufacturing. A transport gap was also demonstrated to appear for specific types of GBs. However, experimentally-observed GBs exhibit more complicated arrangements of non-hexagonal rings [42] than atomic structures of ref. [41]. This necessarily leads to discrepancies between experimental measurements in PG and theoretical predictions of transport properties of "ideal" GBs. Our researches focus on the systematic study of the electronic transport through realistic graphene grain boundaries [43,44] and on the exploration of novel phenomena such as valley-dependent transport, electron optics effects, Aharonov-Bohm interferences and spin-diffusion lengths in polycrystalline graphene [45][46][47].
Graphene GBs atomic structures were obtained using classical molecular dynamics with optimized Tersoff potentials [48]. This empirical potential model has been shown to accurately describe the structural properties of defective graphene [49]. The electronic properties of defective graphene systems under study were calculated using either the first-principles approach [44] or the π-orbital TB technique [43,45,46]. First-principles calculations were performed using DFT within the OPENMX code [50][51][52], a software package based on norm-conserving pseudo-potentials and pseudo-atomic localized basis functions. Parametrized TB models are very effective to model the electronic properties of graphene at low energies, i.e., around the Fermi level. The good agreement between these two approaches was demonstrated in refs. [41,45].
Charge transport through GBs was simulated using the Landauer-Büttiker formalism implemented within the Green's function approach. In this framework, the zero-temperature linear conductance is given by the Landauer formula: where W is the transverse width of the system, T is the transmission probability across the system, a function of carrier energy E and transverse momentum k y . The function T (E, k y ) is calculated through a fully quantum-mechanical approach: where the retarded Green's function G is determined in terms of the device Hamiltonian H D as: In Equation (3), Σ L(R) represents the left (right) self-energy, describing the left (right) lead-to-device coupling. Γ L(R) in Equation (2) is the transfer rate at the left (right) contact determined as Γ L(R) = i(Σ L(R) − Σ † L(R) ). To investigate quantum transport in PG, GBs are modeled as periodic lines of defects with a periodicity that is expressed in terms of two translation vectors d L = (n L , m L ) and d R = (n R , m R ) in the left and right domains, respectively. Yazyev et al. [41] demonstrated that GBs exhibit either metallic or semi-conducting transport regime depending on the interface nature. If n L − m L = 3q but n R − m R = 3q (or vice versa; q, q ∈ Z), a transport gap appears as the allowed momentum-energy manifolds of the two crystalline domains are significantly misaligned. In the case of massless Dirac fermions, this transport gap is given by: where v F is the Fermi velocity, and | d| the common GB period (after applying possibly strain to obtain matching of both domains). Interestingly, the transport gap depends only on the GB periodicity.
GBs are subsequently classified as either transparent (class-I) or reflective (class-II) to electronic transport. First, a symmetric GB with d = (2, 1) and | d| = 6.5 Å is investigated as illustrated in Figure 1. The system is transparent, as expected from its symmetry, with a conductance comparable to pristine graphene despite the large concentration of pentagon-heptagon pairs (see Figure 1a,c). In the asymmetric case (see Figure 1b,d), a transport gap of 0.92 eV appears when d L = (5, 0), d R = (3, 3) and | d| = 12.6 Å, in relative accordance with Equation (4). Outside the transport gap, transmissions are considerably smaller compared to pristine graphene. However, systems considered in Figure 1a,b have short periodicities and are not representative of real GBs which exhibit a more complicated atomic structure. Consequently, more realistic GBs constructed from more complicated arrangements of polygons are studied hereafter. First, an asymmetric class-I GB is examined. This system is defined by d L = (11,3), d R = (10, 5) and | d| = 32.0 Å, as illustrated in Figure1e. The corresponding transmissions are significantly lower than pristine graphene, indicating an increased scattering due to the presence of dangling bonds. An asymmetric class-II GB is then considered, with d L = (24, 0), d R = (16,12) and | d| = 59.6 Å (see Figure 1f). The interface was constructed by creating a vacancy after repeating n times an initial pattern (here n = 4), which corresponds to applying a perturbation approach to the periodicity assumption. The related transmissions depicted in Figure 1d show a transport gap of 0.76 eV. However, a logarithmic scale plot analysis reveals non-zero conductance as large as 10 −2 G 0 nm −1 in the transport-forbidden energy range. This implies that the transport gap depends on the smallest perceptible periodicity while perturbations in periodicity induce leakage currents. In a nutshell, disorder is found to normalize the transport properties of PG (metallic or semi-conducting) [44]. On the one hand, the disordered nature of GBs increases scattering, reducing the overall transmissions in transparent systems, on the other hand, it induces leakage currents which damage the transport gap for reflective configurations.
Owing to the misorientation between grains, the transport through GBs has been also shown to be highly sensitive to the effect of strains [43]. This is essentially due to the fact that even though the metallicity of graphene is still conserved, its Brillouin zone is strongly deformed and Dirac cones can be displaced in the momentum space when strains are applied. Moreover, such displacements of the Dirac cones are strongly dependent of the crystallographic orientation of graphene with respect to the strain direction. On this basis, the strain effects in PG can result in different momentum displacements of Dirac cones in two graphene grains surrounding a given GB, thus strongly modifying the low-energy transport properties through the system. In particular, while the transport gap in the boundary of class-II is modifiable, a uniaxial strain applied in proper directions can always open a gap in the boundary of class-I as illustrated in Figure 2. Interestingly, a large gap of a few hundred meV can be achieved with a small strain of a few percents only. The value of this gap has been shown to depend on three main parameters: the uniaxial strain magnitude, the uniaxial strain direction, and the lattice symmetry of graphene grains [43]. In addition to its dependence on the GBs and crystallographic orientation of grains, the transport through PG also represents interesting magnetotransport phenomena. For instance, it has been shown in ref. [53] that the formation of Landau levels is limited by the grain size and cannot be observed until the magnetic length is smaller than the average grain radius. The electron localization is also found to be unusual. In particular, strongly localized states and extended electronic states are at the center and between Landau levels, respectively. These extended states are developed along the network of GBs, explaining a finite value for the bulk dissipative conductivity and suppression of the quantized Hall conductance.
A novel fascinating magnetotransport feature has been recently predicted in polycrystalline graphene [46]. In particular, the extended states strongly localized along the grain boundaries can be exploited to create tunneling paths connecting quantum Hall edge channels in CVD graphene samples. Consequently, strong Aharonov-Bohm (AB) interferences in the quantum Hall regime can be achieved in the PG systems schematized in Figure 3a. The AB oscillations in the conductance in a typical system are indeed illustrated in Figure 3b. These predicted AB oscillations are clarified by the pictures of the local density of left-injected electronic states in Figure 3c,d and can be essentially explained as follows. First, under a strong magnetic field, electrons in the internal graphene grain have to propagate following quantum Hall edge states. Second, because of defect scattering at the GBs, these electrons can circulate inside such graphene grain (i.e., along its edge). Therefore, depending on the accumulated phase of electron wave in such circulating trajectories, constructive/destructive states are created and tuned when varying the magnetic field. These constructive (see Figure 3c) and destructive (see Figure 3d) states explain the obtained high and low values of conductance, respectively, and the conductance oscillations presented in Figure 3b.

Optical Absorption in Graphene and Borophene
Owing to their extreme thinness, 2D materials are expected to absorb only a small amount of light. Therefore, even metallic 2D layers can be almost transparent in particular in the visible range. Such transparent conducting layers are of broad interest for various applications, e.g., as front current collector in solar panels. To assess the optical properties of 2D materials, a key quantity is the optical conductivity σ, whose real part is related to the imaginary part of the dielectric function ε (Equation (5)). Based on vertical interband electronic transitions, σ can be routinely computed with both TB and DFT approaches using the Kubo formula (Equation (6)).
(σ αβ (hω)) = ωε 0 (ε αβ (hω)), where α, β are the principal in-plane directions (x and y), S is the cell surface, f (E) is the Fermi-Dirac distribution function, E n,m (k), |n and |m are the eigenenergies and eigenstates,v α/β is the velocity operator along a given direction, and η is a small broadening. The longitudinal optical conductivity components are σ xx and σ yy while σ xy is the optical Hall conductivity. Both are given in units of the universal value of σ 0 = e 2 4h (σ 0 = απε 0 c with α, ε 0 and c the fine-structure, the vacuum electric permittivity and the speed of light constants, respectively). As a consequence of intraband transitions in metallic systems, Drude corrections are added to ε (The real part of ε can also be corrected by the where ω p (ω p = ne 2 mε 0 with n the electron density and m an effective mass. In practise, ω p is computed through Kramers-Kronig relation) and γ=τ −1 are the plasmon frequency and the inverse scattering time of the Drude model, respectively. Light propagation is considered orthogonal to the 2D layers and the associated linearly oscillating electric field is thus in-plane along the α direction.
Photon absorption is therefore by essence a 3D phenomenon and to account for the 2D nature of the atomically-thin structure within DFT calculation, the computed dielectric function ε has to be renormalized [58,59] (including both interband transitions and Drude terms), as ε 2D = ε 3D H h with H the total height of the DFT 3D periodic cell which includes the vacuum distance introduced to separate the periodic images of the layer, and h the polarizable electronic thickness which can be estimated e.g., from van der Waals atomic radii [59]. To compute the transmittance T, general Fresnel equations should be used or at least their limit for (ultra-)thin film which include secondary reflections at top and back side interfaces with external media [59][60][61], here vacuum: where n = √ ε is the complex refractive index. The synthesis of monolayers of boron atoms, so-called borophene, has been reported in 2015-2016 [62,63] showing a rich structural polytype. In particular three important structures have been reported, including a buckled triangular lattice and two flat layers named β 12 and χ 3 (Figure 4a). These three borophene layers are (semi-)metallic as clearly shown by their electronic bandstructure and density of states (Figure 4b). Nevertheless and interestingly, the buckled triangular borophene layer displays a zero optical conductivity up to 3 eV (Figure 4c) and, hence, absorbs light only in the UV range or for higher photon energy (with a strong anisotropic character between σ xx and σ yy ) [59]. In contrast, β 12 and χ 3 borophene layers absorb light in the visible range starting from ∼0.5 and ∼1.0 eV, respectively, with a value varying between 0.5 σ 0 and 1.5 σ 0 . For comparison, graphene exhibits a value slowly increasing from 1.0 σ 0 up to roughly 2.0 σ 0 at the end of the visible range. This value of 1.0 σ 0 for graphene corresponds to the well-known absorbance of πα ∼ 2.3% of light, which gives, with reflectance being nearly zero (∼0.05%), the high transmittance displayed in Figure 5. Again, buckled triangular borophene is fully transparent in this visible region neglecting the Drude corrections which, depending on the scattering time τ, can give rise, at most, to a decrease of transmittance down to 96% in the visible range, i.e., comparable to graphene.
Linearly polarized light can be rotated when going through a material in the presence of a magnetic field. This phenomenon is known as Faraday rotation and originates from the optical Hall effect (OHE) for which σ xy (ω) = 0. The OHE has been investigated in thin films and more recently in 2D materials. Strain engineering has been proposed as an alternative to magnetic fields to induce finite optical Hall conductivity [64]. Indeed, without any symmetry breaking mechanisms, σ xy is normally zero after integrating over the whole Brillouin zone (BZ). Uniaxial strain in graphene, but also shear strain, allows to split in energy the van Hove singularities and the related optical absorption peaks (see σ xx in Figure 6). Because of this strain-induced asymmetry, the BZ integration at a given energy yields a positive or negative σ xy (see σ xy in Figure 6 for light linearly polarized along oblique direction φ = 45 • and for two strain directions θ = 0 • and 90 • corresponding to armchair and zigzag respectively) which are responsible for clockwise or anti-clockwise Faraday rotation (Figure 6(inset)). The Faraday rotation angle is however limited to a few degrees for a strained monolayer and its exact value depends on the strain amplitude and the orientation with respect to light polarization. Additionally, the effect can be amplified using a few layers although at the cost of more absorption and, hence, a reduction of the transmitted light intensity.   σ xx graphene (φ=45°) (str=8%,θ=0°) σ xy graphene (φ=45°) (str=8%,θ=0°) σ xx graphene (φ=45°) (str=8%,θ=90°) σ xy graphene (φ=45°) (str=8%,θ=90°)

Electronic Structure and Vibrational Properties of MXenes
Among the recently discovered 2D systems, the emerging family of 2D transition metal carbides and nitrides, known as MXenes [65][66][67], stands out because of the wide chemical diversity allowing for materials property tuning. In contrast to graphite-like layered materials which can be mechanically exfoliated to obtain 2D flakes, MXenes are obtained from the chemical treatment of three-dimensional MAX phases [65]. MAX phases are layered ceramics with the general formula M n+1 AX n , where M represents an early transition metal, A an element from groups 13 to 16, X either a carbon or a nitrogen atom, and n varies from 1 to 3 [68,69]. Due to the use of etching agents, MXenes are always terminated with functional surface groups, such as -F, =O, and -OH, that are randomly distributed [67]. Since the discovery of the first MXene, Ti 3 C 2 T z , at Drexel University in 2011 [65], more than 30 MXenes have been synthesized, and the stability and properties of dozens more have been investigated using ab initio calculations [70]. The exploration of new MAX phases and derivative MXenes is of great interest to further control materials properties and highlight potential applications, such as catalysis, energy storage and related electrochemical applications [70]. Up to now, most of the research has been devoted to the first MXene Ti 3 C 2 T z and consists in experimental works [70]. However, given the wide chemical versatility on both M and T sites, ab initio calculations are powerful to investigate the structural and functional properties of the resulting huge variety of MXenes compositions.
The geometries and electronic properties of the V 2 CT z systems are investigated using DFT calculations. The ground-state structure of the pristine V 2 C shows a hexagonal lattice with the P3m1 space group (Figure 7a). Based on the unit cell of V 2 C, functionalized V 2 CT z MXene structures with terminations are constructed. Two types of hollow sites on the surface of V 2 C can be distinguished (Figure 7a), leading to different functionalization models [71]. The energetically favorable model for functionalized V 2 CT z corresponds to MD2 where the two functional groups are located on fcc sites for which no C atom is present under the V atom (Figure 7b). Similarly to graphene and other 2D materials with hexagonal symmetry, an armchair direction and a zigzag direction can be identified, as depicted in Figure 7c. Moreover, moving to the reciprocal space, the Γ-M and Γ-K high-symmetry paths in the Brillouin zone respectively correspond to the armchair and zigzag direction (Figure 7c. The electronic band structure of pristine V 2 C and fully-terminated V 2 CT 2 with T = -F and/or -OH are presented in Figure 8. The metallic character of V 2 C is ensured by the presence of V-3d orbitals at the Fermi level, and is not affected by any kind of terminations [72]. This is not always verified since Khazaei et al. showed that the metallic pristine Ti 2 C could turn into a semiconductor upon functionalization, depending on the nature of the terminal groups [71]. In general, all pristine MXenes are metallic, while terminated MXenes can be either metallic or semiconducting, depending on the configuration and nature of the terminal groups. This is a powerful property that allows to tune MXenes electronic properties playing on their composition.  The vibrational properties of the V 2 CT z systems are investigated using density functional perturbation theory (DFPT). The three atoms of the primitive cell of the V 2 C monosheet give rise to nine phonon modes at the Γ point of the Brillouin zone. Similarly, the five, seven, and six atoms respectively in the V 2 CF 2 , V 2 C(OH) 2 , and V 2 CF(OH) primitive cells lead to 15, 21, and 18 phonon modes at the Γ point of the Brillouin zone. The phonon spectra of V 2 C-based systems along the armchair (ΓM) and zigzag (ΓK) directions are presented in Figure 9. The stability of all systems is verified since all phonon frequencies are positive in the whole Brillouin zone. The phonon dispersions have three acoustic modes. Two of them exhibit a linear dispersion near Γ and correspond to in-plane rigid-body motions. In contrast, the third acoustic mode corresponding to out-of-plane vibration has a quadratic dispersion close to Γ and a lower energy in the rest of the spectrum. This quadratic dependence is analogous to the one observed in graphene [73] and MoS 2 [74]. In Figure 9a, the first three optical branches in the phonon spectrum of pristine V 2 C demonstrate significantly lower frequencies close to the three acoustic phonon branches, thus inducing a phonon gap of about 250 cm −1 between the lower three and the upper three optical branches. Upon functionalization, this phonon gap between optical phonons is filled by additional optical phonon branches (Figure 9b-d). This makes the most noticeable difference between pristine and terminated MXene systems. Using the group theory, one can have access to the Raman and infrared activity of the phonon modes, hence allowing for an accurate analysis of both Raman and infrared experimental spectra, presented in ref. [72]. Additionally, the knowledge of the phonon spectrum of a material is essential to understand a wide set of macroscopic properties, including the electrical and thermal conductivities.

Magnetoresistance in 2D Tunnel Junctions
Since the isolation of graphene, the library of synthesizable two-dimensional materials has been enriched with several hundreds items featuring versatile electronic properties. This catalog offers unprecedented design opportunities for the rapidly growing field of spintronic where the structural quality of interfaces often plays an instrumental role in determining the functionality and performance of devices. A prototypical spintronic device is the magnetic tunnel junction (MTJ) consisting of two ferromagnetic (FM) layers separated by a thin insulating (I) barrier.
The instrumental property of MTJs is the modulation of their resistance by individually switching the magnetization of the FM electrodes via the application of an external magnetic field. MTJs are thus characterized by two states of electrical resistance associated with the parallel and anti-parallel magnetizations of the electrodes. This is quantified by the tunnel magnetoresistance (TMR) defined as: where R AP,P and G AP,P correspond to the resistance and conductance of the junction in the anti-parallel and parallel configurations. This effect was discovered by Julliere in 1975 [75] who proposed a simple model to rationalize the measured TMR signal based on the assumption that (i) the carriers spin is conserved upon tunneling and (ii) the conductance of a given spin channel is proportional to the product of the density of conducting states in the FM electrodes, Here, ρ ↑ L/R (ρ ↓ L/R ) denotes the density of conducting states of majority (minority) spin in the left and right electrodes. Within this framework, the TMR can be expressed as: where P L,R are the spin polarization of the conducting electrons in the left and right electrodes defined . The spin polarization of the tunneling electrons ultimately depends on the electronic structure of the interfaces. Over the last decade, significant research efforts have been dedicated to introduce 2D materials to control and tailor the magnetoresistive response of the interfaces [76]. Boosted by the theoretical prediction of perfect spin filtering by few of its layers [77][78][79][80], graphene has been the first 2D material to be investigated [81][82][83], rapidly followed by insulating h-BN [84][85][86] and semi-conducting transition metal dichalcogenides (TMDCs) [87][88][89][90]. While very promising experimental results have yet been reported, first-principles simulations provide a complementary tool to unravel the full potential of 2DMs in tunneling junctions. This is illustrated here with the theoretical investigation of a MTJ consisting of a monolayer h-BN sandwiched between two Co ferromagnet electrodes. The in-plane lattice constant of the 0001 hcp Co matches almost perfectly the lattice constant of h-BN. This suggests that epitaxial Co/h-BN interface can be realized experimentally by CVD technique [91,92]. As a consequence, the interface may be perfectly commensurate or not depending on the integration pathway. The epitaxial interface is illustrated in Figure 10 along with a misaligned configuration obtained by twisting the h-BN monolayer by 21.78 • . First-principles calculations reveal strikingly different spin behaviours for these two geometries.
Structural optimization of the aforementioned interfaces indicates a stronger bonding of h-BN with Co in the aligned configuration with respect to the misaligned one. The average distance between h-BN and the FM is of 2.1 Å for the epitaxial structure and of 3.1 Å for the misaligned one. The charge transfer from the ferromagnet to the monolayer and the hybridization of h-BN are consequently affected by the relative crystallographic orientation across the interface. As a result, the computed spin polarization of h-BN (see Figure 10a,c) exhibits opposite sign at the Fermi energy for the two considered geometries. Straight application of the Julliere model, see Equation (10), predicts radically different magnetoresistive behaviour for the aligned and misaligned configurations. Further calculations of the electronic transmission functions by means of the Green's function formalism, see Equation (2), confirm the dependence of the spin-polarized tunneling current on the relative crystallographic orientation of the interfaces. The electronic transmission functions computed for a fully epitaxial Co 1 /h-BN/Co 2 junction (i.e., with Co 1 , h-BN, and Co 2 aligned) is reported in Figure 11a for the parallel and anti-parallel magnetizations of the electrodes. In the anti-parallel configuration the transmission functions associated with majority and minority spins are degenerated as expected for a symmetric junction. At equilibrium, the computed TMR is of 36%. The case of an asymmetric junction, where one electrode has been twisted by 21.78 • , is considered in Figure 11b. In this configuration, the global amplitude of the transmission is reduced with respect to the epitaxial junction and the spin symmetry of the transmission in the anti-parallel configuration is lifted. The computed TMR at equilibrium is reduced to 15%. This illustrates the strong dependence of the magnetoresistive behaviour of 2D-MTJs on the atomic configuration of interfaces. The calculations were performed within the DFT framework as implemented in the SIESTA electronic structure package [93,94]. Converged k-point grids and a basis set of double-ζ plus polarization quality were used.
(a) (b) Figure 11. The spin dependent electronic transmission function of a (a) symmetric and (b) asymmetric Co 1 /h-BN/Co 2 junction with both the parallel and anti-parallel relative orientation of FM magnetization. The transmission functions of the majority and minority spins are depicted in blue and red respectively. In the symmetric junction, Co 1 , h-BN and Co 2 are exactly aligned hence giving rise to two epitaxial interfaces (see Figure 10b). In the asymmetric case, Co 1 and h-BN are aligned but Co 2 has been twisted by an angle of 21.78 • hence corresponding to a less strongly bonded interface (see Figure 10d).

Conclusions
In this review article covering the state-of-the-art of DFT and TB methods devoted to the prediction of 2D materials properties, a particular emphasis was put on the electrical, optical, vibrational, and magnetic properties. Physical quantities such as the charge and spin transmission (across GBs and MJTs), the optical conductivity and normal transmittance, and the influence of chemical terminations on electronic states and phonon modes, are key parameters used to understand the response of 2D metallic materials (i.e., graphene, borophene, MXenes, and h-BN) and hetero-junctions, to fundamental excitations, such as non-oscillating and oscillating electric fields and thermal gradients.
The tremendous properties of 2D materials have proven their outstanding capability to bring new paradigms in the vast field of nanosciences... and the family of 2D crystals is continuously growing, both in terms of variety and number of materials. Indeed, theoretical search for novel 2D materials has been proposed, focusing on how easily the monolayers can be exfoliated from their parent layered compounds [95]. Starting from 108,423 unique, experimentally known 3D compounds, high-throughput calculations using van der Waals DFT allowed the identification of 1825 compounds that are either easily or potentially exfoliable. Since the physical properties of monolayers can often dramatically change from those of their parent 3D materials, these ab initio predicted 2D systems provide a new degree of freedom for applications while also unveiling novel physics. Thanks to this huge number of potentially new 2D layers, computational atomistic modelling will still have a major role to play in the accurate prediction of the corresponding and possibly promising properties using state-of-the-art quantum simulations.
Author Contributions: All the authors contributed equally to this work and should thus be considered as co-first authors. Consequently, the authors are listed in alphabetic order, except for our esteemed supervisor reported as last author. All authors have read and agreed to the published version of the manuscript.