Bulk-edge correspondence of classical diffusion phenomena

We elucidate that diffusive systems, which are widely found in nature, can be a new platform of the bulk-edge correspondence, a representative topological phenomenon. Using a discretized diffusion equation, we demonstrate the emergence of robust edge states protected by the winding number for one- and two-dimensional systems. These topological edge states can be experimentally accessible by measuring diffusive dynamics at edges. Furthermore, we discover a novel diffusion phenomenon by numerically simulating the distribution of temperatures for a honeycomb lattice system; the temperature field with wavenumber \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi $$\end{document}π cannot diffuse to the bulk, which is attributed to the complete localization of the edge state.


Bulk-edge correspondence of classical diffusion phenomena Tsuneya Yoshida * & Yasuhiro Hatsugai
We elucidate that diffusive systems, which are widely found in nature, can be a new platform of the bulk-edge correspondence, a representative topological phenomenon. Using a discretized diffusion equation, we demonstrate the emergence of robust edge states protected by the winding number for one-and two-dimensional systems. These topological edge states can be experimentally accessible by measuring diffusive dynamics at edges. Furthermore, we discover a novel diffusion phenomenon by numerically simulating the distribution of temperatures for a honeycomb lattice system; the temperature field with wavenumber π cannot diffuse to the bulk, which is attributed to the complete localization of the edge state.
In these decades, the notion of topology in condensed matter physics enhances its significance [1][2][3][4][5][6][7] . One of the characteristic topological phenomena is the emergence of robust gapless edge states due to topological properties in the bulk which is known as the bulk-edge correspondence; the chiral edge states emerge 8 corresponding to a finite value of the Chern number in the bulk of systems without symmetry 9 , which is elucidated in Ref. 10 . The topologically protected edge states are sources of novel phenomena, such as the quantized Hall conductance 9,11 , the emergence of Majorana fermions 12-18 , etc. Remarkably, recent works extended the bulk-edge correspondence to several classical systems which are governed by Maxwell equations, Newton equation, etc. [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35] . These progresses beyond quantum systems provide universal understanding from the topology and result in invention of new devises (e.g., the topological laser 36,37 ) thanks to the robust edge states. Therefore, further extending the bulk-edge correspondence beyond quantum systems is considered to be significant in term of both the scientific viewpoint and applications.
In this paper, we point out that classical diffusive systems can be a new platform of the bulk-edge correspondence, which highlights topological aspects of the classical diffusion phenomena; the diffusive systems include a wide variety of systems (e.g., thermal diffusion 38,39 , diffusion of impurities in metals 40 , diffusion of droplets of inks in water, etc.). To this aim, we discretize the diffusion equation based on Fick's law. The discretized diffusion equation allows us to discuss the bulk-edge correspondence of diffusion phenomena for classical systems; the governing equation is expressed in a matrix form that is mathematically equivalent to a tight-binding model of a quantum system. Our numerical data verify the bulk-edge correspondence for diffusion phenomena in the classical systems. Furthermore, our numerical simulation of the temperature distribution elucidates a novel diffusion phenomenon for a honeycomb lattice system; the temperature field with wavenumber k x = π cannot diffuse into the bulk, which is attributed to the complete localization of the edge state with k x = π . Here, we stress that systems we discuss are classical systems in contrast to a previous work 41 analyzing topology of diffusion of electrons.
Before addressing the discretization, let us briefly review Fick's law and the diffusion equation of a continuum scalar field φ(t, x) in one dimension where ∂ t(x) denotes derivative with respect to time t (spatial coordinate x). Here, depending on which system we consider, the scalar field φ(t, x) corresponds to the field of temperatures, the density of the diffusing material, etc.. Fick's law indicates that the corresponding flux J is given by J = −D∂ x φ(t, x) , where D is the diffusion coefficient. By combining this equation and the equation of continuity ∂ t φ(t, x) + ∂ x J(t, x) = 0 , we obtain the diffusion Eq. (1). Now, let us discretize the diffusion Eq. (1) connecting the diffusion phenomena to tight-binding models of quantum systems. In order to show the essential idea, we focus on one-dimensional systems.
Consider a system composed of two sites where the values of the discretized field φ 0 and φ 1 are assigned at each site (see Fig. 1a); for the heat conduction equation, consider two balls (e.g., macroscopic iron balls) where temperatures are T 0 and T 1 . Recalling Fick's law, we can write the flux flowing from site 0 to 1 with φ's, Here, we have chosen the distance between the sites as the unit of length. Thus, the timeevolution of the vector � φ = (φ 0 , φ 1 ) T is described by Therefore, for a one-dimensional chain composed of L x sites (see Fig. 1b), the time-evolution of the vector . To see this, we diagonalize the matrix Ĥ and focus on the long-wavelength limit. By applying the Fourier transformation, x , meaning that the time-evolution is described by Eq. (1) in the long-wavelength limit. Here we have used the correspondence k x ↔ −i∂ x .
In the above, by discretizing the diffusion equation, we have shown that the diffusive dynamics of classical systems can be described by the tight-binding model of quantum systems [see Eq. (3)]. This result implies that the diffusive systems serve as a new platform of topological physics beyond quantum systems.

