Anomalous heat transport in one dimensional systems: a description using non-local fractional-type diffusion equation

It has been observed in many numerical simulations, experiments and from various theoretical treatments that heat transport in one-dimensional systems of interacting particles cannot be described by the phenomenological Fourier's law. The picture that has emerged from studies over the last few years is that Fourier's law gets replaced by a spatially non-local linear equation wherein the current at a point gets contributions from the temperature gradients in other parts of the system. Correspondingly the usual heat diffusion equation gets replaced by a non-local fractional-type diffusion equation. In this review, we describe the various theoretical approaches which lead to this framework and also discuss recent progress on this problem.


I. INTRODUCTION
Transport of heat through materials is a paradigmatic example of non-equilibrium phenomena [1][2][3].When an extended system is attached to two reservoirs of different temperatures at its two ends, an energy current flows through the body from hot region to cold region.At the macroscopic level this phenomena is described by the phenomenological Fourier's law.Considering transport in one dimensional systems, Fourier's law states that the local heat current density j(x, t) inside a system at point x at time t is proportional to the gradient of the local temperature T (x, t): where κ is referred to as the thermal conductivity of the material.This law implies diffusive transfer of energy.To see this we note that the local energy density e(x, t) in a one dimensional system satisfies the continuity equation ∂e(x, t)/∂t = −∂j(x, t)/∂x.Inserting Eq. ( 1) in this continuity equation, and using the relation between the local energy density and the local temperature c v = ∂e/∂T (where c v represents the specific heat per unit volume), one finds the heat diffusion equation where we assume (for simplicity) no variation of κ with temperature.In usual three dimensional systems, the heat diffusion equation takes the form ∂ t T (x, t) = (κ/c v )∇ 2 T (x, t) and describes the evolution of the temperature field in bulk systems.The phenomenological macroscopic description provided by the equations in ( 1) and ( 2) has been used extensively to describe heat transfer phenomena in a wide class of physical systems.A natural question is to ask if it is possible to derive or establish Fourier's phenomenological law theoretically, starting from a complete microscopic description.The issue of deriving Fourier's law has been a long standing question and a very active field of research [1].Several theoretical as well as large scale numerical studies have been performed on different mathematical model systems to understand the necessary and sufficient conditions needed in the microscopic description to validate Fourier's law at the macroscopic level [2][3][4].Surprisingly, these studies suggest that Fourier's law is probably not valid in many one-dimensional systems and one finds that the thermal conductivity κ diverges with system size N as κ ∼ N α where 0 < α < 1 [2][3][4][5][6][7][8][9][10][11][12].This is referred to as anomalous heat transport (AHT).For α = 0, the transport is classified as being diffusive while α = 1 is referred to as ballistic transport [2,3].Recent developments in technology has made it possible to verify some of these theoretical predictions experimentally as well as numerically in real physical systems such as nano-structures, polymers, semiconductor films etc. [13][14][15][16][17][18][19][20], and these have provided further motivation and new directions of study.
Two approaches have mainly been used to look for signatures of anomalous heat transport (AHT): (i) the open system set-up in which a system is connected to heat reservoirs at different temperatures T L and T R at the two ends and (ii) the closed system set-up in which the isolated system is prepared in thermal equilibrium at temperature T and evolves according to Hamiltonian dynamics (or sometimes stochastic dynamics with same conservation laws).In the open system set-up, one usually considers the non-equilibrium steady state (NESS) and measures directly the steady state heat current j and the temperature profile T (x) in a finite system of N particles.For small ∆T = T L −T R , one finds the system size scaling j ∼ N α−1 (implying κ ∼ N α ) and a non-linear temperature profile.These are in contrast with Fourier's law which would predict j ∼ N −1 and a linear temperature profile.In the closed system set-up the idea is to look at the spreading of a heat pulse in a system in equilibrium.From linear response theory we expect that this would evolve in the same way as dynamical correlations of energy fluctuations in equilibrium.Studies on spreading of pulses and energy correlations in systems with AHT show that the process is super-diffusive, with scaling functions described by Lévy distributions [8,21,22].This contrasts systems described by Fourier's law where we expect diffusion and Gaussian propagators.Note that we expect in fact that the thermal conductivity κ obtained in non-equilibrium measurements should be related to equilibrium energy current auto-correlation functions via the Green-Kubo formula [3,23,24].This leads to the understanding of AHT as arising from the fact that the non-integrable long time tails in the auto-correlation function of the total current lead to the divergence of the thermal conductivity.
The natural question that arises for understanding systems with AHT is to find the replacements of Fourier's law in Eq. ( 1) and the heat diffusion equation in Eq. (2).The picture that has emerged from studies over the last few years [4,[25][26][27][28][29][30][31][32][33][34][35][36][37] is that Fourier's law gets replaced by a spatially non-local but linear equation wherein the current at a point gets contributions from temperature gradients in other parts of the system.This has the form j(x, t) = − dx K(x, x ) ∂T (x , t) ∂x , In (a), a localized heat pulse is introduced at some point in a system in thermal equilibrium and its subsequent time-evolution is observed.In (b), the system is attached to two heat reservoirs at different temperatures and the NESS properties such as current and temperature profile are studied.
where now the thermal conductivity is replaced by the non-local kernel K(x, x ).This then leads to a corresponding non-local fractional-type equation for the time evolution of T (x, t).An important difference from the heat diffusion equation is that the fractional-type equation takes different forms in the closed system set-up (infinite domain) and the open system set-up (finite domain).In the infinite domain the evolution of a localized temperature pulse is described by a fractional-type diffusion equation where the fractional operator should be interpreted by its action on plane wave basis states: (−∆) ν/2 e ikx = |k| ν e ikx , with 1 < ν < 2. However it should be noted that the correpsonding Lévy-stable distribution is valid only over the scale x < ∼ t 1/ν .As we will see, the evolution of a heat pulse is restricted to a domain |x| < ct, determined by the sound speed c.For the open system, the precise form of the fractional equation is dependent on the details of boundary conditions.In this review we discuss these developments as well as open questions.
The plan of the review is as follows.In Sec.(II) we discuss the various signatures of AHT in the closed and open set-ups.In Sec.(III) we discuss two theoretical approaches that have been used to understand various aspects of anomalous transport.One of these is a phenomenological approach based on the idea that the heat carriers perform Lévy walks instead of random walk.The second approach is a microscopic one, though still phenomenological, and is based on fluctuating hydrodynamics and applicable to Hamiltonian systems.For a class of stochastic models, it has been possible to provide a complete microscopic derivation of the fractional heat equation in the context of both the closed and open system set-ups.These results are described in Sec.(IV).In the last part of this section we address the difficult issue of treating arbitrary boundary conditions and discuss a heuristic formulation that uses linear response ideas and fluctuating hydrodynamics to arrive at a general form of the kernel K(x, x ) in Eq. ( 3).Finally we conclude in Sec.(V) with a summary of the results presented and some of the outstanding open questions.

