# Superconductivity in striped Hubbard Clusters

## 1 Introduction

Short after the discovery of the high-$`T_c`$ superconductors BED86 , the Hubbard model was introduced as a generic description of the CuO-planes on a microscopic level AND87 . According to the Van Hove scenario we use an extension, the tt’-Hubbard model, to shift the Van Hove singularity in the density of states close to the Fermi energy NEW92 . The experimental result of striped domains CHE89 ; TSU98  in the superconducting CuO-planes inside the high-$`T_c`$ materials has inspired this study of striped Hubbard clusters.

In order to understand superconductivity in the high-$`T_c`$ cuprates on a macroscopic level, the high-$`T_c`$ glass model was introduced in 1987 MOR87 ; MOR89/2 .

It was demonstrated, that the high-$`T_c`$ glass model including the tt’-Hubbard model as a microscopic description of the striped superconducting domains is able to explain several properties of the high-$`T_c`$ cuprates MOR98 , e.g., the d-wave symmetry of the superconducting phase WOL93 ; TSU94  or the pseudogap above $`T_c`$ in the density of states DIN96 ; LOE96 . Furthermore this combined high-$`T_c`$ glass and tt’-Hubbard model picture gives an intuitive description of the experimental puzzle, that different samples of the same material and same doping exhibit a nearly constant superconducting transition temperature $`T_c`$, yet the critical current densities vary from sample to sample MOR99 .

Hubbard clusters were already investigated with numerical algorithms for a number of different geometries and dimensions, e.g. the one-dimensional chains and ladders DAG92 ; FAB92 ; BUL96 ; NOA96 ; DAU99 , two-dimensional (2D) squares SAN89 ; MON94 ; MOR94 ; ZHA97/2 ; FET97/3 , and layered square systems BUL92 ; MOR92 ; MOR92/3 . The stripe instability was found theoretically within Hartree-Fock calculations applied to an extended Hubbard model ZAA89 , and was confirmed by a number of subsequent investigations ZAA99 . But up to now it is not known whether the pure Hubbard model exhibits striping, e.g., in the form of a phase separation. In the closely related 2D $`tJ`$ model the occurrence of stripes is discussed controversially HEL95 ; HEL98 ; WHI99 ; WHI99/2 . Here we study striped clusters of the Hubbard model directly and do not examine the occurrence of stripes per se, but only the existence of superconductivity in striped Hubbard clusters.

In a single striped domain we consider the tt’-Hubbard model, which is described in real space by HUB63 ; GUT63 :

$$=_{kin}+_{pot}$$
(1)   
with the kinetic

$$_{kin}=\underset{i,j,\sigma }{}t_{i,j}(c_{i,\sigma }^{}c_{j,\sigma }+c_{j,\sigma }^{}c_{i,\sigma })$$
(2)   
and the potential part

$$_{pot}=U\underset{i}{}n_{i,}n_{i,}$$
(3)   
of the Hamiltonian. We denote the creation operator for an electron with spin $`\sigma `$ at site $`i`$ with $`c_{i,\sigma }^{}`$, the corresponding annihilation operator with $`c_{i,\sigma }`$, and the number operator at site $`i`$ with $`n_{i,\sigma }c_{i,\sigma }^{}c_{i,\sigma }`$. The hopping $`t_{i,j}`$ is only nonzero for nearest neighbors $`i,j`$ ($`t_{i,j}=t`$) and next nearest neighbors ($`t_{i,j}=t^{}`$). Finally $`U`$ is the interaction. Usually we choose $`t^{}<0`$ to shift the Van Hove singularity in the density of states close to the Fermi energy for less than half filled systems ($`n<1`$, where $`n(n_{e,}+n_{e,})/2`$ and $`n_{e,\sigma }`$ is the number of electrons with spin $`\sigma `$). Throughout this paper we set $`n_{e,}=n_{e,}n_e`$ and the energy unit as $`t=1`$. Additionally we apply periodic boundary conditions both in x- and y-direction.

## 2 Hubbard model and superconducting correlations

Following SCA89 ; FRI91  we use the (vertex) correlation function (instead of the largest eigen value of the reduced two particle density matrix) as an indicator of superconductivity, making use of the standard concept of off-diagonal long range order (ODLRO) YAN62 .

