Effects of heterogeneity in site-site couplings for tight-binding models on scale-invariant structures

TL;DR

This paper investigates how heterogeneity in site couplings affects thermodynamic properties of quantum gases on scale-invariant structures, revealing heterogeneity-dependent oscillations and conditions for Bose-Einstein condensation.

Contribution

It analytically solves the energy spectrum of tight-binding models on scale-invariant structures and uncovers the impact of heterogeneity on thermodynamic behaviors.

Findings

Log-periodic oscillations in specific heat depend on heterogeneity.

Bose-Einstein condensation occurs only in non-homogeneous setups.

Heterogeneity influences low-dimensional quantum phenomena.

Abstract

We studied the thermodynamic behaviors of non-interacting bosons and fermions trapped by a scale-invariant branching structure of adjustable degree of heterogeneity. The full energy spectrum in tight-binding approximation was analytically solved . We found that the log-periodic oscillation of the specific heat for Fermi gas depended on the heterogeneity of hopping. Also, low dimensional Bose-Einstein condensation occurred only for non-homogeneous setup.

Equations

$$N_t=\frac{2}{2+\theta }[(\theta m+2m+1{)}^t+\theta +1].$$

$$\widehat{{\mathcal{H}}_0}=\sum _i{d}_i{\hat{a}}_i^{†}{\hat{a}}_i-\sum _{ij}{a}_{ij}{\hat{a}}_i^{†}{\hat{a}}_j.$$

$$\hat{\mathcal{H}}=-\sum _{ij}{t}_{ij}{\hat{a}}_i^{†}{\hat{a}}_j,$$

$$R(x)=\frac{\theta m+m+1}{\theta m+1}x-\frac{m}{(\theta m+1)x}$$

$$\rho (\epsilon )=\sum _{{\epsilon }^{\prime }\in J_R}\delta (\epsilon -{\epsilon }^{\prime })\underset{t\to \infty }{\lim }\frac{{\text{deg}}^{(t)}({\epsilon }^{\prime })}{N_t}$$

$$I(\epsilon )={\int }_{-1}^{\epsilon }\rho (x)dx.$$

$$I(\epsilon )\sim (\epsilon -{\epsilon }_0{)}^{\frac{d_s}{2}}.$$

$$\underset{\begin{array}{c}t\to \infty \\ \epsilon \to -1\end{array}}{\lim }\frac{N_{t+1}}{N_t}\frac{I(\epsilon )}{I\left(R(\epsilon )\right)}=1.$$

$$d_s=2\frac{\ln (\theta m+2m+1)}{\ln (\theta m+2m+1)-\ln (\theta m+1)}.$$

$$\bar{n}(ϵ)=\frac{1}{e^{\beta (ϵ-\mu )}+1}.$$

$${\int }_{-1}^1\bar{n}(ϵ)\rho (ϵ)dϵ=\gamma$$

$$1-2{\gamma }_c(m,\theta )=\underset{t\to \infty }{\lim }\frac{{\text{deg}}^{(t)}(\epsilon =0)}{N_t}.$$

$$C_{v,f}(T)=\frac{dU}{dT}=\frac{k_B{\beta }^2}{4}\left(<{ϵ}^2\rho (ϵ){>}_{{}_f}-\frac{<ϵ\rho (ϵ){>}_{{}_f}^2}{<\rho (ϵ){>}_{{}_f}}\right)$$

$$<g(ϵ){>}_{{}_f}={\int }_{-1}^1\frac{g(ϵ)dϵ}{{\cosh }^2\frac{\beta (ϵ-\mu )}{2}}.$$

$$C_{v,f}(T)=\frac{k_B{\beta }^2}{4}\left(<(ϵ-E_F{)}^2\rho (ϵ){>}_{{}_f}-\frac{<(ϵ-E_F)\rho (ϵ){>}_{{}_f}^2}{<\rho (ϵ){>}_{{}_f}}\right).$$

