{'fraction_non_alphanumeric': 0.04615598982826488, 'fraction_numerical': 0.011460619361131432, 'mean_word_length': 4.330115113256591, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Fokker-Planck Equation for Lattice Vibration: Stochastic Dynamics and Thermal Conductivity

We propose a Fokker-Planck equation (FPE) theory to describe stochastic fluctuation and relaxation processes of lattice vibration at a wide range of conditions, including those beyond the phonon gas (PG) limit. Using the time-dependent, multiple state-variable probability function of a FPE, we first derive time-correlation functions of lattice heat currents in terms of correlation functions among multiple vibrational modes, and subsequently predict the thermal conductivity based on the Green-Kubo formalism. When the kinetic transport theories are valid, this FPE not only predicts a thermal conductivity that is identical to the one predicted by the phonon Boltzmann transport equation, but also provides additional microscopic details on the multiple-mode correlation functions. More importantly, when the kinetic theories become insufficient due to the breakdown of the PG approximation, this FPE theory remains valid to study the correlation functions among vibrational modes in highly anharmonic lattices with significant mode-mode interactions and/or in disordered lattices with strongly localized modes. At the limit of weak mode-mode interactions, we can adopt quantum perturbation theories to derive the drift/diffusion coefficients based on the anharmonicity data derived from first-principles methods. As temperature elevates to the classical regime, we can perform molecular dynamics simulations to directly compute the coefficients. Because these coefficients are defined as ensemble averages at the short time limit, we can implement massive parallel simulation algorithms to take full advantage of the paralleled high-performance computing platforms. A better understanding of the temperature-dependent coefficients up to melting temperatures will provide new insights on microscopic mechanisms that govern the heat conduction through anharmonic and/or disordered lattices beyond the PG model.


I. INTRODUCTION
The phonon Boltzmann transport equation (BTE) 1,2 has gained some renewed interests as the default choice of transport theory to compute lattice thermal conductivity (κ Latt ) of crystalline solids from first-principles [3][4][5][6][7][8][9] . The theoretical foundation of the phonon BTE is the so-called phonon gas (PG) model [10][11][12][13] , which assumes that interactions among vibrational modes are weak enough that the numbers of phonons of each mode follow the single-particle Bose-Einstein distribution at equilibrium. As a kinetic theory, the phonon BTE further assumes that (1) each quasi-particle phonon travels at a group velocity v g , and (2) the lifetime τ of every phonon is finite because of the scatterings by lattice anharmonicity, lattice defects/disorder, or other particles. For electronic insulators, the necessary inputs for a phonon BTE calculation are the harmonic phonon spectra and the phonon scattering terms, both of which can be numerically calculated using first-principles methods 4, [14][15][16][17][18][19][20][21][22] . Multiple implementations of the phonon BTE methods have been reported in recent years [23][24][25][26][27] , and the calculated results adopting various theoretical and numerical approximations have been systematically bench-marked among themselves and compared with available experimental data. The overall good agreement between the first-principles computational results and available experimental data for a large amount of crystals at moderate temperatures (T ) establishes the phonon BTE as a practical and robust computational tool to design advanced technology materials with optimized thermal transport properties.
Meanwhile, concerns have been raised about the validity of the phonon BTE beyond the PG limit, where interactions among vibrational modes are significant and the weakly interacting quasi-particle approximation becomes insufficient 28 . A schematic plot of a typical temperature dependence of κ Latt in crystals is shown in Fig. 1. Within the PG approximation, the phonon BTE predicts that κ Latt of a crystal decays to zero with increasing T at the rate of 1/T or faster. However, experimental measurements 29,30 reveal that the deviation from the 1/T scaling become noticeable as T approaches the melting temperature (T melt ) of the lattice, with κ Latt eventually reaching a low constant value. The omnipresence of these minimal thermal conductivities (κ min ) 31 in all crystalline lattices suggests that as a lattice approaches its T melt , the increasingly strong anharmonic coupling among vibrational modes causes the breakdown of the PG model. Such breakdown might occur at moderate temperatures in relatively soft solids with large thermal expansion [32][33][34] , or in the high temperature phases of solids whose 0K phonon spectra contain imaginary frequencies 35 . In addition, the phonon BTE incorporates the concept of phonon group velocity, which is not properly defined in non-periodic solids such as alloys, glasses or amorphous semiconductors 36 , even at the conditions where all the vibrational modes remain quasi-harmonic 37 .
When the accuracy of the phonon BTE theory is in question, the statistical linear response transport theory 38 is often combined with equilibrium molecular dynamics (MD) simulations to predict thermal transport properties [39][40][41][42] . For example, the Green-Kubo (GK) formalism states that thermal conductivity is proportional to the time-integral of the auto-correlation function of heat flux 43,44 . Although the GK method is theoretically rigorous and valid beyond the PG approximation, its current implementations, based on the evaluations of atomic trajectories, i.e. displacements and velocities, over a long period of time, usually require much more intensive computational loads. When no reliable empirical force-field interatomic potentials exist, ab initio MD simulations are necessary to simulate the complex lattice vibration. Yet, in practice, typical ab initio MD simulations are often carried out with only relatively short simulation periods (i.e. on the order of a few pico-seconds) and using relatively small super-cell models (i.e. on the order of a couple of hundred atoms) because their computational loads scale as order N 3 , where N is the number of atoms in a supercell model. These numerical finite-size artifacts sometimes impose relatively large uncertainties in the ab initio MD simulation results. Additional approximations are often needed to extract potential energy of each atom from the ab initio total energies of the supercell models in order to evaluate the correlation function of heat currents using the ab initio MD simulation results. [45][46][47][48] .
More importantly, all the atomic trajectories in MD simulations have to be calculated numerically, even at the weak scattering limit of the PG model. This lack of analytical solutions of atomic trajectories in MD simulations hinders the development of quantitative theoretical models to interpret the simulated current-current correlation functions because it provides little insights on improving/correcting the PG model beyond the weak scattering limit. Ladd et al 49 proposed a normal mode analysis (NMA) approach to evaluate the phonon lifetimes τ based on the damped oscillator approximation (DOA). Using the extracted phonon lifetimes, they derived the so-called Peierls phonon-transport expression of κ Latt , which is understood to be only an approximate solution of the phonon BTE theory. Nevertheless, these types of NMA methods have been useful to interpret the phonon scattering in a MD simulation, and these methods have been implemented and further developed in recent years by many groups using both empirical potentials and ab initio methods [50][51][52][53] . However, both the DOA and the concept of phonon lifetime/relaxation-time should be adopted only as semi-quantitative models because the cross-correlations among different vibrational modes can not always be neglected. More robust theoretical models or concepts are needed to quantitatively interpret the NMA results of numerical MD simulations.
In this paper, we present a time-dependent statistical theory to quantitatively describe the thermal fluctuation and correlation properties of vibrational modes using a Fokker-Planck equation 54 for lattice dynamics. First, this vibration FPE theory does not treat the interactions among different vibrational modes as small perturbations. Instead, our theory includes two general sets of parameters, the drift A and the diffusion B coefficients, to explicitly characterize the mode-mode interactions. The results of this vibration FPE, expressed in terms of a time-dependent probability function of multiple-variable vibrational micro-states, provide details of the dynamic relaxation processes of lattice vibration, and are readily used by the linear response transport theory to compute κ Latt beyond the quasiharmonic PG model. Second, this vibration FPE provides detailed information on the time-correlation properties of physical quantities without requirement of long time MD simulations. The proposed vibration FPE derives the correlation functions based on the probability function governed by the drift A and diffusion B coefficients, which are defined in terms of ensemble averages at the δt → 0 limit. It is important to emphasize that no a priori forms of correlation functions are assumed in a FPE calculation of correlation functions. As a result, when implemented with first-principles methods, this vibration FPE is promising to be both accurate and efficient to predict κ Latt of novel and complex solids at wide-ranging conditions. Finally, the κ latt predicted by the vibration FPE converges to the one from the conventional phonon BTE within the PG model. Because the FPE's parameters of a lattice vibration can be evaluated with either perturbative methods or simulation methods at the PG approximation, our vibration FPE theory establishes a systematical computational methodology to analyze errors of the simple PG model and to delineate the breakdown conditions of the PG approximation.

II. STOCHASTIC DYNAMICS OF LATTICE VIBRATION
A. Fokker-Planck equation The first fundamental assumption of this proposed Fokker-Planck equation for lattice vibration is that thermal lattice dynamics is a stochastic process at the microscopic level, and the probabilitic transition dynamics from one vibration micro-state Γ to other thermally accessible micro-states can be modeled with a statistical master equation 38,54 . When a specific micro-state Γ 0 is sampled at time t = 0, the initial probability function is simply: Regardless of the dynamic details of a stochastic process, the equilibrium ensemble theory constrains that at the long time limit of t → ∞, the probability function evolves into the canonical distribution function: where k B is the Boltzmann constant, T represents temperature, E(Γ) denotes the energy of any micro-state Γ, and Z eq (T ) denotes the equilibrium canonical partition function of the lattice vibration. The evolution of this probability function P (Γ, t|Γ 0 ) provides a general and quantitative description of lattice thermal relaxation processes, from a single initially sampled micro-state Γ 0 to a set of the thermally accessible micro-states that correspond to an equilibrium distribution governed by the equilibrium statistics. Here, the ergodic condition in lattice vibration is assumed.
We further adopt the Born-von-Karman periodic boundary condition 10 to specify the vibrational micro-states with total N vibration modes, with N → ∞ for an infinitely large crystal. Using the numbers of phonons at these modes, i.e. n α with α = 1, 2, 3 · · · , N , we specify a vibrational micro-state with a set of N -dimensional state-variables Γ = {n 1 , n 2 , · · · , n N }. Through the Kramers-Moyal expansion of the master equation, the time-evolution of this probability function P (Γ, t|Γ 0 ) = P (n 1 , n 2 , · · · , n N , t|n 0 1 , n 0 2 , · · · , n 0 N ) can be expressed in the form of a FPE 38,54 : The assumption of a FPE is that the third order expansion coefficients are approximately zero. According to the Pawula theorem, all the higher order expansion coefficients are zero if the third order expansion coefficients are zero 54 . Within this theoretical framework, the drift A α (Γ) and diffusion B αβ (Γ) coefficients manifest the interactions among vibrational modes, and they are defined as:  Based on the assumption that a stochastic lattice vibration can be approximated as a random process of transition from one vibrational micro-state Γ to another micro-state Γ with a known rate of transition w Γ→Γ , Eq. 4 can be approximated as: Within the PG model, both initial (Γ) and final (Γ ) quantum vibration states can be represented by the phonon representation |n 1 , n 2 , n 3 , · · · , n N >, and ∆V denotes perturbatively small deviations in the vibration Hamiltonian from that of the ideal phonon gas. We can use Fermi's golden rule to calculate the rate of transition:

B. Thermal relaxation: fluctuation and correlation
At thermal equilibrium, the instantaneous value of a quantity X, either macroscopic or microscopic, fluctuates around its equilibrium value X eq . The dynamical process that brings the fluctuating value of X(t) back toward the X eq is commonly referred as a thermal relaxation process. A self-correlation function of X: is often used to quantify the properties of this thermal relaxation process. When X can be expressed in terms of micro-state variables X(Γ), we can define a time-dependent expectation value X(t|Γ 0 ) based on the probability function P (Γ, t|Γ 0 ) in the vibration FPE, staring with the initial probability function shown in Eq. 1: Clearly, X(t) starts at its initial value of X(Γ 0 ) = dΓδ(Γ−Γ 0 )X(Γ), and eventually relaxes back to its equilibrium value of X eq = dΓX(Γ)P eq (Γ) when P (Γ, t|Γ 0 ) → P eq (Γ) at the limit of t → ∞. Similarly, the corresponding time-dependent statistical variance, defined as By sampling the initial micro-states Γ 0 with the equilibrium probability function P eq (Γ 0 ), we can re-write the time-correlation function of X, defined in Eq. 7, as: where C XX (t = 0) = ∆ X,eq , and based on the approximation that c XX (t) ≈ e −t/τ X .
The dynamical correlation between two different quantities X and Y that fluctuate around their prospective equilibrium values (X eq and Y eq can be quantitatively formulated in terms of a cross-correlation function C XY (t): and this cross-correlation function can be re-written using the probability distribution function P (Γ, t|Γ 0 ) of Eq. 3 : where the ratio c XY ≡ C XY (t = 0)/ ∆ X,eq · ∆ Y,eq is often referred as the correlation ratio, with c XY = 0 being interpreted as that the fluctuations in X and Y are statistically uncorrelated at thermal equilibrium. It is important to emphasize that even at the condition of zero correlation ration, i.e. c XY = 0, a cross-correlation function defined in Eq. 12 in not always zero at t > 0.
Because the self-correlation function formula in Eq. 9 is a special case of the crosscorrelation function formula in Eq. 12 with X = Y , we present only the results of the time derivative of the cross-correlation function here based on Eqs. 8 and 12: where A and B are the parameters (Eq. 4 ) of the vibration FPE (Eq. 3). Using the definitions of y µ ≡ ∂Y ∂nµ · A µ and y µν ≡ ∂ 2 Y ∂nµ∂nν · B µν , we can re-write Eq. 13 in terms of the cross-correlation functions between X and y µ and those between X and y µν : Furthermore, all the higher order time derivatives of C XY (t) functions can also be derived from Eq. 14 in a recursive fashion.
Next, we summarize some key results in the case that X and Y are simply the α-th and β-th state variables n α and n β , with more details on the mathematical derivation given in Appendix A. The commonly adopted concept of phonon occupation number of a vibrational mode can be generalized as the time-dependent expectation value of the state variable n α during a thermal relaxation process, i.e. n α (t|Γ 0 ) ≡ dΓn α P (Γ, t|Γ 0 ), with n α (t|Γ 0 ) → n α,eq and ∆ α (t|Γ 0 ) ≡ n 2 α (t|Γ 0 ) − n α (t|Γ 0 ) 2 → ∆ α,eq at the t → ∞ limit. At the weak phonon scattering limit of the PG model, the thermal equilibrium values of n α,eq follow the Bose-Einstein distribution, and the corresponding statistical variances are ∆ α,eq = n α,eq (n α,eq +1).
Applying the vibration FPE ( Eq. 3) to Eq. 8, we derive the time derivatives of n α (t|Γ 0 ) and ∆ α (t|Γ 0 ) as: Furthermore, using Eqs. 11 and 12, we define the cross-correlation functions between the fluctuating phonon numbers of the α-th mode and the β-th mode (also referred to as two- with C nαn β (t = 0) = δn α · δn β eq = n α · n β eq − n α,eq · n β,eq . We can further define the normalized two-mode correlation functions as: Since X = n α and Y = n β , we have y µ = A β · δ µβ and y µν = 0. Using Eq. 14, we can show that: Multiple-mode correlation functions can be defined in a similar fashion. For example, there is only one type of three-mode correlation function among the α-th, β-th, and γ-th mode: and there are three types of four-mode correlation functions among four (α, β, µ, and ν) modes: Within the PG model, the fluctuations of phonon occupation numbers at two different modes are considered to be statistically independent at a thermal equilibrium, i.e. < n α · n β > eq = n α,eq · n β,eq for α = β. As a result, the values of the normalized time-correlation function at t = 0 are simply c αβ (t = 0) = δ αβ , where δ αβ is the Kronecker-δ symbol. Yet, the PG model does not state the value of a cross-correlation function (Eq. 16) at any other time t = 0, except that c αβ (t) → 0 as t → ∞. Multiple-mode correlation functions remain poorly understood, even within the PG model.

C. Ornstein-Uhlenbeck Processes
The FPE for a well-studied class of stochastic processes, the so-called Ornstein-Uhlenbeck (OU) processes 55 , can be solved analytically. To demonstrate the properties of these OU processes, we start with a new set of zero-mean and unit-variance stochastic variablesΓ = (x 1 , x 2 , x 3 , · · · , x N ), i.e. x λ eq = 0 and x λ 2 eq = 1. The OU processes are defined in terms of their specific form of drift and diffusion coefficients: A λ (Γ) = −γ λ x λ and B λλ (Γ) = 2γ λ δ λ,λ , with γ λ > 0. Consequently, the Fokker-Planck equation for an OU type processes can be re-written in a separable multiple-variable partial differential equation: and its solution can be expresses as: where, x λ (t) = x λ (Γ 0 ) · e −γ λ t . and ∆ λ (t) = 1 − e −2γ λ t . More details on the solution of an OU type FPE can be found in Appendix B. Here we highlight one key result of the time-correlation between any two state variables x λ and x λ of an OU type process: More interesting results on the multiple variable correlation functions, such as the threevariable correlation functions: For a lattice vibration to be classified as an OU process, its set of drift coefficients A(Γ) in the vibration FPE (Eq. 3) must satisfy the following conditions: Here D αβ are matrix elements of the normalized drift matrix D D D, n α,eq and ∆ α,eq are respectively the equilibrium average value of the phonon number at α-th mode and the corresponding statistical variance at the equilibrium with α, β = 1, 2, 3, · · · , N .
The D D D matrix, as defined in Eq. 25, is a positive definite, real, and symmetric N × N matrix with a set of N eigenvalues γ λ and corresponding normalized eigenvectors written as as u λ = (u λ,1 , u λ,2 , u λ,3 , · · · , u λ,N ) for λ = 1, 2, 3, · · · , N . We then can transform the N -dimensional phonon number state variables Γ = {n 1 , n 2 , · · · , n N } into an equivalent set of zero-mean and unit-variance state variablesΓ = (x 1 , x 2 , x 3 , · · · , x N ) using this set of eigenvectors: n α = n α,eq + (∆ α,eq ) 1/2 N λ=1 x λ u λ,α , The linear transformation in Eq. 26 also shows that the diffusion B αβ (Γ) coefficients for an OU type lattice vibration are related to its drift coefficients A α (Γ) through the D D D matrix: In the rest of the paper, the D D D matrix is referred as the normalized drift/diffusion matrix.
We refer to χ αβ (γ) as the Laplace spectral function of c αβ (t). At the N → ∞ limit, a Laplace spectral function χ αβ (γ) converges to a continuous function defined in the spectral regime of [0, γ max ]. The k-th moment of a χ αβ (γ) function, defined as µ αβ (k) ≡ ∞ 0 dγχ αβ (γ) · γ k , is given as: The results in Eqs. 28 and 29 clearly demonstrate that in general the normalized mode self-correlaltion functions of lattice vibration do not decay as an exponetial function of time, and the time-integral of the cross-correlaltion functions are not zero for two different modes.

A. Green-Kubo Theory
The fluctuation-dissipation theorem provides a general statistical theory to connect the equilibrium fluctuation processes of a macroscopic quantity e.g. the total heat current vector J = (J x , J y , J z ) in a solid and the related irreversible transport processes, such as heat conduction at non-equilibrium conditions. Within the statistical linear response transport theory, the thermal conductivity tensor κ IJ , with I, J = x, y, z labeling the Cartesian axes, is expressed in the Green-Kubo formula in terms of the time integral of the current-current correlation functions 43,44 : where Ω and N cell are respectively volume of the unit-cell and total number of cells in a super-cell model with the Born-von Karman periodic boundary.
At the atomistic level, the heat current J is a function of atomic forces, displacements and momenta, and various approximations have been proposed and discussed 56 . Assuming the heat current vector is also a function of phonon numbers of modes, i.e. J = J({Γ}) = J(n 1 , n 2 , n 3 , · · · , n N ), we can use Eq. 9 to evaluate the current-current correlation functions.
Under the condition of small thermal fluctuation, the Cartesian components of the heat current vector can be simplified as: The seminal Peierls formula of the heat current of a phonon gas, J = α ∆n α ω α v α , is an approximation of this class, with Λ Iα = ω α v αI . When the higher order terms(also referred as the non-harmonic terms) in the J formula are included as the corrections to the linear terms formulated in Eq. 36, we can re-write the J I as J I = i Λ Iα ∆n α + δJ I . Consequently, the current-current correlation functions can be expressed as: Wherever the non-harmonic δ J terms in the vibrational heat current in a lattice are not Based on the GK formula, we now express κ Latt in the form of: As shown in Eq. 28 of Sec. II C, when a lattice vibration can be approximated as an Ornstein-Uhlenbeck process, the lattice thermal conductivity is simply:

B. Phonon Boltzmann Transport Equation
As a kinetic transport theory, the phonon BTE theory is valid only within the PG approximation, i.e. at a thermal equilibrium, each mode oscillates at a harmonic frequency ω and the ensemble averaged number of phonons at this mode follows the Bose-Einstein distribution n eq (ω) = n eq = 1 e ( ω/k B T )−1 and ∆ eq = n 2 eq − n 2 eq = n eq · (n eq + 1). In addition, the phonon BTE theory applies only to a crystalline solid, where each vibrational mode of this translation-invariant periodic lattice corresponds to a reciprocal-space k vector and a group velocity v = k (ω).
When a constant temperature gradient r T is imposed on the periodic lattice, the ensemble averaged phonon numbers, n α for α = 1, 2, 3, · · · , N , are no longer able to relax back to their original equilibrium values n eq,α as a result of thermal diffusion. Instead, each n α approaches a space-dependent value when a steady-state is reached: where the diffusion term at the r T → 0 limit is approximated as: A common approximation for the scattering terms in the phonon BTE (Eq. 41) is the so-called linearized approximation: n α,eq · (n α,eq + 1) n β,eq · (n β,eq + 1) · L αβ · (n β − n β,eq ), (43) where L L L is referred as the linear phonon scattering matrix.
By using the results of Eqs. 42 and 43 and the definition of φ α ≡ n α − n α,eq n α,eq · (n α,eq + 1) , the steady-state phonon Boltzmann equation (Eq. 41) can be re-written as a set of linear equations for φ α with α = 1, 2, 3, · · · , N : Similar to what we have derived in Sec. II C, we can solve the set of linear equations using the eigenvectors and the eigenvalues of the matrix L L L: where γ λ and u λ are the λ-th eigenvalue and eigenvector of the matrix L L L, and (L L L) −1 represents the inverse matrix of L L L.
Based on the Peierls formula for the heat current of a phonon gas, the lattice thermal conductivity predicted by the linearized phonon BTE theory can be expressed as: n α,eq (n α,eq + 1) n β,eq (n β,eq + 1) where c = k B · ( ω k B T ) 2 · n α,eq · (n α,eq + 1) is the single mode heat capacity. To compare κ Latt predicted by the phonon BTE (Eq. 46) and the one by the OU type vibration FPE (Eq. 40), we first note that in the limit of weak phonon scattering of the PG model, the variance of the phonon number fluctuation of a mode ∆ α,eq has already been shown to converge to the value of n α,eq · (n α,eq + 1), and the Peierls formula of heat current is

C. Discussions
A comparison chart is shown in Table I  to the scattering L L L matrix in a phonon BTE, we have proved for the first time that the κ Latt derived from the linear response transport theory converges to that from the kinetic transport theory within the PG approximation.
When the interactions among vibrational modes are perturbatively small, the normalized drift/diffusion D D D matrix can be derived by using quantum perturbation theories for lattice vibration at low temperature. As temperature elevates to the semi-classical and classical regime, we can implement numerical algorithms to directly compute normalized drift/diffusion coefficients with first-principles MD simulations. As these coefficients are de-

Molecular Dynamics Simulation
Fokker-Plank Equation
V ∼ many-body inter-atomic potential Quantum, Semi-Classical, and Classical Normal mode analysis: • Long simulation time is necessary • Self-correlation functions c αα (t) often assumed to be e −t/τα • Cross correlations c αβ (t) often assumed to be zero • Three/four-mode correlation functions rarely discussed

OU Approximation of drift/diffusuion coefficients ⇐⇒ Phonon Gas Model
Analytical formulas of multiple-mode time-correlation function: • Two-mode: c αβ = δn α (0) · δn β (t) eq / δn 2 α eq · δn 2 • Four-mode: δn α (0)δn β (0)δn µ (0)δn ν (t) eq = [δ αµ c βν (t) + δ βµ c αν (t) + δ αβ c µν (t)] δn 2 α eq δn 2 β eq δn 2 µ eq δn 2 ν eq fined in the short time limit, high-performance parallel computer platforms can be utilized to distribute the computational loads of such simulations in parallel. When the temperature dependence of the drift/diffusion coefficients are extracted and tested with the OU approximation, we are able to not only quantitatively determine the temperature condition in which the PG model breaks down, but also identify the individual vibrational modes that lead to the breakdown. To study the mechanisms of lattice heat conduction beyond the PG model, it is critical to establish a quantatitive criterion that delineates the breakdown conditions. The theoetical analysis presented in this paper indicates that the OU condition of stochastic lattice dynamics (Eqs. 25   Xu for discussions.
Using the vibration FPE shown in Eq. 3, we then prove that the first-order t-derivatives of these three quantities in Eqs. A1 to A2 have the following forms: As a result, Eq. A4 is now simplified as: We now get: Similarly, we can show that Consequently, we also have: (A8) Appendix B: Analytical solutions of the FPE for an OU process In this appendix we verify the analytical solutions of an OU type FPE (Eq. 22) shown in Sec. II C. For a probability function of one stochastic variable x with zero-mean and unit-variance, the corresponding OU type FPE can is given as: With d dt x = −γx, d dt x 2 = 2γ(1 − x 2 ), and the initial values x and x 2 being x 0 and (x 0 ) 2 respectively, we have: We skip the details of derivation and only show that the analytical solution of Eq. B1 is given as: Using the analytical solution in Eq. B3, we can show that ∂P (x,t) ∂t on the left hand side of Eq. B1 is: ]P (x, t) = γ[1 − ∆ −1 (t) · (1 − xx(t) + x 2 ) + ∆ −2 (t) · (x − x(t)) 2 ]P (x, t).
Plugging these results in the separable multiple-variable formula, we now can verify that the analytical solution of Eq. 22 is indeed the probability function shown in Eq. 23.