Beyond phonon hydrodynamics: Nonlocal phonon heat transport from spatial fractional-order Boltzmann transport equation

Spatially convoluting formulations have been used to describe nonlocal thermal transport, yet there is no related investigation at the microscopic level such as the Boltzmann transport theory. The spatial fractional-order Boltzmann transport equations (BTEs) are first applied to the description of nonlocal phonon heat transport. Constitutive and continuity equations are derived, and two anomalous behaviors are thereafter observed in one-dimensional steady-state heat conduction: one is the power-law length-dependence of the effective thermal conductivity, κeff ∝ Lβ with L as the system length, and the other is the nonlinear temperature profile, [T(x) − T(x = 0)] ∼ x1+η. A connection between the length-dependence and nonlinearity exponents is established, namely, β = −η. Furthermore, we show that the order of these BTEs should be restricted by the ballistic limit. In minimizing problems, the nonlocal models in this work give rise to different results from the case of Fourier heat conduction, namely that the optimized temperature gradient is not uniform. © 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0021058., s


I. INTRODUCTION
In the past decades, temporal fractional-order derivatives have been applied to various physical processes and systems, i.e., anomalous diffusion 1-3 and non-Fourier heat conduction. [4][5][6][7] In recent years, temporal fractional-order modeling of heat conduction attracts increasing interest. 8,9 Generally speaking, heat conduction in three-dimensional bulk materials obeys Fourier's law, In Eq. (1), q = q(x, t) denotes the heat flux, T = T(x, t) is the local temperature, and κ is the thermal conductivity. One unsatisfactory feature of Fourier's law is that it predicts an infinite speed of heat propagation, [10][11][12][13] which is traceable to the parabolic governing equation of the local temperature. The most well-known modification to overcome this unphysical behavior is the Cattaneo model, 14 which adds a relaxation correction to Eq. (1), namely, with τ as the relaxation time. Combining the Cattaneo model with the standard continuity equation yields a hyperbolic governing equation as follows: with c standing for the specific heat capacity per volume. Equation (3) guarantees a finite speed √ κ/cτ yet is still paired with other unsatisfactory problems [15][16][17][18] such as negative temperature and entropy generation.
Many generalizations have been proposed for the Cattaneo model. The Jeffery model 10 is a typical one, which assumes a coexistence of Fourier and hyperbolic transport regimes of heat conduction. Another common example is the phase-lagging class, 19 i.e., q(x, t + τ) = −κ∇T(x, t). (4)

