Electronic structure and magnetic properties of few-layer Cr$_2$Ge$_2$Te$_6$: the key role of nonlocal electron-electron interaction effects

TL;DR

This study investigates the electronic and magnetic properties of few-layer Cr2Ge2Te6, demonstrating the importance of nonlocal electron-electron interactions and hybrid functionals over DFT+U for accurate predictions consistent with experimental data.

Contribution

The paper highlights the failure of DFT+U in predicting properties of Cr2Ge2Te6 and shows that hybrid functionals accurately reproduce experimental observations, revealing new intra-layer exchange couplings.

Findings

DFT+U fails to predict correct ground states, especially in bulk form.

Hybrid functionals match ARPES experimental data.

Discovery of two distinct intra-layer second-neighbor exchange couplings in bilayer Cr2Ge2Te6.

Abstract

Atomically-thin magnetic crystals have been recently isolated experimentally, greatly expanding the family of two-dimensional materials. In this Article we present an extensive comparative analysis of the electronic and magnetic properties of ${\rm Cr}_2{\rm Ge}_2{\rm Te}_6$, based on density functional theory (DFT). We first show that the often-used ${\rm DFT}+U$ approaches fail in predicting the ground-state properties of this material in both its monolayer and bilayer forms, and even more spectacularly in its bulk form. In the latter case, the fundamental gap {\it decreases} by increasing the Hubbard-$U$ parameter, eventually leading to a metallic ground state for physically relevant values of $U$, in stark contrast with experimental data. On the contrary, the use of hybrid functionals, which naturally take into account nonlocal exchange interactions between all orbitals, yields good account of the available ARPES experimental data. We then calculate all the relevant exchange couplings (and the magneto-crystalline anisotropy energy) for monolayer, bilayer, and bulk ${\rm Cr}_2{\rm Ge}_2{\rm Te}_6$ with a hybrid functional, with super-cells containing up to $270$ atoms, commenting on existing calculations with much smaller super-cell sizes. In the case of bilayer ${\rm Cr}_2{\rm Ge}_2{\rm Te}_6$, we show that two distinct intra-layer second-neighbor exchange couplings emerge, a result which, to the best of our knowledge, has not been noticed in the literature.

Figures

Tables

Table 1. Table 1: Irriducible set of atomic coordinates (Wyckoff positions in fractions of the crystallographic cell lattice vectors), main distances (in Å), and angles (in degrees) for bulk Cr 2 ​ Ge 2 ​ Te 6 subscript Cr 2 subscript Ge 2 subscript Te 6 {\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6} . We present both the experimental values (column labelled by “Exp.”) and those obtained by relaxing the atomic coordinates (columns labelled by “spGGA” and “HSE06”, depending on the used approximation)—see Sect.  II.2 . We did not find any appreciable differences between the internal coordinates of monolayer, bilayer, and bulk Cr 2 ​ Ge 2 ​ Te 6 subscript Cr 2 subscript Ge 2 subscript Te 6 {\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6} .

	Exp.Carteaux1995 ; Siberchicot1996			spGGA			HSE06
	$x$	$y$	$z$	$x$	$y$	$z$	$x$	$y$	$z$
Cr-18(f)	0.0	0.0	0.3302	0.0	0.0	0.3343	0.0	0.0	0.3341
Ge-6(c)	0.0	0.0	0.059	0.0	0.0	0.0585	0.0	0.0	0.0577
Te-6(c)	0.663	-0.033	0.2482	0.6692	-0.0386	0.2492	0.6690	-0.0334	0.2487
Cr-Te (Å)	2.749			2.7629			2.7802
${\theta }_{\mathrm{Cr}-\mathrm{Te}-\mathrm{Cr}}$ (deg)	91.6			91.081			90.326
Cr-Cr (Å)	3.942			3.942			3.942
Ge-Ge (Å)	2.459			2.4054			2.3723
Te-Ge (Å)	2.578			2.6022			2.5721
vdW gap (Å)	3.353			3.395			3.374

Table 2. Table 2: Atomic magnetic moments (in units of the Bohr magneton μ B subscript 𝜇 B \mu_{\rm B} ) obtained with the HSE06 functional and the GGA + U GGA 𝑈 {\rm GGA}+U approximation ( U = U sc = 3.9 ​ eV 𝑈 subscript 𝑈 sc 3.9 eV U=U_{\rm sc}=3.9~{}{\rm eV} ) for monolayer, bilayer, and bulk Cr 2 ​ Ge 2 ​ Te 6 subscript Cr 2 subscript Ge 2 subscript Te 6 {\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6} . Bilayer Cr 2 ​ Ge 2 ​ Te 6 subscript Cr 2 subscript Ge 2 subscript Te 6 {\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6} presents two inequivalent Cr atoms, Cr a superscript Cr a {\rm Cr}^{\rm a} and Cr b superscript Cr b {\rm Cr}^{\rm b} —see Fig.  1 g) and Sect.  IV . The magnetic moments lie along the easy-axis.

		HSE06	$\mathrm{GGA}+U_{\mathrm{sc}}$
Monolayer	$\mathrm{Cr}$	3.353	3.645
	$\mathrm{Ge}$	0.028	0.014
	$\mathrm{Te}$	-0.121	-0.167
Bilayer	${\mathrm{Cr}}^a$	3.355	3.662
	${\mathrm{Cr}}^b$	3.361	3.667
	$\mathrm{Ge}$	0.029	0.013
	$\mathrm{Te}$	-0.123	-0.165