We concentrate here on the d$`_{x^2y^2}`$-wave symmetry (abbreviated as d-wave). The full two-particle correlation function is defined for the the d-wave symmetry by

$$C_d(r)=\frac{1}{L}\underset{i,\delta ,\delta ^{}}{}g_\delta g_\delta ^{}c_{i,}^{}c_{i+\delta ,}^{}c_{i+r+\delta ^{},}c_{i+r,}.$$
(4)   
The vertex correlation function $`C_d^V(r)`$ is the two-particle correlation function $`C_d(r)`$ without the contributions of the single-particle correlations of the same symmetry FRI90 . For the d-wave the result is

$$C_d^V(r)=C_d(r)\frac{1}{L}\underset{i,\delta ,\delta ^{}}{}g_\delta g_\delta ^{}C_{}(i,r)C_{}(i+\delta ,i+r+\delta ^{}).$$
(5)   
In equation (5) $`C_\sigma (i,r)c_{i,\sigma }^{}c_{i+r,\sigma }`$ is the single-particle correlation function for spin $`\sigma `$. The phase factors are $`g_\delta `$, $`g_\delta ^{}=\pm 1`$ to model the d-wave symmetry, the number of lattice points is $`L`$ and the sum $`\delta `$ (resp. $`\delta ^{}`$) is over all nearest neighbors.

We averaged the vertex correlation function $`C_d^V(r)`$ only in the large range regime of $`r`$, i. e. for the distances $`|r|>|r_c|`$:

$$\overline{C}_d^{V,P}\frac{1}{L_c}\underset{r,|r|>|r_c|}{}C_d^V(r)$$
(6)   
with the number $`L_c`$ of lattice points with $`|r|>|r_c|`$. The qualitative behavior of our results (concerning the vertex correlation functions) is not influenced by our choice of $`|r_c|`$ as long as we suppress the short range correlations (i. e. $`r_c1.9`$).

Evidence for ODLRO in the d$`_{x^2y^2}`$ channel was already found for the case of the square 2D tt’-Hubbard model MOR94 ; HUS94 ; FET97/3 . We report here the existence of ODLRO in the striped Hubbard model and investigate the influence of shape and interaction strength on the superconducting signal.

Figure 1 shows the d-wave correlation functions (eq. (4) and (5)) as a function of the distance $`|r|`$ between the pair creation and pair annihilation operators of both the vertex and the full correlation function for a striped system. Similar to results for square systems, we obtain huge correlations for the short range part and finite, positive values for $`C_d^V(r)`$ for large distances in the system (inlay of figure 1). Here it becomes obvious that the vertex correlation function is non-negative in the d-wave case.

For a comparison we plot in figure 1 also the correlation function for the extended s-wave (xs-) symmetry for the nearest neighbors. This symmetry obeys the same formulas as equations (4) and (5) only the phase factors $`g_\delta `$ and $`g_\delta ^{}`$ both are set equal to 1. In contrast to the d-wave case, the xs-wave symmetry does not exhibit this long range behavior (inlay of figure 1). This is also in agreement with simulations for the square Hubbard model MOR94 ; HUS96 .

## 3 The PQMC-Method

We calculate the ground state properties of the Hubbard model using the projector quantum Monte Carlo (PQMC) method BLA81 ; KOO82 . In this algorithm the ground state is projected with

$$|\mathrm{\Psi }_0=\frac{1}{𝒩}e^\theta |\mathrm{\Psi }_T$$
(7)   
from a test state $`|\mathrm{\Psi }_T`$ with the projection parameter $`\theta `$ and the normalization factor $`𝒩`$ LIN92 . In order to perform this projection it is necessary to transform the many-particle problem into a single-particle problem. This is done in two steps, first the exponential of the Hamiltonian $``$ is decomposed into two separate parts, $`_{kin}`$ and $`_{pot}`$, using a Trotter Suzuki transformation SUZ76 ; LIN92  and second the interaction term is treated with a discrete Hubbard Stratonovich (HS) transformation, which leads to an effective single-particle problem with additional fluctuating HS fields HIR83 .

We use the second order Trotter Suzuki transformation, which reads as

