Exchange Intervalley Scattering And Magnetic Phase Diagram Of Transition Metal Dichalcogenide Monolayers

TL;DR

This paper investigates the magnetic phase transitions in transition metal dichalcogenide monolayers, revealing a first order ferromagnetic transition influenced by exchange intervalley scattering and spin-orbit coupling, consistent with recent experiments.

Contribution

It introduces a theoretical model predicting a first order ferromagnetic phase transition in monolayer TMDs due to exchange intervalley scattering and spin-orbit effects.

Findings

Prediction of a first order ferromagnetic phase transition.

Identification of non-analytic cubic terms in free energy.

Consistency with recent experimental observations on MoS₂ monolayers.

Abstract

We analyze magnetic phases of monolayers of transition metal dichalcogenides that are two-valley materials with electron-electron interactions. The exchange inter-valley scattering makes two-valley systems less stable to the spin fluctuations but more stable to the valley fluctuations. We predict a first order ferromagnetic phase transition governed by the non-analytic and negative cubic term in the free energy that results in a large spontaneous spin magnetization. Finite spin-orbit interaction leads to the out-of-plane Ising order of the ferromagnetic phase. Our theoretical prediction is consistent with the recent experiment on electron-doped monolayers of MoS$_2$ reported by Roch $\textit{et al.}$ [1]. The proposed first order phase transition can also be tested by measuring the linear magnetic field dependence of the spin susceptibility in the paramagnetic phase which is a direct consequence of the non-analyticity of the free energy.

Equations

$$H=\frac{k^2}{2m}-\alpha {\sigma }_z\tau ,$$

$$v=\nu v_1(0)=\nu v_2(0)>u=\nu v_3(0),$$

$$M_s=\sum _{\sigma ,\tau }\frac{\sigma n_{\tau }^{\sigma }}{4\nu },M_v=\sum _{\sigma ,\tau }\frac{\tau n_{\tau }^{\sigma }}{4\nu },M_{\alpha }=\sum _{\sigma ,\tau }\frac{\sigma \tau n_{\tau }^{\sigma }}{4\nu },$$

$$E_F=\frac{1}{4\nu }\sum _{\sigma ,\tau }{n}_{\tau }^{\sigma }=\frac{n}{4\nu }=const.,$$

$$F=\Omega +\sum _{\sigma ,\tau }{\mu }_{\tau }^{\sigma }{n}_{\tau }^{\sigma }=\Omega -\sum _{\sigma ,\tau }{\mu }_{\tau }^{\sigma }\frac{\partial \Omega }{\partial {\mu }_{\tau }^{\sigma }}.$$

$${\Omega }_0=-\frac{\nu }{2}\sum _{\sigma ,\tau }{\left({\mu }_{\tau }^{\sigma }+\alpha \sigma \tau \right)}^2.$$

$$F_0=2\nu \left(M_s^2+M_v^2+M_{\alpha }^2-2\alpha M_{\alpha }\right).$$

$$F_1=-2\nu \left[v_{+}{M}_s^2+v_{-}\left(M_v^2+M_{\alpha }^2-2\alpha M_{\alpha }\right)\right],$$

$$\delta {\chi }_{zz}\sim {\chi }_0\frac{|{\mu }_B{B}_z|}{E_F},$$

$$\delta F=-\frac{2\nu }{3E_F}\{(v-u{)}^{2}f(M_v,M_{\alpha })+u^2\left(|M_s{|}^3+|M_v{|}^3\right)$$

$$+v^2\left[f(M_s,M_{\alpha })+f(M_s,M_v)\right]\},$$

$$g_{\tau }^{\sigma }(K)=\frac{1}{i\omega -\epsilon (k)+\alpha \sigma \tau +{\mu }_{\tau }^{\sigma }},$$

$${\tilde{\mu }}_{\tau }^{\sigma }={\mu }_{\tau }^{\sigma }+\alpha \sigma \tau .$$

$${\Omega }_0=\sum _{\sigma ,\tau }\sum _K\ln \left[g_{\tau }^{\sigma }(K)\right],$$

$${\Omega }_0=-\frac{\nu }{2}\sum _{\sigma ,\tau }({\tilde{\mu }}_{\tau }^{\sigma }{)}^2,$$

$$F_0=2\nu \left(M_s^2+M_v^2+M_{\alpha }^2-2\alpha M_{\alpha }\right),$$

$${\Omega }_1=-\frac{1}{2\nu }\sum _{\sigma ,\tau ;K,P}\left[v{g}_{\tau }^{\sigma }(K)g_{\tau }^{\sigma }(P)+u{g}_{\tau }^{\sigma }(K)g_{-\tau }^{\sigma }(P)\right],$$

$$\sum _K{g}_{\tau }^{\sigma }(K)=\nu {\tilde{\mu }}_{\tau }^{\sigma }.$$

$${\Omega }_1=-\frac{\nu }{2}\sum _{\sigma ,\tau }\left[v({\tilde{\mu }}_{\tau }^{\sigma }{)}^2+u{\tilde{\mu }}_{\tau }^{\sigma }{\tilde{\mu }}_{-\tau }^{\sigma }\right].$$

