D V Rybkovskiy

A M Prokhorov

I C Gerber

M V Durnev

Atomically inspired k • p approach and valley Zeeman effect in transition metal dichalcogenide monolayers

Keywords: numbers: 73.20.-r, 73.21.Fg, 73.63.Hs, 78.67.De

We developed a six-band k • p model that describes the electronic states of monolayer transition metal dichalcogenides (TMDCs) in K-valleys. The set of parameters for the k • p model is uniquely determined by decomposing tight-binding (TB) models in the vicinity of K ± -points. First, we used TB models existing in literature to derive systematic parametrizations for different materials, including MoS2, WS2, MoSe2 and WSe2. Then, by using the derived six-band k • p Hamiltonian we calculated effective masses, Landau levels, and the effective exciton g-factor g X 0 in different TMDCs. We showed that TB parameterizations existing in literature result in small absolute values of g X 0 , which are far from the experimentally measured g X 0 ≈ -4. To further investigate this issue we derived two additional sets of k • p parameters by developing our own TB parameterizations based on simultaneous fitting of ab-initio calculated, within the density functional (DFT) and GW approaches, energy dispersion and the value of g X 0 . We showed that the change in TB parameters, which only slightly affects the dispersion of higher conduction and deep valence bands, may result in a significant increase of |g X 0 |, yielding close-to-experiment values of g X 0 . Such a high parameter sensitivity of g X 0 opens a way to further improvement of DFT and TB models.

I. INTRODUCTION

Monolayers of transition metal dichalcogenides (TMDCs) have attracted recently much attention due to their exceptional properties, such as coupling of spin and valley degrees of freedom, which allows for the valley polarization with a circularly polarized light in these materials [1][2][3][4][5][6].

Recent magneto-photoluminescence experiments revealed the significant Zeeman splitting of emission lines associated with optical recombination of excitons and trions in different valleys. The effective exciton g-factors of this "valley Zeeman effect" in a magnetic field directed normal to a monolayer plane were found to be g X 0 ≈ -4 for a wide range of investigated materials, including selenides [7][8][9][10][11][12], sulphides [13,14] and tellurides [15].

The multi band k •p method is perfectly suited for theoretical investigation of magneto-optical and transport effects, including Zeeman effect, optical absorption and photogalvanics [START_REF] Ivchenko | Superlattices and other heterostructures[END_REF][START_REF] Voon | The kp Method: Electronic Properties of Semiconductors[END_REF]. The available two-band k•p models describing electronic spectra in TMDCs [START_REF] Kormányos | [END_REF][19][20] account for the bottom conduction and topmost valence bands and are parameterized by density functional theory (DFT) calculations. The simple two-band model is, however, insufficient for calculation of exciton Zeeman effect since the exciton g-factor g X 0 vanishes in the * correspondence address: durnev@mail.ioffe.ru two-band k • p approximation [11]. Therefore to obtain nonzero g X 0 additional terms describing contributions of remote bands should be included in the two-band k • p model [19,21], which make the model less transparent. Moreover parametrization of the k • p model by fitting the DFT band structure only in the vicinity of distinct points of the Brillouin zone is less reliable and allows for a much freedom in the choice of parameters.

In this work we use an alternative approach based on the idea proposed in Ref. [11]. This approach comprises three steps, namely, (i) use DFT-based calculations of electronic states in TMDCs as a starting point, (ii) use atomistic tight-binding (TB) model to fit the electronic spectrum and wave functions along the high-symmetry paths of the Brillouin zone, and (iii) derive k • p Hamiltonian by decomposing the TB model in the vicinity of a given wave vector in the Brillouin zone. As a result, we obtain a multi band k • p Hamiltonian (the number of bands is equal to the number of atomic orbitals included in the TB model) with a set of parameters that are uniquely determined by the TB parametrization. This Hamiltonian can be further used for calculating the valley Zeeman effect, Landau levels, etc.

As a starting point we use several existing eleven-band TB models [22][23][24][25][26]. These models include d-type orbitals of metal atoms and p-type orbitals of chalcogen atoms, and capture all symmetries of the studied system. The resulting k • p Hamiltonian, which describes dispersion of the bottom conduction and topmost va-lence bands, contains six bands that have even parity with respect to mirror reflection in the monolayer plane. We present a systematic parametrization of the six-band k • p Hamiltonian by expansion of TB Hamiltonians [22][23][24][25][26] in the vicinity of K ± points of the Brillouin zone for different materials (MoS 2 , WS 2 , MoSe 2 , WSe 2 ) and different TB models. We also do the whole three-step procedure on our own, i.e., perform post-DFT calculations by applying GW formalism to obtain reliable band gaps and accurate band dispersion, fit it with the TB model and do the k • p decomposition, which allows us, as a result, to obtain an independent k • p parameterization. We want to stress the novelty of using GW set of data to extract TB parameters since in previous studies the main focus was made on DFT calculations using mainly local or semi-local exchange-correlation functionals [22][23][24], whereas hybrid functional was also used [25] to partially correct the severe underestimation of band gap values usually observed at the DFT level for MoS 2 monolayers, see Ref. [27] for a recent review. The use of GW correction strongly affects the effective charge carrier masses [28,29] too.

Then, we use the derived six-band k • p Hamiltonian to calculate effective masses, Landau levels, and the effective exciton g-factor g X 0 in different TMDCs. We analyze the main contributions to g X 0 , which result from the mixing with excited conduction and deep valence bands. We show that the calculated values of g X 0 and effective masses vary in a wide range for different TB parameterizations. The calculated values of g X 0 for available in literature TB models (|g X 0 | 1) are rather far from experimental values. However using our DFT+GW calculations and the TB fitting procedure we were able to find additional k • p parameterization sets, which well describe the experimental exciton g-factor (g X 0 ≈ -4), as well as conduction and valence band effective masses and the wave functions coefficients. This result underlies the importance of the careful choice for the DFT starting point calculations and TB parameterizations: along with effective masses and energy gaps the g-factor value may serve as a test for improving both the DFT calculations and TB models.

II. EFFECTIVE KP-HAMILTONIAN

As a starting point for construction of an effective k•p-Hamiltonian we will use eleven-band tight-binding models developed in Refs. 22-26. These tight-binding models include three p-type orbitals on each of the two chalcogen atoms (X) and five d-type orbitals on a metal atom (M). The electron wave function within the tight-binding approximation is presented as a linear combination of atomic orbitals φ (a) j [11]

Ψ (n) k (r) = a,l,j e ikR a,l C (a) j φ (a) j (r -R a,l ) , (1)
where n is a number of an electronic band, k is a wave vector, a = M, X denotes the type of an atom, l runs through the atoms of a given type, j enumerates the set of orbitals, R a,l gives the position of atoms in a twodimensional lattice, and C

(a) j are coefficients. The basis orbitals of the eleven-band tight-binding model are [22][23][24][25][26]

φ (a) j = {d z 2 , d + , d -, p + , p -, p z,A , d xz , d yz , p x,A , p y,A , p z,S } , (2)
where d ± = d x 2 -y 2 ± 2id xy , p ± = p x,S ± ip y,S , d α denotes the orbital with a d-like symmetry of the M atom and

p β,S = (p β,t + p β,b )/ √ 2, p β,A = (p β,t -p β,b )/ √
2 are the symmetric and asymmetric combinations of the p-type orbitals of the top (t) and bottom (b) X atoms in the unit cell, x and y axes lie in the monolayer plane, and z is the monolayer normal. The spin orbit interaction between electron spin and orbital momenta of atomic orbitals [23,24] is neglected in this work.

If the z → -z mirror symmetry is conserved, i.e. for a free-standing monolayer in the absence of external electric field and strain, the Hamiltonian, which describes the energy spectrum of a monolayer electron with a wave vector k, written in the basis Eq. ( 2) has the form

H(k) = H E 0 0 H O . (3)
Here H E is the 6×6 block acting on the orbitals with even with respect to z → -z symmetry, and H O is the 5×5 block acting on the orbitals with odd symmetry. The exact form of blocks H E and H O depends on a particular tight-binding model [22][23][24][25][26]. It is known that the Bloch functions of the bottom conduction and topmost valence bands, which are of the main interest in this work, are even with respect to z → -z reflection [START_REF] Kormányos | [END_REF]20], and therefore these bands are described by the H E block. We note that magnetic field normal to a monolayer does not break the parity of wave functions, and hence we do not need the H O block in the calculation of g X 0 . We also note, that the mixing of the H E and H O blocks by a perturbation that breaks z → -z symmetry does not affect g X 0 in the first order in this perturbation.

To construct an effective k•p Hamiltonian in the vicinity of K ± = (±4π/3a 0 , 0) points, where a 0 is the lattice constant, we will decompose the tight-binding Hamiltonian H E (k) over a small wave vector q = k -K ± . Up to the second-order terms this decomposition yields

H ± E (q) ≈ H E (K ± ) + α=x,y ∂H E ∂k α (K ± )q α + + 1 2 α,β=x,y ∂ 2 H E ∂k α ∂k β (K ± )q α q β . (4)
Electron wave functions at K ± valleys of MX 2 transform according to irreducible representations (irreps) of the C 3h point group. We denote six wave functions that diagonalize H + E (q) at q = 0 as Ψ

(v-5) E 1 , Ψ (v-4) A , Ψ (v-3) E 2 , Ψ (v) A , Ψ (c) E 1 , Ψ (c+2) E 2
, where a superscript names the electronic band and a subscript denotes the corresponding irreducible representation (A , E 1 and E 2 ), see Tab. I. Note that in addition to two conduction (c and c + 2) and two valence (v and v -3) bands, known from the four-band k • p models [START_REF] Kormányos | [END_REF]20], the six-band model contains two deep valence bands v -4 and v -5, which transform at K + -point via A and E 1 representations, respectively.

In the new basis H ± E (q) can be written as