$$e^{\theta (_{kin}+_{pot})}=\left(e^{\frac{\tau }{2}_{kin}}e^{\tau _{pot}}e^{\frac{\tau }{2}_{kin}}\right)^m+𝒪(\tau ^2),$$
(8)   
where $`m`$ is the number of Trotter slices and $`\tau \frac{\theta }{m}`$. Here a systematic error of order $`𝒪(\tau ^2)`$ enters the calculations for finite $`m`$.

The two parameters $`m`$ and $`\theta `$ influence the correct projection of the ground state $`|\mathrm{\Psi }_0`$ from the test wave function $`|\mathrm{\Psi }_T`$ in the PQMC algorithm FET98/3 ; FET98/4 .

In figure 2 we investigate the dependence of the ground state energy per site $`E_0/L`$ and of the vertex correlation function $`\overline{C}_d^{V,P}`$ on the Trotter parameter $`m`$. Both, $`E_0/L`$ and $`\overline{C}_d^{V,P}`$, level off for large $`m`$, indicating the convergence of the PQMC method. The results resemble similar PQMC simulations for the square 2D-tt’-Hubbard model BOR92 ; ZHA97 ; FET98 ; FET98/4  and the APEX-oxygen model FRI90 . Like in these cases the vertex correlation function is here more sensitive to $`m`$ than the ground state energy per site $`E_0/L`$ (figure 2).

Figure 3 shows the influence of the projection parameter $`\theta `$ on the same observables. The results are again in good agreement with similar simulations for the square Hubbard model FET97/3 ; FET98 ; FET98/4 .

Due to the more rapid convergence of the ground state energy compared to the vertex correlation function the relative changes of $`E_0/L`$ of figures 2 (a) and 3 (a) are significantly different compared to their (b) counterparts showing the vertex correlation function. This is also expressed by the very different scales of the corresponding y-axes.

Figure 2 (b) shows additionally to a $`12\times 4`$ system the results for a twice as large $`24\times 4`$ system. Convergence occurs here at higher values for $`\theta `$, namely $`\theta >16`$. Due to the sign problem (inlay of figure 3 (a)) we were not able to perform simulations for $`\theta >16`$ . Similar effects occur also for PQMC simulations of the square Hubbard model FET98 .

Quantum Monte Carlo simulations are often plagued with the sign problem ASS90 ; LOH90 . The average sign $`sign`$ enters the calculation for the expectation value $`A`$ of an observable $`A`$ by

$$A=\frac{_{\sigma ,\sigma ^{}}w(\sigma ,\sigma ^{})A(\sigma ,\sigma ^{})}{_{\sigma ,\sigma ^{}}w(\sigma ,\sigma ^{})}=\frac{A^+A^{}}{sign}.$$
(9)   
Here $`\sigma `$ and $`\sigma ^{}`$ are configurations of the HS field, $`w(\sigma ,\sigma ^{})`$ is their weight, and $`A(\sigma ,\sigma ^{})`$ is the expectation value of $`A`$ for $`\sigma `$ and $`\sigma ^{}`$ HIR83 . Now $`w(\sigma ,\sigma ^{})`$ can have both positive and negative values, thus when used in a Monte Carlo algorithm for a transition propability one uses the right hand side of equation (9). $`A^+`$ and $`A^{}`$ denote the separate averages of HS-configurations $`\sigma `$ and $`\sigma ^{}`$ with positive resp. negative weights $`w(\sigma ,\sigma ^{})`$. Generally speaking QMC simulations are only meaningful for $`sign`$ close to 1.

The average sign $`sign`$ is known to decrease for increasing system size $`L`$, interaction strength $`U`$ and projection parameter $`\theta `$ (see the inlay of figure 2 resp. 3 and ASS90 ; LOH90 ). But a small average sign leads among others to large statistical errors in the Monte Carlo process and renders the simulation results meaningless.

From the above analysis of the dependence of the ground state energy and the correlation functions on $`m`$ and $`\theta `$ we conclude, that the PQMC simulations are converged for $`U=2`$ when $`\theta 8`$ in smaller systems and $`\theta 16`$ in larger systems. A ratio $`\tau =\frac{\theta }{m}=\frac{1}{8}`$ of the projection parameter and the Trotter parameter was found to be sufficient for a correct decomposition. Due to the sign problem there is only a small range of the parameters, where $`\theta `$ can be chosen sufficiently large, so that the investigation of the vertex correlation function and its long range behavior is possible. This is similar to the case of the square Hubbard model FET97/3 .

