Thermal equilibration in a one-dimensional damped harmonic crystal

TL;DR

This paper investigates how external damping affects thermal equilibration in a one-dimensional harmonic crystal, revealing a two-stage energy decay process with specific asymptotic behaviors for kinetic and potential energies.

Contribution

It provides an analytical description of the damping influence on energy oscillations and decay in a harmonic crystal, supported by numerical verification.

Findings

Energy oscillations decay exponentially in the underdamped case

At large times, potential energy dominates with a $t^{-3/2}$ decay

Kinetic energy decays as $t^{-5/2}$ at large times

Abstract

The features for the unsteady process of thermal equilibration ("the fast motions") in a one-dimensional harmonic crystal lying in a viscous environment (e.g., a gas) are under investigation. It is assumed that initially the displacements of all the particles are zero and the particle velocities are random quantities with zero mean and a constant variance, thus, the system is far away from the thermal equilibrium. It is known that in the framework of the corresponding conservative problem the kinetic and potential energies oscillate and approach the equilibrium value that equals a half of the initial value of the kinetic energy. We show that the presence of the external damping qualitatively changes the features of this process. The unsteady process generally has two stages. At the first stage oscillations of kinetic and potential energies with decreasing amplitude, subjected to exponential decay, can be observed (this stage exists only in the underdamped case). At the second stage (which always exists), the oscillations vanish, and the energies are subjected to a power decay. The large-time asymptotics for the energy is proportional to $t^{-3/2}$ in the case of the potential energy and to $t^{-5/2}$ in the case the kinetic energy. Hence, at large values of time the total energy of the crystal is mostly the potential energy. The obtained analytic results are verified by independent numerical calculations.

Equations

$${\partial }_t{v}_i=F_i,{\partial }_t{u}_i=v_i,$$

$$\displaystyle F_{i}={\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}_{i}u_{i}-\eta v_{i},$$

$${\omega }_0\stackrel{\text{mbox}def}{=}\sqrt{C/m}.$$

$${\mathcal{L}}_i{u}_i=u_{i+1}-2u_i+u_{i-1}.$$

$$u_i(0)=0,v_i(0)={\rho }_i,$$

$$⟨{\rho }_i(0)⟩=0,⟨{\rho }_i(0){\rho }_j(0)⟩={\kappa }_i(0){\delta }_{ij},$$

$$⟨u_i⟩=0,⟨v_i⟩=0.$$

$${\xi }_{p,q}\stackrel{\text{mbox}def}{=}⟨u_p{u}_q⟩,{\nu }_{p,q}\stackrel{\text{mbox}def}{=}⟨u_p{v}_q⟩,{\kappa }_{p,q}\stackrel{\text{mbox}def}{=}⟨v_p{v}_q⟩.$$

$${\mathcal{K}}_{p,q}=\frac{m{\kappa }_{p,q}}{2}$$

$$T_p\stackrel{\text{mbox}def}{=}2k_B^{-1}{\mathcal{K}}_{p,p}=m{k}_B^{-1}{\kappa }_{p,p}$$

$${\partial }_t\xi =\nu +{\nu }^{\top },$$

$$\displaystyle\partial_{t}\nu+\eta\nu={\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}_{q}\xi+\kappa,$$

$$\displaystyle\partial_{t}\kappa+2\eta\kappa={\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}_{p}\nu+{\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}_{q}\nu^{\top},$$

$${\xi }_{pq}(0)=0,{\nu }_{pq}(0)=0,{\kappa }_{pq}(0)={\kappa }_p(0){\delta }_{pq}.$$

$$\partial_{t}\kappa+2\eta\kappa={\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}_{p}\nu+{\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}_{q}\nu^{\top}+\beta,$$

$$\beta ={\kappa }_p(0){\delta }_{pq}\delta (t).$$

$$\xi\big{|}_{t<0}\equiv 0,\qquad\nu\big{|}_{t<0}\equiv 0,\qquad\kappa\big{|}_{t<0}\equiv 0.$$

$$2{\mathcal{L}}^{\mathrm{S}}\stackrel{\text{mbox}def}{=}{\mathcal{L}}_p+{\mathcal{L}}_q,2{\mathcal{L}}^{\mathrm{A}}\stackrel{\text{mbox}def}{=}{\mathcal{L}}_p-{\mathcal{L}}_q,$$

$$2{\nu }^{\mathrm{S}}\stackrel{\text{mbox}def}{=}\nu +{\nu }^{\top },2{\nu }^{\mathrm{A}}\stackrel{\text{mbox}def}{=}\nu -{\nu }^{\top }.$$

$$\displaystyle\partial_{t}\xi=2\nu^{\mathrm{S}},\qquad(\partial_{t}+2\eta)\kappa=2{\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}^{\mathrm{S}}\nu^{\mathrm{S}}+2{\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}^{\mathrm{A}}\nu^{\mathrm{A}}+\beta,$$

$$\displaystyle(\partial_{t}+\eta)\nu^{\mathrm{A}}=-{\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}^{\mathrm{A}}\xi,\qquad(\partial_{t}+\eta)\nu^{\mathrm{S}}={\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}^{\mathrm{S}}\xi+\kappa.$$

$$\left((\partial_{t}+\eta)^{2}(\partial_{t}^{2}+2\eta\partial_{t}-4{\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}^{\mathrm{S}})+4({\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}^{\mathrm{A}})^{2}\right)\kappa\\ =(\partial_{t}+\eta)(\partial_{t}^{2}+\eta\partial_{t}-2{\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}^{\mathrm{S}})\beta,$$

$$\big{(}(\partial_{t}+\eta)^{2}(\partial_{t}^{2}+2\eta\partial_{t}-4{\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}^{\mathrm{S}})+4({\color[rgb]{0,0,0}\omega_{0}^{2}{\mathcal{L}}}^{\mathrm{A}})^{2}\big{)}\xi=2(\partial_{t}+\eta)\beta.$$

