Transition to thermal equilibrium in a deformed crystal

TL;DR

This paper analyzes how an infinite particle chain transitions between equilibrium states with different stiffnesses, revealing oscillatory and non-monotonic temperature evolution through analytical and numerical methods.

Contribution

It provides analytical expressions for temperature evolution during a sudden parameter change in a particle chain, highlighting unexpected oscillatory behavior.

Findings

Transition process is oscillatory and non-monotonic.

Analytical expressions for temperature as a function of time.

Numerical simulations support analytical results.

Abstract

An adiabatic transition between two equilibrium states corresponding to different stiffnesses in an infinite chain of particles is studied. Initially, the chain particles have random displacements and random velocities corresponding to a uniform initial temperature. An instant change of parameters of chain initiates a transient process. Analytical expressions for the chain temperature as a function of time are obtained from the statistical analysis of the dynamics equations. It is shown that the transition process is oscillatory and it converges non-monotonically to a new equilibrium state. Such behavior is usually unexpected for thermal processes. The analytical results are supplemented by numerical simulations.

Equations

$$\begin{array}{r}K\stackrel{\text{mbox}def}{=}\frac{1}{2}m⟨v_n^2⟩,T\stackrel{\text{mbox}def}{=}\frac{2K}{k_B},\end{array}$$

$$\begin{array}{rl}& U(ϵ+{\epsilon }_n)\stackrel{\text{mbox}def}{=}⟨\Pi (ϵ+{\epsilon }_n)⟩,{\epsilon }_n\stackrel{\text{mbox}def}{=}{u}_n-u_{n-1},\\ & \Pi (\epsilon )\stackrel{\text{mbox}def}{=}\frac{C}{2}{\epsilon }^2+\frac{\alpha }{3}{\epsilon }^3,\end{array}$$

