Detecting quantum critical points in the t-t′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t'$$\end{document} Fermi-Hubbard model via complex network theory

A considerable success in phenomenological description of high-Tc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {high-T}_{\text{c}}$$\end{document} superconductors has been achieved within the paradigm of Quantum Critical Point (QCP)—a parental state of a variety of exotic phases that is characterized by dense entanglement and absence of well-defined quasiparticles. However, the microscopic origin of the critical regime in real materials remains an open question. On the other hand, there is a popular view that a single-band t-t′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t'$$\end{document} Hubbard model is the minimal model to catch the main relevant physics of superconducting compounds. Here, we suggest that emergence of the QCP is tightly connected with entanglement in real space and identify its location on the phase diagram of the hole-doped t-t′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t'$$\end{document} Hubbard model. To detect the QCP we study a weighted graph of inter-site quantum mutual information within a four-by-four plaquette that is solved by exact diagonalization. We demonstrate that some quantitative characteristics of such a graph, viewed as a complex network, exhibit peculiar behavior around a certain submanifold in the parametric space of the model. This method allows us to overcome difficulties caused by finite size effects and to identify precursors of the transition point even on a small lattice, where long-range asymptotics of correlation functions cannot be accessed.

Scientific Reports | (2020) 10:20470 | https://doi.org/10.1038/s41598-020-77513-0 www.nature.com/scientificreports/ gas in two-dimensional CuO planes, while the conductivity in the orthogonal direction is suppressed. Other phenomena, such as the formation of Fermi arcs seen in the angle-resolved photoemission spectra of high-Tc compounds 19 , or charge density waves 20 also fit pretty naturally into the context of quantum criticality. The main problem of this approach is its purely phenomenological character. It cannot explain by itself why the high-Tc compounds, contrary to the most of interesting condensed matter systems, do not behave as the Fermi liquid but instead are characterized by minimal quantum viscosity and other fancy properties. Such an explanation requires an analysis of electronic structure of specific materials.
Current understanding of high-T c superconductivity in cuprates assumes a crucial role of strong electron correlations 14,[21][22][23] , which are taken into account within a particular minimal model that was formulated 24 on the basis of the density functional band structure of cuprates,-the single-band t-t ′ Hubbard model on a square lattice given by the Hamiltonian where, the first sum is taken over the pairs i, j of nearest neighbors, the second one-over the pairs l, k of next-to-nearest (diagonal) neighbors, c i,σ is the electron annihilation operator, and the on-site occupation operator is n i,σ = c † i,σ c i,σ . For convenience, in what follows we express all energies in units of t . In an attempt to connect the phenomenological and the microscopic levels of description of HTSC, we shall focus on the t-t ′ Fermi-Hubbard model and try to detect QCP on its phase diagram considering the range of parameters typical for cuprates.
Correlation effects beyond the band structure approximation in this model have been thoroughly analyzed with different methods [25][26][27][28][29] , and there are a number of good indications that it captures all the relevant features of cuprate superconductors. In a series of papers [30][31][32][33][34][35] , perturbative renormalization group studies of the model have been conducted, and the emergence of the superconducting order parameter and the competition between superconductivity and antiferromagnetism were demonstrated. In particular 34 , it was argued that the next-tonearest neighbor hopping t ′ plays a crucial role in the stabilization of superconductivity. A complementary approach is based on the cluster dynamical mean-field studies which consider a 2-by-2 plaquette as an elementary unit 36 . Recently 37 , it was noticed that this plaquette has a very special electronic structure for the parameters and the electron occupation number typical for the the optimal doping regime in YBa 2 Cu 3 O 7 ( t ′ = −0.3 , U ≃ 6 ), with an "accidental" degeneracy of many-electron energy levels and formation of the soft fermion mode due to this degeneracy. The pseudogap forms via this mode by a mechanism of the Fano antiresonance, and the superconducting d-wave susceptibility dominates over other instability channels. This behavior was interpreted in terms of formation of a local plaquette valence bond state. On a larger scale, the ground state of the model has been analyzed by means of density matrix renormalization group (DMRG) [38][39][40] (see also 41 for the related studies of its cousin, t − J-model), and additional arguments in favor of stabilization of superconductivity by the next-to-nearest neighbor hopping were provided (in 40 , the competition between long-range and plaquette d-wave superconducting orders has been demonstrated). In turn, at temperatures above the superconducting phase transition, determinantal Monte Carlo computations 42 demonstrated that the DC resistivity exceeds the Mott-Ioffe-Regel limit and scales linearly with temperature.
The search for the QCP in the t-t ′ Hubbard model has been performed within the dynamical cluster approximation 27 , and its existence has been proven by studying thermodynamics properties of the model at finite temperature and their further extrapolation to T = 0 . However, it is tempting to get a deeper insight into the microscopics of the QCP and demonstrate its emergence due to interactions of electrons at low temperatures.
Since large scale simulations of the fermionic Hubbard model away from half-filling are challenging because of the sign problem, it is natural to ask whether we can extract any information about the tendency to form critical states out of small cluster solutions obtained by means of exact diagonalization. At first, this goal does not seem realistic since studying systems in the critical regime unavoidably requires dealing with long-range correlations, while all the microscopic precursors of the transition on small lattices would be washed out by the finite-size effects. However, it is useful to bear in mind that, in the context of many-body quantum dynamics, the concept of entanglement and the phenomenon of collective emergence go hand in hand. An archetypical example of such relation is the Cooper pairs in the BCS theory of superconductivity: while the ground state wavefunction has a form of a product state of the Cooper pairs, each pair itself is a two-body entangled system. Therefore it is natural to expect that major transitions in phenomenological properties of many-body systems would be reflected in the patterns of entanglement, and quantum criticality should leave its fingerprint on all scales, not only in the deep infrared limit. A nice example of how fruitful this way of thinking can be was given in Refs. 43,44 , where entanglement measures were used to determine universality class of the Mott transition in the 2d Hubbard model.
Recently, a novel approach to phase transitions in quantum lattice models based on complex network theory has been suggested 45,46 . It was noticed that a particular structure that can be computed with relative ease and appears to be very sensitive to reconfigurations of the quantum state is the network of quantum mutual information. The mutual information between two subsystems A and B of a larger systems is defined as where S A = − Tr ρ A log ρ A is the von Neumann entropy, and ρ A = TrĀρ is the density matrix of subsystem A with the trace being taken over the the subsystem Ā which is complementary to A, see section Methods for details. Then we can associate a weighted graph with a state of a quantum lattice system, e.g. the Hubbard model, by considering the lattice sites i = 1 . . . N , where N is the number of sites, as nodes of the graph, and the values of pairwise inter-site mutual information I ij play the role of weights on the graph links (see Fig. 1). This

