Layer-dependent electronic and magnetic properties of Nb 3 I 8

.


I. INTRODUCTION
Graphene synthesis in 2004 inaugurated the twodimensional (2D) materials era, which provided a wide range of physical properties and technological applications in high-performance electronic devices.Silicene, phosphorene, hexagonal boron/gallium nitrides, transition-metal dichalcogenides (TMDs), bismuthene, and MXenes are the most famous examples of the still growing family of synthesized 2D materials, which led to great advances in nanoelectronics and optoelectronics technologies [1][2][3][4].Further impressive results have been shown in the study and synthesis of the van der Waals (vdW) homo-and heterostructures.These systems exhibit novel properties and functionalities unavailable in the single-layer constituents [5,6], tunable due to an unprecedented number of degrees of freedom, such as the single layers' relative orientation (twist angles) [7][8][9][10], their order and number (stacking configuration) [11,12], and the distance between two consecutive single layers (interlayer distance) [13,14].
Although 2D materials and vdW homo-and heterostructures have provided remarkable achievements in valleytronics [15][16][17], straintronics [18], and twistronics [7,19], most of the available pristine 2D materials are intrinsically nonmag-netic, and this limits their practical applications in spintronics.Strain [20][21][22], doping [20,23,24], and defect engineering [25][26][27] can induce magnetism in these materials, but these techniques are still not easily controlled in experiments or, in the last two cases, strongly deteriorate the transport properties because of their tendency to form scattering centers.Many recipes have been conceived to induce magnetism in 2D layered materials, since natively magnetic 2D systems should not be observable due to the Mermin-Wagner theorem.The latter states that thermal fluctuations destroy long-range magnetic ordering at finite temperatures in two dimensions within the isotropic Heisenberg model with continuous spin symmetries [28].However, it has been argued that the magnetic anisotropy can induce the breaking of the Hamiltonian spin symmetry, thus suppressing the effect of thermal fluctuations [29][30][31].This explains the first observed native 2D magnetism in two transition-metal layered systems, CrI 3 [32] and Cr 2 Ge 2 Te 6 [33], both already known for being ferromagnetic (FM) below the Curie temperature T C of ∼ 61 K in the bulk form [34,35]. Experimental measurements show a layer-dependent magnetism for CrI 3 [32], being a ferromagnet with T C ∼ 45 K in monolayer, an antiferromagnet (AFM) in bilayer, and again a ferromagnet in trilayer form.First-principles calculations based upon a quantum anisotropic Heisenberg model confirm the T C measurement for CrI 3 monolayer, predicting a ferromagnet with an easy-axis along the out-of-plane direction and a Curie temperature of ∼42.2 K [36].Furthermore, it has been theoretically demonstrated that the low Curie temperature of CrI 3 monolayer can be remarkably enhanced up to ∼60 K in the CrI 3 /MoTe 2 vdW heterostructure and can be further increased to ∼85 K by applying an out-of-plane pressure of ∼4.2 GPa [37].
Ferromagnetic 2D materials are also explored for biomedical applications, such as magnetic labeling and hyperthermia treatment of tumors, but not all of them are biocompatible.Indeed, since Mn 2+ , Mn 3+ ions [38], and Cr 3+ ions [39] are toxic, CrX 3 (X = Cl, Br, I) monolayers [40] and Mncontaining monolayers (MnX 2 with X = S, Se, Mn-CN, MnO 2 , etc.) [41][42][43] with rather high Curie temperatures (for example, Mn 3 C 12 N 12 H 12 and MnSe 2 show T C of 450 K and 250 K, respectively) could cause permanent neurological damage and therefore are not suitable for biomedical applications.
Recently Nb 3 X 8 (X = Cl, Br, or I) [44] monolayers have been suggested as a new class of intrinsic magnetic 2D materials.They are also biocompatible and stable in atmosphere semiconductors due to the Nb ion content.Firstprinciples calculations predict an FM ground state (GS) with, respectively, 31, 56, and 87 K Curie temperatures, estimated with the mean-field theory (MFT) approximation based on the next-nearest-neighbor Ising model.Moreover, niobium iodide, Nb 3 I 8 , bulk crystal was first synthesized by Magonov et al. in 1993 [45] and very recently cleaved to monolayer and multilayered flakes [46], exhibiting remarkable features for the realization of future high-performance nanodevices.
In this work we use the density functional theory (DFT) to explore the electronic and magnetic properties of Nb 3 I 8 in bulk, monolayer, bilayer, and trilayer forms, which will be referred to as Nb 3 I 8 -bulk, Nb 3 I 8 -1L, Nb 3 I 8 -2L, and Nb 3 I 8 -3L, respectively.In particular, for Nb 3 I 8 -2L and -3L we adopt the bulk stacking configuration, denoted here as original stacking.We demonstrate that the inclusion of magnetism induces very interesting properties, making Nb 3 I 8 a suitable candidate for spintronics applications.In particular, this vdW system shows a layer-dependent magnetism, being FM in monolayer form and AFM in bilayer and trilayer forms.Moreover, we show how the magnetism is an essential factor to reach a satisfactory theoretical description of the measured Nb 3 I 8 -1L semiconducting behavior and of the Nb 3 I 8 few layers' work function.Our analysis of the different magnetic phases also allows us to extract, from the ab initio total energy calculations, the magnetic exchange couplings of a third-nearest-neighbor Ising model Hamiltonian describing Nb 3 I 8 -1L.From this Hamiltonian, by using a statistical Monte Carlo (MC) simulation based on the Metropolis algorithm, we derive the Nb 3 I 8 -1L Curie temperature of about 300 K. Finally, we observe how the choice of the vdW exchange-correlation functional to model the vdW interactions dispersion could significantly affect the predictions about the materials properties.
The paper is organized as follows.In Sec.II we present the computational methods and technical details of our calculations.In Sec.III we show and discuss the electronic and magnetic properties for bulk, monolayer, and multilayer systems, together with the calculation of the Curie temperature for the monolayer.Finally, in Sec.IV, we summarize our results and draw some conclusions.

