Steady-state kinetic temperature distribution in a two-dimensional square harmonic scalar lattice lying in a viscous environment and subjected to a point heat source

TL;DR

This paper models heat transfer in a 2D harmonic lattice in a viscous environment with a point heat source, deriving analytical solutions for the steady-state temperature distribution and validating them against numerical simulations.

Contribution

It introduces a new analytical formula for the steady-state kinetic temperature distribution in a 2D lattice with a point heat source, derived from stochastic differential equations.

Findings

Analytical solution matches numerical results well except near the source.

Derived differential-difference equation describes non-stationary heat propagation.

Singularity occurs in the analytical solution at the source point.

Abstract

We consider heat transfer in an infinite two-dimensional square harmonic scalar lattice lying in a viscous environment and subjected to a heat source. The basic equations for the particles of the lattice are stated in the form of a system of stochastic ordinary differential equations. We perform a continualization procedure and derive an infinite system of linear partial differential equations for covariance variables. The most important results of the paper are the deterministic differential-difference equation describing non-stationary heat propagation in the lattice and the analytical formula in the integral form for its steady-state solution describing kinetic temperature distribution caused by a point heat source of a constant intensity. The comparison between numerical solution of stochastic equations and obtained analytical solution demonstrates a very good agreement everywhere except for the main diagonals of the lattice (with respect to the point source position), where the analytical solution is singular.

Equations

$$\mathrm{d}{v}_{i,j}=F_{i,j}\mathrm{d}t+b_{i,j}\mathrm{d}{W}_{i,j},\mathrm{d}{u}_{i,j}=v_{i,j}\mathrm{d}t,$$

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

$$\mathrm{d}{W}_{i,j}={\rho }_{i,j}\sqrt{\mathrm{d}t},$$

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

$${\mathcal{L}}_i{u}_{i,j}=u_{i+1,j}-2u_{i,j}+u_{i-1,j},$$

$${\mathcal{L}}_{i,j}{u}_{i,j}={\mathcal{L}}_i{u}_{i,j}+{\mathcal{L}}_j{u}_{i,j}$$

$$⟨{\rho }_{i,j}⟩=0,⟨{\rho }_{i,j}{\rho }_{k,l}⟩={\delta }_{i;k}{\delta }_{j;l}$$

$$u_{i,j}(0)=0,v_{i,j}(0)=0.$$

$$b_{i,j}\equiv 0\text{if}\max \{i,j\}>R$$

$$⟨u_{i,j}⟩=0,⟨v_{i,j}⟩=0.$$

$${\xi }_{p,q;r,s}\stackrel{\text{mbox}def}{=}⟨u_{p,q}{u}_{r,s}⟩,{\nu }_{p,q;r,s}\stackrel{\text{mbox}def}{=}⟨u_{p,q}{v}_{r,s}⟩,{\kappa }_{p,q;r,s}\stackrel{\text{mbox}def}{=}⟨v_{p,q}{v}_{r,s}⟩,$$

$${\beta }_{p,q;r,s}\stackrel{\text{mbox}def}{=}{\delta }_{p;r}{\delta }_{q;s}{b}_{p,q}{b}_{r,s}.$$

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

$$2{\mathcal{L}}^{\mathrm{S}}\stackrel{\text{mbox}def}{=}{\mathcal{L}}_{p,q}+{\mathcal{L}}_{r,s},2{\mathcal{L}}^{\mathrm{A}}\stackrel{\text{mbox}def}{=}{\mathcal{L}}_{p,q}-{\mathcal{L}}_{r,s},$$

$$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.$$

$$\kappa\big{|}_{t<0}\equiv 0.$$

$$k\stackrel{\text{mbox}def}{=}p+r,l\stackrel{\text{mbox}def}{=}q+s$$

$$m\stackrel{\text{mbox}def}{=}r-p,n\stackrel{\text{mbox}def}{=}s-q.$$

$$\mathbf{x}\stackrel{\text{mbox}def}{=}\frac{a}{2}(k\mathbf{i}+l\mathbf{j}),$$

$$\mathbf{x}\stackrel{\text{mbox}def}{=}{x}_1\mathbf{i}+x_2\mathbf{j},$$

$$a=ϵ\bar{a},$$

$$c\stackrel{\text{mbox}def}{=}a{\omega }_0$$

$${\omega }_0={ϵ}^{-1}{\bar{\omega }}_0,$$

$${\hat{\zeta }}_{m,n}(\mathbf{x})={\zeta }_{p,q;r,s}.$$

$${\mathcal{L}}_{p,q}{\zeta }_{p,q;r,s}=$$

$$+$$

$${\mathcal{L}}_{r,s}{\zeta }_{p,q;r,s}=$$

$$+$$

