Wiedemann-Franz law and Mott relation for non-linear anomalous transport phenomena

Berry curvature dipole has been shown to produce non-linear anomalous Hall effect which is non-zero even in the presence of time-reversal symmetry. In this paper, within the framework of semiclassical Boltzmann theory, we calculate the analytical expressions for the non-linear anomalous transport coefficients, namely, the non-linear anomalous Hall, Nernst, and thermal Hall coefficients. With these expressions, we predict the fundamental relations between the non-linear anomalous electric and thermo-electric transport coefficients, which are the analog of the Wiedemann-Franz law and the Mott relation in the non-linear regime. We also make experimental predictions for anomalous thermal Hall coefficient for monolayer transition metal dichalcogenide which can be verified in experiments.


I. INTRODUCTION
Conventionally, the intrinsic anomalous Hall effect is triggered by the anomalous velocities induced by a nontrivial Berry curvature Ω k , which requires the breaking of time reversal (TR) symmetry [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] .On the other hand, the Berry curvature dipole, defined as the first order moment of the Berry curvature (∂ k Ω k ) over the occupied states (f k ), can survive under the TR symmetry and is found responsible for the nonlinear anomalous Hall effect in the TR invariant but inversion symmetry broken systems [16][17][18] .Since then, great progress has been made to show the important roles that Berry curvature dipole plays in the non-linear anomalous electrical as well as optical responses [19][20][21][22][23][24] .Specifically, in some recent electric transport experiments 25,26 , the non-linear anomalous Hall effect has been successfully observed in the few-layer WTe 2 , a two-dimensional transition metal dichalcogenide (TMDC) material with broken inversion symmetry and a single mirror symmetry.The electric non-linear anomalous Hall current in the time reversal systems is defined by j a = abc χ abc E b E c with χ abc ∝ D bd , where D bd is the Berry curvature dipole and the indices b, d are associated with its components 17 .Because of their noncentrasymmetric band structures, large spin-orbit coupling and strain-tunability, the 2D TMDCs have draw a great interest in the studies of Berry curvature dipole and the related applications 24,27,28 .
For the non-linear anomalous Hall effect, it is the nonequilibrium Fermi distribution perturbed by the external electric field E (and computed to be proportional to E), that brings the nonlinearity to the Hall response 17,18 .Following a similar analogy, in recent works the non-linear anomalous Nernst effect has been proposed in strained monolayer MoS 2 29 and bilayer WTe 2 30 , where a non-equilibrium Fermi distribution is introduced by the thermal gradient (−∇T ), and a transverse electric current is generated as a second order response to the thermal gradient (which could be defined as j a = abc α cd ∇ bd T ).The non-linear anomalous Hall effect and the non-linear anomalous Nernst effect are both revealed by the charge current, and they are expected to co-exist in the TR invariant systems through satisfying the Onsager reciprocity 31,32 .However, the Bloch electrons can not only carry charges, they also carry energy and heat.Therefore, we expect the possible existence of the non-linear anomalous thermal Hall effect in the presence of TR symmetry.
The thermal Hall effect, describing the generation of a transverse heat current resulting from a longitudinal thermal gradient, has been widely studied in the linear regime for electrons [33][34][35] , and also for charge-free quasi-particles such as phonons [36][37][38] , magnons [39][40][41] and photons 42 .In this paper, we focus only on the anomalous thermal transport carried by electrons, which is shown to be manifested by the intrinsic Berry phase effect through the semiclassical wave packet approach 2,[43][44][45] .In the linear transport, however, the anomalous thermal Hall current is forced to vanish in a TR invariant system because of the symmetry property of Berry curvature Ω k [46][47][48][49][50][51][52][53] .
The non-reciprocal thermal transport has been studied for the 1H monolayer MoS 2 and the polar semiconductor BiTeX (X=I, Br) in a very recent work 54 , where the nonlinearity of the thermal response (∝ (∇T ) 2 ) is induced by the non-linear distribution function that is proportional to the square of the relaxation time (∝ τ 2 ) and is not Berry curvature related.However, the intrinsic Berry curvature along with a linear nonequilibrium distribution function (∝ τ ) could also contribute to the non-linear conductivities, such as the non-linear anomalous Hall effect and non-linear anomalous Nernst effect in previous works 17,18,25,26,29,30 .Based on this, here we propose a new type of non-linear effect, the non-linear anomalous thermal Hall effect (of which the thermal current is defined by (j Q ) a = abc l cd ∇ bd T , see Eq. ( 13)), which is found to be non-zero in the TR invariant systems.To our best knowledge, this non-linear anomalous thermal Hall effect has never been studied before.For measuring the non-linear thermal Hall effect, a schematic experimental setup is given in Fig. (1).
In this work, we systematically derive the expression for the non-linear anomalous thermal Hall coefficient based on the Boltzmann semiclassical approach with a relaxation time approximation.By doing the Sommerfeld expansion in the low temperature limit, we explicitly show a T 2 temperature dependence for the non-linear anomalous thermal Hall coefficient, and relate it to the Berry curvature dipole, which has already been shown playing an important role in the non-linear anomalous Hall effect both in theories and experiments.In addition, based on our results, we show that the Wiedmann-Franz law applicable in the linear transport is violated in the non-linear regime in a way greatly different from previous studies 47,[55][56][57][58][59] .Using a similar approach, the non-linear anomalous Nernst coefficient is found to be related to the nonlinear anomalous Hall coefficient, and remains non-zero even at the limit of zero temperature.Although the non-linear or unconventional thermal Hall transports have been studied 54,57 , so far there is no theory describing the non-linear anomalous thermal Hall coefficient based on which, the thermal gradient induced non-linear conductivities (non-linear anomalous thermal Hall coefficient and non-linear anomalous Nernst coefficient) can be related to the non-linear anomalous Hall conductivity, or to be more specific, the Berry curvature dipole.
The rest of the paper is organized as follows.In Sec.II, after briefly presenting the Berry phase induced anomalous thermal and thermo-electric transport coefficients in the regime of linear response, the non-equilibrium Fermi distribution function under the thermal gradient (−∇T ) is derived by solving the semiclassical Boltzmann equations with a relaxation time approximation and the non-linear anomalous thermal Hall coefficient is formulated.In Sec.III, we derive the temperature (T ) dependence and chemical potential (µ) dependence of the non-linear anomalous thermal Hall coefficient l cd , by using the Sommerfeld expansion in a low temperature limit.We also check the Wiedemann-Franz law and the Mott relation in the non-linear regime, which are found to be different from those in linear transport regime.Based these results, we show that the non-linear anomalous thermal Hall coefficient, non-linear anomalous Nernst coefficient, and the non-linear anomalous Hall coefficient are related to each other through the Berry curvature dipole D bd (µ) and the derivatives of the Berry curvature dipole with respect to the chemical potential ∂ µ D bd (µ).In Sec.IV, we introduce the model Hamiltonian of a real physical system, the strained monolayer MoS 2 , which is a TR invariant 2D TMDC that has been intensively investigated recently.Then in Sec.V, we apply our theoretical expressions to this TR symmetric system.Our results show consistency between the analytical predictions and the numerical results within an appropriate low temperature range, indicating the validity of our theoretical analysis.Here, we also discuss some experimental applications of our results.We end with a conclusion in Sec.VI.