II. METHODS
All calculations were performed using first-principles calculations based on DFT as implemented in the Quantum-ESPRESSO package [47] (version 6.3), based on plane waves and pseudopotentials.The generalized gradient approximation (GGA) is used with projector-augmented wave (PAW) pseudopotentials [48] based on the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [49] to represent the atomic cores [50].The plane wave basis set is truncated using a cutoff energy of 50 Ry for the plane waves and 300 Ry to represent the charge density in all calculations.An adequate vacuum space of 20 Å was set between periodic replicas along the direction orthogonal to the planes (assumed to be the z direction), in order to avoid spurious interactions induced by the periodic boundary conditions.The Brillouin zone (BZ) was sampled using an 8 × 8 × 1 Monkhorst-Pack k-point grid [51] for the normal and FM states, and a 4 × 4 × 1 grid for the 2 × 2 supercell used for the AFM states.These k-point grids have been used for both structural relaxation and selfconsistent runs.As the calculation of the magnetic properties might show poor or slow convergence with respect to the parameters entering the calculation, accurate convergence tests have been carried out.We tested the convergence of the calculated properties (total energy per atom, band structure, Curie temperature) with respect to the wave function cutoff, charge density cutoff and k-point grid.As detailed in the Appendix, the adopted parameters are sufficiently safe to guarantee a very good convergence of the results.
The vdW interaction has been self-consistently accounted for using the vdW-DF2-C09 [52,53] and rev-vdW-DF2 [54] exchange-correlation functionals, which have been compared to show to what extent the treatment of the long-range interaction might influence the predictions of the properties of the hexagonal layered Nb 3 I 8 solid.Both vdW functionals have been proven to be successful in the description of 2D vdW heterostructures with an hexagonal lattice.For example, vdW-DF2-C09 reproduces lattice parameters, interlayer distances, and binding energy in agreement with the experimental data for graphite and pristine bilayer graphene [55].The most recent rev-vdW-DF2 functional provides reasonably small errors for all properties (intralayer and interlayer lattice constants and interlayer binding energy) of different hexagonal layered solids (graphite, h-BN, MoS 2 , MoSe 2 , MoTe 2 , HfTe 2 , WS 2 , WSe 2 ) [56].To describe the on-site Coulomb repulsion of Nb 4d electrons, the LDA + U method has been considered [57][58][59].Based on previous reports for a Nb 3 I 8 monolayer [44], we set U at 2 eV.For all the considered systems, the lattice parameters have been optimized using non-spin-polarized calculations.Using the optimized lattice parameters, the atomic positions are fully optimized (independently for non-spin-polarized and spin-polarized calculations) using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [60][61][62][63], with a convergence criterion such that the total energy difference between consecutive structural optimization steps is less than 10 −4 Ry and that all components of all the forces acting on the atoms are less than 10 −3 Ry/bohr.Small changes (of the order of ∼0.5%) are observed in the optimized atomic positions and interlayer distances, if the optimization is carried out before or after switching on the spin polarization.
The calculation of work function and ionization potential has been carried out starting from the planar average of the electrostatic potential.Let V system ( r) be the total electrostatic potential of a given system and A the area of its unit cell, assumed to be in the xy plane.The planar average of V system ( r) on the surface unit cell is The constant value of the V system (z) in the vacuum sufficiently far from the system is defined as the vacuum level energy E vac .From E vac , the work function (W) of a metallic system is easily calculated as W = E vac − E F , where E F is its Fermi energy.Similarly, the ionization potential (IP) of a semiconducting system can be obtained as IP = E vac − VBM, where VBM is the valence band maximum.
The spin-orbit coupling (SOC) has not been included in the calculations, and we assume that the easy axis of rotation for spin orientation coincides with the z axis, according to Ref. [44].In particular, the Curie temperature for the monolayer system has been calculated starting from the thirdnearest-neighbor Ising model, whose Hamiltonian is where i, j n denotes the pair of nth nearest neighbors i and j sites included only once in the sum, J n is the related magnetic exchange coupling, and s i is the z component of the site i spin.Since the Quantum-ESPRESSO package provides the magnetic moments μ i per site, s i can be obtained by using the relation μ i = g s μ B s i , where g s is the electron spin g-factor that can be approximately set to 2. MC simulations based on the Metropolis algorithm on 2D hexagonal lattices with a size of 200 × 200 unit cells, implementing the periodic boundary conditions, have been performed.Actually, since Nb atoms carry almost all of the magnetic moment, each Nb 3 I 8 single layer can be modeled as a Nb atoms kagome lattice.The exchange couplings have been extracted from the DFT calculations by considering the differences between the FM state and three AFM states, as will be better described in Sec.III C. Then the MC estimation of the Curie temperature, T C,MC , has been compared to the MFT approximation estimation, T C,MFT , which satisfies the equation obtained in the limit of vanishing thermal fluctuations, where S is the total spin, k B is the Boltzmann constant, and q n is the nth (magnetic) nearest-neighbor coordination number at any given site [36].As we shall see, q 1 = q 2 = 2 and q 3 = 4 for Nb 3 I 8 .

