Tuning the critical Li intercalation concentrations for MoX₂ bilayer phase transitions

We calculate total energies as well as the electronic density of states (DOS) of Li ion-intercalated bilayer TMD structures within the framework of ab initio DFT using the VASP code [16]. We use projector augmented wave (PAW) pseudo-potentials [17] and treat the exchange and correlation terms within the generalized gradient approximation with the PBE parameterization [18]. The energy cutoff was set to 400 eV.

We consider bilayer structures where the two monolayers are rotated by 180^∘ with respect to each other. The translational alignment between the two MoX₂ monolayers is such that the Mo atoms are vertically aligned (refer to Fig. 1), since within our DFT parameter settings we find the energy for this configuration to be lower than for other alignments in the presence of Li-ion intercalation. In particular, this translational alignment yields for MoS₂ a total energy (per a rectangular $4\times 2$ supercell that contains $1$ intercalated Li ion) $22$ meV lower than the corresponding lowest energy of the intercalated AA’ stacking (where the top layer’s Mo atoms are aligned with the bottom layer’s S atoms and viceversa [19, 12]).

We construct systems using supercells obtained by laterally replicating the rectangular unit cell of bilayer MoX₂ with X=S, Se or Te. We use $4\times 2$ lateral replicas with lateral area of $\approx$ 1.4$\times$1.2 nm² and vacuum separation along the out-of-plane direction of more than $15$ Å. There are 32 Mo atoms and 64 X atoms in our supercells in the absence of ion intercalation, and $N_{\text{sites}}=16$ lithium atoms at saturation ($\phi=1$) of the low energy sites.

Even with only tens of sites, there is a large configurational space associated with Li ion intercalation in between the MoX₂ monolayers of the bilayer. To explore the configurational space in order to compute ensemble averages, it is sufficient to sample only stable configurations as they carry most of the weight in Eq. (3) given that they correspond to local energy minima. For each phase we identify low energy ion-intercalation sites and select the most stable one, as we find that this crystallographic location has an energy several hundred meV lower than the other, metastable configurations. In the case of the H phase we find that the Li ion is most stable when positioned half-way in between two vertically aligned Mo atoms from different monolayers (see Fig. 2a)) [12]. This results in octahedral coordination to the X atoms (with bond length of about $2.8$ Å). In the T’ case the most stable Li site is not vertically aligned with Mo while coordination with the closest six X neighbors is distorted from the ideal octahedral case with four different bond lengths ranging from $2.72$ Å  to $3.06$ Å (see Fig. 2b). We also checked in the MoS₂ case for both phases that Li ion adsorption yields an energy significantly larger (by more than $100$ meV) than Li ion intercalation, hence the focus of this work is on ion intercalation.

The Li ion-intercalation sites in the 4$\times$2 supercell form a 16-site layer sublattice, i.e.​  there is one possible Li intercalation site for each Mo₂X₄ bilayer formula unit. The Li-concentration $x$ in the stoichiometric formula Li_xMoX₂ ranges from $x=0$ (no intercalation, $N_{Li}=0$, $\phi=0$) to $x=0.5$ at full ($N_{\text{Li}}=16$, $\phi=1$) Li sublattice occupancy. The number of possible configurations for fixed $N_{\text{Li}}$ is $N_{\text{sites}}!/(N_{\text{Li}}!(N_{\text{sites}}-N_{\text{Li}})!)$. Summing over $N_{\text{Li}}$ from $0$ to $16$ leads to 65,536 total configurations per system/phase. We do not consider mixtures of phases (domains) in a monolayer nor the possibility that two monolayers could belong to a different phase.

For each of the three MoX₂ systems and for both H and T’ phases we randomly sample the configurational space such that the configurations with $N_{\text{Li}}$ intercalated Li spanning the range from $0$ to $N_{sites}=16$ are roughly equally represented. Each configuration is optimized by minimizing the stress tensor and relaxing the atomic forces to less than $0.005$ eV/ Å. During optimization the Brillouin zone is sampled with a 4$\times$4$\times$1 k-point grid. Convergence in the electronic density of states of the relaxed structures is achieved with a 10$\times$10$\times$1 k-point grid. We do not include spin-orbit coupling due to the associated computational cost and because it is not expected to affect significantly the structural relaxation; we checked that spin-orbit coupling has an insignificant impact on the density of states near the Fermi level.

