text:
# Superexchange in the quarter- filled two-leg ladder system NaV2O₅

## I Introduction

Since the phase transition into a spin-gapped phase in $`\alpha ^{}`$-NaV2O$`_5`$was reported  this layered insulating compound has attracted much attention. It was found that the phase transition ($`T_c34`$ K) is accompanied both by the lattice distortion  and a charge ordering  in the vanadium layers. Though the details of the low-temperature phase structure in NaV2O$`_5`$remain controversial, some main features of this structure are obtained from experiment. For instance, the same low-T lattice superstructure was measured and reported by three groups of researchers . The experiment  has also ruled out a hypothetical ’in-line’ charge ordering (CO) in favour of a ’zig-zag’ CO. At the same time mechanism of the spin-gap formation is not understood up to now and formulation of adequate spin models for NaV2O$`_5`$is required.

The electronic and magnetic properties of NaV2O$`_5`$are mainly due to the V ions having in average a mixed valence +4.5. In a particular layer V-ions are arranged to give a sequence of two-leg quater-filled ladders (one electron per rung) coupled to each other in a Trellis lattice. It is expected that the main properties of vanadium layers are dominated by the ladder physics and the theoretical consideration should start with an analysis of this one-dimensional structure unit. The most important evidence in favour of this point of view is provided by the band structure calculations  within the density-functional theory (DFT). These calculations have shown that the overall electronic band structure of NaV2O$`_5`$is determined mainly by the large intra-ladder electronic hopping amplitudes while the inter-ladder hopping is rather weak.

Let us now turn to the problem of spin coupling in NaV2O$`_5`$and consider a separate quarter- filled ladder. Since the strongest electronic hopping transfers $``$ ta between the basic vanadium dxy -orbitals are within the rungs  one has to ascribe a spin–$`1/2`$ not to a vanadium site but to the intra-rung V-O-V bonding molecular orbital. With strong Coulomb repulsion $`U`$ on the dxy -orbital the lowest electronic bonding orbital cannot be doubly occupied and the problem is naturally reduced to a half-filled chain in the limit of strong electron correlations. The excitation energy of the lowest 2-particle singlet rung state from now on plays a role of an effective Hubbard repulsion parameter which is roughly estimated to be 2ta in the charge disordered state . Horsch and Mack were the first  who derived a quantum Heisenberg spin–$`1/2`$ chain model to describe a one dimensional behavior of the ladder systems NaV2O$`_5`$in the charge disordered high-T phase.

The superexchange coupling parameters entering into the effective Heisenberg chain model are expected  to be strongly reduced as the system undergoes the phase transition into the low-T charge ordered state. It was predicted, for instance, that for neighbouring spins in a ladder the exchange integral $`J`$, which is antiferromagnetic in high-T phase, may even change the sign to become a ferromagnetic one in the extreme CO phase with a zig-zag charge distribution on the ladder rungs . In the present paper we develop the superexchange theory and obtain estimates for the spin-spin coupling constants in a single ladder both in the charge disordered and in an ordered phase. That allows us to reexamine some earlier statements and estimates . Since the in-line CO is ruled out by experiment the present consideration is restricted to a zig-zag CO state only.

The nature of the CO phase transition in the real NaV2O$`_5`$compound is not well understood up to now. Nevertheless, at the first stage (Section III) the analytical expressions for the superexchange integrals in this system can be derived and analyzed in a rather general context. Next, to reach a quantitative description one has to rely on a particular realistic mechanism driving the charge ordering. We assume that the intersite Coulomb electron repulsion is responsible for the CO phase transition and the lattice degrees of freedom play a secondary role. The inter-ladder short-range Coulomb interaction is treated in a mean field approximation while intra-ladder charge-charge correlations are considered on a more accurate level. To this end considering at some stage (Section IV) the charge degrees of freedom only we transform the electronic Hamiltonian onto the effective Ising model in the transverse magnetic field . The pseudospin variables in this approach are associated with the charge dynamics and the order-disorder transition is determined by the dimensionless interaction parameter $`\lambda =V_b/2t_a`$. Here $`V_b`$ is the characteristic Coulomb repulsion between electrons on neighbouring rungs. We adopt this approach which allows us to examine quantitatively how the superexchange coupling varies as the ladder system evolves from a disordered to the charge ordered phase. We expect that besides the quantitative estimates this particular analysis reveals general, mechanism-independent, trends in the behaviour of the superexchange spin-spin coupling in NaV2O$`_5`$.

## II Electronic model and the eigenstate problem

We start with the electronic Hamiltomian for a quater-filled ladder in the following general form

$$H=\underset{i}{}\left(H_0^{(i)}+H_{t,}^{(i)}\right)+\underset{<ij>}{}\left(H_{t,}^{(ij)}+H_V^{(ij)}\right)$$
(1)
where $`(i)`$ and $`(ij)`$ refer to a particular rung or to a pair of nearest neighbouring ($`nn`$) rungs, respectively. The intra-rung interactions involved into ($`H_0^{(i)}+H_{t,}^{(i)}`$) are described by the electronic hopping amplitude $`t_a`$, the Coulomb on-site, $`U`$, and inter-site, $`V_a`$, parameters. For the $`nn`$ inter-rung coupling ($`H_{t,}^{(ij)}+H_V^{(ij)}`$), there are two, $`t_b`$ and $`t_d`$, electronic hopping amplitudes and the corresponding Coulomb repulsion is given by $`V_b`$ (see Fig.1). By introducing the creation operator $`d_{iL\sigma }^{}`$ ($`d_{iR\sigma }^{}`$) for a spin- $`\sigma `$ electron on the left (right) V–ion within the i-th rung we specify the particular terms in Eq.(1) as follows:

$`H_0^{(i)}`$ $`=`$ $`\epsilon {\displaystyle \underset{\alpha =L,R}{}}n_i^\alpha +V_an_i^Ln_i^R+U{\displaystyle \underset{\alpha =L,R}{}}n_i^\alpha n_i^\alpha `$ (2)
$`H_{t,}^{(i)}`$ $`=`$ $`t_a{\displaystyle \underset{\sigma }{}}(d_{iL\sigma }^{}d_{iR\sigma }+h.c.)`$ (3)
$`H_{t,}^{(ij)}`$ $`=`$ $`t_b{\displaystyle \underset{\alpha =L,R}{}}{\displaystyle \underset{\sigma }{}}(d_{i\alpha \sigma }^{}d_{j\alpha \sigma }+h.c.)`$ (5)
$`t_d{\displaystyle \underset{\sigma }{}}(d_{iL\sigma }^{}d_{jR\sigma }+d_{iR\sigma }^{}d_{jL\sigma }+h.c.)`$
$`H_V^{(ij)}`$ $`=`$ $`V_b{\displaystyle \underset{\alpha =L,R}{}}n_i^\alpha n_j^\alpha `$ (6)
where $`n_i^\alpha =_\sigma n_{i\sigma }^\alpha =_\sigma d_{i\alpha \sigma }^{}d_{i\alpha \sigma }`$ ($`\alpha =L,R`$) and $`\epsilon `$ is a bare vanadium $`d_{xy}`$-orbital energy in an isolated ladder.

