⇦Chinese Physics Letters, 2021, Vol. 38, No. 3, Article code 030501 Advection and Thermal Diode Ying Li (李鹰)1,2,3* and Jiaxin Li (李佳鑫)4 Affiliations 1Interdisciplinary Center for Quantum Information, State Key Laboratory of Modern Optical Instrumentation, College of Information Science and Electronic Engineering, Zhejiang University, Hangzhou 310027, China 2ZJU-Hangzhou Global Science and Technology Innovation Center, Key Lab of Advanced Micro/Nano Electronic Devices & Smart Systems of Zhejiang, Zhejiang University, Hangzhou 310027, China 3International Joint Innovation Center, Zhejiang University, Haining 314400, China 4School of Mechatronics Engineering, Harbin Institute of Technology, Harbin 150001, China Received 28 November 2020; accepted 8 January 2021; published online 2 March 2021 *Corresponding author. Email: eleying@zju.edu.cn Citation Text: Li Y and Li J X 2021 Chin. Phys. Lett. 38 030501    Abstract We prove that under the condition of closed boundary to mass flux, pure advection is not a valid mechanism to make a practical thermal diode. Among the various designs of thermal diodes, many of them involve circulating fluid flow, such as in thermosyphons. However, those designs often employ natural convection, which is basically a nonlinear process. It thus remains unclear how the pure advection of temperature field induced by a decoupled velocity field influences the symmetry of heat transfer. Here we study three typical models with pure advection: one with open boundary, one with closed boundary at unsteady state, and one with closed boundary at steady state. It is shown that only the last model is practical, while it cannot become a thermal diode. Finally, a general proof is given for our claim by analyzing the diffusive reciprocity. DOI:10.1088/0256-307X/38/3/030501 © 2021 Chinese Physics Society Article Text A common symmetry called reciprocity is preserved in many transport processes such as the propagation of electromagnetic,[1] acoustic,[2] and elastic waves,[3] or the diffusion of heat, mass, and charge carriers.[4] It means that the influence of a source on a target should remain the same after their spatial positions are swapped. For systems with two ports, reciprocity implies that the relation between output and input fields is independent of which port receives or sends out signal, energy, or matter. Despite of some issues with lossy materials,[5] reciprocity is believed to be closely related with time-reversal symmetry[6] in which processes in opposite directions of time evolution should be indistinguishable. For macroscopic diffusions, time-reversal symmetry is obviously broken because the amplitude of fields always decays with time. However, reciprocity is still preserved in them, partly because the dynamics of microscopic particles is time-reversal symmetric.[7,8] Time-reversal symmetry can be broken by external fields. For example, a magnetic field can introduce Faraday rotation and thereby achieve non-reciprocal transport of polarized light.[9] It can also break Kirchhoff's law for non-reciprocal thermal radiation.[10,11] The basic mechanism is that the magnetic response is an odd function under time reversal, so the direction of the magnetic field should be reversed to preserve time-reversal symmetry after the positions of source and target are swapped. In actual setups, the magnetic field is treated as an external factor so its direction will not be changed with the position of the source, which breaks the time-reversal symmetry and reciprocity.[12] Non-reciprocity is also an essential property for many thermal metamaterials and devices,[13,14] such as thermal diodes[15] that transport different amounts of heat in opposite directions (for a detailed study on the relation between diffusive non-reciprocity and thermal diode effect, see Ref. [4]). To realize a thermal diode, the thermal convection in moving components is often employed.[16,17] Since the convective term is also an odd function under time reversal, a velocity field should have similar effects as a magnetic field to break the reciprocity of heat transfer. However, it remains a question whether this is the working mechanism of convective thermal diodes, because many other effects such as gravity and thermal expansion are also involved in those devices. To be specific, consider the governing equation for temperature $T({\boldsymbol r},t)$ as a function of time $t$ and position vector ${\boldsymbol r}$: $$ \rho c_{\rm p} \frac{\partial T}{\partial t}=\nabla \cdot ({\kappa \nabla T})-\rho c_{\rm p} \boldsymbol{v}\cdot \nabla T+g,~~ \tag {1} $$ where $\rho$ is density, $c_{\rm p}$ is specific heat capacity at constant pressure, $\kappa$ is thermal conductivity, ${\boldsymbol v}$ is velocity field, and $g$ is density of heat generation.[18] Obviously, if the velocity field ${\boldsymbol v}$ is independent of the temperature field $T$, the term $\rho c_{\rm p}{\boldsymbol v}\cdot\nabla T$ will be an odd function under time reversal that can break the reciprocity. This effect of a decoupled velocity field, such as in forced convection, is simply advection,[19] where the heat transfer process can be regarded as a pure thermal conduction in a moving background, as can be easily seen by making a variable change ${\boldsymbol r}' = {\boldsymbol r} - {\boldsymbol v} t$. However, many convective thermal diodes employ natural convection, where the velocity field is generated by temperature difference in the fluid. In a gravity field, the temperature difference induces thermal expansion and changes the density distribution, which leads to buoyancy driven motions. Moreover, the motion of the fluid will generate viscous heat as a part of $g$, so $g$ is a function of ${\boldsymbol v}$. Therefore, the density $\rho$, velocity ${\boldsymbol v}$, and heat generation $g$ are all temperature dependent, making the heat transfer a complex nonlinear process. It is well known that even without any convection, nonlinear materials such as those with temperature-dependent thermal conductivities[20] can be used to build a thermal diode.[21–24] It is thus unclear how much the advection effect contributes to realization of a thermal diode. On the contrary to the intuition, here we show that a practical thermal diode cannot be realized using pure advection. Our argument is based on an observation that a thermal diode must have a closed boundary to mass flow in practical situations. This means that irrespective of its detailed form, the velocity field ${\boldsymbol v}$ must satisfy a boundary condition ${\boldsymbol v}\cdot{\boldsymbol n} = 0$, where ${\boldsymbol n}$ is the unit normal vector to the boundary of the device. The necessity of this condition is obvious: if there are mass fluxes entering and leaving the device, the device will not have a well-defined boundary to be integrated or applied to a heat transfer system. Nevertheless, let us first examine a system that violates the condition as in Fig. 1(a). Consider a long channel with uniform flow of mass at constant speed $v$ along $x$ direction. Its upper and lower boundaries are thermally insulated. We forcefully select a part in the channel as the system (marked by dashed lines). Pure advection is assumed by ignoring the temperature dependence of all parameters in Eq. (1), as well as the viscous heat generation. At steady state, the governing equation is thus $$ D\frac{\partial^{2}T}{\partial x^{2}}-v\frac{\partial T}{\partial x}=0,~~ \tag {2} $$ where $D=\kappa/\rho c_{\rm p}$ is the diffusivity which is assumed to be uniform and isotropic for simplicity. A general solution to it is $$ T(x)=Ae^{{vx}/D}+B,~~ \tag {3} $$ where $A$ and $B$ are constants. To study the heat transfer through it in forward ($x$) and backward ($-x$) directions, the two ends ($x = 0$, $L$) of the system are maintained at constant temperatures $T_{1}$ and $T_{2}$. In practice, such a boundary condition requires two thermal reservoirs placed inside the channel, which will generally block the mass flow. It supports that the condition of a boundary closed to mass fluxes is reasonable, but here we proceed by assuming that the thermal reservoirs have porous structures[25] that do not alter the mass flow significantly. Combined with the boundary condition, Eq. (3) becomes $$ T(x)=T_{1} +(T_{2} -T_{1})\frac{e^{{vx} / D}-1}{e^{{vL} / D}-1}.~~ \tag {4} $$ The total heat flux through the system is constant along $x$, which is the sum of the conductive and convective parts, $$\begin{alignat}{1} q=-\kappa \frac{\partial T}{\partial x}+\rho c_{\rm p} vT =\rho c_{\rm p}v\Big({T_{1} +\frac{T_{1} -T_{2} }{e^{{vL}/D}-1}}\Big).~~ \tag {5} \end{alignat} $$ We set the parameters as $\kappa = 1$ W/(m$\cdot$K), $\rho c_{\rm p} = 10^{6}$ J/(m$^{3}\cdot$K), and $L = 0.1$ m, and plot the dependence of $q$ on the temperature difference $\Delta T=T_{1} - T_{2}$ in Fig. 1(b), where the average temperature $T_{0} = (T_{1} + T_{2})/2$ is fixed at 300 K. Note that Eq. (5) reduces to the conductive heat flux $\kappa (T_{1} - T_{2})/L$ when $v = 0$. If $T_{1} > T_{2}$, Eq. (5) gives the total heat flux $q_{\rm f}$ in the forward direction. Swapping $T_{1}$ and $T_{2}$ gives the total heat flux $q_{\rm b}$ in the backward direction. From Fig. 1(b), we see that the total heat flux linearly depends on the temperature difference, but is asymmetric around $\Delta T = 0$ for nonzero $v$. Therefore, $q_{\rm f} + q_{\rm b} \ne 0$, which meets the common definition of a thermal diode. Assuming $T_{1} > T_{2}$, the rectification ratio $\gamma$ is defined as $$\begin{align} \gamma&={(q_{\rm f} +q_{\rm b})} / {\max (\vert q_{\rm f} \vert,\vert q_{\rm b} \vert)}\\ &=\frac{{\rm sgn}(v){2T_{0}}}{{{T_{0}+{\rm sgn}(v)[{\frac{1}{2}+\frac{1}{(e^{{vL}/ D}-1)}}]\Delta T}}},~~ \tag {6} \end{align} $$ where sgn($v$) returns $-$1, 0, and 1 for $v < 0$, $v= 0$, and $ v> 0$, respectively. Equation (6) is indefinite for $v = 0$, but it is easy to verify that $\gamma = 0$ at zero velocity. The dependence of $\gamma$ on $vL/D$ is plotted in Fig. 1(c). For common thermal diodes, $\gamma$ is less than 1 as $q_{\rm f} > 0$ and $q_{\rm b} < 0$. However, here their signs depend on the velocity. For example, $q_{\rm b}$ reaches zero at $$ \frac{vL}{D}=\ln\Big({\frac{T_{1}}{T_{2}}}\Big)>0.~~ \tag {7} $$ At this point, the conductive heat in the backward direction is balanced by the convective heat in the forward direction, and the rectification ratio $\gamma$ reaches 1. As $v$ further increases, $q_{\rm b}$ becomes positive despite of the positive temperature gradient along $x$ due to the strong advection. The rectification ratio $\gamma$ becomes greater than 1 and approaches 2 at large $v$, where the conductive heat contributes little to the total heat. It seems that the pure advection system is indeed a thermal diode at nonzero velocity, and even appears to be a very efficient one as the rectification ratio can be greater than 1 (or less than $-$1). However, as mentioned above, the system is actually an inseparable part of the channel. Therefore, it is highly impractical to use it as a real device for heat dissipation, temperature management,[26] heat signal processing,[27,28] etc.
cpl-38-3-030501-fig1.png
Fig. 1. (a) Schematic for an open system with mass fluxes through it (green arrows). Its two ends are maintained at constant temperatures $T_{1}$ (red line) and $T_{2}$ (blue line). (b) Dependence of the total heat flux $q$ on temperature difference $\Delta T$. (c) Dependence of the rectification ratio $\gamma$ on $vL/D$.
Now we turn back to systems that satisfy the condition that no mass flux enters or leaves it. As a closed system, it can be made in direct or indirect contact with the thermal reservoirs for common heat transfer applications. It is immediately noticed that one cannot make such a system by simply enclosing the system in Fig. 1(a), because the mass flux does not conserve at the two ends. One could place an active source (such as a fan) to generate a mass flow near the boundary, such that the velocity changes from 0 to $v$ then back to 0 along $x$, as shown in Fig. 2(a). However, the motion will inevitably modify the density distribution as mass is constantly moved from one end to the other, making it impossible to reach a steady state. More specifically, the continuity equation requires[18] $$ \frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \boldsymbol{v})=0.~~ \tag {8} $$ Therefore, at steady state, the velocity field must satisfy $$ \nabla \cdot (\rho \boldsymbol{v})=0.~~ \tag {9} $$ For the case in Fig. 2(a), we have $$ \frac{d\rho }{dx}v(x)+\rho (x)\frac{dv}{dx}=0.~~ \tag {10} $$ At $x = 0$, the velocity $v(x)$ is zero, but its derivative is not [if the derivative is also zero, one can take a derivative to Eq. (10) and use similar arguments]. Equation (10) thus requires $\rho (0) = 0$ (or $\infty$, as shown later), which is nonphysical.
cpl-38-3-030501-fig2.png
Fig. 2. (a) Schematic for a 1D closed system with interior mass flow (green arrows). Its two ends are maintained at constant temperatures $T_{1}$ (red line) and $T_{2}$ (blue line). The velocity varies along $x$ to ensure no flux at two ends. (b) Distributions of the density $\rho$ and heat flux $q$. A cutoff is made to $\rho$ near the two ends to avoid singularity. (c) Dependence of two rectification ratios on $v_{0}L/D_{0}$.
It is still interesting to study such a system by loosening the steady-state condition a little bit, and requiring that the density distribution changes much slower than the temperature field. As a concrete example, consider the velocity field $$ v(x)=\frac{v_{0}}{2}({1-\cos\alpha}),~~ \tag {11} $$ where $\alpha = 2\pi x/L$, $v_{0}$ is a constant. The velocity vanishes at $x = 0$ and $L$, and reaches maximum $v_{0}$ at $x=L$/2. Substituting Eq. (11) into Eq. (10) gives $$ \rho(x)=\frac{2\rho_{0}}{1-\cos \alpha},~~ \tag {12} $$ where $\rho_{0}$ is a constant. The density approaches infinity at $x = 0$ and $L$, and reaches minimum $\rho_{0}$ at $x=L/2$. To avoid singularity, we can cut the density near $x = 0$ and $L$ within a distance of $d$ to a large finite value. The density distribution $\rho (x)$ is plotted in Fig. 2(b) (green line), where the value of $d$ is chosen to relatively large (5 mm) for a better view. Substituting Eqs. (11) and (12) into the heat transfer equation gives $$ \kappa \frac{\partial^{2}T}{\partial x^{2}}-\rho_{0} c_{\rm p} v_{0} F(x)\frac{\partial T}{\partial x}=0,~~ \tag {13} $$ where $F(x)$ is a continuous function introduced due to the density cutoff. $F(x) = 1$ except for $x$ near 0 and $L$, and $F(0) = F(L) = 0$. Therefore, an integral of $F(x)$ can be approximated with the linear function $x$. Combined with the boundary conditions, the solution has almost the same form as Eq. (4), $$ T(x)=T_{1} +(T_{2} -T_{1})\frac{e^{{v_{0} x} / {D_{0} }}-1}{e^{{v_{0} L} / {D_{0} }}-1},~~ \tag {14} $$ where $D_{0}=\kappa$/($\rho_{0}c_{\rm p}$). However, the heat flux is no longer uniform along $x$, $$\begin{alignat}{1} q(x)={}&\rho_{0} c_{\rm p} v_{0} (T_{1} -T_{2})[1-F(x)]\frac{e^{v_{0} x/D_{0} }}{e^{v_{0} L/D_{0} }-1}\\ &+\rho_{0} c_{\rm p} v_{0} F(x)\Big({T_{1}+\frac{T_{1}-T_{2}}{e^{v_{0}L/D_{0}}-1}}\Big).~~ \tag {15} \end{alignat} $$ The distribution is numerically simulated with COMSOL Multiphysics with $\kappa = 1$ W/(m$\cdot$K), $\rho_{0}c_{\rm p} = 10^{6}$ J/(m$^{3}\cdot$K), $L = 0.1$ m, $T_{0} = 300$ K, $\Delta T = 10$ K, and $d = 0.005$ m. The results are plotted in Fig. 2(b). For most values of $x$, the heat flux is almost the same as in Eq. (5), $$ q(x)=\rho_{0} c_{\rm p} v_{0}\Big({T_{1}+\frac{T_{1}-T_{2}}{e^{v_{0}L/D_{0}}-1}}\Big).~~ \tag {16} $$ However, at the two ends, we have $$\begin{align} &q(0)=\rho_{0} c_{\rm p} v_{0} (T_{1} -T_{2})\frac{1}{e^{v_{0} L/D_{0} }-1},\\ &q(L)=\rho_{0} c_{\rm p} v_{0} (T_{1} -T_{2})\frac{e^{v_{0} L/D_{0} }}{e^{v_{0} L/D_{0} }-1}.~~ \tag {17} \end{align} $$ Obviously, the heat fluxes into and out of the system are unequal. The additional heat comes from the fact that the system is not really at steady state, due to the changes in density distribution. If $T_{1} > T_{2}$, Eq. (17) gives the forward heat fluxes at both ends. The backward heat fluxes are $q_{\rm b}(0) = -q_{\rm f}(0)$ and $q_{\rm b}(L)=-q_{\rm f}(L)$. For a device, we are only interested in the heat flux that enters or leaves it, as the interior heat flux cannot be utilized. Therefore, we determine whether the system is a thermal diode based on $q(0)$ and $q(L)$. Now that $q(0) \ne q(L)$, the rectification ratio cannot be defined using Eq. (6). There are two ways to redefine a thermal diode. One is to calculate their average $q = [q(0)+ q(L)]/2$. Then $q_{\rm f} + q_{\rm b} = 0$, so the system is not a thermal diode in this sense. The other is to compare $q_{\rm f}(0)$ with $q_{\rm b}(L)$ and $q_{\rm f}(L)$ with $q_{\rm b}(0)$, so two rectification ratios can be defined as follows: $$\begin{alignat}{1} \gamma_{1} &=\frac{q_{\rm f} (0)+q_{\rm b} (L)}{\max[{|{q_{\rm f} (0)}|,|{q_{\rm b}(L)}|}]}\\ &={\rm sgn}(v_{0})[{e^{-{\rm sgn}(v_{0})v_{0} L/D_{0}}-1}] \\ \gamma_{2} &=\frac{q_{\rm f} (L)+q_{\rm b} (0)}{\max[{|{q_{\rm f}(L)}|,|{q_{\rm b}(0)}|}]}\\ &={\rm sgn}(v_{0})[{1-e^{-{\rm sgn}(v_{0})v_{0} L/D_{0}}}].~~ \tag {18} \end{alignat} $$ The dependences of them on $v_{0}L/D_{0}$ are plotted in Fig. 2(c), which are in good agreement with the simulated results using $d = 1$ mm. The magnitudes are much smaller than those of an open system in Eq. (6) (note the range of the abscissa), and can never exceed 1, because the direction of the heat flux is always the same as the temperature difference. It is also noted that the rectification ratios are independent of the boundary temperatures. Both $\gamma_{1}$ and $\gamma_{2}$ are nonzero, so the system could be considered as a thermal diode in this sense. However, $\gamma_{1}$ and $\gamma_{2}$ are exactly opposite, so it is still unreasonable to claim that the system favors forward or backward heat transfer. This means that the application of the device is quite limited, not to mention that the system is actually unstable. The above argument shows that there is no stable and closed one-dimensional (1D) configuration with interior mass motion. In higher dimensions, it seems that a practical convective thermal diode system must involve some form of circulating flow like the one in Fig. 3(a). Such a setup is similar to most existing convective thermal diodes, but here we assume no gravity field or thermal expansion. The velocity field is constant and independent of the temperature field. Under pure advection, it is immediately noticed that the configuration in Fig. 3(a) should not be a thermal diode, because of its geometric symmetry. To explain the mechanism, we further simplify the model by assuming a narrow channel such that the velocity is always parallel to the channel with uniform magnitude $v_{0}$, and that the temperature gradient on any cross section of the channel is negligible. We can then study the quasi-1D model as in Fig. 3(b). The parts that are in direct contact with the thermal reservoirs (marked by dashed lines) are assumed to have the same temperature as the thermal reservoirs. The two isothermal regions are connected with two 1D paths. The lower and upper paths are referred to as p1 and p2, respectively. The heat flux has the same form as Eq. (5), $$\begin{align} &q_{{\rm p}1} =\rho c_{\rm p} v_{0}\Big({T_{1} +\frac{T_{1} -T_{2} }{e^{{v_{0} L} / D}-1}}\Big),\\ &q_{{\rm p2}} =-\rho c_{\rm p} v_{0}\Big({T_{1} +\frac{T_{1} -T_{2} }{e^{{-v_{0} L} / D}-1}}\Big).~~ \tag {19} \end{align} $$ Clearly, the difference between the heat fluxes that enters and leaves each isothermal region is balanced by the heat that leaves the system to the adjacent thermal reservoir. For the right end, it is $$ Q=(q_{{\rm p}1} +q_{{\rm p2}})\sigma =\rho c_{\rm p} v_{0} \sigma (T_{1} -T_{2})\frac{e^{v_{0} L/D}+1}{e^{v_{0} L/D}-1},~~ \tag {20} $$ where $\sigma$ is the cross-section area of the channel. When $T_{1} > T_{2}$, Eq. (20) gives the forward heat flux $q_{\rm f}=Q/S$, where $S$ is the contact area between the system and the thermal reservoir. It is noticed that the first terms in the heat fluxes along both paths cancel each other, such that $q_{\rm f} + q_{\rm b} = 0$. The system is indeed not a thermal diode.
cpl-38-3-030501-fig3.png
Fig. 3. (a) Schematic for a 2D closed system with interior circulating mass flow (green arrows). Its two ends are maintained at constant temperatures $T_{1}$ (red line) and $T_{2}$ (blue line). (b) A simplified model for (a), where the velocity is assumed to have constant magnitude $v_{0}$. The regions adjacent to the thermal reservoirs (marked by dashed lines) are assumed to be isothermal. The channel is thereby separated into two paths p1 and p2. (c) The paths can be made asymmetric by placing both thermal reservoirs at the bottom.
The geometric symmetry can also be broken by placing the two thermal reservoirs as in Fig. 3(c). The length of p1 becomes $L_{1}$, while that of p2 becomes $L_{2}$. The same argument gives $$\begin{align} Q=\rho c_{\rm p} v_{0} \sigma (T_{1}-T_{2})\Big({\frac{1}{e^{v_{0} L_{1} /D}-1}+\frac{e^{v_{0} L_{2} /D}}{e^{v_{0} L_{2} /D}-1}} \Big).~~~~~ \tag {21} \end{align} $$ Although the coefficient changes, the heat is still linearly proportional to the temperature difference, so $q_{\rm f} + q_{\rm b} = 0$. The system is not a thermal diode irrespective of the geometry. Our argument is numerically confirmed using simulations on the model in Fig. 4(a). The channel is in rectangular shape of size $8.9{\,\rm cm} \times 1.9{\,\rm cm}$. Its width and thickness are both 3 mm. The corners are rounded with radius 4.5 mm. The velocity field ${\boldsymbol v}$ is generated with an interior fan with flow speed $v_{0}$ placed at the center of the upper side. Using frictionless boundary conditions, the simulated profile of magnitude $|{\boldsymbol v}|$ in Fig. 4(a) shows a uniform velocity magnitude except for the corners, so our assumption of a uniform velocity magnitude is a reasonable approximation. Two thermal reservoirs of size $1{\,\rm cm} \times 0.3{\,\rm cm}$ are attached at the bottom of the system, as shown in Fig. 4(b). The temperature difference between them, $\Delta T=T_{1} - T_{2}$, is set as 10 K for the forward case and $-$10 K for the backward case. The average temperature is fixed at $T_{0} = (T_{1} + T_{2})/2 = 300$ K. The lengths of the two paths separated by the two thermal reservoirs are $L_{1} = 6$ cm and $L_{2} = 12.8$ cm. The temperature profiles in Fig. 4(b) are obtained with $v_{0}$ setting as $v_{0}=D/L_{1}$. Although the forward and backward profiles look asymmetric, the heat fluxes $q_{\rm f}$ and $q_{\rm b}$ are exactly opposite such that $q_{\rm f} + q_{\rm b} = 0$, as shown in Fig. 4(c). It is also confirmed that the theoretical results calculated using the simplified model Eq. (21) meet well with the simulated results.
cpl-38-3-030501-fig4.png
Fig. 4. (a) Profile of the velocity magnitude $|{\boldsymbol v}|$ in a 2D closed system driven by an interior fan (green arrows). (b) Temperature profiles for the forward (upper) and backward (lower) cases. (c) Dependence of the heat fluxes entering (leaving) the system on $v_{0}L_{1}/D$.
We are in a position to provide a general proof for the invalidity of using pure advection to make thermal diodes. The proof is based on an equivalence between thermal diode effect and broken global reciprocity at steady state, which was demonstrated in a recent work.[4] The global reciprocity is defined for a two-port system (${\boldsymbol r} \in V$, ${\boldsymbol r}$ is position vector) in contact with two solid blocks (${\boldsymbol r} \in V_{1}$ and $V_{2}$). The interfaces between the system and the blocks are small enough to ignore temperature variance on them. We can define Green's function $G({\boldsymbol r}|{\boldsymbol r} _{1})$ for the entire region $V\bigcup V_{1}\bigcup V_{2}$, which is the solution to $$ -\nabla \cdot [\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})]+\rho c_{\rm p} \boldsymbol{v}\cdot \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})=\delta (\boldsymbol{r}-\boldsymbol{r}_{1}),~~ \tag {22} $$ where $\delta ({\boldsymbol r})$ is the Dirac delta function and ${\boldsymbol v}$ vanishes on $V_{1}$ and $V_{2}$. The system preserves steady-state global reciprocity if $$ G(\boldsymbol{r}_{2} \vert \boldsymbol{r}_{1})=G(\boldsymbol{r}_{1} \vert \boldsymbol{r}_{2}),~~ \tag {23} $$ where ${\boldsymbol r}_{1} \in V_{1}$ and ${\boldsymbol r}_{2} \in V_{2}$. It was proved that when steady-state global reciprocity is preserved, the system is not a thermal diode, assuming that there is no interior heat generation. In the following we prove Eq. (23). It should be noted that the local reciprocity as well as the global reciprocity at nonzero frequency is generally broken in pure advection. The effects of rotating fluids[29,30] and solid rings[31–34] have been detailedly studied recently. The governing equation for $G({\boldsymbol r}|{\boldsymbol r}_{2})$ is $$ -\nabla \cdot [\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})]+\rho c_{\rm p} \boldsymbol{v}\cdot \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})=\delta (\boldsymbol{r}-\boldsymbol{r}_{2}).~~ \tag {24} $$ We multiply Eq. (22) with $G({\boldsymbol r}|{\boldsymbol r}_{2})$, and Eq. (24) with $G({\boldsymbol r}| {\boldsymbol r}_{1})$. Then perform integral on their difference over the entire region $$\begin{alignat}{1} &\int_{V_{1} \cup V_{2}}\Big\{G(\boldsymbol{r}\vert \boldsymbol{r}_{1})\nabla \cdot [\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})]-G(\boldsymbol{r}\vert \boldsymbol{r}_{2}) \\ &\cdot \nabla \cdot[\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})]\Big\} dV +\int_V \Big\{G(\boldsymbol{r}\vert \boldsymbol{r}_{1})\nabla \cdot [\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})]\\ &-G(\boldsymbol{r}\vert \boldsymbol{r}_{2})\nabla \cdot [\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})]\Big\} dV\\ &+\int_V \rho c_{\rm p} \big[G(\boldsymbol{r}\vert \boldsymbol{r}_{2})\boldsymbol{v}\cdot \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})\\ &-G(\boldsymbol{r}\vert \boldsymbol{r}_{1})\boldsymbol{v}\cdot \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})\big] dV=G(\boldsymbol{r}_{2} \vert \boldsymbol{r}_{1})-G(\boldsymbol{r}_{1} \vert \boldsymbol{r}_{2}).\\~~ \tag {25} \end{alignat} $$ The third term on the left-hand side (LHS) is simplified to an integral on $V$ since ${\boldsymbol v}$ vanishes outside. We can also only integrate over the system as follows: $$\begin{align} &\int_V \Big\{ G(\boldsymbol{r}\vert \boldsymbol{r}_{1})\nabla \cdot [\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})] -G(\boldsymbol{r}\vert \boldsymbol{r}_{2})\\ &\cdot \nabla \cdot[\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})]\Big\}dV +\int_V \rho c_{\rm p} \big[G(\boldsymbol{r}\vert \boldsymbol{r}_{2})\boldsymbol{v}\cdot \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})\\ &-G(\boldsymbol{r}\vert \boldsymbol{r}_{1})\boldsymbol{v}\cdot \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})\big] dV=0.~~ \tag {26} \end{align} $$ The right-hand-side (RHS) is zero since ${\boldsymbol r}_{1}$ and ${\boldsymbol r}_{2}$ are outside $V$. Combining Eqs. (25) and (26) gives $$\begin{alignat}{1} &\int_{V_{1} \cup V_{2}}\Big\{ G(\boldsymbol{r}\vert \boldsymbol{r}_{1})\nabla \cdot [\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})]-G(\boldsymbol{r}\vert \boldsymbol{r}_{2}) \\ &\cdot \nabla \cdot[\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})]\Big\} dV=G(\boldsymbol{r}_{2} \vert \boldsymbol{r}_{1})-G(\boldsymbol{r}_{1} \vert \boldsymbol{r}_{2}).~~~ \tag {27} \end{alignat} $$ Assuming Dirichlet, Neumann, or mixed boundary conditions on the boundaries of $V_{1}$ and $V_{2}$ except for their interfaces with $V$ ($S_{1}$ and $S_{2}$), the LHS becomes $$\begin{alignat}{1} &\int_{V_{1}\cup V_{2}}\{\nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})\cdot [\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})]-\nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})\\ &\cdot [\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})]\} dV +\int_{S_{1} \cup S_{2} } \{G(\boldsymbol{r}\vert \boldsymbol{r}_{2})[\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})]\\ &-G(\boldsymbol{r}\vert \boldsymbol{r}_{1})[\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})]\} \cdot \boldsymbol{n}dS \\ ={}&\int_{S_{1} } {G(\boldsymbol{r}\vert \boldsymbol{r}_{2})[\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})]} \cdot \boldsymbol{n}dS\\ &+\int_{S_{2} } {G(\boldsymbol{r}\vert \boldsymbol{r}_{2})[\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})]} \cdot \boldsymbol{n}dS \\ &-\int_{S_{1} } {G(\boldsymbol{r}\vert \boldsymbol{r}_{1})[\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})]\cdot \boldsymbol{n}dS} \\ &-\int_{S_{2} } {G(\boldsymbol{r}\vert \boldsymbol{r}_{1})[\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})]\cdot \boldsymbol{n}dS},~~ \tag {28} \end{alignat} $$ where ${\boldsymbol n}$ is the unit normal vector on the boundaries of $V_{1}$ and $V_{2}$. At steady state, the heat fluxes entering and leaving the system through $S_{1}$ and $S_{2}$ should be in balance, $$\begin{align} &\int_{S_{1} } {[\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})]} \cdot \boldsymbol{n}dS+\int_{S_{2} } {[\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{1})]} \cdot \boldsymbol{n}dS=0,\\ &\int_{S_{1} } {[\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})]} \cdot \boldsymbol{n}dS+\int_{S_{2} } {[\kappa \nabla G(\boldsymbol{r}\vert \boldsymbol{r}_{2})]} \cdot \boldsymbol{n}dS=0.~~ \tag {29} \end{align} $$ Since the areas of $S_{1}$ and $S_{2}$ are both small, the variance of Green's function on $S_{1}$ and $S_{2}$ is negligible. Therefore, the RHS of Eq. (28) vanishes, so does the RHS of Eq. (27). The steady-state global reciprocity is thus proved. The system is therefore not a thermal diode, meaning that if the two blocks are maintained at constant temperatures $T_{1}$ and $T_{2}$, the heat flux through the system will change sign after the values of $T_{1}$ and $T_{2}$ are swapped. The result can be extended to thermal reservoirs that have finite contact areas with the system, considering the linearity of the system. In conclusion, the effects of pure advection on the symmetry of heat transfer have been detailedly studied. It is proposed that a practical thermal diode should be closed to mass flux. Open and closed systems with pure advection are compared. The former shows a strong rectification effect, but the latter generally does not break the symmetry between forward and backward heat transfer. Therefore, pure advection in closed systems cannot be used to make a practical thermal diode. The argument is proved based on a steady-state global reciprocity for two-port systems. References Reciprocity in opticsExperimental Demonstration of Acoustic Chern InsulatorsNonreciprocity in acoustic and elastic materialsDiffusive nonreciprocity and thermal diodeHidden time-reversal symmetry in dissipative reciprocal systemsElectromagnetic NonreciprocityReciprocal Relations in Irreversible Processes. I.On Onsager's Principle of Microscopic ReversibilityTunable and switchable polarization rotation with non-reciprocal plasmonic thin films at designated wavelengthsNear-complete violation of detailed balance in thermal radiationNear-complete violation of Kirchhoff’s law of thermal radiation with a 03 T magnetic fieldRouting emission with a multi-channel nonreciprocal waveguideTransforming heat transfer with thermal metamaterials and devicesControlling macroscopic heat transfer with thermal metamaterials: Theory, experiment and applicationThermal diodes, regulators, and switches: Physical mechanisms and potential applicationsHeat Transfer in a Liquid Convective DiodeDesign of a Plane-Type Bidirectional Thermal DiodeNegative Thermal Transport in Conduction and AdvectionA Ubiquitous Thermal Conductivity Formula for Liquids, Polymer Glass, and Amorphous SolidsAn oxide thermal rectifierTemperature-Dependent Transformation Thermotics: From Switchable Thermal Cloaks to Macroscopic Thermal DiodesTemperature-dependent transformation thermotics for unsteady states: Switchable concentrator for transient heat flowNonlinear thermal conductivity of periodic compositesThermal Illusion of Porous Media with Convection-Diffusion Process: Transparency, Concentrating, and CloakingTemperature Trapping: Energy-Free Maintenance of Constant Temperatures as Ambient Temperature Gradients ChangeColloquium : Phononics: Manipulating heat flow with electronic analogs and beyondMacroscale Thermal Diode-Like Black Box with High Transient Rectification RatioThermal meta-device in analogue of zero-index photonicsTunable analog thermal materialAnti–parity-time symmetry in diffusive systemsHigh-Order Exceptional Points in Diffusive Systems: Robust APT Symmetry 2 Against Perturbation and Phase Oscillation at APT Symmetry BreakingA Continuously Tunable Solid‐Like Convective Thermal Metadevice on the Reciprocal LineEffective medium theory for thermal scattering off rotating structures
[1] Potton R J 2004 Rep. Prog. Phys. 67 717
[2] Ding Y, Peng Y, Zhu Y, Fan X, Yang J, Liang B, Zhu X, Wan X and Cheng J 2019 Phys. Rev. Lett. 122 014302
[3] Nassar H, Yousefzadeh B, Fleury R, Ruzzene M, Alù A, Daraio C, Norris A N, Huang G and Haberman M R 2020 Nat. Rev. Mater. 5 667
[4] Li Y, Li J, Qi M H, Qiu C W and Chen H S 2021 Phys. Rev. B 103 014307
[5] Silveirinha M G 2019 Opt. Express 27 14328
[6] Caloz C, Alù A, Tretyakov S, Sounas D, Achouri K and Deck-Léger Z L 2018 Phys. Rev. Appl. 10 047001
[7] Onsager L 1931 Phys. Rev. 37 405
[8] Casimir H B G 1945 Rev. Mod. Phys. 17 343
[9] Floess D, Chin J Y, Kawatani A, Dregely D, Habermeier H U, Weiss T and Giessen H 2015 Light: Sci. & Appl. 4 e284
[10] Zhu L and Fan S 2014 Phys. Rev. B 90 220301
[11] Zhao B, Shi Y, Wang J, Zhao Z, Zhao N and Fan S 2019 Opt. Lett. 44 4203
[12] Hu H, Liu L, Hu X, Liu D and Gao D 2019 Photon. Res. 7 642
[13] Li Y, Li W, Han T C, Zheng X, Li J X, Li B W, Fan S H and Qiu C W 2020 arXiv:2008.07964 [physics.app-ph]
[14] Yang S, Wang J, Dai G, Yang F and Huang J 2020 Phys. Rep. (in press)
[15] Wehmeyer G, Yabuki T, Monachon C, Wu J and Dames C 2017 Appl. Phys. Rev. 4 041304
[16] Jones G F 1986 J. Sol. Energy Eng. 108 163
[17] Chen K 1988 J. Sol. Energy Eng. 110 299
[18]Bejan A 2013 Convection Heat Transfer (New York: John Wiley & Sons)
[19] Xu L and Huang J 2020 Chin. Phys. Lett. 37 080502
[20] Qing X, Jinxin Z, Jixiong H, Xiangfan X, Tsuneyoshi N, Yuanyuan W, Jun L, Jun Z and Baowen L 2020 Chin. Phys. Lett. 37 104401
[21] Kobayashi W, Teraoka Y and Terasaki I 2009 Appl. Phys. Lett. 95 171905
[22] Li Y, Shen X, Wu Z, Huang J, Chen Y, Ni Y and Huang J 2015 Phys. Rev. Lett. 115 195503
[23] Li Y, Shen X, Huang J and Ni Y 2016 Phys. Lett. A 380 1641
[24] Dai G and Huang J 2020 Int. J. Heat Mass Transfer 147 118917
[25] Yang F B, Xu L J and Huang J P 2019 ES Energy & Environ. 6 45
[26] Shen X, Li Y, Jiang C and Huang J 2016 Phys. Rev. Lett. 117 055501
[27] Li N, Ren J, Wang L, Zhang G, Hänggi P and Li B 2012 Rev. Mod. Phys. 84 1045
[28] Huang S Y, Zhang J W, Wang M, Lan W, Hu R and Luo X B 2019 ES Energy & Environ. 6 51
[29] Li Y, Zhu K J, Peng Y G, Li W, Yang T, Xu H X, Chen H, Zhu X F, Fan S and Qiu C W 2019 Nat. Mater. 18 48
[30] Xu G, Dong K, Li Y, Li H, Liu K, Li L, Wu J and Qiu C W 2020 Nat. Commun. 11 6028
[31] Li Y, Peng Y G, Han L, Miri M A, Li W, Xiao M, Zhu X F, Zhao J, Alù A, Fan S and Qiu C W 2019 Science 364 170
[32] Cao P C, Li Y, Peng Y G and Qiu C W and Zhu X F 2020 ES Energy & Environ. 7 48
[33] Li J, Li Y, Cao P C, Yang T, Zhu X F, Wang W and Qiu C W 2020 Adv. Mater. 32 2003823
[34] Li J, Li Y, Wang W, Li L and Qiu C W 2020 Opt. Express 28 25894