ARTICLE scitation.org/journal/adv
Cattaneo equation (GCE) by Compte and Metzler,20 which proposed a class of fractional-order generalizations. In this class, temporal fractional-order derivatives are introduced into the constitutive and continuity equations, i.e., In Eq. (5), RL 0 D α t is the temporal Riemann-Liouville operator, namely, wherein n ∈ N ∩ (α, α + 1) and q(x, t) is at least n-order differentiable. The corresponding long-time and short-time asymptotic behaviors of the mean-square displacement (MSD) were thereafter investigated. The above non-Fourier models exhibit memory effects between the heat flux and the temperature gradient, which can be formulated as a convolution as follows: 10 In Eq. (7), the memory kernel K(t − t ′ ) is usually termed the relaxation function, which is an exponential function for the Cattaneo model. For the phase-lagging class, K(t − t ′ ) obeys the Dirac delta function and will become the Mittag-Leffler functions for the fractional-order generalizations.
The formulation in Eq. (7) aims to describe non-Fourier effects when the characteristic time of a thermal transport process is comparable to or even much smaller than the relaxation time. Another natural non-Fourier aspect is that the characteristic size of a thermal transport process is comparable to or even much smaller than the mean free path. Such non-Fourier effects cannot be expected by Eq. (7) because it will reduce to Fourier's law in stationary thermal transport. Therefore, convolution constitutive models including the fractional-order type shall be proposed for the spatial variable as well. For instance, Sapora and co-authors 24 have generalized the Caputo operator and applied it to one-dimensional (1D) steady-state heat conduction in 0 ≤ x ≤ L, namely, with cα as a material constant. A similar nonlocal constitutive relation has been proposed for the silicon thin film 25 as well. Such convolution formulation is from an analogy with the existing elastic theory, namely that the stress at a point x = x 0 in an elastic body depends on not only the strain at x = x 0 but also strains at all other points. Existing studies on the fractional-order nonlocality of thermal transport mainly focus on the macroscopic constitutive models and lack discussion on the underlying microscopic regimes such as the Boltzmann transport theory, which will be addressed in this paper. We mention that the nonlocality from the integer-order Boltzmann transport theory has been studied, [26][27][28][29][30][31] and the nonlocal effects are commonly induced by the second-order derivatives of the heat flux such as (∇ 2 q) and (∇∇ ⋅ q). This class is usually termed phonon hydrodynamics and has been employed to explain non-Fourier behaviors observed in experiments. [32][33][34] Obviously, such nonlocal terms will vanish in 1D heat conduction without the heat source, which is fundamentally different for the fractional-order nonlocal forms such as Eq. (8).
In this work, we consider the phonon thermal transport regime in terms of the phonon Boltzmann transport equation (BTE), i.e., the single mode relaxation time 35,36 and Callaway 37 models. These BTEbased models correspond to a localized physical picture, namely that the derivatives of the distribution function at x = x 0 are only determined by the collision term at x = x 0 . In this work, nonlocal regimes of phonon transport are first investigated via generalizing the phonon BTE into fractional-order forms, wherein the derivatives at x = x 0 depend on the global evolution in 0 ≤ x ≤ L. Based on the generalized BTEs, macroscopic constitutive and continuity equations are then established, which can give rise to anomalous signatures. According to a sound review by Lepri et al., 38 the anomalies in heat conduction can be classified as five signatures including the length-dependent effective thermal conductivity, slow decay of the heat current correlation function, anomalous MSD, nonlinearity of stationary temperature profiles, and fast relaxation of spontaneous fluctuations. In the present work, two anomalous signatures are observed. One is the length-dependent effective thermal conductivity, [39][40][41][42][43] Typical examples 40 are the billiard gas channel models and Fermi-Pasta-Ulam (FPU) chain. For two-dimensional (2D) systems, the lengthdependence commonly becomes logarithmic, κ eff (L) ∼ log L, i.e., 2D momentum conserving systems 44 (we mention that there also exists a numerical calculation 45 that supports the power-law lengthdependence for 2D cases). The other is the nonlinearity of stationary temperature profiles, which is asymptotic to a power-law function at the boundary, Therefore, we provide a BTE-based framework for modeling nonlocal thermal transport with the two anomalies. The two non-Fourier behaviors cannot be expected by previous nonlocal theories such as phonon hydrodynamics, which will generally reproduce Fourier heat conduction in the large size limit. It indicates that the fractionalorder nonlocal theory can extend the existing framework of phonon thermal transport. Furthermore, the fractional-order nonlocal type is able to describe non-Fourier regimes beyond phonon hydrodynamics, i.e., anomalous heat diffusion following the Levy walk. 38 Furthermore, a connection between the divergence and nonlinearity exponents will be established. We mention that our previous papers [46][47][48] have discussed the temporal fractional-order BTEs, which correspond to the relaxation or memory effect rather than the nonlocality. Here, we propose the spatial fractional-order BTEs for the nonlocal effects, which predict non-Fourier anomalies different from the temporal fractional-order situations, i.e., the nonlinear temperature profiles. Another difference from our previous works is ARTICLE scitation.org/journal/adv related to the minimizing problems. The temporal fractional-order models expect the same stationary solution as Fourier's law, and hence, the corresponding optimized results will remain unchanged. In contrast, the nonlocal heat conduction in this work will give rise to non-Fourier optimized results, namely that the optimized temperature gradient is not uniform. These are the main original aspects of this paper.

