Quantum quenches and thermalization on scale-free graphs

Abstract

We show that after a quantum quench of the parameter controlling the number of particles in a Fermi–Hubbard model on scale-free graphs, the distribution of energy modes follows a power-law dependent on the quenched parameter and the connectivity of the graph. This paper contributes to the literature of quantum quenches on lattices, in which, for many integrable lattice models the distribution of modes after a quench thermalizes to a generalized Gibbs ensemble; this paper provides another example of distribution which can arise after relaxation. We argue that the main role is played by the symmetry of the underlying lattice which, in the case we study, is scale free, and to the distortion in the density of modes.

Introduction

There has been recent interest in the effective thermal dynamics following a quantum quench in spin chains [1]. The dynamics out of equilibrium of quantum systems [2] has received a great amount of attention [3, 4]. It became clear that after a quantum quench, many observables at equilibrium after a quantum quench are distributed according to a Generalized Gibbs Ensemble (GGE) [1, 4–8] or a Gibbs ensemble. In addition to these theoretical understandings of the thermalization in quantum systems, these results are supported by recent experiments in trapped cold atomic gases [9–12]. Many of these experiments focused on the role played by dimensionality and conservation laws, which in turn initiated a vigorous effort in understanding the role of the integrability of the system under scrutiny at late times. In fact, the easiest way of driving a quantum system out of equilibrium is indeed a quantum quench (see for instance [13] for a comprehensive summary), i.e. a sudden change in the parameters of a Hamiltonian, and its subsequent relaxation at long times. It has been argued that for integrable models many observable are distributed according to a generalized Gibbs ensemble, meanwhile an effective Gibbs distribution arises in generic systems. Recent studies suggest that the behavior is indeed more complicated, showing a dependence on the initial conditions (state) [14].

In a previous work, we studied quantum quenches in a Fermi–Hubbard model which does not conserve the number of particles [15]. We studied the energy of the excitations, which are invariant under time evolution after the quench, and found that these are distributed according to a GGE. There, the temperature is associated with the gap in the spectrum, which is due to the coupling of non-conserving number of particles term. A similar phenomenon happens in quantum liquids [16].

In this paper we explore a similar approach on a different type of underlying interaction network. Several classical statistical models have been studied on complex networks [17]. Complex networks have become an area of tremendous recent interest since the discoveries of the small world and scale-free properties in many realistic networks. A small-world network is characterized by short network distance a high clustering coefficient. Several reviews of the subjects are now available [18, 19]. Important applications of these techniques are spreading of diseases [20] and syncronization [21] on complex networks. Watts and Strogatz demonstrated that the two small-world characteristics can be obtained from a regular network by rewiring or adding a few long-range links shortcuts, which connect otherwise distant nodes [22]. A regular network intrinsically already has a high clustering coefficient, and has a diameter which is logarithmic in the number of nodes. However, few shortcuts can reduce the distance exponentially unaffecting the clustering coefficient, but dramatically reducing the average diameter of the graph. The scale-free property is characterized by an algebraic degree distribution: where the degree variable k measures the number of links of node in the network, $P(d\gg 1)\approx d^{\mathit{\gamma }}$, and where $\mathit{\gamma }$ is the algebraic scaling exponent. Barabási and Albert discovered the scale-free property and also proposed growth and preferential attachment as the two basic mechanisms responsible for the scale-free property. Here, growth requires that the numbers of nodes and links increase with time and preferential attachment means that when a new node is added to the network, the probability that it connects to an existing node is proportional to the number of links that this node has already had. In this paper, we consider the properties of quantum quenches on Barabási–Albert type of graphs. It is important to understand what kind of distributions can arise in a quantum quench. Although per se the study of quantum phenomena on complex networks might not be physically relevant, it is important from the theoretical point of view. While GGE is a quite common example of distribution arising, here we show that another type of distribution arises. This paper is organized as follows. In second section we introduce Quantum Quenches, Scale-Free graphs and "The model". In next section we describe the "Results", meanwhile "Conclusions" follow in last section.

The model

Quantum quenches and GGE