A. Structural properties
Nb 3 I 8 is a layered transition metal halide belonging to the family of Nb 3 X 8 (X = Cl, Br, I) crystals.The number of Nb 3 X 8 layers per unit cell depends on the X atom, since the Nb 3 Cl 8 crystal structure has two layers and is termed α-Nb 3 Cl 8 (space group D 3 3d − P3m1, No. 164) [64], whereas Nb 3 X 8 with X = Br, I crystal structure has six layers and is termed β-Nb 3 X 8 (space group D 5 3d − R3m, No. 166) [64].In a single Nb 3 X 8 (X = Cl, Br, I) layer, the Nb atoms form triangular Nb 3 clusters, and each Nb atom is under a distorted octahedral environment of X atoms, as shown in Figs.1(a) and 1(b) for Nb 3 I 8 , respectively.In particular, the unit cell contains three Nb atoms and eight X atoms.Furthermore, the single layer side view in Fig. 1(c) shows that a sheet of Nb atoms is sandwiched between two sheets of X atoms.The two sheets of X atoms, each having four atoms, are not equivalent, resulting in an essentially flat (top) and a corrugated (bottom) I sheet.Indeed, in a sheet three X atoms have almost the same z coordinate, and the remaining one is between Nb sheet and three previous X atoms, whereas in the other sheet a X atom is outside of the other three X atoms sheet.Nb 3 I 8 bulk crystal was first synthesized by Magonov et al. in 1993 [45] and very recently cleaved to monolayer and multilayered flakes [46].In the first case single-crystal x-ray diffraction (XRD), atomic force microscopy, and scanning tunneling microscopy (STM) measurements show a lattice parameter a of 7.600 Å, 7.73 ± 0.31 Å, and 7.76 ± 0.43 Å, respectively.The stacking configuration assumed by the layered structure in the Nb 3 I 8 bulk, i.e., the original stacking, is shown in Fig. 1(e).Each layer is obtained by particular transformations of the first one.For example, the second layer results from the inversion operation of the first one with respect to a triangular Nb cluster center (inversion center), whereas the third layer is shifted with respect to the first one of a certain quantity.Instead, first-principles calculations for the monolayer form based on the Grimme's DFT-D dispersion correction [44] predict a lattice parameter a = 7.675 Å.In this work the group theory analysis reveals that the space group D 5 3d of β-Nb 3 I 8 reduces to the space group C 3v in monolayer and multilayer systems.The calculated lattice parameters a for vdW-DF2-C09 and rev-vdW-DF2 exchange-correlation functionals are reported in Table I, together with the lattice parameters c of the bulk crystal.The lattice parameter a of the bulk form is closer to the XRD result, even if atomic force microscopy and STM provide values only 2% greater than that obtained from XRD.A comparison with the literature first-principles calculations can be performed only for the monolayer: the previous DFT-D result gives a lattice parameter a ∼ 1% larger than that of the present work.
Besides the lattice parameters, other structural parameters can be used to better assess the Nb 3 I 8 structure in the different systems under investigation: the distance between the halogens planes in the ith layer (d i ), the distance between two Nb atoms in a triangular cluster in the ith layer (d Nb−Nb,i ), and the distance between Nb atoms of consecutive layers (interlayer distance d i,i+1 ).They can be easily identified in Fig. 1.While the bulk structural parameters do not show any significant dependence on the vdW functional, some differences show for the monolayer and multilayer systems, as can be inferred from a comparison of Tables II(a) and II(b).In particular, the distances between halogens planes in the single layers, the odd-even interlayer distances d 2i+1,2i+2 , and even-odd interlayer distances d 2i,2i+1 calculated with the rev-vdW-DF2 functional are up to ∼1% larger than those calculated with vdW-DF2-C09.On the other hand, Nb-Nb distances are the same for all systems and independent of the vdW functional.Thus, it can be concluded that the rev-vdW-DF2 functional describes more "repulsive" interlayer and intralayer interactions, with the latter conclusion supported by the larger lattice parameter.