A. Anomalous transport equations in linear regime
With a non-trivial Berry curvature, the conventional group velocity of the Bloch electrons acquires an anomalous term 1,2 , which is written as, where v k = −1 ∇ k ε k with ε k the energy dispersion, ṗ = eE + e ṙ × B is the momentum with E, B the electric and magnetic field respectively, and Rather than being a point particle, the Bloch electron in the quantum picture has a finite spread around its mass center (r c , k c ) in the phase space, which could be described by a wave packet |W 2,43 .Because of the finite size, the wave packet generally carries an orbital magnetic moment (m(k)) coming from the self rotation motion around its mass center.Together with the motion of the mass center, a net magnetization (M ) is formed, which in turn produces a magnetization current that should be subtracted from the local current in describing the anomalous transport phenomenon 43 .In the wave packet approach, a general expression for the electric transport current is given by 2,43 , where f k = (e −β(ε k −µ) +1) −1 is the Fermi distribution function of the carrier, ṙ is the electron velocity given by Eq. ( 1).The second term is an integral over the energy (ε = ε k ), where σ 0 (ε) is the zero temperature anomalous Hall conductivity, defined as, Here Θ(x) is the Heaviside step function.Note that in Eq. ( 2), the anomalous Hall response to an electric field, could come from both the ṙ term and the chemical potential gradient term (−∇µ/e).On the other hand, the anomalous thermal Hall response only comes from ∇T in the second term.Similarly, we could find the transport heat current with Berry-phase correction, which is given as, where the first term, with Ω k = 0, is the usual expression for the heat current j Q N in the non-topological systems.We can write the total transverse heat current as j Q trans = j Q E + j Q T , where j Q E is the transverse heat current generated by the electric field, given as the following equation, Note, the electric filed here comes from the electrochemical potential 38,41,44 (μ = −eφ + µ).Similarly, the temperature gradient induced transverse heat current j Q T could be extracted from the second term in Eq. ( 4), written as The linear thermoelectric response in the presence of electric field E and temperature gradient −∇T are given by the transport equations below: where j, j Q are the charge and heat current respectively, σ, α, l are the electrical, Peltier, thermal conductivity tensors respectively.Tensor α is related to α through Onsager's relation (α = T α) 32 , and tensor α could usually be extracted from electrical conductivity and thermopower (S = σ −1 α).
Based on Eqs. ( 2), ( 5) and ( 6), the anomalous Hall conductivity, anomalous Nernst conductivity and anomalous thermal Hall conductivity could be respectively given as, where σ 0 (ε) is given in Eq. ( 3).Note, the finite temperature anomalous Hall conductivity σ ab is found by integrating the zero temperature anomalous Hall conductivity σ 0 along with −∂f /∂ε.Note also, that the integrations for α ab , l ab are weighted by additional factors (ε − µ) and (ε − µ) 2 respectively.Though a Dirac-delta-function (δ-function) like term −∂f /∂ε is present, the anomalous linear transport coefficients are Fermi liquid quantities, since the energy integrals are over the whole energy range.Here, all these anomalous linear responses are scattering time independent, and are expected to be zero for time reversal symmetric systems, because of the symmetry property of Berry curvature i.e., T (Ω k )T † = −Ω −k .The expressions for σ ab , α ab in the nonlinear regime have already been studied 17,18,29,30 .In this paper, we focus on the non-linear anomalous thermal Hall effect (the non-linear contribution for l ab ).In the next part, we derive the equations for the non-linear anomalous thermal Hall effect using the Boltzmann transport fomalism.

