First-principles study of magnetism and electric field effects in 2D systems

Gating a junction with electric fields is a very common experimental method to control functionality and properties of a system. One example with the most significant societal impact is perhaps the field-effect transistor (FET) that, in conjunction with the development of metal-oxide-semiconductors (MOS), because of its high scalability led to a digital revolution in the 1950s. First proposed in the mid-late 1970s,¹ tunneling field effect transistors (TFETs) attracted much attention in mid-2000 because of nanostructures involving carbon nanotubes^2,3 and the intense interest continues. In 2010, low-voltage tunneling FETs based on inter-band tunneling⁴ were fully analyzed, followed by discovery of the single-layer MoS₂ TFET.⁵ Shortly after, junctions with vertical geometry using layered hexagonal boron nitride (h-BN) and MoS₂ 2D materials were built.^6,7 The performance characteristics, for example, the speed, of TFETs are not limited to the Maxwell–Boltzmann tail as in the conventional MOSFET. To date, two-dimensional systems are regarded as promising materials for future low dissipation electronics that are not limited to transistor applications.^8–15 Gating a system provides an easy method to applying an electric field and allows us to investigate the magnetic response of systems and thus magnetoelectric coupling, as well as a method to study charge-doping effects.

Theoretically, modeling of gate field effects was motivated primarily by studies of MOSFETs. In the early days, the circuit model or classical electromagnetic theory (for insulated-gate FETs) was used to model field-effect transistors (FETs).^16,17 Later, studies of the system at the electron level with inclusion of its quantum nature have been carried out within various approximations. The key is to solve the Poisson equation with appropriate boundary conditions.^18,19 In this paper, our focus is on first-principles modeling and simulations of gate effects. Because of the limitations of computational power and algorithms, solutions to the Poisson equation subject to prescribed boundary conditions as commonly used in classical E&M have not been done until very recently, and to our best knowledge, there are, in practice, only a few approaches that faithfully realize the E&M principles. One approach has been developed by our group^20,21 by employing the effective screening medium (ESM) technique as the Poisson solver.²² In the last few years, we have applied this method to study both single- and dual-gate configurations, and a number of two-dimensional junctions and interfaces between bulk systems have been investigated.^{20,21,23–28} In Sec. II, we will highlight some important results from these studies. With this approach, boundary conditions in the z-direction are imposed according to the physical problem, which according to the uniqueness theorem guarantees the correct solution. The advantage of this approach is that it is combined with the non-equilibrium Green's function technique to study electron transport at finite bias, implemented in the TranSiesta package²⁹ in addition to the Quantum ESPRESSO package.³⁰ Besides TranSiesta, QuantumATK³¹ can treat field effects by enabling Dirichlet (potential is held constant) and Neumann (electric field is held constant) boundary conditions, similar to our approach. For a generalized Poisson equation, Bani-Hashemian et al.,³² developed an algorithm to treat Dirichlet-type boundary conditions for first-principles device simulations. In addition to the boundary-condition-driven approach, several other methods should be mentioned because of their impact in current research. One is by Sohier et al.,³³ in which the authors truncate the Coulomb interaction in the direction perpendicular to the slab so that the charging of the slab can be simulated via field effects. An important development in this method is the treatment of flexural phonons in the presence of a field, which breaks the mirror symmetry, thus allowing flexural-phonon-electron coupling. The method is also interfaced with the Quantum ESPRESSO package³⁰ for performing linear response calculations. Closely related to the truncated Coulomb approach is the dipole correction method by Brumme et al.,³⁴ which introduces a way to include a charged plate within a system with periodic boundary conditions. In the first-principles calculations, it is common to simply apply an electric field to understand field effects. Examples include, but are not limited to, graphene-based 2D heterostructures.^35–37 It should be pointed out that applying an electric field to a system is very different from gating the system, since the former can be described as a closed system but the latter is certainly an open system. Experimentally, both gating and electric fields are used, and it is important to choose an appropriate tool to address the right problem. Figure 1 is a sketch of materials and physical properties investigated using various methods by the above-mentioned groups. The remainder of this paper is organized as follows: in Sec. II, we provide a brief description of our approach, models, and Poisson solver, and a review of prior applications; in Sec. III, we present results from recent studies. We conclude our effort of modeling field effects using the first-principles method with an outline of future developments and applications.

Similar to most other recent theoretical studies of field effects, we focus on two-dimensional systems, which comprise a very active research area. In our approach, experimental conditions are first identified, according to which we construct our simulation models. Figure 2 depicts single and dual gate configurations. For both the top and bottom panels, the left part is a sketch of the experiment, and the right is our model for simulation. Two key approximations in the simulation models are: (1) the dielectric layer is kept thin to reduce computational effort; and (2) a vacuum layer is inserted between electrodes (metal) and dielectrics to avoid complications caused by dielectric–metal interfaces, which have very little effect on the physical processes in the 2D systems.

The potential is the solution to the Poisson equation $∇·[ϵ(r→)∇]V(r→)=−4πρtot(r→)$⁠, and the Green's function is defined as $∇·[ϵ(r→)∇]G(r→,r→′)=−4πδ(r→−r→′)$⁠. With the given boundary condition, the Green's function is the potential of a point charge plus its image,

In Eq. (1), G is expanded in the momentum $g→$ within the 2D plane; z is the dimension in which the gate is applied, or the direction perpendicular to the plane, and z₁ is the position where the gate voltage is applied, or the position of the electrode. Integrating G over $z′$ gives the electrostatic potential at any given point z. For the dual-gate configuration, there are two electrodes, which lead to infinite number of image charges; however, we know the analytical expression of the summation. The potential in the x-y directions is easy to obtain in momentum space, since there are no image charges. This is the essence of the effective screening medium technique.²² In the framework of density functional theory,^45,46 the total energy of the system is

where n is the electron density, K is the single-particle, non-interacting electron kinetic energy, E_xc is the exchange-correlation functional, and n_I is nuclear charge distribution. When analyzing interface properties, it is desirable to quantify the electric polarization as a function of the distance measured from the interface. For this purpose, we implemented the so-called hybrid Wannier function in our analysis code. The conventional Wannier function is defined by the Fourier transform of the Bloch wave,