$$M_s=(1+v_{+})\sum _{\sigma ,\tau }\frac{\sigma }{4}{\tilde{\mu }}_{\tau }^{\sigma },M_v=(1+v_{-})\sum _{\sigma ,\tau }\frac{\tau }{4}{\tilde{\mu }}_{\tau }^{\sigma },$$

$$M_{\alpha }=(1+v_{-})\sum _{\sigma ,\tau }\frac{\sigma \tau }{4}{\tilde{\mu }}_{\tau }^{\sigma },$$

$$F_1=-2\nu \left[v_{+}{M}_s^2+v_{-}\left(M_v^2+M_{\alpha }^2-2\alpha M_{\alpha }\right)\right].$$

$$F_0+F_1=a_s{M}_s^2+a_v{M}_v^2+a_{\alpha }(M_{\alpha }^2-2\alpha M_{\alpha }),$$

$${\Omega }_a$$

$${\Omega }_b$$

$${\Omega }_c$$

$${\Omega }_d$$

$${\Pi }_{\tau {\tau }^{\prime }}^{\sigma {\sigma }^{\prime }}(Q)=\frac{1}{\nu }\sum _K{g}_{\tau }^{\sigma }(K)g_{{\tau }^{\prime }}^{{\sigma }^{\prime }}(K+Q)\equiv {\mathcal{P}}_{\tau {\tau }^{\prime }}^{\sigma {\sigma }^{\prime }}(Q)-1,$$

$${\Omega }_a=-\frac{v^2}{4}\sum _{{\tau }_1,{\tau }_2}\sum _{{\sigma }_1,{\sigma }_2}\sum _Q{\mathcal{P}}_{{\tau }_1{\tau }_1}^{{\sigma }_1{\sigma }_1}(Q){\mathcal{P}}_{{\tau }_2{\tau }_2}^{{\sigma }_2{\sigma }_2}(Q),$$

$${\Omega }_b=\frac{v^2}{4}\sum _{\tau ,\sigma }\sum _Q[{\mathcal{P}}_{\tau \tau }^{\sigma \sigma }(Q){]}^2,$$

$${\Omega }_c=-\frac{u^2}{4}\sum _{\tau ,{\sigma }_1,{\sigma }_2}\sum _Q{\mathcal{P}}_{\tau -\tau }^{{\sigma }_1{\sigma }_1}(Q){\mathcal{P}}_{\tau -\tau }^{{\sigma }_2{\sigma }_2}(Q),$$

$${\Omega }_d=\frac{uv}{2}\sum _{\tau ,\sigma }\sum _Q[{\mathcal{P}}_{\tau -\tau }^{\sigma \sigma }(Q){]}^2.$$

$${\mathcal{P}}_{\tau {\tau }^{\prime }}^{\sigma {\sigma }^{\prime }}(\mathbf{q},\epsilon )\approx \frac{|\epsilon |}{\sqrt{{\left(\epsilon -i{\Delta }_{\tau {\tau }^{\prime }}^{\sigma {\sigma }^{\prime }}\right)}^2+{\left(v_{F}q\right)}^2}},$$

$${\Delta }_{\tau {\tau }^{\prime }}^{\sigma {\sigma }^{\prime }}={\tilde{\mu }}_{\tau }^{\sigma }-{\tilde{\mu }}_{{\tau }^{\prime }}^{{\sigma }^{\prime }}={\mu }_{\tau }^{\sigma }-{\mu }_{{\tau }^{\prime }}^{{\sigma }^{\prime }}+\alpha (\sigma \tau -{\sigma }^{\prime }{\tau }^{\prime }).$$

Full text

I Introduction

Multi-valley materials provide an alternative to spintronics where, instead of spin projections, multiple valleys are used as pure states for quantum bits. This additional valley degree of freedom also provides the ground for complex phase diagrams for multi-valley materials kelly ; bloss . Transition metal dichalcogenides (TMDs) are particularly promising two-valley materials for valleytronics wang due to the large spin-orbit interaction (SOI) zhu ; xiao ; chei ; kormanyos ; kormanyos2013 ; klinovaja ; review_1 ; pisoni and the control of the valley polarization by external magnetic fields kormanyos ; cheng ; aivazian . Large SOI also results in the valley-dependent light absorption which yields another way to control the valley polarization at low densities zeng ; mak .

Effects of electron-electron interactions in TMDs were considered previously in Refs. donck ; mukherjee ; braz The results predicted in Refs. donck ; mukherjee are valid if the exchange inter-valley (EIV) interaction is negligible. The limit considered in Ref. braz corresponds to a SOI strength which is much larger than the Fermi energy. In this limit the EIV scattering is forbidden due to the spin-valley locking. Signatures of the valley ferromagnetic phase predicted in Ref. braz were observed in WS2 scrace where the SOI gap is relatively large $\sim 40\,$meV.

