Entanglement scaling in fermion chains with a localization-delocalization transition and inhomogeneous modulations

For this class of models, the on-site energies are all zero $h_{j}=0$, $j=1,\dots,L$, while the hopping amplitudes are position-dependent, either random or follow an aperiodic modulation. In the first case, which we refer to as random model, we assume that the amplitudes $t_{j}$ are independent, identically distributed quenched random variables, drawn from a uniform distribution in the interval $[0,1]$.

In the latter case, we will use two different aperiodic modulations, both defined by an inflation rule. One of them is the Fibonacci modulation, where hopping amplitudes are modulated according to the Fibonacci sequence. It is defined using a two-letter alphabet ($a$ and $b$) by the substitution rule:

The first few realizations of the Fibonacci sequence are $a$, $ab$, $aba$, $abaab$.

Another two-letter sequence, which we refer to as relevant aperiodic modulation (RAM), is obtained by the following inflation rule jz :

Generating a sufficiently long aperiodic sequence by the repeated application of the inflation rule, a modulation pattern can be associated with it, in which letter $a$ ($b$) corresponds to amplitude $t_{a}$ ($t_{b}$). The strength of the modulation can be characterized by the ratio of two types of amplitudes, $r=t_{a}/t_{b}$.

The non-zero energy eigenstates of the randomly disordered model are exponentially localized, and the localization length diverges for zero energy as it was shown with rigorous tools in connection with the off-diagonal Anderson model in deylon1986 ; igloi2017 . According to the strong-disorder renormalization group (SDRG) method mdh ; im , the ground state of the random model is a random-singlet state fisherxx , which is a product of one-particle states $\frac{1}{\sqrt{2}}(|10\rangle-|01\rangle)$ on pairs of sites which can be arbitrarily far away from each other. The method is approximative but is asymptotically exact, giving the low-energy (large-scale) properties of the system correctly. Analogous to this, the ground state of aperiodic models is an aperiodic-singlet state hida , for any $r<1$ in the case of the RAM jz , but, for the Fibonacci modulation, which is a so-called marginal perturbation, only in the limit $r\to 0$.

In such a singlet state, the entanglement entropy in units of $\ln 2$ is simply given by the number of pairs with precisely one constituent in subsystem $A$. The average entanglement entropy of a subsystem of size $\ell$, which is part of an infinite system increases asymptotically (apart from log-periodic oscillations for aperiodic models) as

where the effective central charge depends on the type of modulation. For the random modelrefael , $c^{\rm ran}_{\rm eff}=\ln 2$ , for the RAM, $c^{\rm RAM}_{\rm eff}=\frac{6\ln 2}{\ln\lambda}\frac{(\lambda-3)^{2}}{2+(\lambda-3)^{2}}$, where jz $\lambda=\frac{1}{2}(3+\sqrt{17})$, while for the Fibonacci modulation it varies continuously with $r$, and its limiting value is $\lim_{r\to 0}c^{\rm FM}_{\rm eff}(r)=\frac{2}{(\tau^{2}+1)\log_{2}\tau}$, where $\tau=\frac{1+\sqrt{5}}{2}$ is the golden mean jz ; ijz .

The logarithmic negativity in a singlet state is given by the number of singlets connecting $A$ and $B$, in units of $\ln 2$. As it was shown in Ref. ruggiero2016 for the random singlet state, the average logarithmic negativity of adjacent intervals of size $\ell$ scales as

with the prefactor $\kappa=\frac{\ln 2}{6}$. Recently, a more detailed study has appeared about the negativity spectrum of this model turkeshi2019 . The prefactor $\kappa$ is the half of the prefactor of entanglement entropy, which can be intuitively understood since subsystem $A$ borders to $B$ only on one side, while, in the case of entanglement entropy, $A$ borders to the rest of the system on both sides. According to this, the result in Eq. (7) is expected to hold also for aperiodic singlet states with a prefactor $\kappa=\frac{c_{\rm eff}}{6}$.