For smaller systems we performed also some simulations of the Hubbard model using the stochastic diagonalization RAE92 ; FET97/2 , figure 6. They compare also favorably with their PQMC counterparts, which is an additional indication that the PQMC performs correctly.

## 4 Superconductivity in stripes

We now investigate the dependence of the superconducting properties of the striped Hubbard model on system size, shape and interaction strength.

The geometry has a quite significant effect on the magnitude of the correlation functions. For increasing width $`L_y`$ of the stripes, figure 4, the average long range part of the vertex correlation function $`\overline{C}_d^{V,P}`$ is decreasing significantly. The ratio between $`\overline{C}_d^{V,P}`$ for a rectangular $`12\times 4`$ and a square $`12\times 12`$ system is almost 3. In figure 5 we show $`\overline{C}_d^{V,P}`$ for both, square systems and the rectangular $`12\times L_y`$ systems from figure 4, as a function of the system size $`L=L_x\times L_y`$. Here the rectangular shaped systems always show a higher superconducting signal than square systems of the same size $`L`$.

This is rather surprising when one takes into account that in striped systems, on average, the distances $`|r|`$ between pair creation and pair annihilation operators are larger than in square systems with the same number of sites $`L`$. In our view there are two effects which may increase the superconducting correlations. First the anisotropy of $`L_x`$ and $`L_y`$ which leads to more finite size shells. Finite size shells refer to the energy levels of the free Hubbard clusters ($`U=0`$). It is known, that other ways of introducing additional shells, e.g., anisotropic hopping $`t_xt_y`$ KUR96 ; KUR97 , or additional hopping to next nearest neighbors HUS96  increase the superconducting correlations in the repulsive Hubbard model.

In our view it is a second effect, the squeezing of the system in one dimension, that gives rise to these increased superconducting correlations.

In contrast to the width $`L_y`$, the length $`L_x`$ of the stripes is relatively insensitive to the height of the plateau, figure 6.

Another way to strengthen the superconducting correlations is to increase the (repulsive) Hubbard interaction U, as shown in figure 7. Here we present both common methods for analyzing superconductivity: the full correlation function IMA91 ; SCA89 ; WHI89  and the vertex function SCA89 ; FRI91 ; FET97/2 . The dotted lines in figure 7(a) and (b) indicate the values of the full (vertex) correlation function in the case of no interaction $`U=0`$. Figure 7(a) shows that the full correlation function increases for higher interactions. The vertex correlation function is zero for the case of no interaction ($`U=0`$) and increases also monotonous for increased interaction strength. Due to the sign problem we were not able to perform simulations for an interaction strength $`U>2.5`$, for the system size and filling shown in figure 7.

Thus within the range of parameters accessible by the PQMC method, the superconducting correlations are increasing for an increasing repulsive interaction strength $`U`$. Both full and vertex correlation function show this behavior. We want to note, that the vertex correlation function is much more sensitive to variations of $`U`$ due to the substraction of the background of the single-particle correlation functions (figure 7). Thus the vertex correlation function is the more appropriate observable to analyze superconductivity in small Hubbard clusters. These results are similar to observations made for the BCS reduced Hubbard model FET97/2 .

## 5 Effective interaction in striped Hubbard clusters

Due to the failure of the usual finite size scaling in the square 2D Hubbard model we introduced HUS96 ; FET97/3 ; FET98  an effective model, the BCS-reduced Hubbard model, to compare the superconducting correlations for different system sizes $`L`$. This failure is mainly caused by the underlying shell structure of the free ($`U=0`$) system FET98 ; FET98/3 ; HUS96 .

The BCS-reduced Hubbard model exhibits the same corrections to scaling as the Hubbard model, and has a well chosen interaction term, that produces superconductivity with d-wave symmetry. We calculate for this model the same correlation functions as for the Hubbard model. The effective interaction strength $`J_{eff}`$ is then chosen to give the same values for the correlation functions as the Hubbard model (for details see FET98 ). From this we get a direct estimate of the superconducting interaction strength.

