Entanglement Spectra and Entanglement Thermodynamics of Hofstadter Bilayers

We study Hofstadter bilayers, i.e. coupled hopping models on two-dimensional square lattices in a perpendicular magnetic field. Upon tracing out one of the layers, we find an explicit expression for the resulting entanglement spectrum in terms of the energy eigenvalues of the underlying monolayer system. For strongly coupled layers the entanglement Hamiltonian is proportional to the energetic Hamiltonian of the monolayer system. The proportionality factor, however, cannot be interpreted as the inverse thermodynamic temperature, but represents a phenomenological temperature scale. We derive an explicit relation between both temperature scales which is in close analogy to a standard result of classic thermodynamics. In the limit of vanishing temperature, thermodynamic quantities such as entropy and inner energy approach their ground-state values, but show a fractal structure as a function of magnetic flux.


Introduction
The study of quantum entanglement has by now developed to a mature subfield of many body physics [1,2,3]. Among the recent developments, the concept of the entanglement spectrum [4] has been applied to a plethora of different systems. These comprise quantum Hall monolayers at fractional filling [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] where various spatial geometries and ways of separating the system into subsystems have been investigated. In the case of bilayer quantum Hall systems, a very natural way of defining a subsystem is provided by the double layer structure. Specifically, a numerical investigation of quantum Hall bilayers at total filling factor of unity revealed a striking similarity between the entanglement spectrum of the composite systems and the energy spectrum of a monolayer [20]. A analogous observation was reported slightly earlier in a numerical study of two-leg spin-1/2 ladders [21] which was subsequently understood in terms of perturbation theory in the limit of strong rung coupling [22,23], a result which is remarkably also valid for arbitrary spin length [24]. Similar findings were obtained in the simple case of coupled chains of free fermions where explicit expressions for the entanglement spectrum at arbitrary coupling can be derived [22,25]. The starting point of the present work is to extend these results to Hofstadter bilayers, i.e. coupled hopping models on twodimensional square lattices in a perpendicular magnetic field [26]. The energy spectrum of such systems generates, as a function of the magnetic flux per unit cell, highly selfsimilar and visually appealing structures known as Hofstadter butterflies. As we shall see below, these features immediately translate to the entanglement spectrum via an explicit formula for the entanglement levels in terms of the energy eigenvalues of the monolayer system. Another focus of the present work are thermodynamic properties of the reduced density matrix and its entanglement Hamiltonian. In particular, we derive a thermodynamic entanglement temperature and inner energy starting from an effective coupling parameter of the bilayer system which can be viewed as a phenomenological entanglement scale.
This paper is organized as follows. In section 2 we review and extend previous results on free fermions in two-party lattice systems in the absence of a magnetic field. In the following section 3 we analyze the entanglement spectrum of Hofstadter bilayers and derive an explicit formula for the entanglement levels as a function of the underlying (highly fractal) monolayer system. Section 4 is devoted to a thorough analysis of the thermodynamic properties of the reduced density matrix and its entanglement Hamiltonian. We close with a summary and an outlook in section 5.

Free Fermions without Magnetic Field
Let us first briefly review and generalize results of Ref. [22] on free fermions on coupled lattices. We consider the Hamiltonian where a + k ,a k (b + k ,b k ) generate and annihilate fermions with wave vector k and energy ε A ( k) (ε B ( k)) on a d-dimensional lattice A (B). The two subsystems are coupled by a hopping term proportional to t ⊥ leading, at given wave vector, to a 2 × 2 eigenvalue problem whose elementary solution is and where Let us now focus on the case where the above dispersion bands do not overlap, i.e.
for all wave vectors k 1 , k 2 . This is in particular the case if Indeed previous work [22,25] has concentrated on the situation ε A ( k) = −ε B ( k) where the inequality (8) is obviously valid. For a half-filled system fulfilling (7), only the singleparticle states generated by β + k are occupied in the ground state, and upon tracing out, e.g., subsystem B one obtains, following the very general arguments given in [50,51], a reduced density matrix of the form with Z = tr(exp(−H ent )) and an entanglement Hamiltonian The entanglement levels ξ( k) are given by The exponential form (9) of the reduced density matrix can be derived by demanding that this expression should generate the same one-particle correlations in the monolayer subsystem as the underlying pure bilayer ground-state density operator. Using elementary properties of fermionic systems leads then to the explicit result (11) for the entanglement levels [50]. We note that the entanglement spectrum as discussed is generated by the single-particle operator (10) as appropriate for free fermions. This is different from the situation in interacting systems where the eigenstates of the reduced density matrix are in general nontrivial many-body states. The above considerations extend the results of Peschel and Chung [22] to arbitrary dimension d ≥ 1 and independent dispersions ε A ( k) = −ε B ( k) in the two subsystems. We note that the entanglement spectrum (13) depends, differently from the energy spectrum (3), only on the difference ε A ( k)−ε B ( k), but not on both quantities separately. Moreover, both the energy and the entanglement spectrum are invariant under a change of sign of t ⊥ .
If the two energy bands overlap. i.e. condition (7) is violated, the situation cannot be analyzed in such general terms. Here the ξ( k) will in general not have support in certain parts of the first Brillouin zone, while in other parts two branches of entanglement levels will occur.

