Heat transport via low-dimensional systems with broken time-reversal symmetry
HHeat transport via low-dimensional systems with broken time-reversal symmetry
Shuji Tamaki, Makiko Sasada, and Keiji Saito Department of Physics, Keio University, Yokohama 223-8522, Japan Graduate School of Mathematical Sciences, The University of Tokyo, Komaba, Tokyo 153-8914, Japan (Dated: September 19, 2017)We consider heat transport via systems with broken time-reversal symmetry. We apply magneticfields to the one-dimensional charged particle systems with transverse motions. The standard mo-mentum conservation is not satisfied. To focus on this effect clearly, we introduce a solvable model.We exactly demonstrate that the anomalous transport with a new exponent can appear. We numer-ically show the violation of the standard relation between the power-law decay in the equilibriumcorrelation and the diverging exponent of the thermal conductivity in the open system.
PACS numbers: 05.40.-a, 02.50.-r,63.22.+m, 44.10.+i
Introduction.—
It is generally believed that heat con-duction in low-dimensional nonlinear systems is anoma-lous from many theoretical [1–23] and experimental stud-ies [24–26]. In a one-dimensional system of N particlesconnected at the ends to heat baths with a small temper-ature difference ∆ T , the thermal conductivity is definedas κ = JN/ ∆ T , where J is the steady state current persite. The anomalous heat transport is given by the di-vergence of κ with increasing system size: κ ∼ N α , < α < . (1)The anomalous behavior is related to the equilibrium cur-rent correlation with slow decay in a closed system: C ( t ) = N − (cid:104) J tot ( t ) J tot (cid:105) eq ∼ t − β , < β < , (2)where J tot is the total energy current and (cid:104) ... (cid:105) eq is theequilibrium average. A slow decay leads to the divergingthermal conductivity through the Green-Kubo formula.Generally, in nonlinear chains, there are several con-served quantities in the periodic boundary condition, i.e.,energy, momentum, and the so-called stretch variables[15, 16]. These conserved quantities are key-ingredientsin understanding the anomalous behavior. Recently,there has been significant progress in theories on the equi-librium current correlation by considering the conservedquantities. This remarkable progress included finding anexactly solvable model with anomalous behavior, which isnow called the momentum exchange (ME) model [19, 20].This model contains hybrid dynamics of deterministic dy-namics and stochastic “conservative” noise, which con-serves the three variables. Exact analysis of the currentcorrelation function shows β = 1 / β = 2 / / H = N (cid:88) i =1 | p i − e i A ( q i ) | / V ( r i ) , (3)where we set the masses to unity. The vector q i specifiesthe position of the i th particle, r i is the stretch vector de-fined below in (5), and V is the spring potential betweenthe nearest neighbor sites. The variables p i and A ( q i )are, respectively, the canonical momentum and the gaugepotential, and e i is the charge of the i th particle. Theactual velocity is given by v i ≡ ˙ q i = p i − e i A ( q i ). Weconsider the static magnetic field B , and then the dy-namics are given by the Lorentz force and spring forces˙ v i = e i v i × B − ∂ q i [ V ( r i − ) + V ( r i )] . (4)From these dynamics, one finds that each summation ofthe following local variables is conserved: r i = q i +1 − q i , (5) (cid:15) i = | v i | / V ( r i ) + V ( r i − )] / , (6) P i = v i − e i q i × B = p i − e i A ( q i ) − e i q i × B , (7) a r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p xq x r x q x r x q x r x q x q N − x r N − x q Nx r Nx q N +1 x r N +1 x q N +2 xN P i =1 r i x BBBBB
FIG. 1: Schematic picture of the periodic boundary condition.The x -components of variables are shown. where r i and (cid:15) i are, respectively, the local stretch andenergy variables. The variable P i is a pseudomomentum [28] which is not equivalent to the canonical momentum p i [29]. Hence, the standard momentum conservation isreplaced by the conservation of this variable. From thismodification, the dynamics should be newly categorizedin the context of heat conduction and careful analysis onthe exponent β is required.We note here that for nonlinear systems, it is generallydifficult to obtain accurate values of the exponent even inlarge-scale numerical calculations. Hence, we introducea solvable model by extending the ME model to the caseof finite magnetic fields. Then, we clearly argue that themagnetic fields can generate a new exponent. Velocity exchange models.—
An exactly solvable modelthat we introduce is a harmonic chain with the potential V ( r i ) = | r i | /
2, where the time evolution is composed ofthe deterministic dynamics (4) and conservative noisesthat conserve each summation of (5)-(7). The change ofvariables from time t to t + dt is given by dr ia = ( v i +1 a − v ia ) dt , (8) dv ix = ( r ix − r i − x + e i Bv iy ) dt + dn ix ( v i +1 x − v ix ) + dn i − x ( v i − x − v ix ) , (9) dv iy = ( r iy − r i − y − e i Bv ix ) dt + dn iy ( v i +1 y − v iy ) + dn i − y ( v i − y − v iy ) , (10)where a = x, y . The magnetic field B is applied in the z -direction and we consider only motions of particles inthe xy -plane, which are relevant to the magnetic field.The vector v i = ( v ix , v iy ) is the velocity vector of the i thparticle and r i = ( r ix , r iy ) is the stretch vector definedin Eq.(5). We consider the periodic boundary conditionimposing r i + N = r i and v i + N = v i with an even num-ber N (See Fig.1). The variable dn ia takes the value 0or 1 with the Poisson process satisfying the noise average (cid:104)(cid:104) dn ia (cid:105)(cid:105) = γdt . Hence, the noises stochastically exchangevelocities between the nearest neighbor sites. One caneasily check that each summation of the variables (5)-(7)is conserved. When we switch off the magnetic field, thedynamics for variables of x and y components indepen-dently follow the original ME dynamics.We consider two cases: case (I) with uniform charge e i = 1 and case (II) with alternate charge e i = ( − i .By employing the deterministic dynamics only, one can − −
100 0 100 200 i C (cid:15)(cid:15) ( i , t ) t = 50 t = 100 t = 150 t = 200 C (cid:15)(cid:15) ( i , t ) t = 50 t = 100 t = 150 t = 200case (II)case (I) FIG. 2: Numerical calculation of the space-time correlation C (cid:15)(cid:15) ( i, t ) = (cid:104) δ(cid:15) i +1 ( t ) (cid:15) (cid:105) eq for the dynamics (8)-(10). The fifth-order Runge-Kutta algorithm with dt = 0 .
001 is used for thedeterministic dynamics, and 10 initial states were taken fromthe canonical distribution. Parameters: B = 1, N = 512, T = 1, and γ = 0 .
1. One can clearly see the absence of soundwaves in case (I), while case (II) has finite ballistic peaks. derive the dispersion relation for each case [30] ω I ( k ) = (cid:112) (2 sin( k/ + ( B/ ± B/ , (11) ω II ( k ) = (cid:113) B / ± (cid:112) (2 + B / − (2 sin k ) , (12)where the subscripts I and II indicate the two casesand k is the wave number. From these expressions, thesound velocities are calculated using dω ( k ) /dk (cid:12)(cid:12) k → . Case(I) has zero sound velocity while case (II) has a finitevalue of the velocity. We numerically check this by con-sidering the space-time correlation of the local energy C (cid:15)(cid:15) ( i, t ) := (cid:104) δ(cid:15) i +1 ( t ) δ(cid:15) (0) (cid:105) eq where (cid:15) i is defined in Eq.(6)and δ(cid:15) = (cid:15) − (cid:104) (cid:15) (cid:105) eq . The symbol (cid:104) ... (cid:105) eq is the average overthe canonical ensemble (cid:81) j,a e − ( r j,a + v ja ) / T /Z with tem-perature T and the normalization Z . Here, the Boltz-mann constant is set to unity. In Fig.2, we present nu-merical results for the system size N = 512 with T = 1and B = 1. The figure clearly shows the absence ofsound waves in case (I), while case (II) has finite soundpropagation indicated by ballistic peaks. Thus, cases (I)and (II) have contrasting differences in the dynamics,and hence, we discuss heat conduction with broken time-reversal symmetry, comparing these cases. Methods and main results of equilibrium correlation.—
For zero magnetic field, the exponent β = 1 / (cid:15) i ( t ) − (cid:15) i (0) = − I i [0 , t ]+ I i − [0 , t ], where I i [0 , t ] is the accumulation upto time t of the energy current measured between the i .
01 0 . − − − − − t C ( t ) B = 0case ( I )case (II) ∝ t − / ∝ t − / τ R τ d t C ( t ) B =0case ( I )case (II) ∝ τ / ∝ τ / FIG. 3: Numerical check of the exponent β . Parameters: N = 2048, T = 1, γ = 0 . B = 1 for cases (I) and (II).Numerical procedure is the same as in Fig.1. and ( i + 1)th sites: I i [0 , t ] = (cid:90) t ds (cid:0) J di ( s ) + J si ( s ) (cid:1) + (cid:90) t dJ mi ( s ) , (13) J di ( s ) = − (cid:88) a = x,y r ia ( s ) ( v i +1 a ( s ) + v ia ( s )) / , (14) J si ( s ) = − (cid:88) a = x,y γ ( v i +1 a ( s ) − v ia ( s )) / , (15) dJ mi ( s ) = − (cid:88) a = x,y ( v i +1 a ( s ) / − v ia ( s ) / dm ia ( s ) , (16)where dm ia is the Martingale noise defined as dm ia = dn ia − γdt [31], J d and J s are the instantaneous currentsfrom the deterministic dynamics and average stochasticnoise, respectively. The third current dJ m is a currentfrom the Martingale noise, whose contribution to thethermal conductivity is constant and the correlations be-tween dJ m and J d,s vanish [19, 20]. Since dJ m does not generate power law behavior in the current correlation,we consider only the contribution of J d and J s as C ( t ) ≡ N − (cid:104)(cid:104) J d tot ( t ) J d tot (cid:105)(cid:105) eq = (cid:104)(cid:104) J d tot ( t ) J dc (cid:105)(cid:105) eq , (17)where J d tot = (cid:80) i J di + J si and J dc = ( J dN + J d ) /
2. We used (cid:80) i J si = 0, and the symbol (cid:104)(cid:104) ... (cid:105)(cid:105) eq denotes the averageover the canonical ensemble as well as the average overnoises. We follow the technique developed in Refs.[19,20]. We consider the Laplace transform C ( λ ) = (cid:90) ∞ dte − λt C ( t ) = (cid:90) ∞ dt (cid:104) ( e − ( λ − L ) t J d tot ) J dc (cid:105) eq = (cid:104) (cid:2) ( λ − L ) − J d tot (cid:3) J dc (cid:105) eq , (18)where the operator L is the time evolution operator givenby L = L + γ S , where L and S respectively correspondto the deterministic dynamics and conservative noises: L = (cid:88) i,a ( v i +1 a − v ia ) ∂∂r ia + ( r ia − r i − a ) ∂∂v ia + (cid:88) i e i B (cid:20) v iy ∂∂v ix − v ix ∂∂v iy (cid:21) , (19) S f ( r , v ) = (cid:88) i (cid:88) a f ( r , v i | i +1 ,a ) − f ( r , v ) . (20)Here, the function f is an arbitrary function of r and v and v i | i +1 ,a is obtained from v by exchanging the vari-ables v i a and v i +1 a .The details to derive the function C ( λ ) are providedin the supplementary material [30], and below we discussphysically crucial results. The Laplace transforms in thethermodynamic limit are given as follows: C I ( λ ) = T π (cid:90) π dk cos ( k/
2) ( λ + 4 γ sin ( k/ λ + 8( λγ + 2) sin ( k/ λ + 4 γ sin ( k/ ( λ + 8( λγ + 2) sin ( k/ B λ ( λ + 8 γ sin ( k/ , (21) C II ( λ ) = T π (cid:90) π dk cos ( k/ µ (cid:0) ( B + µ ) a − B (cid:1) + 2 (cid:0) B (2 µ − γa ) − γµ a (cid:1) cos k + a (cid:0) µ cos k + 8 γ cos k (cid:1) ( B + µ ) (( B + µ ) a − B ) + 8 ( B γ a + µ ( a − k − γ a cos k , (22)where the subscripts indicate two cases, and µ = λ + 2 γ , a = 8 − γ + µ , a = − a + 4 and a = 4 + 4 γ − γ (8 + µ ). The asymptotic real-time representation isanalyzed by the inverse Laplace transform, considering asmall wave number for finite B and γ , and one gets C I ( t ) ∼ A t − / + A t − / cos( Bt ) + A t − / , (23) C II ( t ) ∼ A t − / , (24)where A , , , are constant values which depend on B . We now list physically crucial observations for these ex-act results. In both cases, power-law behavior exists.Eq.(23) includes the power law term with oscillation intime which rapidly decays for finite B , and most impor-tantly the new exponent β = 3 / β = 1 /
2, which is the same exponent asfor B = 0. This implies that the universality class de-pends on the charge structure of the system. These exactfindings are the main results in this paper. A numeri-cal evaluation of these observations is presented in Fig.3.Rapid decay with power law behavior are observed forany case. The numerical calculation accurately repro-duces the known exponent β = 1 / β = 1 /
2, and case (I) has β = 3 / Numerical results of the exponent α .— We next con-sider the exponent α in Eq. (1) that is measured inthe nonequilibrium steady state when the system is con-nected to thermal reservoirs. We use a numerical ap-proach here. We attach the Nose-Hoover thermostat tothe end particles [32]. The dynamics for the sites from i = 2 to N − dv (cid:96)a = [ r (cid:96)a − r (cid:96) − a + e (cid:96) B ( δ a,x v (cid:96)y − δ a,y v (cid:96)x )] dt + δ (cid:96), dn a ( v a − v a ) + δ (cid:96),N dn N − a ( v N − a − v Na ) − ξ (cid:96)a v (cid:96)a dt , (25) dξ (cid:96)a = γ (cid:48) ( v (cid:96)a /T (cid:96) − dt , (26)where (cid:96) = 1 or N . T and T N are the reservoir’s tem-peratures at the first and the N th particles, respectively.We show the system-size dependence of the thermal con-ductivity up to N = 8192 in Fig.4. Numerical error wassmaller than the size of the points. The system size issufficiently large to obtain the asymptotic behavior of thepower law divergence. In the figure, the case with zeromagnetic field and case (II) show α = 1 /
2, while the ex-ponent in case (I) is neither 1 / /
3. The best fit is0 . ± . α and β . Toour knowledge, a rigorous derivation of the relationshipbetween α and β has never been made. Thus far, there isonly a phenomenological interpretation of the case whenthe system has finite sound velocity. The argument isbased on the modified Green-Kubo formula κ ∼ (cid:90) τ N dtC ( t ) . (27)When the system has a finite sound velocity, one phe-nomenologically uses τ N ∼ N/c , where c is the soundvelocity, and obtains the relation α = − β + 1. Althoughit is not derived rigorously, thus far, it seems to workwell. In fact, the case of zero magnetic field and case (II)follow this relation. However, in case (I), where no soundwave exists, this relation is not applicable anymore. Nu-merical results indicate τ N ∼ N ν with ν ∼ . ± . Discussion.—
In this paper, we studied the heat trans-fer in one-dimensional systems with broken time-reversal ∝ N / ∝ N / ∝ N . ± . ∝ N / N κ B = 0case ( I )case (II) FIG. 4: Numerical calculation of the system-size dependenceof the thermal conductivity. The fifth-order Runge-Kutta al-gorithm with dt = 0 .
001 is used and the steady state waschecked from the uniformity of local energy current. Param-eters: ( T , T N ) = (2 , γ = 0 . γ (cid:48) = 1, and B = 1 forcases (I) and (II). The exponent α in case (I) is different fromknown exponents. The error bar of the exponent for case (I)is estimated using Gnuplot for the range N = [512 , symmetry for the first time. We considered systems withvery weak charges under a strong magnetic field so thatthe dynamics are dominated by the Lorentz force as wellas the spring forces connecting particles. To clarify theargument on the exponent, we introduced an exactlysolvable model in the spirit of the ME model. Basedon this model, we found that a new power-law decay ex-ponent can appear. We will report elsewhere on severalother results including the effects of higher dimensions[33].In systems without time-reversal symmetry, the stan-dard fluctuating hydrodynamic theory is not applicable,as the Euler equations for conserved quantities are notclosed due to the expression of the pseudomomentum.Physically, the magnetic field induces cyclotron motionand hence, the particles tend to be localized. Based onthis, one might think that the conservation of pseudo-momentum is irrelevant to macroscopic behavior and thesystem may exhibit diffusive heat conduction. We notethat a recent non-acoustic model with momentum con-servation shows diffusive transport [34]. The same con-clusion had been speculated based on the mode-couplingargument for nonlinear systems with zero sound velocityin Ref.[14]. However, our case showed that anomalousheat conduction robustly exists and the new power-lawdecay exponent can appear. In order to understand thesenontrivial results, a precise description of the hydrody-namics is required.The present model clearly shows that the absence ofsound waves cause the violation of the usual relationship α = − β +1 that is satisfied for systems with sound waves.An analytical derivation of α based on the equilibriumcorrelation β = 3 / [1] Thermal Transport in Low Dimensions: From StatisticalPhysics to Nanoscale Heat Transfer , Edited by S. Lepri(Springer, 2016).[2] A. Dhar,
Heat transport in low-dimensional systems , Adv.Phys. , 457 (2008).[3] S. Lepri, R. Livi, and A. Politi, Thermal conductionin classical low-dimensional lattices , Phys. Rep. , 1(2003).[4] S. Lepri, R. Livi, and A. Politi,
Heat conduction in chainsof nonlinear oscillators , Phys. Rev. Lett. , 1896 (1997).[5] P. Grassberger, W. Nadler, and L. Yang, Heat conductionand entropy production in a one-dimensional hard-particlegas , Phys. Rev. Lett. , 180601 (2002).[6] G. Casati and T. Prosen, Anomalous heat conduction ina one-dimensional ideal gas , Phys. Rev. E , 015203(2003).[7] T. Mai, A. Dhar and O. Narayan, Equilibration and uni-versal heat conduction in Fermi-Pasta-Ulam chains , Phys.Rev. Lett. , 184301 (2007).[8] K. Saito and A. Dhar, Heat conduction in a three dimen-sional anharmonic crystal , Phys. Rev. Lett. , 040601(2010).[9] L. Delfini, S. Lepri, R. Livi, and A. Politi,
Self-consistentmode-coupling approach to one-dimensional heat trans-port , Phys. Rev. E , 060201 (2006).[10] A. Pereverzev, Fermi-Pasta-Ulam β lattice: Peierls equa-tion and anomalous heat conductivity , Phys. Rev. E Anomalous heat conduc-tion in one-dimensional momentum-conserving systems ,Phys. Rev. Lett. , 200601 (2002).[12] J. S. Wang and B. Li, Intriguing heat conduction ofa chain with transverse motions , Phys. Rev. Lett. ,074302 (2004).[13] J. Lukkarinen and H. Spohn, Anomalous energy transportin the FPU- β chain , Commun. Pure Appl. Math. Thermalconductivity and bulk viscosity in quartic oscillator chains ,Phys. Rev. E , 031202 (2005).[15] H. van Beijeren, Exact results for anomalous transportin one-dimensional Hamiltonian systems , Phys. Rev. Lett. , 180601 (2012).[16] H. Spohn,
Nonlinear fluctuating hydrodynamics for an-harmonic chains , J. Stat. Phys. , 1191 (2014).[17] C. B. Mendl and H. Spohn,
Dynamic correlators ofFermi-Pasta-Ulam chains and nonlinear fluctuating hy-drodynamics , Phys. Rev. Lett. , 230601 (2013).[18] C. B. Mendl and H. Spohn,
Equilibrium time-correlationfunctions for one-dimensional hard-point systems , Phys.Rev. E , 012147 (2014).[19] G. Basile, C. Bernardin, and S. Olla, Momentum con-serving model with anomalous thermal conductivity in lowdimensional systems , Phys. Rev. Lett., Thermal conductiv-ity for a momentum conservative model , Comm. in Math.Phys.
67 (2009).[21] M. Jara, T. Komorowski and S. Olla,
Superdiffusion ofEnergy in a Chain of Harmonic Oscillators with Noise ,Comm. in Math. Phys. , 407 (2015).[22] S. Lepri, C. Mej´ıa-Monasterio, and A. Politi,
A stochas-tic model of anomalous heat transport: analytical solutionof the steady state , J. Phys. A: Math. Theor. , 025001(2009).[23] L. Delfini, S. Lepri, R. Livi, and A. Politi, NonequilibriumInvariant Measure under Heat Flow , Phys. Rev. Lett. ,120604 (2008).[24] C. W. Chang, D. Okawa, H. Garcia, A. Majumdar, andA. Zettl,
Breakdown of Fourier’s law in nanotube thermalconductors , Phys. Rev. Lett. , 075903 (2008).[25] V. Lee, C. H. Wu, Z. X. Lou, W. L. Lee, and C. W.Chang,
Divergent and ultrahigh thermal conductivity inmillimeter-long nanotubes , Phys. Rev. Lett. , 135901(2017).[26] X. Xu, et al.,
Length dependent thermal conductivity insuspended Graphene , Nat. Comm. , 3869 (2014).[27] A. V. Savin and L. I. Manevitch, Solitons in spiral poly-meric macromolecules , Phys. Rev. E , 7065 (2000).[28] B. R. Johnson, J. O. Hirschfelder and K. Yang, Interac-tion of atoms, molecules, and ions with constant electricand magnetic fields , Rev. Mod. Phys.
109 (1983).[29] One can see the difference by considering an explicitgauge potential, e.g., A ( q i ) = − ( q i × B ) / Stochastic Integration and DifferentialEquations , (2nd Ed. Springer, 2005).[32] J. M. Thijssen,
Computational Physics (Cambridge,2007).[33] K. Saito and M. Sasada, arXiv:1706.09668.[34] T. Komorowski and S. Olla,
Diffusive Propagation of En-ergy in a Non-acoustic Chain , Arch. Rational Mech. Anal. , 95 (2017).
Supplementary Material for“Heat transport via low-dimensional systems with broken time-reversal symmetry”
Shuji Tamaki , Makiko Sasada and Keiji Saito Department of Physics, Keio University, Yokohama 223-8522, Japan Graduate School of Mathematical Sciences, The University of Tokyo, Komaba, Tokyo 153-8914, Japan
NOTATIONS
We here fix several notations to calculate the equilibrium correlation. Let q ia be the a -component of the positionvector of the i th particle q i , and r ia be r ia = q i +1 a − q ia . We impose the periodic boundary condition: r i + N = r i and v i + N = v i . (S.1)We note that for a given initial state, the following quantity is conserved.¯ r a := (1 /N ) N (cid:88) i =1 r ia . (S.2)This conservation law follows the time-evolution ˙ r ia = v i +1 a − v ia . See Fig.5, which schematically depicts thissituation. We define the following new variables ζ ia := q ia − ( i − r a . (S.3)These variables satisfy the following relations: ζ i + N a = ζ ia (S.4) r ia = q i +1 a − q ia = ζ i +1 a − ζ ia + ¯ r a . (S.5) DISPERSION RELATION
The dispersion relation is derived only from the deterministic dynamics. For case (I), the dynamics is given by¨ ζ jx − [∆ ζ ] jx − B ˙ ζ jy = 0 , (S.6)¨ ζ jy − [∆ ζ ] jy + B ˙ ζ jx = 0 , (S.7)where [∆ ζ ] ja = ζ j +1 a + ζ j − a − ζ j a . By the Fourier transform, ζ ja ( t ) = (cid:80) k ˜ ζ ak e − i ( kj − ωt ) /N , where k is the wavenumber k = 2 π/N, π/N, · · · , π ( N − /N, π . Then, we have the equation A ˜ ζ k = , where ˜ ζ k = ( ζ xk , ζ yk ) T , andthe matrix A is given by A = (cid:18) − ω + (2 sin( k/ , − iBωiBω, − ω + (2 sin( k/ (cid:19) . (S.8)The dispersion relation is given by the condition det A = 0: ω I ( k ) = (cid:112) (2 sin( k/ + ( B/ ± B/ . (S.9)The dispersion relation for case (II) is similarly given. The deterministic part in the dynamics is as follows:¨ ζ jx − [∆ ζ ] jx − ( − j B ˙ ζ jy = 0 , (S.10)¨ ζ jy − [∆ ζ ] jy + ( − j B ˙ ζ jx = 0 . (S.11)We define the Fourier transform for even and odd sites as ζ ja = (cid:26) (1 /N ) (cid:80) k ˜ ζ ea,k e − i ( kj − ωt ) j ≡ /N ) (cid:80) k ˜ ζ oa,k e − i ( kj − ωt ) j ≡ . (S.12) xq x r x q x r x q x r x q x q N − x r N − x q Nx r Nx q N +1 x r N +1 x q N +2 x r N +2 x q N +3 x q N − x r N − x q Nx r Nx q N +1 x r N +1 x q N +2 x N ¯ r x N ¯ r x FIG. 5: Schematic picture of the system structure with the periodic boundary condition. The x -components of the variablesare shown. Then, we have A ˜ ζ k = , where ˜ ζ k = (˜ ζ ex,k , ˜ ζ ey,k , ˜ ζ ox,k , ˜ ζ oy,k ) T and A = − ω + 2 , − Biω, − k, Biω, − ω + 2 , , − k − k, , − ω + 2 , Biω , − k, − Biω, − ω + 2 . (S.13)From det A = 0, the following dispersion relation is obtained ω II ( k ) = (cid:113) B / ± (cid:112) (2 + B / − (2 sin k ) . (S.14) CALCULATION OF u λ AND C ( λ ) The correlation function C ( t ) is given by C ( t ) = (cid:104)(cid:104) J d tot ( t ) J dc (cid:105)(cid:105) eq = (cid:104)(cid:104) (cid:88) j,a ( − / v ja ( t ) + v j +1 a ( t )) r ja ( t ) J dc (cid:105)(cid:105) eq , (S.15)where (cid:104)(cid:104) ... (cid:105)(cid:105) eq implies the average over canonical measure and noise average. Based on the relation (S.5), we decomposethe correlation function into two terms C ( t ) = C A ( t ) + C B ( t ) (S.16) C A ( t ) = (cid:88) j,a (cid:104)(cid:104) ( − / v ja ( t ) + v j +1 a ( t ))( ζ j +1 a ( t ) − ζ ja ( t )) J dc (cid:105)(cid:105) eq , (S.17) C B ( t ) = (cid:88) j,a (cid:104)(cid:104) ( − / v ja ( t ) + v j +1 a ( t ))¯ r a J dc (cid:105)(cid:105) eq . (S.18)The Laplace transform of the correlation function is then given by C ( λ ) = C A ( λ ) + C B ( λ ) . (S.19)As we can see below in the detailed calculations for cases (I) and (II), the main contribution is from C A ( λ ), while C B ( λ ) can be neglected. We note that the time evolution of the local current is driven by the systematic part and theMartingale part. The time evolution operator of the systematic part is denoted by L in the main text. In general, thedetailed expression of the Martingale part depends on the quantity that one considers. From the Martingale property,the Martingale term contributes only at t = 0 on average over the noises.We solve the functions u λ,A and u λ,B ( λ − L ) u λ,A := (cid:88) j,a ( − / v ja + v j +1 a )( ζ j +1 a − ζ ja )= (cid:88) i,j ( − / δ i − j, − δ i − j, − )( ζ ix v jx + ζ iy v jy ) , (S.20)( λ − L ) u λ,B := (cid:88) j,a ( − / v ja + v j +1 a )¯ r a = − (cid:88) j,a ¯ r a v ja . (S.21)Then, we obtain the Laplace transform C A ( λ ) = (cid:104) u λ,A J dc (cid:105) eq and C B ( λ ) = (cid:104) u λ,B J dc (cid:105) eq . (S.22)Here (cid:104) ... (cid:105) eq implies an average over the canonical average. Note that we have already taken the noise average usingthe Martingale property. Case (I)
We first explain case (I) following the method in Ref.[1]. Since the energy current J tot is a linear function withrespect to v , we can write the operator L in the following form: L = (cid:88) i,a v ia ∂∂ζ ia + ([∆ ζ ] ia + γ [∆ v ] ia ) ∂∂v ia + (cid:88) i B (cid:20) v iy ∂∂v ix − v ix ∂∂v iy (cid:21) , (S.23)where [∆ f ] := f i +1 a + f i − a − f ia . We search for the solution of u λ,A in the following form: u λ,A = (cid:88) i,j g ( i − j ) ζ ix ζ jy + g ( i − j )( ζ ix v jx + ζ iy v jy ) + g ( i − j )( ζ ix v jy − ζ iy v jx ) + g ( i − j ) v ix v jy , (S.24)where g j ( z ) has the following symmetries for an arbitrary integer zg j ( z ) = − g j ( − z ) and g j ( z ) = g j ( z + N ) . (S.25)We substitute (S.24) into Eq.(S.20). Through straightforward calculation, the following equations are obtained: λg ( z ) − g ( z ) = 0 , (S.26)( λ − γ ∆) g ( z ) + Bg ( z ) = h ( z ) , (S.27) g ( z ) + Bg ( z ) − ( λ − γ ∆) g ( z ) + ∆ g ( z ) = 0 , (S.28) − g ( z ) + ( λ − γ ∆) g ( z ) = 0 , (S.29)where h ( z ) = − (1 / δ z, − δ z, − ) and ∆ g j ( z ) = g j ( z + 1) + g j ( z − − g j ( z ). From these, one can obtain uniquesolutions of the functions g j . We consider the discrete Fourier transform to solve the equations g j ( z ) = (cid:88) k ˜ g j ( k ) e − ikz /N . (S.30)Then Eqs.(S.26)-(S.29) can be solved in Fourier space: A I ˜ g ( k ) = u (S.31)where ˜ g ( k ) and u are vectors given by ˜ g ( k ) = (˜ g ( k ) , ˜ g ( k ) , ˜ g ( k ) , ˜ g ( k )) T and u = (0 , − i sin k, , T , respectively.The matrix A I is given by A I = λ, , k/ , , λ + γ (2 sin( k/ , B, − , − B, λ + γ (2 sin( k/ , (2 sin( k/ , , − , λ + 2 γ (2 sin( k/ . (S.32)Note here that only the function g has a finite contribution to the Laplace transform C A ( λ ): C A ( λ ) = (cid:88) i,j (cid:88) a g ( i − j ) (cid:104) ζ ia v ja J dc (cid:105) eq = ( − T / (cid:88) i,a g ( i ) [ (cid:104) ( ζ i +1 a + ζ ia ) r a (cid:105) eq + (cid:104) ( ζ i +2 a + ζ i +1 a ) r a (cid:105) eq ]= ( − T / (cid:88) i,a g ( i ) [ (cid:104) ( ζ i +1 a + ζ ia − ζ − i +1 a − ζ − ia ) r a (cid:105) eq + (cid:104) ( ζ i +2 a + ζ i +1 a − ζ − i +2 a − ζ − i +1 a ) r a (cid:105) eq ]= ( − T / (cid:88) i,a g ( i ) (cid:104) ( i (cid:88) j = − i ( r j a − ¯ r a ) + i − (cid:88) j = − ( i − ( r j a − ¯ r a )) r a (cid:105) eq + (cid:104) ( i +1 (cid:88) j = − i +1 ( r j a − ¯ r a ) + i (cid:88) j = − i +2 ( r j a − ¯ r a )) r a (cid:105) eq = ( − T / (cid:88) i,a g ( i ) (cid:104) ( i (cid:88) j = − i ( r j a − ¯ r a ) + i − (cid:88) j = − ( i − ( r j a − ¯ r a )) r a (cid:105) eq , (S.33)where we used the symmetry (S.25) to write the expression solely in terms of the stretch variables and the translationalinvariance in the equilibrium correlation between stretches. We define the function F ia as F ia := ( T / (cid:104) ( i (cid:88) j = − i ( r j a − ¯ r a ) + i − (cid:88) j = − ( i − ( r j a − ¯ r a )) r a (cid:105) eq . (S.34)We note here that the following relation from the simple calculation is satisfied, regardless of a = x, y :[∆ F ] ia = T ( − δ i, + δ i, − ) . (S.35)From this, we have the Fourier transform for the function F ia ˜ F ( k ) = iT cos( k/ / sin( k/ . (S.36)Hence, we arrive at the expression for C A ( λ ) C A ( λ ) = − (1 /N ) (cid:88) k (cid:88) j g ( j ) e − ikj ˜ F ( k ) = − (1 /N ) (cid:88) k ˜ g ( − k ) ˜ F ( k ) → ( T /π ) (cid:90) π dk cos ( k/ N ( k ) / D ( k ) , (S.37) N ( k ) = ( λ + 4 γ sin ( k/ λ + 8( λγ + 2) sin ( k/ , (S.38) D ( k ) = ( λ + 4 γ sin ( k/ ( λ + 8( λγ + 2) sin ( k/ B λ ( λ + 8 γ sin ( k/ . (S.39)We consider the inverse Laplace transform for C A ( λ ) to get C A ( t ). C A ( t ) = (cid:90) c + i ∞ c − i ∞ dλ C A ( λ ) e λt / πi . (S.40)We note that the poles λ (cid:63) in the function D are given by λ (cid:63)σσ (cid:48) = − γ (2 sin( k/ + σ [ a ( k ) + σ (cid:48) a ( k )] / / √ , (S.41) a ( k ) = − B − k/ + γ (2 sin( k/ , (S.42) a ( k ) = (cid:2) − γ (2 sin( k/ + ( B + 4(2 sin( k/ + γ (2 sin( k/ ) (cid:3) / , (S.43)where σ, σ (cid:48) take values of ±
1. For finite B and γ , and for a small wave number k , the expansion for the poles can beobtained as λ (cid:63) ++ = − γk /B + · · · , (S.44) λ (cid:63) + − = iB − ( − i/B + γ ) k + · · · , (S.45) λ (cid:63) − + = − γk + · · · , (S.46) λ (cid:63) −− = − iB − (2 i/B + γ ) k + · · · , (S.47)(S.48)The poles λ (cid:63) −− and λ (cid:63) + − provide the oscillation damping t − / cos( Bt ), and the pole λ (cid:63) − + provides the decay t − / .The power law t − / is obtained from the pole λ (cid:63) ++ .We finally consider u λ,B ( λ ). From the simple calculation for Eq.(S.21), we can check that the following expressionis a solution of u λ,B ( λ ): u λ,B = − (cid:88) j (cid:2) v jx (¯ r x λ − ¯ r y B ) / ( λ + B ) + v jy (¯ r x B + ¯ r y λ ) / ( λ + B ) (cid:3) . (S.49)From this expression, C B ( λ ) is proportional to (cid:104) ¯ r x r x (cid:105) eq or (cid:104) ¯ r y r y (cid:105) eq , which is O (1 /N ), and hence it is negligible inthe thermodynamic limit.0 Case (II)
The calculation for case (II) is essentially the same as for case (I). Note here that the operator L is given by L = (cid:88) i,a v ia ∂∂ζ ia + ([∆ ζ ] ia + γ [∆ v ] ia ) ∂∂v ia + (cid:88) i ( − i B (cid:20) v iy ∂∂v ix − v ix ∂∂v iy (cid:21) . (S.50)We formulate an ansatz for u λ,A u λ,A = (cid:88) i ≡ j mod 2 (cid:2) g ( i − j )( − j ζ ix ζ jy + g ( i − j )( − j v ix v jy (cid:3) + (cid:88) i,j (cid:2) g ( i − j )( ζ ix v jx + ζ iy v jy ) + g ( i − j )( − j ( ζ ix v jy − ζ iy v jx ) (cid:3) , (S.51)where we impose the same symmetry as in Eq.(S.25). By direct calculation, the following equations are obtained λg ( z ) + 2 ¯∆ g ( z ) = 0 , for z ≡ , ( λ + 4 γ ) g ( z ) − g ( z ) = 0 , for z ≡ , ( λ − γ ∆) g ( z ) + Bg ( z ) = h ( z ) , for ∀ z , ( λ + γ ¯∆) g ( z ) − g ( z ) + 2 g ( z ) − Bg ( z ) = 0 , for z ≡ , ( λ + γ ¯∆) g ( z ) − g ( z − − g ( z + 1) − Bg ( z ) = 0 , for z ≡ , (S.52)where ¯∆ g j ( z ) := g j ( z + 1) + g j ( z −
1) + 2 g j ( z ), ∆ g j ( z ) = g j ( z + 1) + g j ( z − − g j ( z ) and h ( z ) = ( − / δ z, − δ z, − ).We define the discrete Fourier transform g j ( z ) = (cid:26) (cid:80) k ˜ g je ( k ) e − ikz /N for z ≡ (cid:80) k ˜ g jo ( k ) e − ikz /N for z ≡ . (S.53)Eliminating g and g in Eqs.(S.52), we get the equation A II ˜ g ( k ) = u for the vectors ˜ g =(˜ g o ( k ) , ˜ g e ( k ) , ˜ g o ( k ) , ˜ g e ( k )) T and u = ( − i sin k, , , T , and the matrix A II is given by A II = µ, − γ cos k, B, − γ cos k, µ, , B − B, , µ, γ − / ( µ + 2 γ )) cos k , − B, γ + 2 / ( µ − γ )) cos k, µ (1 + 8 / ( µ − γ )) , (S.54)where µ = λ + 2 γ . We note that the Laplace transform C A ( λ ) is given solely by the function g and obtain thefollowing expression C A ( λ ) = (cid:88) i,j (cid:88) a g ( i − j ) (cid:104) ζ ia v ja J dc (cid:105) eq = ( − / (cid:88) i,a g ( i ) F ia = − (1 /N ) (cid:88) k ˜ g ( − k ) ˜ F ( k ) → ( T /π ) (cid:90) π dk cos ( k/ N ( k ) / D ( k ) , (S.55) N ( k ) = µ (cid:0) ( B + µ ) a − B (cid:1) + 2 (cid:0) B (2 µ − γa ) − γµ a (cid:1) cos k + a (cid:0) µ cos k + 8 γ cos k (cid:1) , (S.56) D ( k ) = ( B + µ ) (cid:0) ( B + µ ) a − B (cid:1) + 8 (cid:0) B γ a + µ ( a − (cid:1) cos k − γ a cos k, (S.57)where the function F ia is the same as that in case (I) and the Fourier transform ˜ F ( k ) is given by Eq.(S.36). We used a = 8 − γ + µ , a = − a + 4 and a = 4 + 4 γ − γ (8 + µ ). We analyzed the poles in the denominator by usingMathematica. The expressions for 6 poles for a small wave number and finite B are given as λ (cid:63) = − B ) γk / (4 + B ) + · · · (S.58) λ (cid:63) = − γ + 8(2 + B ) γk / (4 + B ) + · · · (S.59) λ (cid:63)σσ (cid:48) = − γ + σ (cid:112) − − B + 4 γ + σ (cid:48) ik/ (cid:112) B + · · · , (S.60)where σ, σ (cid:48) = ± . The power law decay is given by the pole λ (cid:63) which gives t − / in C A ( t ).1We finally consider C B ( λ ). One can directly check that the following expression is the solution of u λ,B u λ,B = − λ − / ( λ + 4 γλ + 4 + B ) (cid:88) j (cid:8)(cid:2) ¯ r x ( λ + 4 γλ + 4) − ¯ r y Bλ ( − j (cid:3) v jx + (cid:2) ¯ r x Bλ ( − j + ¯ r y ( λ + 4 γλ + 4) (cid:3) v jy + 2 B ( − j (¯ r x r jy − ¯ r y r jx ) (cid:9) . (S.61)From this expression, C B ( λ ) is proportional to (cid:104) ¯ r x r x (cid:105) eq or (cid:104) ¯ r y r y (cid:105) eq , which is O (1 /N ), and hence it is negligible inthe thermodynamic limit. [1] G. Basile, C. Bernardin, and S. Olla, Thermal conductivity for a momentum conservative model , Comm. in Math. Phys.,287