$${\mathcal{L}}_{p,q}{\zeta }_{p,q;r,s}={\mathcal{L}}_{m,n}{\hat{\zeta }}_{m,n}$$

$$+$$

$${\mathcal{L}}_{r,s}{\zeta }_{p,q;r,s}={\mathcal{L}}_{m,n}{\hat{\zeta }}_{m,n}$$

$$+$$

Full text

††journal: CMAT

1 Introduction

At the macroscale, Fourier’s law of heat conduction is widely and successfully used to describe heat transfer processes. However, recent experimental observations demonstrate that Fourier’s law is violated at the microscale and nanoscale, in particular, in low-dimensional nanostructures chang2008breakdown ; xu2014length ; hsiao2015micron ; cahill2003nanoscale ; liu2012anomalous ; lepri2016thermal8 , where the ballistic heat transfer is realized. The anomalous heat transfer also may be related with the spontaneous emergence of long-range correlations; the latter is typical for momentum-conserving systems lepri2003thermal ; spohn2016fluctuating . The simplest theoretical approach to describe the ballistic heat propagation is to use harmonic lattice models. In some cases such models allow one to obtain the analytical description of thermomechanical processes in solids hoover2015simulation ; daly2002molecular ; krivtsov2003nonlinear ; berinskii2016elastic ; kuzkin2016lattice ; berinskii2015linear ; berinskii2016hyperboloid . In the literature, the problems concerning heat transfer in harmonic lattices are mostly considered in the stationary formulation bonetto2000mathematical ; lepri2003thermal ; 0295-5075-43-3-271 ; dhar2008heat ; rieder1967properties ; allen1969energy ; nakazawa1970lattice ; lee2005heat ; kundu2010heat ; lepri2016thermal2 ; bernardin2012harmonic ; freitas2014analytic ; freitas2014erratum ; hoover2013hamiltonian ; lukkarinen2016harmonic , the non-stationary heat propagation is discussed in le2008molecular ; gendelman2012nonstationary ; tsai1976molecular ; ladd1986lattice ; volz1996transient ; daly2002molecular ; gendelman2010nonstationary ; guzev2018 .

In previous studies krivtsov2014energy ; krivtsov2015heat ; krivtsov-da70 , a new approach was suggested which allows one to solve analytically non-stationary thermal problems for an infinite one-dimensional harmonic crystal — an infinite ordered chain of identical material particles, interacting via linear (harmonic) forces. In particular, a heat transfer equation was obtained that differs from the extended heat transfer equations suggested earlier chandrasekharalah1986thermoelasticity ; tzou2014macro ; cattaneo1958forme ; vernotte1958paradoxes ; however, it is in an excellent agreement with molecular dynamics simulations and previous analytical estimates gendelman2012nonstationary . The properties of the solutions describing heat transfer in a one-dimensional harmonic crystal were discussed in sokolov2017localized ; krivtsov2018one ; PhysRevE.99.042107 . Later this approach was generalized babenkov2016energy ; krivtsov2014energy ; krivtsov2015heat ; kuzkin2017analytical ; kuzkin2017high ; Kuzkin-Krivtsov-accepted ; murachev2018thermal to a number of systems, namely, to an infinite one-dimensional crystal on an elastic substrate babenkov2016energy , to an infinite one-dimensional diatomic harmonic crystal podolskaya2018anomalous , to a finite one-dimensional crystal murachev2018thermal , and to two and three-dimensional infinite harmonic lattices kuzkin2017analytical ; Kuzkin-Krivtsov-accepted ; kuzkin2017high . In the most of the above mentioned papers babenkov2016energy ; krivtsov2014energy ; krivtsov2015heat ; kuzkin2017analytical ; kuzkin2017high ; Kuzkin-Krivtsov-accepted ; murachev2018thermal ; sokolov2017localized only isolated systems were considered.

In recent paper gavrilov2018heat an infinite one-dimensional harmonic crystal that can exchange energy with its surroundings was considered. It was assumed that the crystal lies in a viscous environment (a gas or a liquid, e.g., the air) which causes an additional dissipative term in the equations of stochastic dynamics for the particles. Additionally, sources of heat supply were taken into account as an additional noise term in equations of motions. Unsteady heat conduction regimes as well as the steady-state kinetic temperature distribution caused by a point heat source of a constant intensity were investigated.

In the present paper we generalize to the two-dimensional case the results obtained in paper gavrilov2018heat . We use the model of an infinite two-dimensional harmonic square scalar111This means that the displacements are scalar quantities mielke2006macroscopic ; harris2008energy ; savin2016normal ; nishiguchi1992thermal ; Kuzkin-Krivtsov-accepted lattice performing transverse oscillations. The schematic of the system is shown in Figure 1.

