Reciprocity of thermal diffusion in time-modulated systems

The reciprocity principle governs the symmetry in transmission of electromagnetic and acoustic waves, as well as the diffusion of heat between two points in space, with important consequences for thermal management and energy harvesting. There has been significant recent interest in materials with time-modulated properties, which have been shown to efficiently break reciprocity for light, sound, and even charge diffusion. However, time modulation may not be a plausible approach to break thermal reciprocity, in contrast to the usual perception. We establish a theoretical framework to accurately describe the behavior of diffusive processes under time modulation, and prove that thermal reciprocity in dynamic materials is generally preserved by the continuity equation, unless some external bias or special material is considered. We then experimentally demonstrate reciprocal heat transfer in a time-modulated device. Our findings correct previous misconceptions regarding reciprocity breaking for thermal diffusion, revealing the generality of symmetry constraints in heat transfer, and clarifying its differences from other transport processes in what concerns the principles of reciprocity and microscopic reversibility.

R eciprocity is a fundamental property of wave propagation 1,2 and diffusion 3 , implying symmetric field transport in opposite directions. Breaking reciprocity in energy and information transport 4 is essential in components such as diodes, isolators [5][6][7][8] , rotators 9 , rectifiers 10,11 , and circulators [12][13][14] , spanning electromagnetics, photonics, and acoustic domains. Besides the effects of reciprocal heat transfer in static 15 and moving components [16][17][18][19][20] , breaking the symmetry of heat transfer to achieve thermal non-reciprocity 21 is also of great importance for various applications. Devices like heat pipes and thermosyphon diodes are commonly used for thermal management and energy harvesting 22 . In addition, solid-state thermal diodes 23 and rectifiers 24 are the basic elements for thermal information processing in analogy to electronics 25,26 .
In general, there are three types of approaches that break reciprocity. The first is to apply an external bias that is an odd function of time under time-reversal symmetry, like magnetic fields or mechanical motion 12,27,28 . For heat transfer, a simple external bias can be realized by introducing mass or energy fluxes that enter and leave the system with a preferred directionality. Such straightforward approach is not very practical, because it usually makes the underlying systems hardly integrable. The second is by using nonlinearity 8,[29][30][31] . Asymmetric thermal conduction has been found in nonlinear materials 32 with temperature-dependent properties such as oxides 33 or shapememory alloys 34,35 , but the reliance on exotic materials limits its applicability and working conditions.
The third approach to break reciprocity, inspired by recent efforts in electromagnetics and acoustics 36,37 , has been based on materials with time-varying properties. This scheme has received growing interest, since it is easier to be integrated and broadly applicable compared to the first two approaches. The propagation of electromagnetic waves in coupled waveguides has been shown to be non-reciprocal when the electric permittivity ε is modulated with a traveling wave 36 (Fig. 1a), thanks to asymmetric mode conversions. Similar ideas have been successfully applied to thermal radiation 38 and acoustic waves 39,40 . Interestingly, time modulation can also induce asymmetric transfer of electric charge, which is essentially a diffusive process 41 . Intuitively, this is possible because the governing equation, i.e., Fick's law, contains the same Laplacian term as the wave equation. Different from wave propagation, two material parameters in the diffusion equation-the capacitance C e and electric conductivity σ must be modulated simultaneously (Fig. 1b) to achieve this effect.
Conductive heat transfer in solids is another fundamental diffusive process, whose governing equation (Fourier's law) has the same form as Fick's law. The counterpart of electric conductivity σ is thermal conductivity κ, while the counterpart of capacitance C e is the product of density and specific heat capacity ρc. In practice, the specific heat capacity c is hardly tunable, so we only consider the modulation of density ρ and thermal conductivity κ in the following discussion. It appears quite reasonable to expect thermal non-reciprocity induced by such time modulation 42 , considering the continuous success of this approach in electromagnetic 36,37 and acoustic 39,40 wave propagation, and charge diffusion 41 .
In this work, however, we prove theoretically and present numerical and experimental evidences that it is extremely difficult to break reciprocity in heat transfer using time modulation without resorting to external bias or special materials (Fig. 1c). This result is due to the fact that a time modulation of the density inevitably alters the governing transfer equation by taking into account the necessary mass motion v. Our findings indicate that diffusive heat transfer presents inherent constraints that must be carefully treated in its manipulation to break reciprocity.

Results
Diffusion equation under time modulation. Heat transfer in solids is governed by the diffusion Fourier's law: ∂(ρcT)/ ∂t = ∇·(κ∇T), where T(r,t) is the temperature field, r is position vector, and t is time. If the material is linear and not dynamic, the solutions strictly obey reciprocity 43 . In the case of dynamic materials, the density and thermal conductivity vary with time. If both parameters can be freely modulated without introducing additional effects, the Fourier's law becomes Since Eq. (1) has the same form as the time-modulated Fick's law 41 , it is expected that the solution would be in general nonreciprocal 42 . However, as we will discuss in the following, it is impossible to freely modulate the density, since matter that acts as the carrier of thermal energy cannot be created or destroyed. The variation of density ρ must obey the law of mass conservation, which leads to a different governing equation than Eq. (1). Fig. 1 Transport processes in dynamic materials under time modulation. a Non-reciprocal propagation of electromagnetic wave can be induced by spatial-temporally modulating the electric permittivity ε(x,t) as a traveling wave (with speed v 0 ). b Non-reciprocal diffusion of electric charges can be induced by spatial-temporally modulating the capacitance C e (x,t) and electric conductivity σ(x,t). c The reciprocity of heat transfer cannot be broken by modulating the density ρ(x,t) and thermal conductivity κ(x,t) since it is preserved by the continuity equation. Following the law, mass movements v(x,t) (green arrows) must exist to achieve density modulation. d, e Two types of mass movements (dark gray spheres): d along the heat transfer path (gray region). e moving in and out of the heat transfer path.
In continuous media, mass conservation is preserved by the continuity equation ∂ρ/∂t + ∇·(ρv) = 0 44 . If the density varies with time, we inevitably expect mass movement with velocity v (Fig. 1c). Since thermal energy is inherent in any material, the movement introduces a convective term ∇·(ρcTv), which does not appear in Eq. (1). This implies that Eq. (1) can hardly be realized within a physical system without providing external energy or mass. By adding the convective term into Eq. (1), a massconserving diffusive heat transfer under time modulation becomes the convection-diffusion process where we assume that there is no other thermal effect and that viscous dissipation is negligible. Equation (2) is a correction of Eq. (1) replacing a partial derivative by a material derivative (see Supplementary Note 1 for derivation). The detailed effects of the convective term depend on the velocity field v(r,t). In order to focus on the mechanism of time modulation, it is reasonable to only study systems at time-harmonic steady state without externally applied directional mass or energy flux. To be specific, the boundary conditions are constant, and the modulation of material parameters is periodic with time to ensure a stable (timeharmonic) field. No external mass flux exists: ρv · n = 0, where n is the unit normal vector at the system boundary. In timemodulated systems, it is suitable to require instead that the average external mass flux in a time period t 0 vanishes No accumulated external bias is a central assumption that will be used throughout our analysis. There are only two types of setups that support density modulation without external bias: in the first case, the density is modulated by mass motion along the heat transfer path (Fig. 1d). The spheres illustrated do not represent microscopic particles but macroscopic components. There can be an exchange of mass between the system and two thermal reservoirs at both ends, but the three parts together restore the original state after a period, hence there is no net directional mass flow. In the second one, the density is modulated by mass entering or leaving the heat transfer path cyclically (Fig. 1e).
In order to determine whether a two-port system is reciprocal in diffusion, we define the thermal reciprocity based on the equivalence between steady-state global non-reciprocity and thermal diode effect 21 . Specifically, if the system is reciprocal, the heat transfer is symmetric before and after swapping the boundary conditions at two ports, then it is required that the time-averaged heat flux 〈q〉 at harmonic steady state satisfies where the subscripted b and f represent the cases before and after the exchange of boundary conditions, and 1 and 2 represent the position of two ports. In the following, we prove that the heat transfer in both setups is inherently reciprocal and experimentally build a setup of the second type to validate our prediction.
Density modulations of the first type. The first type of modulation scheme follows the one-dimensional (1D) model shown in Fig. 2a. The density ρ(χ) and thermal conductivity κ(χ) of the material are d-periodic functions of χ = xv 0 t, so their profiles move at constant speed v 0 along x. According to the 1D continuity equation ∂ρ/∂t + ∂(ρv)/∂x = 0, the mass flux ρv satisfies ρv where ρ 0 is the average density and C is a constant (See Supplementary Note 2 for derivation). According to Eq. (3), there is no accumulated mass flux through the system in a time period t 0 = d/v 0 (Fig. 2b). Thus, we have C = 0 (see Supplementary Fig. 1 for the effects of a nonzero C) and can solve for the velocity field as v(χ) = [ρ(χ) − ρ 0 ]v 0 /ρ(χ). The 1D heat transfer then obeys We apply fixed temperature boundary conditions T(0,t) = T cold and T(L,t) = T hot (backward) or T(0,t) = T hot and T(L,t) = T cold (forward) at the two ends, respectively. For Eq. (5), such a symmetry can be proved by comparing the forward and backward heat fluxes. Assuming that T b (x,t) is the solution for the backward case, while T f (x,t) is the solution for the forward case. Given any initial conditions, both solutions at time-harmonic steady state should be unique.
It is easy to check that T s (x,t) = T hot + T cold is a solution, and must be the unique solution thanks to the uniqueness of T f (x,t) and T b (x,t). The heat flux q(x,t) is the sum of conductive and convective heat flux: where the reference temperature T ref is an arbitrary constant that can be selected based on convention (We set T ref equal to the lowest temperature T cold ). Then the forward and backward heat fluxes q f (x,t) and q b (x,t) then satisfy Averaging over time gives 〈q f (x)〉 + 〈q b (x)〉 = 0. Considering that the average heat fluxes in and out of the system should balance, we have 〈q f (0)〉 = 〈q f (L)〉 = −〈q b (0)〉 = −〈q b (L)〉, which meets the condition for a symmetric heat transfer as in Eq. (4), and indicates thermal reciprocity. We can also analytically solve Eq. (5). After a variable change (x,t) to (χ = x − v 0 t, τ = t), it is easy to see that Eq. (5) is periodic on χ, so Floquet-Bloch theorem applies and gives (see Supplementary Note 2 for details): where f(χ) is a periodic function with periodicity d, and K is the Bloch wavenumber. The time-harmonic temperature solution should be a linear combination of Eq. (7). The dissipative nature of heat transfer indicates that K must have nonzero imaginary parts 45 . Substituting Eq. (7) into (5), the temperature field can be analytically solved using the Fourier series of f(χ), based on the periodicity of ρ(χ) and κ(χ).
For comparison, we also plot analytical and numerical solutions to the diffusion equation 42 The solutions to Eq. (5) for the heat transfer are symmetric in the backward (Fig. 2c) and forward (Fig. 2d) directions for all modulation parameters characterized by μ, Δ κ , and Δ ρ . This is in contrast with the solutions to Eq. (8) with concave/convex profiles (see Supplementary Note 3 for accurate solutions and Supplementary Fig. 2c, d for heat flux distributions). The insightful work 42 shows the possibility to generate thermal nonreciprocity at a mathematical level, with the material parameters in Eq. (8) set as space-and time-dependent functions. However, Eq. (8) is a hypothetical one that requires an external energy source, because density modulation ρ(χ) cannot be achieved at no cost, and the implementation of modulation will make the governing equation deviate from the previous theoretical design. In a virtual system following Eq. (8), there should be a difference between the average heat flux entering and leaving the system at the two ends, showing that additional energy input or extraction is required to compensate it.
Density modulations of the second type. Another way to achieve density modulation without net directional flow is to add/remove matter periodically through the heat transfer path. The simplest setup of such type is the two-dimensional (2D) model shown in Fig. 3a, where the heat transfer path under consideration is the transverse narrow region at y = 0. We assume constant temperatures on the left and right sides at x = 0 and L, while periodic boundaries are assumed on the upper and lower sides at y = ±d y / 2. ρ(x,y = 0,t) and κ(x,y = 0,t) are d-periodic functions of x − v 0 t, so we consider 2D distributions that are d-periodic functions of ζ = x + ηy − v 0 t with η = d/d y .
Next, we consider mass motion along y with speed v y (x,y,t) to locally modulate the density. According to the continuity equation ∂ρ(ζ)/∂t + ∂[ρ(ζ)v y ]/∂y = 0, we find ρv y = (ρ − ρ 0 )v 0y + C, where v 0y = v 0 /η. The derivation can be obtained in the same way as in 1D case. When C = 0, the case is realizable with mass oscillations in y direction, which is almost the same as the 1D case. Here, we apply periodic condition to the upper and lower boundaries and set C = ρ 0 v 0y , so that v y (ξ,y,t) = v 0y (Fig. 3a). It is noted that Eq. (3) is still satisfied, and this case can be realized on the side surface of a rotating cylinder. The 2D heat transfer follows in which we consider the general case with a thermal conductivity modeled as an anisotropic tensor. Its xx and yy components are κ x and κ y , while the off-diagonal components are assumed to be zero. Similar to the 1D case, we can prove the thermal reciprocity by analyzing the time-averaged heat fluxes along x direction in forward and backward regimes, which also gives 〈q f 〉 + 〈q b 〉 = 0 and satisfy Eq. (4). The solution to Eq. (9) can be solved with a similar method as in 1D case (see Supplementary Note 4). A practical 3D setup that realizes the time modulation of material parameters is shown in Fig. 3b, which consists of fixed and moving fan-shaped solid plates with density ρ A = 8390 kg m −3 , heat capacity c A = 375 J kg −1 K −1 , and thermal conductivity κ A = 123 W m −1 K −1 . Each plate spans π/2 with inner and exterior radius R 1 = 1 cm and R 2 = 2 cm, and thickness δ = 0.25 cm = d/16. The total length of the system L = 5d = 20 cm. Temperature boundary conditions are T(0,t) = T cold and T(L,t) = T hot (backward) or T(0,t) = T hot and T(L,t) = T cold (forward), with constant temperatures set as T cold = 273 K and T hot = 323 K. All other boundaries are thermally insulated. Naturally, we can regard the heat transfer path as the portion of the system contained in the region ([R 1 ,R 2 ], [π/4,3π/4], [0,L]) of the cylindrical coordinate system (r, θ, x), through which most of the heat flux is conducted. Along the x direction, each moving plate is π/4 ahead of the previous moving one. All of them rotate at angular speed −Ω = −0.06π rad s −1 . As the moving plates enter or leave the heat transfer path at y = 0, the density ρ(x,y = 0,t) and thermal conductivity κ(x,y = 0,t) are effectively modulated (Supplementary Note 5).
We perform numerical simulations on the 2D and 3D models in Fig. 3a, b. As expected, reciprocal heat transfer is confirmed for Fig. 2 Heat transfer under 1D density modulation. a Density ρ and thermal conductivity κ (color maps) move as traveling waves at speed v 0 with wavelength d. To achieve density modulation, actual mass movements (green arrows) at speed v must exist. b The total mass flux in a time period should be zero to keep a cyclic and close setup. c, d Backward (c) and forward (d) temperature distributions of the system (Eq. (5)) at t = Nd/v 0 (N is a large enough integer to achieve time-harmonic steady state), compared with those of a virtual system without mass movements (Eq. (8) all considered temperature distributions along the line at (r = R 2 , θ = π/2) (Fig. 3c, d). The temperature distributions on the plate surfaces are also plotted in the insets (see Supplementary Movie 1 for the evolution with time and Supplementary Fig. 3a for the heat flow in and out of the system). Reciprocity can be easily shown from the symmetric temperature profiles and the identical time-averaged heat fluxes in forward and backward directions. The reason for the absence of non-reciprocity resides in the fact that in these diffusive systems the density variations are averaged out over time due to the conservation of mass and hence there is effectively no density modulation. To further illustrate this point, we performed simulations on a hypothetical model where the density of the moving plates is artificially modulated as ρ A (t) = ρ A {1 + cos[2Ωt − (n − 1)π/4]}/2, for the n-th layer of moving plates. The temperature distributions in Fig. 3e, f show concave and convex profiles, indicating non-reciprocity (see Supplementary Movie 2 for the evolution and Supplementary Fig. 3b for the heat flow). Experimental verification of reciprocal heat transfer in a timevarying material. We have implemented the 3D geometry in Fig. 3b, as shown in Fig. 4a. The system is built with fan-shape plates attached to a fixed beam and a shaft revolved by a lowspeed motor (angular velocity Ω = 0.075π rad s −1 ) to generate a right-handed rotating profile. The plates are made of brass with the same material properties and shape as in the 3D simulations. The supporting beam and shaft are made of nylon with thermal conductivity of 0.3 W m −1 K −1 . Two copper blocks are mounted at the ends of the supporting beam to contact the heat sources directly at 296 and 306 K, generating a temperature difference ΔT = T hot -T cold = 10 K. The blocks are slotted at the top to match the upper four rotating fan-shape plates. Different from the numerical simulations, a 0.2 mm gap is made between adjacent plates, which is filled with thermal grease (thermal conductivity 4.38 W m −1 K −1 ) for conduction and silicone oil (thermal conductivity 0.16 W m −1 K −1 ) for lubrication. Therefore, the interface thermal resistance is nonnegligible.
In addition, the natural convective heat exchange with the ambient air introduces another term h(T ∞ − T), where h is the heat exchange rate and T ∞ is the room temperature. In order to reduce the influence of ambient air convection, we set the system in a vacuum chamber (see Supplementary Fig. 4). Molecular pump (LF-110) and mechanical pump (BSV-16) are employed to reduce the density of air to as low as 10 −3 Pa. Temperature profiles are measured through the inspection window with an infrared camera (Fotric 347). Due to the interface thermal resistance as well as natural convective heat exchange, the temperature gradients at the center are smaller compared to the linear profile (Fig. 4b). However, the two factors do not have any directionality, so the temperature profiles for the backward and forward cases are still symmetric, demonstrating reciprocal heat transfer in the time-modulated system. See Supplementary Movie 3 for videos of the rotating device and the temperature profiles. Reciprocity is further confirmed by the temperature distribution along the heat transfer path of the system (dashed green line in Fig. 4a), where the fixed plates are placed (Fig. 4c). In Fig. 4c, to compare with the experimental results, numerical simulations were performed on the 3D model with interface thermal resistance (0.2 mm thickness and 1 W m −1 K −1 thermal conductivity) and natural air convection (h = 3 W m −2 K −1 , and T ∞ was set as the ambient temperature 300. experimentally measured curves for backward and forward heat transfer are not only symmetric, but also in good agreement with the simulated results. Since the heat flux in x direction is proportional to temperature gradient q x = −κ∂T/∂x, and the equivalent thermal conductivity is homogeneous throughout the system, the reciprocity of heat flow can be reflected from the symmetry of the temperature difference over two intervals of equal length (marked in Fig. 4a). We give the measured ΔT 1 and ΔT 2 in Fig. 4d, and averaging the temperature difference over time gives: 〈ΔT b,1 〉 = −0.874 K, 〈ΔT b,2 〉 = −1.038 K (backward) and 〈ΔT f,1 〉 = 1.046 K, 〈ΔT f,2 〉 = 0.863 K (forward), showing a symmetric relationship: This symmetry is maintained for all angular velocities of the shaft rotated by the motor. As another example in Fig. 4e, the angular velocity Ω is changed to 0.14π rad s −1 and the timeaveraged temperature difference becomes: 〈ΔT b,1 〉 = −0.887 K, 〈ΔT b,2 〉 = −0.992 K (backward) and 〈ΔT f,1 〉 = 0.996 K, 〈ΔT f,2 〉 = 0.872 K (forward). In order to further demonstrate the reciprocity, we give the numerical results of the heat flow Q in and out of the system (see Supplementary Fig. 5). It is easy to check that 〈Q b,1 〉 = −〈Q f,2 〉 and 〈Q b,2 〉 = −〈Q f,1 〉, satisfying the condition for thermal reciprocity in Eq. (4). This reciprocal result exhibited in the experiment is a representative of the general situation where the modulation parameters (e.g., rotating angular velocity) are changed arbitrarily.

Discussion
In this paper, we have shown that physical systems preserve thermal reciprocity under time modulation as a result of a convection correction introduced by density modulation. For other processes governed by the momentum equation (e.g., wave propagation) or the continuity equation (e.g., charge diffusion), the driving approach of time modulation does not alter the form of the governing equation. For heat transport, however, it is a wellknown fact that a material derivative must be used if there is a transport of mass. As such, a convection correction should be added in the governing equation, resulting in the reciprocity in heat transfer. Though only several simulation and experiment examples are given in this work, the conclusion of reciprocity is theoretically rigorous. One evidence is that the proven symmetric relationship between the forward and backward heat flux 〈q f 〉 = −〈q b 〉 always exists, regardless of how the material parameters are modulated. Even if the modulated material parameters are asymmetric in spatial distribution, the system is reciprocal in heat transfer ( Supplementary Fig. 6). While we did not explicitly discuss the case of a varying specific heat capacity c, it is easy to recognize that our arguments equally apply to this scenario. First of all, if the variation is simply a consequence of mechanical motions, the material derivative of ρc must be zero, which gives ∂(ρc)/∂t + ∇·(ρcv) = 0. Therefore, Eq.
(2) and the following analysis still apply. Second, if the specific heat capacity of the material can really be modulated at will, it is possible to avoid the convective term and thereby in principle achieve non-reciprocity. However, this is only possible in very limited scenarios, such as using caloric materials that undergo phase transitions in the presence of electric 46 or magnetic 47 fields. In addition, there are schemes that can break thermal reciprocity beyond our assumptions. Though thermal reciprocity is preserved as the boundary conditions are maintained at two constant temperatures, the non-reciprocity of thermal wave 48,49 under space-time modulation may occur when the boundaries are set as periodically oscillating heat sources 50 . Another trivial way to generate non-reciprocity is introducing directional mass/energy bias while modulating the parameters.
Our work suggests that thermal reciprocity has more fundamental resilience than other transport mechanisms. These findings may have important implications for the design of thermal devices and other dissipative wave propagation systems.