H ± E (q) = H ± 1 (q) + H ± 2 (q) , (5)
where the first term contains linear in q terms:

H + 1 (q) =        E v-5 δ 7 q -δ 6 q + δ 4 q -0 δ 2 q + δ 7 q + E v-4 δ 5 q - 0 δ 3 q + δ 1 q - δ 6 q -δ 5 q + E v-3 γ 2 q + γ 5 q -0 δ 4 q + 0 γ 2 q -E v γ 3 q + γ 4 q - 0 δ 3 q -γ 5 q + γ 3 q -E c γ 6 q + δ 2 q -δ 1 q + 0 γ 4 q + γ 6 q -E c+2        , (6)
and the second one

H ± 2 (q) nl = 2 q 2 2m n δ nl , n, l = 1..6 (7)
is a diagonal matrix with quadratic in q elements. In Eqs. ( 6), (7) E n (n = c + 2, c . . . ) are the band energies at K ± -points, γ j and δ j are parameters, q ± = q x ± iq y , and q 2 = q 2 x + q 2 y . The effective masses m n describe contributions to the band dispersion arising from the mixing with remote bands, which are not present in the k•p model [11]. Note that in the decomposition of the offdiagonal elements of the tight-binding Hamiltonian (4) we retained only the linear in q terms. The phases of wave functions at K + -point (Tab. I) are chosen in such a way that parameters γ j and δ j are real. Parameters of the k • p Hamiltonians ( 6) and ( 7) for different parametrizations, materials and tight-binding models are listed in Tab. III of the Appendix A.

To derive the k • p Hamiltonian H - 1 at the K --point of the Brillouin zone one should replace q + by q -and vice versa in Eq. ( 6) [30]. Note that at k = K -basis wave functions Ψ (c) and Ψ (v-5) transform according to E 2 irreducible representation whereas Ψ (c+2) and Ψ (v-3) transform according to E 1 .

The effective masses of the main conduction and valence bands in the framework of k • p model are

1 m c = 1 m c + 1 m * c , 1 m v = 1 m v + 1 m * v , (8)
where

1 m * c = 2 2 γ 2 5 E c -E v-3 + γ 2 3 E c -E v + (9) + γ 2 6 E c -E c+2 + δ 2 3 E c -E v-4 , 1 m * v = 2 2 γ 2 2 E v -E v-3 + γ 2 3 E v -E c + (10) + γ 2 4 E v -E c+2 + δ 2 4 E v -E v-5
result from the mixing of electronic bands described by Hamiltonian (6), whereas m c and m v account for the mixing with remote bands, see Eq. ( 7).

Figure 1 shows the dispersion of electronic bands in MoS 2 calculated in the framework of TB models of Refs. [24,26], respectively, in panels (b) and (a) and the dispersion ε n = E n + 2 q 2 /2m n , where an effective mass of the n-th band is calculated similar to n = c and n = v, see Eqs. ( 8), ( 9), (10). Figure 1 illustrates that a certain care should be taken when labeling the valence bands: the bands are labeled according to its wave functions representations, see Tab. I, and the order might be different for different TB models. The effective masses m c and m v for different TB models are listed in Tab. II.

III. ZEEMAN EFFECT

In this section we use the developed k • p model to calculate the g-factors of electrons in conduction and valence bands. The main interest, however, is related to the exciton g-factor, which has been measured in a number of recent experiments by optical means. Single carrier Zeeman splittings can be determined, for instance, from the measurements of Shubnikov-de Haas oscillations.

We consider the Zeeman splitting of electrons in K ± valleys in magnetic field B = (0, 0, B z ) directed normal to a monolayer plane. The Zeeman effect contains spin and valley contributions described by g-factors g 0

TABLE I. The nonzero coefficients C

(a) j of wave functions Eq. ( 1) at K + -point and corresponding irreducible representations of the C 3h point group. The phases of the wave functions are chosen in such a way, that αi and βi are real numbers, α 2 i + β 2 i = 1, and αi > 0. We denote pz ≡ pz,A.

Irrep Band Nonzero wave function coefficients

A v, v -4 Ψ (v) : C d + = α1, Cp + = iβ1; Ψ (v-4) : C d + = β1, Cp + = -iα1 E 1 c, v -5 Ψ (c) : C d z 2 = α2, Cp -= iβ2; Ψ (v-5) : C d z 2 = β2, Cp -= -iα2 E 2 c + 2, v -3 Ψ (c+2) : C d -= α3, Cp z = β3; Ψ (v-3) : C d -= β3, Cp z = -α3
Γ M K Γ TB k•p c+2 c v v-5 v-3 v-4 -6 -4 -2 0 2 Γ M K Γ c+2 c v v-5 v-3 v-4
Energy (eV) - Electronic spectra of MoS2 calculated using TB models H. Rostami et al. [26] (a) and S. Fang et al. [24] (b). Red lines are k•p quadratic dispersions at K-point calculated using effective masses for each band (see text for details). Note the different order of deep valence bands in two panels. and g orb , respectively

H B = g 0 2 µ B B z σ z + g orb 2 µ B B z τ z . (11)
Here σ z is a spin operator (σ z = ±1 for spin-up and spindown electrons, respectively), and τ z represents the valley degree of freedom (τ z = ±1 for K + and K -electrons, respectively), µ B is the Bohr magneton. The effective g-factors of K + and K -electrons are defined as [11] g

K+ c,v ≡ g c,v = g 0 + g c,v orb , g K+ c,v = -g K- c,v . (12)
The valley term g c,v orb has an orbital nature and accounts for the mixing of the electronic bands by magnetic field. Within the k • p scheme this mixing is obtained by replacing q in Eq. ( 6) with q -(e/c )A, where e = -|e| is the electron charge, and A is the vector potential of the magnetic field, resulting in (cf. Ref. [11])

g c orb = 4m 0 2 - γ 2 5 E c -E v-3 + γ 2 3 E c -E v - (13)
-

γ 2 6 E c -E c+2 + δ 2 3 E c -E v-4 , g v orb = 4m 0 2 γ 2 2 E v -E v-3 - γ 2 3 E v -E c + (14)
+ γ 2 4 E v -E c+2 - δ 2 4 E v -E v-5 .
The spin g-factor g 0 comprises two contributions, namely, the bare electron g-factor (g 0 = 2) and a small contribution due to the spin-orbit interaction, which is not taken into account in our model. This contribution within the k • p model is of the order of ∼ g orb ∆ so /∆E ij g orb , where ∆ so is the spin-orbit splitting of a given band, and ∆E ij is a characteristic energy distance to other bands.

The σ + and σ -photoluminescence lines observed in experiment originate from the radiative recombination of neutral excitons X 0 with electrons occupying K + and K -valleys, respectively [31]. Therefore the effective Zeeman splitting of X 0 is ∆ Z = g X 0 µ B B z with [11]

g X 0 = g c -g v . (15)
In this difference, according to Eqs. ( 12), ( 13), ( 14), the contribution to g c and g v , which occurs due to the mixing between c and v bands (∝ γ 2 3 ), cancels out, so that nonzero contributions to g X 0 arise due to the mixing of v and c with deep valence and excited conduction bands.

IV. LANDAU LEVELS

The developed k•p model allows for calculation of Landau levels in the system. For this purpose we make replacements q + → √ 2a † /l B and q -→ √ 2a/l B in Eqs. ( 6) and (7), where a † and a are the creation and annihilation operators acting in the space of Landau functions, l B = |e| /|B z |c, and decompose the six-component wave function of the j-th band Ψ (j) in a series of the Landau level functions ϕ n,qy [32]

Ψ (j) = n≥0          a n b n c n d n e n f n          ϕ n,qy , (16)
where n and q y are the quantum numbers, and a n , b n , c n , d n , e n , and f n are coefficients. The numeric diagonalization of obtained Hamiltonian yields the energy position of the j-th band in magnetic field E j (n, B z ). However this energy contains also the valley Zeeman term, described by the second term of Eq. ( 11). To get rid of the Zeeman term we define Landau levels as

ε j (n, B z ) = 1/2[E j (n, B z ) + E j (n, -B z )].
The dispersion of the first four Landau levels for j = c and j = v is presented in Fig. 2. For comparison we also show linear dispersions calculated using the simple formula

ε j (n, B z ) = ω (j) c (n + 1/2), with ω (j) c = |eB z |/m j c
and the effective masses m c and m v given by Eq. ( 8). One can see that the results given by the exact numeric diagonalization of the effective Hamiltonian and the simple analytical formula coincide in the wide range of B z , the discrepancy in the valence band is more noticeable due to more pronounced band non parabolicity.

V. DISCUSSION

Table II presents the values of g-factors and effective masses calculated within the developed k • p model after Eqs. ( 13), (14) and Eqs. ( 8), ( 9), (10). Since the large contributions to g c and g v that originate from the mixing between v and c bands cancel out in the exciton g-factor, the value of g X 0 is defined by the mixing with deep valence and excited conduction bands. The main contributions to g X 0 in the studied parameterizations come from the mixing with v -3 and c+2 bands, i.e. from the terms -γ 2 5 /(E c -E v-3 ) and -γ 2 6 /(E c -E c+2 ) in Eq. ( 13) and terms γ 2 2 /(E v -E v-3 ) and γ 2 4 /(E v -E c+2 ) in Eq. ( 14). As an example, the contribution from c + 2 gives ≈ 15 % of the total g c value, and the contribution from v-3 gives ≈ 30 % of the total g v value for MoS 2 parametrization (a) in Tab. II. As seen from Tab. II, existing in literature TB models result in small absolute values of g X 0 , which are far from the experimentally measured g X 0 ≈ -4.