The aim of the paper is to develop an equation describing non-stationary heat propagation in the lattice and use this equation to obtain analytical formulas describing the steady-state kinetic temperature distribution caused by a point heat source of a constant intensity. The differences in methodology between the present paper and gavrilov2018heat are:

• •

When deriving the differential-difference equation describing non-stationary propagation in the lattice, we avoid to use identities of the calculus of finite differences and instead inverse the corresponding finite difference operators, when necessary, in the Fourier domain. This significantly simplifies the formal procedure. Accordingly, we do not introduce into consideration the quantities, which we call “non-local temperatures” in gavrilov2018heat (Section 4) and deal with covariance variables only.

• •

Since the two-dimensional case is more complicated, the steady-state solution for the kinetic temperature is obtained as a solution of the stationary equations for covariances of the particle velocities (in gavrilov2018heat the steady-state solution was found as a limiting case for the corresponding unsteady solution).

• •

The investigation of unsteady heat conduction regimes is beyond the scope of this paper.

The paper is organized as follows. In Section 2, we consider the formulation of the problem. In Section 2.1, some general notation is introduced. In Section 2.2, we state the basic equations for the lattice particles in the form of a system of stochastic differential equations. In Section 2.3, we introduce and deal with infinite set of covariance variables. These are the mutual covariances of all the particle velocities and all the displacements for all pairs of particles. We use the Itô lemma to derive (see gavrilov2018heat , Appendix A) an infinite deterministic system of ordinary differential equations which follows from the equations of stochastic dynamics. This system can be transformed into an infinite system of differential-difference equations involving only the covariances for the particle velocities. In Section 3, we introduce a vectorial continuous spatial variable and write the finite difference operators involved in the equation for covariances as compositions of finite difference operators and operators of differentiation. In Section 4, we perform an asymptotic uncoupling of the equation for covariances. Provided that the introduced continuous spatial variable can characterize the behavior of the lattice, one can distinguish between slow motions, which are related to the heat propagation, and vanishing fast motions krivtsov2014energy ; babenkov2016energy ; kuzkin2017high ; Kuzkin-Krivtsov-accepted ; kuzkin2017analytical ; kuzkin2019thermal , which are not considered in the paper. Slow motions can be described by a coupled infinite system of second-order hyperbolic partial differential equations for the continualized covariances of velocities. The kinetic temperature is introduced as a quantity proportional to the variance of the particle velocities. In Section 5 we obtain an analytical expression in the integral form for the steady-state solution describing the kinetic temperature distribution caused by a point source of a constant intensity (the corresponding calculations can be found in Appendices A & B). In Section 6, we present the results of the numerical solution of the initial value problem for the system of stochastic differential equations and compare them with the obtained analytical solution. In the conclusion (Section 7), we discuss the basic results of the paper.

2 Mathematical formulation

2.1 Notation

In the paper, we use the following general notation:

$t$

is the time;

$H(\cdot)$

is the Heaviside function;

$\delta(\cdot)$

is the Dirac delta function in a two-dimensional space;

$\delta_{1}(\cdot)$

is the Dirac delta function in a one-dimensional space;

$\langle\cdot\rangle$

is the expected value for a random quantity;

$\delta_{p;q}$

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

$\mathbb{Z}$

is the set of all integers;

$\mathbb{R}$

is the set of all real numbers.

2.2 Stochastic lattice dynamics

Consider the following infinite system of stochastic ordinary differential equations kloeden1999 ; stepanov2013stochastic :

[TABLE]

where

[TABLE]

Here $i,j\in\mathbb{Z}$ are arbitrary integers which describes the position of a particle in the square lattice; the stochastic processes $u_{i,j}(t)$ and $v_{i,j}(t)$ are the displacement and the particle velocity, respectively; $F_{i,j}$ is the force on the particle; $W_{i,j}$ are uncorrelated Wiener processes kloeden1999 ; stepanov2013stochastic ; $b_{i,j}(t)$ is the intensity of the random external excitation; $\eta$ is the specific viscosity for the environment; $C$ is the bond stiffness; $m$ is the mass of a particle; ${\mathcal{L}}_{i,j}$ is the linear finite difference operator describing an infinite two-dimensional stretched square scalar lattice with nearest-neighbor interactions performing out-of-plane vibrations:

[TABLE]

The normal random variables $\rho_{i,j}$ are such that

[TABLE]

and they are assumed to be independent of $u_{i}$ and $v_{i}$. The initial conditions are zero: for all $i$, $j$

[TABLE]

In the case $b_{i,j}\equiv b$, equations (2.1) are the Langevin equations langevin1908theorie ; lemons1997paul for a two-dimensional harmonic stretched square lattice surrounded by a viscous environment (e.g., a gas or a liquid). Assuming that $b_{i,j}$ may depend on ${i,j}$, we introduce a natural generalization of the Langevin equation which allows one to describe the possibility of an external heat excitation (e.g., laser excitation) gavrilov2018heat . This external excitation is assumed to be localized in space:

[TABLE]

for some $R>0$, and much more intensive than the stochastic influence caused by a non-zero temperature of the environment. Therefore, we neglect in (2.1) the stochastic term, which variance does not depend on $i$ and $j$, that models the influence from the environment.222To consider the full problem with a non-zero temperature of the environment we need to take into account the processes of thermal equilibration (the fast motions) being coupled with the heat propagation. The preliminary calculations show that the presence of the viscous environment makes the problem to be much more complicated in comparison with the conservative case considered, e.g., in Kuzkin-Krivtsov-accepted . At the same time, in the framework of this general problem, the processes of heat transfer can be easily separated, and due to linearity of the general problem, they described by exactly the same equations as ones formulated in this paper. Note that zero initial conditions (2.8) and requirement (2.9) guarantee that there are no sources at the infinity, thus we do not need additional boundary conditions at the infinity.

2.3 The dynamics of covariances

According to (2.2), $F_{i,j}$ are linear functions of $u_{i,j},\ v_{i,j}$. Taking this fact into account together with (2.7) and (2.8), we see that for all $t$

[TABLE]

Following krivtsov2016ioffe ; gavrilov2018heat , consider the infinite sets of covariance variables333We use commas to separate the subscripts related with one particle, and semicolon to separate the subscripts related with different particles.

[TABLE]

and the quantities

[TABLE]

In the last formula, we take into account the second formula of (2.7). Thus the variables $\xi_{p,q{;}r,s}$, $\nu_{p,q{;}r,s}$, $\kappa_{p,q{;}r,s}$ are defined for any pair of lattice particles. For simplicity, in what follows we drop all the subscripts, i.e., $\xi\stackrel{{\scriptstyle\mbox{\scriptsize def}}}{{=}}\xi_{p,q{;}r,s}$ etc. By definition, we also put $\xi^{\top}\stackrel{{\scriptstyle\mbox{\scriptsize def}}}{{=}}\xi_{r,s{;}p,q}$ etc. Now we differentiate the variables (2.11) with respect to time taking into account the equations of motion (2.1). This yields the following closed system of differential equations for the covariances (see gavrilov2018heat , Appendix A):

[TABLE]

where $\partial_{t}$ is the operator of differentiation with respect to time; ${\mathcal{L}}_{p,q}$ and ${\mathcal{L}}_{r,s}$ are the linear difference operators defined by (2.6) that act on $\xi_{p,q{;}r,s}$, $\nu_{p,q{;}r,s}$, $\kappa_{p,q{;}r,s}$, $\beta_{p,q{;}r,s}$ with respect to the first pair of subscripts $p,\ q$ and the second one $r,\ s$, respectively. 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 equations (2.13) can be rewritten as follows:

[TABLE]

This system of equations can be reduced (see gavrilov2018heat ) to one equation of the fourth order in time for covariances of the particle velocities $\kappa$:

[TABLE]

In what follows, we deal with Eq. (2.18). This equation can describe a wide class of systems. For example, in the simplest particular case $\mathscr{L}_{i,j}u_{i,j}=-u_{i,j}$ it can be used to obtain the classical results wang1945theory ; stepanov2013stochastic concerning time evolution of the variance of the velocity for a single stochastic oscillator with additive noise.

According to Eqs. (2.8), (2.11), we supplement Eq. (2.18) with zero initial conditions. We state these conditions in the following form, which is conventional for distributions (or generalized functions) Vladimirov1971 :

[TABLE]

To take into account non-zero classical initial conditions, one needs to add the corresponding singular terms (in the form of a linear combination of $\delta(t)$ and its derivatives) to the right-hand sides of the corresponding equations Vladimirov1971 .

Let us note that equation (2.18) is a deterministic equation. What is also important is that (2.18) is a closed equation. Thus, the thermal processes do not depend on any property of the cumulative distribution functions for the displacements and the particle velocities other than the covariance variables used above.

3 Continualization of the finite difference operators

Following krivtsov2015heat ; krivtsov2016ioffe ; gavrilov2018heat , we introduce the discrete spatial variables

[TABLE]

and the discrete correlational variables

[TABLE]

instead of discrete variables $p$, $q$, $r$, $s$. We can also formally introduce the continuous vectorial spatial variable

[TABLE]

where $a$ is the lattice constant (the distance between neighboring particles along a row), $\mathbf{i}$ and $\mathbf{j}$ are orthogonal unit vectors (the directions of the principal axes of the lattice, see Figure 1). In what follows, we represent $\mathbf{x}$ also in the following form:

[TABLE]

where $x_{1},\ x_{2}$ are the spatial coordinates in the lattice.

To perform the continualization, we assume that the lattice constant is a small quantity and introduce a dimensionless formal small parameter $\epsilon$ in the following way:

[TABLE]

where $\bar{a}=O(1)$. To preserve the speed of sound in the lattice

[TABLE]

as a quantity of order $O(1)$, we additionally assume that

[TABLE]

where $\bar{\omega}_{0}=O(1)$. Thus $c=\bar{a}\bar{\omega}_{0}=O(1)$. The basic assumption that allows one to perform the continualization is that any quantity $\zeta_{p,q{;}r,s}$ defined by (2.11) or (2.12) can be calculated as a value of a smooth function $\hat{\zeta}_{m,n}(\mathbf{x})$ of the continuous spatial slowly varying coordinate $\mathbf{x}=\tfrac{1}{2}\epsilon\bar{a}(k\mathbf{i}+l\mathbf{j})$ and the discrete correlational variables $m$ and $n$:

[TABLE]

In accordance with Eqs. (3.1), (3.2), one has

[TABLE]

Applying the Taylor theorem to these formulas yields

[TABLE]

An alternative way of continualization can be realized by letting the number of particles diverge, rather than invoking an increasingly small separation lepri2010nonequilibrium . Despite the algorithmic difference, these approaches apparently lead to the same result, although application of the first approach to an infinite system seems to be more straightforward.

Now we perform the continualization of the operators ${\mathcal{L}}^{\mathrm{S}},\ {\mathcal{L}}^{\mathrm{A}}$. Using (3.10), we obtain

[TABLE]

where

[TABLE]

4 Slow motions

Taking into account assumption (3.7), Eq. (2.18) can be rewritten in the following form:

[TABLE]

Equation (4.1) is a differential equation whose highest derivative with respect to $t$ is multiplied by a small parameter. Therefore, one can expect the existence of two types of solutions, namely, solutions slowly varying in time and fast varying in time nayfeh2008perturbation ; kevorkian2012multiple . The presence of fast and slow motions is a standard property of statistical systems. Fast motions are oscillations of temperature caused by equilibration of kinetic and potential energies. Slow motions are related with macroscopic heat propagation.

Considering slow motions, we assume that

[TABLE]

Vanishing solutions krivtsov2014energy ; babenkov2016energy ; kuzkin2017high ; Kuzkin-Krivtsov-accepted ; kuzkin2017analytical that characterize fast motions, which do not satisfy (4.2), are not considered in this paper.

Now, taking into account Eqs. (3.11), (3.12), we drop the higher order terms and rewrite equation (4.1) in the form of an equation for slow motions:

[TABLE]

where quantity $\hat{\beta}_{m,n}(\mathbf{x})$ is introduced in accordance with Eq. (3.8). Note that according to (2.12), $\hat{\beta}_{m,n}=0$ if $m\neq 0$ or $n\neq 0$, thus

[TABLE]

We identify the following quantities depending on the continuous variable $\mathbf{x}$

[TABLE]

as the kinetic temperature and the heat supply intensity, respectively. Here $k_{B}$ is the Boltzmann constant. We also introduce the following notation:

[TABLE]

Now Eq. (4.3) can be rewritten as the following infinite system of PDE describing non-stationary heat propagation in the lattice:444These equations have a little bit different structure than the corresponding equations obtained in previous paper (gavrilov2018heat , Eq. (4.3)), since the quantities $\theta_{m,n}(x,t)$ ($mn\neq 0$) differ from the quantities, which we call “non-local temperatures” in gavrilov2018heat (see Eq. (4.7)).

[TABLE]

The initial conditions that corresponds to Eq. (2.19) are

[TABLE]

In the particular case where $\eta=0,\ \chi\propto\delta(t)$, system (4.8) was obtained and investigated in Kuzkin-Krivtsov-accepted . In the latter case the derivation of the corresponding system of PDE is much simpler, since it can be based on consideration of system of ordinary differential equations with random initial conditions instead of system of stochastic equations (2.1).

Note that equations (4.8) for slow motions involve only the product $c=\omega_{0}a$ and do not involve the quantities $\omega_{0}$ and $a$ separately, so they do not involve $\epsilon$. Provided that the initial conditions also do not involve $\epsilon$, the solution of the corresponding initial value problem and all its derivatives are quantities of order $O(1)$. The rate of vanishing for fast motions depends on $\epsilon$: the smaller $\epsilon$, the higher the rate. Thus, for sufficiently small $\epsilon$, exact solutions of Eq. (2.18) quickly transform into slow motions.

5 The steady-state solution of equation for slow motions