Bulk	$\mathrm{Cr}$	3.378	3.698
	$\mathrm{Ge}$	0.027	0.013
	$\mathrm{Te}$	-0.124	-0.169

Table 3. Table 3: (Color online) Exchange couplings J i subscript 𝐽 𝑖 J_{i} and magneto-crystalline anisotropy energy (MAE) for monolayer, bilayer, and bulk Cr 2 ​ Ge 2 ​ Te 6 subscript Cr 2 subscript Ge 2 subscript Te 6 {\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6} . In black: results obtained in this work. Different rows refer to different approximations (see main text). In red: results available in the cited literature. (If a cell is empty it means that the corresponding result was not found in the literature.) All the exchange couplings have been calculated with the FSMA approach (see Sect.  IV ). The only exception is the row labeled by “ spGGA ′ superscript spGGA ′ {\rm spGGA}^{\prime} ”, which presents results calculated with the FM-AFM method (see Sect.  I ).

Monolayer
	$J_1\left(\mathrm{meV}\right)$		$J_2\left(\mathrm{meV}\right)$		$J_3\left(\mathrm{meV}\right)$		MAE ($\mu \mathrm{eV}/\mathrm{Cr}$)
spGGA	-5.577		-0.046		-0.484		62
${\mathrm{spGGA}}^{\prime }$	-10.485		0.344		-0.678		62
HSE06	-6.236		-0.083		-0.268		-80
Xu2018Xu2018	-6.64						250
Li2014Li2014b	-5.18
Sivadas2015Sivadas2015	-3.68		0.39		0.43
Zhuang2015Zhuang2015	-3.07						-220
Fang2018Fang2018	-5.662		0.040		-0.129		-107
Dong2019Dong2019	-6.7						-320
Bilayer
	$J_1\left(\mathrm{meV}\right)$	$J_2^a/J_2^b\left(\mathrm{meV}\right)$	$J_3\left(\mathrm{meV}\right)$	$J_{z1}\left(\mathrm{meV}\right)$	$J_{z2}\left(\mathrm{meV}\right)$	$J_{z3}\left(\mathrm{meV}\right)$	MAE ($\mu$eV/ $\mathrm{Cr}$)
spGGA	-5.694	-0.0429/-0.001	-0.629	-0.096	-0.024	-0.365	40
HSE06	-6.399	-0.096/-0.076	-0.432	-0.017	0.196	-0.199	-140
Fang2018Fang2018	-6.338	0.209	-0.475	-0.177			-274
Bulk
	$J_1\left(\mathrm{meV}\right)$	$J_2\left(\mathrm{meV}\right)$	$J_3\left(\mathrm{meV}\right)$	$J_{z1}\left(\mathrm{meV}\right)$	$J_{z2}\left(\mathrm{meV}\right)$	$J_{z3}\left(\mathrm{meV}\right)$	MAE ($\mu$eV/ $\mathrm{Cr}$)
spGGA	-5.755	-0.028	-0.450	-0.1097	-0.043	-0.359	-80
HSE06	-6.607	-0.083	-0.627	-0.010	0.190	-0.225	-260
Gong2017Gong2017 ; Li2018	-3.76	0.08	-0.16	0.05	-0.12	-0.38	-50
Fang2018 Fang2018	-6.995	0.346	-0.773	-0.480			-471

Equations

$${\hat{\mathcal{H}}}_{\mathrm{spin}}=\sum _{i<j}{J}_{ij}{\hat{\mathbf{S}}}_i\cdot {\hat{\mathbf{S}}}_j+\sum _i{A}_i{\hat{S}}_{iz}^2.$$

$$\begin{array}{rl}E=& E_{\mathrm{spin}}+E_0=\\ =& J_{12}{\mathbf{S}}_1\cdot {\mathbf{S}}_2+{\mathbf{S}}_1\cdot {\mathbf{K}}_1+{\mathbf{S}}_2\cdot {\mathbf{K}}_2+E_{\mathrm{others}}+E_0,\end{array}$$

$$E_1=E_0+J_{12}{S}^2+{\mathbf{S}}_1\cdot {\mathbf{K}}_1+{\mathbf{S}}_2\cdot {\mathbf{K}}_2,$$

$$E_2=E_0-J_{12}{S}^2+{\mathbf{S}}_1\cdot {\mathbf{K}}_1-{\mathbf{S}}_2\cdot {\mathbf{K}}_2,$$

$$E_3=E_0-J_{12}{S}^2-{\mathbf{S}}_1\cdot {\mathbf{K}}_1+{\mathbf{S}}_2\cdot {\mathbf{K}}_2,$$

$$E_4=E_0+J_{12}{S}^2+{\mathbf{S}}_1\cdot {\mathbf{K}}_1-{\mathbf{S}}_2\cdot {\mathbf{K}}_2.$$

$$J_{12}=\frac{E_1-E_2-E_3+E_4}{4S^2}.$$

Full text

I Introduction