$$ϵ(t)=\{\begin{array}{ll}0,& t<0\\ ϵ,& t\ge 0,\end{array}$$

$$F_n\stackrel{\text{mbox}def}{=}-F(ϵ+{\epsilon }_n),F(\epsilon )=-{\Pi }^{\prime }(\epsilon )=-C\epsilon -\alpha {\epsilon }^2,$$

$$m{\dot{v}}_n=F_{n+1}-F_n.$$

$$\begin{array}{rl}& K=\frac{1}{2}m⟨v_n^2⟩,U=\frac{C}{2}⟨{\epsilon }_n^2⟩+\frac{\alpha }{3}⟨{\epsilon }_n^3⟩\end{array}$$

$$K=\frac{1}{2}⟨{\epsilon }_n{F}_n⟩.$$

$$m⟨v_n^2⟩=⟨{\epsilon }_n{F}_n⟩=C⟨{\epsilon }_n^2⟩+\alpha ⟨{\epsilon }_n^3⟩.$$

$$k_B{T}_0=m⟨v_n^2⟩=C⟨{\epsilon }_n^2⟩.$$

$$\begin{array}{rl}& K=\frac{1}{2}m⟨v_n^2⟩,U=\Pi (ϵ)+U_T,\\ & U_T=\frac{1}{2}\left(C+2\alpha ϵ\right)⟨{\epsilon }_n^2⟩+\frac{1}{3}\alpha ⟨{\epsilon }_n^3⟩,\end{array}$$

$$\begin{array}{r}{k_{B}T|}_{t\to \infty }=2K=U_T+K.\end{array}$$

$$\begin{array}{r}{k}_{B}T={\frac{1}{2}m⟨v^2⟩|}_{t=0}+{\frac{1}{2}(C+2\alpha ϵ)⟨{\epsilon }_n^2⟩|}_{t=0},\end{array}$$

$$\begin{array}{r}{T|}_{t\to \infty }=T_0\left(1+\frac{\alpha ϵ}{C}\right).\end{array}$$

$$m{\dot{v}}_n=C\Delta {\epsilon }_n+2\alpha ϵ\Delta {\epsilon }_n+\alpha \Delta {\epsilon }_n^2,{\dot{\epsilon }}_n=\Delta v_n.$$

$${\dot{v}}_n={\omega }^2\Delta {\epsilon }_n,{\dot{\epsilon }}_n=\Delta v_n$$

$${v_n|}_{t=0}=\sqrt{\frac{k_B{T}_0}{m}}{\rho }_n,{{\epsilon }_n|}_{t=0}=\sqrt{\frac{k_B{T}_0}{C}}{\varrho }_n,$$

$$\begin{array}{rl}& {\mathcal{K}}_k\stackrel{\text{mbox}def}{=}\frac{1}{2}m⟨v_n{v}_{n+k}⟩,{\mathcal{U}}_k\stackrel{\text{mbox}def}{=}\frac{1}{2}m{\omega }^2⟨{\epsilon }_n{\epsilon }_{n+k}⟩,\\ & {\mathcal{L}}_k\stackrel{\text{mbox}def}{=}{\mathcal{K}}_k-{\mathcal{U}}_k,\end{array}$$

$$\begin{array}{rl}& {\ddot{\mathcal{L}}}_k=4{\omega }^2({\mathcal{L}}_{k-1}-2{\mathcal{L}}_k+{\mathcal{L}}_{k+1}),\\ & t=0:{\mathcal{L}}_k=-\frac{T_0{k}_B\alpha ϵ}{C}{\delta }_k,{\dot{\mathcal{L}}}_k=0,\end{array}$$

$${\mathcal{L}}_k=-\frac{T_0{k}_Bϵ\alpha }{C}{J}_{2k}(4\omega t),$$

$$k_{B}T={\mathcal{L}}_0+K+U_T.$$

$$T=T_{0}+\frac{T_{0}\alpha\epsilon}{C}\Big{(}1-J_{0}(4\omega t)\Big{)}.$$

$$\begin{array}{rl}& x\gg \mu +1:\\ & J_{\mu }(x)=\sqrt{\frac{2}{\pi x}}\cos \left(x-\frac{\pi \mu }{2}-\frac{\pi }{4}\right)+O(x^{-3/2}),\end{array}$$

$$\alpha =-\frac{\beta a{C}^2}{k_B},$$

$$K=\frac{m}{2}⟨u_n{v}_n⟩\dot{}-\frac{1}{2}⟨u_n(F_{n+1}-F_n)⟩,$$

$$⟨u_n(F_{n+1}-F_n)⟩=-⟨{\epsilon }_n{F}_n⟩+⟨F_{n+1}{u}_n⟩-⟨F_n{u}_{n-1}⟩.$$

$$K=\frac{1}{2}⟨{\epsilon }_n{F}_n⟩-\frac{1}{2}{g}^{\prime }+\frac{m}{2}⟨u_n{v}_n⟩\dot{}.$$

$$K=\frac{1}{2}⟨{\epsilon }_n{F}_n⟩.$$

Full text

I Introduction

Nonequilibrium thermal processes in solids at nano- and microlevel are currently a subject of intensive research, partly driven by development of nano technologies [1, 2, 3, 4, 5, 6]. At microlevel, a transition to an equilibrium state for nearly harmonic crystals is a gradual equalization of the kinetic and potential energies of the atomic motion according to the virial theorem [7, 8, 11, 9, 10]. However, this theorem does not describe the processes that occur during the transients. A macroscopic description of such transients is also challenging since it requires application of special constitutive equations for ultrafast atomic processes. Therefore, we develop a simple model that covers analytically both microscopic and macroscopic descriptions of the transient processes.

Crystals with simple lattices are convenient models for studying nonequilibrium processes in solids [12, 13, 14, 15]. Numerical simulations [16] showed that the process of energies equalization in molecular systems was accompanied by high-frequency oscillations. In pioneering paper by Klein and Prigogine [17], the equations of the atomic motion for a one-dimensional harmonic crystal were solved directly and it was shown that the energy oscillations after an instantaneous thermal perturbation were described by the Bessel function of the first kind. In later work [18], this problem was solved by analyzing the dynamics equations for the velocity covariances, which allowed to generalize these results for more complex systems, including multidimensional crystals [20, 19, 21, 22, 23, 24].

In the mentioned works the nanoscale thermal processes were studied separately from the mechanical processes. The advances of the modern technologies bring to attention ultrafast mechanical processes where the speed of mechanical load is comparable or even faster than the speed of the local thermal equilibration in the system. This condition is fulfilled if the material is rapidly deformed by forces uniformly distributed along the length of the sample. Such loads occur in nanoscale electronic components of experimental equipment requiring fast magnetic field switches, for example, necessary for condensed matter physics, plasma physics, or inertial confinement synthesis [25]. Electric pulses can create distributed electromagnetic loads on circle samples for the time up tens nanoseconds [26, 27]. In [28], the concepts of an ultrafast photodetector capable of converting femtosecond light pulses into electric pulses of the same length is proposed. Paper [25] demonstrates an all-optical method to generate magnetic field impulses of the order of several tesla over the course of tens of femtoseconds. In the nearest future it is expected that the experimental electromagnetic impulses can reach a duration less than femtosecond: the femtosecond lasers already exist [29] and attoseconds lasers are underdevelopment [30, 31, 32]. Therefore, an analytical study of the impact of the ultrafast mechanical loads on the thermal equilibration is needed to provide a theoretical basis for the upcoming experimental studies.

In the present paper, we study an adiabatic non-equilibrium process analogous to those considered previously [20, 19, 21, 22, 23, 24]. However, this process is initiated by an instantaneous external load instead of instantaneous heating. The material is deformed by forces uniformly distributed along the length of the sample [26, 27] hence for this case the orders of velocity of mechanical and thermal processes are the same. The interaction between mechanical and thermal processes requires accounting non-linearity, therefore $\alpha$-FPU one-dimensional crystal [33, 34] is considered. The non-linearity is assumed to be large enough to cause an adiabatic heating of the crystal, but, on the other hand, is sufficiently small to analysis the resulting energy oscillations in the framework of the harmonic approximation.

The paper is organized as follows. In Sect.II the mathematical formulation of the problem is presented. In Sect.III the equilibrium value of the kinetic temperature after loading is obtained. In Sect.IV it is shown that the crystal temperature oscillates during the transition process and an analytical expression that describes these oscillations is obtained. In the particular system considered, the kinetic and potential energies gradually become equal over time (in agreement with the virial theorem). The analytical solution is compared with a numerical one. In Sect.V an estimation of the temperature jump for the real crystal is represented.

II Formulation of the problem

A model of an infinite one-dimensional crystal — chain of point masses connected by unmasses springs is considered. It is assumed that the chain particles interact only with their nearest neighbours. In thermodynamic equilibrium state particle velocities $v_{n}$ and bond deformations $\varepsilon_{n}$ are independent stochastic quantities and all statistical characteristics of the crystal are constant over time. To describe the statistical behaviour of the system the following quantities are introduced:

[TABLE]

where $K$ is the mathematical expectations of kinetic energy, $T$ is the kinetic temperature, $n=1,2,3,...$ is the index of the particle, $\langle{...}\rangle$ is the operator of the mathematical expectation, $v_{n}$ is the particle velocity, $k_{B}$ is the Boltzmann constant, and $m$ is the particle mass. The statistical characteristics of this system depend on such crystal parameters as particle mass and bond stiffness. Variation of these parameters transfers the crystal to a non-equilibrium state and the transition process starts, which eventually brings the crystal to a new equilibrium state.

One of the natural way to change the bond stiffness of the crystal is a homogeneous crystal deformation by external loading. For example, such loading can be realized by applying distributed electromagnetic forces to a circle samples [26, 27]. Based on the approach described [18], we investigate the evolution of the kinetic temperature of the crystal during this transition process.

In addition to mathematical expectation of kinetic energy $K$, this paper uses mathematical expectation of potential energy $U$:

[TABLE]

where $\varPi$ is the inter-particle potential, $C$ and $\alpha$ are stiffness of the inter-particle bonds of the first and the second order, $u_{n}$ is the displacement of the particle from its equilibrium position, and deformation is a sum of homogeneous deformation $\epsilon$ and stochastic deformation $\varepsilon_{n}$. Homogeneous deformation $\epsilon(t)$ is applied instantaneously:

[TABLE]

where $\epsilon$ is a constant. The force $F_{n}$ acting on the particle $n-1$ by the particle $n$ is

[TABLE]

and equation of particle dynamics of the crystal is:

[TABLE]

We consider small non-linearity, therefore the second terms in the formulas for $\varPi(\varepsilon)$ (2) and $F(\varepsilon)$ (4) are small compare to the first one. Since the non-linear term is small, a harmonic (linear) approximation is used to describe the thermal processes. At relatively short times, the harmonic approximation is quite accurate for describing thermal processes in a crystal with small non-linearity [20]. However, non-zero stiffness of the second order allows taking into account the influence of homogeneous deformation $\epsilon$ on inter-particle potential.

The analytical expression for kinetic temperature at long times after homogeneous deformation can be obtained using the virial theorem and is given in the following section.

III Crystal in the thermodynamic equilibrium

This section contains expressions for the kinetic temperature of the crystal before loading and at large times after loading.

III.1 Prior to loading

At $t<0$ according to (3), the energies are

[TABLE]

and the crystal is assumed to be in thermodynamic equilibrium. Since only the mathematical expectations of the energies are considered hereinafter, the words “mathematical expectation” are omitted. Following [8] the kinetic energy (1) of the system can be represented by (see appendix A):

[TABLE]

Formula (7) allows one to obtain expressions for the equilibrium kinetic temperature of the crystal before and after the instantaneous deformation, as follows:

[TABLE]

We neglect the small term ${\alpha}\langle{\varepsilon_{n}^{3}}\rangle$ in equation (8) and obtain the kinetic temperature of the crystal prior to loading

[TABLE]

In the next subsection, we consider the thermodynamic quantities after loading.

III.2 After loading

At $t\geqslant 0$, the kinetic and potential energies of the crystal:

[TABLE]

where $U_{T}$ is the thermal part of the potential energy. In the state of thermodynamic equilibrium the kinetic energy $K$ and the thermal part of the potential energy $U_{T}$ are equal, therefore expression for kinetic temperature (1) is:

[TABLE]

Note that $U_{T}+K$ is constant for $t\geq 0$ in the transition process. An expression for this sum can be obtained by substitution to (11) of the initial values of the corresponding quantities:

[TABLE]

where term $\alpha\langle{\varepsilon_{n}^{3}}\rangle/3$ is omitted. Substituting expressions (9) into (12) yields:

[TABLE]

From expression (13) it follows that the change in kinetic temperature $T$, in a first approximation, is proportional to deformation $\epsilon$.

IV Dynamics of the transition process

In this section we obtain an expression for the temperature as a function of time. Substituting expression (4) into (5) we obtain the following equation of motion:

[TABLE]

Term $\alpha\Delta\varepsilon_{n}^{2}$ can be neglected in the case of small deformations. Thus equations (14) become linear:

[TABLE]

where $\omega\stackrel{{\scriptstyle\mbox{\scriptsize def}}}{{=}}\sqrt{\frac{C+2\alpha\epsilon}{m}}$. Note that sum $C+2\alpha\epsilon$ plays a role the first order stiffens after external loading. The initial conditions for the system (15) are determined from equation (9):

[TABLE]

where $\rho_{n}$ and $\varrho_{n}$ are independent random numbers with zero mathematical expectation and unit variance. The initial value problem (15)-(16) describes the stochastic dynamics of the chain particles. Then the kinetic temperature of the crystal as a function of time can be found using the covariance analysis approach [18, 22], by introducing generalized energies

[TABLE]

where ${\cal K}_{k}$ and ${\cal U}_{k}$ are generalized kinetic and potential energies, and ${\cal L}_{k}$ is generalized Lagrangian. Differentiating of (17) and equations of motion (14) lead to the following initial problem for the generalized Lagrangian:

[TABLE]

where $\delta_{k}=1$ for $k=0$ and $\delta_{k}=0$ otherwise. Solution of a similar initial problem is obtained in [18]. According to that solution, generalized Lagrangian oscillates with a monotonically decreasing amplitude:

[TABLE]

where $J_{k}(x)$ is the Bessel function of $k$-th order [35].

Kinetic energy $K$ and thermal part of the potential energy $U_{T}$ are equal to the generalized kinetic and potential energies with zero indexes, therefore kinetic temperature $T$ can be expressed as follows:

[TABLE]

Then the kinetic temperature of the crystal as a function of time can be found using formulas (13) and (19)-(20):

[TABLE]

Note, that according the asymptotic representation for the Bessel function [35]

[TABLE]

the amplitude of the kinetic temperature oscillations decreases as $t^{-1/2}$.

Fig. 1 shows comparison of the analytical solution (21) and the numerical solution obtained by computer simulations of the crystal dynamics consisting of $N=5\cdot 10^{4}$ particles under periodic boundary conditions. In the framework of the numerical experiment, the parameters of the problem under consideration are chosen so that $\alpha\epsilon/C=-0.1$. The simulations use the method of central differences and integration step $0.01/\omega_{e}$. At the initial moment of time, the particle displacements are zero, and the particle velocities are random and correspond to the crystal temperature $2T_{0}$. The process of energies equalization results the crystal temperature oscillates with a decreasing amplitude around $T_{0}$ value. The homogeneous deformation is applied when the temperature oscillations have a negligibly small amplitude. After loading of the crystal, the temperature oscillates around a new equilibrium value.

In [20] it is shown numerically that thermal phenomena for crystals with a sufficiently low nonlinearity do not differ greatly from harmonic crystals. For the simulation, a time span sufficient to describe tens of temperature oscillations after homogeneous deformation has been chosen. However, the time span is chosen not to be large so that temperature oscillations have a significant deviation from the harmonic solution. To calculate the mathematical expectations for the statistical quantities the results are averaged over all particles and $10^{3}$ realizations, which are solutions of the same equations with different random number generations for the initial conditions. According to work [18], crystal temperature before the instantaneous deformation is the Bessel function of zero order. Instant deformation is applied when the amplitude of temperature oscillation is small compared to absolute value of $\alpha\epsilon/C$.

As seen from Fig. 1 analytical solution (21) practically coincides with the results of the numerical integration of the chain dynamics equations (3)–(5) for several tens of oscillation periods.

Formula (13) gives the limiting value for the kinetic temperature. According to expression (20), after the instantaneous deformation, the kinetic temperature oscillates around this limiting value and tends to it for the large times. Thus, for $t\to\infty$ expression (21) coincides with formula (13), obtained from the virial theorem. For arbitrary times this expression gives the desired description of the nonequilibrium transition process.

V Example

In order to estimate the temperature jump in the transition process, we consider one-dimensional ring of carbon atoms, which is in a thermodynamic equilibrium state at initial temperature $T_{0}=300^{o}$ K. We assume that in the result of the homogeneous loading the bond deformation is $1.0$% of the equilibrium interparticle distance $a=0.154$ nm [36]. The mass of the carbon atom $m=1.99\cdot 10^{-26}$ kg and the first-order stiffness coefficient is taken equal to the stiffness of the diamond bond $C=472$ N/m [37]. The second-order stiffness coefficient can be found from the following formula [38]:

[TABLE]

where $\beta=0.7\cdot 10^{-6}~{}{\rm K}^{-1}$ is the coefficient of the thermal expansion of the diamond crystal [39]. The substitution of the stiffness coefficients and deformation into formula (13) gives the new equilibrium temperature value of $298.3$ K. The asymptotic period of the kinetic temperature oscillations in the transient process described by formula (21), is approximately $10.2$ femtoseconds. Thus $1.0$% deformation leads to a change in the crystal temperature by $0.6$%.

VI Conclusions

The paper presents an analytical approach for the analysis of the transition process in one-dimensional crystals (chains) subjected to an instantaneous homogeneous deformation. Such a deformation can be interpreted as an instantaneous change in the stiffness of the inter-particle bonds in a chain. It is obtained that the transition process is accompanied by high-frequency energy oscillations, which have an analytical representation in the terms of the zero Bessel function of the first kind, consequently the amplitude of the transitional oscillations is inversely proportional to the square root of time. After the decay of the oscillations the system reaches a near equilibrium state corresponding to the predictions of equilibrium thermodynamics, however the transition process can be studied in details with the use of the presented approach. The analytical solution is confirmed by numerical simulations. Using technique described in [19] the presented approach can be extended to analyses transition processes in two-dimensional and three-dimensional materials. The obtained results are important for establishing link between mechanical and thermal processes in solids at the femtosecond time scale.

Acknowledgements.

The authors would like to express their gratitude to professor A. V. Porubov, professor M. L. Kachanov for the useful discussions and M. A. Bolshanina and N. D. Mushchak for help in preparing the paper.

This work was supported by the Russian Science Foundation (No. 19-41-04106).

Appendixes

Appendix A The virial relation

Following [8] the kinetic energy (1) of the system can be represented by

[TABLE]

where expression (5) is used. The second term in (24) is

[TABLE]

Defining $g\stackrel{{\scriptstyle\mbox{\scriptsize def}}}{{=}}\langle{F_{n}u_{n-1}}\rangle$, equation (24) takes the form

[TABLE]

Values $g$ and $\langle{u_{n}v_{n}}\rangle$ are constant in the thermodynamic equilibrium state and their derivatives are zero. Therefore the expression for the kinetic energy is

[TABLE]

Formula (27) allows one to obtain expressions for the equilibrium kinetic temperature of the crystal before and after the instantaneous deformation.