Scientific Reports
| (2020) 10:20470 | https://doi.org/10.1038/s41598-020-77513-0 www.nature.com/scientificreports/ representation is appealing for the following reason. Once a wave function on the lattice is known, it is easy to compute the entanglement entropy of a pair of sites and thus the mutual information. At the same time, such a network contains information of quantum correlations which could be very important to understand the dynamics of strongly correlated systems. In the cases of the transverse field Ising and the Bose-Hubbard models in 1d, it was demonstrated that certain characteristics of the mutual information network can be used to detect quantum phase transitions 45,46 . Namely, behavior of the following functions upon changing parameters of the models has been studied: • Clustering of a weighted graph is defined as where N is the total number of sites in the lattice, and I is the N × N matrix of inter-site mutual information.
One can see that this quantity maximizes on graphs with a lot of three-link loops with high weights. For the cases studied in Ref. 45 , it was shown that it serves as sensitive detector that exhibits a clear dip at the phase transition point. A natural explanation of this fact is that, at the criticality, one can expect the corresponding network to be scale-free, and for generic scale-free networks clustering is usually quite low 47 . • Disparity of a single node in a network is defined as a measure to capture how non-uniformly weights on the links attached to this node are distributed: For example, if the node has the same value of mutual information with all the other nodes of the network, its disparity would be Y i = 1/(N − 1) , while if it correlates only with one neighbor, the disparity maximizes as Y i = 1 . Physically speaking, high disparity of a lattice site means that it tends to correlate only with a few other sites, and "factorize out" of the rest of the system. In the context of quantum many-body physics such a behavior would be typical for states that can be nearly decomposed into product states. On the other hand, low disparity means that the site correlates with a large number of degrees of freedom. • Density is an overall characteristic of a network given by i.e. it is the averaged fraction of all the weights (mutual information values) of the network. To gain more intuition on what properties of the many-body quantum state it reflects, we shall estimate an upper bound i.e. the mutual information network is generally sparse even if the system is highly entangled. Note that bound (6) can be saturated in physically very distinct cases. D is maximal if either each single site is maximally entangled with just one partner site, and the state as a whole decomposes into a product of Bell pairs, or if the entanglement between the site and the rest of the system is homogeneously scrambled over all the sites. To distinguish between such configurations one has to refer to the disparity which we defined above. • Pearson correlations measure how much two nodes i and j of a network differ from each other: In Ref. 45 Pearson correlations of neighboring nodes were shown to develop a cusp around the phase transition point. For one-dimensional Ising and Bose-Hubbard models 45 , this approach to detecting quantum phase transitions points was successfully applied for systems of ∼ 10 2 sites, and was demonstrated to be very robust upon finite-size effects. In the two-dimensional case, we are limited by much smaller system sizes (we perform exact diagonalization for a 4-by-4 plaquette), and should not expect our results to be free from finite-size artifacts. Still, as we shall see in the next section, the network measures exhibit clearly distinguishable features on a submanifold of the t-t ′ Hubbard model parametric space. In particular, this submanifold includes the level-crossing point observed in a 2-by-2 plaquette for the choice of parameters corresponding to YBa 2 Cu 3 O 7 superconductor 37 .

