Coulomb electron pairing in a tight-binding model of La-based cuprate superconductors
aa r X i v : . [ c ond - m a t . s up r- c on ] J u l Annalen der Physik, July 27, 2020
Coulomb electron pairing in a tight-binding model ofLa-based cuprate superconductors
K. M. Frahm, and D. L. Shepelyansky ∗ We study the properties of two electrons with Coulombinteractions in a tight-binding model of La-basedcuprate superconductors. This tight-binding model ischaracterized by long-range hopping obtained previ-ously by advanced quantum chemistry computations.We show analytically and numerically that the Coulombrepulsion leads to a formation of compact pairs prop-agating through the whole system. The mechanism ofpair formation is related to the emergence of an effec-tive narrow energy band for Coulomb electron pairswith conserved total pair energy and momentum. Thedependence of the pair formation probability on an ef-fective filling factor is obtained with a maximum arounda filling factor of 20 (or 80) percent. The comparisonwith the case of the nearest neighbor tight-bindingmodel shows that the long-range hopping providesan increase of the phase space volume with high pairformation probability. We conjecture that the Coulombelectron pairs discussed here may play a role in hightemperature superconductivity.
The phenomenon of high temperature superconductiv-ity (HTC), discovered in [1], still requires its detailed phys-ical understanding as discussed by various experts ofthis field (see e.g. [2–4]). The analysis is complicated bythe complexity of the phase diagram and strong interac-tions between electrons (or holes). As a generic model,that can be used for a description of most supercon-ducting cuprates, it was proposed to use a simplifiedone-body Hamiltonian with nearest-neighbor hoppingon a square lattice formed by the Cu ions [5]. In addi-tion the interactions between electrons are consideredas a strongly screened Coulomb interaction that resultsin the 2D Hubbard model [5]. However, a variety of ex-perimental results cannot be described by the 2D Hub- bard model (see e.g. discussion in [6]). Other models oftype Emery [7–10] were developed and extended on thebasis of extensive computations with various numericalmethods of quantum chemistry (see e.g. [6, 11] and Refs.therein). These studies demonstrated the importance ofnext-nearest hopping and allowed to determine reliablythe longer-ranged tight-binding parameters.In this work we use the 2D longer-ranged tight-bindingparameters reported in [6] and study the effects of Coulombinteractions between electrons in the frame work of thistight-binding model. There are different reasons indicat-ing that long-range interactions between electrons maylead to certain new features as compared to the Hub-bard case (see [3, 4, 6]). Recently, we demonstrated thatfor two electrons on a 2D lattice with nearest-neighborhopping the energy and momentum conservation lawsleads to appearance of an effective narrow energy bandfor energy dispersion of two electrons [12]. In such a nar-row band even a repulsive Coulomb interaction leads toelectron pairing and ballistic propagation of such pairsthrough the whole system. The internal classical dynam-ics of electrons inside such a pair is chaotic suggestingnontrivial properties of pair formation in the quantumcase. In this work we extend the investigations of theproperties of such Coulomb electron pairs for a moregeneric longer-ranged tight-binding lattice of one-bodyHamiltonian typical for La-based cuprate superconduc-tors. We find that the long-range hopping leads to newfeatures of Coulomb electron pairs.In Sec. 2 a detailed description of the tight-bindingmodel for two interacting electrons for general latticeswith a particular application to HTC is presented to-gether with an analysis of the effective band width atfixed conserved total pair momentum. Section 3 providesfirst results of the full space time evolution obtained inthe frame work of the Trotter formula approximation. ∗ Corresponding author E-mail: dima(at)irsamc.ups-tlse.fr1 Laboratoire de Physique Théorique, IRSAMC, Université deToulouse, CNRS, UPS, 31062 Toulouse, France
Copyright line will be provided by the publisher . M. Frahm and D. L. Shepelyansky: Coulomb electron pairing in a tight-binding model of La-based cuprate superconductors Section 4 introduces the theoretical basis for the descrip-tion in terms of an effective block Hamiltonian for a givensector of fixed momentum of a pair with technical de-tails provided in Appendix A. In Sec. 5 the phase diagramof the long time average of the pair formation probabil-ity in the plane of total momentum is discussed whileSec. 6 provides some results for the intermediate timeevolution of pair formation. An overview of the results forthe pair formation probability at different filling factorsis given in Section 7. The final discussion is presented inSection 8.
We assume that each electron moves on a square latticeof size N × N with periodic boundary conditions withrespect to the following generalized one-particle tight-binding Hamiltonian: H p = − X r X a ∈ A t a ¡ | r 〉〈 r + a | + | r + a 〉〈 r | ¢ (1)where the first sum is over all discrete lattice points r (measured in units of the lattice constant) and a belongsto a certain set of neighbor vectors A such that for eachlattice state | r 〉 there are non-vanishing hopping matrixelements t a with | r + a 〉 and | r − a 〉 for a ∈ A . To be moreprecise, due to notational reasons, we choose the set A to contain all neighbor vectors a = ( a x , a y ) in one halfplane with either a x > a y > a x = A ′ = A ∪ ( − A ) is the full set of all neighbor vectors. Foreach vector a of the full set A ′ , we require that any othervector ˜ a which can be obtained from a by a reflection ateither the x -axis, y -axis or the x - y diagonal also belongsto the full set A ′ and has the same hopping amplitude t a = t ˜ a .For the usual nearest neighbor tight-binding model(NN-model), already considered in [12], we have the set A NN = {(1,0),(0,1)} with t (1,0) = t (0,1) = t =
1. The nu-merical results presented in this work correspond eitherto the NN-model (for illustration and comparison) orto a longer-ranged tight-binding lattice according to [6]which we denote as the HTC-model. For this case theset of neighbor vectors is A HTC = {(1,0),(0,1),(2,0),(0,2),(1, ± ± ± ± t = t (1,0) = t ′ = t (1,1) = − t ′′ = t (2,0) = t ′′′ = t (2,1 ) = t (4) = t (2,2) = − t = t (1,0) = t (0,1) which is therefore set to unity here; see also Fig. 6a of [6] for the neighbor vectors of the different hopping am-plitudes). The hopping amplitudes for other vectors suchas (0,1), (1, − −
2) etc. are obtained from theabove amplitudes by the appropriate symmetry transfor-mations, e.g. t (1, − = t (1,1) = t ′ = − A and also with a poten-tial generalization to other dimensions.The eigenstates of H p given in (1) are simple planewaves: | p 〉 = N X r e i p · r (2)with energy eigenvalues: E p ( p ) = − X a ∈ A t a cos( p · a ) (3)and momenta p = ( p x , p y ) such that p x and p y are inte-ger multiples of 2 π / N (i.e. p α = π l α / N , l α = N − α = x , y ). For the HTC model, we can give a more explicitexpression of the energy dispersion: E p ( p x , p y ) = − £ cos( p x ) + cos( p y ) ¤ − t ′ cos( p x ) cos( p y ) − t ′′ £ cos(2 p x ) + cos(2 p y ) ¤ − t ′′′ £ cos(2 p x ) cos( p y ) + cos(2 p y ) cos( p x ) ¤ − t (4) cos(2 p x ) cos(2 p y ) (4)which corresponds to eq. (30) of [6] (assuming t = t (5) = t (6) = t (7) = H = H (1)1 p ⊗ (2) + (1) ⊗ H (2)1 p + X r , r ¯ U ( r − r ) | r , r 〉〈 r , r | (5)where H ( j )1 p is the one-particle Hamiltonian (1) of parti-cle j = r j = ( x j , y j ) and ( j ) is the unit operator of particle j . The last term in (5)represents a (regularized) Coulomb type long-range in-teraction ¯ U ( r − r ) = U /[1 + r ( r − r )] with amplitude U and the effective distance r ( r − r ) = p ∆ ¯ x + ∆ ¯ y between the two electrons on the lattice with periodicboundary conditions. (Here ∆ ¯ x = min( ∆ x , N − ∆ x ); ∆ ¯ y = min( ∆ y , N − ∆ y ); ∆ x = x − x ; ∆ y = y − y and the lat-ter differences are taken modulo N , i.e. ∆ x = N + x − x if x − x < ∆ y ). Furthermore, we Copyright line will be provided by the publishernnalen der Physik, July 27, 2020 consider symmetric (spatial) wavefunctions with respectto particle exchange assuming an antisymmetric spin-singlet state (similar results are obtained for antisymmet-ric wavefunctions).In absence of interaction ( U =
0) the energy eigen-values (the classical energy) of the two electron Hamilto-nian (5) (the two electrons) at given momenta p and p are (is) given by: E c ( p , p ) = E p ( p ) + E p ( p ) = − X a ∈ A t a cos( p + · a /2) cos( ∆ p · a ) (6)where p + = p + p is the total momentum and ∆ p = ( p − p )/2 is the momentum associated to the relative co-ordinate ∆ r = r − r . For the NN-model Eq. (6) becomes E c ( p , p ) = − P α = x , y cos( p + α /2) cos( ∆ p α ).Due to the translational invariance the total momen-tum p + is conserved even in the presence of interac-tion ( U
0) and only two-particle plane wave stateswith identical p + are coupled by non-vanishing interac-tion matrix elements. For the case of the NN-model, an-alyzed in [12], the kinetic energy at fixed p + is boundedby ∆ E b = P α = x , y | cos( p + α /2) | . Thus for TIP states with E > ∆ E b the two electrons cannot separate and propa-gate as one pair even if their interaction is repulsive. For p + x = p + y = π + δ being close to π and | δ | ≪ U as soon as ∆ E b ≈ | δ | < U ≪ B with B = + U being the maximal energy bandwidth in 2D. Thusthe conservation of the total momentum of a pair with p + x = p + y ≈ π leads to the appearance of an effectivenarrow energy band with formation of coupled electronpairs propagating through the whole system. However,the results obtained in [12] show that even for other val-ues of p + x , p + y the probability of pair formation is ratherhigh.For the NN-model the effective band width for pairs ∆ E b can be exactly zero for the specific pair momentum p + = ( π , π ). However, this is not the case for the HTC-model where due to the longer-ranged hopping the min-imal width ∆ E b is finite due to the additional terms withfactors cos( p + · a /2) in (6). Therefore, we determined nu-merically for each given value of total momentum p + theeffective bandwidth as: ∆ E b ( p + ) = max ∆ p £ E c ( p , p ) ¤ − min ∆ p £ E c ( p , p ) ¤ (7) B = + U for the band-width of the NN-model. with p = p + /2 − ∆ p and p = p + /2 + ∆ p . Top panels ofFig. 1 show density color plots of ∆ E b ( p + ) for the NN-and the HTC-model. For the HTC-case ∆ E b ( p + ) is max-imal at p + = (0,0) with value ∆ E b ,max = p + = ( π , π ) with value ∆ E b ,min = ∆ E b ,max =
16 at p + = (0,0) and ∆ E b ,min = p + = ( π , π ). The value ∆ E b ,min = ∆ E b ,max ≈
18 and we may expect a somewhatstronger pair formation probability for total momenta p + close to ( π , π ). However, this situation is qualitatively dif-ferent as compared to the NN-model and the HTC-caserequires new careful studies.For comparison, we also show in the lower panels ofFig. 1 the kinetic energy E c at p = p = p + /2 (for thesquare p + ∈ [0, π ] × [0, π ]) corresponding to ∆ p =
0. Whilefor the NN-model this quantity vanishes at p + = ( π , π )there is for the HTC-model a zero-line between the twopoints ( βπ , π ) and ( π , βπ ) where β ≈ ≈ As in [12] the full time evolution of two electrons is com-puted numerically for N =
128 using the Trotter formulaapproximation (see e.g. [12, 13] for computational de-tails). We use the Trotter time step ∆ t = B = + U ) which is the inverse bandwidth for the case of NN-model. A further decrease of the time step does not af-fect the obtained results. At the initial time both elec-trons are localized approximately at ( N /2, N /2) with thedistance ∆ ¯ x = ∆ ¯ y = ∆ x - and ∆ y -axis. The method provides for each time value a wave-function ψ ( x , y , x , y ) from which we extract differentquantities such as the density in x - x plane: ρ X X ( x , x ) = X y , y | ψ ( x , y , x , y ) | (8)or the density ∆ x - ∆ y plane: ρ rel ( ∆ x , ∆ y ) = X x , x | ψ ( x , y , x + ∆ x , y + ∆ y ) | (9)(with position sums taken modulo N ). We also computethe quantity w by summing the latter density (9) overall values such that | ∆ ¯ x | ≤
10 and | ∆ ¯ y | ≤
10 which corre-sponds to a square of size 21 ×
21 in ∆ x - ∆ y plane (dueto negative values of x − x etc.). This quantity gives the Copyright line will be provided by the publisher . M. Frahm and D. L. Shepelyansky: Coulomb electron pairing in a tight-binding model of La-based cuprate superconductors ✲(cid:0)✲✁✂✄✁✁✂✄(cid:0) Figure 1
Top panels show the dependence of the effectiveelectron pair band width ∆ E b ( p + ) on the pair momentum p + = ( p + x , p + y ) . Bottom panels show the kinetic electronpair energy E c ( p , p ) (in absence of interaction) at momenta p = p = p + /2 . Left panels correspond to the NN-modeland right panels to the HTC-model. In all panels the horizon-tal axis corresponds to p + x ∈ [0, π ] and the vertical axis to p + y ∈ [0, π ] . The numbers of the color bar correspond for toppanels to the ratio of the bandwidth over its maximal value andfor lower panels to the quantity sgn( E c ) p | E c | / E c ,max with E c ,max being the maximum of | E c | . In all subsequent color plotfigures the numerical values of the color bar corresponds tothe ratio of the shown quantity over its maximal value. quantum probability to find both electrons at a distance ≤
10 (in each direction) and we will refer to it as the pairformation probability .In Fig. 2 the density ρ X X is shown for U =
2, bothNN- and HTC-models at two time values t = ∆ t and t = ∆ t . These results show that the wavefunction hasa component with electrons separating from each otherand a component where electrons stay close to eachother forming a pair propagating through the whole sys-tem that corresponds to a high density near a diagonalwith x ≈ x . For t = ∆ t the value of w is roughly10% and for t = ∆ t it is roughly 13% for both models.However, the remaining diffusing component of about87-90% probability has a stronger periodic structure forthe NN-model as compared to the HTC-model.Figure 3 shows the density ρ rel ( ∆ x , ∆ y ) for the samecases of Fig. 2. We clearly see a strong enhancement ofthe probability at small values ∆ ¯ x ≈ ∆ ¯ y < < − ✵✵(cid:0)✁✂✵(cid:0)✂✵(cid:0)✄✂✶ Figure 2
2D Wavefunction density ρ X X ( x , x ) in x - x plane (see Eq. (8)) obtained from the time evolution using theTrotter formula approximation for initial electron positions at ≈ ( N /2, N /2) with distance ∆ ¯ x = ∆ ¯ y = for N = , U = and Trotter integration time step ∆ t = B = + U ) .Top (bottom) panels correspond to the time value t = ∆ t ( t = ∆ t ) and left (right) panels correspond to the NN-lattice(HTC-lattice). The corresponding values of the pair formationprobability w are 0.106 (top left), 0.133 (bottom left), 0.0940(top right) and 0.125 (bottom right). Related videos are avail-able at [14, 15]. considerable probability that both electrons stay closeto each other forming a Coulomb electron pair. Further-more, the remaining wavefunction component of inde-pendently propagating electrons, clearly visible in Fig. 2,is not visible in the density shown in Fig. 3 even thoughthis component corresponds to 87-90% probability.The supplementary material contains two videos (for ∼
460 time values in the range ∆ t ≤ t ≤ ∆ t withroughly uniform logarithmic density) of the two densities ρ X X and ρ rel where both models NN and HTC are directlycompared in the same video. The raw-data used for thesevideos is the same as in Figs. 2 and 3. As already mentioned in Sec. 3 the total momentum p + is conserved by the TIP dynamics of the Hamiltonian (5).In order to exploit this more explicitly, we introduce as in Copyright line will be provided by the publishernnalen der Physik, July 27, 2020 ✵✵(cid:0)✁✂✵(cid:0)✂✵(cid:0)✄✂✶
Figure 3
2D Wavefunction density ρ rel ( ∆ x , ∆ y ) in ∆ x - ∆ y plane of relative coordinates (see Eq. (9)) for the same states,cases and parameters of Fig. 2. All panels show the zoomeddensity for ≤ ∆ x , ∆ y < . Related videos are available at[14, 15]. [12], block basis states by: | p + , ∆ r 〉 = N X r e i p + · ( r + ∆ r /2) | r , r + ∆ r 〉 (10)where p + = ( p + x , p + y ) (with p + α = π l + α / N ; l + α = N − α = x , y ) is a fixed value of the total momentum and r , ∆ r are vectors on the square lattice (with position sumsin each spatial direction taken modulo N ). One can show(see Appendix A for details) that the TIP Hamiltonian (5)applied to such state gives a linear combination of suchstates for different ∆ r values but the same total momen-tum value p + which provides for each value or sector of p + an effective block Hamiltonian :¯ h ( p + ) = − X ∆ r X a ∈ A ¯ t ( p + ) a ¡ | ∆ r + a 〉〈 ∆ r | + | ∆ r 〉〈 ∆ r + a | ¢ + X ∆ r ¯ U ( ∆ r ) | ∆ r 〉〈 ∆ r | (11)where ¯ t ( p + ) a = p + · a /2) t a is an effective rescaled hop-ping amplitude depending also on p + and we have forsimplicity omitted the index p + in the block basis states.This effective block Hamiltonian corresponds to a tight-binding model in 2D of similar structure as (1) with mod-ified hopping amplitudes and an additional “potential”¯ U ( ∆ r ). We note that in absence of this external poten-tial ( U =
0) the eigenfunctions of (11) are plane waves and we immediately recover the expression (6) for its en-ergy eigenvalues where ∆ p is the momentum associatedto the relative coordinate ∆ r . For the simple NN-modelthe result for the effective block Hamiltonian was alreadygiven in [12] and the above expression (11) provides thegeneralization to arbitrary tight-binding lattices charac-terized by a certain set of neighbor vectors A and as-sociated hopping amplitudes t a (the generalization toarbitrary spatial dimension is also obvious). As alreadydiscussed in [12] the boundary conditions of (11) in x − ( y − )direction are either periodic if the integer index l + x ( l + y ) of p + x ( p + y ) is even or anti-periodic if this index isodd. This can be understood by the fact that the expres-sion (10) is modified by the factor e ± i p + x N /2 = e ± i π l + x = ( − l + x if ∆ x is replaced by ∆ x ± N and similarly for ∆ y (with ∆ r = ( ∆ x , ∆ y )).Diagonalizing the effective block Hamiltonian (11),we can rather efficiently compute the exact quantumtime evolution | ¯ ψ ( t ) 〉 = e − i ¯ h ( p + ) t | ¯ ψ (0) 〉 inside a given sec-tor of p + . As initial state | ¯ ψ (0) 〉 we choose a state (in thereduced block space) given as the totally symmetric su-perposition of four localized states where ∆ x and ∆ y areeither 1 or N −
1. Such a state corresponds in full spaceto a plane wave in the center of mass direction with totalfixed momentum p + and strongly localized in the relativecoordinate ∆ r . The matrix size of (11) is N which corre-sponds to a complexity of N for the numerical diagonal-ization.However, for a general lattice, such as the HTC-model,one can exploit the particle exchange symmetry to re-duce the effective matrix size to roughly N /2 and for thespecial cases of p + x = p + y or either p + x = p + y =
0a second symmetry allows a further reduction of the ef-fective matrix size to ≈ N /4 (for the NN-model there aretwo or three symmetries for these cases with effective ma-trix sizes of ≈ N /4 or ≈ N /8 respectively; see [12] andAppendix A for details).In view of this, we have been able to compute numer-ically the exact time evolution for the HTC-model in cer-tain p + sectors for a lattice size up to N =
384 for thecase of two symmetries and a limited number of differ-ent other parameters (values of p + and U ). For the caseof one symmetry and the exploration of all possible val-ues of p + x and p + y we used the maximum system size N = N ) that they provide identicalnumerical results.We compute the wavefunction in block representa-tion ¯ ψ ( p + , ∆ r ) for about 700 time values t = − ≤ t / ∆ t ≤ (with a uniform density in logarithmic scale) Copyright line will be provided by the publisher . M. Frahm and D. L. Shepelyansky: Coulomb electron pairing in a tight-binding model of La-based cuprate superconductors where ∆ t = B = + U ) is the time step alreadyused for the Trotter formula approximation given as theinverse bandwidth for the case of the NN-model which isthe smallest time (inverse of the largest energy) scale ofthe system.From the wavefunction we extract in a similar wayas in Sec. 3 the pair formation probability w by sum-ming the (normalized) wavefunction density | ¯ ψ ( p + , ∆ r ) | at fixed p + over the 21 ×
21 square with | ∆ ¯ x | ≤
10 and | ∆ ¯ y | ≤
10. We also compute the inverse participation ra-tio: ξ IPR = ÃX ∆ r | ¯ ψ ( p + , ∆ r ) | ! − (12)which gives roughly the number of lattice sites (in ∆ r space) over which the wavefunction is localized. Bothquantities w and ξ IPR converge typically rather well totheir stationary values at times t > ∆ t with some timedependent fluctuations. Therefore for the cases wherewe are interested in the long time limit we compute thewavefunction only for 70 times values (in the same inter-val as above with uniform logarithmic density) and takethe average over the 21 values with 10 ≤ t / ∆ t ≤ . Wenote that for the case of a uniform wavefunction densitythe ergodic values are w = (21/ N ) and ξ IPR,erg = N .Values of w significantly above w or of ξ IPR below ξ IPR,erg indicate an enhanced probability for the forma-tion of compact electron pairs.We also mention that both quantities w and ξ IPR are invariant with respect to the three transformations p + x ↔ p + y , p + x → − p + x and p + y → − p + y (or p + x → π − p + x and p + y → π − p + y ) corresponding to reflec-tions at the x - y diagonal, the y -axis and the x -axis. Eventhough the effective block Hamiltonian (11) is not (al-ways) invariant with respect to all three of these transfor-mations (see Appendix A for details), the choice of an in-variant initial state ensures that at finite times the wave-function in block space satisfies for example the identity¯ ψ ( p + x , p + y , ∆ x , ∆ y ) = ¯ ψ ( p + y , p + x , ∆ y , ∆ x ) (and similarlyfor the other reflections). In other words a certain reflec-tion transformation for p + results in the equivalent trans-formation for the time dependent block space wavefunc-tion in ∆ r space. Obviously the two quantities w and ξ IPR do not change with respect to these transformations(in ∆ r space) and therefore they are conserved. As a re-sult it is sufficient to compute these quantities only forthe triangle 0 ≤ p + y ≤ p + x ≤ π .In the following sections we present the results forthese quantities and the wavefunction in block represen-tation. ✵✵(cid:0)✁✂✵(cid:0)✂✵(cid:0)✄✂✶ Figure 4
Phase diagram of electron pair formation in theplane of pair momentum p + = ( p + x , p + y ) for the NN-lattice(left panels), the HTC-lattice (right panels) and the interactionvalues U = (top panels), U = (bottom panels). Shownis the pair formation probability w for N = obtainedfrom the exact time evolution for each sector of p + with aninitial electron distance ∆ ¯ x = ∆ ¯ y = and computed from anaverage over 21 time values in the interval ≤ t / ∆ t ≤ . In all panels the horizontal (vertical) axis corresponds to p + x ( p + y ) ∈ [0, π ] and the numerical values of the color barcorrespond to the ratio of w over its maximal value. Themaximum values corresponding to the red region at the topright corner p + = ( π , π ) are w = (both left panels), w = (top right) and w = (bottom right). For com-parison the ergodic value is w = (21/192) = . The phase diagram of the long time average of the pairformation probability w in the p + -plane is shown inFig. 4 for both models and the interaction values U = p + = ( π , π ) and minimal at p + = (0,0).Furthermore, the size of the maximum region is signifi-cantly stronger for U = U = p + values even a rela-tively weak or moderate Coulomb repulsion creates quitestrongly coupled electron pairs.For the NN-model the top ( p + y = π ) or right ( p + x = π )boundary also provide large values with w ≈ U = Copyright line will be provided by the publishernnalen der Physik, July 27, 2020 U = U = U = ≈ U = r g = q p + x + p + y ≈ π for U = w ≈ U = w ≈ r g ≈ π . This circle seems to be less pro-nounced despite its larger value of w as compared to U = U = w ≈ p + = ( π , π )) is roughly twice the maxi-mum value for U = w ≈ w at p + ≈ (0,0) are w ≈ − w ≈ − U = U = p + , for both mod-els and both interaction values U = ξ IPR which is shown in Fig. 5 for the same cases and raw dataof Fig. 4. Large (small) values of ξ IPR corresponds to small(large) values of w and a small (strong) pair forma-tion probability. Here minimum (maximum) values areat p + = ( π , π ) ( p + = (0,0)) as for the effective bandwidthof Fig. 1 (see figure caption for the numerical minimum,maximum and ergodic values). The boundary structureof the NN-model and the circle-structure of the HTC-case are also clearly visible.We have also computed the long time average of thepair formation probability for the HTC-model at largersystem size N =
256 and the special cases of either p + x = p + y or p + y = w for the HTC-model, N = p + = p + x = p + y and both interaction values U = ν = (1 − cos( p + /2))/2. Bothcurves clearly confirm some of the observations of thephase diagrams, i.e. strongest pair formation probabilityat ν = p + x , y = π ) with a somewhat larger maximum ✵✵(cid:0)✁✂✵(cid:0)✂✵(cid:0)✄✂✶ Figure 5
Phase diagram of the inverse participation ratio ξ IPR in the plane of pair momentum p + = ( p + x , p + y ) and computedfrom the same states, data and cases as in Fig. 4. The maxi-mum values corresponding to the red region close to the bot-tom left corner p + = (0,0) are ξ IPR = (top left), ξ IPR = (bottom left), ξ IPR = (top right) and ξ IPR = (bottom right). The minimum values at the top right corner p + = ( π , π ) are ξ IPR = (top left), ξ IPR = (bottom left), ξ IPR = (top right) and ξ IPR = (bottom right). For com-parison the ergodic value is ξ IPR,erg = = and thevalue for the totally symmetrized and localized initial state is ξ IPR,init = . range for U = U = ν = p + x , y =
0) or ν = p + x , y = π ) but still clearly above the ergodic limit forall cases. The precise numerical maximum values of w at ν = w for the HTC model, N =
256 andboth interaction values U = p + y = ν = (1 − cos( p + /2))/2 with p + = p + x . The curve for U = ν ≈ ± p + ≈ π ± π /8) correspond-ing to green-circle with radius r g ≈ π visible in thephase diagram. For U = ν ≈ ± U = w at ν = U = ν = Copyright line will be provided by the publisher . M. Frahm and D. L. Shepelyansky: Coulomb electron pairing in a tight-binding model of La-based cuprate superconductors ✶(cid:0)✲✁✶(cid:0)✲✂✶(cid:0)✵ (cid:0) (cid:0)✄☎ (cid:0)✄✆ (cid:0)✄✝ (cid:0)✄✞ ✶◆✟✠✡☛☞ ✌✰✟✌✰✍✟✌✰✎✇ ✏✑ ♥❂✒✶✓✔✕✖✒✗✘✴☎✙✙✴☎❯❂☎❯❂(cid:0)✄✚❡✛✜✕✢✣✔ Figure 6
Dependence of the electron pair formation proba-bility w on ν = (1 − cos( p + /2))/2 for p + = p + x = p + y andthe HTC-model at U = and N = . w is computedfrom the same long time average as in Fig. 4. The maximumvalue at ν = ν max = is w = ( w = ) for U = ( U = ). See Fig. 5 of [12] for the corresponding fig-ure for the NN-model. For the NN-model the maximum value at ν = ν max = is exactly w = for both interaction values. ✶(cid:0)✲✁✶(cid:0)✲✂✶(cid:0)✵ (cid:0) (cid:0)✄☎ (cid:0)✄✆ (cid:0)✄✝ (cid:0)✄✞ ✶◆✟✠✡☛☞ ✌✰✟✌✰✍☞ ✌✰✎✟❂✇ ✏✑ ♥✒✓✶✔✕✖✗✓✘✙✴☎✚✚✴☎❯✒☎❯✒(cid:0)✄✛❡✜✢✖✣✤✕ Figure 7
Dependence of the electron pair formation probabil-ity w on ν = (1 − cos( p + /2))/2 for p + = p + x , p + y = andthe HTC-model at U = and N = . w is computedfrom the same long time average as in Fig. 4. The value at ν = is w = ( w = ) for U = ( U = ). Figures S1 and S2 of the supplementary material aresimilar to Figs. 6 and 7 respectively but for the inverseparticipation ratio ξ IPR . We also computed a more precise time evolution ofthe pair formation probability w for the larger systemsize N =
384 and certain specific cases p + x = p + y ∈ {0, 2 π /3, π } and p + y = p + x ∈ {7 π /8, π }. The re-sults together with the full space results using the Trot-ter formula approximation at N =
128 are shown in Fig. 8for U = w starts decay-ing from its initial value w = t / ∆ t > t / ∆ t > sometimes with some temporal quasi-periodic fluctua-tions. In most cases the saturation values at U = U = p + y = p + x = π /8 where both saturation values are some-what comparable. In particular, at U = p + y = p + x = π /8 is significantly larger than thevalue for p + y = p + x = π while at U = U = U = p + values, is quitelow if compared to the case p + = U = p + y = p + x = π where the curve is a t ≈ ∆ t evenbelow the ergodic value and saturates later at a value onlyslightly above the ergodic value.Motivated by the observation of the green-circle atradius r g ≈ π in the phase diagram for U = p + x = π /8, p + y = N =
384 and both interaction values U = t / ∆ t = . The first ob-servation is that the diffusive spreading in x -directionis strongly suppressed if compared to the y -directionwhich is expected since p + x is rather close to π while p + y = U = t / ∆ t = , despitea smaller value of w = w = U =
2, has a larger spatial extension of ∼
30 lat-tice sites compared to ∼
12 lattice sites for U =
2. This inrough qualitative agreement with the values ξ IPR = U = ξ IPR =
268 (for U = ξ IPR and the visible Copyright line will be provided by the publishernnalen der Physik, July 27, 2020 ✶(cid:0)✲✁✶(cid:0)✲✂✶(cid:0)✲✄✶(cid:0)✵✶(cid:0)✲✂ ✶(cid:0)✵ ✶(cid:0)✂ ✶(cid:0)✹ ✶(cid:0)✻❯☎(cid:0)✆✝✇ ✞✟ t✠❉t◆✡☛☞✌✍ ✎✰✏✡✎✰✑✡❂◆✡☛☞✌✍ ✎✰✏✡✎✰✑✡✒♣✴☛◆✡☛☞✌✍ ✎✰✏✡✎✰✑✡♣◆✡☛☞✌✍ ✎✰✏✡♣✍ ✎✰✑✡❂◆✡☛☞✌✍ ✎✰✏✡✓♣✴☞✍ ✎✰✑✡❂◆✡✔✒☞✍ ✕✖✗✘✘✙✖✙✖❡✗✚✛✜✍ ◆✡✔✒☞✙✖❡✗✚✛✜✍ ◆✡☛☞✌✢✣✤✥✢✣✤✦✢✣✤✧✢✣★✢✣✤✦ ✢✣★ ✢✣✦ ✢✣✩ ✢✣✪✫✬✭✮ ✯✱ ✳✷✸✳✺✼✽✾✿❀ ❁❃❄✼❁❃❅✼❆✺✼✽✾✿❀ ❁❃❄✼❁❃❅✼❇❈❊✽✺✼✽✾✿❀ ❁❃❄✼❁❃❅✼❈✺✼✽✾✿❀ ❁❃❄✼❈❀ ❁❃❅✼❆✺✼✽✾✿❀ ❁❃❄✼❋❈❊✾❀ ❁❃❅✼❆✺✼●❇✾❀ ❍■❏❑❑▲■▲■▼❏❖P◗❀ ✺✼●❇✾▲■▼❏❖P◗❀ ✺✼✽✾✿
Figure 8
Time dependence of the pair formation probability w for U = (top panel) and U = (bottom panel) anddifferent cases of the exact time evolution in certain p + = ( p + x , p + y ) sectors at N = and the full space Trotter for-mula time evolution at N = . The dashed lines correspondto the ergodic values (21/ N ) = for N = (greydashed) and (21/ N ) = for N = (black dashed). spatial extension in Fig. 9 (for this reason with consider w to be a more suitable quantity than ξ IPR to describethe pair formation probability).
The discussion of the phase diagram given in Fig. 4 hasshown that the pair formation probability is maximal atthe point p + = ( π , π ). However, the surrounding region tothis point is quite small if compared to the green-circlewhere we have a somewhat more modest pair formationprobability. In terms of available values of p + the latterregion is possibly more important. In order to analyzethis point in a more quantitative way, we assume a simple ✵✵(cid:0)✁✂✵(cid:0)✂✵(cid:0)✄✂✶ Figure 9
Color plot of wavefunction amplitude | ¯ ψ ( p + , ∆ r ) | inblock representation in ∆ r = ( ∆ x , ∆ y ) plane obtained from the2D quantum time evolution for the HTC lattice with N = and the sector p + x = π /8 , p + y = . All panels show azoomed region ≤ ∆ x , ∆ y < . Left (right) panels corre-spond to t = ∆ t ( t = ∆ t ; with ∆ t = B = + U ) )and top (bottom) panels correspond to interaction strength U = ( U = ). Related videos are available at [15]. model where both electrons have the same momentum p + /2 (i.e. ∆ p =
0) and where the available states of thistype are filled from smallest to largest energies. We sub-divide these states, ordered in energy, in slices of equalnumber ( ∼ w for each slice which is equivalent to theaverage of w at lines of constant energy. In Fig.10, weshow the dependence of this average on the effective 2D-filling factor ν D which is the weight of slices below a cer-tain energy.For the NN-model we observe a strong peak at ν D = ν D = p + = ( π , π ) and rather strong (top or right)boundary contributions visible in the left panels of Fig. 4.For the HTC-model at U = U = ν D ≈ r g ≈ ν D ≈ p + = ( π , π ). For this partic-ular case, the green circle has a stronger global contribu-tion to the pair formation probability than the maximumregion at p + = ( π , π ). Copyright line will be provided by the publisher . M. Frahm and D. L. Shepelyansky: Coulomb electron pairing in a tight-binding model of La-based cuprate superconductors ✵✵(cid:0)✁✵(cid:0)✂✵(cid:0)✄ ✵ ✵(cid:0)✂ ✵(cid:0)☎ ✵(cid:0)✆ ✵(cid:0)✝ ✁◆◆✇ ✶✞ ♥✷✟❡✠✡☛☞✌✍✎ ◆✏✁✑✂❯✏✂❯✏✵(cid:0)✒✵✵(cid:0)✁✵(cid:0)✂✵(cid:0)✄ ✵ ✵(cid:0)✂ ✵(cid:0)☎ ✵(cid:0)✆ ✵(cid:0)✝ ✁❍✓✔✇ ✶✞ ♥✷✟❡✠✡☛☞✌✍✎ ◆✏✁✑✂❯✏✂❯✏✵(cid:0)✒ Figure 10
Dependence of the electron pair formation prob-ability w on the effective 2D filling factor ν D for the NN-lattice (top) and the HTC-lattice (bottom). The values of w have been obtained from the data of Fig. 4 (for N = ) byan average along lines of constant electron pair energy E c at momenta p = p = p + /2 with p + x , p + y ∈ [0,2 π ] . Low-est (largest) energy corresponds to ν D = ( ν D = ). Thedata points shown correspond to an effective histogram withbin width ∆ ν D ≈ . The red (blue) curve corresponds tothe interaction value U = ( U = ) and the grey dashed linecorresponds to the ergodic value (21/192) = . In our studies we analyzed the electron pair formation ina tight-binding model of La-based cuprate superconduc-tors induced by Coulomb repulsion. Our analytical andnumerical results show that even a repulsive Coulomb in-teraction can form two electron pairs with a high prob-ability. Such pairs have a compact size and propagatethrough the whole system. We expect that such pairs maycontribute to the emergence of superconductivity in La-based cuprates. Of course, our analysis only considers two electronsand in a real system at finite electron density thereis a Fermi sea which can modify electron interactions.However, we expect that electrons significantly belowthe Fermi energy will only create a mean-field potentialwhich will not significantly affect interacting electronswith energies in the vicinity of the Fermi energy. A de-tailed investigation of effects of finite electron densityon the Coulomb pair formation represents an importanttask for future studies.
This work has been partially supported through the grantNANOX N o ANR-17-EURE-0009 in the framework of theProgramme Investissements d’Avenir (project MTDINA).This work was granted access to the HPC resources ofCALMIP (Toulouse) under the allocation 2020-P0110.
A Appendix
In this appendix we present the derivation of the blockHamiltonian (11) and a more detailed discussion aboutits discrete symmetries. In order to simplify the notations,we will use here the full set A ′ = A ∪ ( − A ) of neigh-bor vectors (in the full and not only half plane) for thesummation over the vectors a which allows to reduce thenumber of terms in the following expressions. The TIPHamiltonian (5) can then be written in a more explicitform as: H = − X r , r X a ∈ A ′ t a ¡ | r , r 〉〈 r + a , r | + | r , r 〉〈 r , r − a | ¢ + X r , r ¯ U ( r − r ) | r , r 〉〈 r , r | (13)where for convenience we have written “ r − a ” insteadof “ r + a ” (in the second term of the first line) since for a ∈ A ′ also − a ∈ A ′ . Furthermore, the terms with shiftsof a in the left side have been absorbed by the increasedset A ′ (with respect to A used in (1)) combined with asubsequent shift of the summation index r or r and ex-ploiting the periodic boundary conditions.Applying (13) to a block basis state (10) we find that: H | p + , ∆ r 〉 = − N X r X a ∈ A ′ t a ³ | r − a , r + ∆ r 〉 e i p + · ( r + ∆ r /2) +| r , r + ∆ r + a 〉 e i p + · ( r + ∆ r /2) ´ (14) + ¯ U ( ∆ r ) | p + , ∆ r 〉 . Copyright line will be provided by the publishernnalen der Physik, July 27, 2020
Using the shift r → r + a in the r -sum of the first lineof this expression we obtain: H | p + , ∆ r 〉 = − N X r X a ∈ A ′ t a ³ | r , r + ∆ r + a 〉 e i p + · ( r + a + ∆ r /2) +| r , r + ∆ r + a 〉 e i p + · ( r + ∆ r /2) ´ (15) + ¯ U ( ∆ r ) | p + , ∆ r 〉 which can be rewritten as: H | p + , ∆ r 〉 = − N X r X a ∈ A ′ t a | r , r + ∆ r + a 〉× e i p + · [ r + ( ∆ r + a )/2] ³ e i p + · a /2 + e − i p + · a /2 ´| {z } p + · a /2) + ¯ U ( ∆ r ) | p + , ∆ r 〉 (16) = − X a ∈ A ′ t a cos( p + · a /2) | p + , ∆ r + a 〉 (17) + ¯ U ( ∆ r ) | p + , ∆ r 〉 .The last expression provides exactly the effective blockHamiltonian (11) if we replace the sum over a ∈ A ′ bya sum over a ∈ A with two contributions “ + a ” and “ − a ”and applying for the latter contribution a subsequentshift ∆ r → ∆ r + a in the ∆ r sum. However, there is one ad-ditional complication if ∆ r + a = ( ∆ x + a x , ∆ y + a y ) in (17)leaves the initial square of ∆ x , ∆ y ∈ {0,... N − N to (from) ∆ x + a x and/or ∆ y + a y which provides according to (10) the factor e ± i p + x N /2 = e ± i π l + x = ( − l + x (for ∆ x and similarly for ∆ y ) resultingin either periodic or anti-periodic boundary conditionsin x - ( y -)direction depending on the parity of the integerindex l + x ( l + y ).We close this appendix with a short discussion aboutthe discrete reflection symmetries of the block Hamilto-nian (11) and the possibility to reduce its effective matrixsize N due to such symmetries. For the NN-model, asalready discussed in detail in [12], there are at least twosymmetries with respect to ∆ x → N − ∆ x (reflection at the ∆ y -axis) or ∆ y → N − ∆ y (reflection at the ∆ x -axis) and incase if p + x = p + y there is a third symmetry with respectto ∆ x ↔ ∆ y (reflection at the ∆ x - ∆ y diagonal) which al-lows for an effective matrix size of roughly either N /4 or N /8 (if p + x = p + y ).However, for a more general lattice, such as the HTC-model, or more generally in presence of at least oneneighbor vector a = ( a x , a y ) with both a x a y a = (1,1)) the number of symmetries is reduced.For the most generic case with p + x p + y , p + x p + y ∆ x → N − ∆ x and ∆ y → N − ∆ y which allows for a re-duction of the effective matrix size to ≈ N /2. In thiscase the factors cos( p + · a /2) = cos[( p + x a x + p + y a y )/2] ap-pearing in the effective hopping amplitudes are not mod-ified because the replacement a → − a due the symmetrytransformation only changes the global sign inside thecosine argument. However, this is no longer true if weapply for example the transformation ∆ x → N − ∆ x with-out modifying ∆ y which is equivalent to the replacementof ( a x , a y ) → ( − a x , a y ) of the neighbor vectors. Thereforea single reflection at the ∆ y (or ∆ x ) axis modifies thehopping amplitude (if both a x a y p + x p + y
0) and (11) is (in general) not invariantwith respect to such transformations. However, if either p + x = p + y = ≈ N /4. Also if p + x = p + y ∆ x - ∆ y di-agonal) leading also to an effective matrix size of ≈ N /4.Finally, for the special case p + x = p + y =
0, we have eventhree symmetries (as in the NN-Model for p + x = p + y )with effective matrix size of ≈ N /8. Key words. tight-binding model, interactions, Coulomb elec-tron pairs, cuprates, high-Tc superconductivity
References [1] K. A. Müller, J. G. Bednorz, Z. Phys. B: Condens. Matter , 64, 189.[2] E. Dagotto, Rev. Mod. Phys. , 66, 763.[3] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida,Z. Zaanen, Nature , 518, 179.[4] C. Proust, L. Taillefer, Annu. Rev. Condens. MatterPhys. , 10, 409.[5] P. W. Anderson, Science , 235, 1196.[6] R. Photopoulos, R. Fresard, Ann. Phys. (Berlin) ,1900177.[7] V. J. Emery, Phys. Rev. Lett. , 58, 2794.[8] V. J. Emery, G. Reiter, Phys. Rev. B , 38, 4547.[9] C. M. Varma, Solid State Commun. , 62, 681.[10] Y. B. Gaididei, V. M. Loktev, Phys. Status Solidi ,147, 307.[11] R. S. Markiewicz, S. Sahrakorpi, M. Lindroos, H. Lin,A. Bansil, Phys. Rev. B , 72, 054519.[12] K. M. Frahm, D. L. Shepelyansky, Phys. Rev. Research , 2, 023354.[13] K. M. Frahm, D. L. Shepelyansky, Eur. Phys. J. B ,89, 8.[14] See Supplemental Material at http:XXXX that con-tains extra-figures and videos supporting the mainconclusions and results.
Copyright line will be provided by the publisher nnalen der Physik, July 27, 2020 ✶(cid:0)✵✶(cid:0)✁✶(cid:0)✷✶(cid:0)✸✶(cid:0)✹✶(cid:0)✺ (cid:0) (cid:0)✂✄ (cid:0)✂☎ (cid:0)✂✆ (cid:0)✂✝ ✶◆✞✟✠✡☛ ☞✰✞☞✰✌✞☞✰✍① ■✎✏ ♥❂✑✶✒✓✔✕✑✖✗✴✄✘✘✴✄❯❂✄❯❂(cid:0)✂✙❡✚✛✔✜✢✓ Figure S1
As Fig. 6 but for the inverse participation ratio ξ IPR . ✶(cid:0)✵✶(cid:0)✁✶(cid:0)✷✶(cid:0)✸✶(cid:0)✹✶(cid:0)✺ (cid:0) (cid:0)✂✄ (cid:0)✂☎ (cid:0)✂✆ (cid:0)✂✝ ✶◆✞✟✠✡☛ ☞✰✞☞✰✌☛ ☞✰✍✞❂① ■✎✏ ♥✑✒✶✓✔✕✖✒✗✘✴✄✙✙✴✄❯✑✄❯✑(cid:0)✂✚❡✛✜✕✢✣✔ Figure S2
As Fig. 7 but for the inverse participation ratio ξ IPR . Supplementary Material forCoulomb electron pairing in a tight-binding model ofLa-based cuprate superconductors by K. M. Frahm and D. L. Shepelyansky.Here, we present additional material for the main partof the article.Figure S1 presents data for the inverse participationratio for the case of Fig. 6.Figure S2 presents data for the inverse participationratio for the case of Fig. 7. Two video files for the time evolution obtained bythe Trotter formula approximation corresponding to theparameters of Fig. 2 and Fig. 3 are presented in files videofig2.avi for the density ρ X X ( x , x ) defined inEq. (8) and in videofig3.avi for the density ρ rel ( ∆ x , ∆ y )defined in Eq. (9) (here N = U = t = l j ∆ t (25 valuesper second of video) with integer l =