The family of two-dimensional (2D) materials beyond graphene is now very large Geim2013 ; Mounet2018 and contains also magnetic crystals Shabbir2018 ; Burch2018 ; Gong2019 , notably ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ Gong2017 ; Wang2018 ; Kim2019 , ${\rm Cr}{\rm I}_{3}$ Huang2017 ; Wang2018a ; Klein2018 ; Huang2018 ; Song2018 ; Jiang2018a ; Kim2018 ; Jin2018 ; Jiang2018b ; Song2019 ; Thiel2019 , ${\rm V}{\rm Se}_{2}$ Bonilla2018 , ${\rm Fe}_{3}{\rm Ge}{\rm Te}_{2}$ Fei2018 ; Tan2018 ; Deng2018 ; Wang2019 , ${\rm Fe}{\rm P}{\rm S}_{3}$ Wang2016b ; Lee2016 , ${\rm Mn}{\rm Se}_{2}$ Ohara2018 , ${\rm Cr}{\rm Br}_{3}$ Kim2019 ; Ghazaryan2018 ; Zhang2019 , ${\rm Cr}{\rm Cl}_{3}$ Klein2019 , and ${\rm Gd}{\rm Te}_{3}$ Lei2019b .

The atomic nature of these materials is obviously rather complex and a microscopic understanding of their ground-state properties requires massive use of ab initio theories, especially density-functional theory (DFT) Giuliani2005 . In the context of few-layer magnetic crystals, the greatest deal of efforts has been devoted to the study of ${\rm CrI}_{3}$ Wang2018a ; Zhang2015 ; Soriano2018 ; Webster2018 ; Webster2018a ; Jiang2018 ; Zheng2018 ; McGuire2015 ; Wang2011 ; Liu2016 ; Lado2017 ; Wang2016a ; Gudelli2019 ; Xu2018 ; Lee2019 ; Sivadas2018 ; Jang2018 ; Djurdji2018 ; Abramchuk2018 ; Liu2019 ; Torelli2018 ; Lei2019 , while ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ appears to be less studied Gong2017 ; Wang2018 ; Xu2018 ; Li2014b ; Sivadas2015 ; Zhuang2015 ; Li2018 ; Fang2018 ; Kang2019 ; Dong2019 ; Wang2019b ; Sun2019 and rich in controversial findings, as we discuss below.

In this work we therefore focus our attention on ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$. Bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ is a ferromagnetic semiconductor Carteaux1995 ; Siberchicot1996 below a Curie temperature $T_{\rm C}\sim 60~{}{\rm K}$, with an out-of-plane magnetic easy axis and negligible coercivity. Scanning magneto-optical Kerr microscopy has been carried out Gong2017 in bilayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$, revealing ferromagnetic order below $T_{\rm C}\sim 30~{}{\rm K}$. Electric-field control of the magnetization direction has been demonstrated Wang2018 in field-effect transistors made of ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$, with variable thicknesses in the range $3.5$-$20~{}{\rm nm}$. Ballistic Hall micromagnetometry Geim1997 ; novoselov2003 has finally been applied to six-layer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ to study the temperature dependence of the magnetization—see Supplementary Figure S4 of Ref. Kim2019, .

Widely used approximations for the exchange and correlation (xc) energy functional in DFT, mainly based on parametrizations of the xc energy of the homogeneous electron gas Giuliani2005 , miss important features of the physical properties of these magnetic materials. For instance, both the local spin-density approximation (LSDA) and the spin-polarized generalized gradient approximation (spGGA) fail in describing the insulating behavior of systems with strongly localized and thus, correlated, valence electrons. Leaving aside the well-known “gap problem” Giuliani2005 , these approximations often produce a qualitatively wrong ground state. Due to their moderate computational cost and the ability to capture correlation effects, the Hubbard-corrected LSDA (“${\rm LDA}+U$”) and spGGA (“${\rm GGA}+U$”) approximations anisimov_prb_1991 ; anisimov_prb_1993 ; solovyev_prb_1994 ; Liechtenstein1995 have been extensively used.

Many authors Gong2017 ; Wang2018 ; Xu2018 ; Li2014b ; Sivadas2015 ; Li2018 ; Fang2018 ; Kang2019 ; Dong2019 ; Sun2019 studied the electronic and magnetic properties of ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ using the ${\rm LDA}+U$/${\rm GGA}+U$ approximations in order to deal with electrons belonging to the semi-filled $3d$ orbitals of Cr. The main objection we raise to these approaches is that strong hybridization between Cr $d$ and Te $p$ orbitals pushes ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ away from the atomic limit, invalidating the use of Hubbard-type descriptions. Moreover, we recall that arbitrary choices of the Hubbard-$U$ parameter have been criticized in the literature Cococcioni2005 ; Campo2010 ; Kulik2006 ; Kulik2010 , where it has been shown that $U$ is an intrinsic response property which should be calculated self-consistently.

A few calculations Li2014b ; Zhuang2015 ; Fang2018 ; Kang2019 ; Wang2019b have also been carried out using hybrid functionals, which mix the local (or semi-local) DFT exchange with the exact nonlocal Hartree-Fock (HF) exchange, to treat inaccuracies of the LDA/GGA approximations. However, none of these studies compared results obtained by using different approaches, highlighting the criticality of the ${\rm LDA}+U$/${\rm GGA}+U$ approximations for this specific magnetic material.

In this Article we instead present a DFT study of the ground state properties of ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$, comparing the results obtained in the ${\rm GGA}+U$ approximation (in both its non-self-consistent and self-consistent versions) with those obtained via the use of a hybrid functional.

In order to understand the microscopic properties of magnetic materials it is crucial to derive effective spin models—such as the Heisenberg Hamiltonian Gubanov1992 ; Auerbach1994 —from ab initio calculations. Lichtenstein et al. Liechtenstein1987 (LKAG) introduced a rigorous approach to evaluate the exchange couplings of the Heisenberg model from ab initio calculations. These parameters strongly depend on the local environment around a magnetic ion, such as the spin configuration itself, especially in metals and small-gap semiconductors where the spin density is weakly inhomogeneous. In the LKAG approach Liechtenstein1987 , the exchange parameters are determined by studying the total energy variation with respect to small deviations of the magnetic moments of interest from the ground-state magnetic configuration. It relies on the application of the Andersen’s “local force theorems” Machintosh1980 ; Methfessel1982 , which ensure its validity even in the aforementioned complex cases.