$$k\stackrel{\text{mbox}def}{=}p+q$$

$$n\stackrel{\text{mbox}def}{=}q-p$$

$$q=\frac{k}{2}+\frac{n}{2},p=\frac{k}{2}-\frac{n}{2}.$$

$${\kappa }_p(0)={\bar{\kappa }}_0,$$

$$T_{p}\big{|}_{t=0}=\bar{T}_{0}\equiv 2k_{B}^{-1}\bar{\mathscr{E}}_{0},$$

$${\bar{\mathcal{E}}}_0\equiv \frac{m{\bar{\kappa }}_0}{2}$$

$${\kappa }_{pq}={\hat{\kappa }}_n,{\xi }_{pq}={\hat{\xi }}_n.$$

$$4{\Pi }_{p,p}=m{\omega }_0^2\left(⟨(u_{p-1}-u_p{)}^2⟩+⟨(u_{p+1}-u_p{)}^2⟩\right).$$

$$4\Pi_{p,p}\equiv 2\hat{\Pi}_{0}\\ =m\omega_{0}^{2}\big{(}2\xi_{p,p}+\xi_{p-1,p-1}-2\xi_{p,p-1}+\xi_{p+1,p+1}-2\xi_{p,p+1}\big{)}\\ =-2m\omega_{0}^{2}{\mathcal{L}}_{n}\hat{\xi}_{n}\big{|}_{n=0}.$$

Full text

I Introduction

In this paper, we study the influence of an external viscous environment (e.g., a gas) on the unsteady process of thermal equilibration in an infinite one-dimensional harmonic crystal with nearest-neighbor interactions. The model of a damped harmonic crystal was used in our recent papers Gavrilov et al. (2019); Gavrilov and Krivtsov (2019), where we discuss the ballistic heat propagation in such a structure (“the slow motions”). As opposed to Gavrilov et al. (2019); Gavrilov and Krivtsov (2019) now we consider the fast motions, i.e. the fast vanishing oscillations of the kinetic and potential energy. This oscillations are well known for those who deal with molecular dynamics simulation (see, e.g., Allen and Tildesley (1987), Fig. 5.11).

We assume that initially the displacements of the crystal particles are zero and the particle velocities are random quantities with zero mean and a constant variance. The kinetic energy per particle (as well as the corresponding kinetic temperature) is distributed spatially uniform, whereas the potential energy is zero. Thus, the thermodynamic system is far away from thermal equilibrium. It is well known that in the framework of the corresponding conservative problem (in the absence of external damping) the kinetic and potential energies oscillate and approach the equilibrium value that equals a half of the initial value of the total energy (equipartition of kinetic and potential energy). For a harmonic crystal the process of thermal equilibration was first time investigated by Klein & Prigogine in Klein and Prigogine (1953). Thermal equilibration in harmonic crystals in the conservative case was considered in many studies, e.g., Linn and Robertson (1984); Hemmer (1959); Huerta et al. (1971); Robertson and Huerta (1969, 1970); Kannan et al. (2012); Lepri et al. (2008); Rieder et al. (1967); Lepri et al. (2010); Krivtsov (2014); Babenkov et al. (2016); Kuzkin and Krivtsov (2017a); Guzev (2018); Kuzkin and Krivtsov (2017b, c). The more complete bibliography can be found in recent paper by Kuzkin Kuzkin (2019), where the analytic solution in the integral form is obtained for an infinite harmonic crystal with an arbitrary Bravais lattice and a polyatomic cell with an arbitrary structure. The time evolution the kinetic temperature during the thermal equilibration in a harmonic crystal was considered in Refs. Klein and Prigogine (1953); Hemmer (1959); Krivtsov (2014); Babenkov et al. (2016); Kuzkin and Krivtsov (2017a); Guzev (2018); Kuzkin and Krivtsov (2017b, c); Kuzkin (2019), the entropy was under consideration in Refs. Huerta et al. (1971); Robertson and Huerta (1969, 1970); Sokolov et al. (2019). Thermal equilibration for a system of quantum oscillators is considered in Usha Devi and Rajagopal (2009).

In the paper, we show that the process of thermal equilibration in the presence of a viscous external environment is more complicated than in the conservative case and has two stages in the underdamped case 111In the underdamped case the specific viscosity of the external environment is small enough such that restriction (57) is fulfilled, which is generally assumed in the paper. In the presence of damping, the limiting values for the kinetic and potential energies, clearly, are zero. At the first stage, the qualitative description of the process is as follows: the kinetic and potential energies oscillate approaching an exponentially-decaying curvilinear asymptote. Unexpectedly, for any positive value of the specific viscosity for environment, there is the second stage, which can observed in the underdamped case only for very large values of time. The kinetic and the potential energies at the second stage are subjected to a power decay. Another one unexpected result is as follows: the principal term of the large-time asymptotics is proportional to $t^{-3/2}$ in the case of potential energy and to $t^{-5/2}$ in the case of the kinetic energy. Hence, at very large times the total energy of the harmonic crystal is mostly the potential energy.