$$<g(ϵ){>}_{{}_f}={\int }_{E_F-\delta /\beta }^{E_F+\delta /\beta }\frac{g(ϵ)dϵ}{{\cosh }^2\frac{\beta (ϵ-E_F)}{2}}.$$

$$C_{v,f}(T)=\frac{k_B^{2}T}{4}{\int }_{-\delta }^{\delta }\frac{t^2\rho (t{k}_{B}T+{ϵ}_0+E_F)dt}{{\cosh }^2\frac{t}{2}}.$$

$$\bar{n}(ϵ)=\frac{1}{e^{\beta (ϵ-\mu )}-1}.$$

$$C_{v,b}=\frac{k_B{\beta }^2}{4}\left(<{ϵ}^2\rho (ϵ){>}_{{}_b}-\frac{<ϵ\rho (ϵ){>}_{{}_b}^2}{<\rho (ϵ){>}_{{}_b}}\right)$$

$$<g(ϵ){>}_{{}_b}={\int }_{-1}^1\frac{g(ϵ)dϵ}{{\sinh }^2\frac{\beta (ϵ-\mu )}{2}}.$$

$$T^{(t+1)}=\left(\begin{array}{cc}\frac{\theta m+1}{\theta m+m+1}{T}^{(t)}& S\\ S^T& 0\end{array}\right).$$

$$\frac{\theta m+1}{\theta m+m+1}{T}^{(t)}$$

$$S^T$$

$$\bigg{(}\frac{\lambda(\theta m+1)}{\theta m+m+1}T^{(t)}+SS^{T}-\lambda^{2}\bigg{)}\phi=0$$

$$(S{S}^T{)}_{ij}=\sum _{v_k\in V_{t+1}\backslash V_t}{S}_{ik}{S}_{jk}=\sum _{v_k\in V_{t+1}\backslash V_t}\frac{a_{ik}{a}_{jk}}{d_k\sqrt{d_i{d}_j}}$$

$$(S{S}^T{)}_{ii}=\frac{m}{\theta m+m+1}.$$

$$\left(T^{(t)}-\frac{\theta m+m+1}{\theta m+1}\lambda +\frac{m}{(\theta m+1)\lambda }\right)\varphi =0$$

$$R(\lambda )=\frac{\theta m+m+1}{\theta m+1}\lambda -\frac{m}{(\theta m+1)\lambda }.$$

$$J_R=\{0\}\bigcup \underset{n\to \infty }{\lim }{R}^{-n}(\{1,-1,0\}).$$

Full text

††preprint: APS/123-QED

I Introduction

Tight-binding quantum gases upon quasiperiodic or fractal-like structures with scale symmetry have been studied intensively over the past few decades Cardoso et al. (2008); Nandy et al. (2014); Yamada (2015); Jana et al. (2010); Chakrabarti (2011); van Veen et al. (2016); Moreira et al. (2008); Mauriz et al. (2001); de Oliveira et al. (2004); Coronado and Carpena (2005); Mauriz et al. (2003); Vallejos and Anteneodo (1998); de Oliveira et al. (2009); Aydin and Sisman (2014); Vidal et al. (2011); Ketterle and van Druten (1996); Burioni et al. (2001); de Oliveira et al. (2013); Bagnato and Kleppner (1991); Buonsante et al. (2002); Serva (2014); Brunelli et al. (2004); Lyra et al. (2014). In most cases, the energy spectrum of the ideal gas and corresponding density of states show self-similarity and power-law behaviors at the same time. This is responsible for a sequence of unique behaviors related to localization of states Cardoso et al. (2008); Nandy et al. (2014); Yamada (2015), quantum transport Jana et al. (2010); Chakrabarti (2011); van Veen et al. (2016), specific heat Moreira et al. (2008); Mauriz et al. (2001); de Oliveira et al. (2004); Coronado and Carpena (2005); Mauriz et al. (2003); Vallejos and Anteneodo (1998); de Oliveira et al. (2009); Aydin and Sisman (2014), Bose-Einstein condensation (BEC) Vidal et al. (2011); Ketterle and van Druten (1996); Burioni et al. (2001); de Oliveira et al. (2013); Bagnato and Kleppner (1991); Buonsante et al. (2002); Serva (2014); Brunelli et al. (2004); Lyra et al. (2014), etc. Though without introducing interaction, the simplest model yields lots of interesting phenomena due to the complex topology of fractal-like lattice structures. Different from isotropic models, the hopping of particles is non-trivial in these cases. Naturally one will ask how the heterogeneity of hopping(site-site coupling) influences the model.