SSH model of the heat conduction equation
In order to demonstrate that the diffusive dynamics of classical systems indeed shows topological phenomena we analyze a one-dimensional system with dimerization (see Fig. 2a) which corresponds to the Su-Schrieffer-Heeger (SSH) model 42,43 of quantum systems. In the rest of this paper, we discuss the discretized version of the heat conduction equation for the sake of concreteness.
Let us consider a one-dimensional system illustrated in Fig. 2a. The temperature at each site is described by the following vector, Firstly, let us discuss topological properties in the bulk. Namely, temporarily abandoning the boundary condition illustrated in Fig. 2a, we first impose the periodic boundary condition. In the momentum space, the matrix Ĥ SSH is rewritten as with k x = 2πn x /L x and n x = 0, 1, . . . , L x − 1 . Here, the corresponding vector denoting the temperature is written The Pauli matrices σ 's act on the sublattice degrees of freedom [ σ 1 = 0 1 1 0 , Before analyzing the topological properties, we note that the system shows a gap and preserves the chiral symmetry. Diagonalizing the matrix, we obtain the spectrum This result indicates that the spectrum shows a gap for δ = 1 .
The system also preserves the chiral symmetry; ĥ ′ SSH :=ĥ . Here, we note that the shift described by the identity matrix σ 0 does not affect the eigenvalue problem, meaning that topological properties of the eigenvectors are encoded into ĥ ′ SSH . Because ĥ ′ SSH shows the gap and preserves the chiral symmetry, it may possesses the topologically nontrivial properties which are characterized by the winding number: The winding number counts how many times the off-diagonal element of ĥ SSH winds around the origin of the complex plane, and thus, it takes an integer 44 . Computing the winding number, we can see that the winding number takes the value one ( W = 1 ) for 0 ≤ D ′ < 1 while it takes the value zero ( W = 0 ) for 1 ≤ D ′ .
For one-dimensional quantum systems with chiral symmetry, the winding number predicts the number of the gapless edge modes localized around the edges, which is typical example of the bulk-edge correspondence. We show that the bulk-boundary correspondence can be observed in our classical system. Figure 2b shows the spectrum of Ĥ SSH under the fixed boundary condition. This figure indicates that corresponding the winding number W = 1 ( W = 0 ), there exists an edge state (no edge state) localized at each edge, which is represented as a blue dot for each value of D ′ . Here, the edge state appears at ǫ = D + D ′ because of the term proportional to the identity matrix.
The above results demonstrate that the diffusive dynamics of classical systems exhibit the bulk-edge correspondence which is a unique topological phenomenon. We stress that the edge states appear due to the topological properties in the bulk, which implies the following behaviors as is the case of quantum systems 10 : (1) The edge states survive even in the presence of perturbation (2) A states is localized at the boundary of two SSH models δ > 1 and δ < 1 which would results in the essentially same behavior as the one shown in Fig. 2c. We note that imposing the free boundary condition breaks the chiral symmetry, which shifts the eigenvalue of the edge state away from D(1 + δ).

How to experimentally access the edge states
So far, we have shown that the edge states emerge at ǫ = D + D ′ because of the topological properties in the bulk. In the following, let us discuss how to experimentally access the edge states. and L x = 240 . We set the initial state as � T iα = δ i0 δ αA . These data are obtained by imposing Dirichlet-type boundary conditions. Namely, for the spatial direction, we imposed the fixed boundary condition. For timedirection, we imposed the initial condition T(0) iα = δ i0 δ αA . In either case of parameters, the temperature decays monotonically, which is compatible with the boundary condition of time � T(∞) iα = 0.
Scientific Reports | (2021) 11:888 | https://doi.org/10.1038/s41598-020-80180-w www.nature.com/scientificreports/ One possibility is to observe the time-evolution of the temperature at the edge which is considered to decay exponentially T 0A ∼ e −(D+D ′ )t . In Fig. 2c, the time-evolution of the temperature at edge (i x , α) = (0, A) is plotted. The temperature T 0A shows exponential decay for t 2τ with the half-life τ = 1/(D + D ′ ) = 0.83 for (D, D ′ ) = (1, 0.2) due to the edge state, while it deviates from the line of the exponential decay around t = 0.5 which is shorter than the half-life for (D, D ′ ) = (0.2, 1) . The above behaviors should be observed even in the presence of the disorder preserving the chiral symmetry. This is because regardless of the details of the system, a finite value of the winding number predicts edge states localized around the boundary as is the case of quantum systems. Therefore, we conclude that observing the time-evolution allows us to experimentally access the edge states induced by the bulk topological properties. We note that the time-evolution of the temperature at each site has been measured in Ref. 38 for continuous systems; when the system is composed of aluminum, the half-life τ is estimated to be τ ∼ 1ms (for more details, see Sec. IIB of Supplemental Material).
We also consider that at least in principle, the eigenvectors and eigenvalues of the matrix Ĥ SSH can be extracted from the experimental data in the following procedure. (1) Prepare a set of initial conditions � T (i) (t = 0) l ( l = 0, . . . , L x − 1 ) which are linearly independent each other; for instance, such initial conditions can be prepared by heating at a site. (2) Observe the temperature T l ] i x α = T lα δ li x δ α l α with α l = A, B and T lα taking a random value between 0.5 and 1. (Concerning the temperature, specific choice of the unit does not matter because the ratio of the temperature is discussed throughout this paper.) The eigenvalues e −ǫ ′ t 0 almost reproduce the ones of Ĥ SSH . We note that the deviation for 2.5 ǫ 3 is due to the rounding error; the matrix elements of T (t 0 ) ji exponentially decay. Figure 3b shows the edge state