B. non-linear anomalous thermal Hall effect
In the Boltzmann transport theory, the Fermi distribution function f (k, r, t) in phase space of position and crystal momentum is assumed to satisfy the Boltzmann equation with a collision integral that depends on the detailed collision process.The full Boltzmann equation is given as, In relaxation time approximation, the collision integral could be written as where g k = f (k, r, t) − f 0 is the difference between the perturbed Fermi distribution f k and its value f 0 in equilibrium state, τ is the average scattering time between two successive collisions, which could depend on energy, momentum or impurities potentials.For simplicity we treat the scattering time as independent of energy and momentum in this paper.
To find the response to the thermal gradient −∇T in the absence of the external fields (E = B = 0), we expand g k as where g n k is understood as the n th order response to the thermal gradient, i.e., g n k ∝ (∇T ) n .Thus, in this case, the carriers in a steady state satisfy the Boltzmann equation given by, It is now straightforward to write down the solutions for the components of g k , where the higher order terms g n k could be found by iteration.To derive an expression for the non-linear anomalous thermal Hall conductivity (l ab ), we substitute f k = f 0 + g k in Eq. ( 6).We can see that, besides the contribution from f 0 that gives us the linear component of j Q T and l ab , the perturbation term g k additionally provides the non-linear anomalous thermal responses to higher order of −∇T .In this paper, we only consider the non-linear anomalous thermal effect to the second order of −∇T , such that the non-linear anomalous thermal Hall current is given by, where the prime on (j Q ) indicates the non-linear response.Indices a, b, c, d represent the components x, y, z.Here, only the non-zero component of the non-linear anomalous thermal Hall current under the time reversal symmetry is considered.For the convenience of symmetry analysis, the nonlinear anomalous thermal Hall current in Eq. ( 13) is written as an integral in momentum space, where a summation of all the bands (indexed by n) is introduced.In analogy to the linear case, the non-linear response coefficient could be extracted from (J Q ) a = l cd (∇ bd T ) as, Here the subscripts c, d indicate the index for the Berry curvature Ω n k,c and the momentum k d respectively.In contrast to the anomalous thermal Hall conductivity in linear regime ( Eq. ( 8)), the non-linear anomalous thermal Hall coefficient in the above equation is τ dependent, and the term ∂f0 ∂k d makes the coefficient a Fermi surface quantity with main contribution coming from the occupied states near the Fermi surface, which is different from the linear anomalous transport coefficients.Also, the appearance of the derivative of f 0 with respect to k d renders the integral for l cd non-zero under TR symmetry, i.e., the non-linear anomalous thermal Hall exist in the inversion symmetry broken systems even in the presence of time reversal symmetry.
With the above expression for l cd , in the next section we will extract the parameter dependence of the non-linear anomalous thermal Hall coefficient and find its relations with Berry curvature dipole in a low temperature regime.