The heterogeneity of hopping consists of two aspects: the network topology of lattices and the variation of coupling strength. There have been many results on how the topology of lattice structures gives birth to unusual behaviors of hopping gases. For example, the low dimensional BEC of non-interacting bosons, trapped by diamond hierarchical lattices, only takes place while the branching parameter of the trap structure is lager than 2 Lyra et al. (2014). Recently, the quantum transport on Sierpinski carpets is also found to be determined by structural parameters van Veen et al. (2016). One can guess the topological properties of lattice structures decide how curved the underlying space is for the hopping gas. Though locally similar to an isotropic Euclidean lattice, a fractal-like structure can produce totally different outcomes when serving as traps for hopping gas. To describe those anisotropic structures more quantitatively, some indicators including the fractal dimension and the spectral dimension Cassi and Regina (1993); Burioni and Cassi (1996); Rammal and Toulouse (1983); Rammal (1984); Alexander and Orbach (1982) are introduced. A deterministic relation among them is also provided for some renormalizable structures Burioni et al. (1999).

However, it is rarely reported that how the heterogeneity of the strength of site-site couplings(hopping amplitude) influences the behaviors of quantum gases. The heterogeneity of hopping amplitudes is worth studying since the site-site coupling is suggested to play a important role in other similar models. There are many cases that can not be approached by mean field approximation in real world systems. For example, the heterogeneity in the site-site coupling significantly affects the epidemic spreading  Yang and Zhou (2012); Chu et al. (2011), transportation  Wu et al. (2013); Soh et al. (2010), synchronization  Lu et al. (2006); Huang (2006), random walks Lin and Zhang (2013); Zhang et al. (2013), diffusive processes  Baronchelli and Pastor-Satorras (2010), voter models Baronchelli et al. (2011); Suchecki et al. (2005), etc., on weighted networks. We will fill this gap by a case study regarding the non-interacting Fermi and Bose gases upon a parameterized scale-invariant branching structure. We will show that the heterogeneity of coupling strength has a decisive influence on the thermodynamic behaviors even in the simplest model.

This paper is organized as the following. First we construct a scale-invariant branching structure with two parameters control the heterogeneity of our model. In tight-binding approximation we define the normalized Hamiltonian. By appropriate renormalization the full spectrum is obtained. Then, for Fermi gas, we study its Fermi energy and subsequent log-periodic oscillation of specific heat associated to special weight parameter. As for Bose gas, we investigate its phase transition phenomenon at low temperature and find the relation between the weight parameter and BEC.

II Preliminary

A weighted branching structure is constructed iteratively, see Fig. 1.

$G^{(0)}$ is a chain of length $1$ where two vertices are connected by an edge of unit weight. For $t>0$, $G^{(t)}$ is obtained from $G^{(t-1)}$ by the following transformation. For each edge of weight $w$ in $G^{(t-1)}$, $mw(m>0)$ new vertices are connected to both sides of the edge respectively with unit weight, meanwhile, the weight of the old edge is increased by $m\theta w(\theta\geqslant 0)$. The parameters $m,\theta$ are all integers. Let $G(m,\theta)=\lim\limits_{t\rightarrow\infty}G^{(t)}$. An infinite branching structure forms.