B. Electronic properties of the nonmagnetic ground state
Aiming at bringing out the importance of magnetism in Nb 3 I 8 layered systems, we first discuss the electronic properties of the nonmagnetic (NM) ground state.This will better allow us, for example, to understand the importance of the magnetic interaction to catch the experimental findings regarding the thickness-dependent work function.
The non-spin-polarized band structures of Nb 3 I 8 -bulk, -1L, -2L, and -3L, calculated along the -K-M-line of the BZ [see Fig. 1(d)], are reported in Fig. 2. The comparison between the bands calculated by vdW-DF2-C09 and rev-vdW-DF2 functionals does not show any significant difference, although the geometry has been optimized for each system using the two different functionals (the structural differences, that amount up to 1%, as previously discussed and reported in Table II, are likely too small to induce important modifications on the electronic structure).
NM Nb 3 I 8 -bulk is a semiconductor with a band gap E g = 0.451 eV (0.442 eV) predicted by the vdW-DF2-C09 (rev-vdW-DF2) functional, and a semiconducting behavior also characterizes NM Nb 3 I 8 -2L with E g = 0.477 eV (0.440 eV).In both cases the band gap is indirect, since the VBM is located at the K point of the BZ, whereas the conduction band minimum (CBM) is located at the point.Moreover, TABLE II.Main structural parameters of Nb 3 I 8 -bulk, -1L, -2L, and -3L, as calculated using the vdW-DF2-C09 (a) and rev-vdW-DF2 (b) functionals.The distance between the halogens planes in layers 1 and 3 (d 1 and d 3 , respectively), the Nb-Nb distance within a triangular cluster (d Nb−Nb ), and the odd-even interlayer distances d 2i+1,2i+2 and even-odd interlayer distances d 2i,2i+1 are reported.The results are obtained using non-spin-polarized calculations.The inclusion of the magnetism does not substantially alter these distances, producing variations within ∼1%.Computed band structures of Nb 3 I 8 -bulk, -1L, -2L, and -3L, calculated using both the vdW-DF2-C09 (solid line) and rev-vdW-DF2 (dashed line) vdW functionals in the nonmagnetic ground state.Zero energy is highlighted with a dark green line and corresponds to E 0 , which is E F (VBM) for metallic (semiconducting) systems.Black arrows locate the band gap E g for the Nb 3 I 8 semiconducting systems.
for NM Nb 3 I 8 -bulk, the uppermost three valence bands at the K point are degenerate, and the same is true for the first three conduction bands.On the other hand, due to the odd valence electrons number (95) of a Nb 3 I 8 monolayer, the NM structures with an odd number of layers (Nb 3 I 8 -1L and Nb 3 I 8 -3L) are metallic.
Relevant information can also be inferred from the projected density of states (PDOS) onto atomic orbitals, as shown in Fig. 3. First, the main DOS contributions in Nb 3 I 8 -1L derive from d z 2 (∼33%), d x 2 −y 2 and d xy (∼12%) Nb orbitals [see Fig. 3(a)].A very interesting result is that in Nb 3 I 8 -2L the valence and conduction bands result from a perfect hybridization of the atomic orbitals of the two layers.Instead, Fig. 3(c) shows that in Nb 3 I 8 -3L the addition of a third layer with respect to Nb 3 I 8 -2L dramatically changes this picture, in that at E F the topmost layer provides most of the contribution to the electronic bands.As such, we can conclude that the third layer is responsible for the semiconductor-to-metal transition when moving from Nb 3 I 8 -2L to Nb 3 I 8 -3L.
The description of the electronic properties of the NM ground state can be completed with the discussion of work function (W) or ionization potential (IP), for the metallic or semiconducting systems, respectively.Since, as discussed in Sec.III A, the two sheets of I atoms in a single Nb 3 I 8 layer are not equivalent, these quantities are calculated for both the top and the bottom sides of the selected systems.In Fig. 4(a) the results are reported, together with the experimental data for Nb 3 I 8 monolayer and multilayered flakes [46].Theoretical data show odd-even oscillations due to the semiconducting (metallic) behavior of the systems composed by an even (odd) number of layers.Top and bottom work functions differ by ∼60-70 meV regardless of the vdW functional, whereas the top and bottom ionization potentials are almost identical.Overall, there is a good agreement between the calculated and measured work function.The increasing trend in the experimental results can be identified in the calculated data even though the theoretical estimation seems to underestimate the measured data by ∼0.2 eV.Such a discrepancy might be associated with intrinsic errors induced by the involved approximations (e.g.exchange-correlation functional), as well as the fact that the experimentally determined thickness, for example, is not directly comparable with its theoretical estimation (distance along the z direction between the top and bottom halogen planes).