The recapitulated SDRG method and its modifications has a wide range of applications from the description of entanglement of star-graph like systems juhasz_ober_2018 and disordered surfaces juhasz2017 ; juhasz2017_jstatmech including the dynamics of disordered systems pekker2014 ; vosk2013 ; vosk2014 ; bardason2012 ; roosz2017 and description of their highly excited states pouravarni2013 ; pouravarni2015 ; huang2014 ; vasseur2015 ; storms2014 ; hsin2015 or even a construction of quasi-periodic tensor network model enable us to investigate AdS/CFT correspondence jahn2019 .

Another model, which we will consider is the Harper model. Here, the couplings in the general Hamiltonian in Eq. (5) are constant, $t_{i}=1$, while the local potential $h_{j}$ is a quasiperiodic function of $j$, namely,

where $\tau=\frac{1+\sqrt{5}}{2}$ is the golden mean. The irrationality of $\tau$ makes the model quasiperiodic. This model, when formulated in terms of spin variables is also known as Aubry-André model aubry-andre . It is well-known that this model shows a delocalization transition, which is exactly at $h=1$ due to its self-dual property wilkinson . In the region $h<1$, all one-particle eigenstates are extended, whereas, for $h<1$, they are all localised on a length aubry-andre

At the critical point, $h=1$, the one-particle states show an interesting multi-fractal behaviour evangeleou2000 ; siebesma1987 : they are essentially localised on $\sim L^{D_{2}(n)}$ sites, where the dimension $D_{2}(n)$ of the effective support of the state varies from state to state. Its maximal value is found numericallyevangeleou2000 to be $D_{2}\approx 0.82$ .

To our knowledge the logarithmic negativity in the ground state of the Harper model was not studied so far in the literature, the entanglement entropy has been studied in harper_ent , there the authors focused on the effect of the chemical potential. Here we investigate the size dependence of the entanglement entropy.

First, we consider the upper bound $\mathcal{E}_{u}$. It is formulated in terms of the covariance matrix

where $C_{ij}=\langle c^{\dagger}_{i}c^{\phantom{\dagger}}_{j}\rangle$ is the correlation matrix, with all the other entries of $\gamma$ being zero. For every Hamiltonian in the form (5), the covariance matrix can be obtained following a standard canonical transformation lieb61 . For translational invariant systems, one can easily obtain a closed form of the covariance matrix while, for inhomogeneous systems, it is computable in polynomial time in the number of fermionic modes $L$. The covariance matrix characterizes all correlations in the system, so it also determines the entanglement negativity of any subsystems. However, no simple formula is known for the entanglement negativity in terms of the covariance matrix, making the bounds defined in eisler2018 really valuable. The upper bound we use is given by

where $S_{\alpha}(\rho)$ denotes the Rényi entropy of state $\rho$:

The state $\rho_{\times}$ is defined by its covariance matrix as

where $\gamma_{\pm}=T^{\pm}_{B}\,\gamma\,T^{\pm}_{B}$, and $T^{\pm}_{B}=\bigoplus_{j\in A}\mathds{1}_{2}\bigoplus_{j\in B}({\pm i})\mathds{1}_{2}$.

For constructing the lower bound, the matrix $\Gamma=2C-1$ is divided in the following way:

Using the singular value decomposition of $\Gamma_{AB}$, $\Gamma_{AB}=UDV^{T}$, where $D$ is a diagonal matrix with non-negative elements, whereas $U$ and $V$ are orthogonal matrices, one can transform $\Gamma$ by $U\oplus V$ to the following form

Denoting the diagonal elements of $D$, $U^{T}\Gamma_{AA}U$, and $V^{T}\Gamma_{BB}V$ by $c_{i}$, $a_{i}$, and $b_{i}$, respectively, the lower bound is given by

where

Here, we show that the expression in Eq. (16) can be further simplified making use of the sublattice symmetry, which holds in the absence of an on-site potential for an even $L$. We also assume, that the lattice is half filled. Due to this, the elements of matrix $\Gamma=2C-1$ with indices of the same parity are all zero. By replacing rows and columns of $\Gamma$, let us arrange, separately in block $A$ and $B$, the odd (even) indices to the first (second) $l/2$ places. Then $\Gamma_{AB}=2C_{AB}-1$ will have the form