Obtaining the total energy of all possible distinct initial configurations is computationally too expensive with ab initio DFT even for the relatively small 4$\times$2 supercell. To this end we explore several surrogate models that can sample exhaustively the entire configurational space given DFT training data.

The first method we consider is cluster expansion (CE) [20, 21, 22], a classical technique based on a linear coefficient model with a nearest neighbor expansion

where the features $f_{a}$ are formed from a polynomial expansion of binned nearest neighbors $n_{b}$. Here $n_{b}$ is the total number of neighbors in the $b$-th bin. The bins we consider are annular rings with width small enough such that they contain one type of neighbor. In a $n$-th order expansion $f_{a}$ are complete with powers of $n_{a}$ up to $n_{CE}$. The coefficients $c_{a}$ are found via regression to training data. In the standard form of CE, the polynomial order is set to $n_{CE}=1$; we explored higher order polynomials as well.

We also used several machine techniques such as a basic multilayer perceptron (MLP) neural network (NN) as well as the more complex crystal graph convolutional neural network (CGCNN) [23] and image-based transfer learning with a convolutional neural network [24]. The MLP model is a feed-forward model of densely connected layers. The output of each preceding layer (starting with the selected inputs) is transformed by a linear weight matrix and then a non-linear activation function. The coefficients of the linear transformation of each layer comprise the parameters of the model. We trained the parameters of the MLP using the neighbor counts $n_{b}$ as inputs and the corresponding energies as outputs. The CGCNN formalism is more complex and is based on a universal and interpretable representation of crystalline materials. An input crystallographic structure $\chi_{i}$ is encoded as a graph whose nodes, corresponding to atoms, are connected by edges, corresponding to bonds, with surrounding atoms. The CGCNN architecture then applies a nonlinear graph convolution function that updates the atoms’ feature vectors with information from the surrounding atoms and bonds. A pooling layer then produces an overall feature vector for the entire crystal that satisfies invariance with respect to atom indexing (permutational invariance) and number of atoms in the unit cell (size invariance). Similar to the MLP training, we optimized the parameters of the CGCNN to fit the DFT predicted energy. We also applied a transfer learning technique using a pre-trained convolutional neural network (CNN), googlenet [25], originally used for image recognition. We replaced the classification last layer of this 22-layer deep CNN with a new regression layer and re-trained the network on predicting energies based on input images of the non-self-consistent electrostatic potential projected on a plane that passed through the unrelaxed Li-ion layer. The non-self-consistent electrostatic potential was selected since it is inexpensive to compute and the projection was necessary since googlenet operates on two-dimensional images.

For all these approaches the training set consisted of hundreds to thousands of input-output pairs. The inputs correspond to unrelaxed Li-intercalated MoX₂ structures. Prior to the training models of the energy as a function of configuration, we transformed total energies $E$ obtained from DFT into formation energies $\tilde{E}$ according to a linear transformation Eq. (7). As seen in Fig. 3, this transformation reduces the energy range of the outputs by more than an order of magnitude, from about $50$ eV (the total energy decreases with increasing Li-site occupancy fraction because the Li intercalation energy is positive [12]) to about $3$ eV, facilitating better convergence during the training process.

To make predictions regarding the electronic conductivity of Li ion-intercalated bilayer TMDs as function of ion concentration we consider the electronic density of states (DOS) near the Fermi level (E_F) for each system/phase. Indeed, within band theory, DOS(E_F) is directly proportional to the low voltage bias electronic conductivity.

Fig. 7 and Fig. 8(a) show the DOS for each system/phase in the absence of Li ion intercalation. We observe that the H is semiconducting while the T’ phase is semimetallic. In the H case the band gap decreases with increasing mass of the chalcogen atom X.