The application of Wannier functions in solid state physics is now a common practice because of the availability of the Wannier90 package.⁴⁷ However, the Wannier orbitals obtained via the 3D transformation are not adequate for interfaces between a metal and an insulator because of the delocalized wavefunction in the plane of the interface, and the procedure will not converge. The hybrid Wannier function technique was proposed to overcome this obstacle.⁴⁸ Along the z-axis, we calculate $Mmn(k→)= ⟨um,k→ |un,k→+b→∥ ⟩$⁠, where $um,k→$ is the Bloch wave without the propagation exponential, m and n are band indices, $k→$ is the crystal momentum, and $b→∥$ is the $k→$-spacing in the z-direction. A global matrix $Λ(k→⊥)=∏j=0N∥−1M∼(k→⊥+jb→∥)$ is constructed whose eigenvalues λ_m are related to the center of the Wannier function in the z-direction by $zm=−(L/2π) Im(ln λm)$⁠, and $k⊥$ is the crystal momentum $k→$ in the x-y plane.²⁰ The polarization P in the z-direction for each $k→⊥$ is finally written as,

where the first term is the electron contribution to P and the second term is from ionic displacements, where $Δzm,k→⊥$ is the difference between the ionic coordinate and the center of the Wannier orbital m, $Qα$ is the charge of ion α, and $Δzα$ is the ionic displacement in the z-direction.

When density functional theory (DFT) + U is used in the calculations, the value of U is often taken from the literature or sometimes used as adjustable parameters. If one would like to get a first principles estimate of U for real materials, the constrained random phase approximation (cRPA) method can provide a fully quantum mechanical parameterization of U based on the DFT ground state. The basic idea of cRPA^49,50 is to calculate a partial RPA particle-hole polarization with the constraint of a physically motivated correlation window (e.g., the d-like bands of transition metal atoms). One aims to estimate the screened Coulomb interaction for the correlation window. For this purpose, the particle-hole polarization between all possible pairs of occupied state and unoccupied state is taken into account. Within RPA, the full particle-hole polarization can be written as⁵¹ 

where ψ_i and $εi$ are the single particle eigenfunctions and eigenenergies of DFT. Summations on i and j are restricted such that i must be an occupied state and j must be an unoccupied state.

The selected bands in the correlation window often have a strong orbital character, e.g., d-like in our case. Following the convention in the literature, these bands are called d-space. If both the occupied and unoccupied states are within the d-space, then the polarization contributes to $PdRPA(r,r′;ω)$⁠. All other pairs of occupied and unoccupied states contribute to $PrRPA$⁠, where r stands for the rest of the polarization. Thus, the full polarization is divided into two parts: $PRPA=PdRPA+PrRPA$⁠. Here, the partial polarization $PrRPA$ is the quantity related to the partially screened Coulomb interaction,⁴⁹ 

In this expression, v is the bare Coulomb interaction. According to the Hedin equations and the GW approximation, the full polarization, $PRPA$⁠, screens the bare Coulomb interaction, v, to give the fully screened interaction W. With the same logic, $PdRPA$ screens W_r to give the fully screened interaction W. Thus, W_r is identified as the screened on-site Coulomb interaction for the d-space, i.e., $U(ω)≡Wr(ω)$⁠, that includes the screening effect from the realistic environment of the material. In practice, one calculates the partial polarization $PrRPA$ from the Kohn–Sham susceptibility, which is completely based on the DFT ground state, and then derives $U(ω)$ from $PrRPA$⁠.

The graphene|h-BN|graphene heterostructure²⁰ was the first application of our approach, based on the single-gate system (as shown in Fig. 2, top panel) that was studied experimentally.⁶ In this work, the gate effect at the h-BN and graphene interface was fully analyzed. In order to see whether the immediate contact between graphene and h-BN makes a difference, layer-by-layer electric polarization analysis was performed for h-BN. The electric polarization for each h-BN layer is calculated by $P=−∑ieδzi/Sd$⁠. Here, e is the unit charge, $δzi$ is the change in the hybrid Wannier charge center upon applying a gate voltage in the direction of the gate field, S is the area of the unit cell, and d is the thickness of a h-BN layer which is set to $3.33 Å$⁠. The summation is over all hybrid Wannier functions belonging to this h-BN layer within a unit cell. The layer-by-layer electric polarization analysis shows that the first h-BN layer in direct contact with the graphene sheet has a polarization similar to those further away from the interface (see Fig. 3), and this curve is also similar to that for the bare five-layer h-BN system (not shown). Its inert nature makes h-BN a perfect choice of a supporting material for graphene, let alone that the lattice mismatch is very small. It is clear that compared to h-BN, H-terminated Si shows much stronger interface effects (for simplicity, we do not insert the Si slab between the graphene sheets). We also attempted to compute the transmission function as a function of gate voltage. Due to computational limitations at that time, only monolayer h-BN was considered to illustrate the point (see Fig. 4). However, the model we built is good for other investigations.

As expected from experimental measurements, first-principles calculations also show a Fermi energy shift and gap opening. When there is hole doping, the graphene layer closer to the electrode is doped more than the one further away from the electrode, which is not surprising.

Trilayer graphene was studied experimentally in the dual-gate configuration that inspired us to complete our implementation.²¹ We examined field effects on both ABA stacking and ABC stacking orders, and our calculations reproduce the experimentally observed^52,53 gap opening in ABC stacking and band overlap in ABA stacking. The dual-gate configuration allows one to investigate effects of doping, and our calculations predict possible gap reopening upon doping in the ABA stacking as shown in Fig. 5. We suggest that infrared optical conductivity measurements can confirm the calculated band gap reopening.