Based on our own DFT and post-DFT (GW ) calculations and TB fitting procedure (see computational details and dispersion of energy bands in Appendices B, C) we obtained two additional k•p parameterizations for MoS Solid lines show the results of numeric calculations using decomposition Eq. ( 16), and dashed lines stand for a simple analytical formula εj(n, Bz) = ω we used two different TB models, the eleven-parameters TB model of Ref. [26] and the TB model of Ref. [24] with twenty five independent parameters. Within the fitting procedure we numerically extracted parameters of the k • p Hamiltonian and calculated exciton g-factor using Eqs. ( 13), (14). We then used the value of g X 0 as an extra fitting parameter (we fit it to the experimental value g X 0 ≈ -4) additional to band dispersions and wave function coefficients.

The fitting procedure for the TB model of Ref. [26] results in a good fit of the dispersion of c and v bands across high-symmetry paths of the Brillouin zone and only a slight change of energy position and dispersion of high conduction and deep valence bands compared to the original parameterization of Ref. [26]. However this change is sufficient to obtain a large increase of |g X 0 |, g X 0 ≈ -1.75 (see column (d) of Tab. II). Using the TB model of Ref. [24] we were able to obtain g X 0 ≈ -3.82 as well as a good fit for all six energy bands dispersions, wave function coefficients and effective masses, see column (e) in Tab. II.

The wide spread of calculated g X 0 values underlies the sensitivity of g X 0 to a given parametrization of a DFT or a TB model. Hence, along with effective masses and energy gaps the value of g X 0 may serve as a test tool for a given parametrization of a DFT or a TB model. The TB parameterizations that fit g X 0 can be obtained in principle for other materials (e.g. MoSe 2 , WS 2 , and WSe 2 ), however this is beyond the scope of the present work.

So far in our theory we treated electron and hole in the exciton independently, neglecting the Coulomb interaction between charge carriers. It is well known, however, that the exciton binding is large in TMDCs and plays a significant role in optical experiments [33,34]. The effects of Coulomb interaction between an electron and a hole as well as localization by an in-plane potential [35] can be estimated by introducing the g-factor dependence on a charge carrier energy. Such a dependence is derived by simply replacing E c with E c + ∆E c in Eq. ( 13) and E v with E v -∆E v in Eq. ( 14), where ∆E c and ∆E v are the energy shifts of electron and hole levels [START_REF] Ivchenko | Optical spectroscopy of semiconductor nanostructures[END_REF] due to either localization or Coulomb binding. This dependence is depicted in Fig. 3. Note that negative energies ∆E c , ∆E v reflect the binding of electron and hole in an exciton, whereas positive ∆E c , ∆E v correspond to localization of a charge carrier in a quantum dot. One can see that within a typical scale of exciton binding energy in TMDCs, E B ∼ 400 meV, g c and g v change significantly, resulting in a possible enhancement ∆|g X 0 | ≈ 1.

VI. CONCLUSIONS

To conclude, we developed a six-band k • p model that describes the electronic states of monolayer TMDCs in K-valleys. The set of parameters for the k • p model is uniquely determined by decomposing eleven-band tight-binding models in the vicinity of K ± -points. Using existing in literature TB models we were able to derive systematic parametrizations for different materials (MoS 2 , WS 2 , MoSe 2 , WSe 2 ) and different TB Hamiltonians. Using the derived six-band k • p Hamiltonian we calculated effective masses, Landau levels, and the effective exciton g-factor g X 0 in different TMDCs. We showed that the main contributions to g X 0 result from the mixing with excited conduction band c + 2 and deep valence band v -3. We also obtained two additional sets of k • p parameters by developing our own TB parameterizations based on simultaneous fitting of ab-initio calculated energy dispersion and the value of g X 0 .

The k•p parameterizations extracted from existing TB models result in small absolute values of g X 0 , which are far from the experimentally measured g X 0 ≈ -4. However as we showed using our additional sets of k • p parameters, the change in parameters, which only slightly affects the dispersion of higher conduction and deep valence bands, may result in a significant increase of |g X 0 |. As a result, we obtained g X 0 ≈ -1.75 and g X 0 ≈ -3.82 for the two sets. Such a high parameter sensitivity of g X 0 opens a way to further improvement of DFT and TB models, since g-factor modeling requires at the same time an accurate description of deep valence and high conduction bands. 6), (7). The values of γi and δi are given in eV Å, the units of Ei are eV, the effective masses m n are given in the units of m0. Parameterizations are based on TB models listed in footnotes. f TB model of Ref. [24] for DFT calculations g TB model of Ref. [24] for DFT calculations h TB model of Ref. [24] for DFT calculations software uses the plane-augmented wave scheme [40,41] to treat core electrons. Fourteen electrons for Mo, W atoms and six for S, Se ones are explicitly included in the valence states. All atoms are allowed to relax with a force convergence criterion below 0.005 eV/ Å. A grid of 12×12×1 k-points has been used, in conjunction with a vacuum height of 17 Å, to take benefit of error's cancellation in the band gap estimates [42]. A gaussian smearing with a width of 0.05 eV was used for partial occupancies, when a tight electronic minimization tolerance of 10 -8 eV was set to determine with a good precision the corresponding derivative of the orbitals with respect to k needed in quasi-particle band structure calculations. Spin-orbit coupling was not included to determine eigenvalues and wave functions as input for the full-frequencydependent GW calculations [43] performed at the G 0 W 0 level. The total number of states included in the GW procedure was set to 600, after a careful check of the direct band gap convergence, to be smaller than 0.1 eV.

MoS2

Appendix C: Details of TB fitting procedure and additional TB parameterizations

Most of the modern TB parametrizations are made to reproduce the energy bands of ab-initio calculations. The parameter set is usually found by minimizing the function:

f ({t i }) = i,k w i,k E TB i,k -E GW i,k 2 , (C1)
where {t i } are the TB-parameters, i and k denote the number of the electronic band and the wave vector, respectively, and w i,k are the weight coefficients. In Eq. (C1) E TB are the tight binding energies, which depend on the particular TB Hamiltonian and parameters used, and E GW are the starting-point values, obtained by DFT+GW or another ab-initio method. We found, however, that even if the parameter set reproduces the electronic bands with great accuracy, it does not necessarily give satisfying values of the g-factor. In order to overcome this problem, we included the calculation of the g-factor in our optimization procedure and varied the TB parameters to fit both the GW -energies and g-factor values. To prevent the order change of the energy bands during optimization we also included the error in the eigenvectors at the K-point. The weights were concentrated in the vicinity of the K and Γ points of the hexagonal Brillouin zone and had higher values for v and c bands for better reproduction of the most important electronic states. The fitting was carried out by an adaptive random search algorithm until a compromise between the quality of the band structure and g-factor value for MoS 2 was found.

Resulting parameterizations of TB Hamiltonians of Refs. [24] and [26]

  FIG. 1.Electronic spectra of MoS2 calculated using TB models H. Rostami et al.[26] (a) and S. Fang et al.[24] (b). Red lines are k•p quadratic dispersions at K-point calculated using effective masses for each band (see text for details). Note the different order of deep valence bands in two panels.

FIG. 2 .

 2 FIG. 2. Dispersion of the first fourLandau levels in the bottom conduction (c) and topmost valence (v) bands calculated using parametrization (a) of the k • p model (see Tab. III). Solid lines show the results of numeric calculations using decomposition Eq. (16), and dashed lines stand for a simple analytical formula εj(n, Bz) = ω

4 FIG. 3 .

 43 FIG. 3. Conduction and valence band g-factors, gc and gv, as functions of electron and hole energy, respectively. The energies are counted from the bottom of the conduction band and from the top of the valence band as shown in the inset. The solid and dashed lines show results of calculations for parameterizations (a) and (c) from Tab. III, respectively.

  are presented in Tabs. IV and V. The resulting energy dispersions in comparison with DFT+GW calculations are presented in Fig. 4. Based on these two sets of TB parameters we obtained k •p parameterizations listed in Tab. III, columns (d) and (e), with effective masses and g-factor listed in Tab. II, columns (d) and (e).

FIG. 4 .

 4 FIG. 4. DGT+GW and TB calculations for electronic energy dispersion in MoS2. (a)The fit of TB model of Ref.[24], (b) the fit of TB model of Ref.[26]. The zero energy was aligned to the top of the valence band v.

TABLE II .

 II The values of g-factors and effective masses (in units of m0) of the v and c bands calculated within the k • p model based on different TB parameterizations (listed in footnotes).

		MoS2 a	MoS2 b	MoS2 c	MoS2 d	MoS2 e	MoSe2 f	WS2 g	WSe2 h
	mv	-0.54	-0.72	-0.58	-0.40	-0.56	-0.82	-0.53	-0.57
	mc	0.54	0.86	0.90	0.37	0.37	1.02	0.68	0.76
	gv	8.73	5.57	6.18	11.90	5.59	5.12	6.08	5.64
	gc	7.82	5.41	6.83	10.15	1.77	5.12	6.13	5.79
	g X 0 = gc -gv	-0.91	-0.16	0.65	-1.75	-3.82	0	0.05	0.15