In this section we provide a more detailed introduction to thermalization after a quantum quench, and introduce the generalized Gibbs ensemble which we will later discuss of. Ergodicity in a classical statistical physics setting accounts for the independence of the asymptotic state distribution from the initial condition: that is, at large times thermal equilibrium is approached. Meanwhile in classical mechanics thermalization occurs thanks to ergodicity, in quantum mechanics of isolated systems the unitary dynamics is an obstruction to obtaining an analogous result: if the system is initially prepared in a pure state it will remain in a pure state, as the evolution is period or quasiperiodic, which means that after a sufficient long time it will return to the initial state. However, one can focus on certain expectation values only or trace out some part of the system, i.e. subsystems of the whole are not isolated, and thus the reduced density matrix is not pure anymore. In general, the type of question one asks in thermalization of quantum system is: how close is expectation value in a subsystem, to those of the same degrees of freedom averaged with a Gibbs ensemble? In general, one does not get exactly a Gibbs ensemble, but a generalized Gibbs ensemble:

$$\begin{array}{c}\hfill {\mathit{\rho }}_{\text{GGE}}\approx {\mathrm{e}}^{-{\sum }_j{\mathit{\lambda }}_j{I}_j}\end{array}$$

(1)

in which one introduced a Lagrange multiplier ${\mathit{\lambda }}_j$ and $I_j$ is the integral of motion, which generalizes the standard Gibbs distribution. The GGE-conjecture states that the stationary expectation value of any local observable is equal to the ensemble expectation values or equivalently that of the reduced density matrix of local observables. The subtlety is that one does not have a recipe for choosing the integrals of motion in eq. (1). The general approach to study thermalization is through the device of quantum quenches: you change suddenly a parameter and observe how the system thermalizes. The GGE hypothesis has been tested extensively for many systems, and has turned out to be valid for many quantum quench problems studied recently. In general, this is true for non-interacting models, or models that can be casted into a non-interactive one (quadratic). In general, the mathematical device used to prove that the GGE is valid is the Wick theorem, and that the initial state overlap with the post-quench quasi-particle modes is gaussian. An underlying assumption is thus that different modes are orthogonal, which is crucial to prove Wick’s theorem. In fact, if the underlying lattice is translational invariant, one obtains that two different modes, in general of the form ${\mathrm{e}}^{ikx}$, are orthogonal, e.g. $⟨{\mathit{\varphi }}_{k_1}(x),{\mathit{\varphi }}_{k_2}(x)⟩=\int {\mathrm{e}}^{i(k_1-k_2)x}\mathrm{d}x=\mathit{\delta }(k_1-k_2).$ In this work we consider the case in which the expansion does not have such a clear interpretation in terms of momentum, which occurs in the case in which translational invariance is not explicit in the underlying lattice. However, we consider the case in which the graph has an overall symmetry, i.e. the distribution of the degree is well approximated by a power-law for large values, and thus is scale invariant in its tail, which will be introduced in the following section. In addition to this, the density of energy modes is, differently from the case of translational invariant graphs, distributed according to a power-law.

Scale-free graphs

We now recall the growth algorithm used in the preferential attachment model introduced by Barabási and Albert [23]. The growth algorithm (preferential attachment) is parameterized by a single parameter, $M$. The starting graph is a single node, with no edges. Then, at each step, a new node is added, with $M$ edges. The edges are attached at random to the previous existing nodes, with a probability proportional to the degree of the node. If $d_i$ is the degree of the vertex $i$, at each edge is attached to a node i with probability $p_i=\frac{d_i}{{\sum }_i{d}_i}$. As it is well known, these graphs are scale free, i.e. for $N\gg 1$, the degree distribution is a power-law, $P(d\gg 1)\approx d^{\mathit{\alpha }}$, with $\mathit{\alpha }$ exponent of the power-law (Fig. 1). Another property, here important to mention, is that scale-free graphs are ultra-small: the average distance between two nodes goes as $\sim log(log(N))$, where $N$ is the number of nodes in the network. Notably, Bose–Einstein condensation appears in growing networks if preferential attachment growth is generalized with fitness [24].

Fig. 1

An example of Barabási–Albert graph

Hamiltonian and quantum quench

Here we want to recollect the formalism introduced in [15]. The Hamiltonian we will consider in the present paper is the following Fermi–Hubbard model:

$$\begin{array}{cc}\hfill H({\mathrm{\Gamma }}_M,\mathit{\lambda })=& -J\sum _{i,j=1}^{N_v}{A}_{ij}^{({\mathrm{\Gamma }}_M)}{a}_i^{†}{a}_j\hfill \\ \hfill & +\frac{\mathit{\lambda }}{2}\sum _{i,j=1}^{N_v}{B}_{ij}^{({\mathrm{\Gamma }}_M)}\left(a_i^{†}{a}_j^{†}+\mathrm{h}.c.\right),\hfill \end{array}$$

(2)

where $a_i$$(a_i^{†})$ is the annihilation (creation) fermionic operator that annihilates (create) a particle in the vertex $i$ of the background graph ${\mathrm{\Gamma }}_M$. The matrices $A_{ij}^{({\mathrm{\Gamma }}_M)}$ and $B_{ij}^{({\mathrm{\Gamma }}_M)}$ are, respectively, the adjacency matrix of ${\mathrm{\Gamma }}_M$ and its antisymmetrized form. In the present paper, the adjacency matrix will be the one of a scale-free graph built using the Barabási–Albert growth algorithm. The sum runs over all the $N$ nodes of the graph ${\mathrm{\Gamma }}_M$, where $M$ is the connectivity parameter introduced previously. The coupling $J$ is the tunneling of the particles between two connected sites and $\mathit{\lambda }$ controls the strength of the Hamiltonian terms that do not conserve the number of particles. The physical properties are independent from time scaling if we perform a sudden quantum quench, thus we can measure excitations in units of $J$.¹

In particular, we introduced a notion of particle [with an associated discrete labeling $k$, $\mathit{ϵ}(k)$] given in terms of ladder operators ${\mathit{\eta }}_k$. Once the notion of particle that the detector measures is established, we can determine the energy distribution (number of particles with momentum $k$) of the ground state of the system

$$\begin{array}{c}\hfill n(k)=⟨\text{GS}|{\mathit{\eta }}_k^{†}{\mathit{\eta }}_k|\text{GS}⟩.\end{array}$$

(3)

For graphs with discrete translational invariance, this is associated with the Fourier transform vectors $e^{ikx}$, but for graphs without particular symmetry this identification is lost.

The notion of particle ${\mathit{\eta }}_k$, together with its dispersion relation $\mathit{ϵ}(k)$ will be defined in terms of a test Hamiltonian

$$\begin{array}{c}\hfill H_{\mathrm{test}}=\sum _k\mathit{ϵ}(k){\mathit{\eta }}_k^{†}{\mathit{\eta }}_k.\end{array}$$

(4)

The momentum distribution (3) that the detector measures is given by the overlap between the ground state of the system and the eigen states of the test Hamiltonian $H_{\mathrm{test}}$.

The quantum quench we are going to perform is given by

$$\begin{array}{c}\hfill H({\mathrm{\Gamma }}_M,\mathit{\lambda })\to H({\mathrm{\Gamma }}_M,0),\end{array}$$

from the ground state of $H({\mathrm{\Gamma }}_M,\mathit{\lambda })$, $|\text{GS}⟩$, after the quench. The system $H({\mathrm{\Gamma }}_M,0)$ will see the ground state $|\text{GS}⟩$ as an excited state, and thus it makes sense to calculate the spectrum density $n(\mathit{ϵ}(k))$. The hopping Hamiltonian can be written as

$$\begin{array}{c}\hfill H_{\mathrm{test}}=\sum _{i,j=1}^{N_v}{A}_{ij}^{({\mathrm{\Gamma }}_M)}{a}_i^{†}{a}_j=\sum _{k=1}^{N_v}\mathit{ϵ}(k){\mathit{\eta }}_k^{†}{\mathit{\eta }}_k.\end{array}$$

(5)

The eigenmodes of $H_{\mathrm{test}}$, labeled by an integer $k$, and with energy $\mathit{ϵ}(k)$, define our notion of particle. These are created and annihilated by the operators ${\mathit{\eta }}_k^{†}$ and ${\mathit{\eta }}_k$, and are the excitations that the detector measures, and that we will calculate. Therefore, we need to compute

$$\begin{array}{c}\hfill n(k)=⟨\text{GS}|{\mathit{\eta }}_k^{†}{\mathit{\eta }}_k|\text{GS}⟩,\end{array}$$

(6)

