Diffusive skin effect and topological heat funneling

Non-Hermitian wave system has attracted intense attentions in the past decade since it reveals interesting physics and generates various counterintuitive effects. However, in the diffusive system that is inherently non-Hermitian with natural dissipation, the robust control of heat flow is hitherto still a challenge. Here we introduce the skin effect into diffusive systems. Different from the skin effect in wave systems, where asymmetric couplings were enabled by dynamic modulations or judicious gain/loss engineering, asymmetric couplings of the temperature fields in diffusive systems can be realized by directly contacted metamaterial channels. Topological heat funneling is further presented, where the temperature field automatically concentrates towards a designated position and shows a strong immunity against the defects. Our work indicates that the diffusive system can provide a distinctive platform for exploring non-Hermitian physics as well as thermal topology. Non-Hermitian systems can house a range of unusual physical phenomena such as the skin-effect which has been typically observed in optical lattices and electrical circuits. Here, the authors show the non-Hermitian skin effect in thermal transport, demonstrating that the topologically protected heat flow can be realized by using thermal metamaterials.

N on-Hermitian Hamiltonian plays a core role to describe systems that exchange energy with the environment. In wave systems, the non-Hermitian physics has been explored intensely in the past decade [1][2][3][4][5] . For example, the paritytime symmetric system that operates at an exceptional point (EP) can realize unidirectionally invisible cloak 6 , single mode laser 7 , one-way mode switcher 8,9 , high-sensitivity sensor 10 , and wireless power transfer 11,12 . The dynamically modulated systems with a broken time-reversal symmetry can produce various novel gyrators, isolators and circulators 13 . Besides, an intriguing "skin effect" was recently discovered in non-Hermitian systems, where all the bulk states localize at the edges under the open boundary condition (OBC), which modifies the bulk-boundary correspondence in Hermitian cases [14][15][16][17][18][19][20][21][22] . The non-Hermitian skin effect shows application prospects in the nonreciprocal energy manipulation and has been demonstrated in various fields such as optical lattices 23 , mechanical systems 24 , quantum-walk networks 25 , and electrical circuits 26 , etc. In diffusive systems that are non-Hermitian inherently, the thermal metamaterials advanced rapidly as it enables various agile and flexible manipulation of heat flow 27,28 . Based on the theory of transformation thermotics, fruitful progresses have been achieved, such as thermal cloaks that conceal objects in heat conduction or radiative camouflages against infrared detection [29][30][31] , thermal diodes and concentrators that isolate and harvest heat directionally [32][33][34][35] . Since the diffusive system is dissipative, it provides a natural platform to investigate the non-Hermitian physics, for example, the observation of antiparity-time symmetric phase transition at an EP [36][37][38] . The Hamiltonian respecting anti-parity-time symmetry can be constructed by imposing the convection to diffusive systems, which generates dynamically "stopped" temperature fields 36 . However, the diffusive counterparts of asymmetric coupling and skin effect have not been discovered hitherto. Due to the symmetric energy exchanges in general cases, breaking the time-reversal symmetry with temporal modulations 23 and judicious gain/loss engineering of couplings 26 were considered the few-in-number choices.
In this work, we show that it is possible to realize asymmetric coupling and skin effect in the diffusive systems by directly contacting thermal metamaterials. Different from the cases in wave systems, the effective non-Hermitian Hamiltonians describing asymmetric couplings in diffusive systems are imaginary due to the dissipative nature, as shown in Fig. 1a. Our finding starts from the fact that heat and temperature variations during the energy exchange are not equivalent. Temperature evolution is actually asymmetric between different components. As a result, we realize the diffusive skin effect, which shows that the eigenstates in all bulk bands can become localized edge states under open boundaries. This is fundamentally different from the Hermitian case in which only the edge modes can have boundary localizations. Unlike the non-Hermitian skin effects in the wave systems, where the asymmetric couplings were implemented by temporal modulations or tailored gain/loss couplings, asymmetric coupling of temperature fields can be easily enabled by the directly contacted thermal metamaterials. We further propose an approach to construct a periodic non-Hermitian Hamiltonian in the parameter space from an aperiodic structure. On this basis, the skin-effect-induced heat funneling is showcased, where temperature fields concentrate to the designated position, being regardless of the initial condition and showing a strong immunity against the defects, as schematically shown in Fig. 1b. The robust temperature fields with nontrivial gradients can be useful for thermoelectric power generation or heat harvesting 39 . Moreover, the high sensitivity of diffusive skin effect to the change of boundary conditions makes it possible to achieve topological thermal sensing 40 .

