Rectifying full-counting statistics in a spin Seebeck engine

TL;DR

This paper explores the rectification and negative differential effects in the full-counting statistics of thermal spin transport within a spin Seebeck engine, revealing quantum fluctuation effects at the electron/magnon interface.

Contribution

It formulates the scaled cumulant generating function for the system using nonequilibrium Green's functions, highlighting the role of quantum dot-induced fluctuations.

Findings

Identification of rectification effects in spin transport.

Observation of negative differential conductance phenomena.

Quantum fluctuations significantly influence transport properties.

Abstract

We investigate the rectification and negative differential effects of full-counting statistics of conjugate thermal spin transport in a spin Seebeck engine. The engine is made by an electron/magnon interface diode driven by a temperature bias, of which the scaled cumulant generating function is formulated in the framework of nonequilibrium Green's function. The strongly fluctuating interfacial electron density of states induced by quantum dot is responsible for these intriguing properties. This work is relevant for designing efficient spin caloritronic engine and can broaden our view in nonequilibrium thermodynamics and the nonlinear phenomenon in quantum transport system.

Equations

$$\hat{H}={\hat{H}}_{el}+{\hat{H}}_{mag}+{\hat{H}}_{em},$$

$${\hat{H}}_{el}=\sum _{k\sigma }{ϵ}_{k\sigma }{a}_{k\sigma }^{†}{a}_{k\sigma }.$$

$$f_{L\sigma}(\epsilon_{k\sigma})=\big{\{}\exp[\beta_{L}(\epsilon_{k\sigma}-\mu_{\sigma})]+1\big{\}}^{-1},$$

$${\hat{H}}_{mag}=-J\sum _{⟨i,j⟩}\left[\frac{1}{2}{S}_{i+}{S}_{j-}+\frac{1}{2}{S}_{i-}{S}_{j+}+S_{iz}{S}_{jz}\right],$$

$${\hat{H}}_{mag}\approx \sum _q\hslash {\omega }_q{b}_q^{†}{b}_q+\mathrm{const},$$

$${\hat{H}}_{em}=-\sum _{k{k}^{\prime }q}{g}_q\left[{\gamma }_{k\uparrow }^{\ast }{\gamma }_{k^{\prime }\downarrow }{a}_{k\uparrow }^{†}{a}_{k^{\prime }\downarrow }{b}_q^{†}+H.c.\right]{\delta }_{k-k^{\prime }=q},$$