Let us fix now the range of the model parameters. We use the following set of hopping amplitudes : $`t_a`$=0.38 eV, $`t_bt_d`$ =0.085 eV. The proper signs of these amplitudes are adopted in the definition of $`H`$. Importance of additional hopping processes $`t_d`$, which couple $`d_{xy}`$\- orbitals of two V– ions at the opposite ends of $`nn`$ rungs of a ladder (Fig. 1), was reported in Ref.(). Approximate equality $`t_dt_b`$ allows one to explain the pronounced feature of the band structure in NaV2O$`_5`$. Namely, there is almost no dispersion of the antibonding- type molecular orbital states formed by the $`d_{xy}`$– orbitals on a rung. We will consider also another set: $`t_a`$=0.38 eV, $`t_b`$=0.17 eV and $`t_d`$=0, as a more representative one among the sets with $`t_d`$=0 used by other authors . The value of the on- site Coulomb repulsion $`U`$ is fixed from the supercell LDA+U calculations  to be $`U=4`$ eV. The estimates for the magnitudes of intersite Coulomb repulsion $`V_a`$ and $`V_b`$ (Fig. 1) given in the literature are less accurate. In Ref.() for instance, it is claimed that a screening reduces $`V_a`$ and $`V_b`$ parameters to rather small values $`V_a`$, $`V_b`$ $`<`$ 0.4 eV, and hence the inter- rung Coulomb repulsion cannot stabilize CO. To investigate this point we developed a combined analysis by using the results of DFT calculation  and more generic multi-orbital Hubbard model for a ladder in the Trelllis lattice. We found that inter-orbital vanadium d-d electron correlations lead to an enhancement for the effective $`V_b`$ parameter which enters into the effective single-orbital Hamiltonian Eqs.(16). These results will be presented elsewhere . In the present investigation we take an intermediate value for $`V_a`$0.5 meV and allow $`V_b`$ to vary in a wide range up to 1 eV.
It is worth noting that the validity of the perturbation procedure developed below to the second order in the inter-rung hoppings is based mainly on the smallness of these hoppings, $`t_b`$ and $`t_d`$, compared to $`2t_a`$ provided the on-site Coulomb repulsion U is the largest parameter of the theory. This guarantees also that the quarter- filled ladder is in an insulating state . In contrast to the perturbation approach developed in Ref.() in our case there are no formal restrictions on the values of $`V_a`$ and $`V_b`$. Like in Ref.() these parameters could be arbitrary small giving a meaningfull description of the superexchange in a charge disordered state of NaV2O$`_5`$. Surely, $`V_b`$ should be large enough, $`V_b>V_b^{(\text{cr})}`$, to sustain a zig-zag CO state and the critical value $`V_b^{(\text{cr})}`$ depends on the mechanism chosen for the CO phase transition.

In the perturbation procedure the zeroth order Hamiltonian is defined as $`H_0=H_{<ij>}H_{t,}^{(ij)}`$ where the last term is a perturbation mixing electronic states on $`nn`$ rungs. As the necessary prerequisite of the superexchange theory the 1- and 2- particle electronic states on a particular rung have to be found explicitly. That can be done by considering the inter- rung Coulomb interaction, the term $`_{<ij>}H_V^{(ij)}`$ in $`H_0`$, in a mean field (MF) approximation. In this case the inter- rung Coulomb interactions lead to the following shift of the on- site orbital energies $`(\alpha =L,R)`$:

$$\epsilon _{i\alpha }=\epsilon +V_b\left[n_{i+1}^\alpha _g+n_{i1}^\alpha _g\right]$$
(7)
where $`n_{i\pm 1}^\alpha _g`$ are the averages over the lowest singly occupied electronic states of two rungs neighboring to the i-th rung within the same ladder. Having $`_\alpha n_i^\alpha _g=1`$, one obtains $`(\epsilon _{iL}+\epsilon _{iR})/2=`$$`\epsilon +V_b\overline{\epsilon }`$ and

$$\epsilon _{iL}\epsilon _{iR}=\mathrm{\Delta }\epsilon _i=(1)^i\mathrm{\Delta }\epsilon ;\mathrm{\Delta }\epsilon =2V_b\eta $$
(8)
In Eq.(8) we introduced the parameter $`\eta `$ ($`0\eta 1`$) of charge disproportionation $`n_{iR}_gn_{iL}_g`$$`=(1)^i\eta `$ in a zig-zag CO state of a ladder.

After applying the approximation given in Eq.(7) the zeroth order Hamiltonian $`H_0`$ becomes a sum over individual rungs, $`H_0_iH_0^{(i)}`$. From now on each term $`H_0^{(i)}`$ can be diagonalised independently. It is worth emphasizing that in this approach the local intra- rung correlations due to $`U`$ and $`V_a`$ are treated exactly. To specify explicitly a form of $`H_0^{(i)}`$ in a particular $`n`$ particle sector of the intra- rung electronic states we invoke an extra subscript like in $`H_{0,n}^{(i)}`$ below. For instance, by solving the eigenstate problem in the 1– particle sector one obtains

$$H_{0,1}^{(i)}=\underset{\sigma }{}\left[\epsilon _gg_{i\sigma }^{}g_{i\sigma }+\epsilon _ff_{i\sigma }^{}f_{i\sigma }\right]$$
(9)
where the energy $`\epsilon _{g/f}`$ of the lowest/excited molecular orbital is given by

$$\epsilon _{g/f}=\overline{\epsilon }1/2\sqrt{(\mathrm{\Delta }\epsilon )^2+(2t_a)^2}$$
(10)
and

$`g_{i\sigma }^{}|0_i`$ $`=`$ $`\left[u_id_{iR\sigma }^{}+v_id_{iL\sigma }^{}\right]|0_i|g_{i\sigma }`$ (11)
$`f_{i\sigma }^{}|0_i`$ $`=`$ $`\left[v_id_{iR\sigma }^{}+u_id_{iL\sigma }^{}\right]|0_i|f_{i\sigma }`$ (12)
The $`u_i`$ and $`v_i`$ coefficients in Eq.(12) are defined in the following way

$$u_i=\sqrt{\frac{1}{2}\left[1+\frac{\mathrm{\Delta }\epsilon _i}{\sqrt{(\mathrm{\Delta }\epsilon )^2+(2t_a)^2}}\right]};v_i=\sqrt{\frac{1}{2}\left[1\frac{\mathrm{\Delta }\epsilon _i}{\sqrt{(\mathrm{\Delta }\epsilon )^2+(2t_a)^2}}\right]}$$
(13)
and for $`\mathrm{\Delta }\epsilon _i0`$ the values of $`u_i`$ and $`v_i`$ alternate as i runs along the ladder. The complete basis for the 2- particle states on the i-th rung consists of three singlets $`|S_{i,n}>`$ ($`n=1,2,3`$) and triplet states $`|\tau _i^{(m)}>`$ ($`m=0,\pm 1`$). Explicitly, we choose the following definitions for these states (the rung index i is implied):

