Selectable diffusion direction with topologically protected edge modes

Topological insulators provide great potentials to control diffusion phenomena as well as waves. Here, we show that the direction of thermal diffusion can be selected by the contributions of the topologically protected edge modes via the quantum spin Hall effect in a honeycomb-shaped structure. We demonstrate that when we set our structure to the temperature corresponding to the type of edge mode, the direction of thermal diffusion can be tuned. Moreover, this diffusion system is found to be immune to defects owing to the robustness of topological states. Our work points to exciting new avenues for controlling diffusion phenomena.

Recent theoretical and experimental studies have demonstrated that topological edge modes can be applied to diffusion systems using the SSH model [30][31][32][33].Such edge modes enable temperature localization and robust thermal decay.The pioneering studies implicate the potential for the diffusion control at higher dimensions using alternative topological states.However, topological states with higher DOF have not been investigated sufficiently yet in diffusion systems.
In this paper, we demonstrate that QSHE-based topologically protected edge modes appear in a thermal diffusion system.Our structure consists of honeycomb-shaped unit cells with topological and ordinary states.We consider that the heat transfer in our structure at the topological edge modes appears around the boundary between the topological and ordinary states.From numerical and analytical studies, we show that the temperature corresponding to the edge modes provides the directional heat transfer.The honeycomb-shaped unit cells function as a heat-transfer path while maintaining the edge modes in our structure.Consequently, the diffusion direction can be selected based on the type of excited edge mode.We also verified a well-known unique characteristic of topological edge modes, which is that they are immune to defects.Our results indicate that the use of QSHE-based topological edge modes has the potential to control heat transfer in any direction.Generally, our work should motivate systematic studies to apply topological properties to all diffusion phenomena.

A. Design of a topological diffusion system
Our structure consists of periodically aligned honeycomb-shaped unit cells, as illustrated in Fig. 1(a).The unit cell has six circle sites.The nearest neighboring sites and unit cells are connected by fine beams with the effective diffusivities D 1 and D 2 .As a result, the equation of thermal diffusion in a unit cell is expressed by the 6 × 6 effective diffusivity matrix, where Here, eigenfunctions, T 1 , . . ., and T 6 , denote the temperature at each site as numbered in Fig. 1(a), k is the wavenumber, and a 1(2) is the unit vector as shown in Fig. 1(a).We design the topological and ordinary unit cells by adjusting the ratio r = D 1 /D 2 of the two effective diffusivities.The QSHE-based topological phase transition can be controlled via r [2,16].Diagonalizing the effective diffusivity matrix in Eq. ( 1), we obtain the spectrum of the eigenvalues ϵ.To identify the state of the structure with r > 1 and r < 1, we confirm the eigenfunctions of each doubly degenerated mode specified as i-iv in Fig. 1(b) and i'-iv' in Fig. 1(d).The eigenfunctions in these modes correspond to the site temperature in the unit cell.Figure 1(e) shows the temperature distributions of modes i-iv for r > 1.We observe the dipole p x and p y modes at modes i and ii, and quadrupole d xy and d x 2 −y 2 modes at modes iii and iv, respectively.Thus, the lowly and highly polarized modes appear bellow and above the band gap, respectively.This result indicates that the structure with r > 1 is the ordinary state.In contrast, in Fig. 1(f), for r < 1, modes i' and ii' (iii' and iv') show the quadrupole (dipole) modes.Such an inversion of the order between the dipole and quadrupole modes signifies that the structure with r < 1 is the nontrivial topological state.To demonstrate the QSHE-based topologically protected edge modes in the thermal diffusion system, we consider a supercell with a boundary between the topological and ordinary states, as shown in Fig. 2(a).Figure 2(b) shows the spectra of the supercell.The red and blue lines show the topological edge modes with different polarizations in the unit cells.The gray lines indicate the bulk modes.We focus on the two band-gap-crossing edge modes at k x = 0, indicated by the magenta and cyan arrows in Fig. 2(b) and are specified as modes A and B, respectively.
Figures 2(c) and 2(d) depict the temperature distributions in modes A and B in the supercell.Both temperature distributions are localized around the boundary.The temperatures in modes A and B correspond to the eigenfunction amplitudes obtained by solving the eigenvalue equation, which consists of thermal diffusion equations for all 120 sites in the supercell.Such edge modes have great potential for controlling diffusion phenomena.Indeed, in wave systems, topologically protected edge modes have intriguing characteristics, such as field localization and unidirectional wave propagation [2,8,11,12,16,24].
Based on the eigenfunction amplitudes on each site, we analyze heat transfer in modes A and B. Figures 2(e) and 2(f) visualize the diffusion direction and quantity of heat transferred between the sites in the two ordinary and topological unit cells (specified as O2, O1, T1, and T2).The green arrows denote the flow direction, and their widths are the transferred heat quantity between the nearest neighboring sites.
Focusing on the heat balance in the individual unit cells, we find that modes A and B have a unique diffusion direction.We calculate the heat balances by summarizing the heat quantities along the x-and y-axes [right panels in Figs.2(e) and 2(f)].In Fig. 2(e), the heat balance indicates the thermal flow along the x-axis through the unit cell; there is no heat balance along the y-axis.Thus, mode A transfers heat only along the direction parallel to the boundary.On the other hand, as illustrated in Fig. 2(f), mode B provides thermal flow only along the y-axis.Therefore, our structure has the potential to select the diffusion direction by exciting the different edge modes.