By construction, the total number of the vertices for $G^{(t)}$ is

[TABLE]

Name these vertices by $v_{1},v_{2},\cdots,v_{N_{t}}$. $a_{ij}$ denotes the weight of the edge connecting $v_{i}$ and $v_{j}$. $a_{ij}$ is [math] when $v_{i}$ and $v_{j}$ are not adjacent. Further we define the degree of $v_{i}$ as $d_{i}=\sum_{j}a_{ij}$.

To describe the topological structure of $G^{(t)}$, we introduce the adjacency matrix $(A^{(t)})_{ij}=a_{ij}$ and the degree matrix $(D^{(t)})_{ij}=\delta_{ij}d_{i}$. Let the normalized stochastic matrix Chen and Zhang (2007) for markov chains on $G^{(t)}$ be $T^{(t)}=\sqrt{D^{(t)}}A^{(t)}\sqrt{D^{(t)}}$. Obviously, $t_{ij}=\frac{a_{ij}}{\sqrt{d_{i}d_{j}}}$. For $G(m,\theta)$, define $T=\lim\limits_{t\rightarrow\infty}T^{(t)}$.

III Tight-binding model on $G(m,\theta)$

Suppose the structure we constructed denotes a trapping structure for quantum gases. The edges connecting two vertices represent the correlation of two traps. The tight-binding Hamiltonian describing the system writes Serva (2014); de Oliveira et al. (2013)

[TABLE]

Here $\hat{a}_{i}^{\dagger}$ and $\hat{a}_{i}$ are creation and annihilation operators and $a_{ij}$ denotes the hopping amplitude between two coupled traps. The second summation in Eq. (2) is taken over all neighboring vertices $i$ and $j$. Clearly, when $\theta=0$, $a_{ij}$ is constantly 1 for all existing site-site correlations. This is the most homogeneous case in our model. For non-vanishing $\theta$, the hopping amplitude is heterogeneous.

From Eq. (2), we know the spectrum of $\hat{\mathcal{H}_{0}}$ is unbound for infinite network($t\rightarrow\infty$). However, by rescaling the frequency space(multiplying the Hamiltonian by diagonal operators at both sides), we can normalize $\hat{\mathcal{H}_{0}}$ as

[TABLE]

of which the spectrum lies on $[-1,1]$.

The matrix $T$ we defined previously hence gives a full description of $\hat{\mathcal{H}}$. The allowed energy for $\hat{\mathcal{H}}$ is the eigenvalue spectrum of $-T$. The spectrum is a Julia multiset $J_{R}\subset[-1,1]$ generated by the inverse of the function

[TABLE]

from $\{-1,1\}\bigcup\{0\}$ Teplyaev (1998); Hare et al. (2012). A detailed description of $J_{R}$ is provided in the appendix. Fig. 2 shows how the eigenvalue spectra vary with the weight parameter $\theta$. When $\theta\rightarrow\infty$, the spectrum is dense in $[-1,1]$.

Let $\textrm{deg}^{(t)}(\varepsilon)$ denote the multiplicity of the eigenvalue $\varepsilon$ of $-T^{(t)}$. Take $\textrm{deg}^{(t)}(\varepsilon)=0$ if $\varepsilon$ is not an eigenvalue. The density of states on $[-1,1]$ is

[TABLE]

where $\delta(\varepsilon-\varepsilon^{\prime})$ is the Dirac delta function. Fig. 3 is a schematic representation of $\rho(\epsilon)$ related to different $\theta$. Obviously, $\rho(\epsilon)$ shows self-similar properties.

Fig. 2 and Fig. 3 together show that the spectrum related to $G(m,\theta)$ is highly degenerate and fractal-like. Besides, the spectrum is symmetric with respect to $\epsilon=0$ , which possesses the largest degeneracy. From Fig. 2, we have also found several isolated levels, among which $\epsilon=0$ is fixed for any $\theta$.

The integrated density of states is defined as