$`|S_1=1/\sqrt{2}\left(d_R^{}d_L^{}d_R^{}d_L^{}\right)|0`$ (14)
$`|S_{2,3}=1/\sqrt{2}\left(\pm d_R^{}d_R^{}+d_L^{}d_L^{}\right)|0`$ (15)
$`|\tau ^{(0)}=1/\sqrt{2}\left(d_R^{}d_L^{}+d_R^{}d_L^{}\right)|0`$ (16)
$`|\tau ^{(\pm 1)}=d_R^{}d_L^{}|0,d_R^{}d_L^{}|0`$ (17)
One can check that the triplets are degenerate eigenstates of $`H_{0,2}^{(i)}`$, i.e. $`H_{0,2}^{(i)}|\tau _i^{(m)}=\epsilon _\tau |\tau _i^{(m)}`$, with the energy $`\epsilon _\tau `$$`=2\overline{\epsilon }+V_a`$. In the singlet subspace the Hamiltonian $`H_{0,2}^{(i)}`$ takes the following form

$`H_{0.2}^{(i)}=(|S_{i,1},|S_{i,2},|S_{i,3},\left)\right(\begin{array}{ccc}2\overline{\epsilon }+V_a\hfill & 2t_a\hfill & 0\hfill \\ 2t_a\hfill & 2\overline{\epsilon }+U\hfill & \mathrm{\Delta }\epsilon _i\hfill \\ 0\hfill & \mathrm{\Delta }\epsilon _i\hfill & 2\overline{\epsilon }+U\hfill \end{array}\left)\right(\begin{array}{c}S_{i,1}|,\hfill \\ S_{i,2}|,\hfill \\ S_{i,3}|,\hfill \end{array})`$ (24)
By diagonalizing the symmetric 3$`\times `$3 matrix one obtains the eigenstates

$$|S_i^{(\kappa )}=\underset{n=1}{\overset{3}{}}\gamma _n^{(\kappa )}(i)|S_{i,n}$$
(25)
and the corresponding eigenvalues $`\epsilon _{s\kappa }`$ ($`\kappa `$=1,2,3) in the singlet subspace. More precisely, $`\epsilon _{s\kappa }=2\overline{\epsilon }+V_a+(UV_a)e_\kappa `$, where $`e_\kappa (\kappa =1,2,3)`$, are the solutions of a cubic characteristic equation reading

$`e_\kappa \left[(1e_\kappa )^2\zeta _i^2\right]+\mu ^2[1e_\kappa ]=0`$ (26)
where $`\zeta _i=\mathrm{\Delta }\epsilon _i/(UV_a)`$ and $`\mu =2t_a/(UV_a)`$. Note, that according to Eq.(8) the value of $`\zeta _i^2`$ does not depend on $`i`$. Provided Eq.(26) is solved the 3-component vector $`\stackrel{}{\gamma }^{(\kappa )}(i)`$ in Eq.(25) can be found by using the following expressions

$`\gamma _1^{(\kappa )}`$ $`=`$ $`\pm {\displaystyle \frac{|A_\kappa |}{R_\kappa }}`$ (27)
$`\gamma _2^{(\kappa )}`$ $`=`$ $`\pm sign(A_\kappa ){\displaystyle \frac{\mu (1e_\kappa )}{R_\kappa }}`$ (28)
$`\gamma _3^{(\kappa )}(i)`$ $`=`$ $`sign(A_\kappa ){\displaystyle \frac{\mu \zeta _i}{R_\kappa }}`$ (29)
where $`A_\kappa =(1e_\kappa )^2\zeta _i^2`$ and $`R_\kappa =[A_\kappa ^2+\mu ^2`$$`(1e_\kappa )^2+\mu ^2\zeta _i^2]^{1/2}`$. Note, that the upper/or the lower signs in Eq.(29) should be taken simultaneously. Like $`\zeta _i`$, the third component $`\gamma _3^\kappa (i)`$ is staggered along the chain direction. This particular problem, Eqs.(26, 29), can be solved numerically in a wide range of the underlying model parameters provided the order parameter $`\eta `$ is fixed from some self-consistent procedure. Most representative results will be discussed in Sec.V.

For further purposes it is helpful to use the projection operators $`X_i^{q,q}=|q_iq_i|`$ and the transfer operators $`X_i^{q,p}=|q_ip_i|`$ where $`|q_i`$ and $`|p_i`$ are the states defined above on the i– th rung. More explicitely, for instance, one has

$`X_i^{g\sigma ,g\sigma }`$ $`=`$ $`|g_{i\sigma }g_{i\sigma }|,X_i^{0,g\sigma }=|0_ig_{i\sigma }|,`$ (30)
$`X_i^{S\kappa ,S\kappa }`$ $`=`$ $`|S_i^{(\kappa )}S_i^{(\kappa )}|,X_i^{S\kappa ,g\sigma }=|S_i^{(\kappa )}g_{i\sigma }|\text{etc.}`$ (31)
With these notations the approximate zeroth order Hamiltonian $`H_0_iH_0^{(i)}`$ can be written as

$`H_0^{(i)}`$ $`=`$ $`\epsilon _vX_i^{0,0}+\epsilon _g{\displaystyle \underset{\sigma }{}}X_i^{g\sigma ,g\sigma }+\epsilon _f{\displaystyle \underset{\sigma }{}}X_i^{f\sigma ,f\sigma }+`$ (33)
$`+{\displaystyle \underset{\kappa =1}{\overset{3}{}}}\epsilon _{s\kappa }X_i^{S\kappa ,S\kappa }+\epsilon _\tau {\displaystyle \underset{m=\pm 1,0}{}}X_i^{\tau m,\tau m}`$
In Eq.(33) the vacuum rung state $`|0_i`$ with the reference energy $`\epsilon _v`$ is involved for completeness.

At the same time the 3– and 4– particle states are dropped as these states are irrelevant for the present purposes. At the end of this section the following remark is worth to be made. Actually, while solving the intra- rung eigenstate problems we could not avoid the use of a MF treatment of inter- rung Coulomb interactions, the term $`_{<ij>}H_V^{(ij)}`$ in the original Hamiltonian of Eq.(6). Nevertheless it does not mean the complete neglect of the many body effects due to $`H_V^{(ij)}`$ in our approach. Further analysis in Secs. IV,V will show that charge-charge correlations between singly occupied $`nn`$\- rungs are important to obtain the estimates for the superexchange coupling between spins in a ladder.

## III Perturbation procedure and
superexchange on a ladder

The superexchange coupling between two singly occupied rungs is established due to the electron hopping transitions with the virtual intermediate 2-particle states in one of two rungs. These transitions are due to the kinetic part of the Hamiltonian, Ht,∥, whose action has now to be represented in terms of the molecular orbital states $`|q_i`$ defined in the previous section. In a rather general form one may write:

$`H_{t,}={\displaystyle \underset{<ij>}{}}H_{t,}^{(ij)}={\displaystyle \underset{<ij>,\sigma }{}}{\displaystyle \underset{qq^{}}{}}{\displaystyle \underset{pp^{}}{}}[t_{ij,\sigma }(q^{},q|p^{},p)X_i^{q^{},q}X_j^{p^{},p}+h.c.]`$ (34)
where

$`t_{ij,\sigma }(q^{},q|p^{},p)`$ $`=`$ $`t_b{\displaystyle \underset{\alpha }{}}q_i^{}|d_{i\alpha \sigma }^{}|q_ip_j^{}|d_{j\alpha \sigma }|p_j`$ (36)
$`t_d[q_i^{}|d_{iR\sigma }^{}|q_ip_j^{}|d_{jL\sigma }|p_j+(RL)]`$
and ($`RL`$) means the indices interchange. All the necessary amplitudes $`q_i^{}|d_{i\alpha \sigma }^{}|q_i`$ and their conjugates can be now straightforwardly evaluated to give the relevant part of $`H_{t,}`$ in the following form

$`H_{t,}`$ $`=`$ $`{\displaystyle \underset{ij,\sigma }{}}\xi _{ij}\{{\displaystyle \frac{(2\sigma )}{\sqrt{2}}}{\displaystyle \underset{\kappa =1}{\overset{3}{}}}D_{ij}^{(\kappa )}(X_i^{S\kappa ,g\overline{\sigma }}X_j^{0,g\sigma }+h.c.)+`$ (38)
$`+{\displaystyle \underset{\sigma ^{}}{}}{\displaystyle \underset{m=\pm 1,0}{}}B_{ij}C_{\sigma ^{}\sigma }^{(m)}(X_i^{\tau m,g\sigma ^{}}X_j^{0,g\sigma }+h.c.)\}`$
Here $`\xi _{ij}=1`$ if $`(ij)`$ is a $`nn`$ pair and equals zero otherwise; $`(2\sigma )=+1(\sigma =)`$, $`1(\sigma =)`$ and $`\overline{\sigma }=\sigma `$. The factor $`C_{\sigma ^{}\sigma }^{(m)}`$ is defined as

$$C_{\sigma ^{}\sigma }^{(m)}=\delta _{\sigma ^{},\sigma }\delta _{m,2\sigma }+1/\sqrt{2}\delta _{\sigma ^{},\overline{\sigma }}\delta _{m,0}.$$
(39)
The effective transfer amplitudes $`D_{ij}^{(\kappa )}`$ and $`B_{ij}`$ can be presented in the following form:

$`D_{ij}^{(\kappa )}`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{3}{}}}\gamma _n^{(\kappa )}(i)\left[t_ba_{ij,n}+t_d\overline{a}_{ij,n}\right]`$ (40)
$`B_{ij}`$ $`=`$ $`t_bb_{ij}+t_d\overline{b}_{ij}`$ (41)
where $`a_{ij,n}`$, $`\overline{a}_{ij,n}`$, $`b_{ij}`$ and $`\overline{b}_{ij}`$ are bilinear functions of the $`u`$–, $`v`$– coefficients

$`a_{ij,1}=\overline{a}_{ij,2}=u_iv_j+v_iu_j`$ (42)
$`a_{ij,2}=\overline{a}_{ij,1}=u_iu_j+v_iv_j`$ (43)
$`a_{ij,3}=\overline{b}_{ij}=u_iu_j+v_iv_j`$ (44)
$`b_{ij}=\overline{a}_{ij,3}=u_iv_j+v_iu_j`$ (45)
To derive the superexchange spin coupling to the second order in the effective hopping amplitudes $`D_{ij}^{(\kappa )}`$ and $`B_{ij}`$ we use the Schrieffer-Wolff transformation $`\stackrel{~}{H}`$=$`\mathrm{exp}(\widehat{S})H\mathrm{exp}(\widehat{S})`$ with the generator $`\widehat{S}(=\widehat{S}^{})`$ being determined from the condition $`[H_0,\widehat{S}]=H_{t,}`$. Then the second order correction to $`H_0`$ is given by

$$H_{\text{eff}}=1/2[\widehat{S},H_{t,}]$$
(46)
From Eq.(38) we obtain first the form of the generator

$`\widehat{S}`$ $`=`$ $`{\displaystyle \underset{ij,\sigma }{}}\xi _{ij}\{{\displaystyle \frac{(2\sigma )}{\sqrt{2}}}{\displaystyle \underset{\kappa =1}{\overset{3}{}}}{\displaystyle \frac{D_{ij}^{(\kappa )}}{\epsilon _{s\kappa }+\epsilon _v2\epsilon _g}}X_i^{S\kappa ,g\overline{\sigma }}X_j^{0,g\sigma }+`$ (48)
$`+{\displaystyle \underset{\sigma ^{},m}{}}{\displaystyle \frac{B_{ij}C_{\sigma ^{}\sigma }^{(m)}}{\epsilon _\tau +\epsilon _v2\epsilon _g}}X_i^{\tau m,g\sigma ^{}}X_j^{0,g\sigma }h.c.\}`$
One can see, that there are two perturbation channels, $`(s)`$ and $`(\tau )`$, with the singlet and triplet intermediate virtual states, respectively. Finally we obtain from Eq.(46) the superexchange Hamiltonian

$$H_{\text{eff}}=\underset{<ij>}{}\left(J_{ij}^{(s)}+J_{ij}^{(\tau )}\right)\stackrel{}{S}_i\stackrel{}{S}_j$$
(49)
with the exchange constants in the singlet and triplet channel given by

$`J_{ij}^{(s)}`$ $`=`$ $`{\displaystyle \underset{\kappa =1}{\overset{3}{}}}{\displaystyle \frac{[D_{ij}^{(\kappa )}]^2+[D_{ji}^{(\kappa )}]^2}{\epsilon _{s\kappa }+\epsilon _v2\epsilon _g}}`$ (50)
$`J_{ij}^{(\tau )}`$ $`=`$ $`{\displaystyle \frac{B_{ij}^2+B_{ji}^2}{\epsilon _\tau +\epsilon _v2\epsilon _g}}`$ (51)
In Eq.(49) the spin- $`\frac{1}{2}`$ operators $`S_i^\alpha =1/2_{\sigma \sigma ^{}}|g_{i\sigma }\sigma _{\sigma \sigma ^{}}^\alpha g_{i\sigma ^{}}|`$ and $`\sigma _{\sigma \sigma ^{}}^\alpha `$ are the Pauli matrices $`(\alpha =x,y,z)`$; the notation $`<ij>`$ means that each pair of the $`nn`$ rungs is taken only once.

Note that already at this stage the expressions Eq.(50) and Eq.(51) together with the definitions in Eqs. (41), (43), formulate a complete scheme to estimate $`J_{(ij)}^{(s)}`$ and $`J_{(ij)}^{(\tau )}`$ both in a disordered and a zig-zag CO state. That can be done, at least numerically, by solving first the eigenstate problems of Sec.II for given parameters of the model and with the order parameter $`\eta `$ varying in the range $`0\eta 1`$. In this approach the effects of charge ordering upon spin coupling in a ladder can be examined in its main features without specifying a mechanism driving the CO transition in a quarter- filled ladder system. Preliminary analysis shows that the total superexchange coupling $`J_{(ij)}=J_{(ij)}^{(s)}+J_{(ij)}^{(\tau )}`$ is reduced as the ladder develops from a disorderd to a zig-zag CO state. That is in accordance with a predicition in Ref.(). A more satisfactory quantitative description could be done if one relies on a specific mechanism responsible for CO phase transition. In the next section we pursue one possibility by noting that the underlying model, Eqs.(1 \- 6) contains the expected instability without invoking extra degrees of freedom like charge-phonon coupling. Further analysis is based on the pseudospin description of the charge correlations . Before turning to the pseudospin model itself (Sec.IV) we use this representation to present the main results of this section given in Eqs.(41, 43, 50, 51) in a more condensed form.

In the pseudospin representation the molecular ground state of the i–th rung $`|g_{i\sigma }=|g_i|\sigma _i`$ is given by

$$|g_i=\left[u_i|_i+v_i|_i\right];u_i^2+v_i^2=1$$
(52)
By using the pseudo- spin ($`T=\frac{1}{2}`$) operators $`T_i^\alpha (\alpha =x,y,z)`$ we define now several pair correlation functions and expectation values in the ground state of a ladder as follows

$`Q_{ij}^{()}`$ $`=`$ $`T_i^zT_j^z_g,Q_{ij}^{()}=T_i^xT_j^x+T_i^yT_j^y_g`$ (53)
$`R_{ij}^{(x)}`$ $`=`$ $`T_i^x+T_j^x_g,R_{ij}^{(z)}=T_i^zT_j^z_g`$ (54)
$`R_{ij}^{(xz)}`$ $`=`$ $`T_i^xT_j^zT_i^zT_j^x_g`$ (55)
By using the above definitions and Eqs.(41, 43) one may write the numerators in the expressions of Eqs.(50, 51) in the following form:

$`\left[D_{ij}^{(\kappa )}\right]^2+(ij)=4t_b^2\{{\displaystyle \underset{n=1}{\overset{3}{}}}\left[\gamma _n^{(\kappa )}\right]^2𝒦_{ij}^{(n)}+2\gamma _1^{(\kappa )}\gamma _2^{(\kappa )}𝒦_{ij}^{(12)}+`$ (56)
$`+\gamma _1^{(\kappa )}[\gamma _3^{(\kappa )}(i)\gamma _3^{(\kappa )}(j)]𝒦_{ij}^{(13)}+\gamma _2^{(\kappa )}[\gamma _3^{(\kappa )}(i)\gamma _3^{(\kappa )}(j)]𝒦_{ij}^{(23)}\}`$ (57)
$`B_{ij}^2+(ij)=4t_b^2𝒦_{ij}^{(4)}`$ (58)
where

$`𝒦_{ij}^{(1),(2)}=\left(1+\nu ^2\right)\left[{\displaystyle \frac{1}{4}}+Q_{ij}^{()}\right]\left(1\nu ^2\right)Q_{ij}^{()}+2\nu R_{ij}^{(x)}`$ (59)
$`𝒦_{ij}^{(3),(4)}=\left(1+\nu ^2\right)\left[{\displaystyle \frac{1}{4}}Q_{ij}^{()}\right]\pm \left(1\nu ^2\right)Q_{ij}^{()}`$ (60)
$`𝒦_{ij}^{(12)}=2\nu \left[{\displaystyle \frac{1}{4}}+Q_{ij}^{()}\right]+\left(1+\nu ^2\right)R_{ij}^{(x)}`$ (61)
$`𝒦_{ij}^{(13)}=\nu ^2R_{ij}^{(xz)}{\displaystyle \frac{1}{2}}\nu R_{ij}^{(z)}`$ (62)
$`𝒦_{ij}^{(23)}=\nu R_{ij}^{(xz)}{\displaystyle \frac{1}{2}}\nu ^2R_{ij}^{(z)}`$ (63)
and we defined $`\nu =t_d/t_b`$. It is worth noting that the antisymmetric property, $`𝒦_{ij}^{(13)}=𝒦_{ji}^{(13)}`$ and $`𝒦_{ij}^{(23)}`$$`=𝒦_{ji}^{(23)}`$, makes the complete expression of Eq.(57) to be symmetric.

## IV Pseudo- spin model for the ladder

The effective Hamiltonian in Eq.(49) describes spin-spin interactions in the low- energy subspace spanned by the vector charge- spin manifold $`|G_{\{g\},\{\sigma \}}`$$`=_i[|g_i|\sigma _i]`$. Our aim now is to project the electronic model, Eqs.(1-6) onto the charge sector within this manifold. The projection procedure was discussed several times , it results in the following zeroth- order effective Hamiltonian for a ladder

$`H_{\text{eff}}^{(0)}`$ $`=`$ $`2t_a{\displaystyle \underset{i}{}}T_i^x+2V_b{\displaystyle \underset{<ij>}{}}\left[T_i^zT_j^z+1/4\right]`$ (65)
$`\mathrm{\Delta }\epsilon ^0{\displaystyle \underset{i}{}}(1)^iT_i^z`$
In this one- dimensional (1D) model the effective longitudinal field $`\mathrm{\Delta }\epsilon ^0`$ expresses in a condensed form effects of Coulomb interactions due to higher dimensionality of the real NaV2O$`_5`$compound (see Appendix).

We proceed further in the spirit of the coupled quantum spin chain approach  which involves a MF treatment of the inter-chain coupling and a more sophisticated treatment of the intra-chain, i.e. intra-ladder, interactions. The intra- ladder problem enters into the expressions for the superexchange constants, Eqs.(50, 51) and Eqs.(54-61), in the form of pair correlation functions for the $`nn`$ pseudospins. Having the final aim of calculating these constants we rely on properties of the exact solution for the 1D Ising model in a transverse field.

The term $`\mathrm{\Delta }\epsilon ^0`$ in Eq.(65) breaks explicitly the Z(2) symmetry of the 1D transverse Ising model which results in a non-zero longitudinal magnetisation $`T_i^z_g=(1)^i\eta /2`$ if the coupling constant $`\lambda =V_b/2t_a`$ exceeds some critical value, $`\lambda >\lambda _c`$. One may write

$$\eta =\left[1\left(\frac{\lambda _c}{\lambda }\right)^2\right]^\beta $$
(66)
and choose for the exponent $`\beta =1/8`$ since the presence of a weak longitudinal field does not change the universality class of the 1D Ising model in the transverse field. Note, that for this exactly solvable model $`\lambda _c=1`$ (Ref. ) and we keep this value in the subsequent numerical analysis. The open problem that still remains are unknown dependencies of $`\beta `$ and $`\lambda _c`$ on details of the inter- ladder Coulomb interaction if one goes beyond a mean field treatment of these interaction in the Trellis lattice.

Let us now turn to a discussion of the superexchange constants, Eqs.(50, 51). We estimate the inter- rung pseudospin correlation functions entering into Eq.(54) at two different levels of sophistication. The simplest one is to decouple the pair correlation functions as follows:

$`Q_{ij}^{()}`$ $``$ $`T_i^z^2=\left(\eta /2\right)^2`$ (67)
$`Q_{ij}^{()}`$ $``$ $`T_i^x^2=1/4\left(1\eta ^2\right)`$ (68)
$`R_{ij}^{(xz)}`$ $``$ $`T_i^zT_j^zT_i^x=(1)^{i+1}{\displaystyle \frac{\eta }{2}}\sqrt{1\eta ^2}`$ (69)
In Eq.(68) we take into account that $`i`$ and $`j`$ are refered to $`nn`$ rungs and, hence, $`T_j^z=T_i^z`$ while $`T_j^x=\sqrt{1\eta ^2}/2`$ does not depend ont the rung index $`i`$. Within this first approximation the superexchange constants will be calculated in the next section. The only signatures of the exact solution in Eqs.(66, 68) are due to the choice of the values for $`\beta `$ and $`\lambda _c`$, while the inter- rung fluctuations are ignored. To gain an impression how these fluctuations may change the superexchange constants we adopt some more properties of the exact solution for the 1D transverse Ising model. For this purpose let us consider  the pair correlation functions $`(j=i\pm 1)`$:

$`T_i^zT_j^z_g=k_z;T_i^xT_j^x_g=m_x^2+k_x;T_i^yT_j^y_g=k_y`$ (70)
In particular, one has

$`k_z={\displaystyle \frac{1}{4\pi }}{\displaystyle _0^\pi }𝑑q{\displaystyle \frac{\lambda +\mathrm{cos}(q)}{\mathrm{\Lambda }_q}};m_x={\displaystyle \frac{1}{2\pi }}{\displaystyle _0^\pi }𝑑q{\displaystyle \frac{1+\lambda \mathrm{cos}(q)}{\mathrm{\Lambda }_q}}=T_i^x_g`$ (71)
where $`\mathrm{\Lambda }_q=\sqrt{1+\lambda ^2+2\lambda \mathrm{cos}q}`$. The expressions for $`k_{x,y}`$ can be found in Refs.() as well. We use the property $`|k_x|,|k_y|m_x,|k_z|`$ and drop the weak transverse fluctuations, which means k$`{}_{x,y}{}^{}`$ 0. Breaking of the Z(2) symmetry allows us to write also

$$R_{ij}^{(xz)}(1)^{i+1}\eta m_x$$
(72)
We use the Eqs.(66, 70) and (71) together with the approximation in Eq.(72) as the basis of our second approach to calculate the superexchange constants.

## V Analysis of superexchange integrals,
numerical results

The starting point of our analysis is the exact diagonalisation of the single-rung Hamiltonian in the 1- and 2-particle sectors. The excited singlet and triplet states provide the intermediate states for the superexchange mechanism of magnetic coupling between neighbouring singly occupied rungs. From now on we drop the irrelevant lower indices at $`J`$. The excitation energies $`E_{s\kappa }=\epsilon _{s\kappa }+\epsilon _v2\epsilon _g`$ and $`E_\tau =\epsilon _\tau +\epsilon _v2\epsilon _g`$ which enter as the denominators into the perturbative expressions for $`J^{(s)}`$, $`J^{(\tau )}`$, are determined by the eigenvalues due to Eqs.(10) and (26). These eigenvalues depend on the degree of charge order via the intra-rung d- orbital levels shifts $`\mathrm{\Delta }\epsilon `$ which is proportional to the CO parameter $`\eta `$ in Eq.(8). The symmetry breaking source is assumed to be arbitrary small. The excitation energies in the singlet and triplet channels as function of $`\mathrm{\Delta }\epsilon `$ are shown in Fig.(2), $`\epsilon _v=0`$. There is a crossing of two singlets at very high values of $`\mathrm{\Delta }\epsilon _i`$. For realistic values of Vb however this crossing point is not reached.

While solving the eigenvalue problem in the singlet sector the components of the vectors $`\stackrel{}{\gamma }^{(\kappa )}(i)`$, Eq.(29), are calculated as well. These vector components together with pseudospin correlation functions determine the numerator, Eq.(57), in the expression for $`J^{(s)}`$, Eq.(50). The pseudospin pair correlation functions are calculated either in the modified MF approximation, Eq.(68), or by using a more accurate approach based on Eqs.(70)-(72). In the latter case the most pronounced charge fluctuation contributions to the superexchange are involved. This requires evaluation of only two pair correlation functions, namely, $`T_i^zT_{i+1}^z_g=k_z`$ and $`T_i^xT_{i+1}^x_gm_x^2`$ which are given by Eq.(70),(71). Both of them are shown in Fig.(3) as function of $`\lambda ^1`$=$`2t_a/V_b`$. It is seen that the intra- chain fluctuations of the order parameter described by the ’longitudinal’ correlation function kz extend far above the critical value $`\lambda _c^1`$=1 where the CO vanishes. Since in our approximation k$`{}_{x}{}^{}0`$ the ’transverse’ correlation function is reduced to the square of the ’covalency’ of each rung $`m_x=T_i^x_g`$. The latter approaches one half for $`\lambda ^1\mathrm{}`$. It shows the opposite behaviour compared to the complementary charge order parameter $`m_z=|T_i^z_g|=\eta /2`$. The latter is shown in Fig.(3) as given by Eq.(66) with the exact 1D Ising exponent $`\beta `$=$`1/8`$.

With these necessary ingredients it is now straightforward to calculate the superexchange integrals $`J^{(s)}`$, $`J^{(\tau )}`$ as function of the microscopic model parameters according to Eqs.(50), (51). The hopping and interaction parameters fixed from LDA+U- calculations are discussed in Sec. II where definite values for $`t_a,t_b,t_d`$ and U were given. The intra- rung repulsion $`V_a`$ was assumed to lie close to 0.5 eV and the inter- rung repulsion $`V_b`$ is allowed to vary in a wide range up to 1 eV. The dimensionless control parameter for CO is $`\lambda `$=$`V_b/2t_a`$ with $`\lambda _c`$=1 at the CO transition.
The dependence of $`J^{(s)}`$, $`J^{(\tau )}`$ on $`V_b`$ is shown in Fig.(4) for two cases: The dashed line corresponds to exchange integrals calculated in the modified MF approach, Eq.(68), with the Ising exponent $`\beta `$=$`1/8`$. For $`\lambda <1`$ ($`V_b<V_b^{(\text{cr.})}=2t_a=0.76`$ eV) in the disordered ladder the exchange constants are independent of $`V_b`$ and the total value $`J=J^{(s)}+J^{(\tau )}`$ amounts to $`J_{\text{MF}}=J_{\text{MF}}^{(s)}80`$ meV since $`J_{\text{MF}}^{(\tau )}=0`$ in this region. At $`V_b=V_b^{(\text{cr})}`$ a sharp drop in the singlet exchange constant is observed whose shape is determined by the order parameter $`m_z`$ (Fig.(3)). The most important reason for the reduction of superexchange is the reduced effective hopping between adjacent rung bonding states caused by charge order. The full line in Fig.(4) shows the same exchange integrals with the effect of charge fluctuations included, as described by the pseudo- spin formalism of the previous chapter, Eqs.(70)-(72). In the disordered regime the inclusion of charge fluctuations strongly reduces the exchange constants compared to the MF value. The total superexchange integral falls now into the region 70 meV$`>J>55`$ meV, which is somewhat lower than the theoretical value obtained in Ref.() for the disodered phase while it is somewhat higher than the value $`J^{(\text{exp})}50`$ meV estimated in Ref.(). Fig.4 shows that in the CO regime ($`V_b>V_b^{(\text{cr})}`$), where charge fluctuations are strongly suppressed, both approximations give very similar results. Above the critical region at, for instance, $`V_b1.3V_b^{(\text{cr})}`$ (i.e. at $`V_b=1`$ eV), where CO is nearly complete, $`m_z0.43`$, the superexchange integral reaches the value $`J0.22`$ meV.
In Fig.(5) we show $`J^{(s)}`$ and $`J^{(\tau )}`$as function of Vb, with charge fluctuation effects included, for various intermediate values of the intra-rung Coulomb repulsion $`V_a`$. Moreover, by varying $`V_a`$ in a wider range from zero up to rather high values, $`V_a4t_a`$, a regular behaviour of the exchange constants is found, which is already seen from Fig.(5). First, the triplet contribution $`J^{(\tau )}`$ is almost negligible as compared to the singlet part $`J^{(s)}`$ for the model with $`t_d=t_b`$. In the case $`t_d`$ =0 however $`J^{(\tau )}`$ gives a sizeable (negative) contribution to the total $`J`$. Secondly, with increasing $`V_a`$, the total superexchange integral $`J`$ decreases and remains antiferromagnetic ($`>0`$), it never changes its sign. This is true for both models considered ($`t_d=t_b`$ or $`t_d=0`$). These observations are at strong variance with the results of Ref.()

It is instructive to examine also both the limit of the extreme CO, $`m_z1/2`$, which is reached at $`\lambda ^10`$, i.e. at unphysicaly high values of $`V_b`$, and the the opposite limit $`\lambda ^1\mathrm{}`$, i.e. at $`V_b=0`$. To check these limits a variation of the total superexchange integral $`J`$ is examined in a wide range of $`\lambda ^1`$, Fig.6(a). Two curves there correspond to two sets of the hopping parameters used: (1) $`t_a=0.38`$eV, $`t_b=t_d=0.085`$eV and (2) $`t_a=0.38`$eV, $`t_b=0.17`$eV and $`t_d=0`$. The second set is the most representative one among the sets with $`t_d=0`$ used by other authors . First, at $`\lambda ^10`$ the curve (1) reaches the value $`J=7.2`$meV which coincide exactly with the estimate based on the standard expression $`J=4t_d^2/U`$. The latter estimate is due to the superexchange coupling between two vanadium ions located at the opposite ends of $`nn`$ rungs. From Fig.6 one can see that in the case of finite $`t_d`$ the superexchange is regularly enhanced as compared to the case with $`t_d=0`$. In fact we found numericaly that if $`\lambda ^1\mathrm{}`$, the two dependences, i.e. the curves (1) and (2) in Fig.6, approach the same value $`J=80`$meV which is nothing than the MF value of $`J`$ depicted in Fig.4. In this limit ($`V_b=0`$) the general formulas, Eqs.(50), (51), can be easily treated analytically to give the following expression for the total superexchange integral

$`J=2t^2\left\{{\displaystyle \frac{\left(1+\mu \right)^2}{1+\mu ^2}}{\displaystyle \frac{1}{2t_aJ_a+V_a}}+{\displaystyle \frac{\left(1\mu \right)^2}{1+\mu ^2}}{\displaystyle \frac{1}{2t_a+J_a+U}}\right\}`$ (73)
where $`\mu =2t_a/(UV_a)`$; $`J_a=(2t_a)^2/(UV_a)`$ and $`t=t_b+t_d`$ for the set (1) or $`t=t_b`$ for the set (2). By using the same parameter values as in Fig.6 one finds that Eq.(73) reproduces completely the above mentioned numerical value for $`J`$. A large change in $`J`$ seen in Fig.6(a) at $`\lambda =\lambda _c`$ is provoked by a sharp behaviour of the order parameter $`m_z`$ as the system undergoes the CO phase transition. Considering now the CO regime only ($`\lambda ^1<1`$) we reproduce in Fig.6(b) the same variation of $`J`$, however, with respect to the order parameter $`m_z`$. Fig.6(b) shows a rather smooth behaviour with almost linear decrease of $`J`$ in the regime of weak charge disproportionation, $`m_z<0.2`$.

## VI Discussion and Conclusion

In Sec.II we formulated the underlying electronic model for the compound NaV2O$`_5`$with most of the model paremeters ($`t_a,t_b,t_d`$ and U) fixed from known results of the DFT calculations. By using this model in Sec.III we developed the superexchange theory and derived the second-order perturbation expression (in the effective V-V hopping amplitudes $`t_b`$ and $`t_d`$) for the intra-ladder exchange integral $`J`$ both in the charge disordered and in the zig-zag charge ordered state. Assuming that the Coulomb mechanism drives the charge ordering (Sec.IV) the exchange integral $`J`$ is calculated in a wide range of the variable $`V_b`$ and $`V_a`$ parameters. Since qualitative results are presented in details in Sec.V we summarize here most important qualitative findings.

1. In the intra-ladder superexchange mechanism the triplet channel is of minor importance for the model with $`t_d=t_b`$ and the resulting exchange integral is dominated by the singlet channel, $`JJ^{(s)}`$ (Figs.4,5). For the model with $`t_d=0`$, $`J^{(\tau )}`$ also gives a sizeable (negative) contribution to $`J`$.
2. In the charge disordered state inter-rung charge-charge correlations within a ladder cause a strong reduction of the exchange integral value (Fig.4).
3. With increasing of the zig-zag charge ordering the exchange integral $`J`$ reduces strongly but never changes sign, i.e. $`J`$ remains antiferromagnetic ($`>0`$) for the whole parameter range, Figs.4-6.
4. Variation of the underlying electronic model due to different choices of the hopping amplitude $`t_d`$ (either $`t_d0`$ or $`t_d=0`$) may lead to a considerable quantitative change in the value of $`J`$, especially in the extreme CO limit where $`J`$ stays finite or vanishes respectively (Fig.6).

In common, Figs.4-6 reveal a general trend. That is the decrease of the total exchange integral $`J`$ as the inter-rung charge-charge correlations within a ladder increase. We believe this trend together with less obvious properties (i)-(iv) are of general significance not restricted to the particular mechanism we used to describe the CO phase transition in the ladder system.

In the present paper we concentrated on a study of the the intra-ladder superexchange in NaV2O$`_5`$. To derive a more complete spin model in the Trellis lattice a thorough analysis of the inter-ladder magnetic interactions is required. Such an analysis for the charge disodered phase in NaV2O$`_5`$was done in Ref.(). There the authors claimed that due to the almost perfect cancellation of the triplet- and singlet-interactions in a-direction the ladders are effectively decoupled to give essentially 1D magnetic behavior of the entire spin system. In support of this conjecture we recall also that in the charge disordered Trellis lattice the inter-ladder Heisenberg spin-spin interactions are strongly frustrated. These arguments, however, cannot be transfered directly to a zig-zag CO state. The DFT calculation  of the exchange integrals in a zig-zag CO state have shown, for instance, a rather considerable ferromagnetic (FM) inter-ladder coupling . The opposite sign of this coupling is not surprising in view of different symmetry of V-O-V bonds connecting neighbouring ladders. It is a a challenging problem to reveal the microscopical origins and estimate a value of the resulting FM coupling between neighbouring ladders, which can be done in a properly developed superexchange theory taking into account also effects of direct V-V exchange.

## Acknowledgement

We are grateful to A. Bernert, P. Fulde, A. Yaresko and A.A. Zvyagin for numerous fruitful discussions. One of the autors (V.Yu.) thanks the INTAS organization (INTAS-97-11066) for support.

##

In this Appendix we exploit an extended model including inter-rung Coulomb interactions. Our aim is to show that these interactions lead to the symmetry breaking term $`\mathrm{\Delta }\epsilon ^0`$ in Eq.(65).

We choose the notation $`V_{ab}`$ for the Coulomb repulsion between the electrons on the diagonal V-V bond, i.e. on a pair of $`nn`$ V-ions belonging to adjacent ladders (Fig.(1)). We assume the parameter $`V_{ab}`$ to be not too large otherwise the instability with respect to the in- line charge ordering would be the dominant one. Now the necessary addition to the single ladder model, Eq.(1), can be presented in a form $`H_{ab}=_{\{ij\}}H_{ab}^{\{ij\}}`$, where the summation is over the pairs $`\{ij\}`$ of $`nn`$ rungs belonging to adjacent ladders. For a particular $`\{ij\}`$ pair one has

$$H_{ab}^{\{ij\}}=V_{ab}(i_L,j_R)n_i^Ln_j^R=V_{ab}(i_R,j_L)n_i^Rn_j^L$$
(74)
where the composite index $`i_L`$ ($`i_R`$) means the position of V-ion on the left (right) side of the $`i`$-th rung; the definitions for the electron density operators $`n_i^\alpha (\alpha =L,R)`$ are the same as in the main text. In the pseudospin representation the term $`H_{ab}`$ takes the following form

$`H_{ab}`$ $`=`$ $`{\displaystyle \underset{i}{}}\mathrm{\Delta }\epsilon _i^0T_i^z+{\displaystyle \underset{\{ij\}}{}}V_{ab}(i_L,j_R)\left({\displaystyle \frac{1}{4}}T_i^zT_j^z\right)=`$ (76)
$`=H_{ab,1}+H_{ab,2}`$
where

$`\mathrm{\Delta }\epsilon _i^0={\displaystyle \frac{1}{2}}{\displaystyle \underset{ji}{}}\left[V_{ab}(i_L,j_R)V_{ab}(i_R,j_L)\right]`$ (77)
and summation in (77) runs over the rungs neighbouring to the $`i`$–th rung.

Without distortion in the Trellis lattice all the diagonal bonds are identical, $`V_{ab}(i_L,j_R)=V_{ab}(i_R,j_L)=V_{ab}(|\stackrel{}{r_0}|)`$. In this case $`\mathrm{\Delta }\epsilon _i^0=0`$ and the inter-ladder coupling is given by the bilinear term $`H_{ab,2}`$. Assuming a zig-zag charge ordering in the ladders one may check that in the mean field approximation this term is averaged to zero due to a special geometry of the ladder stacking. It means that on the mean field level the ladders in the non-distorted Trellis lattice are disconnected. Below we argue briefly that a symmetry breaking source appears if the lattice symmetry is lowered as reported, for instance, in Refs.(,) for NaV2O$`_5`$.

Lüdecke et al.  proposed for the low- T lattice superstructure of NaV2O$`_5`$two species of ladders, denoted as A and B in Fig.1. In each A- ladder there is a weak transverse, along the a-axis, modulation of the rung positions while the lattice structure of the B- ladders remain symmetric as in the high-T phase. The doubling of unit cell means, in particular, that in the A-B-A’-B .. array the rung’s position modulation in the A- and A’-ladders are of opposite signs. In the modulated Trellis lattice a shortening/lengthening ($`\pm \stackrel{}{\delta }`$) of the diagonal V-V bonds leads to a special spatial pattern for the intersite Coulomb parameter variation $`V_{ab}(|\stackrel{}{r_0}\stackrel{}{\delta }|)`$ $`V_{ab}(|\stackrel{}{r_0}|)V_{ab}\stackrel{}{\delta }V_{ab}^{()}`$ (see Fig.1). By using Eq.(77) one can easily check that this pattern produces in the A-ladders an effective longitudinal field which is staggered along the ladder b-direction, $`\mathrm{\Delta }\epsilon _{i,A}^0=(1)^i\mathrm{\Delta }\epsilon _A^0`$, with the amplitude $`\mathrm{\Delta }\epsilon _A^02|V_{ab}\stackrel{}{\delta }|`$. Moreover, the field $`\mathrm{\Delta }\epsilon _{i,A}^0`$ is also staggered along a-direction, which sustains the zig-zag CO state with the A- and A’- ladders being ordered in the antiphase manner. Due to smallness of the lattice distortion, $`|\stackrel{}{\delta }|/|\stackrel{}{r_0}|10^2`$, we expect that $`\mathrm{\Delta }\epsilon _A^0V_b`$ for the actual values of $`V_b`$ providing the CO phase transition.
It follows from Eq.(77) that a symmetry breaking field is not present explicitly in the B-subsystem, i.e. $`\mathrm{\Delta }\epsilon _B^0=0`$, for the spatial pattern in Fig.1. Strictly speaking, a presence of non-zero field $`\mathrm{\Delta }\epsilon _B^0`$ is not necessary since the long-range zig-zag CO in B-ladders, if it is present, can be explained due to spontaneous symmetry breaking if $`\lambda >\lambda _c`$.

Note that besides the extracted staggered field term the inter-ladder bilinear interaction, $`H_{ab,2}`$, is still present. With a zig-zag ordering in A-/or both in A- and B- subsystems this term being considered in the MF- approximation is averaged to zero even in the modulated Trellis lattice structure.