Hofstadter Bilayers: Entanglement Spectra
Let us now concentrate on two coupled two-dimensional square lattices with the Hamiltonian where ,n a m,n + a + m,n+1 a m,n + h.c. , and Here a + m,n , b + m,n generate particles at sites m,n in layer A and B, respectively., i.e. compared with equation (1) we have where a is the lattice constant in both layers. In order to implement a perpendicular magnetic field we work, as common, in Landau gauge with the discretized vector potential A = (0, Bma, 0) and concentrate on rational values of the magnetic flux per unit cell Φ = Ba 2 /(h/e) = p/q (in units of h/e) where p and q are integers without any common divisor except unity. Thus, the system is periodic in the x-direction with periodicity qa while the periodicity in the y-direction is still the lattice constant a. The standard Peierls substitution [52] leads the following ansatz for the state vector, with α m ( k) = α m+q ( k) and β m ( k) = β m+q ( k). Here |0 is the fermionic vacuum, and N x , N y are the total numbers of unit cells in x-and y-direction, respectively, assuming periodic boundary conditions commensurate with q. The stationary Schrödinger equation leads to a Harper equation being equivalent to the 2q × 2q eigenvalue problem for the 2q-component spinor and Here the matrix with r m ( k) = 2 cos(2πΦm + k y a) and z( k) = exp(ik x a) corresponds to the classic Hofstadter monolayer problem [26]. The diagonalization of the latter quantity in general requires numerics, and explicit analytical results are possible only for very special values of the magnetic flux Φ = p/q.
Since the off-diagonal elements in (23) are proportional to the q × q unit matrix, they remain unchanged upon a simultaneous diagonalization of both diagonal blocks. Thus, all four blocks are rendered diagonal, and the diagonalization of (23) is reduced to 2 × 2 problems of the form with the obvious eigenvalues and corresponding eigenvectors   Fig. 1 shows the numerically computed classic Hofstadter spectrum of a single square lattice as a function of the magnetic flux Φ per unit cell, along with the energy spectrum of a bilayer system for the particularly simple case t A = t B = t ⊥ . Here the Hamiltonian is invariant under exchange of layers, and the eigensystem consists of states being either symmetric or antisymmetric under this operation. Thus, as seen in the left panel of Fig. 1, the bilayer energy spectrum comprises two monolayer Hofstadter butterflies shifted by ±t ⊥ .
We now concentrate again on the case where both groups of dispersion branches do not overlap, for all k 1 , k 2 and m 1 , m 2 ∈ {0, . . . , q − 1}. Analogously to (8) the stronger condition (i.e. t A and t B need to differ in sign) implies ε + m 1 ( k 1 ) ≥ 0 ≥ ε − m 2 ( k 2 ) and therefore the inequality (30). Thus, under condition (30) only the single-particle states with energies ε − m ( k) are occupied in the ground state of a half-filled system. Tracing out again layer B leads to entanglement levels of the form These entanglement levels enter the entanglement Hamiltonian at given flux Φ = p/q. Here the operators a + km generate Hofstadter monolayer states with eigenvalueε m ( k), i.e. the monolayer Hamiltonian including the magnetic field can be formulated as We note that the matrix (24) has the obvious property with π = (π/a, π/a) which implies where we have the eigenvaluesε m ( k) of h( k) assumed to be given in ascending order. These relations enable to use, with some caveat and specifications, the expression (34) for the entanglement levels even in situations where the condition (30) is not fulfilled: The left panel of Fig. 2 shows the dispersionsε m ( k) for a magnetic flux of Φ = 1/7 and typical parameters. As seen in the figure, the eigenvalues form rather flat bands which do not overlap. The right panel displays the corresponding bilayer energy ε ± m ( k) for parameters violating the inequalities (30), (31). In the ground state of a half-filled system the lowest q = 7 bands are completely occupied such that the highest occupied (valence) band is ε + 1 ( k), which is also the highest band among the branches ε + m ( k) having negative energy. Note that ε + 1 ( k) does also not overlap with any other band ε ± m ( k); otherwise such two bands would only partially be filled each. Let us therefore focus on this "insulator" scenario where the valence band of a half filled system does not overlap with other bands implying that in the ground state all bands are either fully occupied or empty. Let us now assume that some band ε + m ( k) has negative energy and  (30), (31). The situation for other fluxes Φ = p/q is qualitatively similar, but with increasing q more difficult to display due to the larger number of bands.
is therefore fully filled. Thus, according to Eq. (39) ε − q−1−m ( k) has positive energy and is therefore empty. The entanglement branch arising from ε + m ( k) is where we have used Eqs. (29), (39). Thus, up to a rigid shift in the wave vector argument, the entanglement levels arising from the occupied band ε + m ( k) reproduce the missing branch corresponding to ε − q−1−m ( k). In this sense (i.e. with the wave vector dependence being suppressed) the expression (34) for the entanglement spectrum can also be used even if the inequality (30) does not hold, but all bands are, in the ground state, either empty or completely occupied. On the other hand, starting directly from the original 2q × 2q eigenvalue problem (21), the entanglement levels can be expressed as where and the α  In summary, the expressions (27) and (34) provide, in full analogy with Eqs. (3), (13), explicit relations between the energy and the entanglement spectrum, respectively, of the composite system and the energy spectrum of the monolayer. Moreover, the entanglement spectrum ξ m ( k) depends, for a given monolayer spectrum ε m ( k) only on the single quantity which will in the following be referred to as the (effective) coupling parameter. As already discussed above, in the case λ = 0 the layers are maximally entangled with each other, and the all entanglement levels are zero. In Fig. 3 we have plotted the energy spectrum of the bilayer system according to Eq. (27) for t A = t ⊥ and various values of t B . As seen there, if t A and t B differ in sign (here: t B < 0), an energy gap which is independent of the magnetic flux Φ occurs. Fig. 4 displays the corresponding entanglement spectra which form Hofstadter butterflies "deformed" by the inverse hyperbolic sine function occuring in Eq. (34).