TABLE III .

 III Parameters of the k • p model as introduced in Eqs. (

TABLE IV .

 IV Parameters of the TB Hamiltonian (in units of eV) of Ref.[24] for MoS2 obtained after fitting of GW calculations.

			6	7 = 8		9	10 = 11	t	(1) 6,6	t	(1) 7,7	t	(1) 8,8
	-0.913 0.251 -1.538	-2.264	-0.922 0.437 -0.668
	t	(1) 9,9	t	(1) 10,10	t (1) 11,11	t (1) 6,8	t	(1) 9,11	t	(1) 6,7	t	(1) 7,8
	0.240	1.106 -0.003	0.046	-0.041 -0.762 -0.400
	t	(1) 9,10	t	(1) 10,11	t	(5) 9,6	t	(5) 11,6	t	(5) 10,7	t	(5) 9,8	t	(5) 11,8
	-0.168 -0.133 -0.975	0.016	1.829 0.914 -0.045
	t	(6) 9,6	t	(6) 11,6	t (6) 9,8	t	(6) 11,8
	0.935	0.945	0.796	0.449

TABLE V .

 V Parameters of the TB Hamiltonian (in units of eV) of Ref.[26] for MoS2 obtained after fitting of GW calculations.

	0	2	p	z	V pdσ	V pdπ	V ddσ
	-5.707 -5.784 -8.319 -12.171 4.791 -1.606 -1.221
	V ddπ	V ddσ	Vppσ	Vppπ
	0.526	0.359	0.905	-0.396

a TB model of Ref.[26] for DFT calculations b TB model of Ref.[24] for DFT calculations c TB model of Ref.[24] for DFT+GW calculations d DFT+GW and TB model of this work (TB model based on Ref.[26]) e DFT+GW and TB model of this work (TB model based on Ref.[24]) f TB model of Ref.[24] for DFT calculations g TB model of Ref.[24] for DFT calculations h TB model of Ref.[24] for DFT calculations

ACKNOWLEDGMENTS I.C.G. thanks the CNRS and the ANR MoS2ValleyControl project for financial support. He also acknowledges the CALMIP initiative for the generous allocation of computational times, through Project p0812, as well as the GENCI-CINES, GENCI-IDRIS, and GENCI-CCRT for Grant No. x2016096649. M.V.D. was partially supported by RFBR projects No. 14-02-00168 and No. 16-32-60175, the Russian Federation President Grant No. MK-7389.2016.2 and the Dynasty foundation. M.V.D. thanks M. M. Glazov for fruitful discussions.

 6) and (7) for different parametrizations, materials and tight-binding models are listed in Tab. III.

Appendix B: Computational details of DFT+GW method

The atomic structures and the quasi-particle band structures have been obtained from DFT calculations using the VASP package [START_REF] Kresse | [END_REF]38]. The Perdew-Burke-Ernzerhof (PBE) [39] functional was used as approximation of the exchange-correlation electronic term. The


J P Echeverry

I C Gerber

Theoretical investigations of anisotropic optical properties of distorted 1T ReS 2 and ReSe 2 mono-, bi-layers and in the bulk limit

Keywords: numbers: 71.35.-y, 31.15.A, 73.22.-f, 78.67-n

We present a theoretical study of electronic and optical properties of the layered ReX 2 compounds (X = S, Se) upon dimensional reduction. The effect on the band gap character due to interlayer coupling is studied by means of self-energy corrected GW method for optimized and experimental sets of structure's data. Induced changes on the optical properties as well as optical anisotropy are studied through optical spectra as obtained by solving the Bethe-Salpeter equation. At the G 0 W 0 level of theory, when decreasing the thickness of the ReS 2 sample from bulk to bilayer and to a freestanding monolayer the band gap reminds direct, despite a change of the band-gap nature, with values increasing from 1.6, 2.0 and 2.4 eV respectively. For ReSe 2 , the fundamental band gap changes from direct for bulk phase (1.38 eV) to indirect in bilayer (1.73 eV) and becomes again direct for a single layer (2.05 eV). We discuss these results in terms of the renormalization of the band structure. We produce the polarization angular dependent optical response to explore the optical anisotropy present in our results, as well as the fine structure of the lowest excitonic peaks present in absorption spectra.

I. INTRODUCTION

Isolation of a single and few layers from bulk's systems in many Transition Metal Dichalcogenides (TMDCs) has renewed the interest on their structural and optical properties broadly studied in bulk [1][2][3][4][5][6] . Dimensionality reduction's effects have motivated recent works on monolayers (MLs) and few layer's thickness systems, since it has been shown for the first time in MoS 2 that, upon dimensional confinement, the indirect band gap in bulk become direct in a ML 7,8 . Today, several others systems like W-based systems have shown similar behavior [9][10][11] . Such is the case of ReX 2 (X = S, Se) compounds, layered crystals belonging to TMDCs semiconductor's family which crystallize in a distorted 1T diamond-chain structure with triclinic symmetry unit cell usually denoted 1T . The successfull isolation of a single and few layers of ReX 2 , has allowed their structural and optical characterization [12][13][14][15] . However, there are many aspects on ReX 2 optoelectronic properties still remain to be clarified. See Ref. [16] for a recent review on experimental aspects of these TMDCs.

Based on their findings, some groups, have argued that in ReS 2 the interlayer coupling is negligible and the stacked layers in bulk behave as electronically and vibronically decoupled monolayers 12,17 . In contrast, recent Photo-Luminescence (PL) and Reflectance spectroscopy measures in bulk and few layer thickness samples of ReS 2 reveal a dependence of the optical properties on the thickness 15 . Additionally, it was usually claimed that ReS 2 shows a transition from indirect to direct optical band gap when the dimensionality is reduced from bulk to ML as in many other TMDCs 18 . This is in contradiction with the work of Tongay et al.

who have reported an unchanged indirect band gap from bulk to ML. A very recent work on ReS 2 bulk, using gated field-effect transistors, tends to confirm the indirect character of the gap 19 , when recent ARPES experiments tends to prove the contrary 20 .

From the theoretical side, some attempts to clarify these subjects using different methodologies can be found in the literature. Early results on bulk's band-structures based on Density Functional Theory (DFT) have shown ReX 2 compounds to be indirect band gap semiconductors 21 , while more recent theoretical studies tend to prove the contrary 12 . For ML and few-layers systems at standard DFT level, meaning that many-body screening effects are scarcely described, the fundamental gap is found to be direct but largely underestimated for ReS 2 12,[22][23][24] . On contrary, a ReSe 2 ML has an indirect band gap when explicit electronelectron screening is taken into account 14,22 . The only attempts of describing ReX 2 ML with the use of GW self-energy correction also give a direct gap of 2.38 eV at Γ-point, i.e at the center of the first Brillouin Zone (BZ), for ReS 2 as for ReSe 2 ML 22,25 . However for the latter case nothing is said about the difference between the indirect and the direct transition, since it is known that for different computational settings and/or a small strain effect, it is easy to influence the position of valence band maxima and conduction band minima 26 .

The presence of a strong optical anisotropy in ReX 2 compounds [27][28][29][30] , unlike Mo or W based hexagonal TMDCs, is due to their particular crystalline structure (see Fig. 1). It is believed that the structural distortion leads to a much weaker interlayer coupling, so the band renormalization is supposedly absent and so bulk ReX 2 should behave as decoupled monolayers 12,13 . However Aslan et al. 15 , demonstrated that the anisotropic optical properties of single-and few-layer ReS 2 crystals are different through polarization-resolved reflectance and photoluminescence spectroscopy. They have found that the near band-edge excitons in ultra-thin crystals absorb and emit light with preferred linear polarizations. Likewise, they have observed that the transition energies of the excitons can be tuned with layer thickness. They have thus established that ultra-thin ReS 2 has optical transitions with strengths and transition energies that depend on thickness and polarization of the optical radiation.

Here we propose to explore theoretically the optical properties of freestanding ML, bilayers (BL) of ReX 2 systems by including excitonic effects through GW and Bethe-Salpeter Equation (BSE) type calculations and compare them to bulk's spectra. We address the question about the character of the fundamental band gap upon dimensional reduction and we discuss it in terms of screening effects due to the presence of neighboring layers. Likewise, we discuss the dependence of the optical spectra on the thickness and on the polarization angle by rotating the dielectric tensor's components on the layer's plane. The present paper is organized as follows: in the section II we present the computational details for ML, BL and bulk calculations. In the section III the results concerning quasi-particle band-structure of ML, BL and bulk as well as the effects on the nature of the transition band-gaps when increasing the thickness are presented. Optical adsorption spectra for ML, BL and bulk are discussed in Section IV. Finally, we present main conclusions of this work.

II. COMPUTATIONAL DETAILS AND ATOMIC STRUCTURES

The atomic structures, the quasi-particle band structures and optical spectra are obtained from DFT calculations using the VASP package 31,32 . Perdew-Burke-Ernzerhof (PBE) functional 33 is used as approximation of the exchange-correlation electronic term. Core electrons are treated by the projector augmented wave scheme 34,35 . Fifteen electrons for Re atoms and six for S, Se ones are explicitly included in valence states. All atoms are allowed to relax with a force convergence criterion below 0.005 eV/ Å. During bulk and BL geometry optimization calculations, the functional of Langreth and Lundqvist 36,37 is used in the vdW-DF2 form 38,39 in order to describe accurately long-range correlations between layers. A grid of 6 × 6 × 1 k-points has been used, in conjunction with a vacuum height of around 20 Å, to take benefit of error's cancellation in the band gap estimates 40 , and to provide the imaginary part of the dielectric constant, a key quantity in absorption spectra in reasonable agreement with experiments as suggested in different works 41,42 for ML and BL calculations. A grid of 6 × 6 × 6 k-points has been used for bulk calculations. As it can be seen in Appendix A by using a 6 × 6 × 1 grid, it ensures a variation of only 50 meV compared to a 12 × 12 × 1 calculation for the position of the first excitonic peak. We use an energy cutoff of 400 eV and a gaussian smearing with a width of 0.05 eV is used for partial occupancies, when a tight electronic minimization tolerance of 10 -8 eV is set to determine with a good precision the corresponding derivative of the orbitals with respect to k needed in quasi-particle band structure calculations. Spin-orbit coupling was also included non-self-consistently to determine eigenvalues and wave functions as input for the full-frequency-dependent GW calculations 43 performed at the G 0 W 0 level but also at GW 0 level with 2 iterations of the G term, for testing purpose. The total number of states included in the GW procedure is set to 600, in conjunction with an energy cutoff of 100 eV for the response function, after a careful check of the direct band gap convergence, to be smaller than 0.1 eV in function of k-points sampling too. Band structures have been obtained after a Wannier interpolation procedure performed by the WANNIER90 program 44 . Only direct optical excitonic transitions are calculated by solving the BSE 45,46 , with the inclusion of 6 valence and 8 conduction bands in the diagonalization process for ML and bulk calculations. For the BL systems we use 12 valence and 16 conduction bands with for all optical spectra a broadening parameter of 25 meV is set. From these calculations, we report absorbance values defined as proposed in Ref. [47], by using the imaginary part of the complex dielectric function 2 (ω), with the following formula:

A(ω) = ω c 2 (ω)∆z (1)
, where ∆z is the vacuum distance between periodic images of ML and BL. This quantity is measurable and supposedly should not depend on the size of calculation cells in the ML/BL perpendicular direction. As pointed out by Bernardi et al 48 , Eq. ( 1) is a Taylor expansion for ∆z → 0 of the absorbance defined as A = 1 -e -α 2 ∆z for a single or bi-layer of a bulk material with a thickness ∆z, presenting an absorption coefficient

α 2 (ω) = ω 2 (ω) c n(ω)
, where for this case the refractive index is n = 1, since the considered ML/BL is surrounded by vacuum only.

Considering the lattice parameters and atomic positions in the bulk, as well as in ML or BL stacking, few issues need to be addressed. In the bulk's case, as it can be seen in Table I, the lattice parameters and the angles obtained upon geometry optimization using a van der Waals correction are in good agreement with experimental data for both materials. We recall that the 1T structure observed experimentally for bulk ReX 2 phases and used here for atomic configurations is triclinic, the corresponding space group being P1. As a consequence, MLs and bulk present a center of inversion, the symmetry point group is thus C i . For the ML systems, PBE optimization results (denoted D 1 ) are in line with previous theoretical studies 14,22,49 . Additionally we also report the D 2 structure as proposed in Ref. [50], which corresponds to a diamond-flip, or a displacement of Re atoms along the lattice vector b, with a small strain. Interestingly it implies a slight decrease of around 1% of the a parameter.

Since the work of Tongay et al. 12 , it is believed that contrarily to others TMDCs, ReX 2 systems possess very weak interlayer interactions. However when non-local correlation are

properly taken into account it appears that the binding energy between two layers in the simplest AA-stacking (a simple superposition of the two layers), is roughly the same as in the emblematic MoS 2 case, see Appendix B.

In order to investigate the BL stacking, we have kept fixed the lattice parameters from ML D 1 structures and displaced the upper layer in different direction, along a, along b and perpendicular to b, as proposed in previous studies 12,49 . We have taken into account the effects of van der Waals forces in the exchange-correlation term contrarily to previous works, and thus the resulting optimized stacking order is slightly different. For both ReS 2

and ReSe 2 systems the most stable stacking are yielded when the second layer is displaced

a ( Å) b ( Å) c ( Å) α ( • ) β ( • ) γ ( • ) ReS 2
Exp Bulk a 6.510 b Ref. [52] by 0.375 a from the AA-stacking, loosing the presence of a center of of inversion in the calculation's cell.. Mind the energy landscape is very flat with minima only separated by energy differences around few tenth of meV, as shown in Appendix C.

III. BAND STRUCTURES AND QUASIPARTICLE BAND-GAPS A. Monolayer

The G 0 W 0 band structures of distorted 1T -ReS BZ, E Γ g , for ReS 2 ML to be 2.38 eV, in remarkably good agreement with the previous work of Zhong et al. 22 Mind that when reducing lattice parameters to the experimental values, ReS 2 ML remains direct with a slightly smaller gap, 2.23 eV.

Γ K 1 K 2 M 2 Γ K 4 K 3 M 2 Γ -2 -1 0 1 2 3 4
Band Energy (eV)

D 1 D 2 E g Γ =2.38 eV Γ K 1 K 2 M 2 Γ K 4 K 3 M 2 Γ -2 -1 0 1 2 3 4
Band Energy (eV) eV. This result is in line with the reported value of 2.09 eV in Ref. [22], but much smaller than the 2.44 eV proposed in the work of Arora et al. 25 For sake of comparison we report in Figure 2, the behavior of the HVB and LCB around Γ point for the D 1 DFT-optimized lattice parameters and the experimental ones, taken as in the bulk phase. Interestingly when decreasing the lattice parameters, the ReSe 2 fundamental gap becomes indirect, since a maximum (denoted Λ in the HVB along the Γ-K 1 (-K 4 ) direction) pops up, the direct gap remaining 10 meV larger. This observation has been already reported in a previous study of ReX 2 but without any explanation families, for which drastic changes in LCB are observed. 53 This observation tends to confirm the existing electronic decoupling between ReS 2 layers evidenced in experimental works 12,17 .

E g Γ =2.05 eV M 2 Γ K 4 -1 -0.
However for the ReSe 2 BL, it is less obvious. On the one hand, similarly to the ML case, the VB maximum is slightly shifted from Γ-point along the Γ-K 1 (-K 4 ) direction, on the other hand the CB minimum is also slightly shifted away from Γ in Γ-M 2 direction. Thus E g is indirect with a value of 1.73 eV, when the direct gap in Γ is 1.77 eV. We have checked that the influence of the stacking order on the band-gap value is very limited, since for ReS 2 BL in AA-stacking configuration, E Γ g is 2.04 eV.

C. Bulk

In Figure 4 is plotted the self-energy corrected ReS 2 -bulk's band structure using high symmetry directions, as in Figure 1. Our G 0 W 0 calculation predicts bulk ReS 2 to be a semiconductor with a direct fundamental band gap of 1.60 eV in Z (0,0,1/2) coordinates, when the in-plane dispersion does not change that much with a HVB-LCB separation of 1.88 eV at Γ point. Those values are in the same range of the ones reported recently in

Γ K 1 K 2 M 2 Γ K 4 K 3 M 2 Γ -2 -1 0 1 2 3 4
Band Energy (eV) E g Γ =2.00 eV A significant band dispersion of the bands near Fermi level is observed in bulk's band structure compared with ML and BL cases for both ReS 2 and ReSe 2 . It is mainly due to interlayer interaction. Likewise, one can note a blue shift of the band-gap value upon dimensional reduction from bulk to ML as expected.

Γ K 1 K 2 M 2 Γ K 4 K 3 M 2 Γ -0

IV. ABSORBANCE

Absorbance of a single ReS 2 ML, in the D 1 structure, as a function of the frequency for the xx (Γ-M 2 direction) and yy (Γ-K 1 direction) components are given in figure 5. We recall that in ML cases for both ReX 2 systems, the VBM and CBM are doubly degenerated.

This implies, as proposed by Arora et al. 25 , that the first four peaks origins are the four inter-band transitions with different combinations of spin. As a consequence they are split

Γ Z 1 2 3 Z 4 5 6 X Γ Y -2 -1 0 1 2 3 4
Band Energy (eV)

E g Z =1.60 eV Γ Z 1 2 3 Z 4 5 6 X Γ Y -2 -1 0 1 2 3 4
Band Energy (eV) due to exchange interaction. This is clearly shown in the inset of the left part of figure 5. As soon as the active space of the BSE calculations includes at least two occupied and two virtual states, the exchange term split excitons in two different directions. Optical anisotropy is clearly evidenced in the spectrum. Main excitonic features appear in the energy range between 1.5 and 2.4 eV. When the most stable exciton, at 1.64 eV, has contributions along xx and yy directions, the second one, that lies 0.05 eV higher in energy is mainly oriented along Γ-M 2 (xx) direction. A similar behavior is observed for ReSe 2 ML with smaller intensities and a less pronounced anisotropic character in the first two excitonic peaks. We stress the fact that in other TMDCs this splitting does not involve anisotropy but exchange interaction and SOC 54 . Interestingly in Ref. [15], they report an optical gap value of 1.61 eV for ReS 2 ML by PL and reflectance spectroscopy experiments. For ReSe 2 ML the experimental optical gap is 1.50 eV 25 , when it is located at 1.43 eV in our calculations.

The overall agreement between experimental results and our estimates is good, considering the precision of GW +BSE calculations scheme.

The binding energy (BE) of the lowest exciton in ReS 2 ML is 0.74 eV, much larger than other M or W-based TMDCs. For ReSe 2 ML the BE is 0.57 eV when the exchange splitting is only 0.04 eV. The binding energy values are in line with previous reports, they only differ due to different computational details 22,25 . Additionally, we have noticed fine structure of those peaks with the presence of a dark state 7 meV lower in energy than the bright one, and the presence of a grey exciton, possibly coupled to light with z-polarization as recently observed in WSe 2 ML 55 , with a small oscillator strength 4 meV above the bright state.

The next higher energy peak is seen at 1.83 eV with a stronger oscillator strength for yy component than for the xx one.

When a second layer is added to a ReS 2 ML, with two different stacking orders, there is a red-shift of the lowest excitonic peak of 0.04 eV. One can still observe the anisotropic character by decomposing the absorbance function in xx and yy directions. The orientation of the second layer with respect to the first one does not change drastically the spectrum.

The only noticeable difference can be seen when the stacking is highly symmetric (AAstacking), see Fig. 6, the peak at 1.76 eV is much more intense, due to the perfect alignment between the atomic positions. The corresponding BE is 0.40 eV for the optimized stacking order and slightly larger (0.44 eV) for AA-stacking. The screening effect on the first ReS 2 ML due to the presence of a second layer has not significant effect on the energy of the electronic transitions while the band-gap renormalization induces important changes on the BE. Despite the fact that a ReSe 2 BL has an indirect band-gap in our calculations, we have performed BSE calculations, in order to investigate direct transitions. As for ReS 2 BL, those direct transition energies are not altered by the addition of a second layer in ReSe 2 .