In this paper we consider the opposite limit when the SOI gap is small compared to the Fermi energy. This scenario can be realized in the conductance band of MoS2 where the SOI gap is on the order of a few meV. We predict a first-order ferromagnetic phase transition governed by non-analytic terms in the free energy. Such terms for one-valley materials have been calculated in Refs. maslovsaha ; maslovchubukov They originate from the infrared electron-hole fluctuations near the Fermi surface in two-dimensional (2D) interacting systems chubukov . The one-valley calculations maslovsaha ; maslovchubukov ; chubukov ; zakmaslov ; zakQ , however, cannot be applied to two-valley materials if the EIV scattering must be taken into account. Here we argue that the EIV scattering is important as it favors a ferromagnetic instability over valley and spin-valley instabilities. We also argue that the SOI breaks the $O(3)$ symmetry of the ferromagnetic phase leading to the Ising order. Similar scenario is realized in one-valley materials with the Rashba SOI zakmaslov ; zakQ . This work is motivated by the recent experiment in $n$-doped monolayer of MoS2 roch where the ferromagnetic Ising order has been detected. The non-analytic cubic term in the free energy can be directly probed as it results in an unusual linear magnetic field dependence of the spin susceptibility in the paramagnetic phase.

Our paper is organized as follows. In Section II we introduce the matrix elements of the electron-electron interaction. In Section III we derive that the EIV interaction makes possible only ferromagnetic and paramagnetic phases, other phases cannot develop due to the stabilizing effect of the EIV interaction on the valley and spin-valley fluctuations. In Section IV we outline the Ginzburg-Landau theory with the negative cubic non-analyticity in case of a single order parameter. In Section V we present the non-analytic cubic correction to the free energy of two-valley 2D electron gas. Conclusions are given in Section VI. The diagrammatic calculations and other technical details are deferred to Appendices A-C.

II Electron-electron interaction

Monolayers of TMDs are two-valley materials. The electron spectrum of such materials can be effectively described by the following free-particle Hamiltonian kormanyos :

[TABLE]

where $\bm{k}=(k_{x},k_{y})$ is the 2D momentum, $m$ is the effective mass, $2\alpha$ is the SOI gap, $\sigma_{z}$ is the spin Pauli matrix, the index $\tau=\pm 1$ labels the valley. The spectrum of $H$ consists of four bands, see Fig. 1. Each pair of bands is split by the SOI $\alpha$ which also lifts the spin degeneracy and, hence, defines the spin quantization axis.

The presence of two valleys results in three different scattering channels for interacting electrons: the intra-valley interaction $v_{1}(q)$ where interacting electrons are from the same valley, Fig. 2(a); the direct inter-valley interaction $v_{2}(q)$ where electrons from different valleys scatter without changing the valley, Fig. 2(b); the EIV interaction $v_{3}(q)$ where electrons exchange their valleys, Fig. 2(c). Matrix elements $v_{1,2}(q)$ are both described by the screened Coulomb interaction, see e.g. Ref. maslovchubukov , $v_{1}(q)=v_{2}(q)=2\pi e^{2}/[\epsilon(q+k_{TF})]$, where $e$ is the electron charge, $\epsilon$ is the dielectric constant, $k_{TF}=4/a_{B}$ is the Thomas-Fermi screening wave vector, $a_{B}=\epsilon/me^{2}$ is the effective Bohr radius. In all relevant processes the momentum transfer $q\lesssim k_{F}$, where $k_{F}=\sqrt{\pi n}$ is the Fermi momentum, $n$ is the total electron density. At small to moderate densities $k_{F}a_{B}\ll 1$, so the term $qa_{B}$ can be neglected. Neglecting the regular $q$-dependence of the interaction matrix elements does not affect the physics considered here and is done for simplicity.

The EIV matrix element $v_{3}(q)$, Fig. 2(c), is short-range due to the large momentum transfer $2q_{0}+q\approx 2q_{0}\gg k_{F}$, where $2q_{0}$ is the distance between two valleys in the first Brillouin zone. Thus, $v_{3}(q)\approx v_{3}(0)$. One can estimate $v_{3}(0)$ by the bare Coulomb interaction: $v_{3}(0)\sim V_{C}(2q_{0})=\pi e^{2}/\epsilon q_{0}$. However, this estimate is rather crude as it does not take into account the overlap of Bloch amplitudes corresponding to different valleys and the on-site Hubbard interaction, see Ref. roldan As in TMDs $a_{B}\sim 1\,$nm $\sim a_{0}$, with $a_{0}$ being the lattice constant, then $v_{3}(0)\sim v_{1,2}(0)$. For further analysis we introduce the dimensionless interaction matrix elements:

[TABLE]

where $\nu=m/2\pi$ is the 2D density of states per parabolic band. In this work we only consider the case $v>u$. For the opposite case, $v<u$, the system develops a Cooper instability leading to superconducting phases footnote1 .

III Asymmetry between order parameters induced by the EIV scattering

In this section we demonstrate that the EIV interaction narrows down possible magnetic phases to the ferromagnetic and the paramagnetic ones. This is due to stabilizing effect of the EIV interaction on the valley and spin-valley fluctuations. In order to examine magnetic instabilities, we lift the assumption that the chemical potential is identical for all bands. Instead, we introduce four different chemical potentials $\mu^{\sigma}_{\tau}$, $\tau=\pm 1$, $\sigma=\pm 1=\uparrow,\downarrow$, see Fig. 1. The SOI defines the quantization axis for the electron spin and thus leads to the following Ising order parameters:

[TABLE]