Hofstadter Bilayers: Entanglement Thermodynamics
In the spirit of Ref. [20] we now investigate the thermodynamics based on the entanglement Hamiltonian (35) and the pertaining reduced density operator ρ red . From Eq. (34) one finds implying (cf. Eq. (36))) for small coupling parameter λ. This suggests to interpret λ as an inverse "entanglement temperature" and approximate the reduced density matrix as ρ red ≈ ρ 0 with with Z 0 = tr(exp(−λH A /t A )) and H A /t A describing energy of the subsystem. The above findings are completely analogous to the ones in Ref. [20] on quantum Hall bilayers where the above approximate relation (46), valid for strongly coupled layers, was used to numerically analyze a quantum phase transition in the total double layer system. However, having a closed analytical result for the entanglement levels at hand, we shall follow a different route here and employ the full expressions (34),(35) (not necessarily approximated by their first order in λ) to evaluate thermodynamics in a more systematic manner. As we shall explore in more detail below, the coupling parameter λ can be used as a phenomenological temperature scale, but it is not identical to the inverse thermodynamic temperature given by the derivative of the entropy with respect to the (appropriately defined) inner energy. A very simple issue is the average particle number N = N withN = m k a + km a km and · = tr(ρ red ·). Here one always has, independently of coupling parameter and magnetic flux, N (λ.Φ) = N x N y /2 which is clear from the fact that the total bilayer system is half-filled and unbiased with respect to its subsystems. Formally this result can be established, in the thermodynamic limit N x , N y 1, as follows, where we have used Eqs. (33) and (38). This constancy of the averaged particle number can of course also be seen as a consequence of the fact that a reduced density matrix of the form (9) can be viewed as a grand-canonical statistical operator with constant chemical potential µ = 0, which lies exactly in the middle of the symmetric entanglement spectrum. Let us next turn to the other extensive quantities entropy and energy which we define, again in the thermodynamic limit, by The bar at the above quantityĒ is meant to indicate that the definition of energy will be subject to some refinement further below. Moreover, as a consequence of Eq. (40), this quantity is non-positive, for all λ and Φ. Besides, since the entropy S is obviously proportional to the area of the system coinciding with the boundary to the other subsystem, the so-called area law is fulfilled [3].
For small λ 1 one obtains the expansions where we have used the fact that vanishes for odd n, and for low even values one has ω 0 (Φ) ≡ 1, ω 2 (Φ) ≡ 4, and ω 4 (Φ) = 4 (7 + 2 cos(2πΦ)) , as it is readily derived from the explicit form (24) of the matrix h( k). We note that the zeroth-order value s(0, Φ) = (ln 2)/2 of the entropy per unit cell is due to the fact that, at strong coupling, each particle ends up in either layer with equal probability, and the total bilayer system is half-filled. In Fig. 5 we have plotted the entropy s(λ, Φ) per unit cell as a function of the magnetic flux Φ for various coupling strength starting with λ being of order unity, where s(λ, Φ) is a smooth function of Φ. At values of λ ≈ 10 and larger, however, s(λ, Φ) becomes more irregular and shows typical fractal features such as self-similarity. As seen in Fig. 6, this property of s(λ, Φ) becomes more pronounced with increasing λ. A similar behavior is found forē(λ, Φ), as illustrated in Fig. 7. Fig. 8 shows the entropy s(λ, Φ) and the energyē(λ, Φ) as a function of λ for a few representative values of Φ. As seen, both quantities decrease monotonically with λ, and in this sense λ qualifies as a phenomenological (inverse) temperature scale. However, if λ were the true inverse thermodynamic temperature, standard thermodynamics would require it to equal the derivative where the derivatives with respect to λ on the r.h.s can be expressed as In particular, for small λ 1 one finds from Eqs. (55) and (56) ∂s independently of the magnetic flux Φ. As seen in the inset of Fig. 8, and also very explicitly by Eq. (63), the derivative (60) is certainly not equal to λ. The reason for this behavior is that the reduced density matrix formulated as ρ red = exp(−H ent )/Z does not match a canonical equilibrium state characterized by an inverse temperature β. Let us therefore rewrite the entanglement Hamiltonian according to where the inverse thermodynamic temperature β(λ, Φ) is determined as a function of λ as follows: Defining the thermodynamic inner energy and the free energy Here the integration constant ln k, k > 0, reflects the freedom of choosing a unit to measure β, i.e. k plays the same role as Boltzmann's constant in standard thermodynamics. The entropy s and energy e =ē/β per unit cell read as a function of small β Remarkably, the inverse thermodynamic temperature scales for λ 1 as λ 3/2 and is not linear in this parameter, differently what would follow from the ansatz ρ red = ρ 0 , cf. Eq. (47). The technical reason for this observation is that ρ 0 is, while being of canonical form, not the linear approximation to ρ red but contains also arbitrary high powers in λ.
At large λ, i.e. weak coupling between the layers, the entropy s characterizing their mutual entanglement will (along with its derivative) eventually vanish. Thus, for λ 1 Eq. (70) can be simplified as leading to where we have used Eq. (54) and have adjusted the constant k consistent with Eq. (72).
In particular, the above relation implies for large λ. We note that the sign on the r.h.s. of Eqs. (78) depends on the inequality (54) and would be different for positiveē(λ, Φ) at λ 1. The inset of Fig. 9 shows β as a function of λ obtained via numerical integration of Eq. (70) for the same flux values Φ as in Fig. 8. The data depends only weakly on Φ and follows the power law (75) at small λ while the behavior at large λ is well  (72) with an appropriately small but finite initial argument 0 < λ 1.
described by a logarithmic dependence. In the main panels of Fig. 9 we have plotted the entropy s(β, Φ) and the energy e(β, Φ) (as opposed toē) as a function of β. for large β, e(β, Φ) levels off and converges (slowly) to e = −1/k, consistent with Eq. (78). Such a saturation of the thermodynamic energy should be expected since both layers get more and more decoupled in this limit. The originally defined energyē = βe does, in accordance with Eq. (77), not show such a behavior but decreases unboundedly with increasing λ, as shown already in Figs. 7 and 8. While approaching their ground state values in the limit λ → ∞ the inverse thermodynamic temperature β(λ, Φ) as well as the inner energy e(λ, Φ) per unit cell display fractal features, as illustrated in Fig. 10. An analogous behavior was found for the entropy in Fig. 6. In summary, the parameter λ can be viewed as a phenomenological inverse temperature scale having intuitive properties like the decrease of energy and entropy with increasing λ.
The thermodynamic inverse temperature fulfilling standard thermodynamic relations, however, is given by β(λ). Eq. (70) establishes a 1-to-1 mapping between both quantities and is completely universal in the sense that its derivation does not depend on any detail of the underlying system or its entanglement Hamiltonian. Indeed, a relation very similar to Eq. (70) (with in fact an identical l.h.s.) is obtained in standard thermodynamics when connecting the thermodynamic temperature to a phenomenological temperature scale obtained from, say, the change in volume of a given body upon changing its internal energy [53,54]. Also the results (76) to (78) are (up to the sign (54) ofē) very general since they only rely on the vanishing of the derivative of the entanglement entropy s(λ) in the limit of weak coupling λ 1 (⇔ β 1). This statement is of course just the analog of the third law of thermodynamics which requires the entropy to approach a constant in the limit of zero temperature. In the same limit, the thermodynamic energy e per unit cell reaches according to Eq. (78) a saturation corresponding to the ground state energy in classic thermodynamics. Remarkably, within the formalism of entanglement thermodynamics outlined here, this limit is universal and depends only on the constant k, i.e. the unit chosen to measure temperature.