This work was motivated by experimental studies of heterogeneous junctions that consist of two graphene sheets bracketing a monolayer of azobenzene molecules.⁵⁴ The azobenzene molecule has two stable configurations, trans and cis, that can transform from one to another by optical excitation.^55,56 We showed that these two forms of the molecule have different transport properties at zero bias and different I–V characteristics in the one-dimensional configuration.^57,58 Investigations of an azobenzene monolayer between two semi-infinite Au bulk leads indicate that the transport properties of the trans and cis molecules are more complicated.⁵⁹ Our analysis shows that chemisorption of molecules to the top Au lead to a different trans vs cis relation than a physisorption of molecules to the top Au lead. In addition, an ad-atom on an Au surface can change dramatically the I–V characteristics, which explains the experimentally observed I–V curve of an Au–azobenzene–Au break junction.⁶⁰ We further applied azobenzene molecules to modulate the interaction between the two nano-particles and studied the interaction with graphene surfaces in the two configurations.^61,62 These studies demonstrated that one can use the configuration change to manipulate physical properties of a system. Simulations of the graphene|azobenzene|graphene vertical junction are therefore a natural extension of our long-standing interest in the added technique for gating the system.²⁶ We found a rich array of interesting phenomena. The first noticeable finding is that, depending on the sign of the gate field, the trans and the cis become more conducting only in one direction; and second, at some gate voltages, two peaks appear near the Fermi energy. Our analysis shows that gate voltage alters the energy levels in such a way that the interface state (the C–C bond between a molecule and graphene) moves closer to and the Dirac point of the top layer of graphene moves away from the graphene Fermi energy. Even more amazing, one interface state can interfere with another one, resulting in a second, stronger peak, which disappears when we reduce the coverage to 50%. Figure 6 shows the interface state for the trans and the cis molecules in the vertical direction. The transmission functions in the panel (e) clearly show that at full coverage a very strong second peak appears as a result of interference of neighboring C–C interface states.

Transition metal dichalcogenides (TMD) make up one group of 2D semiconductors that have attracted much attention^11,63–65 in the quest of TFETs. Compared to h-BN, the relatively smaller energy gap and strong spin–orbital coupling make them more interesting than h-BN, which has otherwise been a perfect choice for graphene supporting material. We studied field effects of graphene $| n$-layer TMD| graphene (n = 1–5, TM = W, Mo) junctions. Figure 7 compares WS₂ and MoS₂. Two critical quantities, the distance $ΔE$ between the conduction band edge and the Dirac point (the Fermi level) and band splitting $ESO$ due to the spin–orbital coupling, were computed as a function of charge doping. The decrease in $ΔE$ is responsible for the large on/off ratio observed reported by the experimental groups.^5,66 For junctions with thinner WS₂ (n = 1, 2) or MoS₂ (n = 1), $ΔE$ is symmetric between hole doping and electron doping due to the symmetry in the Dirac cone to which electrons/holes are added. For junctions with thicker WS₂ or MoS₂, holes are still added to the graphene only but electrons are added to both the graphene and the TMD, which causes the asymmetry in $ΔE$⁠. The gate field has little effect on $Eso$⁠. This is because both the two valence bands at K are mostly d states from transition metal atoms.⁶⁷ As such, the shift in the two bands is the same when a gate field is applied. In comparison, the junction with n-layer MoS₂ has a much smaller $ΔE$ than the junction with n-layer WS₂. Specifically, $ΔE$ is zero for the graphene $| 4$-layer MoS₂| graphene junction even under zero charge doping.

Bi-Sb-Te-Se₂ (BSTS) is a strong topological insulator with high bulk resistivity and robust surface states. It is promising for applications in high speed and low-energy-consumption spintronic devices. One of the device configurations is the vertical tunneling junction, where two BSTS slabs are stacked together. A question arises for such a configuration that whether topological surface states survive at the interface between two BSTS slabs. First-principles calculations show that topological interface states are absent at the equilibrium inter-slab distance but can be preserved by inserting two or more layers of h-BN between the two BSTS slabs.²⁸ Furthermore, experimental measurements revealed a weak dependence of the electron tunneling current on the gate field at small bias voltages for a BSTS vertical tunneling junction.²⁸ The dependence at high bias voltages is however stronger. Then, one may wonder how topological interface states (when present) respond to a gate electric field. In order to examine this problem, first-principles calculations were performed for the BSTS interface with bilayer h-BN.

Figure 8(a) illustrates the BSTS|bilayer h-BN|BSTS interface under the influence of two gates. A dual gate setup permits independent control of the charge doping concentration n and the average electric field between the two gate electrodes $E$⁠. Figures 8(b)–8(e) show the band structure of the interface under increasing electron doping levels but zero average electric field. There are four species of Dirac states in the system, the top surface species, the bottom surface species, the top interface species, and the bottom interface species, which are represented by red empty circles, blue empty squares, red filled circles, and blue filled squares, respectively. The larger a circle (a square) is, the more localized at the corresponding surface or interface the state is. At zero doping and zero average electric field, all four Dirac points are about the Fermi level, which is set to zero. When electrons are added to the system, the top and the bottom surface bands move downward while the top and the bottom interface bands barely move, which remains true up to a doping concentration of $−3.28×1013 cm−2$⁠. This means that the added electrons mainly go to the top and the bottom surfaces, which helps to understand the weak dependence of the tunneling current on the gate field at small bias voltages. As the doping concentration increases, the conduction bands of bulk BSTS, which are above the Dirac cones in Fig. 8(b), get closer to the Fermi level. Meanwhile, they extend more into the surface regions in real space. Eventually, the otherwise empty conduction bands become partially occupied at a doping concentration of about $−3.28×1013 cm−2$⁠. At a constant charge doping level of $E≥0.5V/Å$⁠, we applied different average electric fields and obtained the band structures shown in Figs. 8(f)–8(i). The major effect of such an electric field is to separate the top surface bands from the bottom surface bands. Specifically, the top (bottom) surface bands are moved upward (downward), meaning that some electrons are transferred from the top surface to the bottom surface. This is consistent with the intuition that a positive electric field along the $E=0.5V/Å$ direction applies a force in the −z direction on electrons. Again, the interface surface bands are not much affected. The conduction bands of bulk BSTS are also shifted downward by a positive electric field, and they cross the Fermi level at $E=0.4 V/Å$ for $0.3V/Å$⁠. Such a doping of the bulk states under a finite average electric field may correlate with the stronger dependence of the tunneling current on the gate field at high bias voltages.

Planar geometry can be viewed as a special case of vertical geometry of minimal thickness, which simplifies the simulation model and reduces computational cost. Here, we highlight a few systems that have been looked at in the last few years.