Results
Directly asymmetric coupling. In order to understand the asymmetric coupling diffusive systems intuitively, we can start from a double-ring toy model that was experimentally demonstrated recently 36 . As shown in Fig. 2a, where two rings are vertically coupled in z direction through an interlayer. According to the Fourier's law in heat conduction, the coupling equations can be written as 36-38 where T 1 (T 2 ) is the temperature field of upper (lower) channel and D 1 ¼ κ 1 ) is the diffusivity of the upper (lower) channel with κ 1 (κ 2 ), ρ 1 (ρ 2 ) and C 1 (C 2 ) being the thermal conductivity, mass density and heat capacity. x is the position along the channel. The heat exchange rate of the upper (lower) channel is where κ i is the thermal conductivity of the coupling layer, b and d are the thicknesses of the ring channel.
As the temperature field in the ring structure is spatially periodic, we can assume that Eq. (1) has a wave-like solution of In the solution, A 1;2 is the amplitude, β is the propagation constant, ω 1;2 is the complex frequency, and T 0 is the reference temperature that is set to be zero for simplicity. Here we employ the temperature gradient and position of maximum temperature to represent the amplitude and phase of the wave-like solution. As all the ring channels have the same size, the propagating constants of circulating "heat wave" in them are equal to be β ¼ m 2π 2πR ¼ m R , where m is the mode order and R is the radius of the ring. Usually, only the first-order mode is taken into consideration, since it is very stable and can be selectively excited. Temperature variations in each channel can be regarded as uniform on condition that the layers are thin enough. As a result, the heat transfer in these channels possess coherent properties. Substituting the wave-like solution into Eq. (1) and considering the boundary continuity condition, we can deduce an effective Hamiltonian about A 1 ; A 2 À Á T 36-38 where S 1 ¼ Àðβ 2 D 1 þ h 21 Þ and S 2 ¼ Àðβ 2 D 2 þ h 12 Þ. In Eq. (2), the Hamiltonian of the diffusive system is imaginary, which is different from the one in the wave system (real), as shown in Fig. 1a. Moreover, when the material parameters ρ 1 C 1 ≠ρ 2 C 2 , the ring coupling becomes asymmetric with ih 21 ≠ih 12 . It needs to be pointed out that the main diagonal elements in Eq. (2) can be unified into iS 0 by adjusting the diffusivities D 1;2 along with h 21;12 , where ih 21;12 are replaced by ic 21;12 in Fig. 1a. Eigenvalues and eigenvectors of Eq. (2) are solved to be ω ± ¼ We can find the amplitudes in channels are asymmetric in the steady state, meanwhile, they are always time dependent with the factor of e Àiω ± t due to the system is dissipative. Therefore, our structure does not require dynamic modulation or gain-loss engineering to achieve asymmetric couplings.
Diffusive skin effect and non-Bloch band theory. For implementation, we can take the densities of directly coupled channels varying with a ratio of a 2 . Figures 2a and 2b display the asymmetric coupling temperature fields with different a, respectively. In order to realize the non-Hermitian diffusive Su-Schrieffer-Heeger (SSH) model, we can connect the asymmetric coupling unit cells with symmetric couplings in Fig. 2c 16 . Each unit cell consists of sublattices (A+B), where the asymmetric intracell couplings and symmetric intercell coupling are i h 1 ± δ À Á and ih 0 , respectively. For the open boundary condition, the Hamiltonian  is 16,20 whereÂ y l andB y l are the creation operators of A and B sites in the lth period. For the periodic boundary condition (PBC), the Blochmode Hamiltonian in momentum space is expressed into which respects the sublattice symmetry ofσ À1 z ðH PBC K ð Þ À iS 0Î Þσ z ¼ ÀðH PBC K ð Þ À iS 0Î Þ withσ x;y;z the Pauli matrices and K the Bloch vector 19 . Eigenvalues of Eq. (4) can be deduced where d x ¼ h 1 þ h 0 cosK and d y ¼ h 0 sinK. When the complex spectrum at the zero-energy level, viz., ω ± K ð Þ À iS 0 ¼ 0, we can obtain EPs at which all eigenstates degenerate (the sign '×' in Fig. 2d). However, the EPs deduced from the PBC case (h 1 ¼ 2:05) contradict the OBC one (h 1 ¼ 1:45), which indicates the bulk boundary correspondence breaks down (the complex spectrum of PBC case can be found in Supplementary Note 1) 16,20 . Back to the OBC case, where the eigen-equation and the Bloch phase factor e iK ! α ¼ re iK 16,20 . According to the non-Bloch band theory, the nearest neighbor coupling of temperature fields can be further transformed into where ϕ A;B ¼ α l ϕ A l ;B l are the eigenstates in lth unit cell and degenerate at ω l À iS 0 ¼ 0 (also at the EPs) with the generalized phase factor Therefore, EPs locate at the hyperbolic curve The corresponding amplitudes of eigenstates for PBC and OBC at the EP are shown in Fig. 2e and Fig. S2 (in Supplementary Note 2).
Here the EP is also a topological phase transition point and the topological phase diagram is shown in the "Methods".
Heat funneling effect. Figure 3a presents the heat funnel model with two mirrored SSH chains splicing at the designated interface. In our case, there are 5 periods for the left chain and 1 period for the right chain. The thicknesses of ring channels and interlayers are b and d. The internal and external radii of each channel are R 1 and R 2 with R 1 % R 2 . We define the propagation constant of temperature fields as β ¼ 1 The thermal conductivities, mass densities and capacities of thermal metamaterials for ring channels and coupling interlayers are ðκ n ; ρ n ; C 0 Þ and ðκ in ; ρ 0 ; C 0 Þ, where the channel number n indicates that the associated material parameters are varied. The coupling coefficients between channels n and n þ 1 can be expressed into ih n;nþ1 ¼ i κ in ρ nþ1 C 0 bd and ih nþ1;n ¼ i κ in ρ n C 0 bd for the forward and backward couplings, respectively. In order to keep the couplings in the non-Hermitian diffusive SSH model respecting the translation symmetry in parameter spaces, the parameters ρ n and κ in of thermal metamaterials should be arranged as ρ n ¼ a nÀ1 ρ 0 ; n ¼ 1; 3; 5; a n ρ 0 ; n ¼ 2; 4; 6; ( ; κ in ¼ a n κ i0 ; n ¼ 1; 2; 3; 4; ; which is shown by the setting of the heat funnel model in Fig. 3b. However, it is difficult to find the natural materials that satisfies these gradient parameters in implementation. According to the effective medium theory, we can realize the effective densities of channels and the effective thermal conductivities of coupling interlayers with composite metamaterials. In this way, the effective material parameters can be satisfied by modulating the doping rates of materials with highly contrasted properties (such as Cu and Polydimethylsiloxane). Here the asymmetric intracell couplings are ih 0 a and ih 0 a, with the symmetric intercell coupling ih 0 , where h 0 ¼ κ i0 ρ 0 C 0 bd . From the semi-infinite model in Fig. 2c, we can easily obtain h 1 ¼ 1 : However, this cannot lead to the equivalent tight-binding model in Fig. 3c, if the diffusivity takes the form of D n ¼ κ 0 ρ n C 0 (the diamond dots in Fig. 3d). The eigenfrequency of iS n ¼ Àiðβ 2 D n þ h n;nþ1 þ h nþ1;n Þ should be uniformed. Here we adjust D n in the way shown by the square dots in Fig. 3d to homogenize iS n into iS 0 ¼ Ài The Hamiltonian of non-Hermitian SSH model can be derived by replacing K with K þ ilna and the eigenvalues are deduced as ω ± K ð Þ À iS 0 ¼ ± i2h 0 cos K 2 . It is interesting to find that the system operates exactly at the topological transition point with a gapless band for any value of a. The generalized Bloch band is plotted in Fig. 3e, where the blue dots denote EPs.
We construct the 2D heat funnel model for simplicity in the theoretical analysis and the full wave simulations here. The structural parameters are chosen as b ¼ 5 mm, d ¼ 1 mm, and R 1 % R 2 ¼ 100 mm. The material parameters are set to be ρ 0 ¼ 1000 kg m À3 , C 0 ¼ 1000 J kg À1 K À1 , κ 0 ¼ 100 Wm À1 K À1 , and κ i0 ¼ 5 W m À1 K À1 . In Figs. 4a, c, we show the imaginary eigenvalues ÀIm ω and the eigenmodes ϕ n À Á at a ¼ 1; 0:4. The result reveals that the temperature gradients of all eigenmodes are almost evenly distributed in the channel array at a = 1. While at a = 0.4, the temperature fields will localize in the interface channels. For experiments, temperature gradient in the channel can be measured by the difference between maximum and minimum temperature values with an infrared camera. The normalized steady-state temperature gradient in each channel corresponds to the eigenmode solved from the Hamiltonian. As heat transfer is inherently dissipative and it normally needs much time for the temperature field evolving into the steady state, the initial temperature gradient should be given larger (for example, 273 K for the cold-side cooling and 320 K for the hot-side heating).
To testify the topological heat funneling, we study the temperature field evolutions in a symmetric coupling structure at a ¼ 1, where κ i0 is set to be 0:3Wm À1 K À1 . Imposing a random stimulation input, the temperature field will naturally concentrate toward the central channels, since the eigenmode of the lowest decay rate is excited and observed. In Fig. 5a, the simulated evolutions of temperature fields agree well with the tight-binding model, since the normalized temperature gradient of each channel is consistent with the theoretical analysis. However, for the asymmetric coupling case at a ¼ 0:4, the temperature field tends to funnel to the designed interface at channels 10 and 11 in Fig. 5b, regardless of initial conditions. Note that the heat funneling is robust against the defects with different coupling strengths, channel numbers, and interface locations (disorder analysis can be found in Supplementary Note 3). For example, it shows the robustness when channels 2 and 3 are set with an asymmetric coupling defect (a ¼ 2:5) in the model, which is different from the case in symmetric coupling systems where the temperature field localizes at the defect in Figs. 5c and 5d. The evolutions of temperature field along z axis are also presented in Figs. 5e and 5f, which validates the concept of heat funneling effect schematically displayed in Fig. 1b. Above discussions are based on the gapless system with equivalent intracell and intercell couplings. Without the loss of generality, we also investigate the gapped diffusive system with the intercell coupling much larger than asymmetric intracell couplings. In this case, the temperature field will concentrate at the channel 12, since the heat flow directions of the inverted structures change to be the same (Supplementary Note 4 and Fig. S4).

Conclusion
In summary, we investigate the physical mechanisms of diffusive skin effect and topological heat funneling. Unlike the cases in wave systems, the diffusive skin effect and heat funneling can be realized in a static framework such as directly contacted thermal metamaterials. The periodically driven Hamiltonian can be constructed in the parameter space with an aperiodic structure. Our results show that the diffusive system provides a distinctive platform to explore the topological physics and non-Hermitian dynamics. Our work is expected to inspire further exploration of other intriguing effects, including the higher-order diffusive skin effect, topological heat flow transfer, high-efficiency thermoelectric effect 39 , heat harvesting, and thermal sensing 40,41 , etc.

Methods
EPs in the topological phase. As the transformation factor and h 0 ¼ h 0 , from which we deduce The topological invariant (winding number W) is calculated by Figure 6 shows the topological phase diagram of the semi-infinite SSH model, where the EPs locate at the hyperbolic boundaries.
From Eq. (6) we also generate As h 1 þ δ ¼ h 0 a and h 1 À δ ¼ h 0 a, Eq. (11) can be simplified as and the coupling equation in real space is According to the diffusive non-Bloch band theory, Eq. (14) can be transformed into a standard form with phase factor e iK ! α ¼ re iK At the zero-mode, Obviously, EPs appear at δ ¼ h 1 where the system is in absolutely asymmetric coupling. The temperature field is exponentially localized at the right boundary with the rate α ± , when 0 < δ < h 1 . In Fig. S6 (Supplementary Note 6), we show the diffusive skin effect in the Hatano-Nelson model with densities of adjacent channels and thermal capacities of interlayers having the same gradient of a 2 .

Data availability
All relevant data are available from the corresponding author X.F.Z. upon request.

Code availability
The code is available from the corresponding author X.F.Z. upon reasonable request.