B. Edge states and heat transport
To verify our theoretical prediction for selecting the diffusion direction, we compute the time evolution of the temperature in modes A and B using COMSOL Multiphysics.Figure 3(a) shows the model used in the calculation.Each of the half structures in Fig. 3(a) consists of topological (T) or ordinary (O) states.For the numerical investigation, we convert the eigenfunction to temperature in Kelvin.Here, T 0 = 293.15K is the reference temperature, α = 100 is the amplification coefficient, T , and s ∈ {O, T} is the ordinary and topological sites.To excite the edge modes to the system, we apply modes A or B to the unit cells of m y = 1 and m y = 2, which are enclosed by the broken black line in Fig. 3(a).The other sites are set as T 0 because the eigenfunction amplitudes are negligibly small in the unit cells of m y ≥ 3, due to the strong localization of the field around the boundary at the edge modes.In fact, we observe that the temperature at the site of m y = 3 (T s,mx,3,n ) is 30% or less of m y = 1 (T s,mx,1,n ) for all m x and n values.
Figures 3(b) and 3(c) show the numerical time evolution of the temperature distributions for modes A and B in the region enclosed by the solid black square in Fig. 3(a).When exciting mode A, we observe high and low temperatures at t = 40 s at left and right edges [regions X left and X right in Fig. 3(a)], see Fig. 3(b).Interestingly, in the topological diffusion system, another edge mode enables heat transfer along the direction perpendicular to the boundary, unlike the topological wave systems.When exciting mode B, we indeed observe high and low temperatures at above and below edges [regions Y above and Y below in Fig. 3(a)], see Fig. 3(c).Thus, modes A and B realizes thermal diffusion only along the x-and y-axis, respectively, as predicted in Fig. 2

(e) and 2(f).
Considering the antiphase of the edge modes, we can further select the diffusion direction.As the phase of the fields is time-invariant in diffusion systems, inphase mode A (B) and antiphase mode A ′ (B ′ ) provides symmetric temperature distributions about the y-(x-)axis.When the temperature distribution of the inphase and antiphase modes coincides with the symmetry of the structure, the eigenequation contains antiphase modes as the eigenvalue solution.Because our structure is axisymmetric about the y-axis, mode A ′ can appear as an edge mode.Figure 4(a) shows the eigenfunctions in cells O1, O2, T1, and T2 of the supercell in mode A (light blue) and A ′ (dark blue).Mode A ′ has an eigenvalue at k x = 0.02, as shown by the blue line in Fig. 2(b).The eigenfunction of mode A ′ exhibits an inverted distribution of mode A at site O2 (also at T1, O2, and T2).Thus, mode A and A ′ are mutually antiphase and this mode property is matched to the symmetry of the structure.
Importantly, the inverted initial field about the y-axis affects the diffusion direction.The inset in Fig. 4(a) shows the snapshot of the temperature distribution at t = 40 s when mode A ′ is excited to the structure in Fig. 3(a).We clearly see high and low temperatures in regions X right and X left , respectively.This proves that mode A ′ inverts the diffusion direction of mode A. However, mode B ′ does not appear in the solutions of the eigenequation based on our supercell because our structure is asymmetry about the x-axis.Thus, we can change the diffusion direction by selecting the initial temperature of the sites based on the three different modes (A, A ′ , and B) in our structure.Furthermore, our structure enables to control the diffusion direction owing to the incoherence of the edge modes.Specifically, the mutually incoherent modes A (A ′ ) and B can be linearly combined.Using the ratio r T of modes A (A ′ ) and B, we can excite our structure in the combined mode as where T A ′ s,v is the eigenfunction of mode A ′ and 0 ≤ r T ≤ 1. Figure 4(b) shows the temperature dependence of r T in regions X left , X right , Y above , and Y below in Fig. 3(a) at t = 40 s.Depending on the type of the combined mode, the temperature in the four regions linearly increases and decreases with an increase in r T .In particular, when r T = 0.5, we obtain a symmetric temperature distribution about the y-axis, see the inset of Fig. 4