The paper is organized as follows. In Section II, we consider the formulation of the problem. In Section II.1, some general notation is introduced. In Section II.2, we state the basic equations for the crystal particles in the form of a system of ordinary differential equations with random initial conditions. In Section II.3, we introduce and deal with infinite set of covariance variables. These are the mutual covariances of the particle velocities and the displacements for all pairs of particles. We obtain two infinite systems of differential-difference equations involving only the covariances for the particle velocities, and only the covariances for the displacements, respectively. The similar approach was used in previous papers Krivtsov (2014, 2015, 2019); Sokolov et al. (2017); Krivtsov et al. (2018); Babenkov et al. (2016); Kuzkin and Krivtsov (2017c, a, b); Murachev et al. (2018); Sokolov et al. (2019); Gavrilov et al. (2019); Gavrilov and Krivtsov (2019). In Section III we simplify the obtained equations using the assumptions of uniformity for the variance of the initial values of the particle velocities. Finally, we obtain four infinite systems of ordinary differential equations for the energetic quantities, which we call the generalized kinetic energy, the generalized potential energy, the generalized total energy, the generalized Lagrangian. The kinetic energy, the potential energy, the total energy, and the Lagrangian, are particular cases of those quantities. It is sufficient to solve the equations for any two of these four energetic quantities to calculate all of them. We choose the generalized Lagrangian and the generalized potential energy as basic variables, since the corresponding equations have a simpler structure. In Section IV we use discrete-time Fourier transform to get the analytical solutions for the Lagrangian and the potential energy in the integral form. In Section V all the energetic quantities are evaluated for the large values of time. To do this we use the method of stationary phase and the Laplace method Fedoryuk (1977). The corresponding calculations are given in Appendices A–D. In Section VI, we present the results of the numerical solution of the initial value problem for the system of ordinary differential equations with random initial conditions and compare them with the obtained analytical solution in the integral and the asymptotic forms. In Section VII we discuss the large-time behavior of energetic quantities in the underdamped case. Finally, in Section VIII, we discuss the basic results of the paper.

II Mathematical formulation

II.1 Notation

In the paper, we use the following general notation:

$t$

the time;

$H(\cdot)$

the Heaviside function;

$\delta(\cdot)$

the Dirac delta function;

$\langle\cdot\rangle$

the expected value for a random quantity;

$\delta_{pq}$

the Kronecker delta ($\delta_{pq}=1$ if $p=q$, and $\delta_{pq}=0$ otherwise);

$\delta_{n}$

$\delta_{n}\stackrel{{\scriptstyle\mbox{\scriptsize def}}}{{=}}\delta_{n0}$;

$J_{0}(\cdot)$

the Bessel function of the first kind of zero order Abramowitz and Stegun (1972);

$\Gamma(\cdot)$

the Euler integral of the second kind (the Gamma function) Abramowitz and Stegun (1972);

$C^{\infty}$

the set of all infinitely differentiable functions;

$\mathbb{Z}$

the set of all integers.

II.2 Dynamic equations for a crystal and random initial conditions

Consider the following system of ordinary differential equations:

[TABLE]

where

[TABLE]

Here $i$ is an arbitrary integer which describes the position of a particle in the chain; $u_{i}(t)$ and $v_{i}(t)$ are the displacement and the particle velocity, respectively; $F_{i}$ is the specific force on the particle; $\eta$ is the specific viscosity for the environment; $C$ is the bond stiffness; $m$ is the mass of a particle; $\partial_{t}$ is the operator of differentiation with respect to time; ${\mathcal{L}}_{i}$ is the linear finite difference operator:

[TABLE]

System of ODE (1) describes the motions of one-dimensional harmonic crystal (an ordered chain of identical interacting material particles, see Fig. 1).

The initial conditions are as follows: for all $i$,

[TABLE]

where the normal random variables $\rho_{i}$ are such that

[TABLE]

where $\kappa_{i}(0)$ is a given function of $i\in\mathbb{Z}$. Later, in Section III, it will be assumed that $\kappa_{i}(0)$ do not depend on $i$. It is useful to proceed with the derivation of basic equations in Section II.3 not taking into account this supposition. Note that in the latter more general case, the boundary conditions at the infinity may be needed. These boundary conditions should guarantee that there are no sources at the infinity, or give a mathematical description for such a source. In the case under consideration in the paper, these special boundary conditions are not necessary, we will use the requirement of spatial uniformity for all physical quantities instead. At the same time, in numerical calculations (see Section VI), where we deal with a model for a finite harmonic crystal, we use periodic boundary conditions (106) to provide the spatial uniformity.

II.3 The dynamics of covariances

According to (2), $F_{i}$ are linear functions of $u_{i},\ v_{i}$. Taking this fact into account together with Eqs. (5), (6), we see that for all $t$

[TABLE]

Following Krivtsov (2016), consider the infinite sets of covariance variables

[TABLE]

Thus, the variables $\xi_{p,q},\ \nu_{p,q},\ \kappa_{p,q}$ are defined for any pair of crystal particles.

We call the quantity

[TABLE]

the generalized kinetic energy. It is clear that quantities $m\kappa_{p,p}$ equal the estimated value for doubled kinetic energy $2\mathscr{K}_{p,p}$ for particle with number $p$. Accordingly, we identify the following quantities

[TABLE]

as the kinetic temperature. Here $k_{B}$ is the Boltzmann constant.

For simplicity, in what follows, we drop the subscripts $p$ and $q$, i.e., $\xi\stackrel{{\scriptstyle\mbox{\scriptsize def}}}{{=}}\xi_{p,q}$ etc. By definition, we also put $\xi^{\top}\stackrel{{\scriptstyle\mbox{\scriptsize def}}}{{=}}\xi_{q,p}$ etc. Now we differentiate variables (8) with respect to time taking into account equations of motion (1). This yields the following closed system of differential equations for covariances:

[TABLE]

where ${\mathcal{L}}_{p}$ and ${\mathcal{L}}_{q}$ are the linear difference operators defined by Eq. (4) that act on $\xi_{p,q},\ \nu_{p,q},\ \kappa_{p,q},\ \beta_{p,q}$ with respect to the first index subscript $p$ and the second one $q$, respectively. The initial conditions that correspond to Eqs. (5), (6) are

[TABLE]

Taking into account these initial conditions, it is useful to rewrite Eq. (13) in the following form