Applying the theory of Ref. Liechtenstein1987, in full glory is however computationally demanding and therefore approximate approaches are often used. (For recent work utilizing the LKAG approach in combination with the LSDA and spGGA on bulk ${\rm Cr}{\rm Cl}_{3}$ and ${\rm Cr}{\rm I}_{3}$ see e.g. Ref. Besbes2019, .) These need to be handled with great care. In particular, the exchange coupling parameters and magneto-crystalline anisotropy of few-layer and bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ have been calculated by DFT methods using two of such approximate methods:

i) One, which will be dubbed below “FM-AFM” approach, relies on calculating the difference between the total energies of the ferromagnetic (FM) and anti-ferromagnetic (AFM) configurations. The FM-AFM approach was used by the authors of Refs. Li2014b, ; Zhuang2015, ; Dong2019, to calculate the intra-layer nearest-neighbor exchange couplings. Similarly, the authors of Refs. Sivadas2015, ; Fang2018, used the same approach but calculated exchange couplings up to third-neighbors. The fact that such couplings are important is deduced from inelastic neutron scattering measurements of spin waves in ${\rm MnPS}_{3}$ Wildes1998 and ${\rm FePS}_{3}$ Wildes2012 . The authors of Ref. Fang2018, evaluated also the nearest-neighbor inter-layer exchange coupling in few-layer and bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$.

The evaluation of the exchange couplings via the FM-AFM approach is not justified Liechtenstein1987 when the ground-state spin density is weakly inhomogeneous and the dependence of the exchange couplings on the magnetic configuration is strong, as in the case of metals and small-gap semiconductors.

ii) Another method, dubbed below “four-state mapping analysis” (FSMA) Whangbo2003 ; Xiang2011 ; Xiang2013a , also relies on total-energy calculations but focuses only on the spin configurations of the two ions “linked” by the exchange coupling of interest, leaving all the other ions in the ground-state configuration. The FSMA approach has been used for ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ by the authors of Refs. Gong2017, ; Li2018, ; Xu2018, . Xu et al. Xu2018 used the FSMA approach to study the intra-layer exchange couplings between nearest-neighbor Cr atoms in monolayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$. The same method has been used by the authors of Refs. Gong2017, ; Li2018, to evaluate the intra- and inter-layer exchange couplings up to third-neighbors for bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$. While the FSMA method is a good compromise between the rigor of the LKAG approach and computational feasibility, it is very sensitive to the super-cell size. In Ref. Xu2018, a $2\times 2$ super-cell was used, while in Refs. Gong2017, ; Li2018, super-cells up to dimension $2\sqrt{3}\times 2\times 1$ were used. Below, we demonstrate that such super-cell dimensions do not give converged values of the exchange coupling parameters.

Finally, we show that bilayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ displays two inequivalent ${\rm Cr}$ atoms, a fact that leads to the existence of two distinct intra-layer second-neighbor exchange couplings, $J^{\rm a}_{2}$ and $J^{\rm b}_{2}$. To the best of our knowledge, this has not been noticed in the literature.

This Article is organized as following. In Sect. II we present the geometrical structure of ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ (Sect. II.2) and the theoretical methods (Sect. II.1) that we have used to compute the electronic properties of this material both in the few-layer and bulk cases. In Sect. III we present a summary of our main results on the electronic structure of few-layer and bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$, comparing the ${\rm GGA}+U$ approach with the hybrid-functional one. Exchange couplings and magneto-crystalline anisotropy are reported in Sect. IV. A summary of our main findings and a brief set of conclusions are finally presented in Sect. V.

II Theoretical approach

In this Article, we study the electronic and magnetic properties of the transition metal trichalcogenide ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ in its monolayer, bilayer, and bulk forms.

In Sect. II.1, which can be skipped by uninterested readers, we discuss the computational tools we have used. In Sect. II.2 we present geometric and structural details.

II.1 Computational details

We carry out DFT calculations by using the Quantum Espresso (QE) QE1 ; QE2 and CRYSTAL14 CRYSTAL14 codes, which use plane waves and atom-centered (Gaussian) basis sets, respectively.

For the calculations with QE we use the SG15 Optimized Norm-Conserving Vanderbilt (ONCV) pseudopotentials Hamann2013 ; Schlipf2015 ; Scherpelz2016 , generated by taking into account both the scalar- and fully-relativistic potentials, in order to turn off/on the contribution from spin-orbit coupling (SOC). The pseudopotentials include also the semi-core states and their valence configurations are: $3s^{2},\ 3p^{6},\ 4s^{2},\ 3d^{4}$ for ${\rm Cr}$; $4d^{10},\ 5s^{2},\ 5p^{4}$ for ${\rm Te}$; and $3d^{10},\ 4s^{2},\ 4p^{2}$ for ${\rm Ge}$. We use an energy cutoff up to $70~{}{\rm Ry}$ for all the calculations. For the BZ integrations we employed a Marzari-Vanderbilt smearing Marzari1999 of $10^{-4}~{}{\rm Ry}$ with a Monkhorst-Pack (MP) MP ${\bm{k}}$-point grid with $6\times 6\times 6$ ($12\times 12\times 12$) points for self-consistent calculations of the charge density (density of states) for the bulk crystal, while for the few-layer crystals we adopt a grid with $8\times 8$ ($16\times 16$) points for self-consistent calculations of the charge density (density of states). For the calculations with CRYSTAL14 we used a 86-411d41 Gaussian all-electron basis set Catti1996 with $24$ valence electrons for ${\rm Cr}$ and a double-zeta basis set with an effective core pseudo-potential (m-cc-pVDZ-PP) Heyd2005 with $24$ valence electrons for ${\rm Te}$ and $22$ valence electrons for ${\rm Ge}$. The charge density integrations over the BZ are performed using a uniform $12\times 12\times 12$ ($12\times 12$) MP ${\bm{k}}$-point grid for bulk (few-layer) ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$, respectively. When we use a super-cell, we scale the MP grid size to assure the same accuracy as in the single-cell calculations.

