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.

Introduction

Reciprocity 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-8], rotators[9], rectifiers[10,11], and circulators[12-14], spanning electromagnetics, photonics, and acoustic domains. Besides the effects of reciprocal heat transfer in static[15] and moving components[16-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-31]. Asymmetric thermal conduction has been found in nonlinear materials[32] with temperature-dependent properties such as oxides[33] or shape-memory 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 and electric conductivity σ must be modulated simultaneously (Fig. 1b) to achieve this effect.

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 v0). b Non-reciprocal diffusion of electric charges can be induced by spatial-temporally modulating the capacitance C(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.

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 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 becomesSince Eq. (1) has the same form as the time-modulated Fick’s law[41], it is expected that the solution would be in general non-reciprocal[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).

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 mass-conserving diffusive heat transfer under time modulation becomes the convection-diffusion processwhere 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 (time-harmonic) field. No external mass flux exists: ρv · n = 0, where n is the unit normal vector at the system boundary. In time-modulated systems, it is suitable to require instead that the average external mass flux in a time period t0 vanishesNo 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 satisfieswhere 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 χ = x – v0t, so their profiles move at constant speed v0 along x. According to the 1D continuity equation ∂ρ/∂t + ∂(ρv)/∂x = 0, the mass flux ρv satisfies ρv = (ρ − ρ0)v0 + C, 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 t0 = d/v0 (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]v0/ρ(χ). The 1D heat transfer then obeysWe apply fixed temperature boundary conditions T(0,t) = Tcold and T(L,t) = Thot (backward) or T(0,t) = Thot and T(L,t) = Tcold (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 Tb(x,t) is the solution for the backward case, while Tf(x,t) is the solution for the forward case. Given any initial conditions, both solutions at time-harmonic steady state should be unique. Their summation Ts(x,t) = Tf(x,t) + Tb(x,t) also satisfies Eq. (5) with boundary conditions Ts(0,t) = Ts(L,t) = Thot + Tcold. It is easy to check that Ts(x,t) = Thot + Tcold is a solution, and must be the unique solution thanks to the uniqueness of Tf(x,t) and Tb(x,t). The heat flux q(x,t) is the sum of conductive and convective heat flux: q(x,t) = −κ∂T/∂x + ρcv(T − Tref), where the reference temperature Tref is an arbitrary constant that can be selected based on convention (We set Tref equal to the lowest temperature Tcold). Then the forward and backward heat fluxes qf(x,t) and qb(x,t) then satisfyAveraging over time gives 〈qf(x)〉 + 〈qb(x)〉 = 0. Considering that the average heat fluxes in and out of the system should balance, we have 〈qf(0)〉 = 〈qf(L)〉 = −〈qb(0)〉 = −〈qb(L)〉, which meets the condition for a symmetric heat transfer as in Eq. (4), and indicates thermal reciprocity.

Fig. 2

Heat transfer under 1D density modulation.

a Density ρ and thermal conductivity κ (color maps) move as traveling waves at speed v0 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/v0 (N is a large enough integer to achieve time-harmonic steady state), compared with those of a virtual system without mass movements (Eq. (8)). Scatter points are simulated results, lines are analytical results, and dashed lines are analytical solutions of the homogenized profiles. For clarity, the results of Eq. (5) at modulating speed μ = 4/d are shifted to have ±0.1 offsets.

We can also analytically solve Eq. (5). After a variable change (x,t) to (χ = x − v0t, τ = 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 κ(χ).

To verify the solution, we build a 1D model with d = 1 cm and L = 10d. The density and thermal conductivity are set to be ρ(χ) = ρ0[1 + Δcos(βχ)] and κ(χ) = κ0[1 + Δcos(βχ)], where β = 2π/d, ρ0 = 2000 kg m−3, Δ = 0.3, κ0 = 100 W m−1 K−1, and Δ = 0.9. The specific heat capacity is c = 1000 J kg−1 K−1. We choose two modulation speeds v0 = μκ0/ρ0c with μ = 1/d and 4/d. Constant temperatures are set as Tcold = 273 K and Thot = 323 K to generate a temperature difference ΔT = Thot – Tcold = 50 K. Our analytical results (solid lines in Fig. 2c, d) are well validated by finite-element simulations (scatter points in Fig. 2c, d) with COMSOL Multiphysics®. The corresponding heat flux distributions are summarized in Supplementary Fig. 2a, b, demonstrating the thermal reciprocity.

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 non-reciprocity 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/2. ρ(x,y = 0,t) and κ(x,y = 0,t) are d-periodic functions of x − v0t, so we consider 2D distributions that are d-periodic functions of ζ = x + ηy − v0t with η = d/d.

Fig. 3

Heat transfer under 2D and 3D density modulation.

a Density ρ and thermal conductivity κ profiles (color maps) of 2D traveling waves such that the profiles move at speed v0 in x direction (marked by dashed lines). The density is modulated by mass movement in y direction at speed v0 (green arrows). The upper and lower boundaries are periodic. b A three-dimensional (3D) model in a (r,θ,x) cylindrical coordinate system. Similar ρ and κ modulations are achieved with fixed (white) and moving (yellow) fan-shaped plates. c–f Backward (c, e) and forward (d, f) temperature distributions on the top line r = R2, θ = π/2 extracted from simulated results for the 2D (lines) and 3D (scatters) models at t = N2π/Ω (N is a large enough integer to achieve time-harmonic steady state). The insets show the entire temperature distributions. The results in e and f use a hypothetical 3D model where the moving plates have time-varying masses, violating the law of mass conservation.

Next, we consider mass motion along y with speed v(x,y,t) to locally modulate the density. According to the continuity equation ∂ρ(ζ)/∂t + ∂[ρ(ζ)v]/∂y = 0, we find ρv = (ρ − ρ0)v0 + C, where v0 = v0/η. 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 = ρ0v0, so that v(ξ,y,t) = v0 (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 followsin which we consider the general case with a thermal conductivity modeled as an anisotropic tensor. Its xx and yy components are κ and κ, 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 〈qf〉 + 〈qb〉 = 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 cA = 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 R1 = 1 cm and R2 = 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) = Tcold and T(L,t) = Thot (backward) or T(0,t) = Thot and T(L,t) = Tcold (forward), with constant temperatures set as Tcold = 273 K and Thot = 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 ([R1,R2], [π/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 all considered temperature distributions along the line at (r = R2, θ = π/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 time-varying 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 low-speed 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 = Thot – Tcold = 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.

Fig. 4

Experimental verification of reciprocal heat transfer under time modulation.

a Photo of the system built according to the 3D modulation model in Fig. 3. The movable brass plates are rotated to generate a right-ward moving profile. b Temperature profiles for backward and forward heat transfer. c Experimentally measured temperature profiles along the heat transfer path of the system (dashed green line in a). The scatter points refer to simulated results considering natural air convection and interface thermal resistance. d, e Experimentally measured temperature difference (between two ends of the red segment line in a) changing with time t. The angular velocity of the shaft rotated by the motor is 0.075π rad s−1 (d) and 0.14π rad s−1 (e). The symmetry of time-averaged temperature difference implies the reciprocity of heat flow in Eq. (4).

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.5 K.) The 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 = −κ∂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 ΔT1 and ΔT2 in Fig. 4d, and averaging the temperature difference over time gives: 〈ΔTb,1〉 = −0.874 K, 〈ΔTb,2〉 = −1.038 K (backward) and 〈ΔTf,1〉 = 1.046 K, 〈ΔTf,2〉 = 0.863 K (forward), showing a symmetric relationship: 〈ΔTb,1〉 ≈ −〈ΔTf,2〉 and 〈ΔTb,2〉 ≈ −〈ΔTf,1〉. 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 time-averaged temperature difference becomes: 〈ΔTb,1〉 = −0.887 K, 〈ΔTb,2〉 = −0.992 K (backward) and 〈ΔTf,1〉 = 0.996 K, 〈ΔTf,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 〈Qb,1〉 = −〈Qf,2〉 and 〈Qb,2〉 = −〈Qf,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 well-known 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 〈qf〉 = −〈qb〉 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.