II. PHONON BOLTZMANN TRANSPORT EQUATIONS AND NONLOCAL GENERALIZATIONS
Usual phonon transport is described by the following BTE with the single mode relaxation time approximation: where h is the reduced Planck constant, ω is the angle frequency, kB is the Boltzmann constant, and vg is the phonon group velocity. For phonon heat transport, the local energy density is given by while the heat flux reads Based on the above expressions, one can obtain the standard energy continuity equation by multiplying Eq. (11) by ̵ hω and integrating it over the wave vector space, namely, When we multiply Eq. (11) by vg ̵ hω and integrate it over the wave vector space, the standard Cattaneo model emerges. It should be mentioned that the standard Cattaneo model is able to be derived from Eq. (11) only in the presence of local equilibrium, which allows the assumption ∇f ≈ ∇f 0 . This condition will be invalidated in the situations that the characteristic size and time are comparable to the mean free path (MFP) l = |vg|τ and relaxation time, respectively. Then, the solution of Eq. (11) will give rise to derivations from the standard Cattaneo model. A typical example is Chen's solution termed ballistic-diffusive heat conduction equations. 36 In Chen's method, the distribution function is divided into two parts, f = fm + f b . fm corresponds to the diffusive transport, which corresponds to the standard Cattaneo model. fm reflects the ballistic contribution beyond the Cattaneo model and leads to the socalled ballistic heat flux. Although the single mode relaxation time approximation is widely applied to heat conduction, its prediction underestimates the thermal conductivity of two-dimensional materials, 49 which are traceable to the effects of the normal scattering. Thus, more refined modifications were applied, i.e., Callaway's dual relaxation approximation, 37 where τR is the relaxation time of resistive process, τN is the relaxation time of normal process, fD = 1 exp( is the displaced Planck distribution, and u is the phonon drift velocity determined by the quasi-momentum conservation. Equations (11) and (15) are commonly used BTEs for phonon heat transport. In this work, we discuss their spatial fractional-order generalizations, namely, where μ ∈ R and D μ x is the spatial fractional-order operator. The above BTEs indicate nonlocal effects in phonon transport, namely that the evolution of the distribution function at x = x 0 depends on the drift and collision terms at not only x = x 0 but also x ≠ x 0 . The corresponding energy continuity equations can be obtained upon multiplying the above BTEs by ̵ hω and integrating them over the wave vector space. When the fractional-order derivative occurs in the collision terms such as Eqs. (16a) and (16c), the energy continuity equation remains Eq. (14). For Eqs. (16b) and (16d), the corresponding energy continuity equation deviates from the standard form, namely, Upon multiplying Eqs. (16a) and (16b) by vg ̵ hω and integrating them over the wave vector space, we can derive the following constitutive equations, respectively: Different from the single mode relaxation time approximation, there are various macroscopic constitutive equations for the Callaway model and its fractional-order generalizations. If we use the eigenvalue analysis method or sub-first-order approximation 34 to solve Eqs. (16c) and (16d), the corresponding constitutive equations are, respectively, given by which will reduce to the well-known Guyer-Krumhansl type as μ = 0. For 1D steady-state problems, we have dq dx = 0. Then,

ARTICLE scitation.org/journal/adv
Eqs. (19a) and (19b) will reduce to Eqs. (18a) and (18b), respectively. Some equations in this section possess similar forms to our previous works. Despite the mathematical similarity, the non-Fourier regime of this work is fundamentally different because the fractional-order operators are here for the spatial variable rather than the temporal variable. The spatial fractional-order operators are parried with the nonlocal effects, while the temporal fractional-order operators yield the memory effects.

III. RESTRICTIONS AND ANOMALIES
In this section, we will show that there exist several physical restrictions on the operator D μ x . We first select the spatial Riemann-Liouville operator as D μ x , with n ′ ∈ N ∩ (μ, μ + 1). In this case, the 1D stationary temperature distribution and heat flux predicted by Eq. (18a) are given by with ΔT = T(x = L) − T(x = 0). Equation (21) implies a power-law length-dependent effective thermal conductivity κ eff , namely, which agrees with the generic behavior of 1D systems. The critical situation is that μ = 1. This is because β = 1 corresponds to the ballistic transport, which is an upper bound of the length-dependence exponent β. Furthermore, μ = 1 will cause a singularity in Eq. (21), which requires dT dx ≡ 0. Thus, the order should be restricted as μ < 1 from either physical or mathematical viewpoint. According to usual understandings of anomalous heat conduction, 38 0 < β < 1 corresponds to the superdiffusive regime, while β < 0 emerges from subdiffusive heat conduction. Therefore, superdiffusive heat conduction can be modeled in terms of Eq. (16a) with 0 < μ < 1, and in the subdiffusive case, μ < 0. We can also observe a nonlinear stationary temperature profile in Eq. (21), namely, Combining β = μ and η = −μ yields a connection between the length-dependence and nonlinearity exponents, This connection means that the transport regime and temperature profile are not independent of each other in nonlocal phonon transport. Subdiffusive heat conduction (β < 0) will correspond to a superlinear temperature profile, while superdiffusive heat conduction (β > 0) must be paired with a sublinear temperature profile. For Eq. (16b), the solutions become Similar to the above case, there also exists a power-law lengthdependent effective thermal conductivity, The ballistic limit β = 1 leads to μ ≥ −1. The subdiffusive and superdiffusive regimes correspond to μ > 0 and −1 < μ < 0, respectively. The stationary temperature profile is nonlinear as well, and the relation between β and η remains Eq. (24). We now consider Eqs. (16a) and (16b) based on the Caputo operator as follows: Noting that the Caputo and Riemann-Liouville operators are identically equal as μ < 0, and hence, we only discuss the case of μ > 0. When μ > 0, there is no near-equilibrium solution for Eq. (16a) unless dT dx ≡ 0, which indicates that the Caputo operator is inapplicable. For Eq. (16b), the solution are still given by Eq. (25).