III. SOMMERFELD EXPANSION IN LOW TEMPERATURE REGIME
In this section we will derive the analytic expression for the temperature dependence of the non-linear anomalous thermal Hall coefficient and show a direct relation with the non-linear anomalous Hall coefficient.This will be the analog of the Wiedemann-Franz law in the non-linear regime.To do this, we use the Sommerfeld expansion 32 , where g(ε) is a general function and First, we recall the expression for non-linear anomalous Hall coefficient given as 17,24 , where the frequency dependence (ωτ 1) could be neglected in experiment 24 .Here, D bd is the Berry curvature dipole defined as, Here, b, d are indices for ∂f k /∂k b and Ω k,b respectively.Comparing Eq. ( 14) with the above equation, it is seen that, the only difference is the additional unitless factor β 3 (ε k − µ) 3 , which modifies the Berry curvature distribution in momentum space.To apply the Sommerfeld expansion, we rewrite the Berry curvature dipole term as, with the Berry curvature dependent term G(ε) defined by, Note that, Berry curvature dipole D bd is a Fermi surface quantity because of the δ-function With the help of Eqs. ( 15) and ( 18), D bd could be expanded as, where the first term G(µ) is the zero-temperature Berry curvature dipole at Fermi energy (chemical potential) µ.This quadratic T -dependence of Berry curvature dipole seems to have been observed in a recent experiment in Ref. [25].Using the expansion in Eq. ( 20), we find the following analytical expression for the non-linear anomalous Hall coefficient, where χ 0 (µ) = e 3 τ G(µ)/2 2 is the zero-temperature nonlinear anomalous Hall coefficient, and χ (2) 0 (µ) denotes the second-order derivative of χ 0 (µ) with respect to chemical potential µ.
Recall that, the non-linear anomalous Nernst coefficient in the inversion symmetry broken but time reversal invariant systems are given by 29,30 , which could be rewritten as below in a similar format as D bd in Eq. ( 18), with G 2 (ε) defined as Here, the indices c, d for G 2 (ε) (and also G 3 (µ) in Eq. ( 27)) are different from the indices b, d for G(ε) given in Eq. (19).For a given time reversal symmetric system, however, it will be shown later that the indices of the non-zero component of these non-linear anomalous transport coefficients (i.e., χ abc , α cd and l cd ) are exactly the same, regardless of the different notations in the proceeding derivations.Similarly, using Eq. ( 15) for the above α cd , we could rewrite the non-linear anomalous Nernst coefficient α cd as, Using Sommerfeld expansion we can also analyze the non-linear anomalous thermal Hall coefficient l cd given in Eq. ( 14).First we rewrite l cd as with G 3 (ε) defined as, Note that, through G 3 (µ) (or G 2 (µ) in Eq. ( 24)), the nonlinear anomalous thermal Hall coefficient (non-linear anomalous Nersnt coefficient) is explicitly related to the non-linear anomalous Hall coefficient in Eq. ( 21), which will be shown in our following discussions.Applying Eq. ( 15) for l cd in the above equation, we find, with the higher order derivatives G (n) (µ) (odd number n ≥ 3) included in O(T 4 ).
So far, we have derived the analytical expressions for the temperature dependence of the non-linear anomalous Hall coefficient χ abc (Eq.( 21)), the non-linear anomalous Nersnt coefficient α cd (Eq.( 25)) and the non-linear anomalous thermal Hall coefficient l cd (Eq.( 28)) based on the Sommerfeld expansion approach.These analytical non-linear anomalous coefficients are the general results that could be applied to any systems with the corresponding non-linear effect present.Now, we can analyze the relations among them.
In the linear transports, it is known that, the amount of charge and heat carried by electrons are closely connected through the Wiedemann-Franz Law 60 .According to this law, the ratio of the the thermal and the electric conductivities is given by the absolute temperature multiplied by a universal constant known as the Lorenz number L, i.e., κ ab σ ab = LT with L = π 2 k 2 B /3e 2 , and σ ab = σ 0 (Eq.( 3)) the zerotemperature anomalous Hall conductivity in the linear regime.Based on Eqs. ( 21) and (28), however, we have the following relation between the anomalous thermal and the anomalous electric conductivies in the non-linear regime, where L 0 = k 2 B π 2 /e 2 = 3L, and χ 0 (µ) is the zerotemperature non-linear anomalous Hall coefficient as given in Eq. ( 21).This equation is different from the conventional Wiedemann-Franz law given by Eq. ( 29).Comparing Eq. ( 30) and ( 29), we see that, the temperature dependence is no longer a linear one (but ∝ T 2 ) in the non-linear regime.Additionally, rather then being proportional to the zero-temperature anomalous Hall coefficient (∝ σ ab ), the non-linear anomalous thermal Hall coefficient is directly proportional to the first order derivative of the zero-temperature non-linear anomalous Hall coefficient with respect to the chemical potential (∝ ∂χ 0 /∂µ).The results in Eq. ( 30) should be taken as a consequence of the intrinsic nonlinearity, rather than a conventional departure from the Wiedemann-Franz law [55][56][57][58][59] .
The differences for Wiedmann-Franz Law in the non-linear regime leads us to examine the Mott relation, which relates the electric and thermo-electric conductivities.In the linear regime, the Mott formula is written as 32 , In the non-linear regime, however, using Eqs.( 25) and (21)  the relation between the non-linear anomalous Nernst coefficient and the non-linear anomalous Hall coefficient is found as This Mott formula in the non-linear regime shows a finite value for the non-linear anomalous Nernst coefficient (α ab ) even at zero temperature.Note that, for both Eqs. ( 30) and (32), only the leading order of their Sommerfeld expansions is considered.Based on the above results, we see that, while the Wiedemann-Franz law in the linear regime does not contain the derivative of the anomalous Hall conductivity (Eq, ( 29)), in the non-linear regime the analog of the Wiedemann-Franz law explicitly involves the first order directive of χ 0 (µ) with respect to the chemical potential (Eq. ( 30)).On the contrary, the Mott relation in linear regime does involve the derivative of the anomalous Hall conductivity (Eq.( 31)); however, the non-linear analog of the Mott formula shows a direct relation with the non-linear anomalous Hall conductivity χ 0 (µ) (Eq.( 32)).