The xc potential is treated in the ${\rm spGGA}$, as parametrized by the Perdew-Burke-Ernzerhof (PBE) formula PBE , with the van der Waals (vdW)-D2 correction proposed by Grimme Grimme2006 . In order to improve the description of the semi-filled $3d$ orbitals of the ${\rm Cr}$ atoms, we adopt both the ${\rm GGA}+U$ scheme Cococcioni2005 ; Kulik2006 ; Kulik2010 and the HSE06 HSE06 hybrid xc functional, which mixes the PBE exchange with $25\%$ of the exact nonlocal HF exchange. This allows us to compare results obtained with different approaches, pointing out an evident criticality of the ${\rm GGA}+U$ method for the material of interest, as discussed below. As anticipated in Sect. I, the Hubbard-$U$ parameter is an intrinsic response property of the material. We calculate it self-consistently in the spin-polarized case, following the linear response approach of Refs. Cococcioni2005, ; Kulik2006, ; Kulik2010, . For bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ we find a self-consistent value of $U$ given by $U_{\rm sc}\equiv 3.9~{}{\rm eV}$. We use this value for all the few-layer crystals under investigation (i.e. monolayer and bilayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$) because we have checked that $U_{\rm sc}$ changes slightly with the number of layers. We also employ $U=U_{\rm sc}$ for the ${\rm GGA}+U$ non-collinear calculations, which include SOC. Finally, we also study the electronic properties, with and without SOC, by varying the value of the Hubbard-$U$ parameter in the range of $0<U\leq U_{\rm sc}=3.9$, as explained below.

In order to evaluate fully-relativistic electronic band structures with the HSE06 hybrid functional we used the Wannier90 W90 code. In fact, it is not possible to calculate them directly from QE due to the heavy computational cost, and SOC is not yet implemented in the CRYSTAL14 code. The Wannier90 code allows us to overcome this problem interpolating the electronic band structure using Maximally-Localized Wannier Functions (MLWFs) Marzari1997 ; Souza2001 extracted from the DFT calculations with QE. As starting guess for the “Wannierisation” procedure, we project the Bloch states onto trial localised atomic-like orbitals: $d$-orbitals for the ${\rm Cr}$ atoms, $p$-orbitals for ${\rm Te}$ atoms, and $p_{z}$-orbital for ${\rm Ge}$ atoms. We made this choice analysing the composition of the density of states around the Fermi energy, as shown in Figs. 2, 3, and 4.

We use CRYSTAL14 to calculate the exchange coupling parameters to reduce the computational cost due to the size of the super-cell and therefore the high number of atoms. This code was especially useful for the calculation of the exchange couplings with the HSE06 hybrid functional. We have been able to switch from QE to CRYSTAL14 after testing the consistency of the results with both codes, comparing electronic band structures and relaxed atomic positions.

We use the VESTA VESTA and Xcrysden Crysden codes to visualize the geometrical structure and BZs, and to produce the plots in Fig. 1.

II.2 Geometrical structure

Bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ forms a layered structure with monolayers separated by a vdW gap—see Fig. 1a). Each monolayer—Fig. 1c)—is formed by edge-sharing ${\rm Cr}{\rm Te}_{6}$ octahedra where the ${\rm Ge}$ pairs are located in the hollow sites formed by the octahedra honeycomb. The layers are ABC-stacked, resulting in a rhombohedral $R\overline{3}$ symmetry (space group number 148), as shown in Fig. 1a), with experimental lattice constants Carteaux1995 ; Siberchicot1996 $a=b=c=7.907$ Å. This structure can also be described as ABC-stacked hexagonal crystal cells with experimental lattice constants $a^{\prime}=b^{\prime}=6.8275$ Å and $c^{\prime}=20.5619$ Å.

Bilayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ is formed by two AB-stacked monolayers—see Fig. 1g).

In Figs 1b) and 1d) we report the bulk and few-layer Brillouin Zones (BZs) Grosso2014 ; Hinuma2017 , respectively.

We adopt the experimental rhombohedral unit cell for the calculation of the bulk properties, while for the monolayer and bilayer we use the hexagonal unit cell shown in Fig. 1c). To simulate the monolayer (bilayer) forms, we consider a super-cell with about $15$ Å ($22$ Å) of vacuum along the $\hat{\bm{z}}$-direction between periodic images.

We optimize the geometrical structures relaxing only the atomic positions until the components of all the forces on the ions are less than $10^{-3}~{}{\rm Ry}/{\rm Bohr}$, while we keep fixed the lattice parameters. Relaxation of the atomic positions has been carried out in the ferromagnetic phase. In Table 1 we report the irriducible set of atomic coordinates and the main atomic distances and angles of both experimental Carteaux1995 ; Siberchicot1996 and optimized structures.