[TABLE]

Since the smallest eigenvalue of $-T$ is always $-1$ through iteration, the lowest allowed energy for $\hat{\mathcal{H}}$ is $\varepsilon_{0}=-1$. Near the lowest energy level, the power-law behavior of the integrated density of states is observed, given by

[TABLE]

$d_{s}$ is called the spectral dimension of $G(m,\theta)$, working as a generalization of the Euclidean dimension defined on homogeneous systems such as regular lattices Rammal (1984); Alexander and Orbach (1982); Cassi and Regina (1993).

Under iteration, the energy level near the bottom is stable ($R(-1)=-1$) and the integrated density of states is fixed, one obtains

[TABLE]

Plugging Eq. (7) into Eq. (8), after some algebra we arrive at

[TABLE]

Obviously, $d_{s}\geqslant 2$, where the equality only holds for $\theta=0$.

IV Fermi gas on $G(m,\theta)$

In this section we will study the behaviors of non-interacting hopping fermions trapped by the structure we have constructed. The Hamiltonian for the system is still given by Eq. (3) while the creation and annihilation operators are fermionic now. We ignore spins for evident reasons.

The Fermi-Dirac statistics give the average occupation number of each state with energy $\epsilon$:

[TABLE]

In Eq. (10), $\beta=\frac{1}{k_{B}T}$ and $\mu$ is the chemical potential. The normalization condition requires

[TABLE]

where $\gamma>0$ is the filling fraction, i.e., the ratio of particles to vertices.

Let $E_{F}$ denote the Fermi energy of this system.

To investigate the influence of the parameter $\theta$ on the Fermi energy while the structural parameter $m$ is fixed, we put $m=2$ without loss of generality. The dependence of $E_{F}$ on $\theta$ is given in Fig. 4. Except for half filling, $E_{F}$ varies with $\theta$ for fixed $\gamma$ and low $\theta$ due to the huge degeneracy of $\epsilon=0$. We can also give the range of $\gamma$ that allows $E_{F}=1$ as $\gamma\in[\gamma_{c}(m,\theta),1-\gamma_{c}(m,\theta)]$. $\gamma_{c}(m,\theta)$ is determined by checking the weight of level $\epsilon=0$:

[TABLE]

Eq. 12 yields $\gamma_{c}(m,\theta)=\frac{1}{\theta m+2m+1}$. As $\theta\rightarrow\infty$, $\gamma\rightarrow 0$ since the level $\epsilon=0$ predominates.

Next we pay attention to the thermodynamic property of the system.

The average internal energy per particle is $U=\frac{1}{\gamma}\int_{-1}^{1}\bar{n}(\epsilon)\rho(\epsilon)\epsilon d\epsilon$. Subsequently the heat capacity per particle or specific heat is

[TABLE]

where $<g(\epsilon)>_{{}_{f}}$ denotes the integral

[TABLE]

A schematic representation of the specific heat as a function of temperature is given in Fig. 5. In low-temperature approximation, Eq. (13) becomes

[TABLE]

The integral Eq. (14) is approximated by

[TABLE]

As a rough estimate we can set $1<\delta<10$. Since $\rho(\epsilon)$ is locally symmetric (as shown in Fig. 3) and $(\epsilon-E_{F})/\cosh^{2}\frac{\beta(\epsilon-E_{F})}{2}$ is strictly odd with respect to $\epsilon=E_{F}$, the second term on the right side of Eq. (15) is negligible. Thus

[TABLE]