[TABLE]

where singular term $\beta$ is

[TABLE]

Equations (11), (12), (15) should be supplemented with initial conditions in the following form, which is conventional for distributions (or generalized functions) Vladimirov (1971):

[TABLE]

Now we introduce the symmetric and antisymmetric difference operators

[TABLE]

and the symmetric and antisymmetric parts of the variable $\nu$:

[TABLE]

Note that $\xi$ and $\kappa$ are symmetric variables. Now Eqs. (11), (11), (15) can be rewritten as follows:

[TABLE]

This system of equations can be reduced (see Gavrilov et al. (2019)) to one equation of the fourth order in time for covariances of the particle velocities $\kappa$

[TABLE]

or, alternatively, to one equation of the fourth order in time for covariances of the displacements $\xi$:

[TABLE]

III The case of a uniform initial kinetic temperature distribution

Following Krivtsov (2015, 2016), we introduce the discrete spatial variable

[TABLE]

and the discrete correlational variable

[TABLE]

instead of discrete variables $p$ and $q$. We have

[TABLE]

In what follows, we consider the case when the initial values $\kappa_{p}(0)\delta_{pq}$ of the covariance variables $\kappa_{pq}$ do not depend on the spatial variable $k$, and depend only on $n$, i.e.

[TABLE]

where $\bar{\kappa}_{0}$ is a given constant. From the physical point of view this means that we have a uniform distribution of the initial value for the kinetic temperature

[TABLE]

where

[TABLE]

is the initial value for both the total and the kinetic energy. In the case of a uniform initial conditions it is natural to assume that for $t>0$ we also have

[TABLE]

In the case (27) of the uniform initial kinetic temperature distribution, in order to measure the estimated value of the (doubled) potential energy, we introduce the following quantity

[TABLE]

Provided that (30) are true, one has

[TABLE]

We call quantities $\Pi_{p,q}\equiv\hat{\Pi}_{n}=-\frac{1}{2}m\omega_{0}^{2}{\mathcal{L}}_{n}\hat{\xi}_{n}$ the generalized potential energy. The following identities are true for any quantity $\zeta_{p,q}$ such that $\zeta_{p,q}=\zeta_{q,p}=\hat{\zeta}_{n}:$

[TABLE]

Now, taking into account Eqs. (17), (16), Eqs. (22), (23), can be rewritten as

[TABLE]

respectively.

It is useful to consider also the following quantities

[TABLE]

and

[TABLE]

We call these quantities the generalized Lagrangian and the generalized total energy, respectively. The corresponding equations for these quantities can be obtained by means of applying of the operator ${\mathcal{L}}_{n}$ to (40) and calculating the sum or difference of both parts of Eqs. (39), (40). Taking into account initial condition in the form of (17), this yields

[TABLE]

respectively. The corresponding initial conditions are

[TABLE]

To calculate every energetic quantity from the set $\hat{\mathscr{K}}_{0}$, $\hat{\Pi}_{0}$, $\hat{\Lambda}_{0}$, $\hat{\mathscr{E}}_{0}$ it is enough to solve any two equations from the set (39), (40), (43), (46). In what follows, we deal with Eqs. (43), (40) which have a simpler structure. The generalized kinetic energy and the generalized total energy in this case can be calculated as follows:

[TABLE]

IV Solution of the equations for energetic quantities

IV.1 The Lagrangian $\hat{\Lambda}_{0}$

We apply the discrete-time Fourier transform Proakis and Manolakis (1996); Slepyan and Yakovlev (1980) with respect to the variable $n$ to Eq. (43). This yields

[TABLE]

where

[TABLE]

Here and in what follows, ${\mathfrak{q}}$ is the wavenumber, $\mathrm{i}$ is the imaginary unit. In order to obtain Eqs. (50), (51) we used the shift property Proakis and Manolakis (1996) of the discrete-time Fourier transform:

[TABLE]

Equation (50) together with initial conditions in the form of Eq. (47) is equivalent Vladimirov (1971) to initial value problem for the corresponding homogeneous equation with the following classical initial conditions:

[TABLE]

The corresponding solution is

[TABLE]

where

[TABLE]

In what follows, we generally assume that the underdamped case

[TABLE]

is under consideration. The critically damped and the overdamped cases are briefly discussed in Section V.6 (see also Figures 4, 5). Now we apply the inverse transform

[TABLE]

and get the solution in the integral form for $\hat{\Lambda}_{n}$. The Lagrangian $\hat{\Lambda}_{0}$ equals

[TABLE]

IV.2 The potential energy $\hat{\Pi}_{0}$

Applying the discrete-time Fourier transform to Eq. (40) yields

[TABLE]

Equation (60) together with initial conditions in the form of Eq. (47) is equivalent Vladimirov (1971) to the initial value problem for the corresponding homogeneous equation with the following classical initial conditions:

[TABLE]

The corresponding solution can be written as follows:

[TABLE]

Now we apply the inverse transform

[TABLE]

and get the solution in the integral form for $\hat{\Pi}_{n}$. The potential energy $\hat{\Pi}_{0}$ equals

[TABLE]

IV.3 The conservative case $\eta=0$

In the conservative case integral (58) can be calculated in the closed form:

[TABLE]

This yields Prudnikov et al. (1986)

[TABLE]

For large time $t$ the asymptotics of the Bessel function $J_{2n}$ is well known Abramowitz and Stegun (1972), thus (69) can be written in the asymptotic form

[TABLE]

In the conservative case it is useful to take $\hat{\mathscr{E}}_{n}$ as the second energetic variable (instead of $\hat{\Pi}_{n}$) since Eq. (46) simplifies to

[TABLE]

The corresponding solution is

[TABLE]

i.e. the total energy $\hat{\mathscr{E}}_{0}$ conserves, keeping the initial value for the total (or kinetic) energy.