We are looking for the singular value decomposition $\Gamma_{AB}=UDV^{T}$, where the diagonal matrix $D$ contains the non-negative singular values of $\Gamma_{AB}$, while the columns of $U$ and $V$ are eigenvectors of $\Gamma_{AB}\Gamma_{AB}^{T}$ and $\Gamma_{AB}^{T}\Gamma_{AB}$, respectively. These matrices are then block diagonal,

and, as a consequence, $U$ and $V$ are block diagonal, as well:

Transforming $\Gamma$ by $U\oplus V$ will bring the diagonal blocks

to the form

The elements $a_{i}$ and $b_{i}$ in Eq. (16) are therefore zero. Using this, the lower bound can be written as

where the prime means that the summation goes over the singular values fulfilling $c_{i}>\sqrt{2}-1$ only.

First, we investigated the behavior of the lower and upper bounds in the homogeneous chain. Here, we assume that the system is infinite ($L\to\infty$), which allows us to use the exact expressions of the correlations $C_{i,j}=\langle c^{\dagger}_{i}c_{j}\rangle$ igloi2020 , while keep the subsystem size $\ell$ finite. The numerical results are shown in Fig. 1.

The upper bound follows the scaling

with $\kappa_{u}=0.25\pm 0.0005$, which is compatible with the behavior of the entanglement negativity known from CFT, see Eq. (4).

The behavior of the lower limit is more complicated: the overall logarithmic trend is decorated with log-periodic oscillations of small amplitude, which is unusual in a homogeneous system. Here, the oscillations are the consequences of the definition of the lower limit [see Eq. (30)], in which only the singular values greater than $\sqrt{2}-1$ contribute to the sum. Putting the singular values in decreasing order, one can observe that the $n$th singular value ($n=1,2,\dots$) slowly increases with increasing $\ell$. When a singular value crosses the threshold $\sqrt{2}-1$, so that the number of terms in the sum increases by one, a new period of the oscillations starts.

To investigate the overall logarithmic trend, one has to fit to identical parts of the periods, like the breaking points indicated by arrows in Fig. 1. To do this, at least a few period long data set is needed. The period of these oscillations is about $3.2$ on a logarithmic scale, which means that one has to handle quite large subsystems of size $10^{5}\dots 10^{6}$. Fortunately, in this size range, only the largest $1\dots 4$ singular values are needed, which can be obtained from the largest eigenvalues of the matrices $QQ^{T}$ and $P^{T}P$. Since we have a closed formula for the matrix elements of $Q$ and $P$, we can use the power method and multiplicate with $QQ^{T}$ ($P^{T}P$) without storing the matrix elements in the memory. The ratios of the consecutive eigenvalues (in decreasing order) of the matrix $QQ^{T}$ ($P^{T}P$) are between $0.25$ and $0.1$, so the power method converges rapidly.

Using the data for the lower bounds at the four breaking points shown in Fig. 1, we calculated effective prefactors from neighboring data points, which are in order, $0.1321$, $0.1298$, and $0.1288$. Assuming corrections of the form $1/\ln\ell$, an extrapolation to $\ell\to\infty$ results in $\kappa_{l}=0.1256$; if the data point obtained from the smallest system sizes is excluded, one obtains $0.1246$, which gives an estimate on the error of the fit. Thus, in the homogeneous chain, the prefactor of the lower bound, $\kappa_{l}=0.125(1)$, is close to the half of the prefactor of entanglement negativity $1/4$.

The dependence of the lower and upper bounds on $L$ for fixed ratios $\ell/L$, as well as the corresponding effective prefactors for the random model are shown in Fig. 2.

As can be seen, both the lower and upper bounds scale logarithmically with $L$ for large $L$. The effective prefactor of the upper bound displays a slow crossover from the clean system’s value $1/4$ toward $\ln(2)/6$ predicted by the SDRG method with increasing $L$, and, at the system sizes available by the numerical method, it is still considerably far from it. We note that, under the same circumstances, the prefactor of the entanglement entropy has similar deviations from the asymptotic value (not shown). In the case of the lower bound, the prefactor in the homogeneous system ($0.13$) is closer to the asymptotic limit $\ln(2)/6$, and, accordingly, it shows a more rapid crossover.