The first two direct transition are still located at 1.41 eV and 1.44 eV for the BL, a similar trend is obtained for ReS 2 . Thus the electronic screening due to the second layer in ReX 2 compounds only implies a reduction of the fundamental gap, possibly altering its direct to indirect character for ReSe 2 only and reducing significantly the BE. Figure 7 shows the absorption coefficient α 2 as a function of the frequency for both bulk of ReS 2 and ReSe 2 systems. Both spectra present a rich excitonic structure in the energy range from 1.3 to 1.6 eV. By comparing in-plane contributions to the dielectric tensor, we clearly evidence the optical anisotropy still present in the bulk phase of Re-based TMDCs. In the mentioned energy range we observe that the excitonic response is slightly stronger for the yy component than for the xx one. For ReS 2 , the lowest-lying feature appears now at 1.40 eV, in good agreement with the experimental determination of the optical gap of Ref. [15], i.e. 1.47 eV. In the same article, it is reported the presence of a feature at 1.4 eV in PL but not in reflectance experiments. Based on this, they have concluded that bulk ReS 2 to be a semiconductor with indirect excitonic band gap. Marzik et al. 56 have concluded from their photo-electrochemical measurements an indirect band gap of 1.4 eV for ReS 2 single crystals, when from low temperature absorption measurements Friemelt et al. 28 have determined a value of 1.56 eV for the direct optical band gap of p-type ReS 2 single crystals and 1.53 eV for n-type ones. Furthermore, they have also found a sharp peak in low temperature absorption spectrum at 1.57 eV which was assigned to an excitonic peak. These values are in reasonable agreement with our findings but the assignment contrast with our conclusion, since our band structure results confirm a direct fundamental gap for ReS 2 in bulk, in Z k-point. Concerning the optical band gap, we safely assign the feature present at 1.40 eV to be a direct excitonic band gap since in any BSE calculation only vertical transitions are taken into account. Two others intense peaks at 1.51 and 1.56 eV are theoretically determined, in very good agreement with values extracted from reflectance contrast measurements, 1.51 and 1.58 eV. The theoretical blue shift when decreasing the thickness of the sample is thus 0.24 eV, when its experimental value is 0.14 eV 15 . The next two main features are shifted in energy by 0.18 and 0.27 eV respectively, when experimentally those values are 0.16 and 0.30 eV. The BE of exciton in ReS 2 bulk is thus estimated to be 0.20 eV, reduced by 75% when compared to the ML limit. For ReSe 2 in the bulk phase, the lowest energy peak in the spectra is composed of three bright excitons, the dominant strength oscillator exciton is located at 1. In order to study the optical anisotropy of ReX From bulk to Ml thickness one can see an increase of the lowest excitonic peak's intensities as well as blue shift of their transition energies. For each excitonic feature in the spectra one can identify the easy polarization axis, a preferential direction where transition probability become maximum. The angle where the intensity reach a minimum corresponds to weakest oscillator strength or a polarization angle where transition are forbidden. We note that oscillator strength of the excitons in the lowest features of ML and BL has different easy polarization direction. As a consequence, in the angular evolution of the features, those excitons seems out of phase. The variation of the intensity of the lowest energy features from xx component to yy component is an evidence of the optical anisotropy of the ReX 2 compounds. Oscillator strength of the excitonic features become stronger over dimensional reduction.

V. CONCLUSION

In the present work, we propose a theoretical study of electronic and optical properties of the layered ReX 2 compounds for ML, BL and bulk phase, by means of GW +BSE scheme.

For ReS 2 the band gap reminds direct, despite a change of the band-gap nature, with values increasing from 1.6, 2.0 and 2.4 eV from bulk to ML respectively. For ReSe

2 FIG. 1

 21 FIG. 1. a) Top view and b) side view of distorted 1T monolayer ReX 2 , Re atoms are in grey, X (S, Se) atoms are in yellow. The cell used during calculations is depicted by black line in a). TheRe chain is along the b direction. c) First Brillouin zone of the ReX 2 monolayer and bulk. For sake of comparison, the notation is the same as the one used in Ref.[25] .

FIG. 2 .

 2 FIG. 2. On the left, the G 0 W 0 band structures of ReS 2 monolayer in D 1 (in red) and D 2 (in black) phases. On the right, the G 0 W 0 band structures of ReSe 2 monolayer in D 1 phase. The inset presents a zoom around the Γ-point for two set of lattice parameters: the DFT-optimized lattice parameters of D 1 phase and the experimental ones. The energy level of the valence band maximum (in Γ) is set to zero.

FIG. 3 .

 3 FIG.3. On the left, the G 0 W 0 band structures of ReS 2 bilayer in the most stable configuration, i.e with a displacement of around 2.4 Å along a of the second layer. On the right, The G 0 W 0 band structures of ReSe 2 bilayer in the most stable configuration, i.e with a displacement of around 2.4 Å along a of the second layer. The energy level of the valence band maximum (in Γ) is set to zero.

FIG. 5 .

 5 FIG. 5. Absorbance of monolayer ReS 2 (left) and ReSe 2 (right) in the Γ-M 2 direction (xx) and Γ-K 1 direction (yy).Vertical dashed line represents the fundamental band gap energy. The inset of the left figure presents the oscillator strengths with respect to the number of occupied and unoccupied states included in the BSE calculations.

FIG. 6 .

 6 FIG. 6. Absorbance of ReS 2 BL in the optimized configuration (left) and in AA-stacking configuration (right) in the Γ-M 2 direction (xx) and Γ-K 1 direction (yy).Vertical dashed line represents the fundamental band gap energy.

FIG. 7 .

 7 FIG. 7. Optical absorption coeeficient of bulk ReS 2 (left) and ReSe 2 (right) in the Γ-M 2 (xx), in the Γ-K 1 (yy) and Γ-Z (zz) directions.Vertical dashed line represents the fundamental band gap energy.

2

 2 compounds and, in particular, to understand the angular dependence of the low-energy features in the spectra, we have rotated the in plane dielectric tensor component on the layer around the c axis (perpendicular to a and b). In figure8, the imaginary part of the dielectric tensor of ReS 2 systems as a function of the polarization angle and the energy of the incoming light is presented. ReSe 2 results are given in Appendix D. The blue line corresponds to yy dielectric tensor component. We note that the optical anisotropy near the excitonic resonances implies that the absorption of these layers varies strongly as a function of the light polarization. For the low lying features below the fundamental gap values, the maximum is always obtained for light polarization aligned with the Re-Re chain, so with an angle of around 120 • , for ML, BL and bulk too.

FIG. 8 .

 8 FIG. 8. Imaginary part of the dielectric tensor of ReS 2 systems, with respect to light polarization angle.

FIG. 10 .

 10 FIG. 10. Comparison between the binding energy of MoS2 and ReS 2 BL for different schemes that include van der Waals forces in the hypothetical AA-stacking.

FIG. 11 .FIG. 12 .

 1112 FIG. 11. Energy differences from AA-stacking for various distances in three directions of displacement.

TABLE I .

 I Calculated and experimental lattice parameters and angles for bulk and monolayer ReX 2 systems.

  values, without changing the positions of the VB maxima and CB minima. The direct gap in Γ values are 2.60 eV and 2.23 eV for ReS 2 and ReSe 2 respectively at this level of theory. Figure 3 provides band structures for both ReS 2 and ReSe 2 BLs. In both cases, no spinsplitting is observed. The interlayer coupling impacts the band structure by renormalizing the band gap by roughly 16%. For ReS 2 BL the positions of the CB minimum and the VB maximum remain the same as in the ML case. As a consequence the ReS 2 BL band gap is still direct in Γ and is decreased to 2.00 eV. This is clearly different for MoX 2 or WX 2

	B. Bilayer

22

. The importance of lattice parameters choices have been pointing out in a previous work on MoX 2 ML systems

26

, especially when Se atoms are involved. As usual, the use of partially self-consistent G 2 W 0 calculations, provides larger band-gaps

  FIG.4. The G 0 W 0 band structures of ReS 2 bulk, on the left and for ReSe 2 bulk on the right.The energy level of the valence band maximum (in Z) is set to zero. High symmetry directions of the 1BZ in bulk crystal are depicted in Figure1.

								E g Z =1.38 eV
		20
	Absorbance [%]	5 10 15	0 10 3	1.6	1.8	2.0	1-2 2-2 4-4 6-8
		1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 Energy [eV] 0	1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 Energy [eV]

  16 eV and it is preceded by a couple of weaker strength oscillator excitons located few meV lower. The BE of bulk ReSe 2 is 0.22 eV considering the direct fundamental band gap also located in Z. Previous experimental works 52,56,57 using different techniques have obtained values of 1.15, 1.17, and 1.22 eV for single crystal optical band gap respectively. More recently PL the first two peaks when reducing layer thickness from bulk to ML is 0.27 and 0.26 eV respectively. This is clearly much larger than in any others TMDCs, but probably overestimated in our calculations, when compared to experiments25 . The corresponding excitonic BE is 0.21 eV.

experiments

58

report a value of 1.26 eV for the optical band gap of bulk ReSe 2 , when a larger value (1.37 eV) is also proposed

25

. Globally in this case, the dispersion of the optical band gap is still relatively large. The next higher energy features are centered around 1.21 and 1.29 eV, the latter being broader for the xx component. As a consequence, the blue shift observed on

Ref.[25]. This is clearly different from others Mo and W based MX 2 families where the HVB maxima always remain at Γ and at the same time the LCB minima still lie in the (k x , k y ) plane. For those systems relatively strong interactions between the chalcogen atoms of the neighboring layers result in the direct-to-indirect band gap crossover with increasing the number of layers, when for Re-systems the nature of band gap itself is changed, favoring

ACKNOWLEDGMENTS

J. P. E. thanks the ANR MoS2ValleyControl project for financial support. The authors also acknowledge the CALMIP initiative for the generous allocation of computational times, through Project p0812, as well as the GENCI-CINES, GENCI-IDRIS, and GENCI-CCRT for Grant No. x2016096649 and A0020906649.