Metal phthalocyanine (MPc) is planar molecules of nanometer size. Abel et al.⁶⁸ attempted to synthesize a covalently bonded 2D framework using MPc molecules. Although it was later proven to be a 2D hydrogen-bonded network, covalently bonded 2D $0.6V/Å$ MnPc networks were realized on Ag surfaces,⁶⁹ and 1D micrometer FePc wires have also been synthesized.⁷⁰ We reported field-effect studies of magnetic order in 2D MPc networks²⁷ and spin-dependent charge transport of MPc junctions.⁴⁴ Figure 9 sketches a general metal phthalocyanine molecule, in which the transition metal ion transfer is at the center. We examined Cr, Mn, and Fe systems and find that each of them loses two electrons to the organic framework. Interestingly, for Cr²⁺- and Fe²⁺-doped Pc, the d-orbital split is large enough that the 2D framework is semiconductors with a band gap, while 2D MnPc is a half metal because the MnPc has d orbitals that are partially occupied (see the middle panel of Kohn–Sham levels in Fig. 9). More importantly, we found that the magnetic order of the FePc, MnPc, and MnTPP [TPP: 5,10,15,20-tetra(phenyl)porphyrin] can be tuned by charge doping (via gating), and the interactions between two spins are mediated by itinerant electrons. Based on these findings, we designed planar MnPc|NiPc|MnPc 2D junctions and investigated (1) gate effects and (2) scattering region length dependence. It is found that this system can be used as a perfect spin filter when it is hole-doped.⁴⁴ Further analysis shows that the role of gating is to align the energy level of the lead (MnPc) with the scattering region (NiPc), such that a conducting channel can appear (see Fig. 10).

The idea of patterning graphene or a graphitic material into functioning circuitry has been a scientific and engineering focus since the time of discovery of carbon nanotubes⁷¹ and single-layer or few-layer graphene.⁷² 

We carried out first-principles investigations of a double barrier in the graphene framework. Figure 11 depicts the model for a double barrier that consists of h-BN. The zigzag edge was chosen and the interface composition was discussed: on one side, we chose the C–B bond and on the other hand, the N–C bond. Vacuum provides another system, in which the zigzag graphene edge was terminated by H atoms. With DFT calculations, details such as inter-edge spin couplings can be considered. The transmission coefficient function (see the left panel of Fig. 12) shows a shift in the rising/falling transmission coefficient between the two spin channels, indicating that spin filtering can be achieved by a bias voltage. Band analysis (see the right panel of Fig. 12) indicates that this originates from the difference between the bands (of the middle graphene ribbon) of the two spins. Only when the graphene ribbon bands intersect with the Dirac cone can the transmission coefficient rise to a non-negligible value. We call this a resonance, which is even more pronounced if the double barrier is made of vacuum. It is interesting to see that the center graphene ribbon is conducting, with sizable contribution from the interface state as the bands cross the Fermi level. Such interface enhancement or induced conducting behavior was observed and characterized in another study from our group where we found that a quasi-1D conducting wire can form at the interface of two semiconductors.²³ 

Other than previously studied materials, many more two-dimensional magnetic materials can be employed in spintronic devices, such as high capacity information storage. Recent discoveries of 2D magnetic material include FePS₃,⁷³ Cr₂Ge₂Te₆,⁷⁴ MnSe₂,⁷⁵ and CrI₃.⁷⁶ Bulk CrI₃ is a layered van der Waals material with a ferromagnetic order at low temperatures (LTs). It has a high temperature (HT) monoclinic phase with space group $E=0.9V/Å$ and a low temperature (LT) rhombohedral phase with space group $0.3V/Å$⁠.⁷⁷ The two phases differ in the interlayer stacking, as shown in Figs. 13(a) and 13(b). Interestingly, bilayer CrI₃ exhibits antiferromagnetic (AFM) interlayer coupling,⁷⁶ which can be tuned by a magnetic field,^76,78,79 gate electric field,^78,79 and pressure.^80,81 In efforts to explain the experimentally observed AFM magnetic order, Sivadas et al. reported stacking-dependent magnetism⁸² and Jang et al. analyzed the interaction between localized e_g and $E>0.3V/Å$ orbitals.⁸³ Among these and other efforts, we will examine the role of the local Coulomb interaction in Subsection III A. In some experimental setups to apply a gate field,^78,79,84 bilayer CrI₃ is in contact with graphene or hexagonal boron nitride (h-BN), which motivated us to model gate field effects on heterostructures such as graphene|bilayer CrI₃|graphene, BN|bilayer CrI₃|BN, and BN|bilayer CrI₃|graphene. These results will be presented in Subsection III B, emphasizing the interface effects on the magnetic phase transition. For these heterostructures, we will denote bilayer CrI₃ by 2-CrI₃ for brevity. We will also discuss how pressure affects structural and magnetic properties of bilayer CrI₃ in Subsection III C.

We relaxed the atomic structure of bilayer CrI₃ based on density functional theory^45,46 as implemented in the Vienna Ab initio Simulation Package (VASP).⁸⁵ The computational details can be found in Ref. 86. Then, we applied the DFT $+U$ method⁸⁷ (the details of estimating the relative U values will be discussed shortly) in VASP to obtain the total energies of bilayer CrI₃ in both the HT and LT stacking. Figure 14(c) shows the energy difference between the HT and LT stackings versus the Hubbard U parameter. The energy difference changes sign as the strength of the local Coulomb interaction increases. If U is small, the LT stacking is energetically preferred, which is the case for bulk CrI₃. If U is larger than $7 eV$⁠, the HT stacking has lower energy than the LT stacking. With such a large value of the U parameter, the HT stacked bilayer CrI₃ energetically prefers the AFM configuration to the ferromagnetic (FM) configuration, as shown in Fig. 14(d). In contrast, the LT stacked bilayer CrI₃ always prefers the FM configuration, no matter how large or small the U parameter is. Therefore, the HT stacking together with a pronounced local Coulomb interaction seems to be responsible for the experimentally observed interlayer AFM magnetism.