$$\mathcal{G}({\lambda }_s,$$

$$\displaystyle\Big{[}$$

$$+$$

$$A(E_{-},E_{+})=\frac{J_{R}(\omega)J_{L\uparrow}(E_{-})J_{L\downarrow}(E_{+})}{\big{|}[\Sigma_{R}^{r}(\omega)+F^{r}(\omega)\big{|}^{2}}.$$

$$F^r(\omega )=i\int \frac{d{\omega }_1}{2\pi }\frac{F^{>}({\omega }_1)-F^{<}({\omega }_1)}{\omega -{\omega }_1+i{0}^{+}}.$$

$$F^{<}(\omega )$$

$$F^{>}(\omega )$$

$$\mathcal{G}({\lambda }_s,{\lambda }_h)=\mathcal{G}(-{\lambda }_s-i{\beta }_L\Delta {\mu }_s,-{\lambda }_h+i({\beta }_L-{\beta }_R)).$$

$$\displaystyle I_{s(h)}=\int\frac{d\omega}{2\pi}\int\frac{dE}{2\pi}\Big{\{}(\hbar\omega)^{\alpha}A(E_{-},E_{+})\times$$

$$\displaystyle\left[f_{\downarrow\rightarrow\uparrow}(E_{+},E-)(N_{R}(\omega)+1)-f_{\uparrow\rightarrow\downarrow}(E_{-},E_{+})N_{R}(\omega)\right]\Big{\}},$$

$$\hat{H}_{el}=\sum_{\sigma}\epsilon_{\sigma}d_{\sigma}^{\dagger}d_{\sigma}+\sum_{k\sigma}\Big{[}\epsilon_{k\sigma}c_{k\sigma}^{\dagger}c_{k\sigma}+(t_{k\sigma}c_{k\sigma}^{\dagger}d_{\sigma}+{\rm H.c.})\Big{]},$$

$$d_{\sigma }=\sum _k{\gamma }_{k\sigma }{a}_{k\sigma },c_{k\sigma }=\sum _{k^{\prime }}{\eta }_{k{k}^{\prime }\sigma }{a}_{k^{\prime }\sigma }$$

$${\gamma }_{k\sigma }=\frac{t_{k\sigma }}{{ϵ}_{k\sigma }-{ϵ}_{\sigma }-{\Sigma }_{k\sigma }},{\eta }_{k{k}^{\prime }\sigma }={\delta }_{k{k}^{\prime }}-\frac{t_{k\sigma }{\gamma }_{k^{\prime }\sigma }}{{ϵ}_{k\sigma }-{ϵ}_{k^{\prime }\sigma }+i{0}^{+}},$$

$$J_{L\sigma }(E)=2\pi g\sum _k|{\gamma }_{k\sigma }{|}^2\delta (E-{ϵ}_{k\sigma })=\frac{g{\Gamma }_{\sigma }}{(E-{ϵ}_{\sigma }{)}^2+{\Gamma }_{\sigma }^2/4},$$

Full text

Introduction– Spin caloritronics has emerged as a new field to utilize the excess heat generated in nano-devices to drive spin currents Bauer2 . Of particular interest is the spin Seebeck effect (SSE) in magnetic insulators, which generates pure spin transfer out of thermal gradient without conducting electron currents, so that opens a new direction to harvest heat in the absence of Joule heating SSE ; SSE1 ; SSE2 ; SSE3 ; SSE4 ; Jie1 ; Jie2 ; Ronetti . The pure spin is carried by magnons (quantized spin waves) that can transfer energy and spin angular momentum. Besides, magnons have long life-time so that they can propagate over long distances in a ballistic way without being dissipated in some materials due to the low damping YIG1 ; YIG2 ; YIG3 .

Many proposals have been made in exploiting the application of magnons, such as energy harvester harvester1 , diodes Jie2 ; diode , transistors transistor and logic deviceslogic . In particular, the two terminal hybrid system involving both metallic and magnonic reservoir has been proposed Jie1 ; Barnas . In such system, the propagation of spin waves can be converted to spin current in the metallic part and vice versa across the junction under a temperature bias. It has been pointed out that strongly fluctuating electron density of states (DOS) in the metallic materials is helpful to obtain the nontrivial rectification and negative differential effects Jie1 ; Jie2 , which suggests that quantum dot (QD) structure is a good candidate in achieving these effects. The reduced dimension of QD can also help to increase the figure of merit (FOM) and harvest from the waste heat more efficiently compared to bulk systems when working as a heat engine QD2 ; QD3 ; QD4 .

So far, the rectification and negative differential effect are only studied in the context of mean values of heat and spin currents. For a nanoscale device, fluctuations and even higher-order cumulants of currents can be significant and play an important roles, which could be well described within the framework of full-counting statistics (FCS)Levitov2 ; Levitov3 ; RMP ; Hanggi ; Ka2 ; gm2 ; gm3 ; gm4 . FCS of a physical quantity is encoded in its cumulant generating function (CGF), which reveals us the higher-order Onsager reciprocity relations Hanggi ; Gaspard1 and the thermodynamic fluctuation symmetry RMP ; Hanggi ; Gaspard2 as well. It would illuminate us in better understanding the nonequilibrium thermodynamics and designing efficient heat harvester to study the rectification of higher-order cumulants, fluctuation symmetry, and heat engine performances for the two terminal hybrid system involving both electronic and bosonic reservoirs.

In this work, we investigate the rectification and negative differential effects of FCS in a spin Seebeck engine sandwiched by one electronic and the other magnonic reservoirs. We have formulated a framework to obtain scaled cumulant generating function (SCGF) of both spin and heat currents using nonequilibrium Green’s function (NEGF), which does not require week couplings to both reservoirs, and treats the interfacial electron-magnon interaction to second order. The fluctuation symmetry is obtained and the condition in realizing rectification and negative differential effect is discussed. The rectification and negative differential effects of cumulants and heat engine performances with respect to temperature reversal are studied in details.

Model and theoretical formalism– The system is schematically illustrated in Fig. 1 wherein the left electronic compartment interacts with the right magnonic reservoir. The spin current can be generated out of a temperature gradient $\Delta T=T_{L}-T_{R}$. The whole system Hamiltonian has contributions from the electrons, magnons and the interaction between them,

[TABLE]

The electron Hamiltonian $H_{el}$ is expressed as:

[TABLE]

$\hat{H}_{el}$ can incorporate the non-interacting electronic scattering region as well, and we take a single two-level QD in the numerical calculation. The electronic spin-$\sigma$ chemical potential is $\mu_{\sigma}$ and their difference is defined as the spin bias with the form $\Delta\mu_{s}=\mu_{\downarrow}-\mu_{\uparrow}$ which could be measured by inverse spin Hall effect SSE2 . The Fermi distribution for the spin-$\sigma$ is

[TABLE]

with inverse temperature $\beta_{L}=1/(k_{B}T_{L})$.

The right insulating magnetic compartment can be described by a Heisenberg lattice

[TABLE]

where $S_{j\pm}$ is the localized raising (lowering) spin operator at site $j$, $S_{jz}$ is the spin operator in the $z$ direction, and $J$ denotes the exchange coupling strength. Using Holstein-Primakoff transformation HP ; Mahan , $S_{j+}=\sqrt{2S_{0}-b_{j}^{\dagger}b_{j}}b_{j},\ S_{j-}=b_{j}^{\dagger}\sqrt{2S_{0}-b_{j}^{\dagger}b_{j}},\ S_{jz}=S_{0}-b_{j}^{\dagger}b_{j}$ with localized spin length $S_{0}$, the spin operators are mapped into bosonic magnons. At large spin limit Jie1 ; harvester1 ; QD4 , i.e., $2S_{0}\gg\langle b_{j}^{\dagger}b_{j}\rangle$, we have the approximation $S_{j}^{-}\approx\sqrt{2S_{0}}b_{j}^{\dagger}$ and $S_{j}^{+}\approx\sqrt{2S_{0}}b_{j}$. The right magnetic insulator can be approximated by the free magnon gas Jie1 :

[TABLE]

after a Fourier transform into the momentum space, where the dispersion of $\omega_{q}$ depends on the material details. The magnonic reservoir obeys the Bose-Einstein distribution $N_{R}(\omega)=1/(e^{\beta_{R}\omega}-1)$, where the chemical potential of magnons is set to zero due to the ferromagnetic insulator phase considered here Magnetic_RC_Circuit . The interfacial electron-magnon interaction Hamiltonian bears the form Jie1 ; Mahan ,

[TABLE]

with the delta function $\delta_{k-k^{\prime}=q}$ ensuring energy and spin angular momentum conservation. The interaction describes magnon assisted spin flip by absorbing or emitting magnons with frequency $\omega_{q}$. The effect from the electronic DOS is incorporated in $\gamma_{k\sigma}$. The electron-magnon interaction strength $g_{q}$ is assumed to be weak and will be treated perturbatively in getting SCGF. For simplicity, one can assume $g_{q}$ is energy independent with $g_{q}=g$.

Spin (s) and heat (h) current operators are defined as the change rate of the electronic spin in the left compartment and heat in the right reservoir, respectively, with the forms $\hat{I}_{s}(t)=-d_{t}\hat{N}_{s}$, $\hat{I}_{h}(t)=-d_{t}\hat{H}_{R}$, and $d_{t}$ being the total differential with respect to time in Heisenberg picture. Here, $\hat{N}_{s}=\frac{1}{2}\sum_{k}(a_{k\downarrow}^{\dagger}a_{k\downarrow}-a_{k\uparrow}^{\dagger}a_{k\uparrow})$ is the operator of number of spin in the left compartment and $\hat{H}_{R}=\sum_{q}\hbar\omega_{q}b_{q}^{\dagger}b_{q}$. Using random phase approximation (RPA), that is treating the interaction strength $g$ between the electrons and magnon modes perturbatively to the second order Ka2 , we could obtain SCGF for spin and heat current (with counting fields $\lambda_{s}$ and $\lambda_{h}$) expressed as (See Supplemental Material supp for details of derivation),

[TABLE]

with

[TABLE]

Here we have defined $f_{\uparrow\rightarrow\downarrow}(E_{-},E_{+})=f_{L\uparrow}(E_{-})[1-f_{L\downarrow}(E_{+})]$ and $f_{\downarrow\rightarrow\uparrow}(E_{+},E_{-})=f_{L\downarrow}(E_{+})[1-f_{L\uparrow}(E_{-})]$ with $E_{\pm}=E\pm\omega/2$. $J_{L\sigma}(E)=2\pi g\sum_{k}|\gamma_{k\sigma}|^{2}\delta(E-\epsilon_{k\sigma})$ is the the spin-$\sigma$ electronic spectral density. The retarded magnonic self-energy due to the dissipation is expressed as $\Sigma_{R}^{r}=-iJ_{R}(\omega)/2$. Without loss of generality, the magnonic reservoir spectral function is considered to be Ohmic and expressed as $J_{R}(\omega)=\pi\alpha\omega e^{-\omega/\omega_{c}}$, where $\alpha$ is the dimensionless dissipation strength and $\omega_{c}$ is the cut-off frequency Weiss ; supp . The expression of retarded electron-hole propagator is,

[TABLE]

with

[TABLE]

The physical picture of SCGF is clear. The first contribution of Eq. (Rectifying full-counting statistics in a spin Seebeck engine) describes that the spin-down electrons with energy $E+\hbar\omega/2$ flip to the spin-up states with energy $E-\hbar\omega/2$ by emitting a magnon with energy $\hbar\omega$ into the magnonic reservoir. The second one describes a reverse process that the spin-up electrons with energy $E-\hbar\omega/2$ flip to the spin-down states with energy $E+\hbar\omega/2$ via absorbing a magnon carrying energy $\hbar\omega$.

Applying the identity $f_{\uparrow\rightarrow\downarrow}(E_{-},E_{+})=e^{\beta_{L}(\omega-\Delta\mu_{s})}f_{\downarrow\rightarrow\uparrow}(E_{+},E_{-})$ in Eq. (Rectifying full-counting statistics in a spin Seebeck engine), we can obtain the following fluctuation symmetry relation RMP ; Hanggi ,

[TABLE]

If $J_{L\sigma}(E)$ were flat near Fermi energy, $A(E_{-},E_{+})$ inside the integral in Eq. (Rectifying full-counting statistics in a spin Seebeck engine) would be treated as a constant. By applying the equality $\int\frac{dE}{2\pi}f_{\downarrow\rightarrow\uparrow}(E_{+},E_{-})=(\hbar\omega-\Delta\mu_{s})N_{L}(\hbar\omega-\Delta\mu_{s})$ and assuming $\Delta\mu_{s}\rightarrow 0$, the SCGF would be almost symmetric by reversing the temperature gradient $\Delta T$ with a constant electron DOS and so are the cumulants of currents, that is, ${\cal G}(\lambda_{s},\lambda_{h})={\cal G}(-\lambda_{s},-\lambda_{h};L\leftrightarrow R)$. This is usually satisfied in the pure fermionic system or bosonic system, but is not satisfied in the electron-magnon hybrid system discussed here. Therefore, strongly fluctuating electron DOS is essential to have asymmetric behavior with respect to reversing the temperature gradient and rectify FCS and heat engine performances of the spin Seebeck engine. The spin current in the electronic reservoir and heat current in the magnonic reservoir could be derived by taking the first order derivative of ${\cal G}(\{\lambda\})$ with respect to $\lambda_{s}$ and $\lambda_{h}$, respectively at $\lambda_{s}=\lambda_{h}=0$, that is, $I_{s/q}=\partial{\cal G}(\{\lambda\})/\partial(i\lambda_{s/q})|_{\lambda_{s}=\lambda_{h}=0}$. Then,

[TABLE]

where $\alpha=0$ for spin current $I_{s}$, and $\alpha=1$ for heat current $I_{h}$.

We now discuss the device working as a thermoelectric engine with the output spin power $P=|I_{s}\Delta\mu_{s}|$. When $\Delta T=T_{L}-T_{R}<0$, heat current flows from the right reservoir to drive the spin current, the heat current calculated from Eq. (Rectifying full-counting statistics in a spin Seebeck engine) is negative and the efficiency should be defined as $\eta=P/|I_{h}|$. One should note that the spin bias which is defined as $\Delta\mu_{s}=\mu_{\downarrow}-\mu_{\uparrow}$ should be chosen as positive $\Delta\mu_{s}>0$ under $\Delta T<0$, otherwise the device works as a unphysical dud engine where output spin power heats both reservoirs. In the case of $\Delta T>0$, i.e., $T_{L}>T_{R}$, the left reservoir is the heat source with spin bias $\Delta\mu_{s}<0$ and heat current is $I_{Lh}=|I_{h}|+P$ using energy conservation law, so that the efficiency has the definition of the form $\eta=P/(|I_{h}|+P)$.

Numerical results– In order to get a fluctuating electron DOS, one can, for example, insert a QD between the electronic and magnonic reservoir. Hamiltonian of the electronic part is then expressed as

[TABLE]

with the first term describing QD, the second term electronic reservoir and the rest for their coupling. $\epsilon_{\sigma}$ are the two energy levels inside QD. We consider a large quantum dot so that the Coulomb interaction effect can be neglected. $\epsilon_{k\sigma}$ and $t_{k\sigma}$ are the energy level of state $k\sigma$ in the electronic reservoir and its coupling strength with the spin-$\sigma$ in QD. The electron Hamiltonian $\hat{H}_{el}$ with QD can be diagonalized and rewritten in the form of Eq. (2) through the linear transformation Bijay ; Segal

[TABLE]

where the dimensionless coefficients are

[TABLE]

with the self-energy function $\Sigma_{k\sigma}=\sum_{k^{\prime}}t_{k^{\prime}\sigma}^{2}/(\epsilon_{k}-\epsilon_{k^{\prime}}+i0^{+})$. The electronic spectral densities in the left compartment (electronic reservoir plus QD) is Lorentzian shaped with the form Bijay

[TABLE]

where the electronic coupling strength $\Gamma_{\sigma}=2\pi\sum_{k}|t_{k\sigma}|^{2}\delta(E-\epsilon_{k\sigma})$ is assumed energy independent.

During the numerical calculation, the coupling strengths for the spin up and down electrons between the electronic reservoir and QD are assumed to be energy independent and equal, i.e., $\Gamma_{\uparrow}=\Gamma_{\downarrow}=\Gamma$. The energy level of spin up is set to be larger than the that of spin down in the QD, $\epsilon_{\downarrow}>\epsilon_{\uparrow}$. The dimensionless dissipation strength and the cut-off frequency of the magnonic reservoir is chosen as $\alpha=0.2$, and $\omega_{c}=20meV$, respectively. The temperatures of the reservoirs are set as $T_{L,R}=T_{0}\pm\Delta T/2$ with $T_{0}=300K$, so that when the normalized temperature gradient $\Delta T/(2T_{0})<0$, the left lead is cooler. We denote the average of spin up and down chemical potential as $\mu_{0}=(\mu_{\uparrow}+\mu_{\downarrow})/2$. The spin currents are displayed in atomic unit.

In Fig. 2, we plot the spin current $I_{s}$ [panel (a)], spin current noise $\langle\langle(I_{s})^{2}\rangle\rangle$ [panel (b)], and Fano factor $F_{2}=\langle\langle(I_{s})^{2}\rangle\rangle/I_{s}$ [panel (c)] as a function of the normalized temperature gradient $\Delta T/(2T_{0})$ by varying $\epsilon_{\sigma}$. The spin up and down chemical potentials are both set to zero, i.e., $\mu_{\uparrow}=\mu_{\downarrow}=0$. The rectification effect where the spin current are asymmetric with respect to temperature gradient reversal can be clearly identified in Fig. 2. In the cases where the two levels of QD are either both above or below the chemical potential ($\epsilon_{\downarrow}>\epsilon_{\uparrow}>0$ in case (i) or $0>\epsilon_{\downarrow}>\epsilon_{\uparrow}$ in case (iv) in Fig. 2), we could observe the negative differential SSE that the spin current decrease or even vanish with the increasing temperature gradient $\Delta T$. This could be understood as follows: when the right reservoir is hotter, i.e., $\Delta T<0$, the process that spin up electron absorbs energy to flip to the spin down state dominates. Near the linear response regime, increasing the magnitude of temperature gradient which corresponds to decreasing $T_{L}$ leads to increased heat current to assist the spin flip process. When the two levels are both above (below) the chemical potential, the electrons on these levels can be depleted (occupied) by decreasing $T_{L}$ further. Compared to the increasing $N_{R}(\omega)$, $f_{\uparrow\rightarrow\downarrow}(E_{-},E_{+})=f_{L\uparrow}(E_{-})[1-f_{L\downarrow}(E_{+})]$ is severely suppressed and so is the flipping process from spin up to down state Jie1 . Thus, the essential ingredient to get negative differential SSE is that the two levels are on one side of the chemical potential. We can also observe from Fig. 2(b) that spin current noise is asymmetric with respect to temperature gradient which signifies rectification as well. If $J_{L\sigma}(E)$ is flat near $\epsilon_{\sigma}$, spin current noise is an even function of $\Delta T$, which is not satisfied in Fig. 2(b). For case (ii) and (iii), noise are monotonic, while for case (i) and (iv), decrease of noise with increasing the magnitude of $\Delta T$ ($\Delta T/(2T_{0})<-0.5$) is due to the severely suppressed number of electrons involving in spin flip. One can observe that the Fano factors shown in Fig. 2(c) are relatively symmetric and almost the same for all cases, except in the region $\Delta T/(2T_{0})\rightarrow-1$.

In Fig. 2(d), spin current $I_{s}$ versus $\Delta T$ by varying $\epsilon_{\uparrow}$, which is above the chemical potential, is plotted. With increasing $\epsilon_{\uparrow}$ away from the electronic chemical potential, negative differential points shift towards $\Delta T=0$, while negative differential and asymmetric behaviors become less obvious. This is due to the fact that the number of electrons near $\epsilon_{\sigma}$ for spin flip is becoming less with increasing $\epsilon_{\uparrow}$. We conclude that the differential points can be adjusted by tuning QD levels.

The spin currents and spin current noises versus $\Delta T/(2T_{0})$ by varying the electronic coupling strength $\Gamma$ with $\epsilon_{\uparrow}=30meV$ and $\epsilon_{\downarrow}=35meV$ are shown in Fig. 3(a) and (b), respectively. With increasing $\Gamma$, the spin current increases first and then decreases. Since when $\Gamma$ is very small, $\Gamma\approx 0.2meV$, the electronic spectral density is quite small so that the spin flip is limited by $\Gamma$ and the spin current will increase with increasing $\Gamma$. However when $\Gamma$ is large, increasing $\Gamma$ reduces the spin current because of the level broadening. The negative differential SSE disappears once $\Gamma$ is too large, since electron DOS becomes flatter near $\epsilon_{\sigma}$ with increasing $\Gamma$. Spin current noise is getting symmetric with increasing $\Gamma$ due to electron DOS flattening, this holds as well for the case of $\epsilon_{\downarrow}>0>\epsilon_{\uparrow}$, which is not shown here.

For a heat engine, maximum output power and maximum efficiency cannot be satisfied at the same time. People usually care about the efficiency $\eta_{\rm MP}$ when the system has a maximum power $P_{\rm MP}$ or the output power under a maximum efficiency at a given $\Delta T$ MP2 ; MP3 ; PED . In Fig. 4, we plot the maximum output spin power $P_{\rm MP}$, and efficiency $\eta_{\rm MP}$ at maximum power as a function of $\Delta T/(2T_{0})$ by varying the average chemical potential $\mu_{0}$. We observe that both $P_{\rm MP}$ and $\eta_{\rm MP}$ are asymmetric with respect to the temperature reversal. When $\mu_{0}$ is above or below both levels $\epsilon_{\sigma}$ of the QD, we can observe the negative differential effect which is also due to the severely decreased electron number near $\epsilon_{\sigma}$ with decreasing $T_{L}$.

Conclusion– In this work we have uncovered the rectification and negative differential effects of FCS and heat engine performances in a electron/magnon diode under a temperature gradient. The SCGF of spin and heat flow of the system were obtained in the framework of NEGF. Spin and heat currents expressions were derived and fluctuation symmetry was verified from the SCGF. We verified that the strongly fluctuating electron DOS induced by QD is crucial to get these intriguing effects. In principle, multiple QDs can be sandwiched between metallic and magnonic reservoirs to form multiple transport channels in parallel to enhance transport. Without loss of generality, we focus on the single two-level QD here, as a demonstration of the intriguing properties studied. The hybrid two terminal setup involving both electronic and bosonic reservoirs exhibits a perfect testing ground for some nontrivial phenomena in nonequilibrium thermodynamics.

G. Tang, X. Chen and J. Wang are financially supported by the General Research Fund (Grant No. 17311116), the University Grant Council (Contract No. AoE/P-04/08) of the Government of HKSAR, and NSF-China (Grant No. 11374246). J. Ren acknowledges support from the National Youth 1000 Talents Program in China, and the 985 startup Grant (205020516074) at Tongji University.