According to Eqs. (41), (42), the generalized kinetic $\hat{\mathscr{K}}_{n}$ and potential $\hat{\Pi}_{n}$ energies equal

[TABLE]

These results are in agreement with ones previously obtained in Klein and Prigogine (1953); Krivtsov (2014).

V Asymptotics for the energetic quantities as $t\to\infty$

Unlike the conservative case the inverse Fourier transforms of quantities (55) and (62) cannot be evaluated in closed forms. Instead of this for large time we can proceed with asymptotic estimation of the corresponding integrals.

In this section we use the following notation

[TABLE]

V.1 The Lagrangian $\hat{\Lambda}_{0}$

Due to (55) the integral in the right-hand side of (58) can be represented as follows:

[TABLE]

It is easy to see that integral $\hat{\Lambda}_{n}^{(1)}$ defined by Eqs. (80), (55) is a sum of two Laplace type integrals and, therefore, it can be estimated by the Laplace method Fedoryuk (1977). On the other hand, integral $\hat{\Lambda}_{n}^{(2)}$ (also defined by Eqs. (80), (55)) is a Fourier type integral and therefore it can be estimated by the method of stationary phase Fedoryuk (1977).

Calculation of the asymptotics for quantities $\hat{\Lambda}_{n}^{(2)}$, $\hat{\Lambda}_{n}^{(1)}$ is presented in Appendices A–B, respectively. In the most interesting particular case $n=0$, the result is 222In what follows, we will see that to calculate the principal term of the asymptotics for $\hat{\mathscr{K}}_{n}^{(1)}$ we need to calculate two first non-zero terms of expansion for $\hat{\Lambda}_{n}^{(1)}$ (see (101) and text below the formula).:

[TABLE]

V.2 The potential energy $\hat{\Pi}_{0}$

Calculation of the asymptotics for generalized potential energy $\hat{\Pi}_{n}$ is presented in Appendices C–D, respectively. We have

[TABLE]

where $\hat{\Pi}_{n}^{(1)}$ is a Laplace type integral, and $\hat{\Pi}_{n}^{(2)}$ is a Fourier type integral. For $n=0$ one gets

[TABLE]

Here symbol $\operatorname{PV}$ means the Cauchy principal value for the corresponding improper integral.

V.3 The total energy $\hat{\mathscr{E}}_{0}$

Calculating the asymptotics by formulas Calculating the asymptotics by Eqs. (41), (42), (80), (87) results in:

[TABLE]

For the non-oscillating term $\hat{\mathscr{E}}_{0}^{(1)}$ due to Eqs. (83), (91) one has

[TABLE]

For the oscillating term $\hat{\mathscr{E}}_{0}^{(2)}$ due to (86), (92) one gets

[TABLE]

The latter term becomes zero as $\eta\to+0$.

V.4 The kinetic energy $\hat{\mathscr{K}}_{0}$

Calculating the asymptotics by Eqs. (41), (80), (87) results in:

[TABLE]

For the non-oscillating term $\hat{\mathscr{K}}_{n}^{(1)}$ due to Eqs. (83), (91) one has

[TABLE]

for any integer $n$. Note that the principal term of expansion for $\hat{\mathscr{K}}_{n}^{(1)}$ is of order $t^{-5/2}$ and does not depend on $n$, whereas expansions for $\hat{\Pi}_{n}^{(1)}$, $\hat{\Lambda}_{n}^{(1)}$, $\hat{\mathscr{E}}_{n}^{(1)}$ have principal terms of order $t^{-3/2}$.

For the oscillating term $\hat{\mathscr{K}}_{0}^{(2)}$ due to (86), (92) one gets

[TABLE]

V.5 The conservative case ($\eta=0$)

In the particular case $\eta=0$ the terms of asymptotic expansions for energetic quantities with superscript “(1)” equal zero since the integration is carried out over the interval of zero length. The corresponding asymptotic formulas for these terms are not valid. At the same time, according to Eqs. (161), (162)

[TABLE]

This yields the same results as ones obtained in Section IV.3 by a different approach.

V.6 The critically damped and the overdamped cases

In the critically damped ($\eta=4\omega_{0}$) and the overdamped cases ($\eta>4\omega_{0}$) the terms of asymptotic expansions for energetic quantities with superscript “(2)” equal zero since the integration is carried out over the interval of zero length. All other formulas remain valid.

VI Numerics

In this section, we present the results of the numerical solution of the system of ordinary differential equations (1) with random initial conditions (5), (6). It is useful to rewrite Eqs. (1) in the dimensionless form

[TABLE]

where

[TABLE]

Here $a$ is a constant with dimension of length, e.g., the lattice constant (the distance between neighboring particles) Gavrilov et al. (2019); Gavrilov and Krivtsov (2019). We consider the chain of $2N+1$ particles and the periodic boundary conditions

[TABLE]

The initial conditions that correspond to (5) are

[TABLE]

where $\rho_{i}$ are generated normal random numbers that satisfy (6), where, without loss of generality, we can take $\kappa_{i}(0)=1$. We use SciPy software Jones et al. : the numerical solutions of system of ODE (104) are found using the standard Python routine scipy.integrate.odeint. We perform a series of $r=1\dots R$ realizations of these calculations (with various independent $\rho_{(r)i}$) and get the corresponding displacements and particle velocities as functions of discrete time $t^{j}$: $\tilde{u}_{(r)i}(t^{j})$ and $\tilde{v}_{(r)i}(t^{j})$, respectively, where $t=\tilde{t}^{j}\stackrel{{\scriptstyle\mbox{\scriptsize def}}}{{=}}j\Delta\tilde{t}$.