where $n^{\sigma}_{\tau}=-\partial\Omega/\partial\mu^{\sigma}_{\tau}$ is the electron density in the band with spin and valley indexes $\sigma$ and $\tau$, $\Omega$ is the grand canonical potential per unit area that depends on $\mu^{\sigma}_{\tau}$ and $\alpha$. The introduced order parameters are independent and have dimension of energy. In this description we keep the total electron density $n=\sum_{\sigma,\tau}n^{\sigma}_{\tau}$ fixed by the condition:

[TABLE]

where $E_{F}$ corresponds to the Fermi energy in the non-interacting case. We emphasize that $E_{F}$ is not a chemical potential in the system.

In order to find which realizations of the chemical potentials $\mu^{\sigma}_{\tau}$ correspond to the ground state of the interacting Fermi gas, we calculate the free energy $F$ per unit area as a function of the magnetizations $M_{i}$. As the Feynman diagram technique is designed for the grand canonical potential $\Omega$ (per unit area), we first calculate interaction corrections to $\Omega$ and then apply the Legendre transform to find the free energy $F$:

[TABLE]

The non-interacting part of the grand canonical potential per unit area is the following, see Appendix A:

[TABLE]

The free energy $F_{0}$ per unit area of the non-interacting electron system is then given by

[TABLE]

Here we have suppressed the constant term $2\nu E_{F}^{2}$. The term $-4\nu\alpha M_{\alpha}$ represents the coupling of the magnetization $M_{\alpha}$ to the SOI $\alpha$ which plays the role of an external field. The EIV interaction leads to an asymmetry between the order parameters already in first order perturbation theory, see Appendix A:

[TABLE]

where $v_{\pm}=v\pm u$. This asymmetry turns out to be very important for the magnetic phase diagram.

In order to illustrate why the asymmetry between different order parameters is important, we first remind the reader of the Ginzburg-Landau (GL) theory that describes second-order phase transitions. This theory assumes that the free energy is analytic with respect to the order parameters. For example, in case of one order parameter $M$ the GL free energy $F$ reads: $F(M)=aM^{2}+bM^{4}+\dots$ The parameters $a$, $b$ depend on the interaction, $b>0$ assures a minimum of $F(M)$ at finite $M$. If the interaction is weak, then $a>0$ and $F(M)$ has its minimum at $M=0$ meaning that there is no net magnetization in the ground state. If the interaction is strong such that $a<0$, then $F(M)$ has its minimum at $|M|=\sqrt{-a/2b}\neq 0$ which manifests a ferromagnetic ground state, see Fig. 3. The phase transition occurs at $a=0$. As there is no magnetization jump at the transition point, this is a second-order phase transition. In our case, $a$ is different for different order parameters: $a_{s}=2\nu(1-v_{+})$ for the spin order parameter $M_{s}$ and $a_{v}=a_{\alpha}=2\nu(1-v_{-})$ for the valley and spin-valley order parameters, see Eqs. (7) and (8). As $v_{+}>v_{-}$ and thus $a_{s}<a_{v}=a_{\alpha}$, the GL theory predicts a second-order ferromagnetic phase transition which occurs at $v_{+}=1$ corresponding to $a_{s}=0$. Therefore, the critical condition is first satisfied for the ferromagnetic phase. Phases other than ferromagnetic or paramagnetic are not possible due to the energy cost coming from the EIV scattering.

In this paper we make use of the language of magnetizations, see Fig. 3. However, the GL theory can be alternatively understood in terms of the linear response (susceptibility) in the paramagnetic phase. In the above example, the susceptibility $\chi\sim 1/\frac{\partial^{2}F(M)}{\partial M^{2}}|_{M=0}\sim 1/a$ represents the Stoner pole for $a\to 0$ and is only well-defined in the paramagnetic phase where $a>0$. At the phase transition $a=a_{cr}=0$ the susceptibility diverges which is the standard signature of a second-order phase transition.

IV First order phase transition due to the non-analytic terms in the free energy

In this section we outline the Ginzburg-Landau theory with the negative cubic non-analyticity in case of a single order parameter belitz . The assumption of analyticity of $F(M)$, which underlies the standard GL theory, is not generally true belitz . The non-analytic corrections coming from the infrared electron-hole fluctuations are especially strong in 2D electron systems where the leading non-analyticity is cubic in $M$: $F(M)=aM^{2}+bM^{4}+c|M|^{3}+\dots$, see Refs. maslovsaha ; maslovchubukov ; belitz The first two terms in $F(M)$ are exactly the same as in the GL theory. If $c>0$, then there is no qualitative difference from the GL theory and in this case the theory describes the second-order phase transition at $a_{cr}=0$. However, if $c<0$, then the second minimum appears even at $a>0$. This second minimum becomes the global one at $a<a_{cr}=c^{2}/4b$ with a finite magnetization jump $|M_{cr}|=-c/2b$, see Fig. 3. Note that in this case $a_{cr}>0$, and the first order phase transition occurs before the Stoner instability of the paramagnetic minimum, corresponding to $a=0$, develops. In fact, at $a<a_{cr}$ the evolution of the paramagnetic point $M=0$ is completely irrelevant as the system is settled in the vicinity of the ferromagnetic minimum. In particular, the Stoner pole in the susceptibility of the paramagnetic phase, $\chi\sim 1/a$, can no longer be reached as the phase transition occurs at $a=a_{cr}>0$. In the ferromagnetic phase the Stoner formula for the susceptibility is no longer valid. The finite magnetization jump and finite susceptibility of the paramagnetic phase at the critical point demonstrate the first order nature of the phase transition described by the theory with negative cubic non-analyticity, see Fig. 3.