It is generally accepted that magnetic moments are associated with localized electrons whose behavior is determined by the competition between the kinetic energy of electrons and the strength of the local Coulomb interaction. Dimensional confinement can dramatically affect the competition, yielding unusual properties. For example, bulk SrVO₃ is known to be a strongly correlated metal. In standard DFT calculations, the isolated three $t2g$ bands of vanadium determine the low-energy properties of this material. In DFT + DMFT dynamical mean-field theory calculations, the Coulomb interactions are explicitly taken into account and the ground state of bulk SrVO₃ is still found to be a metal (with a smaller band width). When the system is under dimensional confinement, orbital re-occupation can happen along with an enhanced Coulomb interaction. In a charge density self-consistent (CSC) DFT + DMFT^88–90 study, 2D monolayer SrVO₃ was found to be insulating,⁹⁰ with a band gap of about $2 eV$⁠. Actually, already at the DFT level, the in-plane d_xy band is no longer degenerate with the other two $t2g$ bands. After DMFT and charge density consistency, the d_xy orbital is found to be half-filled, and the other two out-of-plane $t2g$ bands are almost empty. The change from a bulk correlated metal state to a monolayer Mott insulating state is accompanied by an increase in Coulomb interaction from about 4  to 5.5 eV.

In order to confirm that the local Coulomb interaction increases as the dimension of CrI₃ reduces from 3D to 2D, we calculated the U-matrix within the Kanamori parameterization from first principles using the constrained random phase approximation (cRPA) method as implemented in the FP-LAPW DFT code, a modified version of ELK code.^91,92 The code has been benchmarked^93,94 with other implementations using late transition metal monoxides, and we obtain consistent results. In our calculation, both bulk and monolayer CrI₃ have isolated Cr d-like bands around the Fermi level. We choose the five d-like bands as our correlation window. The ground state includes 100 empty bands to make a reasonable estimation of the partial particle-hole polarization, $PrRPA$⁠. The resulting averaged intra-orbital U increases from about $1.9 eV$ for bulk CrI₃ to about $2.8 eV$ for the monolayer CrI₃. The calculated U parameters are not large enough to bear the AFM ground state of bilayer CrI₃. However, one should keep in mind the cRPA calculation is based on paramagnetic ground states of the two structures, and the Pauli exclusion principle (which is a different mechanism to give rise to on-site repulsion between electrons, especially for magnetic systems) is not taken into account. It is still a non-trivial job to incorporate the Pauli principle with the current cRPA method in one calculation scheme. Here, our observations confirmed an enhancement of about 1 eV in the Coulomb interaction when the structure of CrI₃ reduces from 3D to 2D.

Figure 14(a) illustrates a graphene|2-CrI₃|graphene heterostructure subject to a dual gate setup, which permits a vertical electric field and charge doping. For such a dual gate setup, the average electric field between the gate electrodes will be what we call the electric field $E$⁠,

where $VTG/VBG$ is the electrostatic potential of the top/bottom gate and L is the distance between the two gate electrodes. A metallic part of the heterostructure, which is graphene in the case of Fig. 14(a), is considered to be grounded so that extra charge can be introduced from the environment to the heterostructure. The electric field $E$ and the extra charge density n can be viewed as two independent variables for a dual gate setup. If the heterostructure is insulating, extra charges can hardly be added to the system due to the lack of states around the Fermi energy. Therefore, we will fix the extra charge density n to be zero and tune only the electric field $E$ in our simulations.

We study gate field effects on bare bilayer CrI₃, h-BN|2-CrI₃|h-BN, h-BN|2-CrI₃|graphene, and graphene|2-CrI₃|graphene. Our calculations are based on density functional theory in conjunction with the effective screening medium method as implemented in the SIESTA package.⁹⁵ The computational details of gate calculations with SIESTA can be found in the Ref. 96. The unit cell of the heterostructures is as large as a $3×3$ supercell of the bare bilayer CrI₃, or a 5 × 5 supercell of graphene. Graphene (BN) is compressed by 0.7% (1.1%) to fit with the lattice constant of bilayer CrI₃. The atomic structures of the heterostructures were relaxed by VASP using the same parameters for the bare bilayer CrI₃ except that the force tolerance is set to 0.02 eV/Å and the k-point mesh is $5×5×1$⁠. Here, we consider only the HT stacking bilayer CrI₃, since it is likely what is seen in experiments. First, we consider the effects of electric field under the condition of zero charge doping. Figure 14(b) shows the energy difference between the AFM and the FM magnetic configurations for these systems. At zero electric field, the energy of the AFM state is lower than that of the FM state by 15–$18 meV$ per unit cell (of bilayer CrI₃).⁹⁷ As electric field increases, a magnetic phase transition from the AFM state to the FM state occurs at $∼0.6 V/Å$ for bare bilayer CrI₃. Such an AFM-to-FM magnetic phase transition was also reported in previous experimental⁷⁸ and theoretical⁹⁸ studies. If bilayer CrI₃ is covered by graphene on both the bottom and the top sides, the AFM state is always energetically preferred. Actually, the energy difference $EAFM−EFM$ is less than $−10 meV$ up to an electric field of $0.9 V/Å$⁠. The magnetic phase transition is also absent for the heterostructure h-BN|2-CrI₃|h-BN. However, the AFM and the FM state are quite close in energy at an electric field of $0.4 V/Å$⁠. Immediately after $E=0.4 V/Å, EAFM−EFM$ decreases and reaches about $−6 meV$ at $0.5 V/Å$⁠. If bilayer CrI₃ is covered by BN on the bottom side and graphene on the top side, an AFM-to-FM magnetic phase transition was calculated to occur at around $0.8 V/Å$⁠. Since the heterostructure h-BN|2-CrI₃|graphene is asymmetric in the out-of-plane direction, it is sensitive to the direction of electric field. Based on our calculations, a negative electric field also tends to stabilize the FM state; however, it does not induce any magnetic phase transition down to $−0.9 V/Å$⁠. Second, we consider the effects of charge doping under the condition of zero electric field (⁠$VTG=VBG$⁠). Figure 14(c) shows the energy difference $EAFM−EFM$ versus the extra charge density for h-BN|2-CrI₃|graphene and graphene|2-CrI₃|graphene. These two systems are always in the AFM state under both electron doping and hole doping conditions. The energy difference $EAFM−EFM$ varies between $−20 meV$ and $−15 meV$ for the doping level range of $−1013:+1013 cm−2$⁠.