In what follows, we look for the steady-state solution of initial value problem for system of partial differential equations (4.8) where the heat supply is given in the form of a point source of a constant intensity. Accordingly, we take

[TABLE]

where $\bar{\chi}_{0}=\mathrm{const}$. Looking for the steady-state solution, we substitute $1(t)$ instead of $H(t)$ into Eq. (5.1), replace initial conditions (4.9) by the following boundary condition at infinity

[TABLE]

and drop out the terms involving the time derivatives in (4.8). This yields the following equation

[TABLE]

Now we apply the discrete Fourier transform with respect to the variables $m,\ n$ to Eq. (5.3):

[TABLE]

where subscript $F$ is the symbol of the discrete Fourier transform, $p_{1},\ p_{2}$ are the discrete Fourier transform parameters. Using the shift property for the discrete Fourier transform, one gets

[TABLE]

where

[TABLE]

and dot stands for the scalar product. Multiplying the both sides of the transformed equation (5.3) on \big{(}{\mathcal{L}}_{m\,n}\big{)}_{F}^{-1} results in555In previous paper gavrilov2018heat we inverse the corresponding difference operator directly in space domain (using the identities of the calculus of finite differences), not in the Fourier domain. the following equation:

[TABLE]

wherein

[TABLE]

Vector $\boldsymbol{\mathscr{C}}$ coincides with Kuzkin-Krivtsov-accepted the group velocity in the lattice under condition of zero dissipation $\eta=0$.

Thus, Eq. (5.10) describes the steady-state distribution for the Fourier images of the quantities $\theta_{m,n}$ defined by Eq. (4.7) (in particular, the kinetic temperature $T\equiv\theta_{0,0}$), caused by a point heat source of a constant intensity. According to conditions (5.2), Eq. (5.10) should be supplemented with the following boundary condition:

[TABLE]

Equation (5.10) can be rewritten as

[TABLE]

Here and in what follows,

[TABLE]

for any non-zero vector $\mathbf{a}$,

[TABLE]

for any unit vector $\breve{\mathbf{a}}$, $\times$ is the cross product. Delta-function in the right-hand side of (5.14) can be represented as

[TABLE]

where $\breve{\mathbf{n}}$ is an arbitrary unit vector such that

[TABLE]

We introduce angular variable $\alpha=\arctan\frac{x_{2}}{x_{1}}$ such that

[TABLE]

Due to symmetry, without loss of generality we can assume that (see Figure 1)

[TABLE]

The solution of (5.14) satisfying (5.13) is

[TABLE]

Here we have used Eq. (5.17) wherein $\breve{\mathbf{n}}=\breve{\boldsymbol{\mathscr{C}}}$, and the following formulas:

[TABLE]

Applying the inverse Fourier transform, one get the expression for the kinetic temperature $T=\theta_{0,0}$:

[TABLE]

Since the integrand in the right-hand side of Eq. (5.24) contains the Dirac delta-function $\delta_{1}$ in one-dimensional space, formula (5.24) can be rewritten in a simpler form involving a single integral (see Appendix A):

[TABLE]

where {\mathscr{C}_{0}}\big{(}p_{1},p_{2}^{(j)}\big{)} $(j=1,2)$ are defined by Eqs. (A.13), (A.14). If necessary, the far-field asymptotics of the solution (5.25) can be calculated by means of the Laplace method Fedoruk-Saddle .

5.1 The case $\alpha=0$

In the case $\alpha=0$, i.e. for two orthogonal rows of the lattice that contain the heat source located at $\mathbf{x}=\mathbf{0}$ (see Figure 1), it is possible to obtain the solution in a simpler form. This is done in Appendix B. The final result is $T=T_{1}+T_{2}$, where $T_{1}$ and $T_{2}$ are defined by Eqs.  (B.7), (B.9), respectively.

5.2 The case $\alpha=\pi/4$

Consider the particular case $\alpha=\pi/4$, i.e. the main diagonals of the lattice with respect to a point heat source position (see Figure 1). One has

[TABLE]

Thus, the integrands in the right-hand side of (5.25) have non-integrable singularities at $p_{1}=\pm\pi/2$, and therefore $T=+\infty$. In presence of the dissipation, this seems to be a quite unexpected result, which is discussed in Section 6.

5.3 The solution for a row $x_{2}/a=j$

Consider now the distribution of the kinetic temperature over a row of the lattice $x_{2}/a=j\in\mathbb{Z}$. Provided that assumption (5.20) is true, the analytic expression for the stationary distribution of the kinetic temperature $T$ is (5.25), where

[TABLE]

For $\alpha>\pi/4$, due to symmetry with respect to the main diagonal, one has

[TABLE]

where

[TABLE]

6 Numerics