C. Magnetic properties
The theoretical prediction of a metallic Nb 3 I 8 -1L of the previous calculations disagrees with a recent experimental work [65], where Nb 3 I 8 -1L has been synthesized and has proved to be a semiconductor.Furthermore, Ref. [44] predicted a semiconducting FM ground state.These reasons stimulate the search for a stable magnetic ground state, which can better match the experimental findings.The magnetic properties of Nb 3 I 8 have been studied by considering the FM and three different AFM configurations (referred to as AFM 1 , AFM 2 , and AFM 3 in the following) for Nb 3 I 8 -1L (see also Ref. [44]), Nb 3 I 8 -2L, and Nb 3 I 8 -3L.Only the first has been considered for Nb 3 I 8 -bulk.The total energy of each magnetic ordering state has been accurately evaluated, with the aim of identifying the most stable state.The local spin density approximation (LSDA) has been adopted to compute the spin-polarized band structures, with the magnetization along the z axis (out-of-plane magnetization), as described in Sec.II.
The magnetic configurations of the FM and AFM i , i = 1, 2, 3 states are schematically drawn in Fig. 5 for Nb 3 I 8 -1L.In the case of Nb 3 I 8 -2L and Nb 3 I 8 -3L the second-and thirdlayer magnetic configuration has been obtained from that of the first layer (as in the figure) by considering the symmetry constraints relating the atoms of different layers in the bulk systems (as also briefly discussed in in Sec.III A).
As far as the structural properties are concerned, the magnetic configurations do not significantly alter the interatomic and interlayer distances reported in Table II, with the largest deviations being of the order of ∼1%.As such we do not report again the structural analysis of Table II, which can be considered as a reference for the distribution of the different distances as well as of the performance of the two vdW functionals in catching the geometry of such systems.
It should be pointed out that most of the magnetic moment of the magnetic configurations is carried out by the Nb atoms, as also revealed by the Lowdin analysis for the per atom magnetizations, whereas the adjacent I atoms do not show any significant spin polarization (as also reported in Ref. [44]).In Table III the computed magnetic properties are reported.For all the studied systems a magnetic ground state shows up.However, while the FM magnetic state is the most energetically favored for Nb 3 I 8 -bulk and Nb 3 I 8 -1L, Nb 3 I 8 -2L and Nb 3 I 8 -3L show an AFM ground state, i.e., AFM 3 and AFM 1 , respectively.This result is of utmost importance: Nb 3 I 8 shows a layer-dependent magnetism, being FM in monolayer form and switching on the antiferromagnetism by adding one or two layers.The layer-dependent magnetism demonstrates how the number of layers is a very decisive degree of freedom also for vdW homostructures, in that different compounds, such as Nb 3 I 8 and CrI 3 [32], also may show a significant dependence of their properties on the stacking of the single layers.
The energy differences between the different magnetic orderings do not depend on the vdW functional for Nb 3 I 8 -1L, whereas differences up to few tens of eV show up for the other systems.This is likely due to the different description of the interlayer interaction provided by the functionals.
It should be pointed out that, regardless of the vdW functional used, in the case of Nb 3 I 8 -1L, for which the FM state has the lowest energy, the AFM i energies are the same within less than 0.5 meV, despite the three chosen AFM states being very different from each other.Nb 3 I 8 -2L shows an AFM 3 ground-state total energy differing at most for ∼1 meV by the AFM 2 configuration total energy, whereas the other energy differences have very much larger values.Finally, the total energy of the Nb 3 I 8 -3L AFM 1 ground state differs at most ∼30 meV from the other magnetic configurations, whereas the normal state is the less favored one, also with respect to Nb 3 I 8 -1L and -2L systems.
The choice of three different AFM states allows us to calculate the magnetic exchange couplings J 1 , J 2 , and J 3 in Eq. (2).To this aim, we identify the nearest-, next-nearest-, and third-nearest-neighbor couplings, as shown in Fig. 6.From the latter, we can easily calculate the coordination TABLE III.Magnetic properties of Nb 3 I 8 -bulk, Nb 3 I 8 -1L, Nb 3 I 8 -2L, and Nb 3 I 8 -3L.For each system we set the energy 0 of the most stable magnetic configuration to zero.As such, we report for each magnetic state S = N, FM, AFM 1 , AFM 2 , and AFM 3 (N stands for normal, nonmagnetic state) its energy S referred to 0 .M is the total magnetization in units of Bohr magneton per unit cell.J 1 , J 2 , and J 3 are the extrapolated magnetic exchange couplings of the third-nearest-neighbor Ising model described in Eq. ( 2). , and third-nearest-neighbor couplings (cyan arrows).The coordination numbers are q 1 = q 2 = 2, and q 3 = 4, respectively (see text).
numbers q 1 = q 2 = 2 and q 3 = 4 of each Nb atom up to the third nearest neighbors.
We focus only on the magnetic exchange couplings of the monolayer system, for which both the considered functionals predict the same total magnetization M, as reported in Table III.Since the Lowdin analysis also shows that every Nb atom carries the same magnetic moment per unit cell, the spin z component does not depend on the site i and can be simply denoted by s.Moreover, the s value can be obtained as ∼1/3 of the total magnetization M divided by the g s 2 electron spin g-factor (the factor 1/3 arising from the Nb atoms arrangement in triangular clusters).Using Eq. ( 2), the magnetic energies of FM and AFM states in the 2 × 2 supercell are where E 0 is a reference energy.The energy differences with respect to the FM state energy provide the linear equations system where The resulting magnetic exchange couplings J 1 , J 2 , and J 3 for the monolayer system, as calculated by solving Eqs.(5), are reported in Table III.The sign of J 1 and its magnitude with respect to J 2 and J 3 confirm the FM state as the ground state.The calculation of the magnetic exchange couplings is the starting point for the MC simulations and MFT approximation estimations of the Curie temperature of Nb 3 I 8 -1L.Their values could be further improved by considering a Heisenberg model Hamiltonian containing a magnetic anisotropy.This would require also taking onto account in-plane magnetizations and the spinorbit coupling though.In Fig. 7 the spin-polarized band structures of the ground state for each system are reported for the most stable magnetic configuration.Unlike the nonmagnetic calculations, the spin- polarized band structures depend on the vdW functional.This suggests that the discrepancies originate from the structural differences predicted by the functionals, as will be evident for Nb  IV.
If compared with the nonmagnetic state (see Fig. 2, first panel), the rev-vdW-DF2 FM band structure of Nb 3 I 8 -bulk keeps the semiconducting character, but the magnetic ordering changes its band edges and band gap character (from indirect to direct) depending on the spin channel.The energy gap difference between spin-up and spin-down channel might certainly exploited in optoelectronics applications for the fully spin-polarized excited states induced by suitable visible photons absorption.The Nb 3 I 8 -1L band structure [see third and fourth panels of Fig. 7(a)] is almost independent of the vdW exchangecorrelation functional, as expected for the lack of vdW interlayer interactions.A semiconducting behavior is observed for both the spin-up and spin-down channel, and this is in agreement with the experimental results of Ref. [65].This is the first evidence to our best knowledge that the magnetism on one hand is necessary to correctly describe the Nb 3 I 8 lowest-energy state, and on the other is responsible for the metallic-to-semiconducting "transition" when switched on in the NM state.More importantly, Nb 3 I 8 -1L, as confirmed by the experimental observations and as will be examined in Sec.III D, reveals as a room-temperature ferromagnet.Table IV shows that even in this case the spin-up and spin-down band gaps are very different.
The Nb 3 I 8 -2L and Nb 3 I 8 -3L band structures in the AFM 3 and AFM 1 states, respectively, show semiconducting behaviors, obviously independent of the spin channel [see first and FIG. 7. Computed spin-polarized band structures for Nb 3 I 8 -bulk, -1L (a), -2L, and -3L (b).While for Nb 3 I 8 -1L a FM ground state is predicted, band structures of Nb 3 I 8 -2L and Nb 3 I 8 -3L correspond to AFM 3 and AFM 1 ground states, respectively.Blue (red) bands represent the spin-up (spin-down) channel.Zero energy is highlighted with a dark green line and corresponds to E 0 , which is E F (VBM) for metallic (semiconducting) systems.second panels of Fig. 7(b)].The band gaps are indirect (K-), with the values reported in Table IV.
It should be pointed out that by changing the vdW functional, different gaps and bands dispersions are obtained.It can be argued that the most important reason why this happens is that the two different functionals induce modifications in the structural properties of the systems (see Sec. III A), that in turn are reflected into the band structure.Several tests confirmed such an argument, in that if, at a fixed geometry, the band structure is computed with the two functionals, no significant changes of the electronic properties can be observed.
We should expect that the significant modification to the ground-state electronic properties observed when magnetism is switched on might have a deep impact on the work function/ionization potential too.Such quantities, related to the energy required to remove the most loosely bound electron from the material (either at the Fermi energy or at the VBM) and bring it to a point in the vacuum immediately close to the material surface, can be directly compared, as we have already done for the NM ground state, with the experimental measurements of the thickness-dependent work function of Ref. [46].This comparison is shown in Fig. 4(b).Even if magnetic Nb 3 I 8 -4L, -5L, -6L data have not been examined in this paper (because those systems, mostly in the AFM state, become computationally very demanding and hard to converge), we observe that the magnetism evidently reduces the discrepancy between our results and experimental data, producing a much better agreement.Indeed, as far as Nb 3 I 8 -2L and -3L are concerned, the difference with the experimental data reduces from about 100 meV (200 meV) to about 30 meV (60 meV) for the rev-vdW-DF2 (vdW-DF2-C09) exchange-correlation functional.This analysis demonstrates the need of taking into account the magnetic interaction to describe the ground state of few-layer Nb 3 I 8 .As well, it also gives evidence that the rev-vdW-DF2 exchange-correlation functional is the most appropriate to reproduce the experimental data.

