Finite temperature fluctuation-induced order and responses in magnetic topological insulators

We derive an effective field theory model for magnetic topological insulators and predict that a magnetic electronic gap persists on the surface for temperatures above the ordering temperature of the bulk. Our analysis also applies to interfaces of heterostructures consisting of a ferromagnetic and a topological insulator. In order to make quantitative predictions for MnBi$_2$Te$_4$, and for EuS-Bi$_2$Se$_3$ heterostructures, we combine the effective field theory method with density functional theory and Monte Carlo simulations. For MnBi$_2$Te$_4$ we predict an upwards N\'eel temperature shift at the surface up to $15 \%$, while the EuS-Bi$_2$Se$_3$ interface exhibits a smaller relative shift. The effective theory also predicts induced Dzyaloshinskii-Moriya interactions and a topological magnetoelectric effect, both of which feature a finite temperature and chemical potential dependence.

Conceptually, a perfectly quantized Hall conductivity of e 2 /h arises at zero temperature when the magnetization couples to the topological Dirac-type TI surface states, opening up a gap there in which then the chemical potential µ must lie, see Fig. 1. As a bonus, fermionic quantum fluctuations induce a concomitant linear topological magnetoelectric effect (TME), which couples electric fields directly to the magnetization and, vice versa, magnetization dynamics to electrical polarization [27,28]. At zero temperature for µ outside the gap both Hall conductance and TME fail to be quantized [2,6,29], tending to vanish as µ grows. Nevertheless, this particular situation has the interesting feature, that a Dzyaloshinskii Moriya interaction (DMI) between magnetic degrees of freedom emerges, opening the path towards the formation of various skyrmion-like topolog-FIG. 1. Left: symmetry breaking induced by proximity effect. An exchange coupling is induced across the interface between a FMI and a TI. The FMI polarizes the TI surface by proximity-effect and gaps the surface spectrum like, e.g., heterostructures EuS-Bi2Se3 heterostructues [13]. Right: intrinsic spontaneous symmetry breaking. Here the TI is itself a magnetic insulator like, e.g., MnBi2Te4 [17].
The goal to realize the QAHE, TME and possibly even a DMI at room temperature ties into a number of fundamental and practical questions. The first is how in an ideal situation the surface magnetic ordering temperature T c is affected by coupling to the topological edge states -it has been suggested that this coupling can greatly enhance T c [12,33]. Subsequently the question is how temperature fluctuations affect the conductance quantization, TME and DMI for the different relevant regimes of µ. Particularly interesting would be the existence of a temperature regime in which both DMI and TME are sizeable in which case the TME endows external magnetic and electric fields with novel types of access to DMI-induced skyrmions.
Despite the several theoretical and experimental developments in recent years, a number of fundamental questions remain to be answered. For instance, although arXiv:2101.08645v2 [cond-mat.str-el] 13 Jul 2021 a remarkable enhancement of T c at interfaces of certain ferromagnetic insulators (FMIs) have been reported [12,34], recent works [35,36] question the validity of these findings for the specific case of EuS proximate to either Bi 2 Se 3 or (Bi,Sb) 2 Te 3 . Furthermore, a recent experimental work [35] indicates that the topological electronic states at the interface do not interact strongly with ferromagnetism for the case of a EuS-Bi 2 Se 3 heterostructure. Additionally, in the family of magnetic TIs MnTe(Bi 2 Te 3 ) m , different works that do find a surface spectrum gapped below the Néel temperature also observe the persistence of the gap in the paramagnetic phase [17,22,37], it being unclear whether or not these observations result from the intrinsic magnetism.
Here we develop the finite temperature continuum field theory to address the questions concerning the magnetic phase transitions and dynamics at finite temperature, and apply these results to several experimentally relevant materials systems, using Monte-Carlo simulations (MCS) and density functional theory (DFT) based approaches to obtain quantitative results. Using a minimal model for the coupling of the Dirac fermions to the magnetic Hamiltonian, we show that in a temperature window where bulk magnetism is absent, an out-of-plane surface magnetization can still be nonzero and induce a gap in the Dirac spectrum. As a consequence, the AHE and TME can survive in a certain temperature range above the bulk T c . However, for the experimentally relevant materials and material combinations (EuS-Bi 2 Se 3 , MnBi 2 Te 4 ), the coupling to the topological surface states enhance the bulk T c not more than 15%, even under the most favorable conditions. Significant enhancements would require using TIs with a much lower Fermi velocity of the Dirac cone. In addition, at finite µ we establish the existence of a temperature regime displaying both a substantial fluctuation-induced DMI and TME, even if its Hall conductivity is strongly renormalized, with potentially interesting consequences for skyrmion manipulation and transport.
Saddle-point and induced order -To determine the shift in magnetic ordering temperature due to coupling of the magnetic moments to the fluctuating Dirac fermions we consider the following minimal model Hamiltonian where the Dirac fermions couple to the magnetization via a magnetic exchange interaction J 0 , σ σ σ is the Pauli matrix vector and n(r, t) the unit vector field representing the magnetization direction at r = (x, y). The operator d(−i∇) has the property d 2 = −∇ 2 , where ∇ = (∂ x , ∂ y , 0). Additionally, an electric potential φ has been introduced, which includes contributions of an externally applied electric field and an internal long-range Coulomb interaction as well.
The fermionic quantum fluctuations of the Dirac Hamiltonian (1) are accounted for by the imaginary time path integral, where Ψ = (Ψ ↑ , Ψ ↓ ) T is a spinor of Grassmann fields obtained from the second-quantized Hamiltonian [38]. The above partition function defines a free energy functional F F (n) which provides an additional free energy to the one of the magnetic free energy. As a minimal model leading to the latter, we consider the magnetic Hamiltonian, where J > 0 is the exchange energy and K > 0 is the anisotropy energy density (per unit area). The magnetic partition function is given by the path integral, where S B is the Berry phase that arises in the construction of the spin coherent state path integral [39], and λ is a Lagrange multiplier field enforcing the constraint n 2 = 1. Keeping the magnetic fluctuations classical, we obtain the following effective Hamiltonian after integrating out the Gaussian fluctuations n x and n y , along with the fermions, Variation with respect to n z leads to the saddle-point equation, where E q = ( v F q) 2 + m 2 , ω n = πk B T (2n + 1)/ is a fermionic Matsubara frequency, and we have defined m 2 = J 2 0 n 2 z . Equation (7) is solved together with the saddle-point equation for λ, which occurs at iλ = λ 0 , Setting J 0 = 0 in Eq. (7) reduces the saddle-point equations to one of a classical ferromagnet with easy-axis anisotropy. In this special case the ordered phase immediately implies λ 0 = K and from Eq. (8) it is straightforward to obtain the critical temperature T c by demanding that n z (T c ) = 0, yielding where a cutoff Λ s K/J has been introduced. Note that that the above is consistent with the Mermin-Wagner theorem in the limit K → 0.
Our aim is to calculate the shift of this critical temperature when J 0 = 0, i.e., accounting for the fermionic quantum fluctuations. After explicitly evaluating the Matsubara sum and integral, Eq. (7) becomes where the cutoff Λ F |m|/( v F ). The chemical potential is temperature dependent and can be obtained by fixing the particle density. At the critical temperature T = T c of the interacting system we demand that m = 0 and obtain where µ c = µ( T c ). This finally yields the critical temperature shift relative to the situation where fermions are absent, Since the cutoff is large, it is clear that the argument of the logarithm in Eq. (12) is smaller than unity, and therefore T c > T c in all cases. From Eqs. (11) and (12) we see that smaller values of v F favor larger shifts of the critical temperature.
Material specifics -In order to make quantitative predictions for the material systems of interest, we need to determine the values of the coupling parameters and cutoffs appearing in the continuum theory. We base such values on ab-initio and Monte-Carlo calculations, which we find to be consistent with available experimental data.
Based on DFT for a finite slab, for MnBi 2 Te 4 we find the Fermi velocity as v F = 2.3 ± 0.2 eVÅ and the coupling J 0 ≈ 50 meV, whereas for EuS-Bi 2 Se 3 , we consider v F ≈ 3.3 eVÅ and J 0 ≈ 54 meV [13, 40,41]. Note, the value of J 0 was derived from the gap size in the Dirac dispersion assuming its origin is purely magnetic. In reality higher order effects such as hybridization of the surface fermions with the bulk electronic states of MnBi 2 Te 4 or EuS may also influence the gap such that the actual magnetic gap and correspondingly J 0 might be smaller.
For the fermionic cutoff Λ F , we consider that the average surface density of a completely filled band is 1/A, with A the surface unit cell area. Since our model describes two surface bands, we fix Λ F such that n(µ = 0) = 1/A. Electron-hole symmetry of the model then implies that n(µ → ∞) = 2/A.
To set the anisotropy K and stiffness J of the magnetization field, we first build an anisotropic Heisenberg lattice model which we then map to Eq. (4). As the on-site anisotropy depends crucially on the thickness of the FMI layer in the EuS-Bi 2 Se 3 system [12], while the Mn layers in MnBi 2 Te 4 are well separated with relatively small out-of-plane exchange couplings [17,42], we consider the magnetic subsystem to be monolayer thick. The corresponding lattice in both cases is a two-dimensional triangular lattice spanned by Mn in MnBi 2 Te 4 and Eu on the EuS(111) surface [12]. The considered magnetic interactions comprise of nearest neighbor ferromagnetic exchange couplings J, and an effective on-site out-of-plane anisotropy K. J and K follow from J and K taking the continuum limit [43].
For MnBi 2 Te 4 , we obtain from DFT calculations K ≈ 0.073 meV, and J ≈ 0.18 meV, in good agreement with earlier reported values [44]. For monolayer EuS, on the other hand, K was obtained by extrapolating the data for the 20 nm thick Bi 2 Se 3 layer in the EuS-Bi 2 Se 3 heterostructures in Ref.
[12] to the EuS monolayer thickness. We obtain K ≈ 0.13 meV. The exchange coupling for the monolayer has been previously estimated to be J = 0.017 meV [45]. It is interesting that these two material systems cover a broad range of K/J (from ∼0.4 in the Mn-based compound to ∼7.6 in the EuS-based system).
Last, the spin cutoff was fixed such the critical temperature of the continuum model (without fermions) matches the critical temperature of the corresponding Heisenberg model. To obtain the latter, classical MCS for the lattice model were carried out (see SI [43] for details). For the EuS-Bi 2 Se 3 and MnBi 2 Te 4 system, we obtain T latt c ≈ 5.8 K and ≈ 17.0 K, respectively. Based on these values, we fix Λ s via Eq. (9).
Finally, using Eq. (12), we obtain (T c − T c )/ T c = 10.9% and 14.7% for EuS-Bi 2 Se 3 and MnBi 2 Te 4 , respectively. Note that the compound having smaller Fermi velocity (MnBi 2 Te 4 ) shows indeed a larger shift of the critical temperature (recall that in both cases J 0 is similar).
Fluctuations around the saddle-point -The question on how the DMI, TME and Hall conductivity evolve with temperature and chemical potential requires considering the effect of magnetic fluctuations in the fermionic determinant resulting from integrating out the fermions. Introducing n(r, t) = n zẑ + δn(r, t) we can determine the effective action around the saddle-point approximation up to quadratic order in the fluctuations.
The reason for the appearance of the DMI term [30] is that the magnetization fluctuations effectively break the inversion symmetry of our starting Hamiltonian. As a consequence, this yields a DMI contribution to the magnetic free energy, where, It is important to emphasize that the DMI term is not introduced in an ad hoc way -it is generated by charge fluctuations coupling to the magnetic moments at the interface. The DMI vanishes at the neutrality point and is nonzero away from it. This creates the possibility of manipulating the DMI by controlling the chemical potential, for instance by gating. If we take the mass m to have a mean-field like behavior m(T ) = J 0 1 − T / T c we find the zero temperature value for the DMI D( which demonstrates that the DMI kicks in when the Fermi energy F surpasses a threshold given by the exchange coupling constant J 0 . This feature of the generated DMI can also be seen in Fig. 2 (a) where we show the whole temperature range for different values of the chemical potential with an estimated zero temperature exchange coupling constant of J 0 /k B T c ≈ 29.1 based on our findings for the MnBi 2 Te 4 system. We further see that the lower the temperature gets, the narrower the range for the chemical potential becomes in which there still is a finite DMI. Moreover, we see that the step-function behavior of the generated DMI at zero temperature approximately extends to the whole temperature range as the DMI is only present when the chemical potential exceeds the magnetic gap, meaning it exists only in the metallic regime.
Finally, we determine the fluctuation-induced effective Chern-Simons (CS) action, where we have defined the covariant three-potential A µ = φ c , ± J0 evF d(δn) . The electric potential enters in the time component as usual and the magnetization fluctuations δn act as the vector potential A in the spatial components. Note that the "±" applies to the different choices for the vector d. The coefficient σ arising in Eq. (15) leads to the gap in magnetic susceptibility, in a mechanism closely related to the well known topologically massive photons in a Maxwell-Chern-Simons theory [46]. In our case this topological mass is given by, sinh(βm) cosh(βµ) + cosh(βm) .
We recall here that the topological mass arising in the effective free energy is in general not identical to the Hall conductivity -these quantities differ, for instance, in the metallic regime [6,29], something that is more easily seen in the zero temperature limit. Indeed, for T = 0 the topological mass and Hall conductivity are given by respectively. These zero temperature expressions involving the Heaviside step function are identical only when F < |J 0 |. In fact, σ(T = 0) vanishes in the metallic regime while σ xy is nonzero. This occurs because the Hall conductivity is calculated from the Kubo formula where one first takes the limit q → 0 and then ω → 0, while in case of the topological mass these limits are taken simultaneously. The CS action (15) contains the TME contribution to the free energy, To compare its features to the ones of the DMI we also illustrate its dependency on temperature and chemical potential in Fig. 2 (b). Once more it shows that at zero temperature we have a step function behavior which also approximately extends to finite temperatures. Consequently, the topological mass only exists when the chemical potential lies inside the magnetic gap and is nearly quantized in the bordering regions resulting into plateaus. In comparison to the generated DMI we can see in Fig. 3 that there is a very narrow region where both functions overlap. As a result, the desired simultaneous occurrence of the DMI and the CS action requires fine tuning.
However, upon closer inspection, a different connection between the two terms appears. It turns out that the temperature functions inside both terms complement each other almost perfectly in a temperature and chemical potential plot, as also can be seen for exemplary temperatures in Fig. 3. The two functions are adding up to one creating a plateau that even traverses the chasm that both functions showed individually in the vicinity where the chemical potential crossed the magnetic gap. This means that at the time one of the terms diminishes, the respective other term grows in size equal to the loss of the other, creating a direct correspondence between them. We point out that besides the DMI and CS terms, other interesting terms appear in the effective action, as is shown in explicitly in the SI [43].
Conclusion -We have considered a minimal model for magnetic topological insulators that capture a wealth of interesting properties of the surface of materials like MnTe(Bi 2 Te 3 ) m with m ≥ 1 [17,19,21,22] and MnSb 2 Te 4 [23,24], and also at interfaces of heterostructures EuS-Bi 2 Se 3 [12-14]. An important main result of our analysis is the prediction of the survival of the electronic gap on the surface/interface for temperatures above the bulk ordering temperature. In order to provide quantitative results to be compared with experiment, we have combined the effective field theory analysis with DFT and MCS results applied specifically to MnBi 2 Te 4 and the bilayer system EuS-Bi 2 Se 3 . We also predict that temperature dependent DMI and TME terms are induced by fluctuations. The latter may give rise to new magnetic phenomena at the surface of magnetic TIs, including the interesting possibility of manipulating skyrmions by external or internal electrical fields.

Saddle-point and induced order
The starting point of our field theory is the partition function from the main text, consisting of two parts, namely the magnetic and fermionic contributions. We calculate this partition function in the imaginary time path integral formalism integrating out the fermionic and bosonic fields respectively. Keeping the magnetic fluctuations classical, we look at the magnetic partition function and define a corresponding magnetic action where the magnetic action is given bŷ Here we already split up the different components of the magnetization into in-plane and out-of-plane components, where the sum over a denotes the x and y components, such that in the next step we can perform the path integration over the Gaussian fluctuations n x and n y yielding the effective action Doing the same for the fermionic partition function, just by integration over the fermionic fields we get To note is the additional minus sign in front, resulting from the integration over fermionic fields instead of bosonic ones. Furthermore, the two so defined functional traces differ from each other, which will soon become apparent. Combining the two individual contributions we get the total effective action which implies the effective Hamiltonian given in the main text.
In the next step we apply the saddle-point approximation to the remaining fields to be integrated over n z (r) =ñ z + δn z (r), leading to the saddle point effective action To keep the notation short we redefineλ = λ andñ z = n z by dropping the tilde over both quantities. With this we calculate the functional traces Tr ln −J∇ 2 + iλ = d 2 r r ln −J∇ 2 + iλ r , where which upon reinsertion gives Similarly, the functional trace from the integration over the fermionic fields is given by where with the fermionic Matsubara frequencies ω n = (2n + 1)π/ β. Therefore the functional trace yields where the sum runs over σ ∈ {−1, 1} and we defined E q = 2 v 2 F q 2 + m 2 with the mass m = J 0 n z . Using this we can now formulate the saddle point equations by variation of the effective action with respect to n z and with respect to λ 0 = iλ Setting J 0 = 0 in Eq. (31) reduces the saddle-point equations to one of a classical ferromagnet with easy-axis anisotropy. In this special case the ordered phase immediately implies λ 0 = K and from Eq. (32) it is straightforward to obtain the critical temperature T c by demanding that n z (T c ) = 0, yielding where a cutoff Λ s K/J has been introduced. Our aim is to calculate the shift of this critical temperature when J 0 = 0, i.e. accounting to the fermionic quantum fluctuations. After explicitly evaluating the Matsubara sum and integral, Eq. (31) becomes, where we have assumed that the cutoff Λ F |m|/( v F ). The chemical potential is temperature dependent and can be obtained by fixing the particle density. At the critical temperature T = T c , we demand that m = 0 and obtain where µ c = µ( T c ). Analogously to Eq. (33) we then get This finally yields the critical temperature shift relative to the situation where fermions are absent, Since the cutoff is large, it is clear that the argument of the logarithm in Eq. (37) is smaller than unity, and therefore T c > T c .

Surface particle density and fermionic cutoff
The fermionic cutoff Λ F can be determined from the two dimensional surface particle density given by where in our case the energies are E = ±E q − µ resulting from the Dirac Hamiltonian. By insertion of these energies into Eq. (38) the surface particle density becomes where I is given by with the notation Li n [x] for the polylogarithm. We now consider that the average surface density of a completely filled band is 1/A, with the surface unit cell area A. Since our model describes two surface bands, we fix Λ F such that n(µ = 0) = 1/A. At µ = 0 the Integral I vanishes, giving an expression for the fermionic cutoff Λ F = 4πn(µ = 0).

Mapping between lattice and continuum spin models
We start with an anisotropic Heisenberg model on a two-dimensional triangular lattice: where, J ≥ 0 is the nearest neighbor ferromagnetic Heisenberg exchange coupling and K is the on-site magnetic anisotropy. Introducing n(r i ) = S i /S and ∆R = r j − r i as the distance vector between lattice site i and j the Hamiltonian becomes To estimate the Fermi velocity of the surface state and the surface gap, we performed DFT calculations for a slab structure consisting of six MnBi 2 Te 4 unit cells with a vacuum of 30 Bohr radii [ Fig. 4(a)]. Following Ref. [17,50], we fix U = 5.34 eV and J = 0. The two main surface bands present a gap ∼100 meV, as shown in Fig. 4(b), in good agreement with the value of 88 meV found in [17]. The Fermi velocity found for the upper part of the Dirac cone is approximately 2.3 ± 0.3 eVÅ in good agreement with the experimental results in Ref. [51] and DFT results in Ref. [17].
Lastly, samples of MnBi 2 Te 4 tend to be self-doped, meaning that different kind of defects place the chemical potential µ outside the gap. In particular, the samples tend to be electron-doped. To estimate the value of µ, we consider the estimation of carriers of n c = 2 × 10 19 /cm 3 provided in Ref. [17]. Based on this value and on our slab calculation we estimate µ ∼ 160 meV above the bottom of the conduction band.
Notice that the value of µ = 160 meV with respect to the bottom of the conduction band corresponds to 210 meV with respect to our zero of energies, as we have to add half the size of the gap.

EuS-Bi2Se3 heterostructures
In EuS-Bi 2 Se 3 heterostructures, the interface is typically formed by the (111) surface of the cubic bulk-EuS structure such that the lattice mismatch with the topological Bi 2 Se 3 film is minimal [52]. For this surface, the interfacial layer of Eu atoms span a triangular lattice. The effective lattice constant of this lattice is a = a EuS / √ 2 ≈ 4.22Å, using a EuS ≈ 5.96Å [53].
The Fermi velocity was obtained from [41] to approximately be v F ≈ 3.29 eVÅ for bulk Bi 2 Se 3 . Furthermore, the exchange coupling constant J 0 was estimated in alignment with the magnetic gap reported in [13,40] to have a value of 54 meV.
Regarding magnetism in the EuS-Bi 2 Se 3 , the value of J = 0.017 meV has been reported earlier [45]. The value of the on-site magnetic anisotropy was obtained from the EuS layer thickness dependence of the magnetic anisotropy [12]. We considered the structure with the largest Bi 2 Se 3 layer thickness of 20 nm. Note that these values are available in the continuum limit and were converted to the lattice equivalent values using Eq. (50). The resulting data was modeled with [54]: where K V and K S , respectively denote the bulk and surface magnetic anisotropy contributions and d is the thickness of the EuS layer. For the EuS monolayer along the (111)-direction, with a thickness of d = a EuS / √ 3 ≈ 3.45Å, we obtain K ≈ 0.126 meV, leading to K/J ≈ 7.4.

Monte-Carlo calculations of Tc
Classical Monte-Carlo simulations (MCS) with the Metropolis algorithm were carried out for spins on a twodimensional triangular lattice with 42 × 42 sites. We consider the spin Hamiltonian of Eq. (41). JS 2 = 1 defines the energy scale leaving the ratio K/J as the only free parameter. For each K/J, we started from a high-temperature paramagnetic state, characterized by a random spin configuration, and decreased the temperature in steps of 0.02. At each temperature, the system was allowed to equilibrate over N eq steps and the physical quantities were obtained by averaging over the next N av steps. For K/J < 1, the equilibration was typically reached in 5×10 4 , however, to treat the entire range of K/J < 10 on same footing, we generously consider N eq = 2 × 10 5 and N av = 3 × 10 5 steps. The critical temperature, T latt c , was obtained from peak(s) in the specific heat, which agrees with corresponding values obtained from the magnetization data M vs. T .
To address the materials of interest in this study, we used K/J = 0.40 for MnBi 2 Te 4 and K/J = 7.4 for the EuS-Bi 2 Se 3 heterostructure as discussed earlier. Figure 5 (a) shows the corresponding specific heat data. From the well-defined peaks, we obtain k B T latt c /JS 2 ∼ 1.46 and ∼ 2.40, for the MnBi 2 Te 4 and EuS-Bi 2 Se 3 monolayers, respectively. Considering S = 5/2 for MnBi 2 Te 4 , we obtain T latt c ∼ 16.97 K, while for the EuS-Bi 2 Se 3 heterostructure, S = 7/2 yields ∼ 5.80 K.
For completeness, we also carried out MCS over a wide range of K/J. This allows us to study the evolution of T latt c as a function of K/J and to analyze how close or far are the systems of interest from the Ising limit K/J 1. With increasing K/J, the accessible phase space becomes considerably smaller. Consequently, the equilibration is slower. Therefore, for K/J ≥ 10, MCS were carried out with a total of 1 × 10 6 update steps, out of which the first 5 × 10 5 steps correspond to N eq and were discarded during the averaging and evaluation of the physical properties such as specific heat. Figure 5 (b) shows the evolution of T latt c with K/J. The Ising limit, which corresponds to k B T latt c /JS 2 = 3.642 [55], is approximately reached for K/J 100. where the Greek indices label the spin components resulting from the trace over the Pauli matrices and where we definedĜ k,n = (i ω n + µ) + v F d(k) · σ σ σ − mσ z (i ω n + µ) 2 − 2 v 2 F d 2 (k) − m 2 .
By insertion of the definition of V into Eq. (55) one can now split up the different parts of the action according to their magnetoelectric nature, yielding an electric, magnetic and magnetoelectric contribution.
To arrive at the DMI and CS terms we then trace out the spin components, and calculate the integral over k and the sum over fermionic Matsubara frequencies iω n in a derivative expansion, that means a long wavelength and low frequency expansion in terms of the wavevector q and bosonic Matsubara frequencies iν m . This step is straightforward but very lengthy and shall be omitted here.
Afterwards one transforms back to Euclidean spacetime using the remaining Matsubara sum and momentum integral over iν m and q to find the result for the DMI and CS terms presented in the main text, among other contributions.