Thermal Broadening of Phonon Spectral Function in Classical Lattice Models: Projective Truncation Approximation

Hu-Wei Jia School of Physics, Renmin University of China, Beijing 100872, China Key Laboratory of Quantum State Construction and Manipulation (Ministry of Education), Renmin University of China, Beijing 100872, China Wen-Jun Liu School of Physics, Renmin University of China, Beijing 100872, China Key Laboratory of Quantum State Construction and Manipulation (Ministry of Education), Renmin University of China, Beijing 100872, China Yue-Hong Wu School of Physics, Renmin University of China, Beijing 100872, China Key Laboratory of Quantum State Construction and Manipulation (Ministry of Education), Renmin University of China, Beijing 100872, China Kou-Han Ma International Center for Quantum Materials and School of Physics, Peking University, Beijing 100871, China Lei Wang School of Physics, Renmin University of China, Beijing 100872, China Key Laboratory of Quantum State Construction and Manipulation (Ministry of Education), Renmin University of China, Beijing 100872, China Ning-Hua Tong [email protected] School of Physics, Renmin University of China, Beijing 100872, China Key Laboratory of Quantum State Construction and Manipulation (Ministry of Education), Renmin University of China, Beijing 100872, China

Abstract

Thermal broadening of the quasi-particle peak in the spectral function is an important physical feature in many statistical systems, but it is difficult to calculate. To tackle this problem, we propose the $H$ -expanded basis within the projective truncation approximation (PTA) of the Green’s function equation of motion. A zeros-removing technique is introduced to stabilize the iterative solution of the PTA equations. Benchmarking calculations on the classical one-variable anharmonic oscillator model and the one-dimensional $\phi^{4}$ lattice model show that the thermal broadened quasi-particle peak in the spectral function can be produced on a semi-quantitative level. Using this method, we discuss the low- and high- temperature power-law behaviors of the spectral width $\Gamma_{k}(T)$ of the one-dimensional $\phi^{4}$ model, finding it in contradiction with the assumption of effective phonon theory. A short-chain limit of this model is also discovered. Issues of extending the $H$ -expanded basis to quantum systems and of the applicability of the Debye formula for thermal conductivity are discussed.

I Introduction

The Green’s function (GF) is an important tool for studying the physical features of many-body systems in condensed-matter physics. The idea of projection [1, 2, 3, 4, 5] has been widely used in the calculation of the GF, leading to theories such as the two-pole approximation [6, 7], the composite operator method [8, 9], the self-consistent projection operator approach [10], the operator projection method [11], continued fraction [12, 13], the irreducible Green’s function method [14], mode-coupling theory [15, 16], and the dynamical projective operatorial approach [17], etc.. However, in practical applications, complicated analytical derivation and additional uncontrolled approximations are often required in these theories, which hampers systematic improvement of the results.

Recently, we developed a GF equation of motion (EOM) projective truncation approximation (PTA) method [18, 19, 20, 21]. The PTA overcomes the arbitrariness of traditional Tyablicov type truncation [22] by guaranteeing the analytical structure of the GF, diagonal-positivity of the spectral function, and the equilibrium state conservation of energy [18]. Compared to traditional projection-based methods, the PTA avoids the complicated analytical derivations to a large extent, and minimizes the need for additional uncontrolled approximations. Within the PTA, the study of a given system is reduced to selecting a basis of operators and solving the generalized eigenvalue problem on that basis. The power of modern computers is thus combined with the idea of projection. By a clever selection of the basis and systematically enlarging its size, exact results can be approached efficiently [19].

In the present work, we address the issue of describing the damping of quasi particles in the PTA. The damping (decay, relaxation, etc.) of quasi particles is ubiquitous in the experimental observations of spectra and transport. Modern numerical methods such as quantum Monte Carlo (QMC), density matrix renormalization group (DMRG), and tensor-network renormalization group (TNRG) methods have their own limitations in dealing with the damping. For examples, QMC requires an additional procedure to produce the real frequency spectral function, while DMRG and TNRG are more suitable for low dimensional systems at zero temperature. Theories based on Feynman diagrams usually apply only to weak correlation regime. In the traditional EOM formalism [4, 5] or the generalized Langevin equation approach [1, 2, 3], one needs to evaluate the imaginary part of the self-energy or the memory function to obtain a microscopic description of the damping, which is usually a difficult task. Within the PTA, although the quasi-particle energy can be obtained straightforwardly, the damping is also a challenge. This is because the PTA produces an approximate spectral function composed of $N$ discrete delta peaks, with $N$ being the dimension of the operator basis. Naively increasing $N$ to infinity is shown to be very inefficient for describing the broadening the spectral function and damping of quasi particles. See below.

Here, we propose the PTA on an $H$ -expanded basis to solve a part of this problem, i.e., the damping of quasi particles due to the thermal fluctuations at finite temperature. In the spectral function, it corresponds to the thermal broadening of the quasi-particle peak. The same idea can be extended to tackle the quantum damping and this issue will be addressed in the discussion part. In the framework of the PTA, our idea is first to produce the approximate eigen excitation operator $O$ , with $[O,H]\approx\epsilon O$ , from a low-order PTA. The corresponding spectral function has a delta peak at $\omega=\epsilon$ . Then we construct the new basis of operators by expanding $O$ with (exact or approximate) conserved quantities such as the Hamiltonian $H$ of the studied system or the mode Hamiltonian $H_{\nu k}$ of mode $(\nu k)$ . That is, we choose the new basis operators in the form $\{Of_{i}(H)\}$ , with $f_{i}$ ( $i=1,2,...,N$ ) being linearly independent functions. On the level of low-order PTA, these basis operators are degenerate excitation operators corresponding to the excitation energy $\epsilon$ . A new PTA calculation on the expanded basis will take into account the thermal fluctuation beyond the low-order PTA and hence split the degeneracy, producing refined spectral functions composed of $N$ closely distributed delta peaks around the low-order quasi-particle peak $\delta(\omega-\epsilon)$ . It provides a means to describe the thermal broadening of the quasi particle.

To illustration the implementation and effect of the above idea, we study the phonon damping problem in the classical lattice models which have only thermal damping of quasi particles. We consider two models. One is the one-variable anharmonic oscillator (AHO) model, whose static averages and spectral function are exactly solvable [23]. The other is the one-dimensional (1D) $\phi^{4}$ lattice model which is frequently used in the study of thermal transport in low-dimensional materials. For this model, molecular dynamical (MD) simulation results are available for benchmark.

At finite temperature, the phonon peak in the spectral function of 1D $\phi^{4}$ lattice will be broadened by thermal fluctuations due to the existence of nonlinear interaction. For general models, various methods have been developed to study the phonon decay and its consequence on the thermal transport, including the effective phonon theory [24], the anharmonic phonons [25] and the resonance phonon approaches [26, 27] based on MD simulation, non-equilibrium GF [28], series expansion [29], mode-coupling theory [15, 16], and the Boltzmann equation approach [30]. But a systematic study of the phonon decay for the 1D classical $\phi^{4}$ model, in particular, its low- and high-temperature asymptotic behaviors, is still lacking.

Our results show that PTA on the $H$ -expanded basis can produce semi-quantitatively accurate $\Gamma_{k}$ , the peak width of the phonon spectral function. It works for all temperatures, including the low temperature limit where MD simulation is too time consuming. Using this method, we find the low- and high-temperature asymptotic power laws in the temperature dependence of $\Gamma_{k}(T)$ of the 1D $\phi^{4}$ model.

Notably, in contrast to the assumption of the effective phonon theory $\Gamma_{k}\propto\epsilon\Omega_{k}$ [24], our results suggest a different relation $\Gamma_{k}(T)\approx c\left[\Omega_{k}(T)-\Omega_{k}(0)\right]$ , with $c$ being weakly dependent on $T$ and momentum $k$ . When inserted into the Debye formula [31, 32] for thermal conductivity $\kappa$ , it gives different temperature dependence of $\kappa(T)$ from that of MD. This raised question on the applicability of Debye formula for the 1D $\phi^{4}$ lattice model. We also find a short-chain limit of the $\phi^{4}$ model with distinct temperature dependence of $\Gamma_{k}(T)$ , which is amenable to experimental test in low dimensional materials. The method could be extended to quantum many-body systems where the damping of quasi particles may occur at zero temperature due to quantum fluctuations. We address this issue in the discussion section.

The rest of this paper is arranged as follows. For completeness, in Sec.II, we first review the formalism of GF EOM and PTA for classical systems. In Sec.III, we develop the formalism of the PTA on an $H$ -expanded basis for the AHO model. In Sec.IV, we apply this method to 1D $\phi^{4}$ lattice model, and we analyze the thermal broadening effect in phonon spectral function. Section V is devoted to a discussion and summary.

II Green’s Function EOM and PTA for classical systems

In this section, we review the formalism of the PTA for classical system, which was developed in Ref.[33]. This is to set the stage for developing the PTA on the $H$ -expanded basis.

For two dynamical variables at time $t$ and $t^{\prime}$ , $A[q(t),p(t)]$ and $B[q(t^{\prime}),p(t^{\prime})]$ ( $q$ and $p$ are abbreviation of canonical coordinates $q_{i}$ ’s and momenta $p_{i}$ ’s of the studied system), the retarded GF is defined as [34, 35, 36, 37]

G^{r}\left[A(t)|B(t^{\prime})\right]\equiv\theta(t-t^{\prime})\langle\{A(t),B(t^{\prime})\}\rangle.

(1)

Here $\theta(x)$ is the Heaviside step function. $\langle O\rangle$ is the average of variable $O$ in the equilibrium state of Hamiltonian $H(q,p)$ at temperature $T$ . $\{\,,\,\}$ is the Poisson bracket.

The Fourier transformation of the above GF is

G^{r}(A|B)_{\omega}=\int_{-\infty}^{\infty}G^{r}\left[A(t)|B(t^{\prime})\right]e^{i(t-t^{\prime})(\omega+i\eta)}d(t-t^{\prime}).

(2)

Here, $\eta$ is an infinitesimal positive number. By making derivative to time $t$ and $t^{\prime}$ , we obtain the EOM as [33]

	$\displaystyle\omega G(A\|B)_{\omega}$	$\displaystyle=$	$\displaystyle i\langle\{A,B\}\rangle+iG(\{A,H\}\|B)_{\omega}$		(3)
		$\displaystyle=$	$\displaystyle i\langle\{A,B\}\rangle-iG(A\|\{B,H\})_{\omega}.$		(3)

Here, the Zubarev GF $G(A|B)_{\omega}$ without superscript $r$ is related to the retarded GF $G^{r}(A|B)_{\omega}$ by $G^{r}(A|B)_{\omega}=G(A|B)_{\omega+i\eta}$ .

The average value $\langle BA\rangle$ can be obtained from GF through the fluctuation-dissipation theorem [34, 37]

\langle BA\rangle=\frac{1}{\beta}\int_{-\infty}^{\infty}\frac{\rho_{A,B}(\omega)}{\omega}d\omega+\langle B_{0}A_{0}\rangle.

(4)

Here, $A_{0}$ and $B_{0}$ are the zero-frequency components of $A(t)$ and $B(t)$ , respectively [33]. The spectral function $\rho_{A,B}(\omega)$ in the above equation is defined as

\rho_{A,B}(\omega)=\frac{i}{2\pi}\left[G(A|B)_{\omega+i\eta}-G(A|B)_{\omega-i\eta}\right].

(5)

To solve the EOM of GF by the PTA, we first choose a set of linearly-independent basis variables $A_{1},A_{2},...,A_{N}$ , which span an $N$ -dimensional subspace of variables. Due to the symmetry of the Hamiltonian of the type $H=\sum_{i}p_{i}^{2}/2\mu_{i}+V(q)$ , we can select the basis in one of the subspaces below,

	$\displaystyle\{O_{e}\}=\{f(q)\prod_{i}p_{i}^{m_{i}}\,\,\Big{\|}\,\,\sum_{i}m_{i}=2k,\,\,k\in\mathbb{Z}\},$
	$\displaystyle\{O_{o}\}=\{g(q)\prod_{i}p_{i}^{n_{i}}\,\,\Big{\|}\,\,\sum_{i}n_{i}=2k+1,\,\,k\in\mathbb{Z}\}.$

Here, $f(q)$ and $g(q)$ are arbitrary functions of $q$ . One can prove that the GF matrix $G(\vec{A}|\vec{A})_{\omega}$ , defined for the vector of variables $\vec{A}=\left(A_{1},A_{2},...,A_{N}\right)^{T}$ , fulfils the EOM

\displaystyle\omega^{2}G(\vec{A}\big{|}\vec{A}^{{\dagger}})_{\omega}=-\langle\{\{\vec{A},H\},\vec{A}^{{\dagger}}\}\rangle-G(\{\{\vec{A},H\},H\}\big{|}\vec{A}^{{\dagger}})_{\omega}.

To truncate the EOM in Eq.(II), we introduce the inner product of two classical variables $A$ and $B$ as

(A|B)\equiv\langle\{A^{\ast},\{B,H\}\}\rangle.

(8)

We then truncate the EOM as

\{\{\vec{A},H\},H\}\approx-{\bf M}^{T}\left(\vec{A}-\vec{A}_{0}\right),

(9)

Projecting this equation to each basis operator $A_{k}$ , one obtains the matrix equation ${\bf M}={\bf I}^{-1}{\bf L}$ . Here ${\bf I}$ is the inner product matrix, with element $I_{ij}\equiv(A_{i}|A_{j})$ . ${\bf L}$ is the Liouville matrix, with element $L_{ij}=-(A_{i}|\{\{A_{j},H\},H\})$ . Both are Hermitian and positive semi-definite matrices. ${\bf M}$ is thus guaranteed to have real positive eigenvalues. The GF matrix is formally expressed as

G(\vec{A}\big{|}\vec{A^{{\dagger}}})_{\omega}\approx\left(\omega^{2}-{\bf M}^{T}\right)^{-1}{\bf I}^{T}.

(10)

The spectral function obtained from Eqs.(5) and (10) reads

	$\displaystyle\rho_{A_{i},A_{j}^{\ast}}(\omega)$
	$\displaystyle\approx\sum_{k}\frac{\left({\bf IU}\right)^{\ast}_{ik}\left({\bf IU}\right)_{jk}}{2\sqrt{\epsilon_{k}}}\left[\delta\left(\omega-\sqrt{\epsilon_{k}}\right)-\delta\left(\omega+\sqrt{\epsilon_{k}}\right)\right].$

In this equation, ${\bf U}$ and ${\bf\Lambda}=\text{diag}(\epsilon_{1},\epsilon_{2},...,\epsilon_{N})$ are eigenvector and eigenvalue matrices of ${\bf M}$ , respectively. That is, ${\bf U}^{-1}{\bf M}{\bf U}={\bf\Lambda}$ . The eigenvalues $\epsilon_{k}$ ’s are real and $\epsilon_{k}\geqslant 0$ for all $k$ . They can be obtained by solving the generalized eigen-value problem,

{\bf LU}={\bf IU\Lambda}.

(12)

${\bf U}$ can be made to satisfy the generalized unitary condition ${\bf U}^{\dagger}{\bf I}{\bf U}={\bf 1}$ .

From the above expression for the spectral function, Eq.(4) gives the averages $C_{ij}\equiv\langle A_{i}^{\ast}A_{j}\rangle$ as

{\bf C}\approx\frac{1}{\beta}{\bf I}{\bf L}^{-1}{\bf I}+{\bf C}_{0}.