In order to disentangle the effect of the structure from that of magnetism, all calculations below use the atomic coordinates calculated with the spGGA. As can be seen in Table 1, the effect of nonlocal exchange interactions on the crystal structure is, anyway, small.

III Electronic structure

In this Section we present our main results for the electronic band structures of monolayer, bilayer, and bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$.

In Figs. 2, 3, and 4, we report the relativistic band structure of monolayer, bilayer, and bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$, respectively. Results obtained in the ${\rm GGA}+U$ approximation are compared with those obtained with the HSE06 hybrid functional. As anticipated in Sect. II.1, the ${\rm GGA}+U$ results have been obtained by using the self-consistently-calculated value of the Hubbard-$U$ parameter, i.e. $U_{\rm sc}=3.9~{}{\rm eV}$. The opposite qualitative behavior resulting from the two methods is evident. Indeed, the ${\rm GGA}+U$ approximation—see Figs. 2a), 3a), and 4a)—leads to a metallic ground state. In stark contrast, the hybrid functional returns an insulating ground state with an indirect band gap of $E_{\rm g}\sim 0.69~{}{\rm eV}$ for monolayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$—see Fig. 2b)—$E_{\rm g}\sim 0.57~{}{\rm eV}$ for bilayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$—see Fig. 3b)—and $E_{\rm g}\sim 0.43~{}{\rm eV}$ for bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$—see Fig. 4b). The calculated electronic gap of the bulk crystal is comparable with recent experimental results based on angle-resolved photoemission spectroscopy Li2018a , which seem to indicate that bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ at $50~{}{\rm K}$ is a semiconductor with a gap $E^{\rm exp}_{\rm g}$ of at least $0.38~{}{\rm eV}$. A previous experiment Ji2013 based on infrared transmission spectroscopy reported an optical gap of $\sim 0.74~{}{\rm eV}$ at $T=292~{}{\rm K}$. Since this study was carried out at $T\gg T_{\rm C}\sim 60~{}{\rm K}$, it does not help us in choosing the best xc functional to describe the ferromagnetic phase. To the best of our knowledge, Ref. Li2018a, is the only experimental guidance we have at the present stage, as optical studies (sensitive to the gap size) of few-layer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ in the ferromagnetic phase are currently unavailable.

Looking at the few-layer band structures in Figs. 2 and 3, we observe that the valence band maximum (VBM) is at the $\Gamma$ point and the conduction band minimum (CBM) is located on the path that links the $\Gamma$ and ${\rm K}$ high-symmetry points. In the bulk case, instead, we see from Fig. 4 that the CMB lies on the path that links the ${\rm T}$ and ${\rm H}_{0}$ high-symmetry points. The density of states (DOS) projected onto the atomic orbitals is also plotted in Figs. 2-4. Our results for the DOS demonstrate that the main contribution to the electronic bands near the Fermi energy comes from the ${\rm Cr}$ $3d$ and ${\rm Te}$ $5p$ orbitals, pointing out a strong hybridization between the ${\rm Cr}$ and ${\rm Te}$ orbitals. In fact, the ${\rm Te}$ atoms have a fundamental role in stabilizing the ferromagnetic phase of ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ because they mediate super-exchange interactions in the ${\rm Cr}$-${\rm Te}$-${\rm Cr}$ bonds, as per the Goodenough-Kanamori rule Goodenough1955 ; Kanamori1959 . The DOS analysis shows also a small contribution from the ${\rm Ge}$ $4p$ orbitals near the Fermi energy.

In Fig. 5 we show the relativistic band structure obtained in the ${\rm GGA}+U$ approximation, for different values of $0\leq U<U_{\rm sc}$. We note that the gap decreases with increasing $U$. We ascribe this behavior to the strong hybridization between the ${\rm Cr}$ and ${\rm Te}$ atoms, which leads to delocalized $e_{\rm g}$ states. These are not as sensitive as the other $d$ states from the Cr atoms to the application of the Hubbard-$U$, due to their extended rather than atomic-like nature Wu2019 . This is why we have used the HSE06 hybrid functional, which naturally takes into account nonlocal exchange interactions between all orbitals (and not only between the $3d$ orbitals of Cr atoms). By carrying out spin-polarized non-relativistic calculations, we have checked that HSE06 hybrid functional endows the VBM and CBM with spin-up character, differently from what happens in the ${\rm spGGA}$. We conclude that the application of a scissor operator to open the gap leads to wrong electronic and optical properties Fang2018 . The calculated scalar-relativistic electronic gap for bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ is $E_{\rm g}\sim 0.73~{}{\rm eV}$, which is, as expected, larger than the one calculated with the inclusion of SOC. Finally, we note that the hybridization between the ${\rm Cr}$ $3d$ and ${\rm Te}$ $5p$ orbitals, which leads to the failure of the ${\rm GGA}+U$ approximation, is less evident in few-layer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ than in the bulk crystal.