If $\epsilon=\epsilon_{0}+E_{F}$ is an isolated point in the energy spectrum, the integral in Eq. (17) is 0. Correspondingly, the curves associated to $\theta=0$ $(E_{F}=0.1835)$, $\theta=2$ $(E_{F}=1)$ and $\theta=3$ $(E_{F}=1)$ in Fig. 5 decay rapidly as $T\rightarrow 0$. On the contrary, the curve associated to $\theta=1$ $(E_{F}=0.5)$ displays oscillating decay as $T\rightarrow 0$, because $\epsilon=\epsilon_{0}+E_{F}=-0.5$ is a limit point in the energy spectrum. This means $\rho(\epsilon)$ behaves scale-invariantly near $\epsilon_{0}+E_{F}$. The period of the oscillation, thus, depends on the periodicity of $\log\int_{-\delta}^{\delta}\rho({tk_{B}T}+\epsilon_{0}+E_{F})dt$ with respect to $\log T$. Fig. 6 gives a schematic representation of such relation. One finds $\log(C_{v,f}(T))$ has the same period as $\log\Delta I(k_{B}T)$ with respect to $\log(T)$ at low temperature. This unique effect relates the low-temperature dependence of the specific heat to the self-similarity of the spectrum straightforwardly. Moreover, this effect is only present for certain degree of heterogeneity of the coupling for fixed $\gamma$. And, in large $\theta$, this effect is only allowed for full filling. Other results regarding the log-periodic oscillation of specific heat can be found at de Oliveira et al. (2009); Aydin and Sisman (2014).

Based on above analysis, the value of the structural parameters $m$, $\theta$ and the filling fraction $\gamma$ together determine the low-temperature dependence of the specific heat. The trivial case occurs when $\gamma\in[\gamma_{c}(m,\theta),1-\gamma_{c}(m,\theta)]$, for which $E_{F}=1$ is a isolated energy. Otherwise we need to check whether $E_{F}$ represents an non-isolated level so that the specific heat will display log-periodic oscillation at low temperature.

V Bose gas on $G(m,\theta)$

As for non-interacting hopping bosons on the same structures, we focus on the low-dimensional Bose-Einstein condensation(BEC) induced by the heterogeneity of site-site coupling. BEC purely arising from topological heterogeneity has been studied before, see Lyra et al. (2014). So here we still keep $m=2$.

Ignoring spin, the average occupation number of each state is given by

[TABLE]

At low temperature, $n_{0}(T)=\frac{1}{e^{\beta(\varepsilon_{0}-\mu)}-1}$ is the number of condensed bosons and $f_{0}(T)=n_{0}(T)/N_{t}$ the condensed fraction.

The normalization condition is still Eq. (11). And the specific heat is

[TABLE]

where $<g(\epsilon)>_{{}_{f}}$ denotes the integral

[TABLE]

A schematic representation of the relation between the condensed fraction and the temperature is given by Fig. 7, where we find the value of $\theta$ determines the occurrence of BEC and the transition temperature. The corresponding relation between the specific heat and the temperature is shown in Fig. 8. It is clear that BEC occurs only for non-vanishing $\theta$. In Fig. 8 we also observe that the curves associated to $\theta>0$ have jump discontinuities, indicating the location of the transition temperature. Obviously, the transition temperature $T_{c}(\theta)$ increases monotonically and saturates when $\theta\rightarrow\infty$, as shown in the inset of Fig. 7. As one can predict, the critical relation for condensed fraction is $f_{0}(T)\propto\frac{T_{c}-T}{Tc}$ near the critical temperature.

Through Eq. (9) we find the spectral dimension related to $\theta=0$ is always $2$ no matter how large $m$ is. Since BEC is forbidden for spectral dimension not larger than 2, BEC does not take place upon $G(m,0)$ that is the homogeneous configuration. This conclusion exactly reflects the importance of the heterogeneity of the site-site coupling strength instead of the heterogeneity in topological structures. Because the latter fails to induce a phase transition.

VI Conclusion