So far, we have examined the magnetic phase transition in pure and hybrid bilayer CrI₃ systems via total energy calculations. Next, we will explain the interfacial effects on the energy diagrams of Figs. 14(b) and 14(c) by a detailed electronic structure. Figure 15(a) shows the energy bands of the HT stacking bilayer CrI₃ in both AFM and FM states. The purple circles highlight the conduction band E_c and the valence band E_v at the Γ-point. Figure 15 shows the energy difference between E_c and E_v versus electric field. Under zero electric field, $Ec−Ev$ is $0.45 eV$ for the AFM state and $0.38 eV$ for the FM state. The difference decreases as the electric field increases for both states and reaches zero at an electric field of $∼0.3 V/Å$ for the AFM state. The slope of the $Ec(E)−Ev(E)$ curve for the AFM state changes significantly at $E ∼0.3 V/Å$⁠. This corresponds to the direct band gap closing for the AFM bilayer CrI₃, which is shown in Fig. 15(c). In Fig. 15(c), we see that both the conduction band (the spin-up electrons of the bottom CrI₃ layer) and the valence band (the spin-down electrons of the top CrI₃ layer) touch the Fermi level. Past $E=0.3 V/Å, Ec−Ev$ decreases almost linearly with the electric field. The FM bilayer CrI₃ experiences an indirect band gap closing at $E ∼0.5 V/Å$⁠, where both the conduction and the valence electrons at the Fermi level are spin-up electrons, see Fig. 15(d). The band gap closing seems to be correlated with the plateau of the $EAFM(E)−EFM(E)$ curve between $0.3$ and $0.5 V/Å$⁠. Figure 15(e) [Fig. 15(f)] shows the band structure of h-BN|2-CrI₃|h-BN under an electric field of $0.4 V/Å$ (⁠$0.5 V/Å$⁠). The valence band of the top BN layer crosses the Fermi level at $E=0.5 V/Å$ but not at $E=0.4 V/Å$⁠. Such a band crossing is likely the reason for the significant reduction of the energy difference between the AFM and the FM states, as depicted in Fig. 14(b). Similarly for the heterostructure h-BN|2-CrI₃|graphene, the bottom BN layer starts to lose electrons to the top graphene layer at an electric field around $−0.4 V/Å$⁠. The electron transfer results in an induced electric field which is opposite to the direction of the gate electric field, and thus it weakens the net electric field across the bilayer CrI₃. As a result, the AFM-to-FM magnetic phase transition is hindered.

In order to examine the electron redistribution of the systems in response to a gate field, we plot the plane-averaged electron density difference $ρ(E,n)−ρ0$ in Fig. 16, where ρ₀ is the electron density without any gate field. Since the averaged electron density difference for the AFM and the FM states is quite similar, we will present only that for the AFM state. Figures 16(a)–16(d) show the effects of electric field for bilayer CrI₃, h-BN|2-CrI₃|h-BN, h-BN|2-CrI₃|graphene, and graphene|2-CrI₃|graphene, respectively, under zero charge doping. For bilayer CrI₃, the major electron transfer is from the topmost iodine atomic layer to the bottommost iodine atomic layer. This electron transfer increases gradually from $E=0.1$ to $E=0.9 V/Å$⁠, which does not signal the band gap closing at around $E=0.3 V/Å$⁠. In comparison, the inner chromium and iodine atomic layers experience a relatively small change in electron density. Specifically, the band gap closing is signaled by the change in the electron density between the inner two CrI₃ layers. This electron density change increases gradually at small electric field but saturates after the band gap closing, as shown in the inset of Fig. 16(a). For h-BN|2-CrI₃|h-BN, the electron transfer is similar to that of bare bilayer CrI₃ when the electric field is smaller than $0.4 V/Å$⁠. Within this range of electric field, a local electronic dipole forms for each BN atomic layer without significant electron transfer between BN and CrI₃. Consequently, the energy difference between the AFM and the FM states for h-BN|2-CrI₃|h-BN is quite close to that for bare bilayer CrI₃ when $E<0.4 V/Å$ [see Fig. 14(b)]. However, the major electron transfer is from the top-most BN atomic layer and the bottom-most iodine atomic layer after $E=0.4 V/Å$⁠, which is consistent with the band structure in Fig. 15. The inset of Fig. 16(c) shows the electron density variation for h-BN|2-CrI₃|graphene under $0.1 V/Å$ electric field, a local electronic dipole also forms around the top graphene layer rendering a small amount of electron transfer between graphene and the remaining insulating part of the heterostructure. Again, the energy difference between the AFM and the FM states of h-BN|2-CrI₃|graphene is close to that of bare bilayer CrI₃ at $E=0.1 V/Å$⁠. As the electric field further increases, the electron transfer between the top graphene layer and the bottommost iodine atomic layer gradually becomes dominant. This behavior is correlated with the observation that the magnetic phase transition occurs at a larger electric field for h-BN|2-CrI₃|graphene compared with bare bilayer CrI₃. In contrast to the local dipole formation around graphene in case of $0.1 V/Å$ electric field, a $−0.1 V/Å$ electric field results in a significant amount of electron transfer between the bottom CrI₃ layer and the top graphene layer. When bilayer CrI₃ is covered by graphene on both the top and the bottom sides, the inter-CrI₃-layer electron transfer is greatly reduced due to the electrostatic shielding of the graphene layers. As a result, there is no magnetic phase transition up to an electric field of $0.9 V/Å$⁠. Figures 16(e) and 16(f) show the effects of charge doping for h-BN|2-CrI₃|graphene and graphene|2-CrI₃|graphene, respectively, under the condition of zero electric field (⁠$VTG=VBG$⁠). Upon the addition of electrons or holes, the extra charges go to both the top and the bottom graphene layers resulting in a small or even negligible inter-CrI₃-layer electron transfer. This is the major difference from the case of applying electric field, where there is significant inter-CrI₃-layer electron transfer. The lack of significant inter-CrI₃-layer electron transfer seems the reason why $EAFM(n)−EFM(n)$ does not change much with extra charge density n. Furthermore, there is still some inter-CrI₃-layer electron transfer of h-BN|2-CrI₃|graphene due to the asymmetry between the BN and the graphene layers. This inter-CrI₃-layer electron transfer, although small by itself, is larger than that of graphene|2-CrI₃|graphene. Such a comparison could explain why $EAFM(n)−EFM(n)$ of h-BN|2-CrI₃|graphene is slightly larger than $EAFM(n)−EFM(n)$ of graphene|2-CrI₃|graphene as shown in Fig. 14(c).