and calculate the distribution. As we will see, $n(k)\approx \mathit{ϵ}{(k)}^{\mathit{\gamma }}$, where $\mathit{\gamma }$ is the exponent we will study. The two Hamiltonians have the same number of nodes, thus their Hilbert states overlap (coincident).² If the system was at equilibrium with an external bath at temperature $T$, we would have $n(k)={\mathrm{e}}^{-\frac{\mathit{ϵ}(k)}{T}}$. In the case of an $n$-dimensional tori, we found in [15] for the same model we study here, that $n(k)={\mathrm{e}}^{-\frac{2\mathit{ϵ}(k)}{\mathit{\lambda }}}$, where $\mathit{\lambda }/2$ plays the role of the temperature.

The spectrum of the adjacency matrix of a scale-free network was studied in [25]. It is known that eigenvalues are distributed according to a power-law, while for the eigenvectors an analytical form is still lacking. However, the components of the eigenvectors are strongly localized at the hubs.

As a matter of fact, in our numerical calculation we will neglect the role of the shape of the Fermi surface, approximating it with a sphere. This implies that a small error is made on the temperature, that can become more and more relevant for $\mathit{\lambda }\approx 1$. Thus, our result is valid in the limit $|\mathit{\lambda }|\ll 1$.

Results

Here we show that the distribution of the modes in the ground state follows a power-law. Since analytical techniques are lacking, we evaluated Eq. (6) numerically. All the numerical results are obtained with a number of nodes $N=700$. We can tackle such a big quantum system, due to the Bogoliubov transformation which diagonalizes this model; the quantum dynamics of this system is in fact restricted to a Hilbert space of size $2N$, instead of the $2^N$ of an ordinary spin system. In Fig. 2 we plot $log(n(k))$ against $log(\mathit{ϵ}(k))$. The functional dependence appears, numerically, to be of the form:

$$\begin{array}{c}\hfill n(k)\sim f(M,\mathit{\lambda })\mathit{ϵ}{(k)}^{\mathit{\gamma }}\end{array}$$

(9)

where $\mathit{\gamma }$ is weakly dependent on $\mathit{\lambda }$, in the range $\mathit{\lambda }\in [0,0.2]$ and takes value $\mathit{\gamma }\in [-0.85,-0.95]$. Tables 1 and 2 show the values obtained numerically by fixing $M$ and $\mathit{\lambda }$ and fitting the power-law, together with the error, of the parameters $f$ and $\mathit{\gamma }$. In Figs. 4 and 5, we plot the functional dependence of the constant in front of the Zipf’s law, by fixing the values of $M$ and $\mathit{\lambda }$. These appear to be both convex functions of the parameters. The power-law is rather stable over several order of magnitudes, $log(\mathit{ϵ}(k))\in [-2,2]$, although the density of points is higher in $log(\mathit{ϵ}(k))\gtrsim 0$. For $\mathit{\lambda }\in [0.2,1]$, the distribution changes shape in $log(\mathit{ϵ}(k))<0$. Thus, restoring the units, the power-law for density of excitations is valid only in the limit $\frac{\mathit{\lambda }}{J}\ll 1$. It is clear from this analysis that the power-law exponent, in its domain of validity, is independent from the connectivity parameter $M$ (Fig. 3).

One can give a rough explanation of the results above, by considering the exact expression of $n(\mathit{ϵ})=\frac{{\mathit{\lambda }}^2(1-{\mathit{ϵ}}^2)}{{\mathit{\lambda }}^2(1-{\mathit{ϵ}}^2)+{\mathit{ϵ}}^2+\mathit{ϵ}\sqrt{{\mathit{\lambda }}^2(1-{\mathit{ϵ}}^2)+{\mathit{ϵ}}^2}}$ [15]. For $\mathit{\lambda }\ll 1$, one has that $n(\mathit{ϵ})\approx \frac{\mathit{\lambda }}{2}(\frac{1}{{\mathit{ϵ}}^2}-1)$. However, if the density of eigenvalue $\mathit{\rho }(\mathit{ϵ})$ is distributed according to a power-law, one has to smear the effective distribution considering the deformed density. In fact, if we consider $⟨n(\mathit{ϵ})⟩=\int n(\mathit{ϵ})\mathit{\rho }(\mathit{ϵ})d\mathit{ϵ}\approx \frac{1}{{\mathit{ϵ}}^{1+\mathit{\alpha }}}$, depending on the right functional, asymptotic expression of $\mathit{\rho }(\mathit{ϵ})$, and considering that usually the function $\mathit{\rho }(\mathit{ϵ})$ is peaked at $log(\mathit{ϵ})\approx 0$, which is what we observe numerically. This should explain why we observe a power-law for small values of $\mathit{\lambda }$.