IV. MINIMAL MODEL HAMILTONIAN FOR TMDCS
Based on the derivations in last section, we expect a nonzero net Berry curvature or modulated Berry curvature on the Fermi surface to realize the non-linear anomalous thermal Hall effect in the TR invariant (Ω k = −Ω −k ) systems.Therefore, the inversion symmetry has to be broken (Ω k = Ω −k ) to get a nontrivial Berry curvature.Here we consider the non-centro-symmetric 2D materials TMDCs, which have been extensively studied in the past years because of their many fascinating properties [61][62][63][64][65] .The inversion symmetry is broken when the thickness of the 2D TMDCs is reduced to a monolayer.The monolayer TMDCs essentially show the same physics due to their similar atomic structure.Here, we apply the model Hamiltonian of monolayer TMDCs to demonstrate the non-linear anomalous thermal Hall effect, the Wiedemann-Franz law and the Mott relation in the nonlinear regime.
We consider a model Hamiltonian of tilted 2D Dirac cones, which captures the low energy properties of the strained transition metal dichalcogenides monolayer.The corresponding model Hamiltonian can be written as Here, v F is the Fermi velocity, ∆ is the energy band gap opened at the ±K-valley, α is the tilting parameter and τ x,y,z,0 represent Pauli matrices.The wave vector k is measured from the valley center ±K with index s = ±1 (which also indicates the opposite chirality of the Dirac Fermions).
Note that the two massive Dirac cones H s=±1 are mapped to each other by the TR symmetry.Due to the non-zero tilting along k y -axes, only the single mirror plane M x that takes k x → −k x is preserved, in which case, the non-linear effects are allowed to exist 17,66 .
The low energy dispersion and the corresponding Berry curvature of the above Hamiltonian are respectively given as 34) because of the time reversal symmetry.Note the tilting term does not impact the Berry curvature, though it tilts the energy bands via the term ±αk y in E k .Because of the mismatch of this tilting effect, a non-zero net Berry curvature is possible on the Fermi surface.
With the model Hamiltonian described in this section and the non-linear transport coefficients derived in the previous sections, next we calculate the non-linear anomalous thermal Hall effect for a single layer MoS 2 .Based on these calculations, we discuss the Wiedemann-Franz law and the Mott relation in the non-linear regime.