D. Calculation of the Curie temperature of Nb 3 I 8 -1L
The magnetic exchange couplings obtained in Sec.III C can be used to estimate the Nb 3 I 8 -1L Curie temperature, which is a key parameter for FM spintronics applications.The first estimation results from the MC simulation described in detail in Sec.II, which provides magnetic moment m(T ) and magnetic susceptibility χ (T ) as functions of the temperature.They show the sought magnetic phase transition at T C,MC = 307.5K for both vdW-DF2-C09 and rev-vdW-DF2 functionals.As an example, m(T ) and χ (T ), as computed using the rev-vdW-DF2 functional, are plotted in Fig. 8.This result, which predicts a FM state at a room temperature, certainly improves previous literature results from the MFT [44], since we have accurately taken into account both the lattice type and the number of neighbors of each Nb atom.Of course, although a MC simulation is a reliable method to calculate the magnetic phase transition related to Eq. ( 2), this estimation could be improved by considering a Heisenberg model Hamiltonian containing a magnetic anisotropy, as specified in Sec.III C. It must be pointed out that J 2 and J 3 couplings have a remarkable effect on the Curie temperature, in that including the second-and third-nearest-neighbor interaction produces an approximately 50 K decrease of T C,MC with respect to the firstneighbor approximation, regardless of the vdW functional.Accurate convergence tests, as detailed in Appendix A 2, show a full convergence of the calculated Curie temperature with respect to the parameters entering the calculation (k-point grid size, charge density cutoff).
In Table V we can directly compare Nb 3 I 8 -1L to the most investigated and recent ferromagnetic monolayers predicted by first-principles calculations, in terms of Curie temperature and band gap.Each T C value has been obtained by Monte Carlo simulations except for one case.Nb 3 I 8 -1L is a new competitive 2D magnetic material, due to its Curie temperature and intermediate band gap value.Thus, it already occupies a special position in the 2D magnets scenario, even if a more accurate theoretical treatment, including SOC, might provide even more interesting results for spintronics applications and might be devised.
The Curie temperature has been also estimated within the MFT approximation, as calculated from Eq. ( 3).The values of T C,MFT estimations are ∼908 K (∼948 K) for vdW-DF2-C09 (rev-vdW-DF2).As expected, the MFT approximation overestimates the Curie temperature.The results are very different from the previous literature value 87 K [44], and also if we neglect J 3 to consider the same accuracy of the latter.In any event, the scantiness of MFT approximation with respect to the MC simulation appears evident.