We now discuss about magnetism. With the HSE06 functional and in bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ the ferromagnetic state has lower energy than the paramagnetic one, the energy difference being of $\sim 3.5~{}{\rm eV}/{\rm Cr}$. Moreover, the ferromagnetic state is more stable than all the other antiferromagnetic configurations by at least $\sim 0.5~{}{\rm meV}/{\rm Cr}$ in the same approximation. In the ${\rm GGA}+U$ approximation, the magnetic properties of bilayer and bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$, especially the interlayer exchange couplings, vary strongly with $U$. In fact, for small values of $U$, the ground state remains ferromagnetic, while for $U\gtrsim 1.5~{}{\rm eV}$ bilayer and bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ become antiferromagnetic, with oppositely-oriented magnetic moments in alternating layers Kang2019 . This is in stark contrast with the experimental observations Gong2017 . Monolayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ remains in a ferromagnetic ground state for $1.5~{}{\rm eV}\lesssim U\leq U_{\rm sc}$, despite its exchange couplings greatly vary with $U$. In Fig. 1h) we report a pictorial representation of the magnetic moments, while in Table 2 we summarize the values of the magnetic moments, as obtained by carrying out spin-polarized calculations—with the HSE06 functional and in the ${\rm GGA}+U_{\rm sc}$ approximation—for all the systems under study. The magnetic moments change slightly with the inclusion of SOC. As we will discuss in Sect. IV, in bilayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ we observe two inequivalent ${\rm Cr}$ atoms in the unit cell, in one layer.

IV Exchange couplings and magneto-crystalline anisotropy

According to the Mermin-Wagner theorem Mermin1966 , thermal fluctuations in 2D systems imply the impossibility of long-range magnetic order. A way to bypass the theorem and explain the discovery of 2D magnetic materials Gong2017 ; Huang2017 is to recognize the presence of a large magneto-crystalline anisotropy. The latter lifts the invariance under rotations, an hypothesis of the Mermin-Wagner theorem.

Electron spin resonance and ferromagnetic resonance measurements have been carried out in bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ over a wide frequency and temperature range Zhang2016 ; Zeisner2018 ; Liu2019a ; Khan2019 , revealing the development of 2D spin correlations in the vicinity of $T_{\rm C}$ and proving the intrinsically low-dimensional character of spin dynamics in this material. Such experiments returned a MAE of $\sim 50$-$90~{}\mu{\rm eV}/{\rm Cr}$, which is on the same order of what estimated in Ref. Gong2017, for bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$.

In order to understand the microscopic properties of magnetic materials it is crucial to derive an effective spin model. In the case of localized magnetic moments, the Heisenberg spin Hamiltonian can be used. The total Hamiltonian is the sum a non-magnetic part and a magnetic part, which has the form

[TABLE]

The first term on the right-hand side of Eq. (1) is the rotationally-invariant Heisenberg Hamiltonian with exchange couplings $J_{ij}$ and the sums over $i<j$ run over all the ${\rm Cr}$ pairs without double counting. The second term is the magneto-crystalline anisotropy energy (MAE), with the sum over $i$ running over all the ${\rm Cr}$ sites. The Dzyaloshinski-Moriya interaction term, $\sum_{i<j}{\bm{D}}_{ij}\cdot(\hat{\bm{S}}_{i}\times\hat{\bm{S}}_{j})$, is absent in Eq. (1) due to the existence of an inversion center between the ${\rm Cr}$ ions Moriya1960 .

The MAE originates from SOC Bruno1989 ; Li2014 , and is therefore a quantum effect of relativistic nature. It breaks the rotational invariance with respect to the spin quantization axis and is determined by the interaction between the orbital state of a magnetic ion and the surrounding crystalline field. It can be calculated Gong2017 ; Li2018 ; Li2014 by looking at the total energy difference—obtained through self-consistent calculations in the presence of SOC—between the configuration with all spins perpendicular to the layer (along the $\hat{\bm{z}}$ direction) and that with all spins parallel to the layer (along the $\hat{\bm{x}}$ or $\hat{\bm{y}}$ directions). A positive (negative) sign of the MAE means that the system is an easy-plane (easy-axis) ferromagnet.

As anticipated in Sect. I, an accurate method to calculate the exchange couplings for small-gap semiconductors and metals is the FSMA approach Whangbo2003 ; Xiang2011 ; Xiang2013a . The FSMA approach considers one specific magnetic pair at a time in a super-cell. Without loss of generality, the energy can be written as

[TABLE]

where $E_{\rm spin}$ and $E_{0}$ represent the magnetic and non-magnetic energies, respectively. Here, $J_{12}\bm{S}_{1}\cdot\bm{S}_{2}$ is the exchange coupling between sites $1$ and $2$, while $\bm{S}_{i}\cdot\bm{K}_{i}$ represents the coupling between site $i$ and the all the magnetic sites different from $1$ and $2$. Finally, $E_{\rm others}$ represents the contribution to the energy stemming from the interaction between sites $1$ and $2$ and all the other non-magnetic sites. Consider the following collinear spin states:

[TABLE]

The other spin states are kept fixed, according to the ground-state spin configuration. Using Eq. (2), we can write the energies of the four states as following:

[TABLE]

The solution of this system leads to the equality

[TABLE]

The FSMA method is accurate when the super-cell used for the calculations is large enough. We monitor the accuracy by using a convergence test with respect to the super-cell size. In Fig. 6 we show the convergence of the intra-layer exchange coupling parameters $J_{i}$, depicted in Fig. 1e), with the size of the super-cell for monolayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$. We clearly see that we need at least a $3\times 3$ ($3\times 3\times 1$) supercell for monolayer (bulk) ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ to have converged results. This implies that the results reported in Refs. Gong2017, ; Xu2018, ; Li2018, for the exchange couplings are not converged. In particular, we see that the coupling $J_{2}$ between second neighbors changes character from anti-ferromagnetic to ferromagnetic if we choose a $3\times 3$ supercell instead of a $2\times 2$ or a $2\sqrt{3}\times 2$ supercell, see Fig. 6b).

