Nonintegrability and thermalization of one-dimensional diatomic lattices

TL;DR

This study investigates how breaking integrability in one-dimensional diatomic lattices affects thermalization, revealing a universal inverse-square law for thermalization time regardless of the method used to break integrability.

Contribution

It demonstrates that different methods of destroying integrability lead to distinct thermalization pathways but share a universal inverse-square relationship for thermalization time.

Findings

Thermalization time scales inversely with the square of perturbation strength.

Different methods of breaking integrability produce qualitatively different thermalization routes.

Universal law of thermalization time applies across different 1D diatomic lattice models.

Abstract

Nonintegrability is a necessary condition for the thermalization of a generic Hamiltonian system. In practice, the integrability can be broken in various ways. As illustrating examples, we numerically studied the thermalization behaviors of two types of one-dimensional (1D) diatomic chains in the thermodynamic limit. One chain was the diatomic Toda chain whose nonintegrability was introduced by unequal masses. The other chain was the diatomic Fermi-Pasta-Ulam-Tsingou-$\beta$ chain whose nonintegrability was introduced by quartic nonlinear interaction. We found that these two different methods of destroying the integrability led to qualitatively different routes to thermalization, but the thermalization time, $T_{eq}$, followed the same law; $T_{eq}$ was inversely proportional to the square of the perturbation strength. This law also agreed with the existing results of 1D monatomic lattices. All these results imply that there is a universal law of thermalization that is independent of the method of breaking integrability.

Equations

$$H=\sum _{j=1}^N\left[\frac{p_j^2}{2m_j}+V(q_j-q_{j-1})\right],$$

$$V_T(x)=\frac{e^{2x}-2x-1}{4}.$$

$$V_{\beta }(x)=\frac{x^2}{2}+\frac{\beta x^4}{4},$$

$$H=H_0+H^{\prime },$$

$$H=\sum _{j=1}^N\left[\frac{v_j^2}{2}+\frac{1}{m_j}V(q_j-q_{j-1})\right],$$

$$H_0=\sum _{j=1}^N\left[\frac{v_j^2}{2}+\frac{1}{m_2}{V}_T(q_j-q_{j-1})\right].$$

$$H^{\prime }=\frac{|m_2-m_1|}{m_1{m}_2}\sum _{j=1,3,5,\cdots }^N{V}_T(q_j-q_{j-1}),$$

$$⟨H^{\prime }⟩\propto \frac{|m_1-m_2|}{m_1{m}_2}=\frac{\Delta m}{1-\Delta m^2/4}\sim \Delta m,$$

$$H_0=\sum _{j=1}^N\left[\frac{v_j^2}{2}+\frac{1}{2m_j}(q_j-q_{j-1}{)}^2\right].$$

$$H^{\prime }=\beta \sum _{j=1}^N\frac{(q_j-q_{j-1}{)}^4}{4m_j},$$

$$⟨H^{\prime }⟩\propto \beta \frac{(m_1+m_2)}{m_1{m}_2}=\frac{\beta }{1-\Delta m^2/4}\sim \beta .$$