Figure 17(a) [Fig. 17(c)] shows the magnetic moment M of bare bilayer CrI₃, h-BN|2-CrI₃|h-BN, h-BN|2-CrI₃|graphene, and graphene|2-CrI₃|graphene systems in the AFM (FM) state versus electric field. M is measured per unit cell of bilayer CrI₃ with four chromium atoms. For AFM bare bilayer CrI₃, the magnetic moment remains zero until $E>0.3 V/Å$⁠, where the band gap closes and inter-spin electron transfer transpires, as shown in Fig. 15(c). The magnetic moment M increases linearly with the electric field after $0.3 V/Å$ and reaches $∼0.4 μB$ per unit cell at $E=0.9 V/Å$⁠. An AFM-to-FM magnetic phase transition occurs at around $0.6 V/Å$⁠, where the magnetic moment changes from $∼0.2μB$ to $12μB$⁠. M is insensitive to the electric field for FM bilayer CrI₃ because there is no inter-spin electron transfer even if the band gap closes, as shown in Fig. 15(d). For AFM h-BN|2-CrI₃|h-BN, the magnetic moment also becomes finite after $0.4 V/Å$ for the same reason as with the case of bare bilayer CrI₃, which is evidenced by the band structure of AFM h-BN|2-CrI₃|h-BN under $E=0.4 V/Å$ shown in Fig. 15(e). After $E=0.5 V/Å$⁠, the magnetic moment of AFM h-BN|2-CrI₃|h-BN is smaller than that of AFM bilayer CrI₃. This can be understood since the total inter-layer electron transfer of h-BN|2-CrI₃|h-BN consists of both intra-spin and inter-spin electron transfer, while that of bilayer CrI₃ consists merely of inter-spin electron transfer. The former shows a smaller change in the magnetic moment than the latter under the assumption that the total inter-layer electron transfer is the same, which should be a good approximation for these two systems under the same electric field. The magnetic moment of FM h-BN|2-CrI₃|h-BN does not increase until $E≥0.5 V/Å$⁠, which is due to the electron transfer from the spin down channel of the top BN layer to the spin up channel of the bottom CrI₃ layer, as shown in Fig. 15(f). For h-BN|2-CrI₃|graphene under zero electric field, graphene is slightly doped with electrons regardless of the magnetic state, which can be seen from the band structures in Figs. 17(e) and 17(g). Beginning at $E=0.3 V/Å$⁠, the magnetic moment increases significantly for both the AFM and the FM states. This is because graphene loses electrons from both the spin up and the spin down channels to the spin up channel of CrI₃ as shown in Figs. 17(f) and 17(h). For graphene|2-CrI₃|graphene under zero electric field, both the top and bottom graphene layers gain the same amount of electrons from CrI₃ as evidenced by the same shift of the Dirac cones in Figs. 17(i) and 17(k) for the AFM and the FM states, respectively. The magnetic moment of AFM graphene|2-CrI₃|graphene under zero electric field is exactly zero μ_B due to the spin degeneracy. In contrast, the magnetic moment of FM graphene|2-CrI₃|graphene under zero electric field is not exact $12 μB$⁠, which is the value for bare bilayer CrI₃. This is because the valence bands of the FM bilayer CrI₃ are spin split by $∼0.2 eV$ and only the spin-up energy bands are hole doped. Figure 17(j) [Fig. 17(l)] shows the band structure of the AFM (the FM) graphene|2-CrI₃|graphene under an electric field of $0.3 V/Å$⁠. The major change in the band structure is the shift of graphene bands that corresponds to the process of the top graphene layer losing electrons to the bottom graphene layer. A finite electric field breaks the spin degeneracy of the AFM state leaving the spin-up valence band fully occupied and the spin-down valence band slightly doped with holes. As a result, the magnetic moment of the AFM graphene|2-CrI₃|graphene state under finite electric field is slightly above zero μ_B. In contrast, a finite electric field slightly enhances the hole doping of the spin-up valence band of FM graphene|2-CrI₃|graphene and, in consequence, the magnetic moment reduces a bit. The inset of Fig. 17(e) [Fig. 17(i)] shows a zoomed-in plot of the band structure of h-BN|2-CrI₃|graphene (graphene|2-CrI₃|graphene) in the energy range of $−30:10 meV$ at the Γ point. Note that the energy bands are spin degenerate for graphene|2-CrI₃|graphene but not for h-BN|2-CrI₃|graphene. Such a comparison indicates that the spin degeneracy of bare AFM bilayer CrI₃ is preserved (broken) by the symmetrical (asymmetrical) surrounding chemical environment.

Recently, it has been experimentally observed⁸⁰ that the state of bilayer CrI₃ at zero external magnetic field switches from an interlayer AFM state to an interlayer FM state on increasing the pressure. A corresponding structural transition from HT stacking to LT stacking has also been confirmed by Raman spectroscopy. Interestingly, having experienced an increase in pressure up to $2.70 GPa$ and then taken out of the pressure cell, the bilayer CrI₃ sample remains in the LT structure and the FM state. This implies that the structural transition that occurs during the pressurizing stage is irreversible.