To see the impact of Li ion intercalation on the DOS we consider in Fig. 8(b) the case of MoSe₂ with one Li ion/supercell intercalated in the bilayer. We note that for this doping level the DOS of both phases do not renormalize appreciably. Instead, they can be obtained from the no-intercalation case by a simple shift-up of the Fermi level. Due to the semiconducting band gap in the H case, the shift to the left of the DOS is correspondingly large $>1$ eV. In the semi-metallic T’ case the DOS shift is much smaller. Thus, the change in the DOS(E_F) with low Li ion concentration is expected to be small in the T’ phase case and large for the H phase.

With larger Li ion concentration it is possible that the DOS renormalizes, i.e. cannot be obtained from the no-intercalation ($N_{\text{Li}}=0$) case by a simple shift. This can be seen in Fig. 9(a) which shows the DOS of bilayer MoS₂ (T’ phase) for several intercalation levels. As just mentioned, at low ion concentration the DOS can be obtained from the no-intercalation case by a simple, small shift (compare the black $\phi=0$ and red $\phi=1/16$ lines). However, in the case of half-occupancy of the Li sublattice ($\phi=1/2$, $x=1/4$), the DOS shows very low DOS(E_F), clearly a renormalization feature that is not evident in the absence of Li. At full occupancy of the Li sublattice ($\phi=1$, $x$=1/2), the renormalization is so strong that it induces a small gap in the bandstructure about the Fermi level. The strong renormalization of the electronic structure that takes place at higher ion concentration for all three systems (especially for the T’ phase) is not due to a simple electronic doping effect, but rather by geometrical relaxation/disordering of the bilayer MoX₂ atoms induced by the intercalated Li ions. The disordering effect takes place even at full occupancy of the Li sublattice, where prior to atomic optimization the structure has translational symmetry. This can be seen in Fig. 9(b) for the case of MoS₂ T’ phase at full Li intercalation ($\phi=1$) showing that an opening of the band gap at the Fermi level occurs only upon full atomic relaxation that breaks the translational symmetry and renders the Li ions inequivalent. This disordering effect is subtle, with the mean distortion of atomic positions of $<0.015$ Å. We note that periodic lattice distortions in TMDs have also been connected to a charge-density-wave formation mechanism [26, 27, 26].

The evolution of DOS(E_F) with ion concentration is obtained within the NT ensemble by statistically averaging the DFT calculated DOS(E_F) for each of the $\approx 300$ configurations per bilayer system/phase. Here we do not use a surrogate model to estimate DOS(E_F) in the entire configurational space, but in Appendix B we show how two surrogate models (CE and NN) compare with DFT when representing the DOS.

Fig. 10 illustrates our procedure for obtaining the statistically averaged DOS(E_F) as function of ion concentration from the raw DFT data for the case of the bilayer MoS₂. We note the peculiar feature of very small DOS(E_F) for the T’ phase at half Li sublattice occupancy. It is known that one-dimensional metals are unstable against a distortion of the crystal lattice that lowers the system energy and opens an electronic band gap, the so-called Peierls gap [28]. A similar transformation may occur in low-dimensional systems such as two-dimensional TMDs [29]. We suspect that the band gap opening that we observe at certain Li-ion site occupation fraction in the T’ phase of MoS₂ and MoSe₂ is similar in nature to a Peierls gap opening.

Fig. 11 summarizes our results for the evolution of DOS(E_F) for each of the bilayer system/phase. For MoSe₂ and MoTe₂ systems, DOS(E_F) increases more rapidly with increasing Li concentration for the H phase than for the T’. At relatively low ion concentration, the H phase may be therefore more conductive despite the fact that in the no-intercalation case it is semiconducting. In general, as the concentration increases, the H phase remains the most conductive, with the peculiar situation in MoS₂ case where the T’ phase seems to show very low conductivity near half Li sublattice occupancy ($x\approx 0.25$). Importantly, as opposed to the expectation that a change in phase structure from H to T’ correlates with a (positive) jump in electronic conductivity, we do not find this to hold for any of the three bilayer TMDs. In fact our results suggest that as the system transitions between the two phases, the conductivity suffers a decrease instead of a (positive) jump.