$$u_{j,k}^{\pm }=\{\begin{array}{ll}\frac{\sin \left(\frac{2lk\pi }{2L+1}\right)+\sin \left(\frac{2(l-1)k\pi }{2L+1}\right)}{\sin \left(\frac{2k\pi }{2L+1}\right)},& j=2l-1;\\ \frac{\left[2-m_1{\left({\omega }_k^{\pm }\right)}^2\right]\sin \left(\frac{2lk\pi }{2L+1}\right)}{\sin \left(\frac{2k\pi }{2L+1}\right)},& j=2l,\end{array}$$

$${\omega }_k^{\pm }=\sqrt{\frac{m_1+m_2}{m_1{m}_2}\left[1\pm \sqrt{1-\frac{4m_1{m}_2}{(m_1+m_2{)}^2}{\sin }^2\left(\frac{k\pi }{2L+1}\right)}\right]}.$$

$${\mathbf{U}}_k^{+}$$

$${\mathbf{U}}_k^{-}$$

$$\{\begin{array}{ll}{Q}_k& ={\sum }_{j=1}^N{q}_j{U}_{j,k},\\ P_k& ={\sum }_{j=1}^N{p}_j/m_j{U}_{j,k},\end{array}k=1,2,\cdots ,N.$$

$$E_k=\frac{1}{2}\left(P_k^2+{\omega }_k^2{Q}_k^2\right),$$

$$Q_k=\sqrt{2E_k/{\omega }_k^2}\sin \left({\phi }_k\right),P_k=\sqrt{2E_k}\cos \left({\phi }_k\right).$$

$$\underset{T\to \infty }{\lim }{\bar{E}}_k(T)\simeq \epsilon ,k=1,\cdots ,N,$$

$${\bar{E}}_k(T)=\frac{1}{(1-\mu )T}{\int }_{\mu T}^T{E}_k(P(t),Q(t))dt.$$

$$\xi (t)=\tilde{\xi }(t)\frac{e^{\eta (t)}}{L},\eta (t)=-\sum _{k>L}^N{w}_k(t)\log [w_k(t)]$$

$$\tilde{\xi }(t)=\frac{{\sum }_{k>L}^N{\bar{E}}_k(t)}{\frac{1}{2}{\sum }_{1\le k\le N}{\bar{E}}_k(t)},w_k(t)=\frac{{\bar{E}}_k(t)}{{\sum }_{j>L}^N{\bar{E}}_j(t)}.$$

$$\xi (t{)}_A=\frac{e^{\eta (t)}}{L},\eta (t)=-\sum _{k=1}^L{w}_k(t)\log [w_k(t)],$$

$$\xi (t{)}_O=\frac{e^{\eta (t)}}{L},\eta (t)=-\sum _{k>L}^N{w}_k(t)\log [w_k(t)],$$

Full text

I Introduction

The famous ergodic hypothesis formulated by Boltzmann in the 1870s is at the foundation of statistical physics Khinchin (1949). In the 1950s, Fermi, in collaboration with Pasta, Ulam, and Tsingou (FPUT), conducted the first numerical experiment to verify this hypothesis by observing the rate of mixing and thermalization in microscopic reversible mechanical systems Fermi et al. (1955); Dauxois (2008). However, the result was contrary to general expectations. The system far from equilibrium did not enter the expected thermalized state, but returned to nearly the initial nonequilibrium state. Such a phenomenon is named the “FPUT paradox” Dauxois (2008). This seminal work failed to observe the expected picture but it spurred many great mathematical and physical discoveries such as integrability Zakharov (1991), soliton physics Zabusky and Kruskal (1965), and deterministic chaos Zaslavsky (2005).

More than half a century has passed; the literature on the subject is too varied to summarize here. The state of the art, updated to a few years ago, can be found in a collection of papers Cha (2005) and a status report Gallavotti (2008). Nevertheless, if one tries to draw from these sources any clear conclusion about the mathematical status of the FPUT problem or the physical meaning of the results, one may remain rather confused Gallavotti (2008). The literature is sometimes a bit confusing for various conditions due to the very rich dynamics of one-dimensional (1D) nonlinear chains. For instance, based on Nekhoroshev’s theory on exponential stability Nekhoroshev (1977), equipartition should always occur, and the thermalization time, $T_{eq}$, has a stretched exponential dependence on the energy density, $\varepsilon$; i.e., $T_{eq}\propto\exp(\varepsilon^{-a})$ with some positive $a$ Fucito et al. (1982); Pettini and Landolfi (1990); Galgani et al. (1992); Berchialla et al. (2004). However, some other researchers have shown that the thermalization time follows a power law; i.e., $T_{eq}\propto\varepsilon^{-a}$ with some positive $a$ DeLuca et al. (1995); *PhysRevE.54.2329; *PhysRevE.60.3781. Moreover, some results have shown a crossover from the stretched exponential law to a power law for various conditions Benettin and Ponno (2011). Recently, the wave turbulence theory Zakharov et al. (1992); Majda et al. (1997); Zakharov et al. (2001, 2004); Nazarenko (2011) was applied to attack this problem Onorato et al. (2015); Lvov and Onorato (2018); Pistone et al. (2018a). It was shown analytically and numerically that the exact nontrivial multi-wave resonant interactions were responsible for the thermalization of short FPUT chains in the weakly nonlinear regime and this resulted in $T_{eq}$ following a power law. It was further conjectured that $T_{eq}$ still followed a power law with different exponents in the thermodynamic limit.

More recently, we have shown via extensive numerical calculations that the thermalization follows a universal power law in the thermodynamic limit Fu et al. (2018). This universal law, $T_{eq}\sim\gamma^{-2}$, applies generally to a class of 1D lattices with interaction potential $V(x)=x^{2}/2+\lambda x^{n}/n$, where $n\geq 4$ is an integer and $\gamma=\lambda\varepsilon^{(n-2)/2}$ is the perturbation strength. This scaling law also applies to another class of 1D lattices with symmetric interaction potential $V(x)=x^{2}/2+\lambda|x|^{d}/d$, where $d=m_{1}/m_{2}>2$ with $m_{1}$ and $m_{2}$ being two coprime integers and the perturbation strength $\gamma=\lambda\varepsilon^{(d-2)/2}$. Furthermore, we numerically confirmed that this law of thermalization also held in the perturbed Toda lattices in the thermodynamic limit; i.e., $T_{eq}\propto\epsilon^{-2}$, where the perturbation strength, $\epsilon$, characterizes the distance between the perturbed potential and the Toda potential Fu et al. (2019). It has been shown that the key to identifying the universal exponent $-2$ is to select a suitable reference integrable system so that the perturbation strength can be defined accurately. We noted that the thermalization in the Klein-Gordon lattice Pistone et al. (2018a, b) also follows this law, though this lattice belongs to another class that possesses on-site potential.

In the above studies, the integrability of the system was broken by introducing nonlinearity. However, there are various ways to destroy the integrability Toda (1989), including introducing impurities and raising dimensionality. For example, the nonintegrability of a 1D diatomic Toda lattice was introduced by unequal masses Casati and Ford (1975). This is a somewhat natural perturbation in the sense that isotopic mass impurities do occur in nature. There are many studies on the diatomic Toda chain Dash and Patnaik (1981); Mokross and Büttner (1981); Diederich (1985); Aoki and Takeno (1995); Mokross and Buttner (1983); Jackson and Mistriotis (1989); Hatano (1999); Chen et al. (2014), especially on its heat conduction behavior Mokross and Buttner (1983); Jackson and Mistriotis (1989); Hatano (1999); Chen et al. (2014), but to our knowledge, the thermalization problem of diatomic Toda chains has not been discussed systematically.

In this study, we examined the thermalization rate of 1D diatomic Toda lattices in order to determine whether a universal law of $T_{eq}$ existed for this class of lattices as well, and if the answer was yes, how it differed from $T_{eq}\propto\epsilon^{-2}$ applicable to the 1D lattices in which the nonintegrability is introduced by nonlinearity. For comparison, the diatomic FPUT-$\beta$ lattice was also studied, which is the diatomic harmonic lattice perturbed by the quartic nonlinearity; that is, its nonintegrability was introduced by nonlinearity, though it was a diatomic chain. In the following sections, we first introduce the models in the Sec. II, then give the definitions of perturbation strength for the two models in Sec. III. The physical quantities and the numerical method are provided in Sec. IV. The numerical results are described and presented in Sec. V, followed by the summary and discussions in Sec. VI.

II The models

We considered a 1D diatomic lattice consisting of $L$ unit cells, i.e., the total number of particles $N=2L$. Each cell contained two particles of mass $m_{1}$ and $m_{2}$ ($m_{1}<m_{2}$) situated alternately at the position $2l-1$ and $2l$ in the $l$th unit cell. The physical configuration is illustrated in Fig. 1. Its Hamiltonian is

[TABLE]

where $p_{j}$, $m_{j}$, and $q_{j}$ are the momentum, mass, and displacement from the equilibrium position of the $j$th particle, respectively, and $V$ is the nearest-neighboring inter-particle interaction potential. Without losing generality, we set $m_{1}=1-\Delta m/2$, and $m_{2}=1+\Delta m/2$ ($|\Delta m|<2$ to guarantee positive masses) so that the mean value of the masses was fixed to be a unit; that is, the mass’s density of the diatomic chain remained constant.

In our study, the Toda potential took the form of

[TABLE]

The 1D diatomic Toda lattice increasingly approximates an integrable lattice as $\Delta m$ approaches zero. Therefore, the diatomic Toda chain can be considered the perturbation of the Toda Hamiltonian. In this study, we only focused on the thermalization problem in the near-integrable region; thus, $\Delta m\ll 1$ was required.

For the sake of contrast, we also studied the diatomic FPUT-$\beta$ chain that possessed interaction potential

[TABLE]

where $\beta$ is a positive and free parameter. The diatomic FPUT-$\beta$ chain is always nonintegrable even if $\Delta m=0$. However, it will be the integrable diatomic harmonic lattice when $\beta=0$. Naturally, the diatomic FPUT-$\beta$ chain should be regarded as the integrable diatomic harmonic lattice perturbed by the quartic nonlinearity.

III Definition of perturbation strength

The Hamiltonian of a system can be written as

[TABLE]

where $H_{0}$ and $H^{\prime}$ denote the integrable part and the perturbation, respectively. Intuitively, the larger the perturbation is, the easier the system will be thermalized. To accurately characterize the ability of the system to be thermalized, it was necessary to properly define the perturbation strength, that is, to select the appropriate $H_{0}$. From the Hamiltonian canonical equation, it is easy to prove that the dynamical system described by Eq. (1) is mathematically and strictly equivalent to the homogeneous chain with unit mass described by the following Hamiltonian

[TABLE]

where $v_{j}$ is the velocity of $j$-th particle, and $1/m_{j}$ can be considered the renormalization coefficient of the force constant dependent on the lattice site. Based on the Eq. (5), the definitions of the perturbation strength for different cases are given below.

III.1 The diatomic Toda chains

For this case, the Hamiltonian of the Toda model was adopted as $H_{0}$, and

[TABLE]

Thus, we obtained the perturbation by comparing Eq. (5) and Eq. (6) as

[TABLE]

and the average strength of the perturbation for the fixed energy density

[TABLE]

where $\langle\cdot\rangle$ denotes the ensemble average of the thermodynamic equilibrium state. The perturbation strength being proportional to the mass difference illustrates the fact that the diatomic Toda chain is indeed the Toda Hamiltonian perturbed by unequal masses.

III.2 The diatomic FPUT-$\beta$ chains

For this case, the Hamiltonian of the diatomic harmonic model was adopted as $H_{0}$ and

[TABLE]

Comparing Eq. (5) and Eq. (9), we obtained the perturbation

[TABLE]

and the average strength of the perturbation for the fixed energy density

[TABLE]

The above expression clearly shows that the perturbation strength in the near-integrable region depended only on the nonlinear coefficients for this model.

IV Physical quantities and numerical method

The eigenvector, $\mathbf{u}_{k}$, of 1D diatomic chain for the fixed boundary conditions was derived in Chaturvedi and Baijal (1974). The element of $\mathbf{u}_{k}$ is

[TABLE]

where $\omega_{k}^{\pm}$ is the frequency of the $k$th eigenvector, as shown below:

[TABLE]

In this formula, ‘$-$’ and ‘$+$’ correspond to the acoustic branch and the optical branch, respectively, and $k=1,2,\cdots,L$. It should be noted that the eigenvectors given in formula (12) correspond to two vectors of the same wave number $k$ that are not orthogonally normalized. Therefore, by means of Gram-Schmidt’s orthogonalization technique Leon et al. (2013), the orthogonal and normalized eigenvector of the system can be obtained as follows:

[TABLE]

where $\|\cdot\|$ represents the length of a vector, and $\langle\mathbf{a},\mathbf{b}\rangle$ denotes the inner product of vectors $\mathbf{a}$ and $\mathbf{b}$. To produce a convenience of narration of the results, we sorted $N$ frequencies from small to large and removed the symbols ‘$\pm$’; i.e, $\omega_{k}=\omega_{k}^{-}$, and $\omega_{N-k+1}=\omega_{k}^{+}$ for $1\leq k\leq L$ (see Fig. 2). The superscript ‘$\pm$’ of the corresponding eigenvector was also removed. Namely, $\mathbf{U}_{k}=\mathbf{U}_{k}^{-}$, and $\mathbf{U}_{N-k+1}=\mathbf{U}_{k}^{+}$. Thus, the normal modes of 1D diatomic lattice are defined as

[TABLE]

The energy of the $k$th normal mode is

[TABLE]

and a phase $\varphi_{k}$ is defined via

[TABLE]

Following the definition of equipartition, one expects

[TABLE]

where $\varepsilon$ is the energy per particle and $\bar{E}_{k}(T)$ represents the time average of $E_{k}$ up to time $T$; i.e.,

[TABLE]

In the formula above, $\mu\in[0,1)$ controls the size of the time average window. In our numerical simulations, $\mu=2/3$ was fixed, which could not only speed up the calculations, but also had the advantage of a quicker loss of the memory of the very special initial state, as proposed in Benettin and Ponno (2011).

Based on the defined $\bar{E}_{k}(T)$, we needed to introduce a parameter to measure how close the system was to equipartition. A parameter frequently used for this purpose is the effective relative number of degrees of freedom Livi et al. (1985); Goedde et al. (1992). We employed the quantity $\xi(t)$, as described in Ref. Benettin and Ponno (2011), i.e.,

[TABLE]

and

[TABLE]

When equipartition was approached, $\xi$ will saturated at $1$.

To integrate the motion equations numerically, we used the eighth-order Yoshida method Yoshida (1990). The typical time step was $\Delta t=0.1$. To suppress fluctuations, the ensemble average was done over phases uniformly distributed in $[0,2\pi]$. In all our calculations below, the lowest $10\%$ of the frequency modes were excited. We have checked and verified that no qualitative difference will be resulted in when the percentage of the excited modes was changed.

V Numerical results

In Fig. 3(a), the results of $\langle E_{k}(t)/\varepsilon\rangle$ versus $k/N$ for the Toda chain are presented. It can be seen that only a small portion of the energy spread quickly from the initially excited low-frequency modes to the high-frequency modes. After this, the energy profile kept a stable form. This suggests that for the Toda model, the thermalized state can never be reached Fu et al. (2019). In Fig. 3(b), the results for the diatomic Toda chain with $\Delta m=0.01$ are plotted. It can be seen that the energy of the low-frequency modes (acoustic modes) was shared quickly while the energy of the high-frequency modes (optical modes) increased very slowly. Namely, the acoustic modes entered the equipartition state first, while the optical modes entered the equipartition state last. The system was eventually thermalized. Figure 3(c) shows the results of the diatomic Toda chain with $\Delta m=0.1$. It can be seen that the system reached the thermalized state faster. The system was fully thermalized at time $T\sim 10^{6}$, when $\langle E_{k}/\varepsilon\rangle=1$. As a comparison, Figs. 3(d)-(f) show the results of the diatomic FPUT-$\beta$ chain with $\Delta m=0$, $0.01$, and $0.1$, respectively. It can be seen that the results had a tiny difference, and all the systems were fully thermalized at nearly the same time, e.g., $T\sim 10^{8}$, when $\langle E_{k}/\varepsilon\rangle=1$. However, these results differ qualitatively from those of the diatomic Toda chains. In the FPUT-$\beta$ model, the high frequency modes (optical modes) entered the equipartition state first, while the low frequency modes (acoustic modes) entered the equipartition state last. This phenomenon is also observed in the 1D monatomic chains perturbed by nonlinearity Benettin and Ponno (2011); Fu et al. (2019). We judged that this qualitative difference in the thermalization process was caused by the difference between the mass perturbation and the nonlinearity perturbation.

In order to further clarify the difference of the rate of the energy sharing between the acoustic modes and the optical modes, we defined the quantity $\xi(t)_{A,O}$ for the acoustic modes and the optical modes as

[TABLE]

where $w_{k}(t)=\frac{\bar{E}_{k}(t)}{\sum_{j=1}^{L}\bar{E}_{j}(t)}$, and

[TABLE]

where $w_{k}(t)=\frac{\bar{E}_{k}(t)}{\sum_{j>L}^{N}\bar{E}_{j}(t)}$. Figure 4(a) shows the results for the diatomic Toda chains with fixed energy density $\varepsilon=0.01$, and different $\Delta m$; and Fig. 4(b) shows the results for the diatomic FPUT-$\beta$ chains. It can be seen that $\xi(t)_{A,O}$ quickly reached a constant value less than $1$ for the monatomic Toda chain; that is, it never reached an equipartition state due to its integrability. However, the acoustic and optical modes of the FPUT-$\beta$ model would eventually enter the state of energy sharing. Except for this difference, it can be clearly seen that there was a qualitative difference in the route to thermalization for these two types of diatomic chains. For example, in the diatomic Toda chain, the energy of the acoustic modes entered the equipartition state first, then the energy of the optical modes entered the equipartition state. In contrast, the thermalization process of the FPUT-$\beta$ chain had the opposite order. In the following sections, we will further examine the systematic study of the relationship between the thermalization time and the perturbation strength of the system.

To obtain the thermalization time, we studied the properties of $\langle\xi(t)\rangle$ defined by Eq. (21). Figure 5(a) shows the results for the diatomic Toda chain with the fixed energy density $\varepsilon=0.01$, and different $\Delta m$. It should be noted that on a sufficiently large time scale, all values of $\langle\xi(t)\rangle$ increased from [math] to $1$ with very similar sigmoidal profiles. This suggests that energy equipartition was finally achieved. Additionally, when $\Delta m$ decreased, the time required to reach the thermalized state increased. We adopt the definition of the equipartition time, $T_{eq}$, as the time when $\langle\xi(t)\rangle$ reached the threshold value $0.5$, as described in Benettin and Ponno (2011). By assuming that the threshold value $0.5$ was artificial, it did not influence the scaling law of $T_{eq}$. This can be seen from Fig. 5(b), where the sigmoidal profiles in Fig. 5(a) overlap with each other upon suitable shifts, which suggests that the concrete threshold value did not affect the scaling exponent of $T_{eq}$. With these preparations, we were ready to present the results of $T_{eq}$.

In Fig. 6(a), the numerical results of $T_{eq}$ as a function of $\Delta m$ are shown in log-log scale for the diatomic Toda chain with different energy densities $\varepsilon=0.001$ (red triangles), $\varepsilon=0.005$ (green triangles), and $\varepsilon=0.01$ (blue triangles); system size $N=1024$ (up triangles), and $N=2048$ (down triangles). It should be noted that all the points fall on the lines with a slope of $-2$, suggesting $T_{eq}\propto\Delta m^{-2}$; i.e., $T_{eq}$ was inversely proportional to the square of perturbation strength. In addition, the results of different sizes coincided completely, i.e., the size effect disappeared in our calculated range. For contrast, the numerical results of the diatomic FPUT-$\beta$ chain with the same conditions are presented in Fig. 6(b). It can be clearly seen that all the points fall on the lines with a slope of [math], suggesting $T_{eq}\propto\Delta m^{0}$; i.e., $T_{eq}$ did not depend on $\Delta m$ since the nonintegrability of the diatomic FPUT-$\beta$ chain hardly changed with the variation of $\Delta m$. This means that the results of the diatomic FPUT-$\beta$ chains are completely identical to those of the monatomic FPUT-$\beta$ chains; i.e., $T_{eq}\propto\beta^{-2}$ Fu et al. (2018); Pistone et al. (2018b). In addition, it should be noted that the results vary with the size of system under the same conditions, and the lower the energy density was, the more obvious the difference was; i.e., the stronger the size effect was. The phenomenon in which the symmetric model has a stronger size effect than the asymmetric model has been observed in the literatures Benettin and Ponno (2011); Fu et al. (2018).

VI Summary and discussions

In this study, we examined the thermalization behaviors of two types of 1D diatomic chains; i.e., the mass perturbed Toda chain (the diatomic Toda chain) and the quartic nonlinearity perturbed diatomic harmonic chain (the diatomic FPUT-$\beta$ chain). We showed that the acoustic mode in the diatomic Toda chain entered the state of energy sharing first, and the optical mode entered the equipartition state later. In contrast, for the diatomic FPUT-$\beta$ chain, the optical mode entered the state of energy sharing first, while the acoustic mode entered the equipartition state later. We considered that the qualitative difference in the thermalization process originated from the different methods of perturbation. Although these qualitative differences existed, the scaling law of thermalization time followed the same rule; the thermalization time was inversely proportional to the square of the perturbation strength. This law was also consistent with the thermalization law of 1D monatomic chains with nonlinearity perturbations Fu et al. (2018); Pistone et al. (2018b); Fu et al. (2019). All these results agree with the fact that the thermalization of a 1D chain is universal, and this universality is independent of the way in which nonintegrability is introduced. Specifically, the key to identifying the universal law is the selection of a suitable $H_{0}$ as the reference integrable system, so that the perturbation strength that reflects the ability of the system to be thermalized is defined accurately.

Acknowledgment

This work is supported by NSFC (Grant No. 11335006).