Results.
We have computed the complex network measures discussed above across the space of parameters of a 4-by-4 t-t ′ Hubbard plaquette. Within each fixed particle number sector, from (5, 5) ( 37.5% hole doping) to (7, 7) ( 12.5% hole doping), we scan over t ′ and U, Fig. 2. As an indicative value, we take t ′ = −0.3 , which is estimated to be the next-neighbor hopping in the Hubbard model of YBaCuO compounds, and search for transition points around it. The temperature is fixed to 1/T = β = 100 (all energies are expressed in the units of |t|), and the system is studied in the canonical ensemble.
We assume that a transition point is evident if all the measures exhibit some clear features around the same point. Accepting this criterion, we can claim with a high confidence that, for non-periodic boundary-conditions, there is a family of transition points in each sector (except (7, 7)) forming a nearly perfect straight line in the t ′ -U plane that extends in a certain range of t ′ , Fig. 3 (for a more detailed picture of how complexity measures behave at different values of t ′ (see Supplemental material). For too small |t ′ | , the signs of criticality are faded away from the complexity measures. Moreover, for different values of hole doping, all these lines have very similar slope. This can likely be interpreted as an indication that, in the thermodynamic limit, quantum phase transition occurs on a 2d manifold in the 3d parametric (U/t, t ′ /t , particle number) space of the model.
We stick to non-periodic boundary conditions for the reason that the mutual information network has a richer structure in that case. If periodic boundary conditions are imposed, all lattice sites are identical, and every site has only five inequivalent connections to others, making the I ij matrix highly degenerate. Hence, the corresponding network structures are constrained by symmetries and much less sensitive to variations of the model parameters. Still, we would like to note that, when boundary conditions are changed for periodic ones, all the transition lines are smeared out with the only exception of the (6, 6) sector which corresponds to the hole doping of δ = 25% . For the latter, only the concrete values of Coulomb repulsion U gets shifted (see Supplemental material). While we do not expect the information network constructed for periodic boundary to be sensitive enough to properly detect precursors of phase transitions, it is interesting to note that this single sector where the transition is evident for both choices of b.c. is the same as the one where level-crossing associated with formation of the pseudogap via Fano antiresonance occurs in a 2-by-2 plaquette 37 .
At the same time, in the density of states (d.o.s.) the transition point is (almost) invisible. Some minor peculiarity at the quantum critical point is visible in the density of states at t ′ = −0.3 . Around the transition point identified by means of the complex network theory ( U = 7.5 , sector (6,6)) the ω = 0 peak in the d.o.s. starts splitting and the pseudogap emerges, see Fig. 4. Further decrease of the hole doping leads to enhancement of the gap. The particular role of U in this transition is less clear, as the d.o.s. profile varies very mildly upon changing U, and it would be safer to claim that the spectral properties are not sensitive to the discussed quantum phase transition. Ideologically, this situation is somewhat similar to the Anderson localization in disordered systems which is a clear example of a phenomenon that cannot be detected on the level of the average Green's functions 49 .

Discussion
By associating the quantum state of the t-t ′ Hubbard model with a weighted network of inter-site mutual information, for different values of the next-neighbor hopping t ′ , we have found a set of transition lines in the U-t ′ plane of the model parametric space, where characteristics of the network have a clearly distinguishable feature. Such a behavior was previously shown to be an indication of quantum phase transitions in different one-dimensional models 45,46 . The modern experimental understanding of the putative QCP in cuprates tells that it must be associated with the emergence of the pseudogap phase 4 . For example, for YBaCuO compounds the onset of pseudogap was experimentally demonstrated to happen at hole doping δ ≃ 22% 50 . The hole doping δ = 25% is the closest value one can get for a 4-by-4 cluster (the (6, 6) sector), and, interestingly, it is exactly the sector where the complex network measures demonstrate the most robust transition features. The particular values of the on-site Coulomb repulsion is affected by the finite size effects, and estimated to be about U ≃ 7−8 . At the same time, no peculiarity is seen in the density of states at the transition point, apart from a slight splitting of d.o.s. around the Fermi level, which might be a good indication that the low-order correlation functions that define the spectral and the response properties of the system could be blind to restructuring of many-body quantum states, and does not contain enough information on the role of quantum correlations behind phase transitions in electron systems. Both the strength and the weakness of the employed approach is that it helps to identify any transition point while being ignorant about its nature. Therefore we can claim that we observe clear short-range precursors likely indicating the existence of a manifold of QCP in the thermodynamic limit of the t-t ′ Hubbard model, but www.nature.com/scientificreports/ we cannot deduce what order parameter of the corresponding transition is. Still, we tend to relate the observed transition to the critical point discussed in 37 , where it was associated with emergence of soft fermion modes. Within the exact diagonalization approach, we were able to consider only a small cluster, where one unavoidably has to deal with strong finite size effects. It would be interesting to conduct a similar analysis for a larger system. The state-of-the-art DMRG allows to study quasi-2d stripes as large 4-by-64 atomic sites [38][39][40][41][42] , and we hope to apply the complex network approach to such systems in the future. Still, already at this point it is instructive to compare our results with what has been found by means of other approaches. In the recent paper 39 , a detailed DMRG analysis of the Hubbard model phase diagram with hole doping ≤ 12.5% has been conducted (which would correspond to sectors (7, 7), (7,8), (8,8) in our case), and it was shown that for such hole concentration and negative values of t ′ , the model is in the Luther-Emery liquid phase 51 which is very stable upon varying U and t ′ . This is consistent with our observation that the signs of QCP fade away for all values of the next-nearest hopping and Coulomb repulsion when we approach this level of doping. In 52 , position of the metal-insulator phase transition point on the U-µ plane of the Hubbard model phase diagram has been traced out by means of cellural DMFT. In this study, t ′ was set to 0, so it is clearly a different case (as we mentioned before, we tend to interpret the transition as the onset of pseudogap, and not as the Mott transition). Nevertheless, there is a clear similarity in the trend: when the hole doping is lowered, the transition point is shifted towards smaller U. In diagrammatic Monte Carlo study 53 , the point of the pseudogap onset was associated with t ′ = −0.3 . However, they  . The spectral density of states A(ω) defined as the imaginary part of the retarded Green function, computed for U = 7.5 with non-periodic boundary conditions ( ω = 0 corresponds to the Fermi energy). One can see the pseudogap formation around ω = 0 near the quantum critical point. Its interpretation in terms of the Fano antiresonance due to formation of a "soft fermion" mode was given in Ref. 37  www.nature.com/scientificreports/ considered high temperatures ( β = 5 ) and very low hole doping ( 4% ), which makes it problematic to directly compare the results, although it could be claimed that at this hole concentration (approximately (7,8) sector) we do not see any signs of a transition point for any value of U and t ′ (and, accordingly to 39 , the system should be in the Luther-Emery phase), so this controversy is still to be resolved.

Methods
In this section we give the relevant technical details of the calculation of the entanglement measures defined above. The first step is to diagonalize the Hubbard model (1) for a 4-by-4 cluster. This can be done either for a periodic or a non-periodic model. The diagonalization is performed using the Lanczos algorithm with 200 Krylov basis vectors 54 . The particle number and spin conservation laws are used so that the diagonalization can be restricted to a sector with a fixed number of up-and down-spins. Those eigenstates with the corresponding eigenvectors are then used to calculate the reduced density matrices for each possible pair of sites as well as for each single site. The reduced density matrix is computed using its definition that can be symbolically written as: Here a, a ′ denote the many-particle (Fock) basis states describing the subsystem A we calculate the density matrix for, ā stands for the many-particle basis state of the complementary subsystem Ā , thus a couple of those (a,ā) denotes a basis Fock state for the whole cluster explicitly split into two parts. As before, n stands for a particular eigenvector, the density matrices for given eigenstates are weighted with the Boltzmann factors corresponding to their energies. In a given sector for a given set of parameters we use the Boltzmann factor cut-off of 1% meaning e (E 0 −E k )β > 10 −2 , where E 0 is the ground state energy and E k is the energy of the highest (kth) level taken into account. Note that while performing the partial trace over Ā one has to correctly account for the fermionic commutation relations. To this aim one has to effectively change the numeration of sites so that the sites for which we calculate the density matrix stand first. Explicitly it means that each component of an eigenvector, corresponding to a given basis state of the cluster, gets a factor determined as the parity acquired while "dragging" the occupied sites of A to the beginning past the occupied states of Ā . In other words for each basis vector one takes each occupied site from A and for each occupied spin component counts the number of same spin occupied sites from Ā standing before the considered site in the original numeration. Summing up the parities of those numbers for all occupied sites and spins from A one gets the parity that is assigned to a given basis vector with respect to the subsystem A. Having multiplied the eigenvector components with the acquired parities one finally performs the partial trace over the complementary subset Ā . Given the reduced density matrix, we first calculate the von Neumann entropy of a given subsystem and then, with Eq. (2), the mutual information for each pair of sites, that serves as the basis for our network.
The ω-dependent Green function is given by: Here m, n denote the eigenstates of the system, i and σ denote a given site and spin (in the paramagnetic case the answer is spin-independent), E n is the energy of the n-th state, and Z = m e −βE m is the partition function. Note that m and n necessarily belong to different sectors. The Green function is then used to calculate the spectral density of states, which is defined as To perform numerical computations, delta-peaks in the Green function are broadened with δ = π/β.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.