According to Eqs. (8), (9), (32), and (29) the ratios $\tilde{\mathscr{K}}_{0}{\,}_{i}$ and $\tilde{\Pi}_{0}{\,}_{i}$ of the dimensionless kinetic and potential energies to the initial value of the total energy can be calculated as the following averages:

[TABLE]

respectively. The ratios $\tilde{\Lambda}_{0\,i}$ and $\tilde{\mathscr{E}}_{0}{\,}_{i}$ of the dimensionless Lagrangian and the total energy to the initial value of the dimensionless total energy are

[TABLE]

All calculations were performed for the following values of the problem parameters: $N=20$, $R=10000$. We verify that the numerical results for $\tilde{\mathscr{K}}_{0}{\,}_{i}$ and $\tilde{\Pi}_{0}{\,}_{i}$ actually do not depend on $i$, and use quantities $\tilde{\mathscr{K}}_{0}\equiv\tilde{\mathscr{K}}_{0}{\,}_{0}$ and $\tilde{\Pi}_{0}\equiv\tilde{\Pi}_{0}{\,}_{0}$ to compare the numerical and analytical results.

Numerical results for $\tilde{\mathscr{E}}_{0},$ $\tilde{\mathscr{K}}_{0},$ $\tilde{\Pi}_{0},$ $\tilde{\Lambda}_{0}$ can be compared with the analytical solutions in the integral form given by Eqs. (49), (48), (67), (59), respectively, and corresponding asymptotics (see formulas in Section V). The analytical solutions in the integral form are calculated using the standard Python routine scipy.integrate.quad. A comparison in the case $\tilde{\eta}=0.5$ (the underdamped case) is presented in Figures 2–3. Figure 3 corresponds to the case of large values of time, when exponentially-decaying terms of the asymptotics almost vanish. One can see that the numerical and analytical solutions in the integral form are in a very good agreement. Also, the same graphs confirm the accuracy of formulas for the asymptotic solution for large values of time presented in Section V. The analogous comparisons for the critically damped and overdamped cases are given in Figures 4,5, respectively.

VII Discussion

In this section we discuss the large-time behavior of energetic quantities in the underdamped case.

From the formal point of view, the exponentially decaying terms of the asymptotic expansions (the ones with superscripts “$(0)$” and “$(2)$”) are much less than the power-decaying terms (the ones with superscript “$(1)$”). It seems, therefore, that the exponentially decaying terms can be dropped out. However, the calculations show that actually the terms with different decay approximate the solution at different timescales.

In Figures 2, 3 one can see that in the underdamped case the unsteady process of thermal equilibration has two stages. At the first stage (“the large times”), when the quantity ${{\mathrm{e}}^{-{\eta}\,t}}$ is not yet very small, the qualitative description of the process is as follows: the kinetic and potential energies oscillate approaching the curvilinear asymptote. The asymptote corresponds to the term $\hat{\Pi}_{0}^{(0)}$ described by formula (88), which is equal to $\frac{{\bar{\mathscr{E}}_{0}}\,{{\mathrm{e}}^{-{\eta}\,t}}}{2}$. The total energy also oscillates (with a smaller amplitude than the kinetic and potential energy) approaching the asymptotic level ${\bar{\mathscr{E}}_{0}}\,{{\mathrm{e}}^{-{\eta}\,t}}$. The Lagrangian oscillates around zero and approaches zero. The oscillatory motions are described by the terms $\tilde{\mathscr{E}}_{0}^{(2)},$ $\tilde{\mathscr{K}}_{0}^{(2)},$ $\tilde{\Pi}_{0}^{(2)},$ $\tilde{\Lambda}_{0}^{(2)}$ expressed by formulas (97), (102), (92), (86), respectively. The amplitudes of oscillations for all energetic quantities are of order $O\left(t^{-1/2}{{\mathrm{e}}^{-{\eta}\,t}}\right)$. In the limiting case of zero dissipation the first stage transforms into the solution describing the thermal equilibration in the corresponding conservative system.

Unexpectedly, there is the second stage (see Figure 3), that can be observed only when the quantity ${{\mathrm{e}}^{-{\eta}\,t}}$ becomes very small 333The stage where the asymptotics have a power decay exists for any positive value of specific viscosity $\eta$ (“the very large times”). The expressions for all energetic quantities at the second stage are subjected to a power decay, i.e. from the formal point of view the corresponding terms $\tilde{\mathscr{E}}_{0}^{(1)},$ $\tilde{\mathscr{K}}_{0}^{(1)},$ $\tilde{\Pi}_{0}^{(1)},$ $\tilde{\Lambda}_{0}^{(1)}$ are principal terms of the corresponding asymptotic expansions. The formulas for these terms are (96), (101), (91), (83), respectively. Another one unexpected result is as follows: the principal term of the asymptotic expansion for $\tilde{\mathscr{K}}_{0}$ is proportional to $t^{-5/2}$ in the case of the kinetic energy and to $t^{-3/2}$ for all other energetic quantities. In the limiting case of zero dissipation the second stage disappears.

The calculations show that in the case $\tilde{\eta}\ll 1$ the asymptotic formulas for the power-decaying terms can give wrong results at a timescale that corresponds to the large, but not very large times. Thus, to approximate the solution at such a timescale, it can be preferable to drop out the power-decaying terms and use only the exponentially decaying terms of asymptotics (the ones with superscripts “$(0)$” and “$(2)$”). This fact is illustrated in Figure 6.

VIII Conclusion

It the paper, we have obtained the analytical solutions in the integral form describing the thermal equilibration in a one-dimensional damped harmonic crystal (see Section IV). These solutions describe the time behavior of the energetic quantities (the total energy, the kinetic energy, the potential energy, and the Lagrangian). The analytical solutions are in an excellent agreement with independent numerical calculations (see Section VI).