In this section, we present the results of the numerical solution of the system of stochastic differential equations (2.1)–(2.3) with initial conditions (2.8). It is useful to rewrite Eqs. (2.1)–(2.3) in the dimensionless form

[TABLE]

where

[TABLE]

We consider the square lattice of $(2N+1)^{2}$ particles with the following boundary conditions

[TABLE]

where $i=\overline{-N,N}$. Actually, the specific form of this boundary conditions is not very important in our calculations, since we take large enough $N$ such that the wave reflections from the boundaries do not occur.666In Kuzkin-Krivtsov-accepted it is shown that for the non-dissipative system the continual non-stationary solution caused by a point pulse source is non-zero only in the circle with radius $ct$. It may be shown that the analogous result is true also for the dissipative system. To obtain a numerical solution in the case of the point source of the heat supply located at $i=0$, $j=0$, we assume that $\tilde{b}_{i,j}\rho_{i,j}=\delta_{i;0}\delta_{j;0}\tilde{b}\rho_{i,j}$ and use the scheme

[TABLE]

where $i,j=\overline{-N,N}$. Here the symbols with superscript $k$ denote the corresponding quantities at $\tilde{t}=\tilde{t}^{k}\stackrel{{\scriptstyle\mbox{\scriptsize def}}}{{=}}k\Delta\tilde{t}$: ${\tilde{u}}_{i,j}^{k}={\tilde{u}}_{i,j}(\tilde{t}^{k}),\ \tilde{v}_{i,j}^{k}=\tilde{v}_{i,j}(\tilde{t}^{k})$; $\rho^{k}$ are normal random numbers that satisfy (2.7) generated for all $\tilde{t}^{k}$. Without loss of generality we can take $\tilde{b}=1$. Due to symmetry of the problem we can perform the calculations only for $1/8$ of the whole lattice ($i=\overline{0,N}$, $j=\overline{0,i}$).

We perform a series of $r=1\dots R$ realizations of these calculations (with various independent $\rho^{k}_{(r)}$) and get the corresponding particle velocities $\tilde{v}^{k}_{(r)i,j}$. In accordance with (4.5), in order to obtain the dimensionless kinetic temperature

[TABLE]

we should average the doubled dimensionless kinetic energies:

[TABLE]

Numerical results (6.6) for the kinetic temperature can be compared with the analytical steady-state solution (5.25), expressed in the dimensionless form wherein $a=1,\ c=1$, $x_{1}=i$, $x_{2}=j$, and

[TABLE]

Note that the factor $1/2$ in the right-hand sides of the last formula appears according to Eq. (4.6).

The comparison between analytical and numerical results is presented in Figures 2–6. All calculations were performed for the following values of the problem parameters: $\tilde{\eta}=0.1$, $N=100$, $\tilde{t}=100$, $R=100000$. The time step is $\Delta{\tilde{t}}=0.025$. To perform the numerical calculations we use SciPy software scipy .

In Figure 2 one can see the general view of a central zone of the kinetic temperature distribution pattern in the lattice. One can observe that most of the heat propagates along the main diagonals of the lattice (with respect to a point source position). This result is in a qualitative agreement with results obtained earlier in studies mielke2006macroscopic ; harris2008energy , where energy transport in lattices are discussed in the deterministic case, and also with results of Kuzkin-Krivtsov-accepted .

In Figure 3 we compare the steady-state analytical solution in the form of Eqs. (B.7), (B.9) for the row $x_{2}=0$ (the blue solid line) and the numerical solution (the red crosses).

In Figure 4 we compare the steady-state analytical solution in the form of Eqs. (5.25), (5.29) for the row $x_{2}=7$ (the blue solid line) and the numerical solution (the red crosses). The corresponding view of a central zone of the kinetic temperature distribution is given in Figure 5. One can see that numerical value at the point $|x_{1}|=7$ of the lattice at the main diagonal is finite, whereas the analytic continuum solution is singular at this point.

Finally, in Figure 6 we compare the steady-state analytical solution in the form of Eqs. (5.25), (5.29) for the row $x_{2}=1$ (the blue solid line) and the numerical solution (the red crosses) in a central zone near the source.