The calculation of an effective interaction for the three band Hubbard model was used to identify the pairing mechanism for d-wave superconductivity in this model. The evidence of d-wave pairing in this case is based on symmetry arguments and exact diagonalization results of small clusters CIN97 ; CIN98 ; CIN98/2 ; CIN99 .

In the momentum space the BCS-reduced Hubbard model is described by the Hamiltonian:

$$^{BCS}=_{kin}^{BCS}+_{int}^{BCS}.$$
(10)   
The kinetic part is again equation (2), only transformed to momentum space, i. e. $`_{kin}^{BCS}=_{kin}`$, and the interaction is defined by (for d-wave interaction)

$$_{int}^{BCS}=\frac{J}{L}\underset{\genfrac{}{}{0pt}{}{k,p}{kp}}{}f_kf_pc_{k,}^{}c_{k,}^{}c_{p,}c_{p,}.$$
(11)   
In equation (11) we use the form factors

$$f_k\mathrm{cos}(k_x)\mathrm{cos}(k_y)$$
(12)   
to model the d-wave symmetry in 2D ($`k(k_x,k_y)`$).

We calculate the ground state of this BCS-reduced Hubbard model with the exact and the stochastic diagonalization RAE92 ; RAE92/2 ; FET97/2 .

In figure 8 and 9 we show the effective interaction $`J_{eff}`$ corresponding to the correlation functions and systems shown in figure 4 and 6. Within the error bounds of the simulations we conclude that $`J_{eff}`$ is nearly constant for various geometries of the system. It is not possible to calculate stochastic errors of the physical observables within the SD. But for smaller system sizes our comparison of SD with exact diagonalization results indicates that for the weak interactions $`J`$ used here the errors in the SD are negligible HUS96 ; FET97/2 . The error bars shown in figures 8 and 9 are therefore calculated using only the statistical errors of the PQMC results and fitting these values to the SD results.

In addition to the above mentioned, one has to take into account, that even so we tried to perform the calculations at a constant filling $`n0.8`$ the constraint of closed shells for PQMC simulations leads to different fillings $`n`$ for each of these system sizes. Furthermore, in the case of figure 9 (and 6 respectively) one has to take into account, that all simulations are performed at a constant $`\theta =8`$. Whereas figure 3 indicates that for large system sizes $`L`$ a higher value of $`\theta `$ would lead to slightly higher values of the vertex correlation function in the PQMC runs and thus to a slightly lower effective interactions $`J_{eff}`$.

From figures 8 and 9 we conclude that within the accuracy of the applied methods, the effective interaction strength $`J_{eff}`$ is equal for both square and striped systems. Furthermore $`J_{eff}`$ is insensitive to the length of the striped systems.

## 6 Summary and Conclusions

Here we performed ground state simulations of the 2D tt’-Hubbard model and the BCS-reduced Hubbard model of striped clusters using PQMC and SD techniques. Together with the exact diagonalization these are the most reliable computational tools for this type of calculations.

We concentrated our investigations on the behavior of rectangular striped systems. In agreement with previous calculations of the square Hubbard model we find that these finite systems show evidence for superconductivity in the $`d_{x^2y^2}`$ channel for repulsive interactions $`U`$. Compared to the squared case these correlations are significantly enhanced, and the superconducting signal is nearly insensitive to the length of these stripes.

Using SD-techniques we were capable of estimating the effective superconducting interaction strength $`J_{eff}`$ of a BCS-reduced Hubbard model with the same symmetry of the superconducting correlation functions. Within the accuracy of our methods both square and striped Hubbard model show approximately the same superconducting interaction strength $`J_{eff}`$.

In conclusion, the striped Hubbard model is a promising candidate for the microscopic description of the superconducting striped domains in the high-$`T_c`$ cuprates. Within the larger framework of the high-$`T_c`$ glass model a combined model is able to explain many puzzling properties of the high-T<sub>c</sub> materials.

## 7 Acknowledgment

We want to thank P.C. Pattnaik, D.M. Newns, C.C. Tsuei, T. Doderer, H. Keller, T. Schneider, J.G. Bednorz, and K.A. Müller for very helpful discussions. The LRZ Munich grants us a generous amount of CPU time on their IBM SP2 parallel computer, which is highly appreciated. Finally we acknowledge the financial support of the UniOpt GmbH, Regensburg.