II. SIGNATURES OF ANOMALOUS HEAT TRANSPORT
In the theoretical study of anomalous energy transport in one dimension, one usually considers simple yet non-trivial model systems of interacting particles.Let us consider N particles of unit masses, with positions and momenta given, respectively, by q and p , for = 1, 2, ..., N .One often starts with the following microscopic Hamiltonian: where V (r) is a nearest neighbor interaction potential, and the extra variables q 0 and q N +1 are introduced to incorporate different boundary conditions (BC).For example, fixed BC corresponds to q 0 = 0, q N +1 = 0 while free BC corresponds to setting q 0 = q 1 , q N +1 = q N .The particles in the bulk of the system satisfy Hamiltonian equations of motion One of the well-studied choices for the potential is to take V (r) = k 2 r 2 /2 + k 3 r 3 /3 + k 4 r 4 /4 which leads to the Fermi-Pasta-Ulam-Tsingou (FPUT) model.Another popular choice is the alternate mass hard particle gas which is not in the standard form of Eq. (5).In this model one considers a chain of point particles with masses which alternate between two fixed values, say m 1 , m 2 , and which collide via elastic collisions conserving energy and momentum.For generic interaction potentials V (r) it is expected that the system has three conserved quantities, namely volume of the system (alternatively the total number of particles), total momentum and total energy.Corresponding to each conserved quantity one can write a local continuity equation.For instance, the local energy defined on bulk points as satisfies a continuity equation This equation gives a microscopic definition of the energy current.For quadratic V (r), i.e harmonic chains, there are a macroscopic number of conserved quantities and transport becomes ballistic.In this case a number of studies have considered augmenting the Hamiltonian dynamics with a stochastic component such that the system again has only three conserved quantities [9,[29][30][31].In this case one again recovers the typical features of anomalous transport and several exact results are possible.In this review we will discuss results for both Hamiltonian and stochastic systems.
There are two possible approaches for studying transport properties of a system [3,4].A schematic of the two set-ups is shown in Fig. (1): A. Closed system set-up -in this case, an isolated system is prepared in thermal equilibrium at some temperature T described by the canonical distribution where Z = dqdpe −H/T is the partition function.For any initial condition chosen from this distribution the system evolves according to the pure Hamiltonian dynamics (or the conservative stochastic dynamics).
Transport properties are usually probed by studying the form of spatio-temporal correlation functions of the conserved quantities (volume, momentum, energy) or the decay with time of the energy current auto-correlation function.Another approach that has been used is to study the spreading of an initially localized perturbation in the equilibrated system [see Fig. (1a)].In the closed system set-up one takes the system to be infinite or, in numerical studies, N to be sufficiently large such that the correlations are not affected by the boundaries at the maximum observation times.

B.
Open system set-up -in this case, one considers finite systems attached at the two boundaries to heat reservoirs at different temperatures [see Fig. (1b)].The heat reservoirs are modeled by adding extra force terms to the usual Hamiltonian equations of motion of the boundary particles.One of the standard choices is to consider Langevin type baths, wherein the additional forces consist of a dissipative term and a white noise term, which are related via a fluctuation-dissipation relation.The system is "open" in the sense that energy can flow in and out of the system, though we note that locally in the bulk we still have energy conservation.When the temperatures of the heat reservoirs are different, the system eventually reaches a NESS in which a heat current flows across the system.The main focus of this approach has been to search for anomalous features in the NESS by looking at observables such as the heat current j = j(x, t) neq open and temperature profile obtained from x neq open (the averages are computed in the NESS).There have also been attempts to understand the relaxation to NESS and look at correlations and large deviation properties of the NESS.
In the following sub-sections, we describe various signatures of AHT observed in both the set-ups.• Slow decay of energy current auto-correlations: A commonly followed approach for determining the N dependence of j or equivalently the thermal conductivity κ, is to use the closed system Green-Kubo (GK) formula [23,24]: where J(t) = x j(x, t), with j(x, t) defined in Eq. ( 8), is the total current in the system.The average ... eq closed is evaluated with initial conditions chosen from a thermal distribution and time-evolution given by the closed system dynamics.This formula relates the thermal conductivity κ to the integral of the equilibrium heat current auto-correlation function C J (t) = N −1 J(t)J(0) eq closed .Numerical simulations as well as several theoretical treatments find that C J (t) in a closed system generically decays with time as a power law C J (t) ∼ t α−1 with 0 ≤ α ≤ 1 [2,3,7,9,12,33,[38][39][40][41][42][43][44][45][46][47][48][49][50][51].As an example we show in Fig. (2) data from simulations [7] of the alternate mass hard particle gas, where we see a decay with α ≈ 0.33.Such a power-law time dependence implies, from Eq. 10, a divergent thermal conductivity.To see the dependence on system size one heuristically puts a cutoff t N ∼ N in the upper limit of the time integral, the argument being that this is the time taken by sound modes to explore the full system of size N .Performing the time integral in Eq. (10) with this cut-off, one finally gets κ ∼ N α .An interesting example where this procedure fails has been pointed out in a recent work [52,53].
• Super-diffusive spreading of initially localized energy pulse: Here one looks at the spreading of a localized energy pulse in a thermally equilibrated system.One takes an initial configuration chosen from a thermal distribution with average local energy e 0 = e(x) eq closed , uniform across the system.Imagine putting an extra amount of energy 0 to a few particles in a region inside the bulk to create a pulse of excess energy locally.As the system evolves according to the closed system dynamics, this localized energy perturbation starts spreading across the system.Let (x, t) represent the excess energy density (above e 0 ) at the point x and at time t (averaged over the initial distribution).This quantity starts as a δ-function at t = 0 and then starts spreading with time.Note that dx (x, t) = 0 , the total injected energy is conserved under the closed system dynamics.For a diffusive system, the perturbation would evolve according to the diffusion equation ∂ (x, t)/∂t = D∂ 2 (x, t)/∂x 2 and in macroscopic length-time scales, the perturbation profile at time t would be given by a Gaussian For a system with AHT, one instead finds the following scaling form [4,8] ( with a scaling exponent 1/2 < γ < 1.The two limits γ = 1/2 and 1 correspond respectively to diffusive and ballistic transport.In Fig. 3 we show results for energy pulse spreading obtained in [8] for the alternate mass hard particle gas model.The main plot shows the scaling x ∼ t γ , with γ = 3/5 of the central part of the distribution.The central part of the distribution was found to fit to the Lévy function which is the propagator of Eq. ( 4) with µ = 1/γ.The mean square deviation (MSD) defined as with mean taken as zero, was seen to scale as σ 2 e (t) ∼ t β , with β = 4/3, as opposed to a diffusive system with β = 1.It was also noted that the MSD width exponent, β, is related to the thermal conductivity exponent α as  the Lévy distribution one gets, in the regime t γ << x < ∼ t, the scaling form G(u) ∼ 1/u 1+1/γ .Using these asymptotics and computing σ 2 e (t) = t 0 dx x 2 t −γ G(x/t γ ) gives us the leading behavior σ 2 (t) ∼ t 3−1/γ which then leads to the relation β = 3 − 1/γ.Observations from several other numerical simulations have confirmed the super-diffusive behavior [8,[54][55][56][57][58][59].
• Super-diffusive evolution of density correlations: The anomalous signature discussed in the previous point can also be observed alternatively by looking at the spreading of the equilibrium spatio-temporal correlation function of the energy density e(x, t) defined as C e (x, t) = e(x, t)e(0, 0) − e(x, t) e(0, 0) , where the average is taken over the equilibrium initial conditions.For diffusive systems this correlation has the Gaussian form in Eq. ( 11), while for systems with AHT this has the scaling form in Eq. ( 12) and one again has super-diffusive growth of the MSD [21], now defined as This MSD can be related to σ 2 e (t) defined above, using linear response theory and both have ∼ t β scaling.In the case of AHT, observing the scaling form in Eq. ( 12) usually requires one to subtract contributions of sound modes which travel ballistically.The theory of nonlinear fluctuating hydrodynamics (NFH) provides a framework in which one can systematically describe the super-diffusive scaling of the correlation [22,47,[60][61][62][63].This theory is based on writing hydrodynamic equations for the conserved quantities in the system which for the Hamiltonian in Eq. ( 5) are the total energy, total momentum and the total number of particles (or volume).This framework of NFH is discussed in detail in Sec.(III B).A connection can be made between the superdiffusive scaling (σ 2 c (t) ∼ t β ) of the energy correlations and the power-law decay, ∼ t α−1 , of the current-current correlations [4,58,59], which can be seen as follows.Starting from the continuity equation for energy, one can obtain the relation [61,62] on the infinite line Multiplying by x 2 on both sides and integrating over all the range of x one gets Assuming the expected forms σ 2 (t) ∼ t β and C J (t) ∼ t α−1 we get the relation α = β − 1.

B. Signatures in the open system set-up
• Diverging thermal conductivity: As discussed above, in the open system set-up, one connects the system at the two boundaries to heat reservoirs at unequal temperatures T = T r .A common model for baths is to write system, Fourier's law predicts that the heat current j should be equal to ÿrT, with a that depends on microscopic  Langevin dynamics for the boundary particles involving dissipation and noise term satisfying the fluctuationdissipation relation.For a chain of interacting particles described by the Hamiltonian in Eq. ( 5) the equations of motion for the boundary particles would read where f i = −∂H/∂q i .The noise terms ξ ,r are Gaussian white noise terms, with zero mean and correlations ξ (t)ξ (t ) = 2λT δ(t − t ) and ξ r (t)ξ r (t ) = 2λT r δ(t − t ).The remaining particles evolve according to Eq. (6).After a long time the system reaches a non-equilibrium steady state (NESS) and we can measure the steady state current j as average of the local current j(x, t) defined through Eq. ( 8).In the steady state this will be independent of time as well as the bond where we measure the current.One can then check if the system size N scaling of this steady state current j has the expected form j ∼ N α−1 , where α < 1 for anomalous systems.
Alternatively one can define the κ = jN /(T − T r ) and see how this scales with N .For a large class of nonlinear interaction potentials, it has been observed that the thermal conductivity κ ∼ N α with 0 < α < 1 for large N [6,7,10,11,63,64].As an example, we show in Fig. (4) data from [10] for the FPUT-β chain, where one finds α ≈ 0.33.
• Non-linear temperature profile: The local temperature at a site on the lattice can be defined through the relation T i = p 2 i /m , where the average is taken in the NESS.For diffusive systems, the temperature profile obtained would be linear for small ∆T = T − T r , as expected from solving Fourier's law with a constant κ.It is important to note that non-linear temperature profiles can also be obtained in case of diffusive transport if the thermal conductivity κ is temperature-dependent and ∆T is large.On the other hand, for many systems with AHT, one finds a strongly non-linear temperature profile even when ∆T is made arbitrary small [5,10,11,26,34,36,65].Quite often the profiles are characterized by divergent slopes at the boundaries.In Fig. (5) we show the temperature profile in the FPUT-β model and one can see the characteristic non-linear nature.Note that the definition of local temperature makes sense (and is useful) only if this temperature predicts correctly other local observables, for example higher moments of the velocity.This was also verified in [10] and also shown in Fig. (5).Typically one finds that the temperature difference δT (x) = |T − T (x)| scales as (δx) µ , with distance δx from the boundary, where 0 < µ ≤ 1.The exponent µ has been referred to as the meniscus exponent [66].This exponent is non-universal in the sense that it depends on details of boundary conditions, unlike the conductivity exponent α.
• Green-Kubo-type relation for open systems: Analogous to the Green-Kubo formula in the closed system set-up given by Eq. ( 10), an exact formula exists in the open system set-up that relates the current response to a small temperature difference ∆T = T − T r .This is given by [67] lim The time auto-correlation J(t)J(0) eq open is computed by averaging over equilibrium initial conditions as well as the open system dynamics which includes the stochastic baths (at equal temperatures).This formula is valid for a finite size system.We note that for systems with AHT, unlike with Eq. (10), in the open set-up we do not require the use of an upper cut-off t N ∼ N for estimating the size dependence of conductivity.In this case the linear response current can be evaluated directly from Eq. ( 20) for any finite system of size N and thereby one can verify the form j/∆T ∼ N α−1 .This approach has been discussed for example in [63,64].It was observed in [64] that, for the so-called random collision model studied by them, both J(t)J(0) eq closed and J(t)J( 0) eq open showed a t −2/3 decay at times t < ∼ N .However the exponential decay for the open case begins at t N ∼ N while for the closed system (with periodic boundary conditions) this begins at t N ∼ N 3/2 .This was understood as arising from the time scale associated to the spreading of sound modes.Note that if we put the cut-off t N ∼ N 3/2 as the upper limit in the time-integral of Eq. ( 10) then we would get the wrong conductivity exponent.In order to get the correct exponent in the closed system set-up, one has to by hand set the cut-off at t N ∼ N based on consideration of the practical transport set-up which has baths at the boundaries.Recently, in a model system of AHT the relation in Eq. 3 has been established using the above formula and a heuristic approach based on fluctuating hydrodynamics [36].An explicit expression of the kernel was obtained for a specific model, using which one can understand the divergence of κ as well as the singular features in the temperature profile .A detailed discussion of this method is given later in Sec.(IV B 3).

III. PHENOMENOLOGICAL APPROACHES FOR ANOMALOUS HEAT TRANSPORT
In this section we will discuss two different approaches that have tried to understand the various aspects of AHT mentioned above.The first is a completely heuristic approach where one assumes that the heat carriers perform Lévy walks instead of random walk which is expected for diffusive heat transfer.This method has been used to explain spreading of localized energy pulses, steady state properties and current fluctuations [8,39,57,66,[68][69][70][71].The second approach is a microscopic one where one starts by writing hydrodynamic equations for the conserved quantities of the Hamiltonian dynamics.One then phenomenologically adds noise and dissipation terms satisfying fluctuation dissipation relations and this allows one to study equilibrium fluctuations in the system.In particular, using the formalism of fluctuating hydrodynamics, one can compute dynamical correlation functions which contain information on AHT.

A.
Lévy walk description of anomalous heat transport

Lévy walk description in the closed set-up
In this description one thinks of energy as being carried by Lévy walkers, each of which carry a fixed amount of energy.It follows that the local energy density and energy current at any point can be taken to be directly proportional to, respectively, the particle density and current.Let us also assume that the local temperature is proportional to the local energy density and hence to the density of particles.
Definition of the Lévy walk [72][73][74]: At each step of the walk, a particle chooses a time of flight τ from a specified distribution, φ(τ ), and then moves a distance x = cτ at a fixed speed c, with equal probability in either direction.More generally one can consider the velocity c to be drawn from a distribution.Let P (x, t)dx denote the probability that the particle is in the interval (x, x + dx) at time t.Note that P (x, t) also includes events where the particle is crossing the interval (x, x + dx), in addition to the events in which the particle lands in the interval at time t.If a particle starts at the origin at time t = 0, the probability P (x, t) satisfies where ψ(τ ) = ∞ τ dτ φ(τ ) is the probability of choosing a time of flight ≥ τ .Here we consider Lévy walkers with a time-of-flight distribution which decays, at large times, like a power law φ(t) A t −ν−1 with A = νt ν 0 .For this range of ν the mean flight time t = ∞ 0 dt t φ(t) = t 0 /(ν − 1) is finite but t 2 = ∞.Some properties of the Lévy walk: Taking the Fourier Laplace transform where For asymptotic properties it is useful to find the form of P (k, s) for small k, s.The Laplace transform φ is given by: and Γ(u) is the Gamma-function.Hence we get: where d = bA/(2 t ).Taking the inverse Fourier-Laplace transform of this gives us the propagator of the Lévy walk on the infinite line.This corresponds to a pulse whose central region is a Lévy-stable distribution with a scaling x ∼ t 1/ν .This can be seen by expanding Eq. ( 25) for ck/s << 1 to get The difference with the Lévy-stable distribution is that the Lévy-walk propagator has ballistic peaks of magnitude t 1−ν at x = ±ct and vanishes outside this.The overall behavior of the propagator is as follows [72]: The time evolution of the Lévy-walk propagator, obtained from direct simulations of the Lévy walk, is shown in Fig. (6).
We also plot the Lévy-stable distribution obtained by taking the Fourier transform of P (k, t) = e −c cos(νπ/2)|k| ν t .
Various moments of the distribution can be found using the relation or its Laplace transform given by . Using Eq. ( 25) we get in particular the following leading behavior We see that for 1 < ν < 2 the motion is super-diffusive [73,74].The most interesting characteristics to note about the Lévy walk are the fact that the probability distribution has finite support (|x| ≤ ct), in the bulk it coincides with the Lévy distribution with scaling x ∼ t 1/ν and finally the mean square displacement (MSD) Note that the usual Lévy stable distribution has a diverging second moment, however the Lévy walk has a finite MSD and this follows from the finite support |x| ≤ ct of the corresponding distribution.Indeed, on using this cutoff and the power-law form of the Lévy near the cut-off (see Eq. ( 26)) gives us the expected scaling exponent β = 3 − ν.The inset shows a plot of the mean square displacement and the fourth moment and a comparison with the exact asymptotic forms (dashed lines) given by Eqs.(27,28).In all plots the time to and c are set to one.

Lévy walks and AHT:
The first proposal suggesting the Lévy walk model to describe anomalous heat transport was made in [68].This idea was tested for a microscopic model in [8] where it was shown that the spreading of a heat pulse in a thermally prepared alternate mass hard particle gas was super-diffusive and is well-described by the Lévy walk model.In Fig. (3) we show the evolution of a localized perturbation.The main plot shows the x ∼ t 3/5 scaling of the central part of the distribution while the inset shows a fit to the expected Lévy distribution (for a LW with ν = 5/3) with a single fitting parameter.It was also shown that the MSD of the energy has the scaling ∼ t 4/3 as expected from the relation β = 3 − ν for LW.Finally it was proposed using linear response ideas that the exponent β and the conductivity exponent α should be related as α = β − 1 which gives α = 1/3 for the present system.This agrees with known results for the alternate mass hard particle gas.The validity of the Lévy walk description of pulse propagation was further verified in [39] for a hard particle gas interacting via a square well potential and in [57] for the FPUT chain.All these cases were described by the same Lévy-walk exponent ν = 5/3.

Lévy walk description of the open set-up
We now discuss the case of the open system consisting of a finite segment (0, L) that is connected to two reservoirs at the ends.The use of the Lévy walk model to study NESS properties in AHT was first proposed in [66] where the authors considered a finite lattice of N sites containing a collection of Lévy walkers.The system was connected at it's two ends to infinite reservoirs that contained sources emitting Lévy walkers at fixed constant rates.A Lévy walker crosses from the reservoir into system with probability one, but while exiting from system into reservoir, it can get reflected with probability R. A particle exiting the reservoir is eliminated.The authors in [66] considered a discrete version and studied this problem numerically.The strategy was to write appropriate master equations for the probability evolution and obtain the steady state solution numerically.One of the main observations in the paper was that the NESS profile for P (x) was non-linear and was singular at the boundaries.In Fig. ( It was noted in [66] that the value R = −0.1 was unphysical but made mathematical sense in the master equation (see [66] for further discussions on this point) and gave the best agreement with the momentum exchange simulation profile.Some exact results were obtained for the Lévy walk model with particle reservoirs, for the special case of perfectly transmitting boundary walls (i.e R = 0) [69] which we now describe.We note that for the Lévy walker, at any given time, a particle could either be passing over a point x or could have landed precisely at the point.Hence, in addition to the probability density P (x, t), it is convenient to define the probability Q(x, t)dxdt that a particle lands precisely between x to x + dx in the time interval (t, t + dt).We now specify the boundary conditions required to set up a non-equilibrium current carrying steady state.For this, we consider the region x ≤ 0 as the left reservoir with Q(x, t) = Q l for all points in this region.Similarly, we set Q(x, t) = Q r in the region x ≥ L corresponding to the right reservoir.In the steady state, the distributions become time-independent and Q(x, t) = Q(x), P (x, t) = P (x) satisfy [69] where as mentioned earlier after Eq. ( 21).The terms on the right hand side of the above equation for Q(x) represent different contributions.The first term represents the contributions from walkers that start from various points y and land at x.The second and the third term represent contributions from walkers starting, respectively, from left and right reservoirs and landing at x. Similarly, the terms on the right hand side of Eq. (31) for P (x) can also be interpreted in the same way except now, events in which the walkers are passing over x, in addition to the events in which they land at x at a given time, also contribute.
Interestingly, it turns out that the problem of finding Q(x) can be related to the problem of the escape probability [75] of a Lévy walker on the interval (0, L).Let H(x) denote the probability with which a Lévy walker, starting at position x, arrives at the left reservoir (region x < 0) before arriving at the right reservoir (region x > L).It can be shown that H(x) satisfies [69] The probability Q(x) can now be expressed in terms of H(x) as which can be checked easily to satisfy Eq. (30).
If one considers a Lévy flight with distribution ρ(z) = [φ(z/c) + φ(−z/c]/(2c) of steps z, the probability H(x) that starting at x, the flight hits first the left bath satisfies exactly Eq. (32).Hence by following the same mathematical steps as in [75] to study equations such as (30) or (32), one can show that, in the large L limit, the solution Q(x) of (30) (and H(x) of ( 32)) satisfies with and H(0) = 1 and H(1) = 0 for (32)] with a solution of (33), for a φ(τ ) decaying as in (22), which satisfies We can integrate this equation to get Q(x), with the integration constant and B being then determined from the boundary conditions where 2 F 1 (a, b, c, z) is the hypergeometric function.For large L, the right hand side of Eq. ( 31) is dominated by the range |y − x| L and therefore The exact results of Eqs.(34) have been verified in [69] from direct numerical solution of Eqs.(30,31) and it was noted that density profiles were similar to the temperature profiles seen in AHT.
Next we discuss the steady state current j(x) which is given by This can be seen to be the difference between the flow from left to right and from right to left.The contribution from 0 < y < ∞ to the integral corresponds to particles crossing the point x from left to right which is obtained by multiplying the density of particles at x − y with the probability ψ(y/c) that they have a flight time larger than y/c.The contribution from −∞ < y < 0 to the integral corresponds to a similar right-to-left current.After performing a partial integration and using the boundary conditions Using Eq. ( 33) it is easy to see that dj/dx = 0 which implies that the current in the steady state is independent of x, as expected.Hence, evaluating the current at x = 0 and using Eq. ( 34), we get for large L From Eq. ( 27) we then get the relation α = β − 1, between the conductivity exponent of AHT and the MSD exponent for Lévy-walk diffusion.This relation for Lévy diffusion was pointed out in [68] and verified in simulations in 1D heat conduction models [8,21].A derivation based on linear response theory has been given in [59].Finite size corrections to the results in Eqs.(34,39) were recently obtained in [76].
In the large L limit by using Eq. ( 36) in Eq. ( 38) we obtain Above equation is the analogue of Fourier's law Eq.( 1) with the important difference that in the linear response regime the current at a point gets contributions from the temperature gradients at other parts of the system as well.
The above treatment can be generalized for arbitrary values of the reflection probability R [37] and this leads to the following general non-local form of the current where Remarkably we note that for ν = 3/2 (α = 2 − ν = 1/2), the expression above is identical to the expression for K R (v, v ) with v = x/L, v = y/L, given later in Eq. (112).

B. Nonlinear fluctuating hydrodynamics description of anomalous heat transport
We now discuss a completely different approach for understanding AHT.In this approach the starting point is the Hamiltonian dynamics of the system.The idea is to consider the effective dynamics of the slow conserved fields using some coarse graining.One finds that the evolution of small fluctuations around equilibrium can be described by fluctuating hydrodynamics [22,47,[60][61][62].Solving these equations using mode coupling theory, detailed predictions can be made on the form of equilibrium spatio-temporal correlation functions of the conserved fields.In particular, we will see that it predicts the super-diffusive spreading of energy perturbations with Lévy-law scaling, and the slow decay of energy current auto-correlation functions.We will here describe the theory for generic anharmonic systems with three conserved quantities, namely volume, momentum, energy [61] and present some numerical results which verify the predictions of the theory.
Let us consider N particles of unit masses with positions and momenta denoted by {q( ), p( )}, for = 1, . . ., N .The particles move on a ring of size L so that we have the boundary conditions q(N + 1) = q(1) + L and p(N + 1) = p(1).The Hamiltonian is taken to be where we have defined the stretch variables r( ) = q( + 1) − q( ).It is easy to see from the Hamiltonian equations of motion that stretch r( ), momentum p( ), and energy ( ) are locally conserved and hence satisfy corresponding continuity equations.In the continuum limit, these equations take the form where the label index has been denoted by the corresponding continuous variable x and P (x) = −V (x) is the local force.Assume that the system starts in a state of thermal equilibrium at zero total average momentum characterized by the temperature (T = β −1 ) and pressure (P ), which fix the the average energy and average stretch of the chain.
The distribution corresponding to this ensemble is Since the fields r(x, t), p(x, t) and e(x, t) satisfy continuity equations, they evolve slowly suggesting a slowly evolving local equilibrium picture.We consider small fluctuations of the conserved quantities about their equilibrium values, u 1 (x, t) = r(x, t) − r eq , u 2 (x, t) = p(x, t) and u 3 (x, t) = (x, t) − eq .Inserting these into Eqs.(44) one obtains ∂ t u α = −∂ x j α , where j α are the corresponding Euler currents which are functions of u α s.Expanding these currents to second order in the fields as j α = β A αβ u β + β,γ H α βγ u β u γ , and then adding dissipation and noise terms (to ensure thermal equilibration) one arrives at the following noisy hydrodynamic equations where repeated indices are summed over.The noise and the dissipation matrices, B, D, are related to each other by the fluctuation-dissipation relation DC + C D = B B T , where the matrix C corresponds to equilibrium correlations and its elements are C αβ (x) = u α (x, 0)u β (0, 0) .It is useful to define normal modes of the linearized equations (dropping u 2 terms in Eq. ( 46)) through the transformation (φ − , φ 0 , φ + ) = φ = R u, where the matrix R acts only on the component index and diagonalizes A, i.e.RAR −1 = diag(−c, 0, c).The diagonal form implies that there are two sound modes, φ ± , traveling at speed c in opposite directions and one stationary but decaying heat mode, φ 0 .The quantities of interest are the equilibrium spatio-temporal correlation functions C ss (x, t) = φ s (x, t)φ s (0, 0) , where s, s = −, 0, +.Because the modes separate linearly in time, one argues that at large times the off-diagonal matrix elements of the correlator are small compared to the diagonal ones and that the dynamics of the diagonal terms decouples into three single component equations.After including the non-linearity it is seen that to leading order the equations for sound modes have self-coupling terms of the form φ 2 ± .These then have the structure of the noisy Burgers equation, for which the exact scaling function, denoted by f KPZ , are known.For the heat peak the self-coupling coefficient vanishes for any interaction potential.Thus one has to study the sub-leading corrections, and calculations using the mode-coupling approximation result in the symmetric Lévy walk distribution, with a cut-off at x = ct.While this is an approximation, it seems to be very accurate.For the generic case of non-zero pressure, i.e.P = 0, which corresponds either to asymmetric inter-particle potentials or to an externally applied stress, the prediction for the left moving, resp.right moving, sound peaks and the heat mode are where f KPZ (x) is the KPZ scaling function discussed in [61,77], and tabulated in [78].The scaling function f ν LW (x) is given by the Fourier transform of the Lévy characteristic function e −|k| ν , with a cut-off at x = ct.The scaling parameters λ s and λ e are known explicitly.On the other hand for an even potential at zero pressure i.e.P = 0, all self-coupling coefficients vanish.As a result the scaling solutions within mode-coupling approximation change and one obtains where f G (x) is the unit Gaussian with zero mean.The scaling parameters λ 0 s is not known from microscopics while λ 0 e is known explicitly in terms of λ 0 s .Here we present molecular dynamics simulation results for the FPUT chain that were obtained in [63] which verify the predictions of NFH.In Fig. (8, top panel), the two-point correlation functions C 00 (x, t), C ++ (x, t) and C −− (x, t) are plotted as a function of x for three values of time t = 800, 2400 and 3200.The parameters used in this plot are k 2 = 1.0, k 3 = 2.0, k 4 = 1.0, T = 5.0, P = 1.0 for which one gets c = 1.80293 and we also see there is a good separation of the heat and sound modes.In Fig. (8, bottom panel) we also find an excellent collapse of the heat mode and the sound mode data with the expected scalings .The scaled data for the heat mode fits very well to the Lévy-scaling function whereas the same for the sound-mode still shows some asymmetry but is quite close to the KPZ function.
The numerically estimated values of the constants λ s,e are λ s = 0.46 and λ e = 5.86.These are in close agreement to the theoretically obtained values λ s = 0.396 and λ e = 5.89.

IV. STOCHASTIC MODELS: EXACT RESULTS ON FRACTIONAL EQUATION DESCRIPTION
It is now well understood that conservation laws play an important role in observation of super-diffusive transport in one-dimensional systems.As we saw in the previous section, NFH provides some understanding of the emergence of Lévy-walk behaviour, which seems to capture several aspects of anomalous transport.However, providing a completely rigorous microscopic derivation of the Lévy-walk picture in a Hamiltonian model has been difficult, though there have been some attempts [79].While generic non-linear Hamiltonian models are difficult to analyze, analytical results have been obtained for harmonic chains whose Hamiltonian dynamics is perturbed by stochastic noise that breaks integrability of the system [9,30,52].These stochastic models attempt to mimic nonlinear chains and for these models, several exact results both in the closed system set-up and the open system set-up have been obtained.In particular one can rigorously establish non-local response relation Eq. ( 3) and the fractional diffusion equation Eq. ( 4).There are two widely studied stochastic models which we discuss below.

A. Harmonic chain with volume exchange
This model is defined on a one dimensional lattice where each site carries a 'stretch' variable η i , i ∈ Z and the energy of the system is x/(λ e t)  8. Top Panel: Plots of the heat mode correlation C00(x, t) (central peaks) and the sound mode correlations C++(x, t) and C−−(x, t) (right and left moving peaks) in the FPUT chain, at three different times, for the parameter set with k2 = 1.0, k3 = 2.0, k4 = 1.0, T = 5.0, P = 1.0 and system size 16384.The speed of sound was c = 1.80293.We see that the heat and sound modes are well-separated.The numerical data in this plot were obtained by averaging over around 10 6 − 10 7 initial conditions.Bottom Panel: The heat mode (a) and the left moving sound mode (b) correlations, respectively, C00(x, t) and C−−(x, t) are plotted at different times, using a Lévy-type-scaling for the heat mode and KPZ-type scaling for the sound mode.Here we observe a very good collapse of the data at different times.Moreover, we observe a good fit to the Lévy-stable distribution with λe = 5.86 and a reasonable fit to the KPZ scaling function, with λs = 0.46.The parameters used in this plot are k2 = 1.0, k3 = 2.0, k4 = 1.0, T = 5.0, P = 1.0.(Adapted from Das et.al.with permission from [63] Copyright (2014) by American Physical Society).
a stochastic exchange part where ηs from any two randomly chosen neighboring sites, are exchanged at a constant rate γ.We refer to this model as Harmonic chain with volume exchange (HCVE).This model was introduced in [30] where it was shown that the energy current auto-correlation decays as ∼ 1/ √ t, implying super-diffusive transport.It is easy to see that this system has only two conserved quantities namely, the total "volume" i η i and the total energy i η 2 i .The evolution of the density fields corresponding to these conserved quantities at the macroscopic length and time scales was studied in [62] using NFH, where it has been shown that this model has two normal modes -one diffusive sound mode and a 3  2 -asymmetric Lévy heat mode.Subsequently, it was rigorously shown that the local energy density e(x, t) satisfies a (3/4)-skew-fractional equation [31] ∂ with ∆ as the usual Laplacian operator.In the Fourier domain, defined by e(k, t) = ∞ −∞ e(x, t)e ixk dx, the above equation reads as Note that for the diffusive case the analogue of the above equation would be ∂ t e(k, t) = −Dk 2 e(k, t).The above results suggest that, in the open set-up where the system is connected to two reservoirs at different temperatures, this model would exhibit anomalous scaling of the steady state current j with system size N .In [30], it has been numerically shown that indeed j ∼ 1/ √ N .Recently, an understanding of the open system was achieved using the fractional equation description, which we now discuss [34].An aspect that we will point out here is that the fractional-equation-type description in the open-set up is strongly dependent on boundary conditions (fixed or free or mixed).
For the open system case, we consider a finite lattice of size N , connected to two thermal reservoirs at temperatures T and T r on the left and right boundaries.The dynamics of the η i , i = 1, 2, ..., N now gets modified to + stochastic exchange at rate γ. ( The Langevin terms at the boundaries i = 1 and i = N appear due to the baths and ξ ,r (t) are independent Gaussian white noises with mean zero and unit variance.We consider fixed boundary conditions η 0 = η N +1 = 0. Our main interest is to obtain an equation in this finite system, analogous to Eq. ( 51), to describe the evolution equation of the temperature profile .To do this we first define the local temperature T i (t) = η 2 i (t) and the off-diagonal correlations C i,j (t) = η i (t)η j (t) , i = j, which characterize the non-equilibrium state of the system.Interestingly, it turns out that the equations for two point correlations do not depend on higher order correlations and this property leads to the model's solvability.The evolution of these quantities in the bulk (2 < i, j < N − 1) can be obtained from Eq. 53 as: The equations involving the boundary terms are given in [34].Note that in an infinite system, we get the same set of equations with i, j ∈ Z.For the finite open system, solving the above equations exactly seems to be difficult.However, it was observed numerically [34] that for large N the temperature field T i (t) scales as T i (t) = T i N , t N 3/2 and the correlation field C i,j (t) scales as Inserting these into (54), and expanding in powers of 1/ √ N , we find at leading order the following equations where the scaling variables Note that for the isolated infinite system, one can follow the same procedure as above, but now replacing the scale parameter 1/N → a where a is the lattice spacing, to obtain the same set of Eqs.55-57 with a different domain {−∞ ≤ u ≤ ∞ ; − ∞ ≤ v ≤ ∞}.These equations can be solved by Fourier transforms to get a skew fractional evolution equation for T (v, τ ) of the same form as Eq. ( 52).Defining Fourier transforms Solving the first equation ( 58), with the condition that correlations vanish at u = ±∞, we get The equation ( 59) relates the constant A k to Tk :
Using Eqs.(61,62) in Eq. ( 60) we get the infinite line equation in Eq. (52).We now go back to the open system case where the solution is more non-trivial.To solve these equations in the open set-up, we proceed as done for the regular diffusive heat equation, and write the solution as sum of a steady state part and a relaxation part where we have defined z = 1 − v.We note that under this transformation, the "anti-diffusion" Eq. ( 55), becomes a diffusion equation, with v as the time variable and z the space variable.The relaxation part satisfies the equations given in Eqs.(55,56,57), while the steady state part satisfies these equations but with ∂ τ T ss (v) = 0.The boundary conditions for the steady state part are given by [34] C ss (u, where we have used Eq. ( 57) to identify J = 2C ss (u = 0, z) as the NESS current which gets determined by the boundary conditions for T ss (v).In terms of the original unscaled variables, the true current is given by j ss = J/ √ N .The solution of the steady state equations is given by [34] T  70).The first 4 eigenvalues are completely real and distinct.The higher eigenvalue comes in complex conjugate pairs.For large µn ∼ (nπ) 3/2 (1 ± i).For smaller n, there is a deviation from asymptotic scaling due to finite definition of the operator.(Adapted from Kundu et.al.with permission from [34] Copyright (2018) by American Physical Society).
In Fig.
(9a), we show a comparison of the above result for steady state temperature profile with those obtained from direct simulations of the microscopic model, and we see very good agreement.It is interesting to note that the temperature profile is non-symmetric under space reversal as the microscopic model itself does not have such symmetry.This fact is also reflected in hydrodynamics where this shows in the existence of a single sound mode.
For the relaxation part we look for solutions which satisfy the initial condition T r (z, 0), C r (u, z, τ = 0) = 0 and boundary conditions C r (u, z, τ )| u→∞ = 0, T r (0, τ ) = T r (1, τ ) = 0.The solution of the "anti-diffusion" Eq. ( 55), with z as time variable, with the boundary condition (56) can be obtained as [34] C r (u, z, τ Using this in (56) then gives finally the evolution equation for the temperature field inside the domain 0 ≤ z ≤ 1.This is a non-local equation which can be recognized as a continuity equation where the current j is precisely in the form stated in Eq. ( 3).This is the open-system analogue of Eq. ( 51).For the infinite system, a similar computation leads to Eq. (68) but with the lower limit of integration replaced by z = −∞ and, by taking Fourier transforms, this can be shown to reduce to Eq. ( 52).We now proceed to solve Eq. ( 68) to find the temperature evolution.It is natural to expand the temperature profile T r (z, τ ) in a basis set satisfying Dirichlet boundary conditions, and we choose the set α n (z) = √ 2 sin(nπz), n = 1, 2, 3.... Substituting T r (z, τ ) = n Tn (τ )α n (z) in Eq. ( 68), we get Further we expand the function where L v nl = (nπ)ζ nl and the column vector | T has elements Tn = α n | T .The above equation is an infinitedimensional matrix representation of the non-local Eq. ( 68).To solve this, we diagonalize the matrix L v as R −1 L v R = µ, which gives the time dependent solution as | T (τ ) = Re κµτ R −1 | T (0) where R n,l = α n |ψ l denotes the nth element of the l-th right-eigenvector of the matrix L v and the diagonal matrix µ contains the corresponding eigenvalue µ l .The matrix L v is real but non-symmetric and it has left eigenvectors χ l | whose elements are given by χ l |α n = R −1 l,n .The formal solution for the temperature field T r (z, τ ) can then be written as where Finding the eigenspectrum of the matrix L v is a difficult problem as the matrix is infinite-dimensional and non-symmetric.However, one can truncate the matrix at some order and diagonalize it numerically, assuming that the spectrum converges with increasing truncation order.In [34] the authors used this apprach to compute the eigenspectrum and thereby study the time evolution of the temperature profile.This is shown in Fig. (9b).The spectrum is shown in Fig. (10) where it is seen that for large n, µ n ∼ π 2 |nπ| 3/2 (1 ± i) which is similar to the spectrum of the non-local operator L v in Eq. ( 52) describing the evolution in infinite system.In Fig. (11) we show the left and right eigenvectors χ n (z) = l=1 R −1 nl α l (z) and ψ n (z) = l=1 R ln α l (z) respectively, corresponding to the first eight eigenvalues.One observes that the eigenvectors corresponding to the first four eigenvalues are real whereas the eigenvectors corresponding to the eigenvalues with n > 4 are complex and come in conjugate pairs.

B. Harmonic chain with momentum exchange
In the previous section we discussed transport in the HCVE model which has two conserved quantities, namely volume and energy.In this section, we discuss heat transport in the harmonic chain momentum exchange (HCME) model which has three conserved quantities, namely volume, momentum and energy, that are the same as the ones in usual anharmonic chains with Hamiltonian dynamics [3,4].The model consists of a harmonic chain of particles each of unit mass and described by the degrees of freedom q i , p i , with i ∈ Z, corresponding respectively to position and momentum.As for the HCVE system, the dynamics of the HCME model also has two parts: (i) the usual deterministic part given by the Hamiltonian equations qi = p i , ṗi = ω 2 (q i+1 − 2q i + q i−1 ), i ∈ Z, where ω is the strength of the harmonic interaction and (ii) a stochastic part consisting of exchanges of momenta between neighboring particles (chosen at random) occurring with rate γ.In the absence of the stochastic exchange, the underlying Hamiltonian dynamics is integrable and the transport in this system is ballistic due to the absence of any scattering mechanism.The stochastic exchange introduces a momentum conserving scattering mechanism, which should make the transport behavior nonballistic.However, it turns out that the stochastic mixing is not sufficient to make the transport behavior diffusive.It has been shown rigorously that the energy current correlation in equilibrium of an infinite chain decays as t −1/2 similar to that in the HCVE model [9].This, through the closed system GK formula in Eq. ( 10), implies the anomalous system size scaling of the steady state current as j ∼ N −1/2 .
The HCME dynamics conserves the following three quantities: (a) total stretch i r i where r i = q i+1 − q i , (b) total momentum i p i and (c) the total energy i e i with e i = p 2 i /2 + ω 2 r 2 i /2.As a consequence, the corresponding local densities evolve slowly in the macroscopic length and time scales.In [29], it has been analytically shown that the local energy density e(x, t) in the isolated system evolves according to the following fractional diffusion equation and the fractional operator in Fourier space is represented by (−∆) 3/4 e ikx = |k| 3/2 e ikx .The NESS of this system was analyzed in detail in [26][27][28] where the scaling j ∼ N −1/2 and a closed form for the nonlinear temperature profile

Left-eigenvectors
Right-eigenvectors The real parts are indicated by blue lines while orange denotes the imaginary part.Note that the eigenvectors corresponding to the real eigenvalues (n = 1, 2, 3, 4) are also real where the eigenvectors corresponding to complex eigenvalues (n = 5, 6, 7, 8....) are complex.(Adapted from Kundu et.al.with permission from [34] Copyright (2018) by American Physical Society).
were established.More recently the fractional-equation-type description of this system in the open set-up was further discussed in [37].We summarize below some of these results for the open system.We first discuss the steady state and relaxation properties which is followed by the discussion on the evolution of the fluctuations and in the end we discuss the role of boundary conditions.

Typical behaviour of temperature, current and other correlations
In the open system HCME set-up, the two ends are attached to two reservoirs at temperatures T and T r .The equations of motion are now modified by adding Langevin forces to the 1st and the N th particles: +stochastic exchange of momenta at rate γ, for i = 1, 2 . . ., N , where ξ 1,N are independent Gaussian white noises with mean zero and unit variance, λ is the friction coefficient, and the parameter ζ has been introduced to describe different boundary conditions.Free boundary conditions correspond to ζ = 1 while fixed boundary conditions are given by ζ = 2.We will first discuss the fixed boundary case i.e. ζ = 2.We will be interested not only in NESS properties such as the form of the temperature profile and the current scaling with system size but also in the temporal evolution of the temperature from some arbitrary initial profile to the steady state form.As in the case of the HCVE model, the analytical tractability of the HCME system comes from the fact that the evolution of the two-point correlations is given by a closed set of equations.The two point correlations include U i,j = q i q j ,V i,j = p i p j , and Z i,j = q i p j and the local temperature defined as T i (t) = p 2 i consists of the diagonal elements of V.For these, one obtains a set of coupled linear equations, similar in form to Eqs. (54), which one needs to solve with appropriate boundary and initial conditions.The number of equations in this case is much larger than the HCVE case and hence it is even more difficult to solve them analytically for finite N .Observations from numerical solutions of these equations reveal [27] that for large N , the temperature field T i and the correlations z show the following scaling behaviors: . Hence, for large N it is instructive to construct solutions of these scaling forms.Inserting these scaling forms in the discrete equations of the two point correlations and taking the large N limit one finds, at leading order in 1/ √ N , the following partial differential equations [27] where the scaling variables We again note that for the isolated infinite system, one can follow the same procedure as above, but now replacing the scale parameter 1/N → a where a is the lattice spacing, to obtain the same set of Eqs.74-76 with ik Solving the first equation ( 77), with the condition that correlations vanish at u = ±∞, The equation ( 78) relates the constant A k to Tk : Using Eqs.(80,81) in Eq. ( 79) we get which is the Fourier representation of Eq. ( 72), with κ = ω 3/2 /2 √ 2γ.

(a)
< l a t e x i t s h a 1 _ b a s e 6 4 = " s O j Z U o N K 3 3 N c 8 n q m G B 8 y g q 4 L + + F q 0 5 J 5 s 5 h j 9 w P n 8 A i f m N T g = = < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 t  83) (solid black line) with the same obtained from direct numerical simulations of microscopic system for N = 128, 256, 512.The agreement between theory and numerics becomes better for larger N as can be seen in the inset where the difference between theoretical curve (Eq.( 83) and simulation data are plotted for various system sizes.(b) Time evolution of temperature starting from an initial step profile.The function Θ(v, τ ) = [T (v, τ ) − T ]/∆T , with T (v, τ ) given by Eq. ( 93), is plotted and compared with direct numerical simulations.The dashed lines indicate simulation results for the time-evolution in HCME at different scaled times (τ ), for system sizes N = 128 (red), N = 256 (blue), N = 512 (magenta), while the solid lines are obtained from the theory.The boundary temperatures were fixed at T = 1.5 and Tr = 0.5.(Adapted from Kundu et.al.with permission from [37] Copyright SISSA Medialab Srl, IOP Publishing).
In the steady state, the analytical solutions of these equations [with ∂ τ T (v, τ ) = 0] were obtained in [26] and are given by where T = (T + T r )/2, ∆T = T − T r and φ n (v) = δ n,0 + (1 − δ n,0 ) √ 2 cos(nπv) for n ≥ 0. From Eq. (76) we see that the current J = −ω 2 C(u, v, τ )| u→0 is given by Note that both the temperature profile and the current are independent of the friction coefficient λ.This is true only for the special case of fixed boundary conditions.Note also that the temperature profile in the steady state is intrinsically non-linear as can be seen in Fig. (12a) where one observes excellent agreement with data from simulations of the microscopic dynamics in Eq. ( 73).It can be shown that the temperature profile at both boundaries scales as ∼ (δv) µ with µ = 1/2 where δv is the distance from the boundary [26].This singular behavior of T ss (v) is a common signature of anomalous transport and it is characterized by the meniscus exponent µ.The value of µ however is non-universal and depends strongly on the boundary conditions.We will discuss this in Sec.(IV B 3).
To solve for the approach toward the above steady state results, we proceed as for the HCVE model.Separating the relaxation part and the steady state part we write Since the relaxation parts satisfy Dirichlet boundary conditions C r (u, 0, τ ) = C r (u, 1, τ ) = 0 and T r (0, τ ) = T r (1, τ ) = 0, we expand them in the Dirichlet basis α n (v) = √ 2 sin(nπv) for n = 1, 2, 3, ... as After inserting these expansions in Eqs.(74)(75)(76) and using the orthogonality property of the α n (v) functions, one gets the following (infinite order) matrix equation for the evolution of the components Tn : L p nl Tl , n = 1, 2, . . ., ∞, where: with S nl = α n |φ l = 1 0 dzα n (z)φ l (z), Λ nl = λ n δ nl is a diagonal matrix with λ n = (nπ) 2 and the constant κ = ω 3/2 /(2 √ 2γ).In the position basis, the above equation can be written as where the operator L p is represented as From this representation one can identify the action of L p on the set of basis functions φ m (which satisfy Neumann boundary conditions) [4,37] For the time evolution we need the eigenspectrum of L p with Dirichlet boundary conditions.The eigenstates ψ n (y) and eigenvalues µ n can be obtained by diagonalizing the matrix L p nm in Eq.(90).In [27] the spectrum was obtained numerically by diagonalizing truncated form of the infinite-dimensional matrix L p .An alternate method was recently proposed in [37] which gives the spectrum directly as roots of a transcendental equation and explicit series form expressions for the wave functions in the φ n basis.The numerical values of the computed eigenvalues are plotted in Fig. (13a), where we see that for large n the eigenvalues scale as µ n ≈ (nπ) 3/2 .At smaller values n there is a systematic deviation from the Neumann spectrum, λ n , for example the first three eigenvalues (µ n ) are given by µ 1 ≈ 2.75, µ 2 ≈ 12.02, µ 3 ≈ 24.22.As shown in the inset of Fig. (13a) the relative difference between µ n and λ n decreases as 1/n.The first few numerically computed eigenvectors are shown in Fig. (13b) where they are compared with the basis functions α n which are the Dirichlet eigenfunctions of the usual Laplacian.We observe that they are different and in particular show a non-analytic behavior at the boundaries.For example near the boundaries one finds ψ n (δv) ∼ √ δv, where δv is the distance from the boundaries.The eigenspectrum of fractional operator in a bounded domain, with different boundary conditions, has been discussed earlier in the literature, using somewhat heuristic approaches [75,[80][81][82].However their connection to the spectrum of L p defined here is unclear.
Using these Dirichlet eigenvalues and eigenfunctions, we follow the steps leading to Eq. ( 71) and obtain the following for the time evolution of an arbitrary initial profile: In Fig. (12b), a numerical verification of the above time evolution is shown.We note that Eq. ( 91) can be cast in the form of a continuity equation ∂ τ T r (v, τ ) = −∂ v j(v, τ ) with j in the form [37] j(v, τ ) = −κ where the kernel K is defined through it's action on a test function The operator L p can be expressed in terms of K as

Characterization of fluctuations
The discussions till now describe only the average or typical behavior of the conserved density fields and the associated current fields.The equation (91) describes the evolution of the average temperature profile as well as the evolution of a localized energy pulse in a thermally equilibrated system.However, other interesting aspects that characterize the NESS are the distributions of density and current fluctuations, long range correlations and the large deviations.To study these aspects, one requires to have a stochastic description of the evolution at the macroscopic length and time scales.
In the context of diffusive transport, a general framework called the macroscopic fluctuation theory has been developed in the last decade which allows to provide such a description for fluctuations [83][84][85].Starting from the microscopic description of the system one can show that in the diffusive scaling limit, the fluctuating energy density field e(x, t) and the corresponding fluctuating current J e (x, t) still satisfy the continuity equation but now, in addition to the regular diffusive part of the current, there is a fluctuating part J e (x, t) = −D(e) ∂e(x,t) ∂x + χ(e) η(x, t), where χ(e(x, t)) is the mobility of the system which is related to the diffusivity D(e(x, t)) through the fluctuation dissipation relation and η(x, t) is a mean zero white Gaussian noise with the properties η(x, t) = 0 and η(x, t)η(x , t ) = δ(x − x )δ(t − t ).The evolution equation for the energy density is given by ∂e Starting from this stochastic equation one can compute various moments, fluctuations and correlations of e(x, t) and j(x, t) both in stationary and non-stationary regime.This description also allows one to compute the probabilities of observing atypical density and current profiles which are characterized by large deviation functions.The whole program has been established and applied in several microscopic systems which show diffusive behavior at macroscopic scales.We ask if a similar procedure works for our system, displaying anomalous transport, and described by the fractional diffusion equation.Recently such an extension has been proposed in [37] which we now describe.The approach in [37] is to include a noise part in the current expression in such a way that the fluctuation-dissipation theorem is satisfied.For a system in equilibrium at temperature T this leads to the unique choice where η(v, τ ) is white Gaussian noise with η(v, τ ) = 0, η(v, τ )η(v , τ ) = δ(v − v )δ(τ − τ ) and the fluctuationdissipation theorem implies the relation with B † defined as the adjoint of B. It was verified in [37] that Eq. ( 98) reproduces correctly results on energy correlations and current fluctuations in equilibrium.Extending this approach to the non-equilibrium situation was also attempted in [37] and a conjecture for long-range correlations in the NESS was proposed.

Role of boundary conditions: Hydrodynamic theory
In the previous section we have mainly discussed the fixed boundary condition, in which case we have learned that the transport behavior in HCME model is anomalous with exponent α = 1/2 and the Fourier's law gets modified to a non-local linear response relation as in the form of Eq. ( 94) with an explicit form for the kernel K(v, v ) given in Eq. ( 95).Also in this case the evolution of the temperature profile is given by a non-local equation (91) with L p defined through Eqs. ( 95) and (96).In this section we would like to understand the dependence of these results on the choice of boundary conditions.In particular we focus on the case of free boundary conditions i.e. for ζ = 1 in Eq. (73).
Energy transport in HCME with free boundary condition was studied numerically in [28] where it was observed that the system size scaling of the current j in the steady state is again proportional to 1/ √ N , as for fixed BC.However, in contrast to the fixed BC case, the proportionality constant depends on the friction coefficient λ.It was also observed that the temperature profile in this case is non-linear but the associated meniscus exponent µ depends strongly on the relative values of λ and ω.For this case finding the appropriate boundary conditions for Eqs.(74,75,76) is a difficult problem [28] and has so far not been possible.A different approach, based on linear response theory and NFH was proposed in [36] and we present some details here.
This approach starts with the following non-local linear response result which is based on a linear response calculation as done in [67] but around a local equilibrium state characterized by a temperature profile.According to this calculation the Kernel is related to the equilibrium current-current correlation [36] K where j(x, t) is the local current and a is a constant.For systems with AHT we expect N j(x, t)j(y, 0) eq ∼ t 1−α which means that K N (x, y) should scale as N α−1 .Hence we expect that the limit exists, which implies also that j = J/N 1−α with J given by where the temperature profile T (x) is assumed to have the scaling form T (x) = T + ∆T Θ(x/N ).This equation can then be used to compute the NESS temperature profile and also the current.The remaining task now is to compute the kernel K(v, v ).
For HCME, the kernel K(u, v) has recently been computed in [36] using the techniques of NFH as introduced in Sec.(III B).Following this procedure for the HCME model, one finds that on hydrodynamic length and time scales, a random fluctuation created inside the system decomposes into two ballistically moving but diffusively spreading sound modes φ ± and a stationary heat mode φ 0 .In terms of the local stretch r i = q i+1 − q i and energy e i = p 2 i /2 + ω 2 r 2 i /2, the sound modes and the heat mode are expressed as φ ± = ωr ∓ p and φ 0 = e, respectively.The evolution of these modes are given by [4] ∂ where c s = ω is the speed of sound, η + , η 0 and η − are uncorrelated Gaussian white noises, G = ω 4 and D and D 0 are phenomenological diffusion coefficients.
The instantaneous energy current can be read from (104), neglecting the sub-dominant terms arising from the momentum exchange and the noises η ± [62].The stochastic momentum exchange process generate a diffusive contribution [see Eq. ( 104)] which becomes sub-leading at large N and the noises η ± also do not contribute since their time averages vanish.
In order to compute the the kernel in (101) using the form of j(x, t) in (105), one needs to solve the equations of φ ± in (104) inside a finite domain with suitable BCs.At this point we would like to mention that originally the NFH theory was formulated for an infinite domain [62].The work in [36] provides an extension to incorporate boundary conditions for a finite domain, in the context of the HCME model.As the equations for φ + and φ − are independent of φ 0 , it is straightforward to write the solution in terms of the appropriate Green's function, as shown later.
We now discuss how to get the boundary conditions of fields φ ± .The strategy that has been followed in [36] is to introduce extra stretch and momentum variables in such a way that the equations at the boundary points (i = 1, N ) are also included into the structure of the bulk equations.This can be acheived by introducing additional conditions, which after appropriate coarse-graining become the hydrodynamic BCs .To explain the procedure let us consider the free BC case as an example.We first introduce an extra dynamical variable r 0 in such a way that the form of the equation satisfied by p 1 becomes same as that of the bulk evolution equations with the condition where we have neglected the noise terms in (73).This provides one BC.We need another BC as the Eq. ( 104) is of second order in space.As before, introducing p 0 in such a way that one can make r 0 to satisfy a regular equation of motion in the bulk at the cost of an extra condition, provides the second BC.Taking single derivative with respect to time on both sides of the first condition yields One can get two other boundary conditions by applying similar procedure to the equations of the last (N th) particle.Finally, coarse-graining over space and expressing the stretch r and momenta p in terms of the sound modes φ ± , we obtain the following BCs for free boundaries: where These BCs can be interpreted physically as some sort of partially 'reflecting' boundaries.The BCs on the first (second) line of Eq. (108) mean that when a φ + (resp.φ − ) Gaussian peak hits the right (resp.left) boundary, it gets reflected as a φ − (resp.φ + ) Gaussian peak with area under the peak reduced by a factor w.This feature has been observed in numerical simulations and the validity of (108) has been confirmed [36].There are two interesting cases w = 0 and w → 1.In case of resonance (also called impedance matching) λ = ω i.e. w = 0 [66], once a φ ± peak hits the boundary nothing gets reflected because everything gets absorbed at the boundary reservoirs.On the other hand, w → 1 corresponds to almost perfectly reflecting case.This situation arises for the fixed BCs in the microscopic dynamics.Following a similar procedure as done for free BCs, it is possible to show that one arrives at the same hydrodynamic BCs Eq. ( 108) except now w = 1.From Eq. (109), one can easily see that the w → 1 limit is achieved for λ → ∞.In this limit, the 1st and the N th particles hardly move i.e. their positions q 1 and q N stay very close to 0 for all times due to infinite dissipation and therefore mimic the fixed BCs for the microscopic dynamics.So for fixed BCs we have the hydrodynamic BCs Eq. ( 108) with w = 1.Since the hydrodynamic equations (104) for φ + and φ − along with the BCs (108) are linear, it is easy to solve them for arbitrary initial condition.The solutions are best expressed in terms of the four Green's functions f σ,τ (x, y, t) for σ, τ = ±, as where, with w = 1 for fixed BCs and w = λ−ω λ+ω for free BCs.Using this expression in Eq. ( 105) one finally gets from Eqs. (101,102) the following expression for the kernel: where the constant A = G 2 S 2 T 2 √ Dcs with S = φ + (x, 0) 2 eq = φ − (x, 0) 2 eq = 2 T and R = w 2 .The diffusion constant D appearing in the equation for φ ± arises from the exchange mechanism and it can be shown from a microscopic calculation that D = γ/2.This then gives A = ω 3/2 /(2 √ 2γ) which we note coincides with the expression for κ in Eq. ( 72), and so we identify A = κ.One can use this kernel in Eq. ( 103) to compute the current and the temperature profile Θ(v).
Let us define the Greens function, G R , corresponding to the kernel K R through the equation Then Eq.( 103) can be inverted to give Solving this equation with the boundary conditions Θ(0) = 1/2, Θ(1) = −1/2 gives us the expressions for the current and temperature profile One uses this in Eq. ( 103) to solve for the temperature profile Θ(v).The above analysis, based on linear response calculation, assumes |∆T | << T .However for HCME, one observes that the quadratic correlations satisfy a closed set of linear equations with a source term proportional to ∆T [26].Hence the temperature profile Θ(v) in ( 116) is also valid for any ∆T .It turns out the equation ( 114) can be solved analytically and exact expressions of the temperature profile Θ(v) can be obtained in the following two limiting cases -(i) Free resonant case R = 0: In this case the kernel is simply given K 0 = 1/ 2π |v − v | which is same as that of an infinite system.For this kernel, the solution of Eq. ( 103) can be directly written using standard results on solution of integral equations [86] as This can be solved with the boundary conditions to give the temperature profile This profile is verified numerically in Fig. 14 (left panel), where we observe diverging derivatives at the boundaries.
From the above expression it is possible to show that the meniscus exponent is µ = 3/4.(ii) Perfectly reflecting case R → 1: As mentioned above this is equivalent to fixed BC for which the temperature profile, given in Eq. ( 83), was computed from microscopic calculation in the previous section.In this case it is known [37] that the eigenfunctions of the operator K R are precisely the sine-functions α n (v), i.e which is consistent with Eq. ( 95).This then gives us the corresponding Green's function Using this and Eqs.(114,115) we recover the exact expressions for the temperature profile and current given in Eqs.(83,85) [37].
For free BCs with λ = ω we have 0 < R < 1.In this case it is difficult to solve Eqs.(103,112) analytically but numerical solutions have been obtained.In Fig. 14 (right panel) a comparison of the temperature profile obtained from the numerical solution and from direct microscopic simulations for R = 1/2 and one can observe excellent agreement.Note again that the temperature profile is singular at the boundaries.It turns out that the exponent µ characterizing this singularity depends on not only on α but also on R [66].To determine this dependence we take a derivative with respect to v of Eq. ( 103) and get Although the integral is identically zero for all v, the individual terms in the integrand have divergences.For example, the kernel diverges as K R ∼ |v − v | −1/2 while system, beyond which time the correlations decay exponentially.The scaling function is given by the Lévy-stable distribution in the bulk and the finite cut-off leads to the width of the pulse scaling as σ(t) ∼ t β/2 .
In this review we discussed these signatures of anomalous transport and showed how they can be understood within three different but related frameworks -(a) a phenomenological model where the heat carriers are taken to be Lévy walkers, (b) a microscopic phenomenological approach based on nonlinear fluctuating hydrodynamics and (c) exact results obtained for certain stochastic models.The main picture that emerges is that anomalous heat transport can be understood by replacing Fourier's law in Eq. ( 1) by a non-local fractional-type diffusion equation given in Eq. ( 3), where the precise form of the kernel K R (x, y) depends on the specific set-up and boundary conditions.For the stochastic models the form of the kernel is known explicitly both for the closed system (infinite line) and the open system.In the Lévy walk picture, where the distribution of flight times has a power-law dependence ∼ 1/t ν+1 , the kernel has the asymptotic form K R (x, y) ∼ 1/|x − y| ν−1 .We saw from the various approaches, that all the different exponents mentioned above are related to each other and in fact can be expressed in terms of the Lévy walk exponent as For the Hamiltonian models that we discussed, namely the alternate mass hard-particle gas and the FPUT model, the various exponents are given by α = 1/3, β = 4/3, γ = 3/5 and correspond to a Lévy-walk exponent ν = 5/3.For the stochastic momentum exchange model we have α = 1/2, β = 3/2, γ = 2/3 which corresponds to ν = 3/2.The meniscus exponent µ is non-universal and depends on ν and on boundary conditions through a single dimensionless number R, which can be interpreted as the reflection coefficient of the Lévy walkers at the boundaries.In the context of the exactly solvable stochastic models, we discussed the spectrum of the fractional-type Laplacian operator [specified by the kernel K R (x, y)] in the open set-up, and pointed out important differences with the spectrum of the usual Laplacian for diffusive processes.We conclude by mentioning some outstanding open questions in the field.
• Hamiltonian systems -The Lévy walk behaviour has been clearly observed in large number of simulations.The formalism of NFH gives a microscopic justification of the Lévy walk model and the fractional-diffusion type description of the heat mode.Some open questions include: 1.A more rigorous microscopic derivation of the evolution equation of a localized heat pulse in an equilibrium system, to show that the central peak satisfies a fractional-diffusion type equation of a form similar to that in Eq. (3).
2. Extension of the NFH formalism to the non-equilibrium case to study transport in finite open system and understand the role of BCs.Detailed simulations are also required to understand the effect of BCs.
3. Establishing the Lévy walk picture from a microscopic viewpoint ?
• Stochastic systems -For the HCME model, the non-local version of Fourier's law has been established and the response kernel K R computed so far using two methods: (i) exact microscopic method for the BC corresponding to R = 1 and (ii) using NFH for arbitrary R. Is it possible to extend the exact microscopic approach to find the non-local kernel K R for general boundary conditions.Similarly for the HCVE it would be interesting to explore the role of BCs.
• For the HCME model, it has been possible to find the eigenspectrum of the non-local kernel K R for the case R = 1 and it was observed that the eigenvalues for Dirichlet and Neumann boundary conditions differ (unlike for the usual Laplacian).Finding the spectrum of the non-local kernel K R for general R, for Dirichlet and Neumann boundary conditions, is an interesting mathematical problem.The knowledge of the spectrum, namely eigenvectors and eigenvalues, enables one to study the time-evolution.
• For the HCME model we showed that it is possible to write a stochastic non-local equation [Eq.(98)] to describe equilibrium fluctuations.An open problem is to write such an equation in the non-equilibrium set-up.For diffusive systems this is given by Eq. (97) and this equation enables one to compute long-range correlations in the NESS and large deviation functions.

FIG. 1 .
FIG.1.Schematic illustration of the (a) closed system set-up and the (b) open system set-up, commonly used to study heat transport.In (a), a localized heat pulse is introduced at some point in a system in thermal equilibrium and its subsequent time-evolution is observed.In (b), the system is attached to two heat reservoirs at different temperatures and the NESS properties such as current and temperature profile are studied.

FIG. 3 .
FIG. 3. Scaled perturbation profiles at times t = 40, 80, 160, 320, 640, 1280, 2560, and 3840 , with γ = 3/5.The profiles have been obtained by averaging over large number of realizations.In the inset, the profile at t = 640 (solid line) is compared with the propagators of a Lévy walk with an exponent ν = 5/3 with a fixed velocity v = 1 (dotted line) or with velocity chosen from a Gaussian distribution with mean 1 and variance 0.036 (dashed line).(Adapted from Cipriani et.al.with permission from [8] Copyright (2005) by American Physical Society)

FIG. 3 (
FIG. 3 (color online).Heat current as a function of N. Root mean square errors from O10 5 ÿ 10 6 measurements are shown except when they are smaller than the points.(Top) Conductivity versus N. The last five points fit to a slope of 0:333 0:004.The baths are described in Fig.2.(Bottom) jN 1ÿ versus N for 1=3 and 2=5.In the large N regime, is definitely less than 2=5 and appears to agree quite well with the 1=3 prediction for all data sets.Langevin baths with 0:4, 2, and 10, and one data set with Nose-Hoover baths, are shown.

FIG. 6 .
FIG.6.Plot of the scaled distribution t 2/3 P (x, t) versus x/t 2/3 of the Lévy walk on the open line for ν = 3/2 at three different times.Also shown is a plot of the Lévy-stable distribution.The inset shows a plot of the mean square displacement and the fourth moment and a comparison with the exact asymptotic forms (dashed lines) given by Eqs.(27,28).In all plots the time to and c are set to one.

FIG. 7 .
FIG. 7. (a) Temperature profile of the oscillator chain with conservative noise with free boundary condition and λ = γ = 1 (solid line) and solution of the master equation with reflection coefficient R = −0.1 (dashed line).(Reprinted from[66]).(b) Dependence of the meniscus exponent µ on the reflection probability r for ν = 3/2.Note that r in this figure is denoted by R in this review.Full circles and stars are measures from fitting of P (x) and Q(x), respectively (see text).The dashed line is given by formula(29).(Adapted from Lepri et.al.with permission from[66] Copyright (2011) by American Physical Society).
The dynamics has two parts: (a) a deterministic part given by dηi dt = η i+1 − η i−1 and (b)

2 FIG. 10 .
FIG.10.The real and imaginary part of the alternate eigenvalues for the matrix L v in Eq.(70).The first 4 eigenvalues are completely real and distinct.The higher eigenvalue comes in complex conjugate pairs.For large µn ∼ (nπ) 3/2 (1 ± i).For smaller n, there is a deviation from asymptotic scaling due to finite definition of the operator.(Adapted from Kundu et.al.with permission from[34] Copyright (2018) by American Physical Society).
FIG.12.(a) Comparison of temperature profiles obtained theoretically from Eq. (83) (solid black line) with the same obtained from direct numerical simulations of microscopic system for N = 128, 256, 512.The agreement between theory and numerics becomes better for larger N as can be seen in the inset where the difference between theoretical curve (Eq.(83) and simulation data are plotted for various system sizes.(b) Time evolution of temperature starting from an initial step profile.The function Θ(v, τ ) = [T (v, τ ) − T ]/∆T , with T (v, τ ) given by Eq. (93), is plotted and compared with direct numerical simulations.The dashed lines indicate simulation results for the time-evolution in HCME at different scaled times (τ ), for system sizes N = 128 (red), N = 256 (blue), N = 512 (magenta), while the solid lines are obtained from the theory.The boundary temperatures were fixed at T = 1.5 and Tr = 0.5.(Adapted from Kundu et.al.with permission from[37] Copyright SISSA Medialab Srl, IOP Publishing).

FIG. 14 .
FIG. 14. Rescaled temperature profile for resonant BCs R = 0 (left panel) and free BCs with R = 1 2 (right panel).In the main plots results of Monte-Carlo simulations for increasing system sizes N = 100, 200 and 400 are compared to the theoretical predictions given by Eq. (118) for the R = 0 case (left panel) and the numerical solution of Eq. (103,112) with R = 1 2 for the plots in the right panel.In the insets the differences between measurements and theory are shown.The other parameter values are T+ = 1.5, T− = 0.5, and ω = γ = 1.(Adapted from Cividini et.al.with permission from [36] Copyright SISSA Medialab Srl, IOP Publishing.) To compute the MSD and relate the exponents β and γ is a somewhat subtle issue and requires one to note that the scaling function is valid in the bulk region |x| < ∼ t, beyond which (x, t) decays rapidly [see discussion in Sec.(III A 1) in the context of Lévy-walk model].From properties of [37] 13.(a)Eigenvalues of the fractional operator in Eq. (90) corresponding to Dirichlet boundary conditions.For large n, the slope is seen to approach that of n 3/2 (black dashed line).For small n there is a systematic difference between the Dirichlet and Neumann eigenvalue snd the inset plots the difference between the two.For large n the difference between the two decreases inversely with n.(b)The first six eigenvectors, ψn(v) (black lines), are compared to the corresponding Dirichlet eigenfunctions of Laplacian i.e. sin-functions (blue dashed lines).The eigenstates are different from sin-functions, especially near the boundaries, even for large n.(Adapted from Kundu et.al.withpermissionfrom[37]Copyright SISSA Medialab Srl, IOP Publishing.) s