In Table 3 we present a summary of our results (in black) and the comparison with existing calculations (in red). All the exchange couplings have been calculated with the FSMA approach. The only exception is the row in Table 3 labeled by “${\rm spGGA}^{\prime}$”. Results in this row have been calculated with the FM-AFM method (see Sect. I) with the only aim of pointing out how much these differ from the ones calculated with the FSMA, in the same approximation (spGGA) for the xc functional. Note that, for monolayer and bilayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$, the sign of the MAE calculated in the spGGA (${\rm MAE}>0$) is opposite to that calculated with the HSE06 hybrid functional (${\rm MAE}<0$). For all the calculations and the comparisons with the literature we use $S=3/2$, as justified by the atomic magnetic moments reported in Table 2.

Following Ref. Torelli2018, and taking into account only the largest exchange coupling (i.e. $J_{1}$) and the MAE, we estimate the Curie temperature $T_{\rm C}$. We find: $T_{\rm C}\sim 34~{}{\rm K}$ for monolayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$; $T_{\rm C}\sim 40~{}{\rm K}$ for bilayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$; and $T_{\rm C}\sim 48~{}{\rm K}$ for bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$. The trend of $T_{\rm C}$ with respect to the number of layer is in qualitative agreement with the experimental results of Ref. Gong2017, . The results for $J_{i}$, MAE, and $T_{\rm C}$ emphasize a strong dimensionality effect, in agreement with the experimental behaviors Gong2017 .

We finally note that bilayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ presents a particular magnetic structure resulting from the lack of the ABC-stacking configuration and periodicity along the $\hat{\bm{z}}$ direction. Indeed, in AB-stacked bilayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ we can identify two inequivalent Cr atoms in the unit cell—see Fig. 1g). With reference to this figure, we clearly see that in the top layer (orange) there are Cr atoms which have another Cr atom underneath in the bottom layer (black)—these Cr atoms are label by the letter “a” in blue in Fig. 1g)—while there are other Cr atoms, which do not have a Cr partner underneath in the bottom layer—these Cr atoms are label by the letter “b” in green in Fig. 1g). In Table 3 we report the calculated values of the two associated second-neighbor exchange couplings, $J^{\rm a}_{2}\neq J^{\rm b}_{2}$. Note that the inter-layer couplings affect the interaction between these two inequivalent Cr atoms. This is evident if we compare the results obtained with the spGGA and HSE06 functionals. In fact, the spGGA overestimates the inter-layer couplings, resulting in a vanishing $J^{\rm b}_{2}$. In contrast, smaller inter-layer couplings, extracted from the HSE06 functional, lead to two different but comparable second-neighbor exchange couplings $J^{\rm a}_{2}$ and $J^{\rm b}_{2}$. In the bulk crystal the equivalence between the two Cr atoms is restored and we find $J^{\rm a}_{2}=J^{\rm b}_{2}=J_{2}$.

V Summary and conclusions

In summary, in this Article we have reported on a comparative density functional theory study of the electronic and magnetic properties of the recently discovered atomically-thin magnetic material ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$.

We have clearly shown that the often-used ${\rm LDA}+U$ and ${\rm GGA}+U$ approaches fail in predicting the ground-state properties of this material in both its monolayer and bilayer forms, and even more spectacularly in its bulk form. According to these approaches, ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ should have a metallic ground state, with its magnetic properties strongly depending on the number of layers and the values of the Hubbard-$U$ parameter used, including the self-consistently-calculated one, $U_{\rm sc}=3.9~{}{\rm eV}$ for the bulk crystal. On the contrary, the use of the HSE06 hybrid functional, which naturally takes into account nonlocal exchange interactions between all orbitals (and not only between the $d$ orbitals of Cr), yields an insulating ferromagnetic ground state—see Figs. 2-4—and an indirect gap, which, for the case of bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$, is comparable with the one recently measured Li2018a via angle-resolved photo-emission spectroscopy.

Finally, we have calculated all the relevant exchange couplings for monolayer, bilayer, and bulk ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$, comparing results obtained via different methods and/or super-cell size—see Table 3. Using the four-state mapping analysis approach Whangbo2003 ; Xiang2011 ; Xiang2013a , we have shown that the slow convergence of the exchange couplings with the super-cell size may be responsible for discrepancies in the literature Gong2017 ; Li2018 ; Xu2018 . Such a convergence check has been reported in Fig. 6.

In the case of bilayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$, we have shown that two distinct intra-layer second-neighbor exchange couplings emerge—see the portion of Table 3 that concerns bilayer ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$.

In the future, it would be interesting to study the dependence of these results on doping, obtained via the usual electrical field effect, plasmons, and plasmon-magnon coupling in ${\rm Cr}_{2}{\rm Ge}_{2}{\rm Te}_{6}$ and other magnetic 2D materials.

Acknowledgements.

This work was supported by the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 785219 - GrapheneCore2. G.M. was supported by the EDISON-Volta prize (2018). We thank L. Paulatto for help with the open_grid package of Quantum Espresso and C. Tresca and M.I. Katsnelson for useful discussions. Finally, we acknowledge the allocation of computer resources from CINECA, through the “ISCRA C” projects “HP10CXQECI” and “HP10CXRX1O”, PRACE (Project No. 2017174186), and EDARI (Grant A0050901202).