(13)

Here $(C_{0})_{ij}=\langle A_{i0}^{\ast}A_{j0}\rangle$ is the correlation of the static component of variables. This equation, once being used to calculate the averages appearing in ${\bf I}$ and ${\bf L}$ , closes the set of non-linear equations for the averages. The PTA solution to the averages can be obtained by self-consistently solving this set of equations. PTA solution to GF and spectral function are then obtained from Eqs.(10) and (II).

III One-Variable Anharmonic Oscillator Model

In this section, the idea of expanding the basis with $H$ to describe the thermal broadening of quasi-particle peak is applied to a simple classical statistical model, the one-variable AHO model. Its Hamiltonian reads

	$\displaystyle H=\frac{1}{2\mu}p^{2}+V(x),$
	$\displaystyle V(x)=\frac{1}{2}\mu\omega_{0}^{2}x^{2}+\alpha x^{4}.$		(14)

This model describes the motion of a single oscillator of mass $\mu$ in a quadratic potential with the basic frequency $\omega_{0}$ and an additional quartic potential of strength $\alpha$ . The spectral function $\rho_{x,x}(\omega)$ of variable $x$ in the equilibrium state was solved exactly in the canonical ensemble [23]. This is an ideal model for demonstrating the applicability of the PTA on the $H$ -expanded basis. At finite temperature, the spectral function $\rho_{x,x}(\omega)$ has a peak around the effective oscillating frequency $\Omega_{p}$ with width $\Gamma$ . The broadening of the spectral peak can be understood as follows. When the nonlinear potential is present, the oscillating frequency of an oscillator depends on its initial energy $E_{0}$ . In the equilibrium state of a finite temperature $T$ , $E_{0}$ has a thermal fluctuation of order $T$ in the canonical ensemble. This leads to an extended distribution of the oscillating frequencies, and subsequently to the broadening of the spectral peak.

III.1 basis $B_{1}\equiv\{x^{2m-1}p^{2n-2}\}$

Before we present the main subject of this paper, i.e., PTA calculation of the spectral function on the $H$ -expanded basis, let us first explore the convergence behavior of the spectral function from an ordinary basis, i.e., the basis composed of powers of basic variables $x$ and $p$ . The results will be contrasted with the PTA results from the $H$ -expanded basis.

If we start from the GF $G_{x,x}(\omega)$ and carry out the EOM an even number of times, we will get the higher order GFs $G_{x^{2m-1}p^{2n-2},x}(\omega)$ ( $m,n=1,2,...$ ). The specific form of the generated variable $x^{2m-1}p^{2n-2}$ is due to the parity symmetry of the potential $V(-x)=V(x)$ . We therefore define our basis variable as

A_{mn}=x^{2m-1}p^{2n-2},\,\,\,\,(m=1,2,...,M;\,n=1,2,...,N)

(15)

They form the basis $B_{1}=\{A_{mn}\}$ ( $m=1,2,...,M;\,n=1,2,...,N$ ), which span a $D=M\times N$ dimensional subspace. Due to the symmetry $V(-x)=V(x)$ , the static components of these variables are zero, $\left(A_{mn}\right)_{0}=0$ .

Following the standard formalism of PTA summarized in Sec.II, on this basis we have the inner product matrix

$\displaystyle I_{rs,mn}$	$\displaystyle\equiv$	$\displaystyle\left(A_{rs}\|A_{mn}\right)$	(16)
	$\displaystyle=$	$\displaystyle c_{1}\langle x^{2r+2m-4}\rangle\langle p^{2s+2n-4}\rangle$
		$\displaystyle+c_{2}\langle x^{2r+2m-3}V^{\prime}(x)\rangle\langle p^{2s+2n-6}\rangle$
		$\displaystyle+c_{3}\langle x^{2r+2m-2}V^{\prime\prime}(x)\rangle\langle p^{2s+2n-6}\rangle,$

with $c_{1}=[(2m-1)/\mu]\left[(2r-1)(2n-1)-(2s-2)(2m-2)\right]$ , $c_{2}=-(2n-2)\left[(2r-1)(2n-3)-(2s-2)(2m-1)\right]$ , and $c_{3}=(2n-2)(2s-2)$ . The Liouville matrix is obtained as

$\displaystyle L_{rs,mn}$	$\displaystyle=$	$\displaystyle-\frac{1}{\mu^{2}}(2m-1)(2m-2)\,I_{rs,m-1\,n+1}$	(17)
		$\displaystyle+\omega_{0}^{2}\left[(2m-1)(4n-3)+2n-2\right]\,I_{rs,mn}$
		$\displaystyle+\frac{4\alpha}{\mu}\left[(2m-1)(4n-3)+3(2n-2)\right]\,I_{rs,m+1\,n}$
		$\displaystyle-\mu^{2}\omega_{0}^{4}(2n-2)(2n-3)\,I_{rs,m+1\,n-1}$
		$\displaystyle-16\alpha^{2}(2n-2)(2n-3)\,I_{rs,m+3\,n-1}$
		$\displaystyle-8\alpha\mu\omega_{0}^{2}(2n-2)(2n-3)\,I_{rs,m+2\,n-1}.$

The calculation of matrices ${\bf I}$ and ${\bf L}$ requires evaluating the averages of the type $\langle x^{2k}\rangle$ and $\langle p^{2k}\rangle$ . The latter can be obtained exactly from the Gaussian integral, which gives $\langle p^{2k}\rangle=(2k-1)!!(\mu T)^{k}$ ( $k=1,2,...$ ). $\langle x^{2k}\rangle$ needs to be computed self-consistently from GFs. However, since our purpose is to check the convergence behavior of the spectral function on this basis, we use the numerically exact value of $\langle x^{2k}\rangle$ .

From ${\bf I}$ and ${\bf L}$ , we can calculate the spectral function $\rho_{x,x}(\omega)$ from the following formalism,

\rho_{x,x}(\omega)=\sum_{k}\frac{(\bf{IU})_{1k}^{2}}{2\sqrt{\epsilon_{k}}}\left[\delta(\omega-\sqrt{\epsilon_{k}})-\delta(\omega+\sqrt{\epsilon_{k}})\right].

(18)

In the above equation, the matrix $\bf{U}$ and the excitation energy $\epsilon_{k}$ ’s are determined by the generalized eigen-value problem

\bf{LU}=\bf{IU\Lambda},

(19)

where $\bf{U}$ is the generalized eigen vector matrix, and ${\bf\Lambda}=\text{diag}(\epsilon_{1},\epsilon_{2},...,\epsilon_{D})$ is the generalized eigen value matrix.

III.2 $H$ -expanded basis $B_{2}\equiv\{xe^{\lambda_{i}H}\}$ and $B_{3}\equiv\{x,x^{3}\}\otimes\{e^{\lambda_{i}H}\}$

In this subsection, we implement PTA on two successively larger $H$ -expanded bases, $B_{2}\equiv\{xe^{\lambda_{i}H}\}$ and $B_{3}\equiv\{x,x^{3}\}\otimes\{e^{\lambda_{i}H}\}$ ( $i=1,2,...,N$ ). Here $\{\lambda_{i}\}$ are $N$ different real numbers. We choose the exponential function $e^{\lambda_{i}H}$ to expand the low order basis variables $\{x\}$ and $\{x,x^{3}\}$ , because with this functional form, the calculation of ${\bf I}$ and ${\bf L}$ can be simplified significantly. For details, see Appendix A. The PTA results from bases $B_{2}$ and $B_{3}$ will be compared with the exact ones.

III.2.1 basis $B_{2}\equiv\{xe^{\lambda_{i}H}\}$

We first derive PTA equations based on the $H$ -expanded basis $B_{2}$ composed of variables

A_{i}\equiv xf_{i},\,\,\,\,\,(i=1,2,...,N).

(20)

Here, $f_{i}=e^{\lambda_{i}H}$ . To put the physically interesting variable $x$ into the basis, we assign $\lambda_{1}\equiv 0$ , meaning $A_{1}=x$ . $\{A_{i}\}$ form a $N$ -dimensional basis, with zero static components $A_{i0}=0$ . The parametrization of $\{\lambda_{i}\}$ will be specified in Sec.III.D.

The inner product matrix $\bf{I}$ and Liouville matrix $\bf{L}$ are derived as

	$\displaystyle I_{ij}=\frac{\beta}{\mu(\beta-\lambda_{i}-\lambda_{j})}\langle f_{i}f_{j}\rangle,$
	$\displaystyle L_{ij}=\frac{\beta}{\mu^{2}(\beta-\lambda_{i}-\lambda_{j})}\left[\mu\omega_{0}^{2}\langle f_{i}f_{j}\rangle+12\alpha\langle x^{2}f_{i}f_{j}\rangle\right].$

Both ${\bf I}$ and ${\bf L}$ contain the correlation among the expanding factors $\{e^{\lambda_{i}H}\}$ that takes into account the energy fluctuation effect. Details of the derivation are given in Appendix A.

To evaluate the averages appearing in $I_{ij}$ and $L_{ij}$ , we first note that $\langle x^{2}f_{i}f_{j}\rangle$ in $L_{ij}$ can be expressed in terms of $C_{ij}\equiv\langle A_{i}A_{j}\rangle$ and thus can be calculated self-consistently from the fluctuation-dissipation theorem Eq.(13) as

	$\displaystyle C_{ij}$	$\displaystyle=$	$\displaystyle\langle x^{2}f_{i}f_{j}\rangle$		(22)
		$\displaystyle=$	$\displaystyle\frac{1}{\beta}\left({\bf I}{\bf L}^{-1}{\bf I}\right)_{ij}.$		(22)

Here, ${\bf C}_{0}=0$ has been used. In particular, we have $\langle x^{2}\rangle=C_{11}$ . The other unknown average $\langle f_{i}f_{j}\rangle$ cannot be written in the form $\langle A_{i}A_{j}\rangle$ . For this specific model, it can be evaluated by numerical integration, or by a quadratic variational approximation [38]. We have compared the two methods and find little difference in the final results. Here we use the latter method,

	$\displaystyle\langle f_{i}f_{j}\rangle$	$\displaystyle\approx$	$\displaystyle\langle f_{i}f_{j}\rangle_{vari}$		(23)
		$\displaystyle=$	$\displaystyle\frac{\int\int_{-\infty}^{\infty}dxdp\,e^{-(\beta-\lambda_{i}-\lambda_{j})H_{vari}}}{\int\int_{-\infty}^{\infty}dxdp\,e^{-\beta H_{vari}}}.$		(23)

The variational Hamiltonian $H_{vari}$ is obtained from a quadratic variational approximation as

H_{vari}=\frac{1}{2\mu}p^{2}+\theta x^{2},

(24)

with $\theta=\frac{1}{4}\mu\omega_{0}^{2}+\sqrt{\left(\frac{1}{4}\mu\omega_{0}^{2}\right)^{2}+3\alpha T}$ . Using this variational Hamiltonian, it is easy to obtain the expression for $\langle f_{i}f_{j}\rangle$ as

\langle f_{i}f_{j}\rangle\approx\frac{\beta}{\beta-\lambda_{i}-\lambda_{j}},

(25)

being independent of the parameters in $H$ .

The spectral function $\rho_{x,x}(\omega)$ is then calculated from Eqs.(18) and (19).

III.2.2 basis $B_{3}\equiv\{x,x^{3}\}\otimes\{e^{\lambda_{i}H}\}$

Here, we derive the ${\bf I}$ and ${\bf L}$ matrices on basis $B_{3}$ . We denote the basis operator as

A_{\mu i}=a_{\mu}f_{i},\,\,\,\,(\mu=1,2;\,\,i=1,2,...,N).

(26)

where $a_{1}=x$ and $a_{2}=x^{3}$ , and $f_{i}=e^{\lambda_{i}H}$ ( $i=1,2,...,N$ ). These variables also have zero static components $\left(A_{\mu i}\right)_{0}=0$ . The basis $B_{3}\equiv\{a_{\mu}f_{i}\}$ ( $\mu=1,2;\,\,i=1,2,...,N$ ) is a $2N$ dimensional subspace, being larger than $B_{2}$ .

We denote the matrix elements of ${\bf I}$ and ${\bf L}$ as $I_{\mu\nu}^{ij}=\left(A_{\mu i}|A_{\nu j}\right)$ , and $L_{\mu\nu}^{ij}=-\left(A_{\mu i}|\{\{A_{\nu j},H\},H\}\right)$ , respectively. Using the same method for deriving ${\bf I}$ and ${\bf L}$ on the $B_{2}$ basis and some additional simplifications, we obtain

\displaystyle{\bf I}^{ij}=\frac{\beta}{\mu\left(\beta-\lambda_{i}-\lambda_{j}\right)}\left(\begin{array}[]{cc}\langle f_{i}f_{j}\rangle&3\langle x^{2}f_{i}f_{j}\rangle\\ 3\langle x^{2}f_{i}f_{j}\rangle&9\langle x^{4}f_{i}f_{j}\rangle\\ \end{array}\right).

(29)

The matrix elements of ${\bf L}$ read

$\displaystyle L_{11}^{ij}$	$\displaystyle=$	$\displaystyle\frac{\beta}{\mu^{2}\left(\beta-\lambda_{i}-\lambda_{j}\right)}\left[\mu\omega_{0}^{2}\langle f_{i}f_{j}\rangle+12\alpha\langle x^{2}f_{i}f_{j}\rangle\right],$
$\displaystyle L_{12}^{ij}$	$\displaystyle=$	$\displaystyle L_{21}^{ij}$
	$\displaystyle=$	$\displaystyle\frac{3\beta}{\mu^{2}\left(\beta-\lambda_{i}-\lambda_{j}\right)}\left[\mu\omega_{0}^{2}\langle x^{2}f_{i}f_{j}\rangle+12\alpha\langle x^{4}f_{i}f_{j}\rangle\right],$
$\displaystyle L_{22}^{ij}$	$\displaystyle=$	$\displaystyle-\frac{54\beta}{\mu^{2}\left(\beta-\lambda_{i}-\lambda_{j}\right)^{2}}\langle x^{2}f_{i}f_{j}\rangle$
		$\displaystyle+\frac{\beta}{\beta-\lambda_{i}-\lambda_{j}}\left[\frac{63\omega_{0}^{2}}{\mu}\langle x^{4}f_{i}f_{j}\rangle+\frac{324\alpha}{\mu^{2}}\langle x^{6}f_{i}f_{j}\rangle\right].$

Details of the derivation are given in Appendix A.

$C_{\mu\nu}^{ij}=\langle a_{\mu}a_{\nu}f_{i}f_{j}\rangle$ is still given by Eq.(13), with ${\bf C}_{0}=0$ . Some unknown averages in ${\bf I}$ and ${\bf L}$ can be calculated self-consistently through

	$\displaystyle\langle x^{2}f_{i}f_{j}\rangle=C_{11}^{ij},$
	$\displaystyle\langle x^{4}f_{i}f_{j}\rangle=C_{12}^{ij},$
	$\displaystyle\langle x^{6}f_{i}f_{j}\rangle=C_{22}^{ij}.$		(32)

The other quantities, including $\langle f_{i}f_{j}\rangle$ , GF, and the spectral function $\rho_{x,x}(\omega)$ , can be calculated similarly as for basis $B_{2}$ .

III.3 removing the redundant basis variables

In our numerical solution of the above PTA equations, we supplemented a technique that we called the zeros-removing method to stabilize the iteration process. It will be useful in general PTA calculation. Therefore, in this subsection, we introduce this method in the general framework of PTA.