V. RESULTS AND DISCUSSIONS
In this section, we would calculate the non-linear anomalous thermal Hall coefficient for single layer MoS 2 , with the parameters taken from the experiment systems 62,67,68 .We consider only the contribution to the non-linear anomalous thermal Hall effect from the conduction band.Therefore, the spin-orbit coupling effect could be neglected since it is small for the conduction band.
Based on Eq. ( 34), a colormap plot in k x -k y momentum space of Ω k is shown in Fig. 2(a), where the colors indicate the magnitude of the Berry curvature, and a Fermi surface at constant energy µ = 1.5∆ is plotted as the black dash line.Now, the Berry curvature possesses a uniform azimuthal distribution around the valley center ±K, while the Fermi sur- face at a constant energy µ is shifted along the tilt axes (k yaxes).Comparing with the Berry curvature dipole D bd given by Eq. ( 16), we denote the term β 3 (E k −µ) 3 Ω k in Eq. ( 14) the modulated Berry curvature for the non-linear anomalous thermal Hall effect, and its projection in k x -k y momentum space is shown in Fig. 2(b).Because of the factor β 3 (E k − µ) 3 , the modulated Berry curvature near the Fermi energy (black dash) is greatly reduced and a downwards shift of the Berry curvature is observed.Unlike the Berry curvature (Fig. 2(a)), the modulated Berry curvature (Fig. 2(b)) is not symmetric along k y -direction due to the tilting, with only the x-direction mirror symmetry (M x ) that takes k x → −k x preserved.Based on Eqs. ( 16) and ( 22), we could also define another modulated Berry curvature for the non-linear anomalous Nernst effect as β 2 (E k − µ) 2 Ω k , which could be discussed similarly as above.
As has been pointed out earlier in this paper, the non-linear anomalous coefficients namely χ abc , α bd and l cd , are Fermi surface quantities, due to the δ-function-like term (∂f k /∂k b ) in their expressions.As is shown in Fig. (2), the derivative of the Fermi distribution function f k with respect to k x (in (c)) and k y (in (d)) at Fermi energy µ = 1.5∆ is a ring in k x -k y space for the minimal model (Eq.( 34)).The magnitudes of  38) and (39).The constant 7π 4 /15 is included into G (µ).The inset is a zoom-in of the plot around µ = ∆ (black dash line), where we could see that, the numerical data points are consistent with the red line (i.e., the analytical result) at low temperatures (T = 5K, 20K, 50K) for all chemical potentials µ > ∆; for higher temperatures (T = 100K ∼ 300K), the numerical and analytical results differ from each other because of the higher order (O(T 4 )) temperature dependence ignored in the analytical calculations (Eq.( 39)).Here the unit for the y-axes is τ k 2 B / 2 , the other parameters are the same as in Fig.
(  19)) through G 2 (µ) in Eq. ( 24) and G 3 (µ) in Eq. ( 27), respectively.Next, we show the analytical expressions for the non-linear anomalous coefficients based on the model Hamiltonian given in Eq. (33).We find that the analytical results based on Sommerfeld expansion in Sec.(III) are well-matched with the numerical results based on Sec.(II B).
To derive the analytical expressions for the non-linear anomalous Hall coefficient χ abc , the non-linear anomalous Nernst coefficient α ab and the non-linear anomalous thermal Hall coefficient l cd in the low temperature regime, we first show the analytical expression for the Berry curvature dipole.Based on Eqs. ( 16)∼( 21), the zero temperature Berry curvature dipole G(µ) for a two-band system with Ω k and E k given in Eq. ( 34) is computed as 17,66 , Here n 0 = 3vα 0 ∆/4π, n 1 = (1 + α 2 0 )/(1 − α 2 0 ) with α 0 = α/v describing the tilting strength, and µ > 0 for the conduction band is considered.The result based on Eq. ( 35) is valid for chemical potentials µ 2 ≥ ∆ 2 (1 − α 2 0 ) in a tilted system with α 0 .This holds for the non-linear anomalous Nernst coefficient α cd but not for the non-linear anomalous thermal coefficient l cd , since the leading term for l cd is temperature dependent.
As shown in the Sommerfeld expansion for D bd in Eq. ( 20), the temperature-dependent terms (G (2) (µ), and the higher order G (n) (µ)) are relatively small in comparison to the first term G(µ).Therefore, we will not consider the temperature dependence of the Berry curvature dipole.The non-linear anomalous Hall coefficient given in Eq. ( 21) is computed as, Based on the derivations in Sec.(III), especially Eq. ( 25), the analytical expression for the non-linear anomalous Nernst coefficient α cd could be written as, ) where we have ignored the temperature-dependent terms (namely G (n) (µ) with n ≥ 2) as well.
In contrast to χ abc and α cd , we are interested in the temperature dependence of the non-linear thermal Hall coefficient l cd .For small tilting strength i.e., α 0 1 and α µ, the first order derivative with respect to µ of Eq. ( 35) could be approximately written as Substituting the above equation into Eq.( 28), the analytical expression for non-linear thermal Hall coefficient l cd is found as, Through the above equation, we could find the chemical potential dependence of l cd as µ −5 (∆ 2 − µ 2 /2), as well as a (k B T ) 2 temperature dependence for the non-linear anomalous thermal Hall coefficient.It has been shown in Ref. [30], that the non-linear anomalous Nernst coefficient shows a similar chemical potential dependence as that of the non-linear anomalous Hall coefficient in Ref. [24].Now it is clearly seen the underlying reason for it from Eqs. ( 36) and (37).The comparison between the nonlinear anomalous Hall coefficient and the non-linear anomalous Nernst coefficient indicates that, our analytical result for the Mott relation given by Eq. ( 32) is valid in the non-linear regime.In the following, we discuss the Wiedemann-Franz law in the non-linear regime given by Eq. (30).
To verify the chemical potential dependence for the nonlinear anomalous thermal Hall coefficient given by Eq. ( 39), in Fig. (3) we show both the numerical results calculated from Eqs. ( 14) and (34), and anayltical results based on the Sommerdfeld expansion in Eq. (39).As is shown in Fig. (3), the analytical results (the red line) based on Sommerfeld expansion are well matched with the numerical results (the rest of the data besides the red line) at low temperatures with T = 5 ∼ 50K.The numerical and analytical results differ from each other at higher temperatures (T = 100 ∼ 300K) near the chemical potential µ = ∆ as is shown by the inset in Fig. (3).This deviation is caused by the omitted higher orders terms in temperature (O(T 4 )) in the analytical expression given by Eq. (39).Because the magnitude of G (µ) (= G (1) (µ)) around µ = ∆ is large while it is much smaller at the higher chemical potentials, as is shown in Fig. (3), the deviations appears around µ = ∆ and disappear for higher chemical potentials.In Fig. (3), we can also see the sign of l cd changing around µ = √ 2∆, which has been indicated by the chemical potential dependence in Eq. (39).
To verify the temperature dependence, we plot the nonlinear anomalous thermal Hall coefficient l cd as a function of (k B T ) 2 in Fig. (4), where the circles are the numerical results at different chemical potentials based on Eq. ( 14), and the black lines attached to these circles are the corresponding theoretical predictions based on Eq. (39).At low temperatures (T ≤ 50K), the numerical and analytical results are consistent with each other at different chemical potentials, while the numerical result with µ = 1.05∆ (blue circles) start to deviate from the analytical black line around T = 60K.From Fig. (3)  we could see that, the magnitude of G (µ) at µ = 1.05∆ is the maximum.Thus, the contribution from the omitted terms (∝ O(T 4 )) in Eq. ( 39) for µ = 1.05∆ become more evident when the temperature is increasing, as expected.This result is also consistent with that in Fig. (3).The T 2 temperature dependence of the non-linear anomalous thermal Hall coefficient is robust against all the monotonous modulation of the band gap ∆ such as tuning effect by external field 69 , finite temperature effect due to electron-phonon coupling 70 , doping effect through the mixing of chalcogens in MoX 2 (X=S,Se,or Te) 71 , etc., as well as strength of the tilting parameter due to strain [72][73][74][75] .