Though it may be satisfactory to study homogeneous model for many real-world solids when disorder is introduced, the non-uniform model is necessary for plenty of other applications as we have mentioned in Sec. I. In this spirit, we investigated trapped non-interacting bosons and fermions upon a scale-invariant trap structure, with the branching parameter $m$ and the weight parameter $\theta$ controlling the heterogeneity of the site-site coupling strength. Through exact renormalization, we analytically determined the entire energy spectrum of the normalized Hamiltonian with tight-binding approximation. The spectrum was shaped by the value of $\theta$. Then we studied the unique thermodynamic behaviors of our model at low temperature with different degree of heterogeneity. To fermions, log-periodic oscillation of specific heat only occurs for certain $\theta$ provided fixed filling fraction. To bosons, the Bose-Einstein condensation is forbidden for the homogeneous setup($\theta=0$) regardless of the value of $m$.

The scale-invariant branching structure in this study is exactly renormalizable hence allows efficient analytical study. At the same time other irregular structures can not enjoy such convenience. Nevertheless, the deterministic renormalizable structures give us precious insight into the heterogeneity of network-like structures. Be more specific, our model serves as a good prototype for plenty of scale-invariant models all displaying fractal-like spectra Cardoso et al. (2008); Nandy et al. (2014). In fact one can distort a solved structure while keeping the emergent behaviors unchanged. For example, the branching parameter $m$ in our model can actually vary with different branches while keeping BEC forbidden for $\theta=0$.

As a final comment, if we regard the topological aspect of structural heterogeneity as the shape of a circuit, the heterogeneity in site-site couplings alters the distribution of resistance. The latter puts the emphasis on some parts of a lattice and weakens the rest. In this circumstance, more features are added to a model that is previously uniform. This is the reason why we can find lots of interesting phenomena in the simplest model we consider here.

Acknowledgements.

This work was supported by the National Natural Science Foundation of China under Grant No. 11275049. P. Xie and B. Yang were also supported by Fudan’s Undergraduate Research Opportunities Program, FDUROP. Also, P. Xie appreciated the fruitful discussion with Jin-Yue Su.

Appendix

For clarity, name the vertex set of $G^{(t)}$ by $V_{t}$. The set of the new nodes generated at iteration $t+1$ is denoted by $V_{t+1}\backslash V_{t}$. By proper arranging $T^{(t+1)}$ can be expressed as

[TABLE]

In Eq. (21) $S$ reflects the coupling between the new vertices generated at iteration $t+1$ and the rest. Correspondingly partition the eigenvector of $T^{(t+1)}$ as $\psi=(\phi,\varphi)$ with respect to eigenvalue $\lambda$. Hence

[TABLE]

which yields

[TABLE]

provided $\lambda\neq 0$.

By the construction procedure,

[TABLE]

is nonzero if and only if $i=j$ since each $v_{k}$ is adjacent to one vertex only. By simple algebra one arrives at

[TABLE]

Due to Eq. (26), Eq. (24) is simplified as

[TABLE]

which entails

[TABLE]

$R(\lambda)$ is the eigenvalue of $T^{(t)}$ associated to eigenvector $\phi$. Further one will find the multiplicity of $R(\lambda)$ to $T^{(t)}$ is the same as the multiplicity of $\lambda$ to $T^{(t+1)}$.

Without proof we claim that the preimage of $\lambda$ under function $R$ is a subset of the eigenvalue spectrum of $T^{(t+1)}$. Any eigenvalue in the preimage of $\lambda$ carries on the multiplicity of $\lambda$. Notice that by tracking the preimage of the spectrum $T^{(t)}$ under function $R$ one can obtain only $2N_{t}$ eigenvalues of $T^{(t+1)}$. Fortunately, the left $N_{t+1}-2N_{t}$ eigenvalues are uniformly 0 that is the only singularity of $R$. More detailed discussion about the iterative derivation of spectrum can be found at Xie et al. (2016, 2015).

So far we are able to obtain the full spectrum of $T^{(t)}$ by tracking the flow generated by $R$. As $t\rightarrow\infty$, the spectrum grows into a Julia set $J_{R}$ given as

[TABLE]

As for $-T$, the result is the same since the spectrum of $T$ is symmetric about [math].