Appendix A: Convergence issue of the imaginary part of the dielectric constant with respect to k-points sampling 1. 4 1 To investigate convergence issues with respect to the k-points sampling, we have neglected SOC term but kept all the others calculation parameters as they are given in the computational details section. If one checks the first excitonic peak in the imaginary part of the dielectric constant for 3 different grids, the difference is around 0.05 eV. The position of higher energy excitons also follow the same trend, but with 12 × 12 × 1 grid new features appear, probably due to a better description of excited states of the excitons. Mind that by adding SOC term, symmetry is removed in the calculation, thus increasing significantly the computational cost, that makes a 12x12x1 calculation intractable in a reasonable time. In figure 10, the binding energy of MoS 2 and ReS 2 bilayers in AA-stacking are given. PBE results are compared to the curves obtained by either the Tkatchenko-Scheffler method 59 , the DFT-D3 method 60 or the vdW-DF2 scheme 38,39 too. Globally similar range of energy The stability of the stacking with respect to translational displacement, along a, along as well as perpendicular to b have been investigated using vdW-DF2 functional. In figure 11, the primitive cell and the lattice vectors directions are recalled, as well as the energy differences from the hypothetical AA-stacking, i.e the simplest superposition of the two layers. The atomic positions were allowed to relax too.


Prof. Dr Ralf Seppelt
email: ralf.seppelt@ufz.de

Peter H Verburg

Albert Norström

Wolfgang Cramer

Tomáš Václavík

Focus on cross-scale feedbacks in global sustainable land management

Keywords: place-based research, transferability, land use, tele-coupling, scenario analysis

Introduction

Human land-use activities have transformed most of the Earth's land surface (Foley et al 30 2005[START_REF] Ellis | Anthropogenic transformation of the terrestrial biosphere[END_REF][START_REF] Gauthier | Boreal forest health and global change[END_REF]. While land-use activities differ in many ways across the world, their combined impact is becoming a force of global importance.

Consequently, sustainable land management (SLM) has been identified as a key lever for achieving global sustainability. For example, six out of 17 sustainable development goals (SDGs), adopted in the United Nations 2030 Agenda for Sustainable Development, relate 35 directly to land management: (i) land management is key for providing goods and services for human kind relating to "zero hunger" (SDG 2); (ii) land management is responsible for 20-40% of GHG emissions relating to "climate action" (SDG 13), and (iii) land management is a major driver of biodiversity loss relating to "life of land" (SDG 15).

Indirectly land management affects "affordable clean energy" (SDG 7) and "sustainable 40 cities and communities" (SDG 11).

A growing human population, associated with increasing consumption rates and demands on commodities, requires a true paradigm shift regarding the management of the land for long-term sustainability. At the same time we are witnessing a progressive scarcity of available productive land, and the production peak of many renewable 45 resources has already been passed [START_REF] Lambin | Inaugural Article: Global land use change, economic globalization, and the looming land scarcity[END_REF]Meyfroidt 2011, Seppelt et al 2014).

These demands and the limits to supply underlie the many linkages between the different social, economic and ecological goals and targets that are being charted out by Agenda 2030 [START_REF] Geijzendorffer | Ecosystem services in global sustainability policies[END_REF]. Pathways to some goals are synergistic, while others present trade-offs for their mutual achievement. For example, an increase of 50 provisioning goods and services from ecosystems, such as food and fibre, could be achieved through further intensification of land use [START_REF] Mauser | Global biomass production potentials exceed expected future demand without the need for cropland expansion[END_REF], which might lead to an increasing loss of biodiversity [START_REF] Newbold | Global effects of 345 land use on local terrestrial biodiversity[END_REF], Gerstner et al 2014, Stein et al 2014) but also results in higher greenhouse gas (GHG) emission through higher energy use and fertilizer application. At the same time, biodiversity has to be maintained 55 for many societal objectives including its potential to support ecosystem functions such as pollination [START_REF] Cardinale | Biodiversity loss and its impact on humanity[END_REF], Seppelt et al 2016). A significant reduction of GHG emissions through large-scale deployment of new biofuels is in conflict with the production of food or conservation of natural habitats for biodiversity or carbon objectives strongly depends on the local land system and the environmental, socioeconomic and cultural context in which this land system is operating. Thus, achieving one SDG might compromise others [START_REF] Pradhan | A Systematic Study of 350 Sustainable Development Goal (SDG) Interactions Earth's Future[END_REF] and trade-offs on various scales need to be expected, which can be moderated by appropriate land management.

Whilst place-based research provides essential knowledge on the biophysical and socio-65 economic boundaries of land use, its findings are naturally contingent upon the specific geographical context and rarely account for off-site effects. The conversion of a conventional agricultural system at one location to organic farming may have positive impacts on local sustainability, but it may, due to lower production, displace some impacts to other locations that need to compensate for the loss of production. On the 70 other hand, many land use drivers such as climate change, population growth or consumption patterns are well captured at the global scale, but there are significant uncertainties about how they interact with local conditions. Both regional and global studies on food production rarely account for these tight links and interactions between socio-economic and biophysical processes. These uncertainties and incongruences in 75 spatial scales prevent effective integration, synthesis and transferability of findings from research to sustainable land management.

Against this background, we here review and synthesize the contributions of the focus collection on "Cross-scale Feedbacks in Global Sustainable Land Management", which collates papers that investigate the links between global change processes and local 80 realities through, e.g. integration of local and global drivers impacting economic and biophysical processes or assessing the transferability or up-scaling of findings from place-based research.

Synthesis: Emerging Topics in Sustainable Land Management and Land

Systems Research

85

The articles in this focus issue illustrate new approaches to investigate global and regional land systems, and identify key research frontiers important for sustainable land management across scales to achieve the SDGs. Three major clusters of research frontiers have been identified: Improving our understanding of these different cross-scale dynamics, in diverse land-95 use systems is critical. [START_REF] Dorninger | Assessing sustainable biophysical human-nature connectedness at regional scales[END_REF] conceptualize the "human-nature connectedness" as a new methodological framework that can be applied in any region of the world to assess how closely connected people are to their regional ecosystems. The authors identify two key mechanisms that disconnect humans from nature on a regional scale: (1) the flow of external non-renewable inputs into the land-use system and ( 2) 100 teleconnections with distant systems. While these mechanisms allow for greater regional resource use, they pose challenges for sustainability through waste generation, depletion of non-renewable resources and environmental burdens shifted to distant regions.

The topic of environmental burdens is elaborated by [START_REF] Pascual | Off-stage ecosystem service burdens: A blind spot for global sustainability[END_REF], who argue for 105 a better recognition of the distant, diffuse and delayed impacts that land management often has on biodiversity and ecosystem services. They define these impacts as "offstage ecosystem service burdens" and identify four typical pathways based on biodiversity conservation policies, and the management of provisioning, regulating and cultural services. The authors advocate for their incorporation in land management 110 decisions and ecosystem service assessments such as those conducted by the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES).

Finally, [START_REF] Sietz | Nested archetypes of vulnerability in african drylands: Where lies potential 375 for sustainable agricultural intensification[END_REF] apply the emerging approach of "archetype analysis" and assess vulnerability in African drylands to environmental change. A cluster analysis 115 reveals archetypical patterns of how vulnerable farming systems are to land degradation and support understanding the heterogeneity of vulnerability determinants (e.g. water availability, agro-ecological potential or population pressure) across sub-Saharan Africa. Importantly, their spatially explicit framework offers the opportunity of evaluating a specific region's potentials and challenges in its wider context across 120 nested scales.

Synthesis of Place-Based Research Results

While new conceptual and analytical frameworks such as the above provide guidance in designing specific analysis and provide suggestions for similarities between case A c c e p t e d M a n u s c r i p t unresolved and defines the second emerging topic. Three studies focus on synthesis of data and local case studies on sustainable land management, highlighting the need for generalization and transferability of findings. [START_REF] Hermans-Neumann | Why do forest products become less available? A pan-tropical comparison of drivers of forest-resource degradation[END_REF] analyse the drivers of changes in tropical forest products using a standardized, pan-tropical dataset of more than 200 villages with forest access. Their analysis show that forest 130 resources (e.g. timber, fuel wood and food) declined over the last five years, though with marked differences across continents. The strongest degradation of forest resources occurred in places with both growing resource use and immigration.

Similarly, [START_REF] Carter | 280 Large scale land acquisitions and REDD+: a synthesis of conflicts and opportunities[END_REF] synthesize comprehensive data on large scale land acquisition (LSLA) to reveal that land available for agriculture, accessibility and political 135 stability are the main factors that explain whether a country will be targeted for LSLA.

The synthesis of such comprehensive datasets allows for globally comparative analyses that go beyond case studies in terms of generalizable conclusions and transferability of findings.

The issue of transferability is specifically addressed by [START_REF] Václavík | Investigating potential transferability of placebased research in land system science[END_REF] who build 140 on a previously developed concept of land system archetypes [START_REF] Václavík | Mapping global land system 380 archetypes[END_REF][START_REF] Levers | Archetypical patterns and trajectories of land systems in Europe Reg[END_REF] to investigate potential transferability of regional case studies that focus on land management and ecosystem services across four continents. The proposed method is offered as a blueprint for large research frameworks that need to assess the relevance and representativeness of place-based research for other geographical areas 145 and to identify possible gaps in research efforts.

Future Perspectives and New Integrated Scenarios

Finally, the third emerging topic deals with potential future developments of land management, based on scenarios which capture the diversity of land systems. Scenarios of land management need to consider that land is a limited resource which can be used 150 to produce food and fibre or to maintain non-provisioning ecosystem services but that trade-offs exist in the decisions how to manage land [START_REF] Titeux | Global scenarios for biodiversity need to better integrate climate and land use change[END_REF]. Modelling provides the base for future forest management optimization at the regional and national scale.

Outlook: Guiding questions of the emerging topics