In principle, the general PTA formalism works for arbitrary basis dimension $N$ . In practice, however, if the PTA calculation is implemented naively, it is limited to some dimension $N$ by numerical instability. This is because the orthogonality of the basis operators are difficult to control and they easily become linearly dependent when $N$ is large. As a result, the inner product matrix $\bf{I}$ , which should be positive-definite in theory, often becomes nearly singular for large $N$ . This makes the iterative process unstable. This problem, arising from the non-orthogonality of the basis variables, is quite common in PTA calculations. For example, in our study of Anderson impurity model with the PTA [18], the singularity of $\bf{I}$ also appears.

To partly overcome this difficulty, here we propose a method commonly used in the ab initio calculation with unorthogonal basis [39]. The idea of the method is to first figure out the redundant basis variables by diagonalizing ${\bf I}$ and find certain smallest eigen values and the associated eigen vectors. Then one removes these vectors from the basis and redoes the PTA calculation on the remaining basis. Our calculation shows that this method can significantly stabilize the iteration and makes it possible to do PTA calculation on a much larger basis. Below, we present the formalism of this method.

Suppose we have a set of basis variables $\{A_{i}\}$ ( $i=1,2,...,N$ ). When $\bf{I}$ and $\bf{L}$ matrices are obtained, we first diagonalize $\bf{I}$ using the unitary transformation,

({\bf U}_{I})^{-1}{\bf I}{\bf U}_{I}={\bf\Lambda}_{I}.

(33)

Here ${\bf U}_{I}$ is an unitary matrix and ${\bf\Lambda}_{I}$ is diagonal. Among the $N$ eigen values of ${\bf I}$ , suppose the smallest $N_{1}$ of them are close to zero, with their absolute value below a criterion, say, $10^{-12}$ . The other $N_{2}=N-N_{1}$ eigenvalues $\{i_{1},i_{2},...,i_{N_{2}}\}$ are larger than this criterion. Then we can approximately write the matrix ${\bf\Lambda}_{I}$ in the following block form,

{\bf\Lambda}_{I}\approx\left(\begin{array}[]{cc}{\bf 0}&{\bf 0}\\ {\bf 0}&{\bf\Lambda}_{I2}\\ \end{array}\right),

(34)

with ${\bf\Lambda}_{I2}={\text{diag}}(i_{1},i_{2},...,i_{N_{2}})$ and $i_{k}>0$ ( $k=1,2,...,N_{2}$ ). We can also write ${\bf U}_{I}$ as

{\bf U}_{I}=\left(\begin{array}[]{cc}{\bf U}_{I1}&{\bf U}_{I2}\end{array}\right),

(35)

where $\left({\bf U}_{I1}\right)_{N\times N_{1}}$ contains the $N_{1}$ column vectors of the zero-length variables, while $\left({\bf U}_{I2}\right)_{N\times N_{2}}$ contains the $N_{2}$ column vectors of finite-length variables. From the orthonormal and complete conditions ${\bf U}_{I}^{\dagger}{\bf U}_{I}={\bf U}_{I}{\bf U}_{I}^{\dagger}={\bf 1}$ , we have the following properties for ${\bf U}_{I1}$ and ${\bf U}_{I2}$ ,

	$\displaystyle{\bf U}_{I1}{\bf U}_{I1}^{\dagger}+{\bf U}_{I2}{\bf U}_{I2}^{\dagger}={\bf 1}_{N\times N},$
	$\displaystyle{\bf U}_{I1}^{\dagger}{\bf U}_{I1}={\bf 1}_{N_{1}\times N_{1}},\,\,\,\,\,\,\,\,{\bf U}_{I2}^{\dagger}{\bf U}_{I2}={\bf 1}_{N_{2}\times N_{2}},$
	$\displaystyle{\bf U}_{I1}^{\dagger}{\bf U}_{I2}={\bf 0}_{N_{1}\times N_{2}},\,\,\,\,\,\,\,\,{\bf U}_{I2}^{\dagger}{\bf U}_{I1}={\bf 0}_{N_{2}\times N_{1}}.$		(36)

We normalize the finite-length variables by defining matrix $\tilde{\bf U}_{2}={\bf U}_{I2}\frac{1}{\sqrt{{\bf\Lambda}_{I2}}}$ . Now we can write down the orthogonalized variables. The zero-length variables are denoted as $O_{1k}$ and the normalized finite-length variables as $O_{2k}$ . They read

	$\displaystyle O_{1k}=\sum_{i=1}^{N}A_{i}\left(U_{I1}\right)_{ik},$
	$\displaystyle O_{2k}=\sum_{i=1}^{N}A_{i}\left(\tilde{U}_{2}\right)_{ik}.$		(37)

Since $O_{1k}$ has zero length, one can prove that the GF involving $O_{1k}$ are zero, i.e., $G(O_{1k}|O_{1p})_{\omega}=G(O_{1k}|O_{2p})_{\omega}=0$ . Its static correlation function is also zero, i.e., $\langle O_{1k}X\rangle=0$ for arbitrary variable $X$ . We can therefore focus only on the remaining $N_{2}$ -dimensional subspace. The original basis variable $A_{i}$ is now expressed in terms of $O_{1k}$ and $O_{2k}$ as

A_{i}=\sum_{k=1}^{N_{1}}O_{1k}\left({\bf U}_{I1}^{\dagger}\right)_{ki}+\sum_{k=1}^{N_{2}}O_{2k}\left({\bf\Lambda}_{I2}\tilde{\bf U}_{2}^{\dagger}\right)_{ki}.

(38)

The Liouville matrix on the normalized basis $O_{2k}$ reads

\tilde{\bf L}_{22}=\tilde{\bf U}_{2}^{\dagger}{\bf L}\tilde{\bf U}_{2}.

(39)

Here $\tilde{\bf L}_{22}$ is a $N_{2}\times N_{2}$ Hermitian matrix. We diagonalize it by a unitary transformation

[{\bf V}_{22}]^{-1}\tilde{\bf L}_{22}{\bf V}_{22}={\bf\Lambda}_{2},

(40)

with ${\bf\Lambda}_{2}=\text{diag}(\epsilon_{1},\epsilon_{2},...,\epsilon_{N_{2}})$ being the square of the $N_{2}$ eigen excitation energies.

In the $N_{2}$ -dimensional space of the normalized basis variables, we employ the standard PTA to write down the GF matrix as

	$\displaystyle\tilde{\bf G}_{22}(\omega)$	$\displaystyle=$	$\displaystyle\left(\omega^{2}{\bf 1}-\tilde{\bf L}_{22}^{T}\right)^{-1}$		(41)
		$\displaystyle=$	$\displaystyle\left({\bf V}_{22}^{-1}\right)^{T}\left(\omega^{2}{\bf 1}-{\bf\Lambda}_{2}\right)^{-1}{\bf V}_{22}^{T}.$		(41)

In the above equation, $(\tilde{G}_{22})_{kp}(\omega)=G(O_{2k}|O_{2p}^{\ast})_{\omega}$ . Using Eq.(38), we obtain the GF of original basis variables $A_{i}$ as

G(A_{i}|A_{j}^{\ast})_{\omega}=\sum_{k,p=1}^{N_{2}}\sqrt{i_{k}i_{p}}\left(U_{I2}\right)^{\ast}_{ik}\left(U_{I2}\right)_{jp}(\tilde{G}_{22})_{kp}(\omega).

(42)

The spectral function is obtained similarly as

	$\displaystyle\tilde{\rho}_{O_{2k},O_{2p}^{\ast}}(\omega)$
	$\displaystyle=\sum_{m=1}^{N_{2}}\frac{\left(V_{22}\right)^{\ast}_{km}\left(V_{22}\right)_{pm}}{2\sqrt{\epsilon_{m}}}\left[\delta(\omega-\sqrt{\epsilon_{m}})-\delta(\omega+\sqrt{\epsilon_{m}})\right],$

and

\rho_{A_{i},A_{j}^{\ast}}(\omega)=\sum_{k,p=1}^{N_{2}}\sqrt{i_{k}i_{p}}\left(U_{I2}\right)^{\ast}_{ik}\left(U_{I2}\right)_{jp}\tilde{\rho}_{O_{2k},O_{2p}^{\ast}}(\omega).

(44)

From the spectral theorem, the averages $(\tilde{C}_{22})_{kp}=\langle O_{2k}^{\ast}O_{2p}\rangle$ are calculated as

\tilde{\bf C}_{22}=\frac{1}{\beta}\left(\tilde{\bf L}_{22}\right)^{-1}.

(45)

The correlation function between the original basis variables reads

\langle A_{j}^{\ast}A_{i}\rangle=\sum_{k,p=1}^{N_{2}}\sqrt{i_{k}i_{p}}\left(U_{I2}\right)^{\ast}_{ik}\left(U_{I2}\right)_{jp}\left[\tilde{C}_{22}\right]_{pk}.

(46)

The closed equations for $\langle A_{j}^{\ast}A_{i}\rangle$ are obtained by expressing ${\bf I}$ and ${\bf L}$ in terms of $\langle A_{j}^{\ast}A_{i}\rangle$ . They can be solved iteratively. Once $\langle A_{j}^{\ast}A_{i}\rangle$ are obtained, the GF and spectral functions can be calculated from Eqs.(42) and (44).

III.4 results for one-variable AHO model

In this subsection, we present the results for the one-variable AHO model, obtained from the basis $B_{1}\equiv\{x^{2m-1}p^{2n-2}\}$ ( $1\leq m\leq M$ , $1\leq n\leq N$ ), $B_{2}\equiv\{xe^{\lambda_{i}H}\}$ ( $1\leq i\leq N$ ), and $B_{3}\equiv\{x,x^{3}\}\otimes\{e^{\lambda_{i}H}\}$ ( $1\leq i\leq N$ ). We compare these results with those from exact numerical integration (for averages) and exact formula [23] (for the spectral function). We want to show that, although the $B_{1}$ basis becomes complete in the limit $M=N=\infty$ , the shape of the primary peak in the spectral function converges very slowly. The $B_{2}$ basis quite faithfully produces the shape of the main peak with moderate dimension $N$ . The $B_{3}$ basis produces both the main peak and the secondary peak. These results demonstrate that the $H$ -expanded basis can efficiently describe the thermal broadening effect in the spectral function. For the data shown below, we typically use the smallness criterion $10^{-14}\sim 10^{-9}$ for truncating the spectrum of ${\bf I}$ and the convergence criterion $10^{-9}\sim 10^{-7}$ for iterative solution of $\langle A_{i}^{\ast}A_{j}\rangle$ . All our results in this paper are in the unit $k_{B}=1$ .

The $H$ -expanded bases $B_{2}$ and $B_{3}$ are characterized by $N$ real parameters $\{\lambda_{i}\}$ ( $i=1,2,...,N$ ). We find that the value of these parameters will influence the convergence speed of PTA calculation and the quality of the spectral function. Since in $\langle f_{i}f_{j}\rangle\sim\text{Tr}[e^{-(\beta-\lambda_{i}-\lambda_{j})}]$ , $\lambda_{i}+\lambda_{j}<\beta$ is required, we let $\lambda_{i}\in(-\infty,\beta/2)$ . The lower and higher $\lambda_{i}$ ’s are responsible for describing low and high energy thermal fluctuations, respectively. We use the following scheme for assigning them,