The negative cubic non-analyticity cannot be ignored as it defines the nature of magnetic phase transitions. In our case we have three independent order parameters and the previous theories that considered a single order parameter cannot be applied to our problem. In the next section we calculate the non-analytic cubic correction to the free energy in the leading order in the interaction.

The statement that the EIV scattering narrows down all non-trivial magnetic phases to the ferromagnetic one stays true even when the cubic non-analyticity is taken into account. Indeed, in the previous section we showed that the smallest $a$ corresponds to the spin order parameter $a_{s}$, so $a_{s}$ is always closer than $a_{v}=a_{\alpha}$ to the critical value $a_{cr}$ at which the first order phase transition occurs. Therefore, due to the EIV scattering the ferromagnetic instability always wins over the valley and spin-valley instabilities.

The considered mechanism of the first order phase transition relies on the negative cubic non-analyticity in the free energy. It can be independently probed in the paramagnetic phase where this non-analyticity results in an unusual linear magnetic field dependence of the spin susceptibility, see Appendix C:

[TABLE]

where $\chi_{0}>0$ is the spin susceptibility of the non-interacting 2DEG, $\mu_{B}$ is the Bohr magneton, and $B_{z}$ is the perpendicular magnetic field, and the proportionality coefficient is given in the appendix C. We emphasize that $\delta\chi_{zz}$ must be positive in order to produce the negative cubic term in the free energy. Similar non-analyticity of the susceptibility in the paramagnetic phase was predicted theoretically for one-valley 2D electron systems maslovsaha ; maslovchubukov ; chubukov ; zakmaslov ; zakQ . In our case the prefactor is different due to the valley contribution, see Appendix C.

V Non-analytic cubic correction to the electron free energy

Here we calculate the second order Feynman diagrams that contribute to the non-analytic terms in the free energy. The second order diagrams that appear in the RPA series are not included here as they are fully absorbed in the quadratic part of the free energy, see Appendix A. The diagrams in Fig. 4(a),(b) do not involve the EIV scattering and thus correspond to the one-valley calculations maslovchubukov ; zakmaslov . In particular, the non-analytic part of the diagram in Fig. 4(b) vanishes. Going beyond this, we now also consider contributions of the processes involving one or two EIV scattering events that are represented by diagrams in Fig. 4(c),(d). The diagram in Fig. 4(d) is especially interesting because its sign is opposite to the diagrams in Fig. 4(a),(c). We outline the calculations of these diagrams in Appendix B.

The cubic correction to $\delta\Omega$ given by the diagrams in Fig. 4 leads to corresponding correction to the free energy $\delta F$, see Eq. (5):

[TABLE]

where $f(x,y)=(|x+y|^{3}+|x-y|^{3})/2$. This is the central result of the paper. Contributions from different diagrams can be recognized by the interaction factors $v^{2}$, $u^{2}$, and $uv$ corresponding to the diagrams in Fig. 4(a), Fig. 4(c) and Fig. 4(d), respectively (the diagram in Fig. 4(b) does not contribute to the cubic non-analyticity). Notice that Eq. (10) contains mixed non-analytic terms like $|M_{s}\pm M_{v}|$ which makes it qualitatively different from the case of a single order parameter. Eq. (10) is valid under the condition that all order parameters $M_{i}$ are much smaller than the Fermi energy $E_{F}$ and corresponds to the limit of zero temperature $T=0$. The temperature effects become important at $T\gtrsim|M^{*}|$, where $M^{*}$ is the spontaneous magnetization of the ground state (from one of the order parameters). At this temperature regime the non-analyticity disappears as $|M|^{3}\to M^{4}/T$. Therefore, one can expect the critical temperature $T_{c}\sim|M^{*}|$.

In order to analyze Eq. (10), it is useful to consider three limiting cases when only one of the order parameters is non-zero. If $M_{v}=M_{\alpha}=0$, then $\delta F=-2\nu|M_{s}|^{3}(2v^{2}+u^{2})/3E_{F}$. For $u=0$, it reproduces results obtained for the one-valley case of Ref. maslovchubukov If $M_{s}=M_{\alpha}=0$, then $\delta F=-2\nu|M_{v}|^{3}(2v^{2}+2u^{2}-2uv)/3E_{F}$. If $M_{s}=M_{v}=0$, then $\delta F=-2\nu|M_{\alpha}|^{3}(2v^{2}+u^{2}-2uv)/3E_{F}$. We notice that the spin term has the largest prefactor $2v^{2}+u^{2}$. The difference in prefactors is again due to the EIV scattering. Here we point out that the asymmetry mostly arises from the diagram in Fig. 4(d) which reduces prefactors for the valley and spin-valley cubic terms. The sum of $\delta F$ and the quadratic part $F_{0}+F_{1}$ [see Eqs. (7) and (8)] is minimal for the spin term, this allows us to conclude once again that the ferromagnetic instability $M_{s}\neq 0$, $M_{v}=M_{\alpha}=0$ wins over other possible magnetic instabilities. We stress that the ferromagnetic instability is favored over others because of the EIV scattering. The phase diagram of monolayers of TMDs in absence of the EIV scattering is predicted to be much richer, see Ref. donck This is why in this work we argue that the EIV scattering is an important interaction channel that cannot be neglected.