Honeycomb lattice system
Topological phenomena of the diffusive dynamics can also be found for two-dimensional systems. To show this, we analyze a honeycomb lattice system illustrated in Fig. 4a. We have supposed that the sites are coupled with the diffusion coefficient D. The dynamics of the temperature at each site T is described by ∂ t � T = −Ĥ honey � T . As is the case of the SSH model, Ĥ honey corresponds to the honeycomb lattice of the tight-binding model; Ĥ ′ honey :=Ĥ honey − 3D1l preserves the chiral symmetry. Under the periodic (fixed) boundary condition for the x-(y-) direction, the system can be regarded as a set of one-dimensional systems aligned along the momentum space −π ≤ k x < π . Noting that the one-dimensional system specified by k x preserves the chiral symmetry, we can compute the winding number; the winding number takes the value one ( W = 1 ) for 2π/3 < |k x | < π , while it takes the value zero ( W = 0 ) for 0 ≤ |k x | < 2π/3 . Correspondingly, only for 2π/3 < |k x | < π , the edge state appears 13,45 . We note that along the armchair edge, no edge states can be observed. For more details of the spectrum, see Sec. IIIA of Supplemental Material.
The presence or absence of the edge state can affect the diffusive dynamics. Figure 4b,c shows the timeevolution at site i cz ( i az ). Figure 4b shows the dynamics obtained for the two cases of the initial condition spatially modulating either k x = 0 or k x = π . The data of k x = π are obtained by subtracting data obtained with the initial condition 2 T iz1 from the ones obtained with T iz2 (For specific form of T iz1 and T iz2 , see Sec. IIIB of Supplemental Material). Figure 4b indicates that at the zigzag edge, the temperature field with k x = π exponentially decays T icz ∼ e −3Dt while the data of the temperature field with k x = 0 deviates from e −3Dt . The above time-evolution is consistent with the presence of the edge state for 2π/3 < |k x | < π whose eigenvalue is 3D. We note that the time-evolution of the armchair edge deviates from e −3Dt for either initial condition (see Fig. 4c).
Furthermore, the edge state at k x = π results in counter intuitive dynamics; for the zigzag edge, the initial state with k x = π cannot diffuse to the bulk (see Fig. 5a) while for the armchair edge, the initial state diffuses to the bulk (see Fig. 5b). This intriguing behavior is due to the complete localization of the edge state with k x = π . The above counterintuitive behavior is due to the complete localization of the state around the zigzag edge.

Summary
In this paper, we have elucidated the topological aspect of the diffusive dynamics, providing a new platform of the bulk-edge correspondence.
Specifically, based on Fick's law, we have introduced the discretized form of the diffusion equation, bridging the diffusive dynamics of classical systems and a tight-binding model discussed for quantum systems. The correspondence between the classical and quantum systems allows us to discuss the topological phenomena (e.g., the bulk-edge correspondence) for the diffusive dynamics of classical systems; we have numerically elucidated that topological properties characterized by the winding number in the bulk induces the edge states for the onedimensional system and the honeycomb lattice system. Furthermore, our numerical simulation has revealed a novel diffusion phenomenon for the honeycomb lattice system; at zigzag edges, the temperature field with spatial modulation k x = π cannot diffuse to the bulk.
Our results provide topological insights into diffusion phenomena, indicating the potential existence of diffusion phenomena analog of topological insulators for other symmetry classes and higher-order insulators. Their realization is left as future works to be addressed.  T(t = 1)/T 0 for a zigzag (an armchair) edge. The initial condition is chosen as T iz3 ( T ia3 ) for data of zigzag (armchair) edges, which allows us to observe the mode with k x = π ( k y = π ) for the zigzag (armchair) edge. For more details of the initial condition, see Fig. 2 and Sec. IIIB of Supplemental Material. The data shown in panel (a), (b) are obtained under the periodic and fixed (fixed and periodic) boundary conditions along the x-and y-directions. Here, we have taken T 0 = 0.0497 (0.0292) for data of zigzag (armchair) edges. The data are obtained for L x = L y = 40 and D = 1.