Fig. 2

Log–Log plot of the distribution $n(\mathit{ϵ}(k))$, evaluated for $N=700$, $M=10$, $\mathit{\gamma }=0.1$. $log(n)$ on the y-axis, $log(\mathit{ϵ})$ on the x-axis

Fig. 3

Log–Log plot of the distribution $n(k)$ for $M=10$., $\mathit{\lambda }=0.8$ and $N=700$

Conclusions

In this paper we discussed quantum quenches of a Fermi–Hubbard model on scale-free graphs, motivated by the search to alternative distributions from those of the generalized Gibbs ensemble hypothesis, relevant in particular in the case of integrable lattice models. A previous analysis of the quench protocol discussed in the present paper was done in [15], and solved analytically for the case of $n$-dimensional torii, showing that the generalized Gibbs ensemble emerges (GGE). The quenched parameter, $\mathit{\lambda }$, controls the conservation of the particle number and, in particular, introduces a gap in the spectrum. For non-integrable systems, it is a known fact that the expectation value of several observables after a quantum quench, is similar to those calculated on a thermal state, i.e. a GGE. The observables we considered were the density of eigenmodes calculated over the ground state of the unquenched Hamiltonian. In this work, we have shown that the spectrum of the excitations for scale-free graphs can be well approximated by a power-law for the case in which the quenched parameter is small. In particular, we have analyzed the functional dependence of the two parameters of the power-law on the topological properties of the scale-free graph, the connectivity, and the quenched parameter. We found that while for small values of $\mathit{\lambda }$ this distribution is scale free, meanwhile for higher values of the quenched parameter ($\mathit{\lambda }>0.2$), the distribution changes shape and the power-law behavior is lost.

The outcome of our analysis shows that the underlying symmetry of the graph does indeed contribute to the shape of the distribution of modes. While GGE naturally emerges in the case of translational invariant lattices, we have given evidences for a counter example, in which a power-law distribution for the modes spectral density arises. Due to the lack of analytical understanding of the eigenmodes for scale-free graphs, we have reached these conclusions using a numerical approach. We have focused on a Fermi–Hubbard model for a quite simple, technical reason. The scale-free properties of graphs generated using a preferential attachment become evident for large graphs, e.g. for the case $N\to \infty$. We thus had to focus on a model which can be diagonalized numerically for large graphs. For the Fermi–Hubbard model considered in the present paper, we could take advantage of a Bogoliubov transformation which allows to tackle the diagonalization process on a matrix of size proportional to the number of nodes, and not exponentially increasing, as standard in quantum mechanics.

We have tried to show in a numerically treatable model the relevance of the symmetry of the underlying lattice for the GGE hypothesis. Although striking evidences have been put forward, both analytically and numerically, of the universality of the GGE for integrable models [4], these results rely on the underlying symmetry of the graph, in particular on the properties of the eigenvectors of the adjacency matrix of the lattice. These properties are indeed lost when translational invariance is not present, which is the case for the scale-free graphs. However, we have shown that the distribution does indeed retain some properties of the underlying lattice, which for the case of scale-free graphs is the scale-freeness of the spectrum distribution. We have provided a rough, analytical explanation of why, in the case in which the density of modes is not uniform, one has to consider a smeared version of the spectral density.

Appendix A

See Tables 1, 2; Figs. 4, 5.

Table 1 Table of $\mathit{\gamma }$ and $f$ as function of $M$, for fixed $\mathit{\lambda }=0.05$ and $N=700$

Table 2 Table of $\mathit{\gamma }$ and $f$ as function of $\mathit{\lambda }$, for fixed $M=10$ and $N=700$

Fig. 4

Plot of the $f$ as a function of $\mathit{\gamma }$ for $M=10$ and $N=700$

Fig. 5

Plot of $f$ as a function of $M$, for $\mathit{\lambda }=0.05$ and $N=700$