The Ising ferromagnetic phase was recently observed experimentally in monolayers of MoS2 roch . The reported spin polarization is close to 100%, i.e. $|M_{s}|/E_{F}\approx 1$. Our theory relies on the expansion over small parameters $M_{i}/E_{F}$, where $M_{i}$ are magnetizations defined in Eq. (3), and cannot be extended to the case when $M_{i}/E_{F}\sim 1$. However, it still allows us to grasp the origin of the observed instability. In particular, we predict the first order phase transition that is in agreement with the large spin magnetization observed in Ref. roch Next, we argue that the ferromagnetic instability is favored due to the EIV scattering. In fact, the limit of total spin magnetization with no valley and spin-valley magnetization corresponds to the case when the EIV scattering contribution to the interaction is maximized. Thus, our theory provides us with the qualitative argument why other phases were not observed.

VI Conclusions

In conclusion, we predict the Ising ferromagnetic phase with large spontaneous magnetization in monolayers of TMDs in case of $E_{F}\gg\alpha$. The first-order nature of the transition comes from the special role of electron-electron interactions in 2D which gives rise to unusual cubic magnetization terms with a negative sign. The non-analyticity can be potentially probed by measuring the linear magnetic field dependence of the spin susceptibility in the paramagnetic phase, see Eq. (9). The Ising order is due to the finite SOI that breaks the $O(3)$ rotational symmetry and defines the quantization axis for spin. The phases with non-zero spontaneous valley polarization cannot develop for $E_{F}\gg\alpha$ due to the EIV scattering. While our theoretical prediction agrees with recent experiments on MoS2 monolayers based on exciton spectra roch , it might be interesting to look for the spontaneous spin polarization directly by e.g. magnetization measurements.

Acknowledgments. We would like to thank J. Roch and R. Warburton for very stimulating discussions about their experimental findings which motivated this work. We acknowledge support by the Swiss National Science Foundation and NCCR QSIT as well as the Georg H. Endress foundation. This project received funding from the European Union’s Horizon 2020 research and innovation program (ERC Starting Grant, grant agreement No 757725).

Appendix A QUADRATIC PART OF THE FREE ENERGY

The electron Matsubara Green function corresponding to the Hamiltonian $H$ Eq. (1) is given by

[TABLE]

where $K=(\bm{k},\omega)$ is the “relativistic” notation for the momentum $\bm{k}=(k_{x},k_{y})$ and the Matsubara frequency $\omega$, $\varepsilon(k)=k^{2}/2m$. It is convenient to introduce the notation:

[TABLE]

The non-interacting part of the grand canonical potential (per unit area) is given by the following expression:

[TABLE]

where $\sum\limits_{K}=T\sum\limits_{\omega}\int d\bm{k}/(2\pi)^{2}$. Evaluating the sum in Eq. (13), we get

[TABLE]

where $\nu=m/2\pi$ is the 2D density of states. Using the Legendre transform Eq. (5) and the definition of magnetic order parameters Eq. (3), we find the non-interacting part of the free energy (per unit area):

[TABLE]

where we suppressed the constant term $2\nu E_{F}^{2}$ which does not depend on the order parameters. The Fermi energy $E_{F}$ is defined in Eq. (4).

The first order interaction correction to the grand canonical potential is given by the diagrams in Fig. 5 and reads explicitly:

[TABLE]

where the first and second terms correspond to the diagrams in Fig. 5(a) and Fig. 5(b), respectively. Here we just have to calculate $\sum_{K}g^{\sigma}_{\tau}(K)$ and find

[TABLE]

Substituting it back to Eq. (16), we get

[TABLE]

In order to calculate the correction to the free energy, we have to take into account the renormalization of the order parameters due to $\Omega_{1}$:

[TABLE]

where $v_{\pm}=v\pm u$. The SOI strength $\alpha$ is also renormalized as $\tilde{\alpha}=(1+v_{-})\alpha$. In what follows we use the notation $\alpha$ for the renormalized SOI. Applying the Legendre transform, we get the first order interaction correction to the free energy (per unit area)

[TABLE]

Together with $F_{0}$ we get the quadratic part of free energy accounting for the first order interaction correction:

[TABLE]

where $a_{s}=2\nu(1-v_{+})$, $a_{v}=a_{\alpha}=2\nu(1-v_{-})$. The higher order diagrams of the random phase approximation (RPA) series [see Fig. 6] lead to the renormalization of order parameters. In Fig. 6, wavy lines represent both components $v$ and $u$ of the interaction for brevity. It turns out that Eq. (21) also holds in this case if parameters are replaced by renormalized ones.

Appendix B NON-ANALYTIC CORRECTION TO THE FREE ENERGY