C. Robustness of the edge modes
The edge modes gradually collapse with time after those are applied to the structure.As directional heat transfer of the unit cells continues as long as the edge modes are maintained, the temporal robustness of the edge modes is crucial for designing topological diffusion systems.Focusing mode A, we evaluate the temperature variation with time for the topological unit cell.Figure 5(a) shows the relative temperature R n (t) = TT,n (t)/ TT,2 (t), considering the unbiased temperature TT,n (t) = (T T,7,1,n (t) − T 0 )/α of the unit cell near the boundary, (m x , m y ) = (7, 1) in Fig. 3(a).For t < 20 s, R n (t) on each site has a constant value, that is, each temperature uniformly decays at the topologically protected decay rate.For t ≥ 20 s, mode A collapses because the temperatures of the unit cell almost converge to T 0 , see inset of Fig. 5(a); thus mode A cannot be maintained.The temporal mode robustness contributes to the duration of the directional heat transfer.Figure 5(b) shows the temperature in regions X left (circles) and X right (crosses) as a function of time.For t < 20 s, in region X left (X right ), the temperature increases (decreases) with time due to the directional heat transfer provided by mode A. For t ≥ 20 s, the temperature in X left (X right ) decays (saturates) because the edge mode plays a minor role in the structure.
Finally, we evaluate a unique characteristic of the topologically protected edge modes, which is their immunity to defects.To introduce the defects into our structure, we remove the beams from sites 1 and 2, see Fig. 5(c), for two unit cells (s, m x , m y ) = (T, 6, 1) and (T, 8, 1).The time dependence of the temperature in regions X left and X right is plotted with red symbols in Fig. 5(b).We observe small variations in the temperatures with and without defects.Hence, the edge states are robust against disorders due to topological protection as well as other topological systems.

III. DISCUSSION
We have shown that QSHE-based topologically protected edge modes appear in the diffusion system consisting of the honeycomb-shaped structure.While previous topological edge modes in diffusion systems could only design the decay rate [30][31][32], the QSHE based edge states in our structure realize that the diffusion direction is selected based on the type of edge mode.In addition, it is found that the edge modes are immune to defects.The robust modes stably lead to directional heat transfer in realistic diffusion systems.The topological diffusion systems provide a fruitful avenue for temperature controls, thermal managements, and control for other diffusion phenomena with the topological insulators.
Again, our approach taken here exploits the topological edge modes in the honeycomb structure for the directional heat transfer.This can potentially be extended to other diffusion phenomena, e.g., ion transport.An interesting point is that the ionic transport is described by internal and external factors unlike the thermal diffusion.Specifically, ionic diffusivity is affected by the external factor (e.g., ambient temperature) as well as the internal material properties such as number of the mobile ions, charge, and activation energy.The ambient temperature can be adjusted by applying external signals.Hence, it will be possible actively to control the topological and ordinary states via the ambient temperature in ion transport systems.Design for topological systems depending on both internal and external factors is a crucial future work.

A. Numerical simulation of time evolving temperature
The time evolution of the temperature is calculated by the MEMS module in COMSOL Multiphysics.We consider the finite structure consisting of 13 × 12 unit cells to obtain the edge states between the ordinary and topological unit cells, see Fig. 3(a).Each unit cell has six disc-shaped sites (thickness of 10 mm and diameter of 20 mm).The cells have intra-(inter-)cell connection using the beam of L = 35 mm with the thermal diffusivity D 1 (D 2 ), see Fig. 1(a).We set the initial temperature to the unit cells, corresponding to the edge modes calculated from the eigen functions of the effective diffusivity matrix.We apply adiabatic boundary on whole surface of the structure to eliminate the thermal radiation and convection.In the numerical model, we design the thermal diffusivities of the beams by directly adjusting the value of thermal conductivity in the original material parameters.The temperature is calculated in the time duration from 0 to 100 s with the step of 1 s.
Figures 1(b)-(d) show the spectra for an infinite periodic honeycomb lattice with r > 1, r = 1, and r < 1.When all diffusivities in the unit cell have the same value (i.e., r = 1), we observe the Dirac cone at Γ point [Fig.1(c)].For r > 1 and r < 1, the spectra have a bandgap and two doubly degenerated modes [Figs.1(b) and 1(d)].

FIG. 3 .
FIG. 3. (a) Structure with the boundary between the topological and ordinary unit cells.The edge modes are excited in the region enclosed by the broken black line.Initial temperature is T0 = 293.15K in the area except the enclosed region.Regions X left (red line), X right (blue line), Y above (magenta line), and Y below (cyan line) are the nearest neighboring sites from the heated and cooled regions.The thickness of the simulated model is 10 mm, the diameter of the site is d = 20 mm, and the distance between the nearest neighboring sites is the same as the value for the model in Fig. 1 (L = 35 mm).(b), (c) Snapshots of thermal diffusion at t = 0, 20 s, and 40 s for modes (b) A and (c) B in the area enclosed by the solid black square in (a).

FIG. 5 .
FIG. 5. (a) Time dependency of the relative temperature Rn (t) on the six sites in the topological unit cell of (mx, my) = (7, 1) at mode A. The inset shows the temperatures in time on each site.(b) Time dependency of the temperatures in regions X left (circles) and X right (crosses).Red and black symbols show the results with and without the defect, respectively.(c) Schematic of the unit cell including the defects.Sites 1 and 2 lack beams.We use the same simulation settings as that in Fig. 3 except for the defects in (c).