\lambda_{i}=\left\{\begin{array}[]{lll}\frac{1}{2}\left(\frac{1}{T}-\frac{1}{T_{i}}\right)&\,\,\,(2\leq i\leq N),\\ \\ 0&\,\,\,(i=1).\end{array}\right.

(47)

Here, we set $\lambda_{1}=0$ to include $x$ in the basis. We let $\ln{T_{i}}$ distribute equal-distantly in the interval $[\ln{T}-\Delta,\ln{T}+\Delta]$ . That is, for $N>1$ , we use

\ln{T_{i}}=\ln{T}-\Delta+(i-1)\frac{2\Delta}{N-1},\,\,\,\,\,\,\,\,(2\leq i\leq N).

(48)

Typically $\Delta=4.0\sim 9.0$ is used in our calculation.

Our test shows that the above scheme works quite well for the AHO model. The peaks in the spectral function $\rho_{x,x}(\omega)$ shift slightly with $\Delta$ in the range $\Delta=1.0\sim 10.0$ . The key features extracted from $\rho_{x,x}(\omega)$ , such as the peak position $\Omega_{p}$ , peak width $\Gamma$ , and the moments of the spectral function, change very little. As $\Delta$ decreases, $\Gamma$ stays qualitatively the same until $\Delta\sim 10^{-6}$ , at which point $\Gamma$ jumps to zero, recovering the expected result for $\Delta=0$ . For tiny but finite $\Delta$ , $\rho_{x,x}(\omega)$ has two delta peaks, since the basis is effectively reduced to a two-dimensional basis $\{x,xH\}$ in the small $\Delta$ limit. The same holds for the 1D $\phi^{4}$ lattice model. When we need to average the spectral function $\rho_{x,x}(\omega)$ over $N_{r}$ random assignments of $\{\lambda_{i}\}$ , we let $\ln{T_{i}}$ take random numbers from the uniform distribution within $[\ln{T}-\Delta,\ln{T}+\Delta]$ .

In this work, when presenting the spectral function, we always broaden the delta peaks in the spectral function by a Gaussian function with standard variance $\sigma$ . Depending on the occasion, the value of $\sigma$ is either fixed to be $0.2\sim 0.05$ , or it is taken to be proportional to the width of the main peak $\Gamma$ , $\sigma=r\Gamma$ with $r=0.2\sim 0.3$ .

Considering that the spectral function in the present study always has a single primary peak in the $\omega>0$ regime, and it is almost symmetric in shape, we extract the peak position $\Omega_{p}$ of the primary peak through

\Omega_{p}=m_{1},

(49)

and the width $\Gamma$ of the primary peak through

\Gamma=\sqrt{m_{2}-m_{1}^{2}}.

(50)

Here, $m_{k}$ is the $k$ -th order moment of the spectral function on the positive frequency axis. This method was used in, e.g., defining the spin wave line width [40]. $m_{k}$ is evaluated from the weight $w_{i}$ and pole $p_{i}$ of the raw data of spectral function $\rho_{x,x}(\omega)=\sum_{i}w_{i}\delta(\omega-p_{i})$ as $m_{k}=(\sum_{i,p_{i}>0}w_{i}p_{i}^{k})/(\sum_{i,p_{i}>0}w_{i})$ . In this way, both $\Omega_{p}$ and $\Gamma$ are independent of the broadening procedure used for plotting the curve of $\rho_{x,x}(\omega)$ . Note that in the literature, there are various ways for extracting the width of phonon peak from the MD simulation data [25, 26, 41, 42].

Refer to caption — Figure 1: (color online) $\langle x^{2}\rangle$ as functions of $T$ for AHO model. Solid lines are obtained from $B_{1}$ -PTA with $M=1$ , $2$ , $3$ , and $6$ from bottom to top ( $N=M$ is always used). The dashed line is for the exact solution. Inset: The errors with respect to the exact results as functions of $M$ ( $N=M$ ) in $\langle x^{2}\rangle$ (red squares), $\langle x^{4}\rangle$ (green dots), and $\langle x^{6}\rangle$ (blue triangles) at $T=0.3$ . The dashed line is $M^{-2}$ for guiding eyes. Model parameters are $\mu=1.0$ , $\omega_{0}=0.3$ , and $\alpha=0.25$ .

Below, we first examine the convergence behavior of the basis $B_{1}$ . In the limit $M=N=\infty$ , $B_{1}$ becomes complete and the results from PTA are expected to become exact. In Fig.1, we show the average value $\langle x^{2}\rangle$ as a functions of $T$ for various $M$ (fixing $N=M$ ). It clearly shows that as $M=N$ increases, the curve uniformly approaches the exact dashed line. The inset of Fig.1 shows that the errors of $\langle x^{2}\rangle$ , $\langle x^{4}\rangle$ , and $\langle x^{6}\rangle$ at $T=0.3$ decrease with increasing $M$ as $M^{-2}$ (fixing $N=M$ ). It is also observed that at about $M=N=10$ , i.e., at the dimension of basis $D=M\times N=100$ , the error from PTA abruptly increases. We have observed that the largest (smallest) eigenvalue of ${\bf I}$ increases (decreases) exponentially with increasing $M=N$ . At about $M=N=10$ , the condition number of ${\bf I}$ deteriorates to such an extent that the iterative solution of PTA equations can no longer be carried out accurately. In our previous application of PTA [18], we also observed similar phenomena, hinting that this problem may be universal for PTA. Note that this problem can not be solved by the use of zeros-removing method discussed in Sec.III.C.

In Fig.2, the spectral function $\rho_{x,x}(\omega)$ of AHO model at $T=0.3$ is shown. In each panel, the red curve (from PTA on basis $B_{1}$ ) is compared to the black dashed line (from exact solution [23]). The exact $\rho_{x,x}(\omega)$ has a series of broad peaks at successively higher frequencies, whose weights decay exponentially with frequency. These multiple peaks are the overtone phenomenon due to the non-linear confinement, broadened by thermal fluctuations. At very low temperature, the peak positions tend to frequencies $\omega_{0}$ , $3\omega_{0}$ , $5\omega_{0}$ , etc.. The exact spectral function vanishes abruptly for $\omega<\omega_{0}$ , in agreement with the fact that $\omega_{0}$ is the lowest oscillation frequency for all initial conditions.

It is seen that as $M=N$ increases up to $M=N=8$ , the PTA-produced spectral function contains more and more peaks (we broaden the delta peaks a little for presentation purpose). The peak positions and weights are distributed in such a way that the envelope curve uniformly approaches the exact one, including the primary peak and the higher order ones. Fig.1 and Fig.2 show that both the thermodynamical and the dynamical quantities indeed converge towards the exact value as the dimension of basis increases.

However, if we look at the primary peak alone in Fig.3 on the linear scale, we will find that with increasing $M=N$ , the primary peak of $\rho_{x,x}(\omega)$ from PTA receives little improvement. Up to $M=N=8$ , the primary peak has got only $4$ discernible peaks on the linear scale. Clearly, on the naive basis $B_{1}$ , the poles of GF distribute uniformly in the whole frequency axis. Among the $M\times N$ delta peaks, most of them have high frequencies and tiny weights, such that they contribute little to the broad primary peak. This leads to a very slow convergence in the shape of the primary peak with increasing basis dimension, although the extracted peak position $\Omega_{p}$ and width $\Gamma$ converge fast with increasing $M=N$ , as shown in the inset of Fig.3.

Fig.4 shows the spectral function $\rho_{x,x}(\omega)$ obtained from PTA on the $H$ -expanded basis $B_{2}$ (red solid lines), together with the exact curve (black dashed lines). They are shown in log-linear axes. In contrast to the case of basis $B_{1}$ (shown in Fig.3), here as the dimension of basis increases, in the frequency regime of the primary peak, one gets significantly more delta peaks. The envelope tends to the exact broad primary peak much faster than that of $B_{1}$ basis. A basis with dimension $N=20$ already produces about $10$ visible peaks densely located within the frequency regime of primary peak, leading to qualitatively accurate description of the primary peak, if each delta peak of PTA is properly broadened with a Gaussian function. However, $B_{2}$ basis also meets its limitation at about $N=30$ , beyond which the spectral function can hardly receive any improvement. This is due to the deterioration of the condition number of ${\bf I}$ matrix as the basis dimension increases.

In Fig.5, we compare the spectral function $\rho_{x,x}(\omega)$ averaged over $500$ random $\lambda_{i}$ realizations (solid lines) with the exact curve (dashed lines) for different temperatures. Only the primary peak is shown in this plot with a linear-linear axis. The agreement is quite good in a wide temperature range. The exact $\rho_{x,x}(\omega)$ curve has a natural frequency cut off at $\omega=\omega_{0}$ , below which there is no spectral weight. This feature is nicely captured by PTA on basis $B_{2}$ with averaging on $500$ random $\{\lambda_{i}\}$ ’s (not shown in Fig.5). The inset shows that the relative peak position $\Omega_{p}-\omega_{0}$ and the peak width $\Gamma$ from Eqs.(49) and (50) are in qualitative agreement with the exact one from low- to high-temperature limits. The quantitative deviation in $\Gamma$ becomes apparent only in the high $T$ regime, with PTA value lower by a factor $0.77$ . $\Omega_{p}-\omega_{0}$ and $\Gamma$ have the same power law of $T$ , being $T^{1}$ and $T^{1/4}$ in the low and high temperature limits, respectively.

These power laws can be understood as follows. The peak position $\Omega_{p}$ is qualitatively described by the $N=1$ PTA expression

\omega_{p}=\sqrt{\omega_{0}+(12\alpha/\mu)\langle x^{2}\rangle}.

(51)

In the low temperature limit, the system tends to be harmonic, $\langle x^{2}\rangle\sim T\ll\omega_{0}$ due to equipartition theorem. In the high temperature limit, the particle’s movement is dominated by the $x^{4}$ potential. Exact analysis gives $\langle x^{2}\rangle\sim T^{1/2}\gg\omega_{0}$ [33]. Equation (51) then gives the asymptotic behavior of $\Omega_{p}(T)$ observed in the inset of Fig.5,

\Omega_{p}(T)\sim\left\{\begin{array}[]{lll}\omega_{0}+c\,T^{1},&\,\,\,(T\ll T_{cr})\\ \\ c^{\prime}\,T^{\frac{1}{4}},&\,\,\,(T\gg T_{cr}).\end{array}\right.

(52)

Here, $c$ and $c^{\prime}$ are $T$ -independent coefficients. The crossover temperature $T_{cr}\approx 1.0$ .

The spectral width $\Gamma(T)$ is generated by the thermal fluctuation of oscillating frequency which, due to the nonlinear potential, comes from the thermal fluctuation of the initial energy $E_{0}$ of the system. Guided by Eq.(51), we assume that the particle oscillates with a frequency $\omega\sim\sqrt{\omega_{0}^{2}+12\alpha x_{0}^{2}}$ , with $x_{0}$ being the initial coordinate. $x_{0}$ depends on $E_{0}$ through $E_{0}\sim\frac{1}{2}\mu\omega_{0}x_{0}^{2}$ at low $T$ , and $E_{0}\sim\alpha x_{0}^{4}$ at high $T$ . Inserting them into the expression for $\omega$ and considering $E_{0}\sim\delta E_{0}\sim T$ due to thermal fluctuation, we obtain the fluctuation of frequency $\delta\omega\sim\delta E_{0}\sim T$ at low $T$ , and $\delta\omega\sim E_{0}^{-3/4}\delta E_{0}\sim T^{1/4}$ . This explains the observed asymptotic behavior of $\Gamma(T)$ in the inset of Fig.5,

\Gamma(T)\sim\left\{\begin{array}[]{lll}d\,T^{1},&\,\,\,(T\ll T_{cr})\\ \\ d^{\prime}\,T^{\frac{1}{4}},&\,\,\,(T\gg T_{cr}).\end{array}\right.

(53)

If we further enlarge the basis $B_{2}$ by adding $\{x^{3}e^{\lambda_{i}H}\}$ , we obtain basis $B_{3}$ . Fig.6 shows the spectral function $\rho_{x,x}(\omega)$ from PTA on basis $B_{3}$ . The exact curve is shown as the black dashed line. The basis dimension of $B_{3}$ is $D=2N$ . At $N=1$ , the two delta peaks are located at the central positions of the primary and the secondary peaks of the exact curve. With increasing $N$ , PTA produces more and more delta peaks in the frequency regime of the primary and the secondary peaks, forming two broadened curves. Combining the results for $B_{2}$ and $B_{3}$ bases, we conclude that the $H$ -expanding of basis can lift the degeneracy of the quasi-particle excitation energies obtained from a lower order PTA, and split the corresponding delta peak into a bunch of densely located peaks, which effectively describe a thermally broadened peak.

IV $\mathcal{H}_{k}$ -expanded basis for the 1D nonlinear $\phi^{4}$ lattice model

In this section, we apply the $H$ -expanded basis formalism to the classical 1D nonlinear $\phi^{4}$ lattice model, and we calculate the phonon spectral function $\rho_{Q_{k},Q_{k}^{\ast}}$ , with a focus on the thermal broadening effect of the phonon peak at $T>0$ . Here, $Q_{k}=1/\sqrt{L}\sum_{j}e^{-ijk}x_{j}$ is the Fourier component of coordinates of particles. The Hamiltonian of the $\phi^{4}$ model reads

H=\sum_{i=1}^{L}\Big{[}\frac{p_{i}^{2}}{2m}+\frac{K}{2}(x_{i}-x_{i-1})^{2}+\frac{\gamma}{4}x_{i}^{4}\Big{]},

(54)

Here, $L$ is the total number of particles, $x_{i}$ represents the deviation of the $i$ -th particle from its equilibrium position, $K$ is the nearest-neighbor coupling strength, and $\gamma$ is the coefficient of the on-site potential. We use periodic boundary condition $x_{0}=x_{L}$ . This model is obtained by discretizing the classical $\phi^{4}$ field [43]. It has been widely used in the study of 1D heat transport [44, 45, 46, 47, 48, 26, 27, 49]. Similar to the FPU- $\beta$ model [50], the $\phi^{4}$ model has a scaling property which makes it sufficient to study only the parameter point $K=\gamma=1$ [33]. Below, we will focus on this parameter point.

For this model, the previous PTA study employs the basis $Q_{k}$ [33]. To describe the thermal fluctuation effect, it is tempting to expand it into the basis $A^{k}_{i}\equiv Q_{k}e^{\lambda_{i}H}$ ( $i=1,2,...,N$ ), similar to the case of one-variable AHO model. However, we find that this scheme has two problems. First, the peak width $\Gamma_{k}$ in the phonon dispersion function $\rho_{Q_{k},Q_{k}^{\ast}}(\omega)$ obtained from this basis is proportional to $1/\sqrt{L}$ , which vanishes in the large size limit $L\rightarrow\infty$ . This is unphysical. The reason is that, unlike in the one-variable AHO model, the thermal broadening of the phonon peak of $\rho_{Q_{k},Q_{k}^{\ast}}(\omega)$ for a long chain of atoms is no longer contributed by the fluctuation of total energy $H$ , but mainly by the thermal fluctuation of energy $\mathcal{H}_{k}$ of mode $k$ . Indeed, in the microscopic canonical ensemble where the total energy $H$ is fixed, the $H$ -expanded basis becomes identical to the single variable $Q_{k}$ and hence cannot produce a broadened phonon peak, while the MD simulation in the microscopic canonical ensemble can correctly produce the phonon broadened effect.

In fact, in the canonical ensemble in which our formalism is written, the thermal fluctuation $\delta H$ of the total energy $H$ , when distributed equally to each mode, gives a magnitude $\delta H/L$ to the fluctuations of each mode $H_{k}$ . Due to the nonlinear coupling, it further gives an uncertainty of the frequency $\delta\omega_{k}\sim\sqrt{\delta H/L}$ . This explains the $1/\sqrt{L}$ dependence of the incorrect results obtained from the basis $Q_{k}e^{\lambda_{i}H}$ .

The second problem of the basis $\{Q_{k}e^{\lambda_{i}H}\}$ is that the complex variable $Q_{k}$ becomes real at momenta $k=0$ and $k=\pi$ , which have distinct dynamical and statistical properties from the complex variables at $k\neq 0,\pi$ . As a result, special care must be taken in PTA formalism to avoid the possible discontinuity in physical results at $k=0$ and $\pi$ . This will complicate the formalism.

Due to the considerations above, in the following, we will work on the real basis $\text{Re}Q_{k}$ and $\text{Im}Q_{k}$ , and we expand the basis using the effective local energy $\mathcal{H}_{k}$ of the mode $k$ , instead of the total Hamiltonian $H$ . Although $\mathcal{H}_{k}$ is no longer a strictly conserving quantity like $H$ , it is approximately conserving on the level of $Q_{k}$ -PTA.

IV.1 $\{Q_{k}\}$ -PTA reformulated on a real basis

Since we will construct the $\mathcal{H}_{k}$ -expanded real basis, in this subsection, we first derive the PTA formalism on the real basis $\{\text{Re}{Q_{k}},\text{Im}{Q_{k}}\}$ . The obtained equations are equivalent to those of PTA on the complex basis $\{Q_{k}\}$ discussed in Ref.[33]. Then we expand it with $e^{\lambda_{i}\mathcal{H}_{k}}$ , where $\mathcal{H}_{k}$ is the effective energy of the real mode $k$ (see below). In this part, we reformulate this PTA in terms of the real basis, because we will need the results of this lower order PTA to identify the approximate eigen-exciatation variable in order to introduce the proper mode energy $\mathcal{H}_{k}$ and to expand the former with the latter.

We first define the following dynamical variables,

	$\displaystyle Q_{k}\equiv\frac{1}{\sqrt{L}}\sum_{j}x_{j}e^{-ijk},$
	$\displaystyle P_{k}\equiv\frac{1}{\sqrt{L}}\sum_{j}p_{j}e^{ijk},$
	$\displaystyle R_{k}\equiv\frac{1}{\sqrt{L}}\sum_{j}x^{3}_{j}e^{-ijk}.$		(55)

The lattice momentum $k$ takes $L$ different values $k=(2\pi/L)n_{k}$ ( $n_{k}=0,1,...,L-1$ ). Here and below, we assume that $L$ is an even number so that $k=\pi$ is possible. We take the real and imaginary part of $Q_{k}$ and $P_{k}$ as our dynamical variables. Considering the independence of variables, we will confine the momentum $k$ to $[0,\pi]$ , and we take the following definition for the dynamical variables

	$\displaystyle Q_{1k}\equiv\alpha_{k}\text{Re}{Q_{k}},\,\,\,\,\,\,\,k\in[0,\pi],$
	$\displaystyle Q_{2k}\equiv-\alpha_{k}\text{Im}{Q_{k}},\,\,\,\,\,\,\,k\in(0,\pi),$
	$\displaystyle P_{1k}\equiv\alpha_{k}\text{Re}{P_{k}},\,\,\,\,\,\,\,k\in[0,\pi],$
	$\displaystyle P_{2k}\equiv\alpha_{k}\text{Im}{P_{k}},\,\,\,\,\,\,\,k\in(0,\pi).$		(56)

In the above equations, the factor $\alpha_{k}$ is for normalizing the variables. $\alpha_{k}=1$ for $k=0$ and $\pi$ , and $\alpha_{k}=\sqrt{2}$ for $k\neq 0$ or $\pi$ . Note that since $Q_{k=0}$ and $Q_{k=\pi}$ are real variables, for $Q_{1k}$ and $P_{1k}$ , $k$ takes $L/2+1$ values, $k=(2\pi/L)n_{k}$ ( $n_{k}=0,1,...,L/2$ ). For $Q_{2k}$ and $P_{2k}$ , $k$ takes $L/2-1$ values, $k=(2\pi/L)n_{k}$ ( $n_{k}=1,...,L/2-1$ ). Altogether we have $L$ independent real variables.

One can prove that the above-defined variables $Q_{\nu k}$ and $P_{\nu k}$ ( $\nu=1,2$ ) fulfill the properties of canonical variables. That is, $\{Q_{\nu k},Q_{\mu p}\}=\{P_{\nu k},P_{\mu p}\}=0$ and $\{Q_{\nu k},P_{\mu p}\}=\delta_{\nu\mu}\delta_{kp}$ . They are a set of complete and independent real dynamical variables in place of $\{x_{i},p_{i}\}$ or $\{Q_{k},P_{k}\}$ .

The PTA formalism on the basis $Q_{\nu k}$ ( $k\in[0,\pi]$ for $\nu=1$ , and $k\in(0,\pi)$ for $\nu=2$ ) can be obtained from the previous PTA equations on basis $Q_{k}$ [33]. We obtain the inner product matrix $I_{\nu k,\mu p}$ and the Liouville matrix $L_{\nu k,\mu p}$ as

	$\displaystyle I_{\nu k,\mu p}=\frac{1}{m}\delta_{\nu\mu}\delta_{kp},$
	$\displaystyle L_{\nu k,\mu p}=\frac{1}{m}\omega_{k}^{2}\delta_{\nu\mu}\delta_{kp}.$		(57)

The phonon excitation energy $\omega_{k}$ is given by

\omega_{k}=\sqrt{\frac{1}{m}\left[2K\left(1-\cos{k}\right)+3\gamma\langle x^{2}\rangle\right]}.

(58)

From these expressions, one obtains the spectral functions $\rho_{Q_{\nu k},Q_{\mu p}}(\omega)$ as

\rho_{Q_{\nu k},Q_{\mu p}}(\omega)=\delta_{\nu\mu}\delta_{kp}\frac{1}{2m\omega_{k}}\left[\delta(\omega-\omega_{k})-\delta(\omega+\omega_{k})\right].

(59)

The PTA self-consistent equation for $\langle x^{2}\rangle$ is obtained from the spectral theorem Eq.(4) as

	$\displaystyle\langle x^{2}\rangle$	$\displaystyle=$	$\displaystyle\frac{1}{L}\sum_{\nu=1}^{2}\sum_{k_{\nu}}\langle Q_{\nu k_{\nu}}^{2}\rangle$		(60)
		$\displaystyle=$	$\displaystyle\frac{1}{L}\sum_{k\in[0,2\pi)}\frac{1}{\beta m\omega_{k}^{2}}.$		(60)

In the last line above, we have used the fact that $\langle Q_{1k}^{2}\rangle=\langle Q_{2k}^{2}\rangle=\langle Q_{1\,2\pi-k}^{2}\rangle$ for $k\in(0,\pi)$ . The final expression is same as that of the $\{Q_{k}\}$ basis PTA.

IV.2 effective Hamiltonian $\mathcal{H}_{\nu k}$ of the mode $(\nu k)$

In order to expand the basis $Q_{\nu k}$ using the mode energy, we first have to construct an effective energy $\mathcal{H}_{\nu k}$ for mode $(\nu k)$ . We will do this construction approximately based on the real basis $\{Q_{1k},Q_{2k}\}$ PTA equations discussed above. The effective energy $\mathcal{H}_{\nu k}$ is constructed by the following principles. First, it is a real variable with the unit of energy; second, it generates the expected dynamics, i.e., $\{Q_{\nu k},\mathcal{H}_{\nu k}\}\approx\{Q_{\nu k},H\}$ and $\{P_{\nu k},\mathcal{H}_{\nu k}\}\approx\{P_{\nu k},H\}$ ; and third, it is a zero-th order conserving quantity, i.e., $\{\mathcal{H}_{\nu k},H\}\approx 0$ . The third principle is required because the broadening of the phonon peak will be regarded in our theory as splitting the degeneracy of the excitation energy $\omega_{k}$ by high order fluctuations beyond the zero-th order approximation, as discussed in the case of the one-variable AHO model. In this work, the two last principles are fulfilled at the level of PTA on the basis $\{Q_{1k},Q_{2k}\}$ .

From Eq.(IV.1) of $\{Q_{1k},Q_{2k}\}$ -PTA and the general principle of PTA, we have the approximate EOM as

	$\displaystyle\{Q_{\nu k},H\}=\frac{1}{m}P_{\nu k},$
	$\displaystyle\{P_{\nu k},H\}\approx-m\omega_{k}^{2}\,Q_{\nu k}.\,\,\,\,\,\,\,(\nu=1,2)$		(61)

It is then easy to find the two eigen excitation variables associated to mode $(\nu k)$ by diagonalizing the EOM matrix. We obtain

	$\displaystyle a_{\nu k}^{(1)}=c_{k}\left(m\omega_{k}^{2}\,Q_{\nu k}-i\omega_{k}P_{\nu k}\right),$
	$\displaystyle a_{\nu k}^{(2)}=\tilde{c}_{k}\left(m\omega_{k}^{2}\,Q_{\nu k}+i\omega_{k}P_{\nu k}\right).$		(62)

They obey $\{a_{\nu k}^{(1)},H\}\approx i\omega_{k}\,a_{\nu k}^{(1)}$ , and $\{a_{\nu k}^{(2)},H\}\approx-i\omega_{k}\,a_{\nu k}^{(1)}$ , respectively.

One can construct the approximately conserving quantity $\mathcal{H}_{\nu k}$ as $\mathcal{H}_{\nu k}=f_{k}a_{\nu k}^{(1)\ast}a_{\nu k}^{(1)}$ . It satisfies the approximate conservation $\{\mathcal{H}_{\nu k},H\}\approx 0$ . Fixing the coefficients $c_{k}$ and $\tilde{c}_{k}$ from the second principle, i.e., $\{Q_{\nu k},\mathcal{H}_{\nu k}\}\approx\{Q_{\nu k},H\}$ , we obtain the effective Hamiltonian of mode $(\nu k)$ as

\mathcal{H}_{\nu k}=\frac{1}{2}m\omega_{k}^{2}\,Q_{\nu k}^{2}+\frac{1}{2m}P_{\nu k}^{2},

(63)

which has the desired properties

	$\displaystyle\{\mathcal{H}_{\nu k},H\}\approx 0,$
	$\displaystyle\{Q_{\nu k},\mathcal{H}_{\nu k}\}\approx\{Q_{\nu k},H\}$
	$\displaystyle\{P_{\nu k},\mathcal{H}_{\nu k}\}\approx\{P_{\nu k},H\}.$		(64)

IV.3 PTA on $\mathcal{H}_{\nu k}$ -expanded basis

Following the same practice of the PTA on the $H$ -expanded basis for the one-variable AHO model, we construct the expanded basis for the 1D $\phi^{4}$ lattice model in the subspace $(\nu k)$ as

A^{\nu k}_{i}=Q_{\nu k}e^{\lambda_{i}\mathcal{H}_{\nu k}}.\,\,\,\,\,(\lambda_{i}\,\,\,\text{real},\,\,\,\,i=1,2,...,N)

(65)

Here, $(\nu k)=(1,k_{1})$ and $(2,k_{2})$ , with $k_{1}\in[0,\pi]$ and $k_{2}\in(0,\pi)$ , respectively. We assign $\lambda_{1}=0$ to make sure $A^{\nu k}_{1}=Q_{\nu k}$ .

The inner product matrix ${\bf I}^{\nu k}$ and the Liouville matrix ${\bf L}^{\nu k}$ are derived from straightforward calculation, along the line for the AHO model. In the derivation, we have used the properties Eq.(IV.2) of $\mathcal{H}_{\nu k}$ . We obtain

$\displaystyle I^{\nu k}_{ij}$	$\displaystyle=$	$\displaystyle\frac{1}{m}\frac{\beta}{\beta-\lambda_{i}-\lambda_{j}}\langle e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\nu k}}\rangle,$
$\displaystyle L^{\nu k}_{ij}$	$\displaystyle=$	$\displaystyle\frac{1}{m}2K(1-\cos{\omega})I^{\nu k}_{ij}$	(66)
		$\displaystyle+\frac{1}{m^{2}}\frac{\beta}{\beta-\lambda_{i}-\lambda_{j}}3\gamma\langle x^{2}e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\nu k}}\rangle.$	(66)