Here we outline the calculation of second-order diagrams that first contribute to the non-analytic cubic correction to the grand canonical potential $\Omega$, see Fig. 4, and thus to the free energy $F$. The derivation of these diagrams follows the procedure introduced in Refs. maslovsaha ; maslovchubukov ; zakmaslov We do not include the second order diagram from the RPA series [see Fig. 6], as it already has been taken into account. The algebraic representation for all four diagrams in Fig. 4 is the following:

[TABLE]

where $\nu=m/2\pi$ is the 2D density of states, sums are taken over all spin and valley indexes and over all momenta $K$, $Q$, and $P$. Notice the additional factor of $2$ in $\Omega_{d}$ which comes from the fact that $\Omega_{d}$ is represented by two diagrams: one is shown in Fig. 4(d) and the other is different by swapping dashed and wavy interaction lines.

Next, we perform the summations over $K$ and $P$ merging the Green functions into a particle-hole bubble:

[TABLE]

where $\mathcal{P}^{\sigma\sigma^{\prime}}_{\tau\tau^{\prime}}(Q)$ is the so-called dynamic part of the particle-hole bubble maslovsaha ; maslovchubukov ; zakmaslov . In fact, we only have to account for the dynamic part of the particle-hole bubbles because the static part is taken into account via the Thomas-Fermi screening of the Coulomb interaction. With this, we get the following expressions for the diagrams:

[TABLE]

The non-analyticity comes from the vicinities of two special points of the particle-hole bubble $\mathcal{P}^{\sigma\sigma^{\prime}}_{\tau\tau^{\prime}}(Q)$, $Q=(\bm{q},\varepsilon)$: (i) the Landau damping at $q=0$, $\varepsilon=0$ and (ii) the Kohn anomaly at $q=2k_{F}$, $\varepsilon=0$, $k_{F}$ is the Fermi momentum. The Landau damping only contributes if at least one of the bubbles $\mathcal{P}^{\sigma\sigma^{\prime}}_{\tau\tau^{\prime}}(Q)$ does not have both $\sigma=\sigma^{\prime}$ and $\tau=\tau^{\prime}$, see Refs. maslovchubukov ; zakmaslov Therefore, the Landau damping gives contributions to $\Omega_{c}$ and $\Omega_{d}$ only. For the Kohn anomaly only the difference in chemical potentials of two particle-hole bubbles is important maslovchubukov . Therefore, the Kohn anomaly does not contribute to $\Omega_{b}$ and $\Omega_{d}$ because the two particle-hole bubbles are identical there. This simple analysis readily shows that neither Landau damping nor Kohn anomaly contribute to $\Omega_{b}$, meaning that $\Omega_{b}$ is analytic.

First, we calculate the Landau damping contributions using the following approximation of the particle-hole bubble at small $Q=(\bm{q},\varepsilon)$ maslovsaha ; maslovchubukov :

[TABLE]

where $v_{F}=k_{F}/m$ is the Fermi velocity and

[TABLE]

Integration over $q$ can be done with the help of the following identity:

[TABLE]

where we have regularized the ultraviolet divergence by subtracting the integral at $\Delta_{1}=\Delta_{2}=0$. Therefore, the Landau damping contribution is given by

[TABLE]

where the index $L$ of the sum indicates that $Q$ is in the vicinity of the Landau damping point, $\mathcal{F}(z)$ represents the sum over $\varepsilon=2\pi nT$, with $T$ being the temperature (Boltzman constant $k_{B}=1$) and $n$ being an integer:

[TABLE]

This sum is ultraviolet divergent. However, we only need the essential dependence on $z$, meaning that one can ignore the divergent part by rewriting the sum in terms of a contour integral. The integration over the logarithm branch cut readily yields the following representation maslovchubukov :

[TABLE]

Therefore, the Landau damping contribution is given by Eq. (34) with the function $\mathcal{F}(z)$ given by Eq. (36). At small temperature $T\to 0$ the argument of $\mathcal{F}$ in Eq. (34) is large and one can use the asymptotic expansion:

[TABLE]

We note that Eq. (37) establishes the non-analytic behavior of the Landau damping contribution.

In order to calculate the Kohn anomaly contribution, we will use the trick from Ref. maslovchubukov First of all, for the Kohn anomaly only the difference $\delta E_{F}$ in chemical potential matters:

[TABLE]

where the index $K$ of the sum indicates that $Q$ is in the vicinity of the Kohn anomaly. The difference $\delta E_{F}$ in chemical potential is the following:

[TABLE]

where $\tilde{\mu}^{\sigma}_{\tau}$ is given by Eq. (12). Next, we notice that Eq. (38) is very similar to the $Q$-integral of $\mathcal{P}_{\uparrow}\mathcal{P}_{\downarrow}$ in Ref. maslovchubukov where $\mathcal{P}_{\uparrow}$ ($\mathcal{P}_{\downarrow}$) are particle-hole bubbles corresponding to the spin up (down) in a one-valley Fermi gas with Zeeman field. In Ref. maslovchubukov it is shown that one can rewrite the Kohn anomaly contribution of $\mathcal{P}_{\uparrow}\mathcal{P}_{\downarrow}$ as the Landau damping contribution of $\mathcal{P}_{\uparrow\downarrow}^{2}$, where $\mathcal{P}_{\uparrow\downarrow}$ is the particle-hole bubble constructed from up and down components of the Green function. Applying this trick to Eq. (38), we get:

[TABLE]

where the sum on the right is taken in the vicinity of Landau damping point. Applying Eq. (34), we obtain the Kohn anomaly contribution as

[TABLE]

where $\delta E_{F}$ is given by Eq. (39).

Using Eqs. (34) and (41) we get the non-analytic parts of the diagrams in Fig. 4:

[TABLE]

As $\delta\Omega_{b}=0$, the non-analytic corrections coming from other diagrams can be recognized by the interaction factors $v^{2}$, $u^{2}$ and $uv$. The non-analyticity is especially pronounced at $T=0$ where the asymptotic expansion Eq. (37) of the $\mathcal{F}$-function is valid. Performing the Legendre transform and taking into account the renormalization of order parameters due to the non-analytic correction $\delta\Omega$, we get Eq. (10).

Appendix C MAGNETIC SUSCEPTIBILITY IN THE PARAMAGNETIC PHASE

In this section, we calculate the non-analytic part of the magnetic field dependence of the magnetic susceptibility in the paramagnetic phase. In this paragraph we only account for the Zeeman coupling and neglect the orbital effects related to the vector-potential. The one-particle Hamiltonian in the paramagnetic phase with the magnetic field $\bm{B}=(B_{x},B_{y},B_{z})$ is the following kormanyos :

[TABLE]

where \mbox{\boldmath{\sigma}}=(\sigma_{x},\sigma_{y},\sigma_{z}) are the spin Pauli matrices, $g_{s}$ and $g_{v}$ are the spin and valley $g$-factors, $\mu_{B}$ is the Bohr magneton.

We define the magnetic susceptibility as follows:

[TABLE]

Here we are interested in $\chi_{zz}$ and $\chi_{xx}=\chi_{yy}$. Note that for the susceptibility $\chi_{zz}$ defined this way the formulation in terms of the linear response functions should be done using $\langle(g_{s}\sigma_{z}+g_{v}\tau)(g_{s}\sigma_{z}+g_{v}\tau)\rangle$ correlation function instead of the standard spin-spin correlation function $\langle g_{s}\sigma_{z}g_{s}\sigma_{z}\rangle$. This is because the magnetic field $B_{z}$ couples both to the spin and the valley, see Eq. (45).

If the magnetic field is perpendicular to the 2DEG plane, then we consider eigenstates of $\sigma_{z}=\sigma=\pm 1$. For the non-analytic correction to the grand canonical potential $\Omega$ we can readily use Eqs. (42)–(44) thanks to the correspondence

[TABLE]

As the SOI in TMDs is not smaller than a few meV, we always assume here that the magnetic field splitting is always smaller than the SOI. The corresponding non-analytic correction $\delta\chi_{zz}$ to the susceptibility at zero temperature is then linear in the magnetic field:

[TABLE]

where $f(x,y)=(|x+y|^{3}+|x-y|^{3})/2$, and $\chi_{0}=\nu\mu_{B}^{2}(g_{s}^{2}+g_{v}^{2})$ is the magnetic susceptibility of the non-interacting 2DEG. The linear non-analyticity in $\chi_{zz}$ is due to the cubic non-analyticity in the grand canonical potential. The same non-analyticity results in the first order phase transition. Therefore, measuring the unusual linear dependence of the magnetic susceptibility on the magnetic field allows to probe directly the non-analytic cubic correction to the grand canonical potential $\Omega$.

In case of an in-plane magnetic field the quantization axes in different valleys do not coincide. This leads to non-trivial spin traces in the diagrams in Fig. 4. Apart from this complication, the rest is the same. For example, instead of Eqs. (27)–(30), we get similar expressions:

[TABLE]

where the spin tensors $\mathcal{B}^{\sigma_{1}\sigma_{2}}$ are due to the non-collinear spin quantization axes in different valleys and given by

[TABLE]

Here, \Delta_{\tau}=|\mbox{\boldmath{\Delta}}_{\tau}| is the spin gap in valley $\tau$ with

[TABLE]

Using Eqs. (34) and (41) we get the non-analytic parts of the diagrams defined in Eqs. (49)–(52),

[TABLE]

In case of $B_{x}=B_{y}=0$, Eqs. (55)–(57) coincide with Eqs. (42)–(44) upon the correspondence outlined in Eq. (47).

In case if $B_{z}=0$ and $T=0$, Eqs. (55)–(57) can be further simplified:

[TABLE]

where $\delta\Omega$ is the total non-analytic correction, $B_{\parallel}=\sqrt{B_{x}^{2}+B_{y}^{2}}$ is the in-plane magnetic field, $\Delta=\sqrt{\alpha^{2}+(g_{s}\mu_{B}B_{\parallel}/2)^{2}}$. It follows from Eq. (58) that there is no linear dependence of the in-plane susceptibility on the magnetic field in case $g_{s}\mu_{B}B_{\parallel}\ll\alpha$ as the SOI vector $(0,0,\alpha)$ is orthogonal to the in-plane magnetic field. Therefore, only the $\chi_{zz}$ component of the magnetic susceptibility is useful to detect the cubic non-analyticity.