Conclusions and Outlook
We have derived an explicit expression for the entanglement levels of Hofstadter bilayers in terms of the energy eigenvalues of the underlying monolayer system. For strongly coupled layers we find the (expected) proportionality between the entanglement Hamiltonian and the energetic Hamiltonian of the monolayer system with the proportionality factor given by an effective coupling parameter. This parameter, however, is not identical to the inverse thermodynamic temperature, but represents a phenomenological temperature scale. We have devised an explicit relation between both temperature scales which is in close analogy to a standard result of classic thermodynamics. The introduction of the thermodynamic temperature also implies a redefinition of the inner energy. In the limit of vanishing temperature, thermodynamic quantities such as entropy and inner energy approach their ground-state values, but show a fractal structure as a function of magnetic flux.
The relation between the phenomenological temperature scale (given by an appropriate coupling parameter of the total system) and the thermodynamic entanglement temperature applies certainly also to other and more general situations. For instance, in Ref. [20] an entanglement temperature was obtained for bilayer quantum Hall systems by a numerical fit to exact-diagonalization data as a function of layer separation. In the light of the present work, this temperature scale should be seen as a phenomenological one. The analysis of the quantum phase transition performed there, however, should qualitatively not be affected by this issue, since the relation between a phenomenological and the thermodynamic temperature is expected to be a smooth function. On the other and, very explicit investigations as done here for non-interacting particles are obviously more difficult for an interacting system since clearly less analytical tools are available [20].
Similar comments apply to the situation of spin ladders in the limit of strong rung coupling [22,23,24]: The proportionality factor found there between entanglement spectrum and energy spectrum should also be viewed as a phenomenological temperature scale, but not the thermodynamic one. Moreover, Ref. [55] also introduces, studying block entanglement in a spin-1/2 XX-chain, an effective inverse temperature which grows monotonically with increasing block size. Although the situation there and the present one show obvious differences their possible interrelations are of interest.
Finally, another possible extension of the present study is to consider other types of lattices. Given the large deal of interest devoted presently to graphene, the hexagonal geometry [56] and its bilayer versions [57] are obvious candidates . Indeed, it is an interesting speculation whether, for example, the exponent 3/2 occuring in Eq. (75) depends on the lattice geometry.