In the above equations, like in the study of AHO model, we evaluate $\langle e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\nu k}}\rangle$ using the quadratic variational approximation

	$\displaystyle\langle e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\nu k}}\rangle$	$\displaystyle\approx$	$\displaystyle\langle e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\nu k}}\rangle_{vari}$		(67)
		$\displaystyle=$	$\displaystyle\frac{\beta}{\beta-\lambda_{i}-\lambda_{j}},$		(67)

with the variational Hamiltonian

	$\displaystyle H_{vari}=\sum_{k}H_{k},$
	$\displaystyle H_{k}=\frac{1}{2m}P_{k}^{\ast}P_{k}+\frac{1}{2}m\omega_{k}^{2}Q_{k}^{\ast}Q_{k}.$		(68)

To calculate $\langle x^{2}e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\nu k}}\rangle$ self-consistently, we write it as

	$\displaystyle\langle x^{2}e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\nu k}}\rangle$	(69)
$\displaystyle=$	$\displaystyle\frac{1}{L}\sum_{\mu p}\langle Q_{\mu p}^{2}e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\nu k}}\rangle$
$\displaystyle=$	$\displaystyle\frac{1}{L}\sum_{\mu p}\langle Q_{\mu p}^{2}e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\mu p}}\frac{e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\nu k}}}{e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\mu p}}}\rangle$
$\displaystyle\approx$	$\displaystyle\frac{1}{L}\sum_{\mu p}\left[\langle Q_{\mu p}^{2}e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\mu p}}\rangle\frac{\langle e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\nu k}}\rangle}{\langle e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\mu p}}\rangle}\right]$
$\displaystyle\approx$	$\displaystyle\frac{1}{L}\sum_{\mu p}\langle A^{\mu p}_{i}A^{\mu p}_{j}\rangle.$

In the above approximation, the momentum dependence of $\langle x^{2}e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\nu k}}\rangle$ is neglected. This makes the interaction part of the $\bf{L}$ matrix $k$ -independent, reminiscent of the local self-energy in the dynamical mean-field theory [51, 52]. Note that the additional approximation used in Eq.(69) exaggerates the coupling between different modes and tends to overestimate the broadening. See Fig.11. The final expression can be evaluated from ${\bf I}^{\nu k}$ and ${\bf L}^{\nu k}$ via

\langle A^{\nu k}_{i}A^{\nu k}_{j}\rangle\approx\frac{1}{\beta}\left[{\bf I}^{\nu k}\left({\bf L}^{\nu k}\right)^{-1}{\bf I}^{\nu k}\right]_{ij}.

(70)

From the self-consistent solution of the above equations, one can finally obtain the physical quantities of interest, ı.e., the spectral function $\rho_{Q_{k},Q_{k}^{\ast}}(\omega)$ . It is expressed by the spectral functions of real variables $\rho_{Q_{1k},Q_{1k}}(\omega)$ as