The most important result of the paper is large-time asymptotic formulas for energetic quantities presented in Section V (see the derivation in the Appendices A–D).

The main conclusions of the paper can be formulated as follows:

• •

In presence of small viscous external damping (i.e., in the underdamped case) the process of thermal equilibration is more complicated than in the corresponding conservative system and has two stages that correspond to large and very large times (see Figures 2, 3).

• •

In the critically damped (see Figure 4) and the overdamped (see Figure 5) cases the process of thermal equilibration has only one stage, which is similar to the second stage in the underdamped case.

• •

At very large times (i.e., during the second stage) the total energy of an underdamped harmonic crystal is mostly the potential energy. The same conclusion is true for large time behavior of a critically damped and an overdamped harmonic crystal.

Acknowledgments

The authors are grateful to V.A. Kuzkin and E.V. Shishkina for useful and stimulating discussions. This work is supported by Russian Science Foundation (Grant No. 19-41-04106).

Appendix A Calculation of the asymptotics for $\hat{\Lambda}_{n}^{(2)}$

Following to the general procedure of the method of stationary phase Fedoryuk (1977), according to the localization principle, we claim that for $t\to\infty$ the integral $\hat{\Lambda}_{n}^{(2)}$ defined by Eqs. (80), (55) equals the sum

[TABLE]

of contributions $\tilde{\Lambda}_{n}^{(2)}({\mathfrak{q}}_{i})$ from the critical points ${\mathfrak{q}}={\mathfrak{q}}_{i}$ in the interval of integration $[\bar{\mathfrak{q}};\pi]$. In the case under consideration, there are two critical points: the stationary point ${\mathfrak{q}}_{1}=\pi$ of the phase function

[TABLE]

and the boundary of integration interval ${\mathfrak{q}}_{2}=\bar{\mathfrak{q}}$. The contribution from the stationary point ${\mathfrak{q}}_{1}=\pi$ can be calculated as

[TABLE]

Here and in what follows, $\chi({\mathfrak{q}})\in C^{\infty}$ is a non-negative even function such that

[TABLE]

$\bar{\epsilon}>0$ is a small enough number to get the integrand with a unique isolated critical point. The expression for the principal term of the asymptotics for the contribution from an isolated boundary stationary point ${\mathfrak{q}}={\mathfrak{q}}_{\ast}$ is Fedoryuk (1977)

[TABLE]

One has

[TABLE]

Taking the real part of formula (116), wherein $f({\mathfrak{q}})=\lambda({\mathfrak{q}})$,

[TABLE]

to calculate the integral in the right-hand side of Eq. (114) results in

[TABLE]

To finalize the calculation of asymptotics for $\hat{\Lambda}_{n}^{(2)}$ we need to estimate the contribution from the critical point ${\mathfrak{q}}_{2}=\bar{\mathfrak{q}}$. One has

[TABLE]

as ${\mathfrak{q}}\to\bar{\mathfrak{q}}+0$. Here $\phi_{1/2}$ is a positive constant. It follows from the Erdeliy lemma Fedoryuk (1977) that

[TABLE]

where $A$ is a non-zero constant. Applying Eq. (123) in the cases $\beta=1,\ \alpha=1/2$ and $\beta=1/2,\ \alpha=1/2$ yields

[TABLE]

Thus, the asymptotics of $\hat{\Lambda}_{n}^{(2)}$ equals to the right-hand side of Eq. (120).

Appendix B Calculation of the asymptotics for $\hat{\Lambda}_{n}^{(1)}$

One has

[TABLE]

where

[TABLE]

Following to the general procedure of the Laplace method Fedoryuk (1977), we claim that for $t\to\infty$ the each integral $\hat{\Lambda}_{n}^{(1)(\pm)}$ can be asymptotically approximated by contribution

[TABLE]

from the global maximum point ${\mathfrak{q}}_{\pm}$ for the functions $\pm\phi({\mathfrak{q}})$ (defined by (113)) lying in the interval of integration $[-\bar{\mathfrak{q}};\bar{\mathfrak{q}}]$.

At first, consider the integral $\Lambda_{n}^{(1)(+)}$. The maximum point for $\phi({\mathfrak{q}})$ is the internal point ${\mathfrak{q}}=0$, and $\phi(0)=\eta$. The expression for the corresponding contribution in the case of an internal isolated maximum point ${\mathfrak{q}}={\mathfrak{q}}_{\star}$ is (see Fedoryuk (1977), formula (1.25) 444 Note that there are misprints in (Fedoryuk (1977), formula (1.25)) that were corrected in (129). ):

[TABLE]

where

[TABLE]

Now we take $f({\mathfrak{q}})=L_{+}({\mathfrak{q}})$, where

[TABLE]

to calculate the coefficients $c_{k}$. To do this we use Maple symbolic calculation software. In this way, we find

[TABLE]

i.e. we deal with a degenerate case, and

[TABLE]

Thus,

[TABLE]

To finalize the calculation of asymptotics for $\hat{\Lambda}_{n}^{(1)}$ we need to estimate the integral $\tilde{\Lambda}_{n}^{(1)(-)}$. The maximum points for $-\phi({\mathfrak{q}})$ are boundary points $\pm\bar{\mathfrak{q}}$, and $-\phi(\bar{\mathfrak{q}})=0$. Due to the symmetry, these two points bring equal contributions, therefore we can estimate only one of them at ${\mathfrak{q}}=\bar{\mathfrak{q}}$. One has

[TABLE]

as ${\mathfrak{q}}\to\bar{\mathfrak{q}}-0$. Here $L_{1/2}$ is a non-zero constant. It follows from the Watson lemma Fedoryuk (1977) that

[TABLE]

where $B$ is a non-zero constant. Taking into account (121) and using (136) in the case $\beta=1/2,\ \alpha=1/2$ one gets

