A Renormalization group approach for highly anisotropic 2D Fermion systems:
application to coupled Hubbard chains
S. Moukouri
Department of Physics and Michigan Center for Theoretical Physics
University of Michigan, 2477 Randall Laboratory, Ann Arbor MI 48109
I apply a two-step density-matrix renormalization group method to the anisotropic two-
dimensional Hubbard model. As a prelude to this study, I compare the numerical results to the
exact one for the tight-binding model. I find a ground-state energy which agrees with the exact
value up to four digits for systems as large as 24 × 25. I then apply the method to the interact-
ing case. I find that for strong Hubbard interaction, the ground-state is dominated by magnetic
correlations. These correlations are robust even in the presence of strong frustration. Interchain
pair tunneling is negligible in the singlet and triplet channels and it is not enhanced by frustration.
For weak Hubbard couplings, interchain non-local singlet pair tunneling is enhanced and magnetic
correlations are strongly reduced. This suggests a possible superconductive ground state.
I. INTRODUCTION
Quasi-one dimensional organic1 and inorganic2 mate-
rials have been the object of an important theoretical in-
terest for the last three decades. The essential features of
their phase diagram may be captured by the anisotropic
Hubbard model (AHM),
H = −t‖
i,l,σ
i,l,σci+1,l,σ + h.c.) + U
ni,l,↑ni,l,↓
i,l,σ
ni,l,σni+1,l,σ − µ
i,l,σ
ni,l,σ
i,l,σ
i,l,σci,l+1,σ + h.c.). (1)
or a more general Hubbard-like model including longer
range Coulomb interactions. The indices i and l label
the sites and the chains respectively. For these highly
anisotropic materials, t⊥ ≪ t‖. Over the years, the
AHM has remained a formidable challenge to condensed-
matter theorists. Some important insights on this model
or its low energy version, the g-ology model, have been
obtained through the work of Bourbonnais and Caron3,4
and others. They used a perturbative renormalization
group approach to analyze the crossover from 1D to 2D
at low temperatures. More recently, Biermann et al.5 ap-
plied the chain dynamical mean-field approach to study
the crossover from Luttinger liquid to Fermi liquid in this
model. Despite this important progress, crucial informa-
tion such as the ground-state phase diagram, or most no-
tably, whether the AHM displays superconductivity, are
still unknown. So far it has remained beyond the reach
of numerical methods such as the exact diagonalization
(ED) or the quantum Monte Carlo (QMC) methods. ED
cannot exceed lattices of about 4 × 5. It is likely to re-
main so for many years unless there is a breakthrough
in quantum computations. The QMC method is plagued
by the minus sign problem and will not be helpful at low
temperatures. The small value of t⊥ implies that, in or-
der to see the 2D behavior, it will be necessary to reach
lower temperatures than those usually studied for the
isotropic 2D Hubbard model. Hence, even in the absence
of the minus sign problem, in order to work in this low
temperature regime, the QMC algorithm requires spe-
cial stabilization schemes which lead to prohibitive cpu
time.6
II. TWO-STEP DMRG
I have shown in Ref. 7 that this class of anisotropic
models may be studied using a two-step density-matrix
renormalization group (TSDMRG) method. The TS-
DMRG method is a perturbative approach in which the
standard 1D DMRG is applied twice. In the first step,
the usual 1D DMRG method9 is applied to find a set
of low lying eigenvalues ǫn and eigenfunctions |φn〉 of a
single chain. In the second step, the 2D Hamiltonian is
then projected onto the basis constructed from the tensor
product of the |φn〉’s. This projection yields an effective
one-dimensional Hamiltonian for the 2D lattice,
E‖[n]|Φ‖[n]〉〈Φ‖[n]| − t⊥
i,l,σ
i,l,σ c̃i,l+1,σ + h.c.)(2)
where E‖[n] is the sum of eigenvalues of the different
chains, E‖[n] =
l ǫnl ; |Φ‖[n]〉 are the corresponding
eigenstates, |Φ‖[n]〉 = |φn1〉|φn2〉...|φnL〉; c̃
i,l,σ, c̃i,l,σ, and
ñi,l,σ are the renormalized matrix elements in the single
chain basis. They are given by
i,l,σ)
nl,ml = (−1)ni〈φnl |c
i,l,σ|φml〉, (3)
(c̃i,l,σ)
nl,ml = (−1)ni〈φnl |ci,l,σ|φml〉, (4)
(ñi,l,σ)
nl,ml = 〈φnl |ni,l,σ|φml〉, (5)
where ni represents the total number of fermions from
sites 1 to i− 1. For each chain, operators for all the sites
are stored in a single matrix
http://arxiv.org/abs/0704.1618v1
l,σ = (c̃
1,l,σ, ..., c̃
L,l,σ), (6)
c̃l,σ = (c̃1,l,σ, ..., c̃L,l,σ), (7)
ñl,σ = (ñ1,l,σ, ..., ñL,l,σ). (8)
Since the in-chain degrees of freedom have been inte-
grated out, the interchain couplings are between the
block matrix operators in Eq.( 6, 7) which depend only
on the chain index l. In this matrix notation, the effec-
tive Hamiltonian is one-dimensional and it is also studied
by the DMRG method. The only difference compared to
a normal 1D situation is that the local operators are now
ms2 ×ms2 matrices, where ms2 is the number of states
kept during the second step.
The two-step method has previously been applied
to anisotropic two-dimensional Heisenberg models.7 In
Ref. 8, it was applied to the t− J model but due to the
absence of an exact result in certain limits, it was tested
against ED results on small ladders only. A systematic
analysis of its performance on a fermionic model on 2D
lattices of various size has not been done. In this paper,
as a prelude to the study of the AHM, I will apply the
TSDMRG to the anisotropic tight-binding model on a
2D lattice, i.e., model (1) with U = V = 0. I perform
a comparison with the exact result of the tight-binding
model. I was able to obtain agreement for the ground-
state energies on the order of 10−4 for lattices of up to
24 × 25. I then discuss how these calculations may be
extended to the interacting case, before presenting the
U 6= 0 results.
III. WARM UP: THE TIGHT-BINDING MODEL
The tight-binding Hamiltonian is diagonal in the mo-
mentum space, the single particle energies are,
ǫk = −2t‖coskx − 2t⊥cosky − µ, (9)
with k = (kx, ky), kx = nxπ/(Lx+1) and ky = nyπ/(Ly+
1) for open boundary conditions (OBC); Lx, Ly are re-
spectively the linear dimensions of the lattice in the par-
allel and transverse directions. The ground-state energy
of an N electron system is obtained by filling the low-
est states up to the Fermi level, E[0](N) =
However in real space, this problem is not trivial and it
constitutes, for any real space method such as the TS-
DMRG, a test having the same level of difficulty as the
case with U 6= 0. This is because the term involving U is
diagonal in real space and the challenge of diagonalizing
the AHM arises from the hopping term.
I will study the tight-binding model at quarter filling,
N/LxLy = 1/2, the nominal density of the organic con-
ductors known as the Bechgaard salts. Systems of up to
Lx×Ly = L×(L+1) = 24×25 will be studied. During the
first step, I keep enough states (ms1 is a few hundred) so
0 20 40 60 80 100
FIG. 1: Low-lying states of the 1D tight-binding model (full
line) and of the 1D Heisenberg spin chain (dotted line) for
L = 16 and ms2 = 96.
that the truncation error ρ1 is less than 10
−6. I target the
lowest state in each charge-spin sectorsNx±2, Nx±1, Nx
and Sz ± 1, Sz ± 2, Nx is the number of electrons within
the chain. It is fixed such that Nx/Lx = 1/2. There is a
total of 22 charge-spin states targeted at each iteration.
For the tight-binding model, the chains remain discon-
nected if t⊥ < ǫ0(Nx + 1) − ǫ0(Nx) or t⊥ < ǫ0(Nx) −
ǫ0(Nx − 1), where Nx is the number of electons on sin-
gle chain. In order to observe transverse motion, it is
necessary that at least t⊥ >∼ ǫ0(Nx + 1) − ǫ0(Nx) and
t⊥ >∼ ǫ0(Nx)− ǫ0(Nx − 1). These two conditions are sat-
isfied only if µ is appropietly chosen. The values listed in
Table (I) corresponds to µ = (ǫ0(Nx+1)−ǫ0(Nx−1))/2.
This treshold varies with L. I give in Table (I) the val-
ues of t⊥ chosen for different lattice sizes. In principle,
for the TSDMRG to be accurate, it is necessary that
∆ǫ = ǫnc − ǫ0, where ǫnc is the cut-off, be such that
∆ǫ/t⊥ ≫ 1. But in practice, I find that I can achieve
accuracy up to the fourth digit even if ∆ǫ/t⊥ ≈ 5 using
the finite system method. Five sweeps were necessary
to reach convergence. Note that this conclusion is some-
what different from my earlier estimate of ∆ǫ/t⊥ ≈ 10
for spin systems.8 This is because in Ref. 8, I used the
infinite system method during the second step.
The ultimate success of the TSDMRG depends on the
density of the low-lying states in the 1D model. For
fixed ms2 and L, it is, for instance, easier to reach
larger ∆ǫ/J⊥ in the anisotropic spin one-half Heisen-
berg model, studied in Ref. 7, than ∆ǫ/t⊥ for the tight-
binding model as shown in Fig.1. For L = 16, ms2 = 96,
and J⊥ = t⊥ = 0.15, I find that ∆ǫ/J⊥ ≈ 10, while
∆ǫ/t⊥ ≈ 5. Hence, the TSDMRG method will be more
accurate for a spin model than for the tight-binding
model. Using the infinite system method during the
second step on the anisotropic Heisenberg model with
J⊥ = 0.1, I can now reach an agreement of about 10
with the stochastic QMC method.
Two possible sources of error can contribute to reduce
the accuracy in the TSDMRG with respect to the conven-
8× 9 16× 17 24× 25
t⊥ 0.28 0.15 0.1
µ -1.2660 -1.3411 -1.3657
∆ǫ/t⊥ 6.42 5.40 5.78
TABLE I: Transverse hopping and chemical potential used in
the simulations for different lattice sizes
ms2 8× 3 16× 3 24× 3
64 -0.241524 -0.211929 0.204040
Exact -0.241524 -0.211931 0.204049
TABLE II: Ground-state energies of three-leg ladders.
tional DMRG. They are the truncation of the superblock
from 4×ms1 states to onlyms2 states and the use of three
blocks instead of four during the second step. In Table
(II) I analyze the impact of the reduction of the number of
states to ms2 for three-leg ladders. The choice of three-
leg ladders is motivated by the fact that at this point,
the TSDMRG is equivalent to the exact diagonalization
of three reduced superblocks. It can be seen that as far as
t⊥ >∼ ǫ0(Nx+1)−ǫ0(Nx) and t⊥ >∼ ǫ0(Nx)−ǫ0(Nx−1), the
TSDMRG at this point is as accurate as the 1D DMRG.
Note that the accuracy remains nearly the same irrespec-
tive of L as far as the ratio ∆ǫ/t⊥ remains nearly con-
stant. Since ∆ǫ decreases when L increases, t⊥ must be
decreased in order to keep the same level of accuracy for
fixed ms2. In principle, following this prescription, much
larger systems may be studied. ∆ǫ/t⊥ does not have to
be very large, in this case it is about 5, to obtain very
good agreement with the exact result.
The second source of error is related to the fact that
the effective single site during the second step is now a
chain having ms2 states, I am thus forced to use three
blocks instead of four to reduce the computational bur-
den. In Table (III), it can be seen that this results in a
reduction in accuracy of about two orders of magnitude
with respect to those of three leg-ladders. These results
are nevertheless very good given the relatively modest
computer power involved. All calculations were done on
a workstation.
The DMRG is less accurate when three blocks are used
instead of four. This can be understood by applying the
following view on the formation of the reduced density
matrix. The construction of the reduced density matrix
ms2 8× 9 16× 17 24× 25
64 -0.24761 -0.21401 0.20504
100 -0.24819 -0.21414 0.20509
120 -0.24832 -0.21419
Exact -0.24857 -0.21432 0.20519
TABLE III: Ground-state energies for different lattice sizes;
a single state was targeted in the second step.
may be regarded as a linear mapping uΨ : F
∗ → E,
where E is the system, F is the environment and, F∗
is the dual space of F. Using the decomposition of the
superblock wave function Ψ[0] =
i ⊗ φRi , with φLi ∈
E and ΦRi ∈ F, for any φ∗ ∈ F ∗,
〈φ∗|φRi 〉φLi . (10)
Let |k〉, k = 1, ...dimE and |l〉, l = 1, ...dimF be the
basis of E and F respectively. Then, |l〉 has a dual basis
〈l∗| such that 〈l∗|l〉 = δl,l∗ . The matrix elements of uΨ
in this basis are just the coordinates of the superblock
wave function Φ[0]k,l . The rank r of this mapping, which
is also the rank of the reduced density matrix is always
smaller or equal to the smallest dimension of E or F,
r < Min(dimE, dimF). Hence, if ms2 states are kept in
the two external blocks, the number of non-zero eigenval-
ues of ρ cannot be larger than ms2. Consequently, some
states which have non-zero eigenvalues in the normal four
block configuration will be missing. A possible cure to
this problem is to target additional low-lying states above
Ψ[0](N). The weight of these states in ρ must be small
so that their role is simply to add the missing states not
to be described accurately themselves. A larger weight
on these additional states would lead to the reduction
of the accuracy for a fixed ms2. In table (IV), I show
the improved energies when, besides the ground state,
I target the lowest states of the spin sectors Sz = −1
and Sz = +1 with N electrons. The weights were re-
spectively 0.995, 0.0025, and 0.0025 for the three states.
This lowers E[0](N) in all cases, but the gain does not
appear to be spectacular. But I do not know whether
this is due to my choice of perturbation of ρ or whether
even the algorithm with four blocks would not yield bet-
ter E[0](N). If the lowest sectors with N + 1 and N − 1
electrons which have Sz = ±0.5 are projected instead, I
find that the results are similar to those with Sz = ±1
sectors, there are possibly many ways to add the missing
states. A more systematic approach to this problem has
recently been suggested.10 It is based on using a local
perturbation to build a correction to the density matrix
from the site at the edge of the system. Here, such a
perturbation would be ∆ρ = αc
l ρcl, where α is a con-
stant, α ≈ 10−3 − 10−2, and c†l , cl are the creation and
annihilation operators of the chain at the edge of the sys-
tem. This type of perturbation resulted in an accuracy
gain of more than an order of magnitude in the case of a
spin chain.10 The three block method was found to be on
par with the four block method. It will be interesting to
see in a future study how this type of local perturbation
performs within the TSDMRG.
To conclude this section, as a first step to the investi-
gation of interacting electron models, I have shown that
the TSDMRG can successfuly be applied to the tight-
binding model. The agreement with the exact result is
very good and can be improved since the computational
power involved in this study was modest. The extension
ms2 8× 9 16× 17
64 -0.24803 -0.21401
100 -0.24828 -0.21417
Exact -0.24857 -0.21432
TABLE IV: Ground-state energies for different lattice sizes;
three states were targeted in the second step: the ground state
itself and the lowest states of Sz = 0 and Sz = 1 sectors.
0 1 2 3 4
FIG. 2: Width ∆ǫ for the low-lying states of the 1D Hubbard
chain as function of U for L = 16 and ms2 = 128.
to the AHM with U 6= 0 is straightforward. There is no
additional change in the algorithm since the term involv-
ing U is local and thus treated during the 1D part of the
TSDMRG. The role of U is to reduce ∆ǫ as shown in
Fig.2. For fixed L and ms2, ∆ǫ decreases linearly with
increasing U . For L = 16 and ms2 = 128, I anticipate
that for U <∼ 3 the interacting system results will be on
the same level or better than those of the non-interacting
case with ms2 = 100 for the same value of L.
IV. GROUND-STATE PROPERTIES OF
COUPLED HUBBARD CHAINS
I now proceed to the study of U 6= 0. One of the
main motivations for such a study is the possibility to
gain insight into the mechanism of superconductivity in
quasi 1D systems. The mechanism of superconductivity
in the quasi 1D organic materials Bechgaard and Fabre
salts, is still an open issue.13 Since these materials are 1D
above a crossover temperature Tx ≈ t⊥/π, it is broadly
accepted that the starting point for the the understand-
ing of their low T behavior should be pure 1D physics.
The occurence of the low T ordered phases is driven by
the interchain hopping t⊥. Two main hypotheses have
been suggested concerning superconductivity. The first
hypothesis (see a recent review in Ref. 13) relies on a
more conventional physics: t⊥ drives the system to a 2D
electron gas which is an anisotropic Fermi liquid which
becomes superconductive through a conventional BCS
0 2 4 6 8
-0.004
-0.003
-0.002
-0.001
FIG. 3: Transverse Green’s function G(y) for td = 0 (circles),
td = 0.1 (squares).
mechanism. However, it has been argued11 that given
the smallness of t⊥, the resulting electron-phonon cou-
pling would not be enough to account for the observed
Tc. The second hypothesis, which has gained strength
over the years given the absence of a clear phonon sig-
nature, is that the pairing mechanism originates from an
exchange of spin fluctuation.11
Interest in this issue was recently revived by the NMR
Knight shift experimental finding that the symmetry of
the Cooper pairs is triplet12 in (TMTSF )2(PF )6. No
shift was found in the magnetic susceptibility at the tran-
sition for measurement made under a magnetic field of
about 1.4 Tesla. A triplet pairing scenario was subse-
quently supported by the persistence of superconduc-
tivity under fields far exceeding the Pauli breaking-pair
limit19. However there is no simple explanation of this
scenario. Triplet pairing would be unfavorable in a BCS
like scenario for which a singlet s-wave is most likely.
Triplet pairing is also less likely in the spin fluctuation
mechanism for which a singlet d-wave is predicted by an-
lytical RG13 or by perturbative approaches17. It has be
argued that these difficulties in both mechanisms can be
circumvented. In the BCS case, the association of AFM
fluctuations with an open Fermi surface to the electron-
phonon mechanism may lead to a triplet pairing18. In the
spin fluctuation case, the addition of interchain Coulomb
interactions may favor a triplet f-wave in lieu of the sin-
glet d-wave13,17. The more exotic Fulde-Ferrel-Larkin-
Ovchinnikov phase can also been invoked to account for
the large paramagnetic limit. However, the Knight shift
result which was thought to bring a conclusion to this
long standing issue has only revived the old controversy.
The conclusion of this experiment itself has been recently
challenged. In Ref. 15, it was pointed out that the ob-
servation of triplet superconductivity claimed in Ref. 12
could be a spurious effect due to the lack of thermaliza-
tion of the samples. A recent Knight shift experimenent
performed at lower fields reveals a decrease in the spin
susceptibility. This is consistent with singlet pairing.16
0 2 4 6 8
-0.002
-0.001
0.001
FIG. 4: Transverse spin-spin correlation C(y) for td = 0 (cir-
cles), td = 0.1 (squares).
The 1D interacting electron gas is now fairly well
understood.3 There is no phase with long range order.
There are essentially four regions in the phase diagram,
characterized by the dominant correlations i.e., SDW,
charge density wave (CDW), singlet superconductivity
(SS) and triplet superconductivity (TS). The essential
question is whether the interchain hopping will simply
freeze the dominant 1D fluctuation into long-range order
(LRO) or create new 2D physics. The estimated values
of U and V for the Bechgaard salts suggest that they
are in the SDW region in their 1D regime. This suggests
that superconductivity in these materials is a 2D phe-
nomenon. Interchain pair tunneling was suggested soon
after the discovery of superconductivity in an organic
compound.14 Emery argued instead that a mechanism
similar to the Kohn-Luttinger mechanism might be re-
sponsible for superconductivity in the organic materials.
When t⊥ is turned on, pairing can arise from exchange of
short-range SDW fluctuations. The reason is that the os-
cillating SDW susceptibility atQ = (2kF , k⊥) would have
an attractive region if k⊥ 6= 0. In particular if k⊥ = π
as I found, then the interaction would be attractive be-
tween particles in neighboring chains. In this study, I
will restrain myself to the study of interchain pair tun-
neling. I was unable to compute correlation functions of
pairs in which each electron belongs to a different chain.
The reason is that in the DMRG method, for the correla-
tion functions to be accurate, at least two different blocks
should be involved. This means that for pair correlation
for which each electron of the pair is on a different chain,
at least four blocks are needed. However, the introduc-
tion of four blocks in the second step of the TSDMRG
leads to a prohibitive CPU time.
With the hope of frustrating an SDW ordering which
is usually expected, I will add an extra terms to model
(1) . These are the diagonal interchain hopping,
Hd = −td
i,l,σ
i,l,σci+1,l+1,σ +
0 2 4 6 8
-0.002
-0.001
0.000
0.001
0 1 2 3 4 5 6 7 8
-0.0001
0.0000
0.0001
FIG. 5: Transverse local singlet correlation SS(y) for td = 0
(circles), td = 0.1 (squares).
h.c) + (c
i+1,l,σci,l−1,σ + h.c), (11)
and the next-nearrest neighbor interchain hopping,
H ′⊥ = −t′⊥
i,l,σ
i,l,σci,l+2,σ + h.c).
I will also add the interchain Coulomb interaction,
HV = V⊥
i,l,σ
ni,l,σ.ni,l+1,σ (12)
I set t⊥ = 0.2, ms1 = 256, ms2 = 128, and (L × (L +
1) = 16×17. A second set of calculations with t⊥ = 0.15,
same values of ms1 and ms2, and (L× (L+1) = 24× 25
lead to the same conclusions. Therefore, they will not be
shown here. In order to analyze the physics induced by
the transverse couplings, I compute the following inter-
chain correlations: the transverse single-particle Green’s
function, shown in Fig.3,
G(y) = 〈cL/2,L/2+yc†L/2,L/2+1〉, (13)
the transverse spin-spin correlation function, shown in
Fig.4,
C(y) =
〈SL/2,L/2+ySL/2,L/2+1〉, (14)
the transverse local pairs singlet superconductive corre-
lation, shown in Fig.5,
SS(y) = 〈ΣL/2,L/2+yΣ†L/2,L/2+1〉, (15)
where
Σi,l = ci,l↑ci,l↓, (16)
the transverse triplet superconductive correlation, shown
in Fig.6,
ST (y) = 2〈ΘL/2,L/2+yΘ†L/2,L/2+1〉, (17)
where
Θi,l =
(ci,l↑ci+1,l↓ + ci,l↓ci+1,l↑), (18)
and the transverse non-local singlet pair superconductive
correlation function, shown in Fig.7,
SD(y) = 2〈∆L/2,L/2+y∆†L/2,L/2+1〉, (19)
where
∆i,l =
(ci,l↑ci+1,l↓ − ci,l↓ci+1,l↑). (20)
A. Strong-coupling regime
Let us first consider, the regime U >∼ 4, I choose
for instance U = 4, V = 0.85, µ = 0, and
td = t
⊥ = V⊥ = 0; besides single-particle hop-
ping, t⊥ also generates two-particle hopping both in
the particle-hole and particle-particle channels. These
two-particle correlation functions are roughly given
by the average values t2⊥〈c
i,lσci,l−σc
i,l+j−σci,l+jσ〉 and
t2⊥〈c
i,lσc
i,l−σci,l+jσci,l+j−σ〉 for an on-site pair created
at (i, l) and then destroyed at (i, l + j). It is expected
that the dominant two-particle correlation are SDW with
k⊥ = π. This is seen in Fig.(4-7). The transverse pair-
ing correlations are all found to be small with respect
to C(y). Among the pairing correlations, SS(y) decays
faster then ST (y) and SD(y). These results are consis-
tent with the view that the role of t⊥ is to freeze the
dominant 1D correlations into LRO.
When td 6= 0, it is expected that for a strong enough td,
the magnetic order will vanish because of the frustration
induced by td. A simple argument is that td induces an
AFM exchange between next-nearest neigbhors on chains
l and l + 1 which compete with the AFM exchange be-
tween nearest neigbhors. The hope is that there could be
a region of the phase diagram where superconductivity
could ultimately win either by pair tunneling between
the chains or by the Emery’s mechanism. However, in
Fig.(4-7) it can be seen that, while td slightly reduces
C(y), the dominant correlations are still SDW even for a
strong td/t⊥ = 0.5. SS(y), ST (y) and SD(y) are barely
affected by td. The fact that td does not strongly af-
fect the SDW order can be understood in the light of
recent study of coupled t− J chains8. It was shown that
the frustration strongly suppresses magnetic LRO only
0 1 2 3 4 5 6 7
-0.002
-0.001
0.000
0.001
0 1 2 3 4 5 6 7 8
-0.0002
-0.0001
0.0000
0.0001
FIG. 6: Transverse triplet superconductive correlation ST (y)
for td = 0 (circles), td = 0.1 (squares).
0 1 2 3 4 5 6 7
-0.002
-0.001
0.000
0.001
0 1 2 3 4 5 6 7 8
-0.0001
0.0000
0.0001
0.0002
FIG. 7: Transverse singlet non-local superconductive correla-
tion SD(y) for td = 0 (circles), td = 0.1 (squares).
close to half-filling. For large dopings, two neighboring
spins in a chain do not always points to opposite direc-
tion as the consequence, td does not necessarily frustrate
the magnetic order. This is illustrated in a simple sketch
in Fig.(8). td could even enhance it as seen in the study
of t− J chains. In Fig.3, it can be seen that td enhances
G(y). This enhancement, together with the decrease of
C(y), suggests a possible widening of an eventual Fermi
liquid region at finite T above the ordered phase. When
t⊥ 6= 0, I also found (not shown) that magnetic correla-
tion are not effectively suppressed even when t′⊥ = t⊥/2.
For this value, it would be expected that the ratio of the
effective exchange term generated by t′⊥ to that generated
by t⊥ is about one quarter. In the frustrated J1−J2 spin
chain, a spin gap opens around this ratio. This simple
picture does not seem to work here.
FIG. 8: sketch of the spin texture (arrows) in two consecutive
chains in an SDW. The bold horizontal lines represent the
chains. The full diagonal lines show bonds for which td tends
to increase the SDW order. The diagonal dotted lines show
bonds for which td frustrates the magnetic order.
B. Weak-coupling regime
I now turn in to the regime where U <∼ 4. I set U = 2,
V = 0, µ = −0.9271, t⊥ = 0.2, td = 0, and V⊥ = 0.4,
where V⊥ is the interchain Coulomb interaction between
nearest neighbors. It can be seen in Fig.9 that C(y) is
now strongly reduced with respect to its strong coupling
values. It is already within our numerical error for the
next-nearest neighbor in the transverse direction. This
is an indication that the ground state is probably not an
SDW. It is to be noted that this occurs even in the ab-
sence of td or t
⊥. This seems to be at variance with the
RG analysis which requires t′⊥ to destroy the magnetic
order. A possible explanation of this is that at half-filling
the perfect nesting occurs at the wave vector Q = (π, π)
for the spectrum of equation (9). Away from half-filling
the nesting is no longer perfect this leads to the reduc-
tion of magnetic correlations. The first correction to the
nesting is an effective frustration term which is roughly
t2⊥cos2k⊥. This expression is identical to a term that
could be generated by an explicit frustration t′⊥ = t
The discrepancy between the TSDMRG and the RG re-
sults could be that this nesting deviation is undereval-
uated in the RG analysis. This mechanism cannot be
invoked in the strong coupling regime where band effects
are small.
The suppression of magnetism is concommitant to a
strong enhancement of the singlet pairing correlations as
seen in Fig. 11. Triplet correlations, shown in Fig. 10,
remain very small. However, while it is clear from the
behavior of C(y) that the ground state is non magnetic.
This result strongly suggests that the ground state is a
superconductor in this regime. A finite size analysis is,
however, necessary to conclude whether this persists to
the thermodynamic limit. I cannot rule out the possibil-
ity of a Fermi liquid ground state, which is implied by
0 2 4 6 8
-0.002
-0.001
0.001
FIG. 9: Transverse spin-spin correlation C(y) for U = 4 (cir-
cles), U = 2 and V⊥ = 0.4 (squares).
0 1 2 3 4 5 6 7
-0.001
0.001
0.002
FIG. 10: Transverse triplet superconductive correlation
ST (y) for U = 4 (circles), U = 2 and V⊥ = 0.4 (squares).
strong single particle correlations.
V. CONCLUSION
In this paper, I have presented a TSDMRG study of
the competition between magnetism and superconductiv-
ity in an anisotropic Hubbard model. I have analyzed the
effect of the interchain hopping in the strong and weak
U regimes. In the strong-coupling regime, the results are
consistent with earlier predictions that the role of t⊥ is
to freeze the dominant 1D SDW correlations into a 2D
ordered state. But at variance with analytical predic-
tions, this is only true in the strong U regime. In this
regime, I find that even the introduction of frustration
does not disrupt the SDW order which remain robust
up to large values of the frustration. In the weak cou-
pling regime singlet pair correlations are dominant. The
ground state seems to be a superconductor. This behav-
ior is somewhat in agreement with experiments in the
Bechgaard or Fabre salts. The phase diagram is domi-
0 1 2 3 4 5 6 7
-0.001
0.001
0.002
FIG. 11: Transverse singlet non-local superconductive cor-
relation SD(y) for U = 4 (circles), U = 2 and V⊥ = 0.4
(squares).
nated by magnetism at low pressure (strong U) and by
superconductivity at high pressure (weak U). Because of
experimental relevance, I restricted myself to the com-
petion between magnetism and superconductivity. I did
not analyze CDW correlations. These are likely to be
important given that I applied open boundary conditions
which are known to generate Friedel oscillations20 that
very decay slowly from the boundaries. They may also
genuinely generated by V⊥, leading to a CDW ground
state instead of a superconductor.
Acknowledgments
I am very grateful to C. Bourbonnais for very helpful
exchanges. I wish to thank A.M.-S. Tremblay for helpful
discussions. This work was supported by the NSF Grant
No. DMR-0426775.
1 C. Bourbonnais and D. Jérome in ”Advance in Synthetic
Metals” Eds. P. Bernier, S. Lefrant and G. Bidan (Elsevier,
New York), 206 (1999).
2 J.W. Allen, Sol. St. Comm. 123, 469 (2002).
3 C. Bourbonnais and L.G. Caron, Int. J. Mod. Phys. B 5,
1033 (1991).
4 T. Giamarchi in ”Quantum Physics in One Dimension”,
Clarendon Press Eds, P. 254-269 (2004).
5 S. Biermann, A. Georges, A. Lichtenstein, and T. Gia-
marchi, Phys. Rev. Lett. 87, 276405 (2001).
6 S.R. White, D.J. Scalapino, R.L. Sugar, E.Y. Loh, J.E. Gu-
bernatis, and R.T. Scalettar, Phys. Rev. B 40, 506 (1989).
7 S. Moukouri, Phys. Rev. B 70, 014403 (2004).
8 S. Moukouri, J. Stat. Mech. P02002 (2006)
9 S.R. White, Phys. Rev. Lett. 69, 2863 (1992). Phys. Rev.
B 48, 10 345 (1993).
10 S.R. White, Phys. Rev. B 72, 180403 (2005).
11 V.J. Emery, Synthetic Metals, 13, 21 (1986).
12 I.J. Lee et al., Phys. Rev. Lett. 88, 017004 (2002).
13 N. Dupuis, C. Bourbonnais and J.C. Nickel,
cond-mat/0510544.
14 D. Jerome and H.J. Schulz, Adv. Phys. 31, 299 (1982).
15 D. Jerome, Chem. Rev. 104, 5565 (2004). D. Jerome
and C.R. Pasquier, in Superconductors, edited by A.V.
Narlikar (Springer Verlag, Berlin, 2005).
16 Y. Shinagawa, et al., cond-mat/0701566 (2007).
17 Y. Tanaka and K. Kuroki, Phys. Rev.B 70, 060502 (2004).
18 Mahito Kohmoto and Masatoshi Sato, cond-mat/0001331
(2000).
19 I.J. Lee, M.J. Naughton, P.M. Chaikin, Physica B 294-
295, 413 (2001).
20 S. R. White, Ian Affleck, and D. J. Scalapino, Phys. Rev.
B 65, 165122 (2002).
http://arxiv.org/abs/cond-mat/0510544
http://arxiv.org/abs/cond-mat/0701566
http://arxiv.org/abs/cond-mat/0001331

  I apply a two-step density-matrix renormalization group method to the
anisotropic two-dimensional Hubbard model. As a prelude to this study, I
compare the numerical results to the exact one for the tight-binding model. I
find a ground-state energy which agrees with the exact value up to four digits
for systems as large as $24 \times 25$. I then apply the method to the
interacting case. I find that for strong Hubbard interaction, the ground-state
is dominated by magnetic correlations.
  These correlations are robust even in the presence of strong frustration.
Interchain pair tunneling is negligible in the singlet and triplet channels and
it is not enhanced by frustration. For weak Hubbard couplings, interchain
non-local singlet pair tunneling is enhanced and magnetic correlations are
strongly reduced. This suggests a possible superconductive ground state.

A Renormalization group approach for highly anisotropic 2D Fermion systems:
application to coupled Hubbard chains
S. Moukouri
Department of Physics and Michigan Center for Theoretical Physics
University of Michigan, 2477 Randall Laboratory, Ann Arbor MI 48109
I apply a two-step density-matrix renormalization group method to the anisotropic two-
dimensional Hubbard model. As a prelude to this study, I compare the numerical results to the
exact one for the tight-binding model. I find a ground-state energy which agrees with the exact
value up to four digits for systems as large as 24 × 25. I then apply the method to the interact-
ing case. I find that for strong Hubbard interaction, the ground-state is dominated by magnetic
correlations. These correlations are robust even in the presence of strong frustration. Interchain
pair tunneling is negligible in the singlet and triplet channels and it is not enhanced by frustration.
For weak Hubbard couplings, interchain non-local singlet pair tunneling is enhanced and magnetic
correlations are strongly reduced. This suggests a possible superconductive ground state.
I. INTRODUCTION
Quasi-one dimensional organic1 and inorganic2 mate-
rials have been the object of an important theoretical in-
terest for the last three decades. The essential features of
their phase diagram may be captured by the anisotropic
Hubbard model (AHM),
H = −t‖
i,l,σ
i,l,σci+1,l,σ + h.c.) + U
ni,l,↑ni,l,↓
i,l,σ
ni,l,σni+1,l,σ − µ
i,l,σ
ni,l,σ
i,l,σ
i,l,σci,l+1,σ + h.c.). (1)
or a more general Hubbard-like model including longer
range Coulomb interactions. The indices i and l label
the sites and the chains respectively. For these highly
anisotropic materials, t⊥ ≪ t‖. Over the years, the
AHM has remained a formidable challenge to condensed-
matter theorists. Some important insights on this model
or its low energy version, the g-ology model, have been
obtained through the work of Bourbonnais and Caron3,4
and others. They used a perturbative renormalization
group approach to analyze the crossover from 1D to 2D
at low temperatures. More recently, Biermann et al.5 ap-
plied the chain dynamical mean-field approach to study
the crossover from Luttinger liquid to Fermi liquid in this
model. Despite this important progress, crucial informa-
tion such as the ground-state phase diagram, or most no-
tably, whether the AHM displays superconductivity, are
still unknown. So far it has remained beyond the reach
of numerical methods such as the exact diagonalization
(ED) or the quantum Monte Carlo (QMC) methods. ED
cannot exceed lattices of about 4 × 5. It is likely to re-
main so for many years unless there is a breakthrough
in quantum computations. The QMC method is plagued
by the minus sign problem and will not be helpful at low
temperatures. The small value of t⊥ implies that, in or-
der to see the 2D behavior, it will be necessary to reach
lower temperatures than those usually studied for the
isotropic 2D Hubbard model. Hence, even in the absence
of the minus sign problem, in order to work in this low
temperature regime, the QMC algorithm requires spe-
cial stabilization schemes which lead to prohibitive cpu
time.6
II. TWO-STEP DMRG
I have shown in Ref. 7 that this class of anisotropic
models may be studied using a two-step density-matrix
renormalization group (TSDMRG) method. The TS-
DMRG method is a perturbative approach in which the
standard 1D DMRG is applied twice. In the first step,
the usual 1D DMRG method9 is applied to find a set
of low lying eigenvalues ǫn and eigenfunctions |φn〉 of a
single chain. In the second step, the 2D Hamiltonian is
then projected onto the basis constructed from the tensor
product of the |φn〉’s. This projection yields an effective
one-dimensional Hamiltonian for the 2D lattice,
E‖[n]|Φ‖[n]〉〈Φ‖[n]| − t⊥
i,l,σ
i,l,σ c̃i,l+1,σ + h.c.)(2)
where E‖[n] is the sum of eigenvalues of the different
chains, E‖[n] =
l ǫnl ; |Φ‖[n]〉 are the corresponding
eigenstates, |Φ‖[n]〉 = |φn1〉|φn2〉...|φnL〉; c̃
i,l,σ, c̃i,l,σ, and
ñi,l,σ are the renormalized matrix elements in the single
chain basis. They are given by
i,l,σ)
nl,ml = (−1)ni〈φnl |c
i,l,σ|φml〉, (3)
(c̃i,l,σ)
nl,ml = (−1)ni〈φnl |ci,l,σ|φml〉, (4)
(ñi,l,σ)
nl,ml = 〈φnl |ni,l,σ|φml〉, (5)
where ni represents the total number of fermions from
sites 1 to i− 1. For each chain, operators for all the sites
are stored in a single matrix
http://arxiv.org/abs/0704.1618v1
l,σ = (c̃
1,l,σ, ..., c̃
L,l,σ), (6)
c̃l,σ = (c̃1,l,σ, ..., c̃L,l,σ), (7)
ñl,σ = (ñ1,l,σ, ..., ñL,l,σ). (8)
Since the in-chain degrees of freedom have been inte-
grated out, the interchain couplings are between the
block matrix operators in Eq.( 6, 7) which depend only
on the chain index l. In this matrix notation, the effec-
tive Hamiltonian is one-dimensional and it is also studied
by the DMRG method. The only difference compared to
a normal 1D situation is that the local operators are now
ms2 ×ms2 matrices, where ms2 is the number of states
kept during the second step.
The two-step method has previously been applied
to anisotropic two-dimensional Heisenberg models.7 In
Ref. 8, it was applied to the t− J model but due to the
absence of an exact result in certain limits, it was tested
against ED results on small ladders only. A systematic
analysis of its performance on a fermionic model on 2D
lattices of various size has not been done. In this paper,
as a prelude to the study of the AHM, I will apply the
TSDMRG to the anisotropic tight-binding model on a
2D lattice, i.e., model (1) with U = V = 0. I perform
a comparison with the exact result of the tight-binding
model. I was able to obtain agreement for the ground-
state energies on the order of 10−4 for lattices of up to
24 × 25. I then discuss how these calculations may be
extended to the interacting case, before presenting the
U 6= 0 results.
III. WARM UP: THE TIGHT-BINDING MODEL
The tight-binding Hamiltonian is diagonal in the mo-
mentum space, the single particle energies are,
ǫk = −2t‖coskx − 2t⊥cosky − µ, (9)
with k = (kx, ky), kx = nxπ/(Lx+1) and ky = nyπ/(Ly+
1) for open boundary conditions (OBC); Lx, Ly are re-
spectively the linear dimensions of the lattice in the par-
allel and transverse directions. The ground-state energy
of an N electron system is obtained by filling the low-
est states up to the Fermi level, E[0](N) =
However in real space, this problem is not trivial and it
constitutes, for any real space method such as the TS-
DMRG, a test having the same level of difficulty as the
case with U 6= 0. This is because the term involving U is
diagonal in real space and the challenge of diagonalizing
the AHM arises from the hopping term.
I will study the tight-binding model at quarter filling,
N/LxLy = 1/2, the nominal density of the organic con-
ductors known as the Bechgaard salts. Systems of up to
Lx×Ly = L×(L+1) = 24×25 will be studied. During the
first step, I keep enough states (ms1 is a few hundred) so
0 20 40 60 80 100
FIG. 1: Low-lying states of the 1D tight-binding model (full
line) and of the 1D Heisenberg spin chain (dotted line) for
L = 16 and ms2 = 96.
that the truncation error ρ1 is less than 10
−6. I target the
lowest state in each charge-spin sectorsNx±2, Nx±1, Nx
and Sz ± 1, Sz ± 2, Nx is the number of electrons within
the chain. It is fixed such that Nx/Lx = 1/2. There is a
total of 22 charge-spin states targeted at each iteration.
For the tight-binding model, the chains remain discon-
nected if t⊥ < ǫ0(Nx + 1) − ǫ0(Nx) or t⊥ < ǫ0(Nx) −
ǫ0(Nx − 1), where Nx is the number of electons on sin-
gle chain. In order to observe transverse motion, it is
necessary that at least t⊥ >∼ ǫ0(Nx + 1) − ǫ0(Nx) and
t⊥ >∼ ǫ0(Nx)− ǫ0(Nx − 1). These two conditions are sat-
isfied only if µ is appropietly chosen. The values listed in
Table (I) corresponds to µ = (ǫ0(Nx+1)−ǫ0(Nx−1))/2.
This treshold varies with L. I give in Table (I) the val-
ues of t⊥ chosen for different lattice sizes. In principle,
for the TSDMRG to be accurate, it is necessary that
∆ǫ = ǫnc − ǫ0, where ǫnc is the cut-off, be such that
∆ǫ/t⊥ ≫ 1. But in practice, I find that I can achieve
accuracy up to the fourth digit even if ∆ǫ/t⊥ ≈ 5 using
the finite system method. Five sweeps were necessary
to reach convergence. Note that this conclusion is some-
what different from my earlier estimate of ∆ǫ/t⊥ ≈ 10
for spin systems.8 This is because in Ref. 8, I used the
infinite system method during the second step.
The ultimate success of the TSDMRG depends on the
density of the low-lying states in the 1D model. For
fixed ms2 and L, it is, for instance, easier to reach
larger ∆ǫ/J⊥ in the anisotropic spin one-half Heisen-
berg model, studied in Ref. 7, than ∆ǫ/t⊥ for the tight-
binding model as shown in Fig.1. For L = 16, ms2 = 96,
and J⊥ = t⊥ = 0.15, I find that ∆ǫ/J⊥ ≈ 10, while
∆ǫ/t⊥ ≈ 5. Hence, the TSDMRG method will be more
accurate for a spin model than for the tight-binding
model. Using the infinite system method during the
second step on the anisotropic Heisenberg model with
J⊥ = 0.1, I can now reach an agreement of about 10
with the stochastic QMC method.
Two possible sources of error can contribute to reduce
the accuracy in the TSDMRG with respect to the conven-
8× 9 16× 17 24× 25
t⊥ 0.28 0.15 0.1
µ -1.2660 -1.3411 -1.3657
∆ǫ/t⊥ 6.42 5.40 5.78
TABLE I: Transverse hopping and chemical potential used in
the simulations for different lattice sizes
ms2 8× 3 16× 3 24× 3
64 -0.241524 -0.211929 0.204040
Exact -0.241524 -0.211931 0.204049
TABLE II: Ground-state energies of three-leg ladders.
tional DMRG. They are the truncation of the superblock
from 4×ms1 states to onlyms2 states and the use of three
blocks instead of four during the second step. In Table
(II) I analyze the impact of the reduction of the number of
states to ms2 for three-leg ladders. The choice of three-
leg ladders is motivated by the fact that at this point,
the TSDMRG is equivalent to the exact diagonalization
of three reduced superblocks. It can be seen that as far as
t⊥ >∼ ǫ0(Nx+1)−ǫ0(Nx) and t⊥ >∼ ǫ0(Nx)−ǫ0(Nx−1), the
TSDMRG at this point is as accurate as the 1D DMRG.
Note that the accuracy remains nearly the same irrespec-
tive of L as far as the ratio ∆ǫ/t⊥ remains nearly con-
stant. Since ∆ǫ decreases when L increases, t⊥ must be
decreased in order to keep the same level of accuracy for
fixed ms2. In principle, following this prescription, much
larger systems may be studied. ∆ǫ/t⊥ does not have to
be very large, in this case it is about 5, to obtain very
good agreement with the exact result.
The second source of error is related to the fact that
the effective single site during the second step is now a
chain having ms2 states, I am thus forced to use three
blocks instead of four to reduce the computational bur-
den. In Table (III), it can be seen that this results in a
reduction in accuracy of about two orders of magnitude
with respect to those of three leg-ladders. These results
are nevertheless very good given the relatively modest
computer power involved. All calculations were done on
a workstation.
The DMRG is less accurate when three blocks are used
instead of four. This can be understood by applying the
following view on the formation of the reduced density
matrix. The construction of the reduced density matrix
ms2 8× 9 16× 17 24× 25
64 -0.24761 -0.21401 0.20504
100 -0.24819 -0.21414 0.20509
120 -0.24832 -0.21419
Exact -0.24857 -0.21432 0.20519
TABLE III: Ground-state energies for different lattice sizes;
a single state was targeted in the second step.
may be regarded as a linear mapping uΨ : F
∗ → E,
where E is the system, F is the environment and, F∗
is the dual space of F. Using the decomposition of the
superblock wave function Ψ[0] =
i ⊗ φRi , with φLi ∈
E and ΦRi ∈ F, for any φ∗ ∈ F ∗,
〈φ∗|φRi 〉φLi . (10)
Let |k〉, k = 1, ...dimE and |l〉, l = 1, ...dimF be the
basis of E and F respectively. Then, |l〉 has a dual basis
〈l∗| such that 〈l∗|l〉 = δl,l∗ . The matrix elements of uΨ
in this basis are just the coordinates of the superblock
wave function Φ[0]k,l . The rank r of this mapping, which
is also the rank of the reduced density matrix is always
smaller or equal to the smallest dimension of E or F,
r < Min(dimE, dimF). Hence, if ms2 states are kept in
the two external blocks, the number of non-zero eigenval-
ues of ρ cannot be larger than ms2. Consequently, some
states which have non-zero eigenvalues in the normal four
block configuration will be missing. A possible cure to
this problem is to target additional low-lying states above
Ψ[0](N). The weight of these states in ρ must be small
so that their role is simply to add the missing states not
to be described accurately themselves. A larger weight
on these additional states would lead to the reduction
of the accuracy for a fixed ms2. In table (IV), I show
the improved energies when, besides the ground state,
I target the lowest states of the spin sectors Sz = −1
and Sz = +1 with N electrons. The weights were re-
spectively 0.995, 0.0025, and 0.0025 for the three states.
This lowers E[0](N) in all cases, but the gain does not
appear to be spectacular. But I do not know whether
this is due to my choice of perturbation of ρ or whether
even the algorithm with four blocks would not yield bet-
ter E[0](N). If the lowest sectors with N + 1 and N − 1
electrons which have Sz = ±0.5 are projected instead, I
find that the results are similar to those with Sz = ±1
sectors, there are possibly many ways to add the missing
states. A more systematic approach to this problem has
recently been suggested.10 It is based on using a local
perturbation to build a correction to the density matrix
from the site at the edge of the system. Here, such a
perturbation would be ∆ρ = αc
l ρcl, where α is a con-
stant, α ≈ 10−3 − 10−2, and c†l , cl are the creation and
annihilation operators of the chain at the edge of the sys-
tem. This type of perturbation resulted in an accuracy
gain of more than an order of magnitude in the case of a
spin chain.10 The three block method was found to be on
par with the four block method. It will be interesting to
see in a future study how this type of local perturbation
performs within the TSDMRG.
To conclude this section, as a first step to the investi-
gation of interacting electron models, I have shown that
the TSDMRG can successfuly be applied to the tight-
binding model. The agreement with the exact result is
very good and can be improved since the computational
power involved in this study was modest. The extension
ms2 8× 9 16× 17
64 -0.24803 -0.21401
100 -0.24828 -0.21417
Exact -0.24857 -0.21432
TABLE IV: Ground-state energies for different lattice sizes;
three states were targeted in the second step: the ground state
itself and the lowest states of Sz = 0 and Sz = 1 sectors.
0 1 2 3 4
FIG. 2: Width ∆ǫ for the low-lying states of the 1D Hubbard
chain as function of U for L = 16 and ms2 = 128.
to the AHM with U 6= 0 is straightforward. There is no
additional change in the algorithm since the term involv-
ing U is local and thus treated during the 1D part of the
TSDMRG. The role of U is to reduce ∆ǫ as shown in
Fig.2. For fixed L and ms2, ∆ǫ decreases linearly with
increasing U . For L = 16 and ms2 = 128, I anticipate
that for U <∼ 3 the interacting system results will be on
the same level or better than those of the non-interacting
case with ms2 = 100 for the same value of L.
IV. GROUND-STATE PROPERTIES OF
COUPLED HUBBARD CHAINS
I now proceed to the study of U 6= 0. One of the
main motivations for such a study is the possibility to
gain insight into the mechanism of superconductivity in
quasi 1D systems. The mechanism of superconductivity
in the quasi 1D organic materials Bechgaard and Fabre
salts, is still an open issue.13 Since these materials are 1D
above a crossover temperature Tx ≈ t⊥/π, it is broadly
accepted that the starting point for the the understand-
ing of their low T behavior should be pure 1D physics.
The occurence of the low T ordered phases is driven by
the interchain hopping t⊥. Two main hypotheses have
been suggested concerning superconductivity. The first
hypothesis (see a recent review in Ref. 13) relies on a
more conventional physics: t⊥ drives the system to a 2D
electron gas which is an anisotropic Fermi liquid which
becomes superconductive through a conventional BCS
0 2 4 6 8
-0.004
-0.003
-0.002
-0.001
FIG. 3: Transverse Green’s function G(y) for td = 0 (circles),
td = 0.1 (squares).
mechanism. However, it has been argued11 that given
the smallness of t⊥, the resulting electron-phonon cou-
pling would not be enough to account for the observed
Tc. The second hypothesis, which has gained strength
over the years given the absence of a clear phonon sig-
nature, is that the pairing mechanism originates from an
exchange of spin fluctuation.11
Interest in this issue was recently revived by the NMR
Knight shift experimental finding that the symmetry of
the Cooper pairs is triplet12 in (TMTSF )2(PF )6. No
shift was found in the magnetic susceptibility at the tran-
sition for measurement made under a magnetic field of
about 1.4 Tesla. A triplet pairing scenario was subse-
quently supported by the persistence of superconduc-
tivity under fields far exceeding the Pauli breaking-pair
limit19. However there is no simple explanation of this
scenario. Triplet pairing would be unfavorable in a BCS
like scenario for which a singlet s-wave is most likely.
Triplet pairing is also less likely in the spin fluctuation
mechanism for which a singlet d-wave is predicted by an-
lytical RG13 or by perturbative approaches17. It has be
argued that these difficulties in both mechanisms can be
circumvented. In the BCS case, the association of AFM
fluctuations with an open Fermi surface to the electron-
phonon mechanism may lead to a triplet pairing18. In the
spin fluctuation case, the addition of interchain Coulomb
interactions may favor a triplet f-wave in lieu of the sin-
glet d-wave13,17. The more exotic Fulde-Ferrel-Larkin-
Ovchinnikov phase can also been invoked to account for
the large paramagnetic limit. However, the Knight shift
result which was thought to bring a conclusion to this
long standing issue has only revived the old controversy.
The conclusion of this experiment itself has been recently
challenged. In Ref. 15, it was pointed out that the ob-
servation of triplet superconductivity claimed in Ref. 12
could be a spurious effect due to the lack of thermaliza-
tion of the samples. A recent Knight shift experimenent
performed at lower fields reveals a decrease in the spin
susceptibility. This is consistent with singlet pairing.16
0 2 4 6 8
-0.002
-0.001
0.001
FIG. 4: Transverse spin-spin correlation C(y) for td = 0 (cir-
cles), td = 0.1 (squares).
The 1D interacting electron gas is now fairly well
understood.3 There is no phase with long range order.
There are essentially four regions in the phase diagram,
characterized by the dominant correlations i.e., SDW,
charge density wave (CDW), singlet superconductivity
(SS) and triplet superconductivity (TS). The essential
question is whether the interchain hopping will simply
freeze the dominant 1D fluctuation into long-range order
(LRO) or create new 2D physics. The estimated values
of U and V for the Bechgaard salts suggest that they
are in the SDW region in their 1D regime. This suggests
that superconductivity in these materials is a 2D phe-
nomenon. Interchain pair tunneling was suggested soon
after the discovery of superconductivity in an organic
compound.14 Emery argued instead that a mechanism
similar to the Kohn-Luttinger mechanism might be re-
sponsible for superconductivity in the organic materials.
When t⊥ is turned on, pairing can arise from exchange of
short-range SDW fluctuations. The reason is that the os-
cillating SDW susceptibility atQ = (2kF , k⊥) would have
an attractive region if k⊥ 6= 0. In particular if k⊥ = π
as I found, then the interaction would be attractive be-
tween particles in neighboring chains. In this study, I
will restrain myself to the study of interchain pair tun-
neling. I was unable to compute correlation functions of
pairs in which each electron belongs to a different chain.
The reason is that in the DMRG method, for the correla-
tion functions to be accurate, at least two different blocks
should be involved. This means that for pair correlation
for which each electron of the pair is on a different chain,
at least four blocks are needed. However, the introduc-
tion of four blocks in the second step of the TSDMRG
leads to a prohibitive CPU time.
With the hope of frustrating an SDW ordering which
is usually expected, I will add an extra terms to model
(1) . These are the diagonal interchain hopping,
Hd = −td
i,l,σ
i,l,σci+1,l+1,σ +
0 2 4 6 8
-0.002
-0.001
0.000
0.001
0 1 2 3 4 5 6 7 8
-0.0001
0.0000
0.0001
FIG. 5: Transverse local singlet correlation SS(y) for td = 0
(circles), td = 0.1 (squares).
h.c) + (c
i+1,l,σci,l−1,σ + h.c), (11)
and the next-nearrest neighbor interchain hopping,
H ′⊥ = −t′⊥
i,l,σ
i,l,σci,l+2,σ + h.c).
I will also add the interchain Coulomb interaction,
HV = V⊥
i,l,σ
ni,l,σ.ni,l+1,σ (12)
I set t⊥ = 0.2, ms1 = 256, ms2 = 128, and (L × (L +
1) = 16×17. A second set of calculations with t⊥ = 0.15,
same values of ms1 and ms2, and (L× (L+1) = 24× 25
lead to the same conclusions. Therefore, they will not be
shown here. In order to analyze the physics induced by
the transverse couplings, I compute the following inter-
chain correlations: the transverse single-particle Green’s
function, shown in Fig.3,
G(y) = 〈cL/2,L/2+yc†L/2,L/2+1〉, (13)
the transverse spin-spin correlation function, shown in
Fig.4,
C(y) =
〈SL/2,L/2+ySL/2,L/2+1〉, (14)
the transverse local pairs singlet superconductive corre-
lation, shown in Fig.5,
SS(y) = 〈ΣL/2,L/2+yΣ†L/2,L/2+1〉, (15)
where
Σi,l = ci,l↑ci,l↓, (16)
the transverse triplet superconductive correlation, shown
in Fig.6,
ST (y) = 2〈ΘL/2,L/2+yΘ†L/2,L/2+1〉, (17)
where
Θi,l =
(ci,l↑ci+1,l↓ + ci,l↓ci+1,l↑), (18)
and the transverse non-local singlet pair superconductive
correlation function, shown in Fig.7,
SD(y) = 2〈∆L/2,L/2+y∆†L/2,L/2+1〉, (19)
where
∆i,l =
(ci,l↑ci+1,l↓ − ci,l↓ci+1,l↑). (20)
A. Strong-coupling regime
Let us first consider, the regime U >∼ 4, I choose
for instance U = 4, V = 0.85, µ = 0, and
td = t
⊥ = V⊥ = 0; besides single-particle hop-
ping, t⊥ also generates two-particle hopping both in
the particle-hole and particle-particle channels. These
two-particle correlation functions are roughly given
by the average values t2⊥〈c
i,lσci,l−σc
i,l+j−σci,l+jσ〉 and
t2⊥〈c
i,lσc
i,l−σci,l+jσci,l+j−σ〉 for an on-site pair created
at (i, l) and then destroyed at (i, l + j). It is expected
that the dominant two-particle correlation are SDW with
k⊥ = π. This is seen in Fig.(4-7). The transverse pair-
ing correlations are all found to be small with respect
to C(y). Among the pairing correlations, SS(y) decays
faster then ST (y) and SD(y). These results are consis-
tent with the view that the role of t⊥ is to freeze the
dominant 1D correlations into LRO.
When td 6= 0, it is expected that for a strong enough td,
the magnetic order will vanish because of the frustration
induced by td. A simple argument is that td induces an
AFM exchange between next-nearest neigbhors on chains
l and l + 1 which compete with the AFM exchange be-
tween nearest neigbhors. The hope is that there could be
a region of the phase diagram where superconductivity
could ultimately win either by pair tunneling between
the chains or by the Emery’s mechanism. However, in
Fig.(4-7) it can be seen that, while td slightly reduces
C(y), the dominant correlations are still SDW even for a
strong td/t⊥ = 0.5. SS(y), ST (y) and SD(y) are barely
affected by td. The fact that td does not strongly af-
fect the SDW order can be understood in the light of
recent study of coupled t− J chains8. It was shown that
the frustration strongly suppresses magnetic LRO only
0 1 2 3 4 5 6 7
-0.002
-0.001
0.000
0.001
0 1 2 3 4 5 6 7 8
-0.0002
-0.0001
0.0000
0.0001
FIG. 6: Transverse triplet superconductive correlation ST (y)
for td = 0 (circles), td = 0.1 (squares).
0 1 2 3 4 5 6 7
-0.002
-0.001
0.000
0.001
0 1 2 3 4 5 6 7 8
-0.0001
0.0000
0.0001
0.0002
FIG. 7: Transverse singlet non-local superconductive correla-
tion SD(y) for td = 0 (circles), td = 0.1 (squares).
close to half-filling. For large dopings, two neighboring
spins in a chain do not always points to opposite direc-
tion as the consequence, td does not necessarily frustrate
the magnetic order. This is illustrated in a simple sketch
in Fig.(8). td could even enhance it as seen in the study
of t− J chains. In Fig.3, it can be seen that td enhances
G(y). This enhancement, together with the decrease of
C(y), suggests a possible widening of an eventual Fermi
liquid region at finite T above the ordered phase. When
t⊥ 6= 0, I also found (not shown) that magnetic correla-
tion are not effectively suppressed even when t′⊥ = t⊥/2.
For this value, it would be expected that the ratio of the
effective exchange term generated by t′⊥ to that generated
by t⊥ is about one quarter. In the frustrated J1−J2 spin
chain, a spin gap opens around this ratio. This simple
picture does not seem to work here.
FIG. 8: sketch of the spin texture (arrows) in two consecutive
chains in an SDW. The bold horizontal lines represent the
chains. The full diagonal lines show bonds for which td tends
to increase the SDW order. The diagonal dotted lines show
bonds for which td frustrates the magnetic order.
B. Weak-coupling regime
I now turn in to the regime where U <∼ 4. I set U = 2,
V = 0, µ = −0.9271, t⊥ = 0.2, td = 0, and V⊥ = 0.4,
where V⊥ is the interchain Coulomb interaction between
nearest neighbors. It can be seen in Fig.9 that C(y) is
now strongly reduced with respect to its strong coupling
values. It is already within our numerical error for the
next-nearest neighbor in the transverse direction. This
is an indication that the ground state is probably not an
SDW. It is to be noted that this occurs even in the ab-
sence of td or t
⊥. This seems to be at variance with the
RG analysis which requires t′⊥ to destroy the magnetic
order. A possible explanation of this is that at half-filling
the perfect nesting occurs at the wave vector Q = (π, π)
for the spectrum of equation (9). Away from half-filling
the nesting is no longer perfect this leads to the reduc-
tion of magnetic correlations. The first correction to the
nesting is an effective frustration term which is roughly
t2⊥cos2k⊥. This expression is identical to a term that
could be generated by an explicit frustration t′⊥ = t
The discrepancy between the TSDMRG and the RG re-
sults could be that this nesting deviation is undereval-
uated in the RG analysis. This mechanism cannot be
invoked in the strong coupling regime where band effects
are small.
The suppression of magnetism is concommitant to a
strong enhancement of the singlet pairing correlations as
seen in Fig. 11. Triplet correlations, shown in Fig. 10,
remain very small. However, while it is clear from the
behavior of C(y) that the ground state is non magnetic.
This result strongly suggests that the ground state is a
superconductor in this regime. A finite size analysis is,
however, necessary to conclude whether this persists to
the thermodynamic limit. I cannot rule out the possibil-
ity of a Fermi liquid ground state, which is implied by
0 2 4 6 8
-0.002
-0.001
0.001
FIG. 9: Transverse spin-spin correlation C(y) for U = 4 (cir-
cles), U = 2 and V⊥ = 0.4 (squares).
0 1 2 3 4 5 6 7
-0.001
0.001
0.002
FIG. 10: Transverse triplet superconductive correlation
ST (y) for U = 4 (circles), U = 2 and V⊥ = 0.4 (squares).
strong single particle correlations.
V. CONCLUSION
In this paper, I have presented a TSDMRG study of
the competition between magnetism and superconductiv-
ity in an anisotropic Hubbard model. I have analyzed the
effect of the interchain hopping in the strong and weak
U regimes. In the strong-coupling regime, the results are
consistent with earlier predictions that the role of t⊥ is
to freeze the dominant 1D SDW correlations into a 2D
ordered state. But at variance with analytical predic-
tions, this is only true in the strong U regime. In this
regime, I find that even the introduction of frustration
does not disrupt the SDW order which remain robust
up to large values of the frustration. In the weak cou-
pling regime singlet pair correlations are dominant. The
ground state seems to be a superconductor. This behav-
ior is somewhat in agreement with experiments in the
Bechgaard or Fabre salts. The phase diagram is domi-
0 1 2 3 4 5 6 7
-0.001
0.001
0.002
FIG. 11: Transverse singlet non-local superconductive cor-
relation SD(y) for U = 4 (circles), U = 2 and V⊥ = 0.4
(squares).
nated by magnetism at low pressure (strong U) and by
superconductivity at high pressure (weak U). Because of
experimental relevance, I restricted myself to the com-
petion between magnetism and superconductivity. I did
not analyze CDW correlations. These are likely to be
important given that I applied open boundary conditions
which are known to generate Friedel oscillations20 that
very decay slowly from the boundaries. They may also
genuinely generated by V⊥, leading to a CDW ground
state instead of a superconductor.
Acknowledgments
I am very grateful to C. Bourbonnais for very helpful
exchanges. I wish to thank A.M.-S. Tremblay for helpful
discussions. This work was supported by the NSF Grant
No. DMR-0426775.
1 C. Bourbonnais and D. Jérome in ”Advance in Synthetic
Metals” Eds. P. Bernier, S. Lefrant and G. Bidan (Elsevier,
New York), 206 (1999).
2 J.W. Allen, Sol. St. Comm. 123, 469 (2002).
3 C. Bourbonnais and L.G. Caron, Int. J. Mod. Phys. B 5,
1033 (1991).
4 T. Giamarchi in ”Quantum Physics in One Dimension”,
Clarendon Press Eds, P. 254-269 (2004).
5 S. Biermann, A. Georges, A. Lichtenstein, and T. Gia-
marchi, Phys. Rev. Lett. 87, 276405 (2001).
6 S.R. White, D.J. Scalapino, R.L. Sugar, E.Y. Loh, J.E. Gu-
bernatis, and R.T. Scalettar, Phys. Rev. B 40, 506 (1989).
7 S. Moukouri, Phys. Rev. B 70, 014403 (2004).
8 S. Moukouri, J. Stat. Mech. P02002 (2006)
9 S.R. White, Phys. Rev. Lett. 69, 2863 (1992). Phys. Rev.
B 48, 10 345 (1993).
10 S.R. White, Phys. Rev. B 72, 180403 (2005).
11 V.J. Emery, Synthetic Metals, 13, 21 (1986).
12 I.J. Lee et al., Phys. Rev. Lett. 88, 017004 (2002).
13 N. Dupuis, C. Bourbonnais and J.C. Nickel,
cond-mat/0510544.
14 D. Jerome and H.J. Schulz, Adv. Phys. 31, 299 (1982).
15 D. Jerome, Chem. Rev. 104, 5565 (2004). D. Jerome
and C.R. Pasquier, in Superconductors, edited by A.V.
Narlikar (Springer Verlag, Berlin, 2005).
16 Y. Shinagawa, et al., cond-mat/0701566 (2007).
17 Y. Tanaka and K. Kuroki, Phys. Rev.B 70, 060502 (2004).
18 Mahito Kohmoto and Masatoshi Sato, cond-mat/0001331
(2000).
19 I.J. Lee, M.J. Naughton, P.M. Chaikin, Physica B 294-
295, 413 (2001).
20 S. R. White, Ian Affleck, and D. J. Scalapino, Phys. Rev.
B 65, 165122 (2002).
http://arxiv.org/abs/cond-mat/0510544
http://arxiv.org/abs/cond-mat/0701566
http://arxiv.org/abs/cond-mat/0001331