A. Symmetrization
The above fractional-order generalizations reflect asymmetric nonlocality, namely that phonon transport at x depends on the evolution of the distribution function in [0, x], but the evolution in (x, L] has no contribution. In order to acquire a symmetrical nonlocality, the above fractional-order operators can be symmetrized as follows: ARTICLE scitation.org/journal/adv For Eq. (29a), the stationary solutions are given by The effective thermal conductivity is still Eq. (22), and the nonlinearity of the temperature profile is symmetric, The stationary solutions for Eq. (29b) are written as The effective thermal conductivity remains as in Eq. (26), and the temperature profile possesses the following asymptotics: which is also symmetric.

B. Minimizing problem
In optimization design of heat transfer, minimizations of physical quantities have been widely used as optimizing strategies. [50][51][52] The entransy dissipation minimization 52 is a typical strategy. For 1D steady-state problems in 0 ≤ x ≤ L, the total entransy dissipation ratė G is given byĠ In the following, we will show the derivations between the nonlocal models in this work and Fourier's law in the entransy dissipation minimization. We consider a minimizing problem that is constrained by the boundary heat flux q(x = 0) = q 0 and total amount of the thermal conductivity, For Fourier heat conduction, minimization ofĠ will lead to a constant temperature gradient, Owning to the existence of fractional-order derivatives, the usual Lagrange multiplier approach is invalid for nonlocal heat transport. In order to minimize the entransy dissipation, we introduce the following Cauchy-Bunyakowsky-Schwarz inequality: (38) and thereupon,Ġ fulfills the following inequality: which becomes equality if and only if Different from the optimized result in Eq. (36), the optimized temperature gradient in Eq. (40) is not uniform. It implies that nonlocal heat transport based on the fractional-order derivatives will perform non-Fourier behaviors in minimizing problems as well. For Eq. (18b), the Cauchy-Bunyakowsky-Schwarz inequality yieldṡ G ≥ q 2 0 L 2+μ κ 0 l μ Γ(1 + μ) ( and the corresponding optimized temperature gradient satisfies which still differs from the optimized result in Fourier heat conduction.

V. CONCLUDING REMARKS
1. The spatial fractional-order derivatives are first introduced into the phonon BTEs, which reflects nonlocal effects in heat transport. The corresponding constitutive and continuity equations are thereupon derived, which also contain fractional-order derivatives. Based on the constitutive and continuity equations, the temperature distribution and heat flux are derived for 1D steady-state heat conduction. 2. Anomalies are thereafter observed, including the power-law length-dependence of the effective thermal conductivity and nonlinear temperature profile. We also establish a connection between the length-dependence κ eff ∝ L β and nonlinearity [T(x) − T(x = 0)] ∼ x 1+η , namely, β = −η. Furthermore, we ARTICLE scitation.org/journal/adv show that the order μ should be restricted due to the ballistic limit. 3. In order to reflect a symmetric nonlocality in phonon transport, the Caputo and Riemann-Liouville operators are symmetrized. The symmetrized modifications predict the same length-dependence as the original models. The corresponding temperature profiles are nonlinear as well, which are asymptotic to power-law functions at the two boundaries. 4. Minimizing problems in steady-state heat conduction are discussed. For the entransy dissipation minimization, the nonlocal models in this work will lead to different optimized results from the result of Fourier heat conduction, namely, the non-uniform optimized temperature gradient.