In order to study a possible structural transition and to investigate the variation of magnetism in each state, we have calculated the Gibbs free energies of the FM state in the LT structure (LT-FM) and the AFM state in the HT structure (HT-AFM) as functions of pressure. Our calculations are based on density functional theory as implemented in the VASP package with the same computational details as in Ref. 84. To simulate the effect of pressure, the distance d between the topmost monoatomic and bottommost monoatomic layers of I atoms [see the inset of Fig. 18(a)] has been tuned, while the in-plane coordinates of these I atoms, all coordinates of the other atoms, and the lattice constant a of the hexagonal unit cell are fully relaxed. In the end, at each d the pressure is calculated by summing up the atomic forces on the topmost or bottommost monoatomic layer of I atoms and then dividing the result by the in-plane area of the unit cell. In Figs. 18(a) and 18(b), lattice-related parameters d and a are plotted, respectively. As expected, d, as a measure of interlayer distance, decreases with increasing pressure, while the lattice constant a increases. The Gibbs free energy difference between the LT-FM and HT-AFM states, $ΔG=GHT−AFM−GLT−FM$⁠, obtained from the linearly interpolated Gibbs free energy $G=E+pV$ as a function of p in each state, is shown in Fig. 18(c). Here, E is the total energy of the state, p is the pressure, and the volume is V = Ad, where $A=3 a2/2$ is the in-plane area of the unit cell. The free energy shows that the HT-AFM state is less stable than the LT-FM state at all experimental pressures and becomes increasingly unstable compared to the LT-FM state with increasing pressure. It worth mentioning that the Gibbs free energy difference $ΔG$ depends on the choice of the Hubbard U parameter of the DFT $+U$ method, and U is $5 eV$ here. However, We have checked that the trend of $ΔG$ as a function of pressure p does not change for $U=8 eV$⁠. We conjecture that this increasing relative instability leads to a transition from an initial metastable HT-AFM state to the more stable LT-AFM state at some intermediate pressure, so that after the sample is taken out of the pressure cell, the LT-FM state persists. Spin–orbit coupling has been taken into account in the calculation, and the result shows that the magnetizations of all atoms are in the out-of-plane direction. The magnetizations of several single atoms in the bottom CrI₃ layer, which also range over all the different atoms in the sense of symmetry, are shown in Figs. 18(d) and 18(e). Notice that in both the HT-AFM and LT-FM states, the magnetizations of I atoms are in the opposite direction from that of the Cr atoms in the same CrI₃ layer. Finally the magnitude of total magnetization per unit cell on each CrI₃ layer as a sum of the magnetizations of individual atoms is plotted in Fig. 18(f), from which we see that only the strength of the AFM order is enhanced by pressure, while there is essentially no change in that of the FM order.

So far, we have demonstrated that a first-principles description of the gate field effect provides physical insight that cannot be revealed by empirical methods. As an experimental tool, gating is a convenient knob for controlling the electronic, magnetic, and electron transport properties of two-dimensional materials. As we reduce the size of systems or devices to the nanometer scale and if quantum mechanical laws are governing the physical processes, it is inevitable to deal with system performance at the electron level and its response to gate fields in various configurations. Much of the scientific outcomes were discussed in previous publications; here, we add a few lines to conclude the study of bilayer CrI₃ and the heterogeneous systems between bilayer CrI₃ and graphene/BN. We showed that the local Coulomb interaction (the Hubbard U parameter specifically) increases by about 1 eV as the dimension of CrI₃ reduces from 3D to 2D, which stabilizes the AFM state of the high temperature stacking. Our calculations also show that the FM state of the low temperature stacking can be stabilized by increasing pressure, which agrees with experiments. Inserting bilayer CrI₃ between two graphene sheets prohibits a magnetic phase transition in electric fields up to at least $0.9 V/Å$ due to electrostatic shielding. Note that the current experimental limit is only about 0.1 $V/Å$⁠. If bilayer CrI₃ is placed between a h-BN sheet (the bottom) and a graphene sheet (the top), a magnetic phase transition can be driven by an electric field but the field strength required for the phase transition is larger than that for bare bilayer CrI₃. According to our simulations, charge doping does not induce a magnetic phase transition, at least for doping concentrations below $1×1013 cm−2$⁠. For the h-BN|2-CrI₃|h-BN system, the top BN layer becomes hole doped at $∼0.4 V/Å$⁠, which inhibits the magnetic phase transition. The magnetic moment of the bare bilayer in the AFM state does not increase until the band gap closes at ∼0.3 $V/Å$⁠. In contrast, the magnetic moment of h-BN|2-CrI₃|graphene in the AFM state increases at a smaller electric field due to charge transfer between graphene and bilayer CrI₃.

Before the end of the paper, we present the outlook for ongoing and near future projects. The tunability of the interlayer magnetic order by a magnetic field permits giant tunneling magnetoresistance.⁸⁴ One of our ongoing projects is exploration of gate field effects on spin-dependent electron tunneling properties through bilayer CrI₃ in a graphene|bilayer CrI₃|graphene vertical tunneling junction. Inspired by experimental activities in the Center for Molecular Magnetic Quantum Materials (M²QM), we plan to study adsorption of Mn₁₂ molecules and other single-molecule magnets onto monolayer MoS₂, which has not yet been reported. It was shown that charging is an effective way to tune the magnetic anisotropy of a magnetic molecule;⁹⁹ and we have previously simulated the electron transfer between graphene and Mn₁₂ molecules without a gate field.¹⁰⁰ Following this thread, we started investigations to address how charge doping affects the magnetic anisotropy of single-molecule magnets such as Mn₁₂. A glance in the field of first-principles transport studies shows there is very little theoretical work^101,102 addressing the role of spin–orbit coupling of molecules or 2D junction under finite gate or bias voltages. One of the reasons for such studies being scarce is the lack of computational tools. Another one is the computational cost. At the frontier of methodology development in our group, a few things are on the horizon: first-principles modeling and algorithm development for computing the Schottky barrier, implementing spin–orbital coupling in ESM and transport calculations, and adding phonon–electron coupling in non-equilibrium Green's function and the ESM framework. In particular, the capability of including the spin–orbital coupling allows us to look at Janus monolayer MoSSe, which was recently synthesized in experiments.^103,104 An intrinsic out-of-plane potential buildup exists in MoSSe since S and Se atoms have different electron affinity. Our plan of making shared tools for the community will be based on the packages TranSIESTA and QuantumEspresso. Looking forward, we will then test and apply these tools to study in depth the roles of spin–orbit coupling, gate effects on magnetoelectric coupling and magnetorestriction, spin-phonon coupling, and manifestations of these couplings and effects in transport measurements.

This work was supported by the U.S. Department of Energy (DOE), Office of Basic Energy Sciences (BES), under Contract No. DE-FG02–02ER45995. Computations were done using the utilities of the National Energy Research Scientific Computing Center and University of Florida Research Computing.

The data that support the findings of this study are available within the article.