One can see that in all cases, the analytical and numerical solutions are in a very good agreement everywhere except the main diagonals of the lattice $|x_{1}|=|x_{2}|$ (see Figures 3–6). The continuous analytical solution predicts singularities in the stationary solution at the main diagonals (see Section 5.2). Usually, such a result means that a corresponding non-stationary solution grows with time, and a steady-state solution at singular points does not exist. However, the numerical calculations based on the original infinite system of stochastic ODE contradict with such a hypothesis and predict that despite most of the heat propagates along the main diagonals of the lattice (see Figure 2), the kinetic temperature at the main diagonals apparently converges to finite values (see Figure 7, where plots of the numerical solution for the kinetic temperature at several fixed positions versus the time are given). In particular, this is true for the point $x_{1}=x_{2}=0$ where the heat source is located (see Figure 7). Note that in one-dimensional case considered in gavrilov2018heat we also observe the singularity in the continuous solution at the point where heat source is applied, whereas the numerics predicts the finite value of the kinetic temperature at that point777This is not discussed in gavrilov2018heat . In our opinion, such paradoxical result is caused by the continualization procedure (in particular, by the choice of heat supply intensity in singular form (5.1)), and it definitely needs an additional investigation. Note that in the conservative non-stochastic case the singularities at main diagonals were discovered in mielke2006macroscopic ; harris2008energy ; giannoulis2006continuum , where they are associated with non-smoothness of the dispersion relation.

7 Conclusion

In the paper, we started with equations (2.1) for stochastic dynamics of a two-dimensional harmonic square scalar lattice in a viscous environment. We introduced in the standard way the kinetic temperature in the lattice as a quantity proportional to the variance of the particle velocities. The most important results of the paper are the differential-difference equation (4.8) describing non-stationary heat propagation in the lattice and the analytical formula in the integral form (5.25) describing the steady-state kinetic temperature distribution in the lattice caused by a point heat source of a constant intensity.

The comparison between numerical solution of equations (2.1) and analytic steady-state solution (5.25) of differential-difference equation (4.8) demonstrates a a very good agreement everywhere except the main diagonals of the lattice with respect to the point source position (see Figures 3–6). The continuous analytical solution predicts singularities in the stationary solution at the main diagonals (see Section 5.2), whereas the numerical solution predicts that despite most of the heat propagates along the main diagonals of the lattice (see Figure 2), the kinetic temperature at the main diagonals apparently converges to finite values (see Figure 7). In our opinion, such a paradoxical result can be caused by the continualization procedure (in particular, by the choice of heat supply intensity in singular form (5.1)). Another possible oversimplification in physical modeling that makes the main diagonals to be preferable directions for heat propagation, apparently, is the choice of the potential of interaction involving for a given particle only four neighbours (as it is stated by Eqs. (2.5), (2.6)). This situation needs an additional investigation.

We expect that the results obtained in the paper can be used to describe the heat transfer in low-dimensional nanostructures and ultra-pure materials chang2008breakdown ; xu2014length ; goldstein2007mechanics . The subject of our future work is to generalize the theoretical results to the case of the graphene lattice and to verify the results by means of experiments with the laser heating of graphene hwang2016measuring ; indeitsev2017two .

Acknowledgements

The authors are grateful to V.A. Kuzkin, A.S. Murachev, and E.V. Shishkina for useful and stimulating discussions.

Appendix A Calculation of the steady-state solution in the integral form

To calculate the right-hand side of Eq. (5.24), one needs to use the formula (see G-Sh-1 )

[TABLE]

where $y_{j}$ are the roots of $f(y)$ lying inside the interval $I$.

One has

[TABLE]

[TABLE]

where

[TABLE]

such that $p_{2}^{(j,k)}\in[-\pi,\pi]$. The typical structure of roots $p_{2}^{(j,k)}$ in the case $0<\alpha<\pi/4$ is presented in Figure 8.

For $p_{1}\neq 0$ there exist exactly two roots lying in the interval $[-\pi,\pi]$. The first one corresponds to the choice $j=1,\ k=0$, and the second one corresponds to the choice $j=2,\ k=0$ for $p_{1}>0$ and $j=2,\ k=-1$ for $p_{1}<0$. Thus, for $0<\alpha<\pi/4$ we put

[TABLE]

Applying now formula (A.1), one gets

[TABLE]

According to (A.5) one has:

[TABLE]

Using (5.12), and taking into account identities

[TABLE]

one obtains

[TABLE]

According to (3.4), (5.11), one gets

[TABLE]

where $j=1,2$. Finally, substituting of these expressions into Eq. (A.8) and simplifying of expression obtained, results in the formula for the solution in the solution in the integral form:

[TABLE]

Appendix B Calculation of the steady-state solution in the particular case

$\alpha=0$

In the case under consideration, due to (A.5) one has

[TABLE]

Since in the case $p=p_{2}$ the roots lie at the boundaries of the interval $[-\pi,\pi]$, we take the corresponding contributions multiplied by $1/2$. Since these contributions are equal to each other, we can put

[TABLE]

and use the formulas obtained above. According to Eqs. (A.13), (A.14) one gets

[TABLE]

Thus,

[TABLE]

Let $\gamma=\frac{1}{\cos y}$. One has

[TABLE]

and, therefore,

[TABLE]

where

[TABLE]

is the exponential integral abramowitz1972handbook . For $T_{2}$ one gets the solution in the integral form

[TABLE]