Obviously, a focus collection of publications as this cannot fully cover the complex topic of sustainable land management in a comprehensive or exhaustive manner. Global land 175 management is characterized by a diverse set of key challenges. These range from sustainable resource appropriation, the preservation of biodiversity and ecosystem functioning, human well-being, equity and so forth, as comprehensively covered by the SDGs, or as illustrated by Figure 1. Based on the emerging topics identified by and discussed in the contributions of this focus collection we can collate key questions, 180 which serve for stimulation as well as for guidance of future research directions, c.f. Box 1.

A c c e p t e d M a n u s c r i p t

190

There is a predominant paradigm that an improved understanding of the system of interest, here land systems, is key to improve decision making. However, despite a significant increase in our understanding of land system dynamics over the past two decades, the uptake and integration of scientific knowledge into decision-making processes remains limited [START_REF] Kirchhoff | Actionable Knowledge for Environmental Decision Making: Broadening the Usability of Climate Science[END_REF]. This is because decision-making only 195 partly relies on well-established scientific knowledge. Of equal importance are the underlying value systems of the involved decision-makers, beneficiary and stakeholders as well as the governance system and power structures in which decision can be taken, namely rules, values and knowledge (rvk), c.f. Goddard et al (2016). Consequently, a simplistic focus on generating more understanding of land system dynamics, will likely 200 in itself be insufficient to foster sustainable land management.

A c c e p t e d M

a n u s c r i p t production (Mauser et al 2013), could increase the relevance and usability of land-use science for society and decision-making. With respect to global land governance, new emerging processes such as large-scale land acquisition or spill-off and offsite-effects 205 [START_REF] Carter | 280 Large scale land acquisitions and REDD+: a synthesis of conflicts and opportunities[END_REF][START_REF] Pascual | Off-stage ecosystem service burdens: A blind spot for global sustainability[END_REF][START_REF] Seppelt | A quantitative review of ecosystem service studies: approaches, shortcomings and the road ahead[END_REF] pose challenges to land management which is mostly implemented through law, rule or incentives at the local to regional scale. Surprisingly, large-scale land acquisitions are not an issue in global scale agro-economic models [START_REF] Debonne | Representing large-scale land acquisitions in land use change scenarios for the Lao PDR Reg[END_REF]. A proper representation of changes in farming structure, including their underlying social, economic and political drivers, is 210 important to be able to analyse the environmental, economic and social impacts of such changes and the ways in which these new modes of land governance impact on the relations between global and local processes. The limitation of the available land surface and the limitation of its goods and services produced simply suggests that novel ideas to govern land as global commons are required (Seppelt et al 2014[START_REF] Creutzig | Govern land as a global commons[END_REF]. 215

This focus collection also contrasts two different conceptual approaches to the synthesis of place-based research results: (a) global-scale analysis and modelling, that builds on the basic assumption to fully capture global processes related to land use and (b) linkage of a variety of locations studied as different case studies. While global trade models are limited with respect to spatial scale but also with respect to the commodities captured, a 220 similar limitation holds for the synthesis of place-based results. For the latter, concepts like tele-coupling, off-site effects or spill-over emerge quite logically [START_REF] Pascual | Off-stage ecosystem service burdens: A blind spot for global sustainability[END_REF][START_REF] Liu | Systems integration for global sustainability Science 347 1258832-1258832 sustainability[END_REF][START_REF] Seppelt | A quantitative review of ecosystem service studies: approaches, shortcomings and the road ahead[END_REF]. Although tele-coupling is well conceptualized, operationalization in research methods is still challenging. The literature is full of local case study evidence of potential impacts of emerging value chains and the role of 225 market-based commodities and tele-coupled land management [START_REF] Lenzen | International trade drives biodiversity threats in developing nations[END_REF].

However, these local insights are poorly coupled to larger scale assessments and life cycle analysis (LCA) where impacts are only considered 'on average' ignoring the importance of local land systems as determinants of the impacts of these global

relations. 230

For developing future perspectives on sustainable land management based on scenario approaches a better integration of feedbacks is needed, c.f. [START_REF] Delzeit | Global economic-biophysical assessment of midterm scenarios for agricultural markets-biofuel policies, dietary patterns, cropland expansion, and productivity growth[END_REF]. Gaps still relate to understating land use intensity, landscape homogenization and the feedback between landscapes, agricultural production and biodiversity of managed landscapes (Seppelt et al 2016, Verburg et al 2016), specifically as humans shape emerging or novel 235 A c c e p t e d M a n u s c r i p t ecosystems. Two understudied feedbacks pose major challenges for future global land systems research. First, the mutual dependence between biodiversity and agricultural production is understudied in global studies and models: Biodiversity is negatively affected through land-use intensification, which is mostly applied to boost yields.

Maintaining yields on a high level, however, requires various facets of biodiversity for 240 support of important ecosystem functions such as nutrient cycling, biocontrol or pollination (Seppelt et al 2016). Research for embedding these feedbacks quantitatively in global scale models is advancing and could become crucial for global assessments in the near future (Rosa et al 2016).

Second, the feedback between commodity production and consumption are also 245 understudied. Usually demand trajectories are predefined, such as by predefined scenarions, e.g. [START_REF] Delzeit | Global economic-biophysical assessment of midterm scenarios for agricultural markets-biofuel policies, dietary patterns, cropland expansion, and productivity growth[END_REF]. Jevons paradox-like, rebound phenomena, which denote the increasing demand for a resource after establishing a more resource efficient production method [START_REF] Alcott | Jevons' paradox[END_REF], are mostly neglected in today's scenarios analysis. This might hamper understanding of rebound-effects and probably leads to 250 overly positive estimates on certain scenarios.

Given the limited land resources available and multiple competing claims on these resources sustainable land management should also include sustainable consumption [START_REF] Scherer | Mapping and linking supply-and demand-side measures in 360 climate-smart agriculture[END_REF]. While sustainable consumption is also one of the SDGs it is hardly related to land management. Recent work of [START_REF] Alexander | Human appropriation of land for food: The role of diet Glob[END_REF] shows the 255 strong impacts of consumption and value-chain losses on agricultural production, indicating the potential reduction on land resources that can be achieved through improved consumption and value-chain management.

This focus collection provides several contributions in the fields of land system science to the development of concepts, models and tools for sustainable land management. To 260 advance beyond the current state of the art, future research directions need to address a diversity of topical challenges such as poverty reduction, large scale land acquisition, global feedbacks of agricultural production and biodiversity. While research questions can be developed easily, we acknowledge that further research needs shall not hamper action with respect to lowering pressure on the environment by all possible means. 265

Research needs are no excuse for inaction [START_REF] Voinov | Values in socioenvironmental modelling: Persuasion for action or excuse for inaction[END_REF].

  (1) new frameworks to understand cross-scale dynamics of land-use systems, (2) synthesis of place-based research, and (3) addressing future 90 perspectives of land use by development of consistent scenarios.

  frameworks to understand cross-scale dynamics of land-use systems Land management dynamics are seldom just local or place-specific anymore, but are influenced by multiple global drivers with complex connections to other places.

  approaches that integrate the interplay of biophysical and socioeconomic factors in scenarios of global or regional change are promising tools to study future land-use impacts and trade-offs.[START_REF] Delzeit | Global economic-biophysical assessment of midterm scenarios for agricultural markets-biofuel policies, dietary patterns, cropland expansion, and productivity growth[END_REF] provide a set of scenarios of global drivers 155 until 2030 that can be used consistently in a range of regional and local case studies of land use. The impacts of biofuel policies, dietary patterns, cropland expansion and A c c e p t e d M a n u s c r i p t productivity changes on agricultural markets are investigated in a modelling framework that couples an economic model with a crop growth model. One of these global scenarios together with regionally-tailored land-use and climate 160 change scenarios is applied by Langerwisch et al (2018) who quantify the combined effects of land-use and climate change on four ecosystem services in rice-production regions in Southeast Asia. Here, the vegetation and hydrology model LPJmL shows clear trade-offs in the future provision of ecosystem service, but also the potential of land management to partially offset the negative impacts of climate change on rice 165 production, carbon storage and sequestration. Following a similar framework, Gutsch et al (2018) quantify the effects of alternative land management scenarios and climate impacts on forest variables indicating ecosystem services related to timber, habitat, water and carbon. Again, the combination of modelling tools applied under scenarios of future change allows to better balance the trade-offs between ecosystem services and 170

Figure 1 :

 1 Figure 1: The diverse facets of global sustainable management of land systems (Artist: M. Volk). The main task of balancing different trade-offs, such as between various SDGs, by the artistic guy in the center of the

-

  Outstanding questions of the emerging question in sustainable land management Enhancing understanding of the diversity of land systems -Are the correct drivers addressed for investigating solutions on sustainable land management, considering the knowledge, values and rules define the decision context? -How to capture countries activities and characteristics properly to account for emerging issues such as large-scale land acquisition, or long-distance externalization of effects within global agro-economic models? -What are the options to govern land as global commons? Synthesis of place-based research results -Which are the next steps to enable global agro-economic models to address a larger set of commodities, different land-holding systems, capture nutrient cycling and provide sufficient information on food security question on a finger spatial resolution? -How can tele-coupling concept be operationalized in research to better underpin and embed life cycle analysis in global relations? -Which data gaps should be closed to better account for local variations in the socio-economic context of sustainable land management? Future perspectives and new scenarios -How to implement the mutual feedback of biodiversity and agricultural production in today's global model system estimation global agricultural yields and estimate optimum intensification levels? To what extent do concepts like sustainable intensification that claim to have synergies between SDGs really have potential, what are the trade-offs hidden in these systems and in what local context are such concepts applicable? -How can integrated scenarios capture the links between production and consumption, rebound effects and Jevons paradox? 270 A c c e p t e d M a n u s c r i p t A c c e p t e d M a n u s c r i p t

Acknowledgements

This study was funded by the German Federal Ministry of Education and Research (BMBF) as part of the GLUES project (01LL0901A). We are grateful to M. Volk and his creativeness providing us Figure 1 with a cartoon-like illustration of the topics in this focus collection. 275