VI. CONCLUSION
In conclusion, we study the non-linear anomalous thermal Hall effect of the time reversal invariant but inversion symmetry broken systems, which is attributed to a combined effect from the modulated Berry curvature (β 3 (ε k − µ) 3 Ω k ) and the non-equilibrium Fermi distribution (f k ) purturbed by a thermal gradient (−∇T ).We have systematically derived the transverse non-linear thermal current as a second order response to the longitudinal thermal gradient using the Boltzmann semiclassical approach with a relaxation time approximation.The non-linear anomalous thermal Hall coefficient l cd is a pseudotensorial quantity, which plays a similar role as the non-linear anomalous Hall coefficient and non-linear anomalous Nersnt coefficient respectively 17,24,29,30 .By Sommerfeld expansion in low temperature regime, we analytically found the non-linear anomalous thermal Hall coefficient, which is shown to have an explicit relation with the zero temperature Berry curvature dipole and a T 2 temperature dependence.We also found the analog of the Wiedemann-Franz law and the Mott relation relating the anomalous transverse transport coefficients in the non-linear regime.An important byproduct of these calculations is the realization that the non-linear anomalous Nernst coefficient remains non-zero even in the limit of zero temperature.We then apply our theoretical derivations to the monolayer MoS 2 , a TR invariant but inversion symmetry broken TMDC that has been intensively studied recently.Our predictions for the analog of the Wiedeman-Franz law and the Mott formula relating the anomalous transport coefficients in the non-linear regime can be experimentally verified on inversion asymmetric but time reversal symmetric topological systems where the linear anomalous transport coefficients are identically zero by symmetry.In addition, we make specific predictions for the non-linear anomalous themal Hall coefficient in monolayer TMDCs that can be probed in experiments.