\rho_{Q_{k},Q_{k}^{\ast}}(\omega)=\left\{\begin{array}[]{lll}\rho_{Q_{1k},Q_{1k}}(\omega),&\,\,\,k\in[0,\pi]\\ \\ \rho_{Q_{1\,2\pi-k},Q_{1\,2\pi-k}}(\omega),&\,\,\,k\in(\pi,2\pi).\end{array}\right.

(71)

The spectral function $\rho_{Q_{1k},Q_{1k}}(\omega)$ is then calculated from the standard PTA formalism Eqs.(II) and (12). Note that in the case of $N=1$ , the $\mathcal{H}_{\nu k}$ -expanded basis recovers the basis $\{Q_{\nu k}\}$ , and the excitation energy is given by the $\omega_{k}$ in Eq.(58).

In the practical calculation, the above equations are combined with the zeros-removing method described in section III.C.

IV.4 spectral function $\rho_{Q_{k},Q_{k}^{\ast}}(\omega)$ for 1D $\phi^{4}$ lattice model

Here, we present the numerical results obtained from the above formula for 1D $\phi^{4}$ lattice model. We use the scheme Eqs.(47) and (48) for assigning $\{\lambda_{i}\}$ values in our calculation.

First, in Fig.7, we benchmark PTA on the $\mathcal{H}_{\nu k}$ -expanded basis by checking the convergence of its results with increasing basis dimension $N$ , taking other results as references. In Fig.7(a), we plot the dispersion $\Omega_{k}$ curve obtained from successively larger bases. The curve from the lower bound harmonic variational approximation (LH) [53] (dashed line) is shown for comparison. It is seen that as $N$ increases, the long-wavelength part of $\Omega_{k}$ shifts towards the LH curve which is expected to be more accurate, but the short-wavelength part shifts only slightly. The $\Omega_{k}$ curve converges at about $N=10$ . Since our basis is not complete in the limit $N=\infty$ , we do not expect the dispersion to converge to the exact curve. We see that the PTA with $N=1$ already produces quite accurate $\Omega_{k}$ and for this quantity only limited improvement can be gained by enlarging $N$ .

Fig.7(b) shows the relative error in $\langle x^{2}\rangle$ calculated from PTA, taking the value from a two-site cluster variational method (CVM) [54, 55] as a reference. The CVM was originally designed for Ising-like lattice models with discrete degrees of freedom and was applied to the Klein-Gordan lattice model with continuous degrees of freedom [56]. It is known that the two-site CVM produces exact static averages for 1D models with nearest-neighbour interaction [57]. So here we take the CVM results as a reference. As $N$ increases from $N=1$ , the error of PTA first decreases and then saturates at about $N=10$ . The saturated error increases with increasing temperature. This shows both the effect and the limitation of the $\mathcal{H}_{\nu k}$ -expanded basis.

In Fig.7(c), the spectral widths $\Gamma_{k=0}(T)$ and $\Gamma_{k=0.8\pi}(T)$ extracted from $\rho_{Q_{k},Q_{k}^{\ast}}(\omega)$ of PTA are compared with those from MD simulations. In MD simulations, $\Gamma_{k}(T)$ is extracted from the power spectrum using the half-area method [41]. For this purpose, we first define the complex normal variables

\displaystyle a_{k}=i\tilde{P}_{k}+\tilde{Q}_{k},

(72)

where

	$\displaystyle\tilde{P}_{k}$	$\displaystyle=\sqrt{\frac{2}{L}}\,\sum_{j=1}^{L}p_{j}\sin(jk),$		(73)
	$\displaystyle\tilde{Q}_{k}$	$\displaystyle=\sqrt{\frac{2}{L}}\,\sum_{j=1}^{L}x_{j}\sin(jk).$		(74)

are momentum and coordinate variables. We calculate $|\hat{a}_{k}(\omega)|^{2}$ , the Fourier transform of the complex normal variable $a_{k}(t)$ in Eq. (72). $\Gamma_{k}$ is then defined as the width of the frequency window around the center $\tilde{\omega}_{k}$ of the power spectrum $|\hat{a}_{k}(\omega)|^{2}$ , in which the power is equal to half of the total power. That is,

\displaystyle\frac{\displaystyle\int_{\tilde{\omega}_{k}-\Gamma_{k}/2}^{\tilde{\omega}_{k}+\Gamma_{k}/2}|\hat{a}_{k}(\omega)|^{2}\,d\omega}{\displaystyle\int_{0}^{\infty}|\hat{a}_{k}(\omega)|^{2}\,d\omega}=\frac{1}{2}.

(75)

Here, $\tilde{\omega}_{k}$ is the center frequency of the peak in $|\hat{a}_{k}(\omega)|^{2}$ . In the case of large fluctuations, this definition gives a much more robust result in the numerical simulation than the standard method does [41].

In Fig.7(c), for both $k=0$ and $k=0.8\pi$ , we see a qualitative agreement between PTA and MD results. Both curves for $k=0$ (squares) and $k=0.8\pi$ (up triangles) show qualitative agreement in all temperature regime. In particular, the asymptotic power law behaviors, i.e., $\Omega_{k}(T\rightarrow\infty)\sim T^{1/4}$ , $\Omega_{k=0}(T\rightarrow 0)\sim T^{1/3}$ , and $\Omega_{k\neq 0}(T\rightarrow 0)\sim T^{2/3}$ are observed in both results. The source of these asymptotic power laws will be analysed in Fig.12. Part of the quantitative deviation between PTA and MD results may be traced back to the truncation error and additional errors introduced in Eq.(69). We find that the present PTA on the $H_{\nu k}$ -expanded basis slightly underestimates $\Gamma_{k}$ in the high-temperature regime, and overestimate it by a factor of $5.0$ in the low- $T$ regime. This is similar to the case of the AHO model, as shown in the inset of Fig.5. Another source of the discrepancy may lie in the different definition and extraction method of $\Gamma_{k}$ . Other spectral width data from MD simulation in the literatures [25, 27, 42] are also compared with our results. We find qualitative agreement in both the $k$ and $T$ dependences. See Fig.11.

In Fig.8, we show the evolution of spectral function $\rho_{Q_{k},Q_{k}^{\ast}}(\omega)$ of 1D $\phi^{4}$ lattice model with momentum $k$ at $T=1.0$ , obtained from $\mathcal{H}_{\nu k}$ -expanded basis PTA and averaging over $200$ random realization of $\{\lambda_{i}\}$ ’s. Panels (a) and (b) show the curve on the linear and log y axis, respectively. As $k$ increases, the peak position increases and the width of the peak shrinks. The heights of the peaks do not change much, since the total weight of $\rho_{Q_{k},Q_{k}^{\ast}}(\omega)$ in $\omega>0$ is not conserved. In (b), the spectral function shows a feature similar to that of AHO model: a low frequency cut-off at frequency $\omega_{c}(k)$ . In AHO model, $\rho_{x,x}(\omega)$ has a low frequency cut-off at $\omega_{c}=\omega_{0}$ , which is well described by the $H$ -expanded PTA with random averaging. The similarity between PTA equations of the AHO model, Eqs.(III.2.1), (22), and (25), and those of the 1D $\phi^{4}$ lattice model, Eqs.(IV.3), (67), (69), and (70) can explain the low frequency cut-off in the $\phi^{4}$ lattice model and gives the cut-off frequency

\omega_{c}(k)=\sqrt{(2K/m)(1-\cos{k})}.

(76)

The inset of Fig.8(b) compares the observed cut off frequency (symbols) with this equation (line) and perfect agreement is obtained.

Fig.9 shows the evolution of the phonon spectral function $\rho_{Q_{k=0},Q_{k=0}^{\ast}}(\omega)$ with temperature. We used basis dimension $N=20$ . For a given $T$ , the spectral function curve has a single, almost symmetric peak in the positive frequency regime, located around $\omega=\Omega_{k=0}$ , with width $\Gamma_{k=0}$ . Both $\Omega_{k=0}$ and $\Gamma_{k=0}$ increases with increasing temperature. We will see from Fig.12 that both increase as $T^{1/3}$ in the low-temperature regime, while they increase as $T^{1/4}$ in the high-temperature limit. Note that the frequency with the largest $\rho_{Q_{k},Q_{k}^{\ast}}(\omega)$ is numerically very close to the $\Omega_{k}$ calculated from Eq.(49) from one-shot raw $\rho_{Q_{k}^{\ast},Q_{k}}(\omega)$ data. Therefore, in this paper, we always use the latter definition for the dispersion.

Comparing $\Omega_{k}$ to the dispersion $\omega_{k}$ Eq.(58) obtained from a single $Q_{k}$ -basis calculation (equivalent to the $\mathcal{H}_{\nu k}$ -expanded basis at $N=1$ as well as to the quadratic variation [33]), we find the same qualitative behavior and small quantitative deviation. We conclude that the phonon dispersion from $N=1$ PTA is already quite accurate and quantitative improvement of $\Omega_{k}$ with increasing $N$ is limited.

Fig.10 shows the dispersion $\Omega_{k}$ and width $\Gamma_{k}$ , extracted from $\rho_{Q_{k},Q_{k}^{\ast}}(\omega)$ using Eqs.(49) and (50), respectively. We present them as functions of $k$ for various temperatures. The dispersion $\Omega_{k}$ shown in Fig.10(a) is quantitatively close to Eq.(58), as checked in Fig.7(a). The gap at $k=0$ increases with temperature in a two-section power-law fashion, as summarized in Fig.12(a). The spectral width shown in Fig.10(b) is a decreasing function of $k$ for all $T$ , in agreement with the MD result [42]. The long wave length phonons that are important to the thermal transport are strongly damped by thermal fluctuations. This leads to a finite thermal conductivity in the thermodynamic limit at $T>0$ , in contrast to the lattice models with continuous spatial translational invariance. Note that the data in Fig.10 converge with respect to increasing chain length $L$ . To obtain reliable data for the thermodynamic limit, we used a chain length up to $L=4000$ , which is larger than the correlation length of the lowest studied temperature.

In Fig.11, we make quantitative comparisons between the results from the PTA and various MD simulations. Fig.11(a) shows the phonon life time $\tau_{k}\equiv 1/(2\Gamma_{k})$ and Fig.11(b) shows the mean free path $l_{k}\equiv\tau_{k}v_{k}$ at $T=1.0$ . Note that the extraction of $\tau_{k}$ and $l_{k}$ from MD simulation is a nontrivial task and there are various ways of doing it [25, 27, 41, 42]. This leads to quantitative deviations among different MD results. The PTA results for $\tau_{k}$ and $l_{k}$ are smaller than those of MD by a factor of $5$ or so. When multiplied by a factor $5$ , PTA result has qualitatively the same $k$ -dependence with the MD curves. We believe that, apart from the difference in the definition of $\Gamma_{k}$ , the approximation in Eq.(69) is a major source for the underestimation of $\tau_{k}$ and $l_{k}$ in PTA. In Eq.(69), the correlation between $Q_{\mu p}^{2}$ and $e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\nu k}}$ for modes $(\mu p)\neq(\nu k)$ is replaced with those of the same mode. This tends to exaggerate $\langle x^{2}e^{(\lambda_{i}+\lambda_{j})\mathcal{H}_{\nu k}}\rangle$ and hence overestimates the broadening $\Gamma_{k}$ at low temperatures (see Fig.7(c)). Further improvement of our theory in this aspect is in progress.

From the benchmark comparison with MD in Fig.7 and Fig.11, we pinpoint certain aspects where the present method of PTA can compete with and has some advantages over MD. For the low temperature limit where the equilibrium state needs to be built in a much longer time, MD is time-consuming. Both the low and high frequency limits of the dynamical correlation function are resource-demanding for MD. The methods of extracting relevant physical quantities are yet to be unified in MD, leading to quantitative deviations among different MD studies, as shown in Fig.11. In contrast, PTA produces results for well-defined quantities within the same order of time for arbitrary temperature and frequency with well-controlled precision.

Having presented the benchmarking results, below we focus on new physical results obtained by PTA for the 1D $\phi^{4}$ lattice model in Figs.12, 13, and 14. Fig.12 shows the temperature dependence of $\Omega_{k}$ and $\Gamma_{k}$ , for a series $k$ from $0$ to $\pi$ . The dispersion $\Omega_{k}(T)$ shown in Fig.12(a) is a three-section power law curve of temperature except at $k=0$ . In the high temperature regime $T\gg T_{\text{cr}}\sim 1.0$ , $\Omega_{k}(T)\sim T^{1/4}$ and the curves for all momentum merge into the curve $\Omega(T)$ of the AHO model (stars in the figure). This shows that $T\gg T_{\text{cr}}$ is a local-oscillator regime where oscillations are dominated by local $x^{4}$ potential and are dispersionless. As temperature decreases below $T_{cr}$ , $\Omega_{k}(T)\sim T^{1/3}$ , before it finally enter a plateau at very low temperature $T\ll T_{k}$ . The crossover temperature $T_{k}$ decreases when $k$ decreases. In particular, $T_{k=0}=0$ and the intermediate temperature power law $\Omega_{k=0}(T)\sim T^{1/3}$ is thus maintained to $T=0$ .

Fig.12(b) shows the spectral width $\Gamma_{k}(T)$ . It is also an increasing function with three-section power-law behavior, except at $k=0$ . Being different from $\Omega_{k}$ , $\Gamma_{k}$ decreases with increasing $k$ . It tends to zero at $T=0$ for all momentum. At high temperatures $T\gg T_{cr}\sim 1.0$ , $\Gamma_{k}(T)\sim T^{1/4}$ for all momenta and they merge into the $\Gamma(T)$ curve of AHO model (stars). As temperature lowers, $\Gamma_{k}(T)\sim T^{1/3}$ in a temperature window $T_{k}\ll T\ll T_{cr}$ . The $T_{k}$ here is the same $k$ -dependent crossover temperature in the $\Omega_{k}(T)$ curve. Finally, when $T\ll T_{k}$ , we have the third section with power law $\Gamma_{k}(T)\sim T^{2/3}$ . For $k=0$ , $T_{k=0}=0$ . So $\Gamma_{k=0}(T)\sim T^{1/3}$ is maintained to the low temperature limit. As shown in Fig.7(c), these power laws are confirmed by the MD results.

The inset of Fig.12(b) shows the ratio $\Gamma_{k}(T)/\delta\Omega_{k}(T)$ as functions of $T$ for the same $k$ ’s as in the main figure. Here $\delta\Omega_{k}(T)\equiv\Omega_{k}(T)-\Omega_{k}(T=0)$ , with $\Omega_{k}(T=0)=\sqrt{(2K/m)(1-\cos{k})}$ being the dispersion at zero temperature. It is seen that the power law $T$ -dependence in the numerator and denominator cancels, giving a ratio that is weakly dependent on $T$ and $k$ . This is consistent with the asymptotic expressions for $\Omega_{k}(T)$ and $\Gamma_{k}(T)$ in Eqs.(77) and (78).

The peculiar temperature dependence of $\Omega_{k}(T)$ and $\Gamma_{k}(T)$ can be understood in terms of the quadratic variation result for $\omega_{k}$ in Eq.(58), by a similar analysis as for the AHO model. It is well known [43, 33, 56] that for $\phi^{4}$ model, $\langle x^{2}\rangle\sim T^{2/3}$ for $T\ll T_{cr}$ and $\langle x^{2}\rangle\sim T^{1/2}$ for $T\gg T_{cr}$ . $T_{cr}$ depends on the ratio between $K$ and $\gamma$ , and $T_{cr}\sim 1.0$ in our calculation. Putting this into Eq.(58), one obtains the asymptotic $T$ dependence of dispersion as