IV. CONCLUSIONS
In summary, we investigated the structural, electronic, and magnetic properties of Nb 3 I 8 in bulk, monolayer, and some multilayer (bilayer and trilayer) forms.Magnetism has turned out to be necessary to reproduce the experimental data on Nb 3 I 8 -1L semiconducting behavior and Nb 3 I 8 few layers work function.Nb 3 I 8 -1L shows a room Curie temperature of 307.5 K, as predicted by Monte Carlo simulations for both vdW-DF2-C09 and rev-vdW-DF2 exchangecorrelation functionals (the mean-field theory approximation estimation 900 K unambiguously overestimates its value).We observed how the number of layers significantly affects the magnetic ordering of the lowest energy state, since it is possible to switch on the antiferromagnetism by adding one or two layers to the FM monolayer.Our results for the electronic and magnetic properties of Nb 3 I 8 render it an ideal and promising candidate for future spintronics research and applications.
Our first-principles calculations also provide a further benchmark of vdW-DF2-C09 and rev-vdW-DF2 vdW exchange-correlation functionals for hexagonal van der Waals nanostructures, with the result that rev-vdW-DF2 seems to produce more reliable results as far as the comparison with the available experimental data is concerned.For all the calculations the plane wave basis set was truncated using a cutoff energy of 50 Ry for the plane waves and 300 Ry to represent the charge density in all calculations.The Brillouin zone (BZ) was sampled using an 8 × 8 × 1 Monkhorst-Pack k-point grid for the normal and FM states, and a 4 × 4 × 1 grid for the 2 × 2 supercell used for the AFM states.These k-point grids have been used for both structural relaxation and self-consistent runs.As detailed in the following, thorough tests have been carried out to show that our choice can be considered safe to ensure the convergence of all the calculated properties, from the structural and electronic ones to the relative stability of the different magnetic phases to the Curie temperature.The ferromagnetic (FM) monolayer Nb 3 I 8 (Nb 3 I 8 -1L) using the rev-vdW-DF2 exchange-correlation functional was used as test case.
The adopted cutoff energy of 50 Ry for the plane waves basis set gives a convergence of the total energy within 15 meV per atom (that reduces to ∼5 meV at 60 Ry and 1 meV at 70 Ry); however, much better convergence (within very few meVs) is obtained for the energy differences (e.g., those used to evaluate the relative stability of the different magnetic phases).We considered 50 Ry as the best compromise between computational resources needed to carry out all the needed calculations and the accuracy of the calculated   properties.At a fixed wave-function cutoff value of 50 Ry, we considered nks × nks × 1 k-point grids with nks = 4, 6, 8, 10, 12 and increased the ecutrho (charge density cutoff) starting from 225 Ry to 475 Ry in steps of 25 Ry.The results are shown in Fig. 9. Interestingly, we notice no dependence on the k-point grid, showing that the 8 × 8 × 1 we used produces fully converged results.On the other hand, as far as ecutrho in concerned, while more significant variations are observed for the lowest considered values, starting from 300 Ry the changes become as small as 1 meV.A more insightful analysis with respect to ecutrho is depicted in Fig. 10, where we fix the 4 × 4 × 1 k-point grid and increase ecutrho up to 1000 Ry.While we might consider the total energy with ecutrho = 500 Ry as fully converged, that with ecutrho = 300 Ry differs from the latter by only 0.75 meV per atom and was chosen for all the other calculations.
Since the Curie temperature, as predicted by the Monte Carlo simulations (T C,MC ), can significantly depend of the energy differences between the different magnetic phases, it is worth carrying out the convergence analysis also on the energy differences E AFM i − E FM .The results are reported in Fig. 11, for nks × nks × 1 k-points grids with nks = 8, 12, 16 and ecutrho = 300 Ry, 375 Ry, and 400 Ry.The energy differences are reported in meV.We observe that all the calculated values are within an interval of a less than 0.3 meV amplitude, with largest percentage variations with respect to the ecutrho = 300 Ry case being 0.16% for E AFM 1 − E FM , 0.46% for E AFM 2 − E FM , 0.25% for E AFM 3 − E FM and 8 × 8 × 1 k-point grid, and 0.06% for E AFM 1,2 − E FM and 0.21% for E AFM 3 − E FM with a 12 × 12 × 1 k-point grid.It is worth pointing out that, while only for the 8 × 8 × 1 k-point grid the relative stability of the three AFM phases slightly changes (since they are degenerate within 0.3 meV), the FM phase is  (K ), for all the combinations of parameters considered in Fig. 11.
definitely the most stable for all the values of the considered parameters (ecutwfc, ecutrho, k-point grid).We are going to show that such small variations cannot affect the Curie temperature estimation.
Finally, to conclude this convergence analysis, we performed a spin-polarized band structure calculation for FM trilayer Nb 3 I 8 (Nb 3 I 8 -3L) using the rev-vdW-DF2 exchangecorrelation functional, taken as an example, and a 16 × 16 × 1 k-point grid with the ecutrho values of 400 Ry and 600 Ry, to be compared with the results reported in the paper, obtained with a 8 × 8 × 1 k-point grid and ecutrho = 300 Ry.The results are shown in Fig. 12.In agreement with the conclusions drawn so far, the spin-polarized band structures show only negligible changes at the largest considered size of k-point grid and value of ecutrho.Indeed, a substantially perfect superposition of bands is shown for both the spin-up and spin-down cases.
Based upon these calculations and the subsequent ones, we retain that the 8 × 8 × 1 k-point grid and an ecutrho value of 300 Ry are good choices to describe this system, being a good compromise between a safe convergence of the calculated results and the computational resources needed to carry out spin-polarized calculations also containing the Hubbard U parameter.