FIG. 1 .
FIG. 1. (Color online) Schematic experimental setup for measuring the non-linear thermal Hall effect.A transverse thermal current (∆T ) could be measured as a second-order response of the thermal gradient (−∇T ) even in the presence of time-reversal symmetry.Here we consider the thermal Hall contribution from the electrons (e − ) only.

FIG. 2 .
FIG. 2. (Color online) (a) Berry curvature Ω n=1,s=1 k and (b) modified Berry curvature β 3 (E k − µ) 3 Ω n=1,s=1 k projected in k space for Hamiltonian given in Eq. (33) respectively.Here, panel (a) is azimuthally symmetric whereas (b) is only symmetric with respect to kx.The black dash lines indicate the Fermi surface at µ = 1.5∆.The center of the Fermi surface is shifted upwards along ky and the factor β 3 (E k − µ) 3 greatly reduces the weight of the Berry curvature near the Fermi surface.Panel (c) and (d) show the derivative of Fermi distribution function at Fermi energy µ = 1.5∆ for ∂xf k and ∂yf k respectively.For a massive Dirac cone tilted along y-axis, ∂xf k is symmetrically opposite (red and blue in (c)) along kx on the Fermi surface, while the value of ∂yf k is either positive (blue in (d)) or negative (blue in (d)) along kx on the Fermi surface.The parameters used here are t = 1.1eV, a = 3.19 Å, v = at, α = 0.1v, ∆ = 1.8eV , kx,y ∈ [−0.5π, 0.5π], β = 1/eV is considered for (a), (b) and temperature T = 100K is applied for (c), (d).

FIG. 3 .
FIG. 3. (Color online) Plot of l cd /(kBT ) 2 versus chemical potential µ at different temperatures T .All the data at different temperatures are from numerical calculations, while the red line represents the analytical result with G (µ) given in Eqs.(38) and(39).The constant 7π 4 /15 is included into G (µ).The inset is a zoom-in of the plot around µ = ∆ (black dash line), where we could see that, the numerical data points are consistent with the red line (i.e., the analytical result) at low temperatures (T = 5K, 20K, 50K) for all chemical potentials µ > ∆; for higher temperatures (T = 100K ∼ 300K), the numerical and analytical results differ from each other because of the higher order (O(T 4 )) temperature dependence ignored in the analytical calculations (Eq.(39)).Here the unit for the y-axes is τ k 2 B / 2 , the other parameters are the same as in Fig.(2).
FIG. 3. (Color online) Plot of l cd /(kBT ) 2 versus chemical potential µ at different temperatures T .All the data at different temperatures are from numerical calculations, while the red line represents the analytical result with G (µ) given in Eqs.(38) and(39).The constant 7π 4 /15 is included into G (µ).The inset is a zoom-in of the plot around µ = ∆ (black dash line), where we could see that, the numerical data points are consistent with the red line (i.e., the analytical result) at low temperatures (T = 5K, 20K, 50K) for all chemical potentials µ > ∆; for higher temperatures (T = 100K ∼ 300K), the numerical and analytical results differ from each other because of the higher order (O(T 4 )) temperature dependence ignored in the analytical calculations (Eq.(39)).Here the unit for the y-axes is τ k 2 B / 2 , the other parameters are the same as in Fig.(2).

FIG. 4 .
FIG. 4. (Color online) Similar as Fig. (3) but is a plot of l cd versus (kBT ) 2 with different chemical potentials.The circles are from numerical calculations, while the black lines corresponding to each chemical potential are the analytical results based on Eq. (39).Here the units for y-axes is τ k 2 B / 2 , the applied temperature T ∈ [5K, 100K] with a unit step (1K), and all other parameters are the same as in Fig. (2).