\Omega_{k}(T)\sim\left\{\begin{array}[]{lllll}\Omega_{k}(T=0)+c\,T^{\frac{2}{3}},&\,\,\,(T\ll T_{k})\\ \\ c^{\prime}\,T^{\frac{1}{3}},&\,\,\,(T_{k}\ll T\ll T_{cr})\\ \\ c^{\prime\prime}\,T^{\frac{1}{4}},&\,\,\,(T\gg T_{cr}).\end{array}\right.

(77)

Here, $T_{k}$ is obtained by equating $2K(1-\cos{k})$ with $3\gamma\langle x^{2}\rangle$ in Eq.(58), which gives $T_{k}\sim\left[(K/\gamma)(1-\cos{k})\right]^{3/2}$ . It tends to zero in the limit $k=0$ . Therefore, we call the three regimes of temperature the local oscillator regime ( $T\gg T_{cr}$ ), the pure harmonic regime ( $T\ll T_{k}$ ), where anharmonic interaction does not appear in the dispersion $\Omega_{k}\approx\omega_{k}(\gamma=0)$ , and the renormalized harmonic regime $T_{k}\ll T\ll T_{cr}$ , where the oscillators still form collective waves but with renormalized dispersion since $3\gamma\langle x^{2}\rangle\gg\omega_{k}$ in Eq.(58).

For the thermal broadening $\Gamma_{k}(T)$ , we invoke the expression $\omega_{k}(E_{0})\sim\sqrt{\frac{1}{m}\left[2K(1-\cos{k})+3\gamma x_{0}^{2}\right]}$ , according to Eq.(58). Here $E_{0}$ is the initial energy of the mode $k$ , which has a fluctuation magnitude of $T$ . In the low-temperature limit, $x_{0}^{2}\sim\langle x^{2}\rangle\sim T^{2/3}\sim E_{0}^{2/3}$ . In the high-temperature limit, the energy is dominated by the $x^{4}$ local potential and we have $x_{0}^{4}\sim E_{0}$ . Inputting these relations into $\omega_{k}(E_{0})$ and taking $E_{0}\sim\delta E_{0}\sim T$ , we obtain the expression

\Gamma_{k}(T)\sim\left\{\begin{array}[]{lllll}d\,T^{\frac{2}{3}},&\,\,\,(T\ll T_{k})\\ \\ d^{\prime}\,T^{\frac{1}{3}},&\,\,\,(T_{k}\ll T\ll T_{cr})\\ \\ d^{\prime\prime}\,T^{\frac{1}{4}},&\,\,\,(T\gg T_{cr}).\end{array}\right.

(78)

Here, we have assumed $T_{k}\ll T_{cr}$ . Eqs.(77) and (78) summarize the observed results in Fig.12.

Eqs.(77) and (78) suggests a relation between the width $\Gamma_{k}(T)$ and dispersion $\Omega_{k}(T)$ ,

\Gamma_{k}(T)\approx c\left[\Omega_{k}(T)-\Omega_{k}(0)\right].

(79)

Here, $c$ is weakly dependent on $T$ and $k$ . See the inset of Fig.12(b) for its behavior. This relation is in contrast to the assumption of the effective phonon theory [24] that $\Gamma_{k}(T)\approx\epsilon(T)\Omega_{k}(T)$ . Here, $\epsilon(T)=\langle E_{n}\rangle/\langle E_{n}+E_{l}\rangle$ , with $E_{n}$ and $E_{l}$ being the nonlinear and linear potential energies, respectively. The results $\Gamma_{k}(T)$ and $\Omega_{k}(T)$ shown in Fig.12 do not agree with this assumption both in the momentum and in the temperature dependences.

IV.5 short-chain limit of 1D $\phi^{4}$ lattice model

So far, we have focused solely on the thermodynamical limit where the chain length $L$ is much larger than the correlation length $\xi(T)$ . Our results presented above are obtained with sufficiently large $L$ to avoid any finite-size effect. Here, we report that the $\phi^{4}$ model has also a short-chain limit characterized by $L\ll\xi(T)$ , where physical quantities have distinct power laws from those in the thermodynamical limit. This regime is pertinent to the thermal transport experiments with small size low-dimensional systems, but to our knowledge, has not been discussed in the literature.

In Fig.13, $\langle x^{2}\rangle(T)$ is plotted for chain size ranging from $L=2$ to $L=2000$ . For each $L$ , in the high temperature limit $T\gg T_{cr}\approx 1.0$ , we have $\langle x^{2}\rangle(T)\sim T^{1/2}$ , the power law of the atomic limit discussed in previous subsections. As $T$ decreases below $T_{cr}$ , $\langle x^{2}\rangle\sim T^{2/3}$ . Here $T_{cr}$ is the $L$ -independent local-nonlocal crossover temperature. Since the correlation length $\xi$ increases with decreasing $T$ as $\xi\sim T^{-1/3}$ in the low temperature limit [56] (see the inset of Fig.13), when $T$ decreases further and is below a short-chain crossover temperature $T_{sh}$ where $\xi=L$ (marked by the arrows in the figure), $\xi>L$ occurs and $\langle x^{2}\rangle(T)$ begins to deviate from the $T^{2/3}$ law. In the low $T$ limit where $L\ll\xi(T)$ , $\langle x^{2}\rangle\sim T^{1/2}$ , being different from the $T^{2/3}$ behavior of the thermodynamical limit. We have checked that the dependence of $T_{sh}$ on $L$ is determined by $\xi(T_{sh})=L$ , giving $T_{sh}(L)\sim L^{-3}$ , with a coefficient about $1.0$ . In the thermodynamical limit, $T_{sh}(L=\infty)=0$ and the $T^{2/3}$ behavior extends to zero temperature. Note that $T_{sh}(L)\leq T_{cr}\approx 1.0$ holds for its largest value $L=2$ . This behavior is summarized as

\langle x^{2}\rangle(T)\sim\left\{\begin{array}[]{lllll}T^{\frac{1}{2}},&\,\,\,(T\ll T_{sh}(L))\\ \\ T^{\frac{2}{3}},&\,\,\,(T_{sh}(L)\ll T\ll T_{cr})\\ \\ T^{\frac{1}{2}},&\,\,\,(T\gg T_{cr}).\end{array}\right.

(80)

This behavior can be understood from the $N=1$ solution, where $\langle x^{2}\rangle$ is determined by the self-consistent equations Eqs.(58) and (60). In the limit of small $L$ and low $T$ , the summation of $k$ in Eq.(60) is dominated by the $k\sim 0$ terms. The result is therefore similar to that of the atomic limit at $K=0$ or $T\gg T_{cr}$ .

A similar short-chain limit exists in the dynamical quantities $\Omega_{k}(T)$ and $\Gamma_{k}(T)$ . In Fig.14(a) and (b), we show the curves of $\Omega_{k=0}(T)$ and $\Gamma_{k=0}(T)$ for a series $L$ . For momentum $k=0$ , both quantities have a common three-section power-law behavior. As $T$ decreases, they change from the high temperature atomic-limit behavior $T^{1/4}$ , to the non-local long-chain behavior $T^{1/3}$ in $T_{sh}(L)\ll T\ll T_{cr}$ , and to the short-chain behavior $T^{1/4}$ in the low temperature limit $T\ll T_{sh}(L)\leq T_{cr}$ . The behavior is summarized as

\Omega_{k=0}(T),\,\,\,\Gamma_{k=0}(T)\sim\left\{\begin{array}[]{lllll}T^{\frac{1}{4}},&\,\,\,(T\ll T_{sh}(L))\\ \\ T^{\frac{1}{3}},&\,\,\,(T_{sh}(L)\ll T\ll T_{cr})\\ \\ T^{\frac{1}{4}},&\,\,\,(T\gg T_{cr}).\end{array}\right.

(81)

In Fig.14(c) and (d), we show the quantities for $k=\pi$ . In the low $T$ limit, $\Omega_{k=\pi}(T)$ is a finite constant, not showing the short-chain behaviors. However, the quantities $\delta\Omega_{k=\pi}(T)\equiv\Omega_{k=\pi}(T)-\Omega_{k=\pi}(0)$ and $\Gamma_{k=\pi}(T)$ do have short-chain limit behavior. For other $k>0$ , similar behaviors are found and they are summarized as

\delta\Omega_{k>0}(T),\,\,\,\Gamma_{k>0}(T)\sim\left\{\begin{array}[]{lllll}T^{\frac{1}{2}},&\,\,\,(T\ll T_{sh}(L))\\ \\ T^{\frac{2}{3}},&\,\,\,(T_{sh}(L)\ll T\ll T_{cr})\\ \\ T^{\frac{1}{4}},&\,\,\,(T\gg T_{cr}).\end{array}\right.

(82)

These results are in principle amenable to experimental tests on low dimensional materials with varying length and temperature. They also show that one must be careful when interpreting the temperature behavior of the experimental or MD data related to thermal transport in low dimensional systems, since at least two energy scales are present to separate the scaling regimes. One is related to the length of chain, and the other is related to the competition between local and non-local interactions. In the present work, we fix parameters $m=K=\gamma=1$ and hence have a definite $T_{cr}\approx 1.0$ and $T_{sh}(L)\approx 1.0\times L^{-3}$ . For a system with arbitrary parameters, the concrete values of the two crossover temperatures will shift according to the exact or approximate scaling properties of the model. Therefore before one can carry out reliable power-law fitting, one must collect data in a sufficiently wide range of temperature to discern the possible crossovers.

V Discussions and Summary

We first discuss the implication of our results for the temperature dependence of thermal conductivity $\kappa(T)$ . In the Debye formula $\kappa=\sum_{k}C_{k}v_{k}l_{k}$ [31, 32], $\kappa$ is expressed in terms of the specific heat capacity $C_{k}=(1/L)\partial\langle H_{k}\rangle/\partial T$ , the phonon group velocity $v_{k}=\partial\Omega_{k}/\partial k$ , and the mean free path $l_{k}=v_{k}/(2\Gamma_{k})$ .

PTA results for $\Omega_{k}$ and $\Gamma_{k}$ can be used to calculate $\kappa$ through the Debye formula. It is known that $C_{k>0}$ ranges between $1.0$ (low $T$ ) and $0.75$ (high $T$ ) and is weakly $k$ dependent, and $C_{k=0}=0.75$ [33]. Approximating $C_{k}$ as a constant, and inputting Eqs.(77) and (78) into the Debye formula, we obtain the asymptotic $T$ dependence of $\kappa$ as

\kappa(T)\sim\left\{\begin{array}[]{lll}T^{-\frac{2}{3}},&\,\,\,(T\ll T_{cr})\\ \\ T^{-\frac{3}{4}},&\,\,\,(T\gg T_{cr}),\end{array}\right.

(83)

with $T_{cr}\approx 1.0$ . Direct numerical evaluation of the Debye formula using PTA input confirms this result. However, this behavior does not agree with the high temperature result $\kappa(T)\sim T^{-1.35}$ from MD simulation in Ref.[44] and the expression $\kappa(T)\sim T^{-4/3}$ from the effective phonon theory [48]. This poses question on the applicability of the Debye formula for describing the thermal conductivity of the $\phi^{4}$ lattice model. A detailed answer to this question is beyond the scope of this paper and we leave it for a future study.

In this work, we developed the idea of expanding a low level eigen-operator $Q_{\nu k}$ by multiplying it with functions of approximate mode energy $\mathcal{H}_{\nu k}$ , to describe the thermal-broadened quasi-particle peak in the spectral function. Here, we have used the basis variable $Q_{\nu k}e^{\lambda_{i}\mathcal{H}_{\nu k}}$ , but it is interesting to know the importance of the inter-mode coupling for the effectiveness of basis. For example, PTA on the more complete basis $Q_{\nu k}e^{\sum_{i\mu p}\lambda^{i}_{\mu p}\mathcal{H}_{\mu p}}$ needs to be investigated to answer this question. We expect that it can give quantitative improvement of $\Gamma_{k}$ . It is also noted that the present $H$ -expanded basis does not cover all types of basis operators or scattering processes. For example, the variable $Q_{\nu_{1}k_{1}}Q_{\nu_{2}k_{2}}Q_{\nu_{3}k_{3}}$ ( $k_{1}+k_{2}+k_{3}=k\mod{2\pi})$ with different $k_{1}\sim k_{3}$ is not covered. Scrutinization is therefore needed for applying the present $H$ -expansion of the basis to general situations.

For quantum systems, quantum fluctuation is an important source of quasi-particle broadening. In particular, in systems such as frustrated quantum antiferromagnetism or system close to quantum criticality, quantum damping plays a dominant role. The present scheme can be generalized to study the quantum damping of quasi particles or quantum broadening of the quasi-particle peak in the spectral function. A direct application of the present formalism needs additional care, because for quantum systems the conservation identity Eqs.(84) and (85), used to simplify the matrices ${\bf I}$ and ${\bf L}$ , must be modified. The GF and the whole formalism of the PTA needs to be replaced with the quantum version. Other than this, there seem to be no other constraints. In case the scattering processes other than those contained in the $H$ -expanded basis play a vital role, one must consider the expanding schemes with other (approximately) conserved quantities other than energy, such as the particle number, momentum, spin, etc., to describe the decay of quasi-particles due to various fluctuations. What we have in mind is the spin-wave damping in frustrated antiferromagnets [40]. Further study of this interesting subject is in progress.

In summary, we developed the PTA on the $H$ -expanded basis to describe the thermal broadening of quasi particles. A zeros-removing technique is used to stabilize the iterative solution of PTA equations. Benchmarking calculations on the classical AHO model and the 1D $\phi^{4}$ model show that this method can produce semi-quantitative thermal broadening of a quasi-particle peak in the spectral function. Using this method, we disclose the low- and high-temperature power laws of the phonon spectral width $\Gamma_{k}$ for the 1D $\phi^{4}$ lattice model, contradicting the assumption in effective phonon theory and raising questions on the applicability of the Debye formula for this model. A short-chain limit of this model is also discovered. Future development and possible extension of this method are discussed.

VI Acknowledgments

We acknowledge helpful discussions with Y.J. Wu, C.L. Ji, X.G. Ren, L. Wan, Y. Wan, and T. Li. We are grateful to Professor Ren for showing us the method to treat the redundant basis variables. This work is supported by NSFC (Grant No.11974420 and No.12075316).

Appendix A Derivation of ${\bf I}$ and ${\bf L}$ matrices for AHO model on bases $B_{2}$ and $B_{3}$

In this Appendix, we present the derivation of ${\bf I}$ and ${\bf L}$ matrices for AHO on the bases $B_{2}$ and $B_{3}$ .

Before presenting the details, we present the identity to be used for simplifying the expressions [33]. For arbitrary variables $A$ and $B$ , we have

\langle\{A,B\}\rangle=\beta\langle\{A,H\}B\rangle,

(84)

and for arbitrary number $\theta$ ,

\langle\{A,H\}Be^{\theta H}\rangle=\frac{1}{\beta-\theta}\langle\{A,B\}e^{\theta H}\rangle.

(85)

Eq.(84) is proved in Appendix B. Eq.(85) is obtained by letting $K=1$ , $X_{1}=H$ , and $g(X_{1})=e^{\theta H}$ in the identity Eq.(99).

A.1 Basis $B_{2}$

For basis $B_{2}=\{A_{i}=xf_{i}\}$ ( $f_{i}=e^{\lambda_{i}H}$ ), the matrix elements $I_{ij}=(xf_{i}|xf_{j})$ is calculated through Eq.(84) as

	$\displaystyle I_{ij}$	$\displaystyle=$	$\displaystyle\langle\{xf_{i},\{xf_{j},H\}\}\rangle$		(86)
		$\displaystyle=$	$\displaystyle\beta\langle\{x,H\}\{x,H\}f_{i}f_{j}\rangle.$		(86)

Using Eq.(85) (letting $A=x$ , $B=\{x,H\}$ , and $\theta=\lambda_{i}+\lambda_{j}$ ), and $\{x,\{x,H\}\}=1/\mu$ , we obtain

	$\displaystyle I_{ij}$	$\displaystyle=$	$\displaystyle\langle\{x,\{x,H\}\}f_{i}f_{j}\rangle+(\lambda_{i}+\lambda_{j})\langle\{x,H\}\{x,H\}f_{i}f_{j}\rangle$		(87)
		$\displaystyle=$	$\displaystyle\frac{\beta}{\mu\left(\beta-\lambda_{i}-\lambda_{j}\right)}\langle f_{i}f_{j}\rangle.$		(87)

The matrix element $L_{ij}=-(xf_{i}|\{\{xf_{j},H\},H\})$ is calculated similarly as

	$\displaystyle L_{ij}$	$\displaystyle=$	$\displaystyle-\langle\{xf_{i},\{\{\{x,H\},H\},H\}f_{j}\}\rangle$		(88)
		$\displaystyle=$	$\displaystyle-\beta\langle\{x,H\}\{\{\{x,H\},H\},H\}f_{i}f_{j}\rangle.$		(88)

Using Eq.(85) (letting $A=x$ , $B=\{\{\{x,H\},H\},H\}$ , and $\theta=\lambda_{i}+\lambda_{j}$ ) and $\{x,\{\{\{x,H\},H\},H\}\}=-(1/\mu^{2})\left(\mu\omega_{0}^{2}+12\alpha x^{2}\right)$ , we obtain

\displaystyle L_{ij}

\displaystyle=

\displaystyle\frac{\beta}{\mu^{2}\left(\beta-\lambda_{i}-\lambda_{j}\right)}\left[\mu\omega_{0}^{2}\langle f_{i}f_{j}\rangle+12\alpha\langle x^{2}f_{i}f_{j}\rangle\right].

Eqs.(87) and (A.1) are Eq.(III.2.1) in the main text.

A.2 Basis $B_{3}$

For basis $B_{3}$ , the basis variables are denoted in Eq.(26) as $A_{\mu i}=a_{\mu}f_{i}$ , with $a_{1}=x$ , $a_{2}=x^{3}$ , and $f_{i}=e^{\lambda_{i}H}$ . The matrix elements of ${\bf I}$ and ${\bf L}$ read

I_{\mu\nu}^{ij}=\langle\{A_{\mu i},\{A_{\nu j},H\}\}\rangle,

(90)

and

L_{\mu\nu}^{ij}=-\langle\{A_{\mu i},\{\{\{A_{\nu j},H\},H\},H\}\}\rangle,

(91)

respectively. Using the same method used for basis $B_{2}$ above, we obtain

	$\displaystyle I_{\mu\nu}^{ij}$	$\displaystyle=$	$\displaystyle\frac{\beta}{\beta-\lambda_{i}-\lambda_{j}}\langle I_{\mu\nu}^{(0)}f_{i}(H)f_{j}(H)\rangle,$
	$\displaystyle L_{\mu\nu}^{ij}$	$\displaystyle=$	$\displaystyle\frac{\beta}{\beta-\lambda_{i}-\lambda_{j}}\langle L_{\mu\nu}^{(0)}f_{i}(H)f_{j}(H)\rangle.$		(92)

Here, the variables $I^{(0)}_{\mu\nu}=\{a_{\mu},\{a_{\nu},H\}\}$ and $L^{(0)}_{\mu\nu}=-\{a_{\mu},\{\{\{a_{\nu},H\},H\},H\}\}$ . Calculation gives

	$\displaystyle I_{11}^{(0)}=\frac{1}{\mu},$
	$\displaystyle I_{12}^{(0)}=I^{(0)}_{21}=\frac{3}{\mu}x^{2},$
	$\displaystyle I_{22}^{(0)}=\frac{9}{\mu}x^{4},$		(93)

and

	$\displaystyle L_{11}^{(0)}=\frac{1}{\mu^{2}}V^{{\prime}{\prime}}(x),$
	$\displaystyle L_{12}^{(0)}=-\frac{18}{\mu^{3}}p^{2}+\frac{18}{\mu^{2}}xV^{\prime}(x)+\frac{3}{\mu^{2}}x^{2}V^{\prime\prime}(x),$
	$\displaystyle L_{21}^{(0)}=\frac{3}{\mu^{2}}x^{2}V^{\prime\prime}(x),$
	$\displaystyle L_{22}^{(0)}=-\frac{54}{\mu^{3}}x^{2}p^{2}+\frac{54}{\mu^{2}}x^{3}V^{\prime}(x)+\frac{9}{\mu^{2}}x^{4}V^{\prime\prime}(x).$

Here $V(x)$ is the potential function Eq.(III). Inserting these equations into Eq.(A.2), we get the expressions for $I_{\mu\nu}^{ij}$ and $L_{\mu\nu}^{ij}$ . They contain the averages $\langle f_{i}f_{j}\rangle$ , $\langle x^{2}f_{i}f_{j}\rangle$ , $\langle x^{4}f_{i}f_{j}\rangle$ , $\langle p^{2}f_{i}f_{j}\rangle$ , and $\langle x^{2}p^{2}f_{i}f_{j}\rangle$ . The first three can be calculated self-consistently. The last two can be simplified by the conservation identity Eq.(85).

In Eq.(85), we let $A=x$ , $B=p$ , and $\theta=\lambda_{i}+\lambda_{j}$ , we obtain

\langle p^{2}f_{i}f_{j}\rangle=\frac{\mu}{\beta-\lambda_{i}-\lambda_{j}}\langle f_{i}f_{j}\rangle.

(95)

Similarly, letting $A=x^{2}$ and $B=xp$ , we obtain

\langle x^{2}p^{2}f_{i}f_{j}\rangle=\frac{\mu}{\beta-\lambda_{i}-\lambda_{j}}\langle x^{2}f_{i}f_{j}\rangle.

(96)

Using Eqs.(95) and (96), we finally obtain Eqs.(29) and (III.2.2). Note that the conservation relation Eq.(85) is also used to reduce $L^{ij}_{12}$ to $L^{ij}_{21}$ . The final result for ${\bf L}$ is Hermitian.

Appendix B Proof of a conservation identity

In this Appendix, we derive an identity to be used in the derivation of expressions for ${\bf I}$ and ${\bf L}$ in Appendix A.

First, we cite the conservation identity Eq.(15) of Ref.[33], i.e., for arbitrary variables $A$ and $B$ , we have

	$\displaystyle\langle\{A(t),B(t^{\prime})\}\rangle$	$\displaystyle=$	$\displaystyle\beta\langle\{A,H\}(t)B(t^{\prime})\rangle$		(97)
		$\displaystyle=$	$\displaystyle-\beta\langle A(t)\{B,H\}(t^{\prime})\rangle.$		(97)

In the above equation, the average is defined as the canonical ensemble average in an equilibrium state at temperature $T$ . $\beta=1/k_{B}T$ is the inverse temperature. In the equal time situation $t=t^{\prime}$ , we get

	$\displaystyle\langle\{A,B\}\rangle$	$\displaystyle=$	$\displaystyle\beta\langle\{A,H\}B\rangle$		(98)
		$\displaystyle=$	$\displaystyle-\beta\langle A\{B,H\}\rangle.$		(98)

In the first equation of Eq.(98), we replace $B$ with $Bg(X_{1},X_{2},...,X_{K})$ . Here $X_{i}$ is a conserved quantity satisfying $\{X_{i},H\}=0$ ( $i=1,2,...,K$ ) and $g$ is an arbitrary $K$ -variable function. We get the following conservation identity,

	$\displaystyle\langle\{A,B\}g(X_{1},X_{2},...,X_{K})\rangle$	(99)
$\displaystyle=$	$\displaystyle\beta\langle\{A,H\}Bg(X_{1},X_{2},...,X_{K})\rangle$
	$\displaystyle-\sum_{i=1}^{K}\langle\{A,X_{i}\}B\frac{\partial}{\partial X_{i}}g(X_{1},X_{2},...,X_{K})\rangle.$

Eq.(85) is a special case of Eq.(99), with $K=1$ , $X_{1}=H$ , and $g(X_{1})=e^{\theta H}$ .

References

[1] H. Mori, Transport, collective motion, and Brownian motion, Prog. Theor. Phys. 33, 423 (1965).
[2] H. Mori, A continued-fraction representation of the time-correlation functions, Prog. Theor. Phys. 34, 399 (1965).
[3] R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford University Press, New York, 2001).
[4] Yu. A. Tserkovnikov, A method of solving infinite systems of equations for two-time thermal Green’s functions, Theor. Math. Phys. 49, 993 (1981).
[5] Yu. A. Tserkovnikov, Two-time temperature Green’s functions in kinetic theory and molecular hydrodynamics: I. The chain of equations for the irreducible functions, Theor. Math. Phys. 118, 85 (1999).
[6] L. M. Roth, New method for linearizing many-body equations of motion in statistical mechanics, Phys. Rev. Lett. 20, 1431 (1968).
[7] L. M. Roth, Electron correlation in narrow energy bands. I. The two-pole approximation in a narrow s band, Phys. Rev. 184, 451 (1969).
[8] For a review, see A. Avella, Composite operator method analysis of the underdoped cuprates puzzle, Adv. Condens. Matter Phys. 2014, 515698 (2014).
[9] L. Haurie et al., Bands renormalization and superconductivity in the strongly correlated Hubbard model using composite operators method, J. Phys.: Condens. Matter 36, 255601 (2024).
[10] Y. Kakehashi and P. Fulde, Self-consistent projection operator approach to nonlocal excitation spectra, Phys. Rev. B 70, 195102 (2004).
[11] S. Onoda and M. Imada, Operator projection method applied to the single-particle Green’s function in the Hubbard model, J. Phys. Soc. Jpn 70, 632 (2001).
[12] M. H. Lee, Orthogonalization process by recurrence relations, Phys. Rev. Lett. 49, 1072 (1982).
[13] M. H. Lee, Solutions of the generalized Langevin equation by a method of recurrence relations, Phys. Rev. B 26, 2547 (1982).
[14] A. L. Kuzemsky, Irreducible Green functions method and many-particle interacting systems on a lattice, Riv. Nuovo Cimento 25, 1 (2002).
[15] S. Lepri, Relaxation of classical many-body Hamiltonians in one dimension, Phys. Rev. E 58, 7165 (1998).
[16] A. Castellano, J. P. A. Batista, and M. J. Verstraete, Mode-coupling theory of lattice dynamics for classical and quantum crystals, J. Chem. Phys. 159, 234501 (2023).
[17] A. Eskandari-asl and A. Avella, Time-resolved ARPES signal in pumped semiconductors within the dynamical projective operatorial approach, Phys. Rev. B 110, 094309 (2024).
[18] P. Fan, K. Yang, K. H. Ma, and N. H. Tong, Projective truncation approximation for equations of motion of two-time Green’s functions, Phys. Rev. B 97, 165140 (2018).
[19] P. Fan and N. H. Tong, Controllable precision of the projective truncation approximation for Green’s functions, Chin. Phys. B 28, 047102 (2019).
[20] K. H. Ma and N. H. Tong, Interacting spinless fermions on the square lattice: Charge order, phase separation, and superconductivity, Phys. Rev. B 104, 155116 (2021).
[21] P. Fan, N. H. Tong, and Z. G. Zhu, Spatial spin-spin correlations of the single-impurity Anderson model with a ferromagnetic bath, Phys. Rev. B 106, 155130 (2022).
[22] For a review, see Y. G. Rudoy, The Bogoliubov–Tyablikov Green’s function method in the quantum theory of magnetism, Theor. Math. Phys. 168, 1318 (2011).
[23] Y. Onodera, Dynamic susceptibility of classical anharmonic oscillator - A unified oscillator model for order-disorder and displacive ferroelectrics, Prog. Theor. Phys. 44, 1477 (1970).
[24] N. Li, P. Tong, and B. Li, Effective phonons in anharmonic lattices: Anomalous vs. normal heat conduction, Europhys. Lett. 75, 49 (2006).
[25] S. Liu, J. Liu, P. H $\ddot{a}$ nggi, C. Wu, and B. Li, Triggering waves in nonlinear lattices: Quest for anharmonic phonons and corresponding mean-free paths, Phys. Rev. B 90, 174304 (2014).
[26] L. Xu and L. Wang, Dispersion and absorption in one-dimensional nonlinear lattices: A resonance phonon approach, Phys. Rev. E 94, 030101(R) (2016).
[27] L. Xu and L. Wang, Resonance phonon approach to phonon relaxation time and mean free path in one-dimensional nonlinear lattices, Phys. Rev. E 95, 042138 (2017).
[28] J. S. Wang, J. Wang, and N. Zeng, Nonequilibrium Green’s function approach to mesoscopic thermal transport, Phys. Rev. B 74, 033408 (2006).
[29] R. G. Della Valle and P. Procacci, Computer-aided series expansion for phonon self-energy, J. Comp. Phys. 165, 428 (2000).
[30] B. Nickel, The solution to the 4-phonon Boltzmann equation for a 1D chain in a thermal gradient, J. Phys. A 40, 1219 (2007).
[31] J. M. Ziman, Electrons and Phonons (Oxford at the Clarendon Press, 1960), Chaps. 7 and 8.
[32] P. B. Allen and J. L. Feldman, Thermal conductivity of disordered harmonic solids, Phys. Rev. B 48, 12581 (1993).
[33] K. H. Ma, Y. J. Guo, L. Wang, and N. H. Tong, Projective-truncation-approximation study of the one-dimensional $\phi^{4}$ lattice model, Phys. Rev. E 106, 014110 (2022).
[34] J. C. Herzel, Green’s functions and double-time distribution functions in classical statistical mechanics, J. Math. Phys. 8, 1650 (1967).
[35] M. J. Smith, Classical Green’s functions for the ideal gas, Phys. Rev. Lett. 24, 1398 (1970).
[36] A. Cavallo, F. Cosenza, and L. De. Cesare, Classical Heisenberg ferromagnetic chain with long-range interactions: A spectral density approach, Phys. Rev. B 66, 174439 (2002).
[37] For a review, see A. Cavallo, F. Cosenza, and L. De. Cesare, in New Developements in Ferromagnetism Research, edited by V. N. Murray (Nova Science 2005), Chap 6.
[38] T. Dauxois, M. Peyrard, and A. R. Bishop, Dynamics and thermodynamics of a nonlinear model for DNA denaturation, Phys. Rev. E 47, 684 (1993).
[39] X. G. Ren (private discussion).
[40] S. Winterfeldt and D. Ihle, Spin dynamics and antiferromagnetic short-range order in the two-dimensional Heisenberg model, Phys. Rev. B 59, 6010 (1999).
[41] Y. J. Guo, Y. C. Sun, and L. Wang, Thermalization process in a one-dimensional lattice with two-dimensional motions: The role of angular momentum conservation, Phys. Rev. E 108, 014127 (2023).
[42] Y. Liu and D. He, Approach to phonon relaxation time and mean free path in nonlinear lattices, Chin. Phys. Lett. 38, 044401 (2021).
[43] D. Boyanovsky, C. Destri, and H. J. de Vega, Approach to thermalization in the classical $\phi^{4}$ theory in $1+1$ dimensions: Energy cascades and universal scaling, Phys. Rev. D 69, 045003 (2004).
[44] K. Aoki and D. Kusnezov, Bulk properties of anharmonic chains in strong thermal gradients: non-equilibrium $\phi^{4}$ theory, Phys. Lett. A 265, 250 (2000).
[45] B. Hu, B. Li, and H. Zhao, Heat conduction in one-dimensional nonintegrable systems, Phys. Rev. E 61, 3828 (2000).
[46] K. Aoki, J. Lukkarinen, and H. Spohn, Energy transport in weakly anharmonic chains, J. Stat. Phys. 124, 1105 (2006).
[47] A. Dhar, Heat transport in low-dimensional systems, Adv Phys. 57, 457 (2008).
[48] N. Li and B. Li, Scaling of temperature-dependent thermal conductivities for one-dimensional nonlinear lattices, Phys. Rev. E 87, 042125 (2013).
[49] N. Li, J. Liu, C. Wu, and B. Li, Temperature and frequency dependentmean free paths of renormalized phonons in nonlinear lattices, New J. Phys. 20, 023006 (2018).
[50] K. Aoki and D. Kusnezov, Fermi-Pasta-Ulam $\beta$ Model: Boundary jumps, Fourier’s law, and scaling, Phys. Rev. Lett. 86, 4029 (2001).
[51] W. Metzner and D. Vollhardt, Correlated lattice fermions in $d=\infty$ dimensions, Phys. Rev. Lett. 62, 324 (1989).
[52] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-ﬁeld theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
[53] J. Liu, S. Liu, N. Li, B. Li, and C. Wu, Renormalized phonons in nonlinear lattices: A variational approach, Phys. Rev. E 91, 042910 (2015).
[54] R. Kikuchi, A theory of cooperative phenomena, Phys. Rev. 81, 988 (1951).
[55] For a review, see Theory and Applications of the Cluster Variation and Path Probability Methods, edited by J. L. Morán-López and J. M. Sanchez (Plenum, New York 1996).
[56] H. W. Jia and N. H. Tong, Thermodynamics of classical one-dimensional Klein-Gordon lattice model, arXiv:2410.19249 (unpublished).
[57] A. Pelizzolla, Cluster variation method in statistical physics and probabilistic graphical models, J. Phys. A 38, R309 (2005).

Thermal Broadening of Phonon Spectral Function in Classical Lattice Models: Projective Truncation Approximation

Abstract

I Introduction

II Green’s Function EOM and PTA for classical systems

III One-Variable Anharmonic Oscillator Model

III.1 basis B1≡{x2​m−1​p2​n−2}B_{1}\equiv\{x^{2m-1}p^{2n-2}\}

III.2 HH-expanded basis B2≡{x​eλi​H}B_{2}\equiv\{xe^{\lambda_{i}H}\} and B3≡{x,x3}⊗{eλi​H}B_{3}\equiv\{x,x^{3}\}\otimes\{e^{\lambda_{i}H}\}

III.2.1 basis B2≡{x​eλi​H}B_{2}\equiv\{xe^{\lambda_{i}H}\}

III.2.2 basis B3≡{x,x3}⊗{eλi​H}B_{3}\equiv\{x,x^{3}\}\otimes\{e^{\lambda_{i}H}\}

III.3 removing the redundant basis variables

III.4 results for one-variable AHO model

IV ℋk\mathcal{H}_{k}-expanded basis for the 1D nonlinear ϕ4\phi^{4} lattice model

IV.1 {Qk}\{Q_{k}\}-PTA reformulated on a real basis

IV.2 effective Hamiltonian ℋν​k\mathcal{H}_{\nu k} of the mode (ν​k)(\nu k)

IV.3 PTA on ℋν​k\mathcal{H}_{\nu k}-expanded basis

IV.4 spectral function ρQk,Qk∗​(ω)\rho_{Q_{k},Q_{k}^{\ast}}(\omega) for 1D ϕ4\phi^{4} lattice model

IV.5 short-chain limit of 1D ϕ4\phi^{4} lattice model