Curie temperature
Based upon the calculations shown in Fig. 11, we computed the Curie temperature T C,MC for different values of the parameters, and we report the results in Table VI.
As already pointed out before, the very small absolute (within ∼0.3 meV) and percentage variations on the energy differences E AFM i − E FM do not significantly affect the Curie estimation, even if they lead to differences up to 20 K on T C,MC (∼6%).It is remarkable that the value calculated with the 8 × 8 × 1 k-point grid and ecutrho 300 Ry is just the same as the most converged result, calculated with the 16 × 16 × 1 k-point grid and ecutrho 400 Ry.

FIG. 1 .
FIG. 1. Nb 3 I 8 structure.(a) Top view of the Nb 3 I 8 -1L unit cell.(b) Distorted octahedral environment of I atoms around each Nb atom.(c) Side view of a Nb 3 I 8 -1L sheet.(d) 2D BZ. b 1 and b 2 are the reciprocal lattice primitive vectors.The high-symmetry points , K, M for the band structure calculation along the -K-M-path are shown.(e) Side view of the bulk crystal unit cell (containing six Nb 3 I 8 monolayers).

FIG. 3 .
FIG. 3. (a) PDOS onto the valence atomic orbitals (p and d for Nb, p for I) of NM Nb 3 I 8 -1L and layer-resolved DOS of (b) NM Nb 3 I 8 -2L and (c) NM Nb 3 I 8 -3L.Only the data calculated by the rev-vdW-DF2 exchange-correlation functional are reported, as an example.

FIG. 4 .
FIG. 4. Comparison between experimental measurements of the work function (W) in Nb 3 I 8 multilayered flakes and the theoretical estimations of the work function, W (ionization potential, IP), of NM Nb 3 I 8 multilayers with odd (even) number of layers.The results obtained for both the NM (a) and the lowest energy magnetic ground state (b) are shown.
3 I 8 -3L and Nb 3 I 8 -bulk.First of all, for Nb 3 I 8 -bulk vdW-DF2-C09 predicts a metallic FM ground state, whereas rev-vdW-DF2 predicts a semiconducting ground state with the same magnetic ordering.Spin-up and spin-down channels show the expected six-band bundle close to the Fermi energy or the VBM.However, in the case of the majority spin and for Nb 3 I 8 -1L and Nb 3 I 8 -3L, another band bundle can be clearly identified, centered at ∼0.5 eV [see first and second panels of Fig.7(a)].As far as the semiconducting rev-vdW-DF2 band structure is concerned, the spin-up (spin-down) VBM and CBM both are located at the K ( ) point of the BZ, with a band gap of 0.254 eV (0.983 eV) reported in Table

FIG. 8 .
FIG.8.Magnetic moment m(T ) and magnetic susceptibility χ (T ) (normalized to 1) for Nb 3 I 8 -1L as results of Monte Carlo simulations.Data refer to rev-vdW-DF2 functional, as an example.The vertical line highlights the Curie temperature, whereas the cyan solid line is a fit intended as a guide for the eye.

FIG. 9 .
FIG. 9. Total energy per atom (meV) as a function of ecutrho (Ry) at a fixed wave-function cutoff (50 Ry) for Nb 3 I 8 -1L using the rev-vdW-DF2 exchange-correlation functional.We considered nks × nks × 1 k-point grids with nks = 4, 6, 8, 10, 12 for ecutrho starting from 225 Ry to 475 Ry in steps of 25 Ry.The energies are referred to the total energy E 0 of the run with a 4 × 4 × 1 k-point grid and ecutrho = 225 Ry.The number of atoms in the unit cell is N at = 11.

FIG. 10 .
FIG. 10.Total energy per atom (meV) as a function of ecutrho (Ry) for Nb 3 I 8 -1L using the rev-vdW-DF2 exchange-correlation functional and a 4 × 4 × 1 k-point grid.The energies are referred to the total energy E 0 of the run with a 4 × 4 × 1 k-point grid and ecutrho = 225 Ry.The number of atoms in the unit cell is N at = 11.

FIG. 12 .
FIG.12.Spin-polarized band structure calculation for FM Nb 3 I 8 -3L using the rev-vdW-DF2 exchange-correlation functional, 8 × 8 × 1 k-point grid and ecutrho = 300 Ry, 16 × 16 × 1 k-point grid and ecutrho value of 400 Ry and 600 Ry.Spin-up bands are reported in panel (a), spin-down bands in panel (b).Zero energy is highlighted with a dark green line and corresponds to E 0 , which is the valence band maximum of the spin-up channel.

TABLE I .
Lattice parameters a and c for Nb 3 I 8 bulk, monolayer, and multilayer.Both experimental (XRD, atomic force microscopy, and STM) and first-principles theoretical results are reported.The present work results are obtained using non-spinpolarized calculations.
a This work.

TABLE IV .
Spin-up (spin-down) band gap energies E g,↑ (E g,↓ ) in eV of the studied systems in their most stable magnetic configuration and for the two different vdW exchange-correlation functionals.GS stands for the (magnetic) ground state.Zero gap means a metallic band structure.

TABLE V .
Ferromagnetic semiconducting monolayers predicted by first-principles calculations based on both density functional theory and Monte Carlo simulations (except for one case).For every material, together with the calculation method, we point out whether the SOC is included and report the Curie temperature T C and band gap E g .
a Nonlinear spin wave theory.b This work.

TABLE VI .
Estimation of the Curie temperature predicted by the Monte Carlo simulations, T C,MC