[TABLE]

Thus, the asymptotics of $\hat{\Lambda}_{n}^{(1)}$ equals to the right-hand side of Eq. (134).

Appendix C Calculation of the asymptotics for $\hat{\Pi}_{n}$

The integral in the right-hand side of Eq. (66) can be represented as follows:

[TABLE]

where ${\epsilon}>0$ is a small enough number. Again, the integral $\hat{\Pi}_{n}^{(1)}$ is a Laplace type integral, whereas the integral $\hat{\Pi}_{n}^{(2)}$ is a Fourier type integral.

The asymptotic expansion for large time of the integral $\Pi_{n}^{(2)}$ is the sum of the doubled contribution from the boundary stationary point ${\mathfrak{q}}={\mathfrak{q}}_{\ast}=\pi$ (where $\phi^{\prime}({\mathfrak{q}}_{\ast})=0$):

[TABLE]

and the doubled contribution from the boundary point ${\mathfrak{q}}=\bar{\mathfrak{q}}+{\epsilon}$. The latter term is discussed at the end of this Appendix (after formula (157)), where we deal with the estimation of the reminder $R$. Applying formula (116), wherein $f({\mathfrak{q}})=P({\mathfrak{q}})$,

[TABLE]

to calculate the integral in the right-hand side of Eq. (145), and taking into account Eqs. (117), (118), results in

[TABLE]

We represent integral $\hat{\Pi}_{n}^{(1)}$ as follows:

[TABLE]

where

[TABLE]

At first, consider the integral $\Pi_{n}^{(1)(+)}$. Again (see Appendix  B), the maximum point for $\phi({\mathfrak{q}})$ defined by Eq. (113) is ${\mathfrak{q}}=0$, and $\phi(0)=\eta$. Applying formulas (128), (129) (wherein $f({\mathfrak{q}})=P({\mathfrak{q}})$, and $P({\mathfrak{q}})$ is defined by (146)), to calculate the corresponding contribution, one can get

[TABLE]

i.e. we deal with a degenerate case, and

[TABLE]

Here, to calculate coefficients $c_{k}$, we again used Maple symbolic calculation software. In this way we find

[TABLE]

Now we need to consider the integral $\Pi_{n}^{(1)(-)}$. The maximum points for $-\phi({\mathfrak{q}})$ defined by Eq. (113) are boundary points ${\mathfrak{q}}=\pm(\bar{\mathfrak{q}}\mp{\epsilon})$:

[TABLE]

The corresponding contribution is discussed in what follows (after formula (157)), where we deal with the estimation of the reminder term $R$.

Now consider the reminder term $R$ defined by Eq. (144). We need to estimate the contribution to integral (138) from a neighbourhood of the point ${\mathfrak{q}}=\bar{\mathfrak{q}}$. We asymptotically approximate integrand in the neighbourhood of ${\mathfrak{q}}=\bar{\mathfrak{q}}$ using Eq. (121), and

[TABLE]

Accordingly, one has

[TABLE]

where

[TABLE]

Here $\phi_{1/2}$ is defined by (121), $\operatorname{Ci}(\cdot)$ is the integral cosine, $\operatorname{Chi}(\cdot)$ is the integral hyperbolic cosine Abramowitz and Stegun (1972). Thus, for large $t$ the reminder $R$ is the sum of contributions from the boundary points $\bar{\mathfrak{q}}\pm{\epsilon}$, which must be totally compensated in the sum with the corresponding contributions from the boundary points for integrals $\hat{\Pi}_{n}^{(1)(\pm)}$ and $\hat{\Pi}_{n}^{(2)(\pm)}$. Hence, to estimate the contribution to integral (138) from a neighbourhood of the point ${\mathfrak{q}}=\bar{\mathfrak{q}}$ we need to take into account the correction terms in expansions (121), (155).

We introduce the substitution

[TABLE]

where $\phi({\mathfrak{q}})$ is defined by Eq. (113) (and expansion (122)). One has

[TABLE]

Analogously, the next term $R_{0}$ in (156) is the sum of contribution of the boundary points (which again are totally compensated) and a contribution from ${\mathfrak{q}}=\bar{\mathfrak{q}}$. Taking into account (155) and (159), one can estimate the latter term as follows:

[TABLE]

The last formula is obtained using the Erdeliy lemma and Watson lemma in the case $\beta=1,\ \alpha=1/2$ and taking $\bar{\epsilon}<\epsilon$ in Eq. (115).

Thus, the asymptotics of $\hat{\Pi}_{n}$ equals the sum of the right-hand side of Eq. (147), the right-hand side of Eq. (153), and the term $\hat{\Pi}_{n}^{(0)}$ calculated in Appendix D in the particular case $n=0$.

Appendix D Calculation of the integral $\hat{\Pi}_{0}^{(0)}$

Since $\bar{\epsilon}$ introduced by Eqs. (138)–(144) is an arbitrary positive number, the integral $\hat{\Pi}_{0}^{(0)}$ can be calculated as

[TABLE]

where symbol $\operatorname{PV}$ means the Cauchy principal value for the corresponding improper integral.

In what follows, we consider only the case $n=0$. According to Eqs. (63), (51) one gets:

[TABLE]

Thus,

[TABLE]

where $\mathscr{J}({\mathfrak{q}})$ is the following integral that can be calculated in the closed form Prudnikov et al. (1986)

[TABLE]

One can see that

[TABLE]

whereas at ${\mathfrak{q}}=\bar{\mathfrak{q}}$ the argument \big{|}\mathscr{Z}({\mathfrak{q}})\big{|} of logarithmic function has indeterminate form $0/0$. One can obtain the following expansions ($\epsilon\to+0$):

[TABLE]

Now it follows from Eqs. (170), (167) that

[TABLE]

and

[TABLE]