In the case of the relevant aperiodic modulation, the lower and upper bounds, as well as the effective prefactors as a function of $L$ are exemplified in Fig. 3 for $r=0.5$. The system sizes were $L=17,61,217,773,2753$, which are the lengths of words obtained by the repeated application of the inflation rule starting with letter $a$, while the subsystem size was $\ell=[L/8]$, where $[\cdot]$ stands for the integer part. By this choice of $L$, log-periodic oscillations in the data can be avoided jz . As can be seen, the bounds follow a logarithmic scaling and their asymptotic prefactors shown in the inset for different disorder strengths $r$ are in agreement with the expectation $\kappa_{l,u}=\kappa=c_{\rm eff}/6$ valid for singlet states.

For the Fibonacci modulation, which is a marginal one, we calculated the size dependence of lower and upper bounds for different modulation strengths $r$. Here, the system sizes and the subsystem sizes were chosen to be $L=F_{n}$ and $\ell=F_{n-4}$, respectively, where $F_{n}$ denotes the $n$th term of the Fibonacci sequence. As can be seen in Fig. 4, the bounds plotted against $\ln L$ show log-periodic oscillations with a period of three data points. This is in accordance with that the self-similarity of the aperiodic singlet state is achieved by applying the third power of the inflation transformation jz . Therefore, in estimating the prefactors, we used every third data points. For both bounds, we find a logarithmic dependence on $L$ as given in Eq. (31). The prefactor of the upper bound shows a variation with $r$ qualitatively similar to that of the prefactor of the entanglement entropy jz : For moderate strengths of the modulation, it has a slight deviation from the clean system’s value, then it tends rapidly to the singlet-state limit. In the case of the lower bound, the prefactor in the clean system ($0.125(1)$) and in the singlet-state limit ($0.1327\dots$) hardly differ, therefore the variation with $r$ is very weak.

In the Harper model, we investigated first the scaling of the entanglement entropy in the ground state. In these calculations, system sizes were chosen to be twice of Fibonacci numbers $L=2F_{n}$, and the size of the subsystem was $F_{n-4}$.

The results are shown in Fig. 5. We find that, in the extended phase, the entanglement entropy scales logarithmically as $S=0.33\ln L+\textnormal{const}$, and the measured prefactor ($0.33$) is compatible with that of the homogeneous system $1/3$ up to the precision of the estimation. The reason of this agreement is the extended nature of the eigenstates, which is similar to the eigenstates of a homogeneous chain.

In the critical point, the entanglement entropy is still found to scale logarithmically as

however, the effective central charge $c_{\rm eff}=0.78$ differs significantly from the central charge of the homogeneous system.

In the localized phase, the entanglement entropy saturates to finite values in the limit $L\to\infty$. As it is demonstrated in Fig. 5, in large systems $L\gg l_{\rm loc}$, in which the length scale is set by the localization length $l_{\rm loc}$ rather than the system size, the entanglement entropy follows the law

with the same prefactor as found in Eq. (34). Approaching the critical point, this law is, however, deteriorates when the diverging localisation length becomes comparable with the system size.

Numerical results on the entanglement negativity are shown in Fig. 6. In the delocalized phase and in the critical point, the upper bound increases logarithmically with the system size. In the delocalized phase, just like for the entranglement entropy, the prefactor is found to agree with the that of the conformally invariant homogeneous system. In the critical point, the prefactor of the upper bound is reduced to $\kappa_{u}=0.212$. As can be seen in Fig. 6, the size dependence of the lower bound shows oscillations which may originate both in the discreteness in its definition (see section IV.1) and in the quasiperiodic modulation. For this reason, a completely reliable estimate of the prefactor cannot be obtained.

As can seen in Fig. 6b, the upper and lower bounds become nearly equal in the localized phase. Their ratio approaches to one, if $h\to\infty$, and the ratio is $1.05$ already for $h=1.5$. This makes possible to investigate the behavior of the negativity in the localized phase. Although the difference of the upper and lower bounds becomes larger as one approaches the transition point, there is an intermediate regime, where the transition point is not too far, but the difference of the bounds is still moderate. Similarly to the behavior of the entanglement entropy in the localized phase, the upper bound of entanglement negativity increases according to

in the vicinity of the transition point, where $\kappa_{u}=0.212$ is the prefactor characteristic of the critical point.