format: text

text:
# Staggered flux and stripes in doped antiferromagnets

## I Introduction

It is well known that topological spin-textures are important in quantum Hall ferromagnets.  It is also known that an electromagnetic flux through a system of tight binding electrons on a lattice can lower the electronic kinetic energy.  Topological arguments also suggest that a possible source of this flux could be the formation of a spin-texture, which would suggest a relationship between doping and magnetic flux.  However, for the two dimensional antiferromagnetic Heisenberg model it has been argued  that, taking the continuum limit and looking at long wavelength fluctuations about the Néel state, there is no topological term in the effective action. Although this argument is correct, it does not, however, answer the question whether or not spin-textures can be important on a length scale comparable with the lattice. Furthermore it does not address the issue of second-neighbor hopping. These ideas will be explored in the present paper which extends ideas presented previously.  In particular, we generalize the topological arguments given in Ref. \[\] and present a numerical comparison of the energy of flux generating spin textures and flux free spin configurations for uniform systems. Besides uniform systems, we also consider if spin- and charge-stripes arise naturally as a topological fictitious flux generating spin texture.

The paper is organized as follows: In Section II we introduce the $`tJ`$-model and provide some background material on its properties. We proceed in Section III by deriving the effective model which we will work with. This model turns out to include topological fluxes which are discussed in Section IV and in Section V we review the effect of such a flux on a system of free electrons. Section VI contains a numerical mean-field investigation of the energetics of the system, comparing flux generating states and more regular spiralling states where there is no flux generated. In Section VII we extend the discussion of uniform systems to include stripes forming antiphase domain walls between Néel ordered regions. We describe of stripe model and present data from numerical calculations comparing different stripes. Finally, we summarize our results in Section VIII.

## II The $`tJ`$ model

In order to explore these ideas we employ the $`tJ`$ model which most simply captures the competition between Heisenberg exchange and kinetic energy:

$$H=\underset{\mathrm{𝐫𝐫}^{}}{}\left[t\left(c_{𝐫\sigma }^{}c_{𝐫^{}\sigma }+\text{h.c.}\right)+J\left(𝐒_𝐫𝐒_𝐫^{}\frac{1}{4}n_𝐫n_𝐫^{}\right)\right].$$
(1)
The summation is restricted to nearest-neighbor pairs on the square lattice and the spin operator is given by $`𝐒_𝐫=\frac{1}{2}c_{𝐫\alpha }^{}\stackrel{}{\sigma }^{\alpha \beta }c_{𝐫\beta }`$, where $`\stackrel{}{\sigma }=(\sigma _x,\sigma _y,\sigma _z)`$ is the vector of Pauli matrices. All states containing doubly occupied sites have been excluded from the Hilbert space, leaving three states per site: $`|`$, $`|`$, and $`|`$. A natural generalization, deferred to later in this paper, is to add Coulomb-repulsion between particles occupying nearest-neighbor sites.

Striped phases have been found experimentally in high-$`T_c`$ materials to which the present model has been applied. There is an ongoing debate regarding the existence of stripes in the $`tJ`$\- and Hubbard models. Stripes were first found in Hartree-Fock solutions of the Hubbard model , but the stripes found in these calculations had one hole per unit length of the stripe, in contrast to the results from experiments where half a hole per unit stripe length is found. From DMRG calculations on finite systems (of the order $`20\times 10`$). White and Scalapino  find stripes in a wide range of dopings. For instance, using $`J=0.35t`$ they find stripes for dopings in the interval $`0<x<0.3`$. For $`x<0.125`$ the stripes have half a hole per unit cell of stripe, in agreement with experiments, and the distance between two consecutive stripes is $`d=(2x)^1`$. For higher dopings they find that there is one hole per unit cell of the stripe and that the inter-distance between the stripes is $`d=x^1`$. On the other hand, using quantum Monte-Carlo calculations Hellberg and Manousakis  find that uniform or phase-separated states are energetically favorable. In this case, the formation of stripes would also require the existence of a long-range Coulomb interaction preventing an ordinary phase-separation .

Incommensurate states were discussed in connection to the $`tJ`$ model before the notion of stripes was introduced. It was found by Shraiman and Siggia  using a continuum limit of the model, that the antiferromagnetic order of the undoped $`tJ`$ model is unstable against the formation of a spiral state for small dopings. Using various mean-field approaches other authors  came to similar conclusions. We will return to the spiral instability in Section VI using the effective model to be described in the next section. It will be shown how a small twist in the spin-order leads to a reduction of the kinetic energy of the order $`tx`$ while the loss in exchange energy is of the order $`Jx^2`$; showing that for small enough dopings there is energy to be gained by twisting the spin-order.

In addition to spin textures, Affleck and Marston  discussed the possibility of flux textures. They replaced the two spin components of the electrons by $`N`$ different flavors, extending the SU(2) spin symmetry to SU($`N`$). Taking the limit $`N\mathrm{}`$, they obtained an essentially exact mean-field model to which they numerically looked for solutions. In particular they found a phase, called the flux phase, where the sum of the phases of the link-operators $`\chi _{\mathrm{𝐫𝐫}^{}}=c_𝐫^{}c_𝐫^{}`$ around a plaquette equals $`\pm \pi `$. This is interpreted as half a flux quantum penetrating each plaquette. These phases do not come from a real electromagnetic field and are therefore referred to as fictitious.

Work by Hasegawa et al.  showed that the energy of non-interacting spinless fermions has a minimum when a uniform flux, corresponding to one flux quantum per particle, threads the system. States having $`\mathrm{\Phi }=n`$ (in units of the flux quantum) are referred to as commensurate flux states. This shows that a fictitious flux can lower the kinetic energy of the particles. The commensurate flux states have also been considered in connection with the $`tJ`$ model.

Another possibility for a flux state is to have a staggered flux through the system. In the case of half a flux quantum per plaquette there is no difference between uniform- and staggered fluxes, so the Affleck-Marston state can be thought to belong to this category as well. Inspired by the work of Shraiman and Siggia, Kane et al.  suggested a double spiral state showing a staggered chiral spin order, and hence also, according to a result of Wen et al. , a staggered fictitious flux. In a staggered flux state, the time-reversal symmetry can be broken locally but not globally, as the system is invariant under a time-reversal operation followed by a lattice translation, just like a Néel state. Staggered flux phases have also been investigated by other groups . The effective model used by Barford and Lu  coincides with the model derived in the next section.

Our paper expands on results found by previous authors; we find that certain spin textures and charged stripes in particular are coupled by the creation of a fictitious $`\pi `$ flux which we show is a natural consequence of a stripe with broken rotational symmetry.

## III Deriving an effective model

In order to make progress, we make certain simplifications of the $`tJ`$ model. First, following Schulz , we introduce a local quantization axis $`\widehat{\mathrm{\Omega }}_𝐫`$ at site $`𝐫`$. In terms of spherical coordinates we write $`\widehat{\mathrm{\Omega }}_𝐫=(\mathrm{sin}\theta _𝐫\mathrm{cos}\varphi _𝐫,\mathrm{sin}\theta _𝐫\mathrm{sin}\varphi _𝐫,\mathrm{cos}\theta _𝐫)`$. This local SU(2) transformation on $`c_𝐫`$, denoted by $`U_{\widehat{\mathrm{\Omega }}_𝐫}`$, must fulfill the equation

$$U_{\widehat{\mathrm{\Omega }}_𝐫}\sigma _zU_{\widehat{\mathrm{\Omega }}_𝐫}^{}=\widehat{\mathrm{\Omega }}_𝐫\stackrel{}{\sigma }.$$
(2)
As can be seen from the above equation, specifying $`\widehat{\mathrm{\Omega }}_𝐫`$ determines $`U_{\widehat{\mathrm{\Omega }}_𝐫}`$ only up to a rotation about the new local $`z`$-axis. For example, we may choose our SU(2) transformation according to

$$U_{\widehat{\mathrm{\Omega }}_𝐫}=\mathrm{exp}\left[i\frac{\theta _𝐫}{2}\widehat{\omega }_𝐫\stackrel{}{\sigma }\right],$$
(3)
where $`\widehat{\omega }_𝐫=\widehat{(\widehat{z}\times \widehat{\mathrm{\Omega }}_𝐫)}=(\mathrm{sin}\varphi _𝐫,\mathrm{cos}\varphi _𝐫,0)`$.

Expressing the $`tJ`$ Hamiltonian in terms of this new spin-coordinate system we find

$`H={\displaystyle \underset{\mathrm{𝐫𝐫}^{}}{}}`$ $`[`$ $`t\left(c_{𝐫\alpha }^{}M_{\alpha \beta }^{\mathrm{𝐫𝐫}^{}}c_{𝐫^{}\beta }+\text{h.c.}\right)`$ (5)
$`+J(S_𝐫^\alpha S_𝐫^{}^\beta Q_{\alpha \beta }^{\mathrm{𝐫𝐫}^{}}{\displaystyle \frac{1}{4}}n_𝐫n_𝐫^{})],`$
with

$`M^{\mathrm{𝐫𝐫}^{}}`$ $`=`$ $`(U_{\widehat{\mathrm{\Omega }}_𝐫})^{}U_{\widehat{\mathrm{\Omega }}_𝐫^{}}`$ (6)
$`Q^{\mathrm{𝐫𝐫}^{}}`$ $`=`$ $`R_{\widehat{\mathrm{\Omega }}_𝐫}^1R_{\widehat{\mathrm{\Omega }}_𝐫^{}}`$ (7)
$`(R_{\widehat{\mathrm{\Omega }}})_{ij}`$ $`=`$ $`\mathrm{cos}\theta \delta _{ij}+(1\mathrm{cos}\theta )\omega _i\omega _j\mathrm{sin}\theta ϵ_{ijk}\omega ^k.`$ (8)
We note that $`R_{\widehat{\mathrm{\Omega }}}`$ is the SO(3) rotation-operator induced by the SU(2)-transformation $`U_{\widehat{\mathrm{\Omega }}}`$.

Thinking of the $`tJ`$ model as the large-$`U`$ limit of the Hubbard model, we know that there is a gap between the Hubbard bands scaling as $`U`$, corresponding to the energy cost for a double occupancy. Following Schulz , who neglected holes in the lower Hubbard band, we will throw away the upper Hubbard band because of this large gap when we consider hole-doping, corresponding to the removal of states containing double occupancies from the Hilbert space of the Hubbard model. Since the quantization axis at a site is locally determined by $`\widehat{\mathrm{\Omega }}_𝐫`$, we can arbitrarily assume that the upper Hubbard band is associated with spin down relative $`\widehat{\mathrm{\Omega }}`$. Hence, our effective model is obtained by keeping only the terms in Eq. (5) associated with spin-up particles. The spin of an electron at site $`𝐫`$ will now be determined by the field $`\widehat{\mathrm{\Omega }}_𝐫`$. As the simplest approximation, we will consider the $`\widehat{\mathrm{\Omega }}`$-field as a classical field, neglecting spin-fluctuations in the system. Keeping only terms containing particles aligned with the positive local $`z`$-axis we obtain an effective Hamiltonian

$$H_{\mathrm{eff}}=\underset{\mathrm{𝐫𝐫}^{}}{}\left[(\tau ^{\mathrm{𝐫𝐫}^{}}c_𝐫^{}c_𝐫^{}+\text{h.c.})+K^{\mathrm{𝐫𝐫}^{}}n_𝐫n_𝐫^{}\right]$$
(9)
with $`\tau ^{\mathrm{𝐫𝐫}^{}}=tM_{11}^{\mathrm{𝐫𝐫}^{}}`$, $`K^{\mathrm{𝐫𝐫}^{}}=\frac{1}{4}J(\widehat{\mathrm{\Omega }}_𝐫\widehat{\mathrm{\Omega }}_𝐫^{}1)`$, and $`c_𝐫=c_𝐫`$. This Hamiltonian describes a system of spinless fermions moving in a lattice with position dependent hopping amplitudes and interaction strengths.

The effects of the spin texture on the charge-motion is therefore to generate an effective hopping amplitude, $`tt\mathrm{cos}\frac{\theta }{2}`$ and the appearance of a fictitious magnetic field. We also note that there is a Coulomb like nearest-neighbor interaction of strength $`\frac{1}{4}J(\mathrm{cos}\theta 1)`$. When $`\mathrm{cos}\theta <1`$ this leads to an effective attraction between particles, which hints that the system may favor a phase-separation when being doped.

We have already mentioned that there is a degree of freedom in $`U_{\widehat{\mathrm{\Omega }}}`$ not being fixed by $`\widehat{\mathrm{\Omega }}`$. Since the effect of a spin-rotation about the local $`z`$-axis on the “up” spin only introduces a phase factor, it will be indistinguishable from a local electromagnetic gauge transformation in our approximation. Hence, the set of physically inequivalent choices of $`U_{\widehat{\mathrm{\Omega }}}`$ are determined by $`\widehat{\mathrm{\Omega }}`$, i.e. they belong to SU(2)/U(1)$`S^2`$.

The physical degrees of freedom of the hopping are contained within the size of the hopping amplitude and the gauge invariant parts of the complex phases of the hopping elements. In case of nearest-neighbor hopping only, the smallest closed loop that can be formed is around a plaquette in the lattice, see Fig. 1, and the flux enclosed by such a counter-clockwise path $`01230`$ is given by

$$\mathrm{\Phi }_{0123}=\mathrm{Im}\mathrm{ln}(\tau ^{01}\tau ^{12}\tau ^{23}\tau ^{30}).$$
(10)
One may show that this flux $`\mathrm{\Phi }_{0123}`$ is equal to half the solid angle enclosed by the shortest path on the spin sphere connecting the points $`\{\widehat{\mathrm{\Omega }}_i\}_{i=0}^3`$. Thus the flux is equal to $`2\pi Q`$, where $`Q`$ is the topological charge represented by the plaquette.

In the following discussion we will refer to the flux as a fictitious flux, in contrast to a “real” electromagnetic flux that would come from an applied magnetic field. The reason for this distinction is that the fictitious spin generated flux is only seen by the spinless fermions in the system and is unrelated to a physical electromagnetic field. Furthermore, the fictitious flux does not couple to the charge of the fermions, but rather to the $`z`$-component of the spin measured in the local spin-coordinate system. Since all particles in our system are polarized along the positive $`z`$-axis, they will appear as having the same fictitious charge. However, the flux can still drive currents through the system and in principle it is possible for the fictitious flux to cancel the effect of an external electromagnetic flux on the particles in the system.

A physical magnetic field giving one flux quantum per plaquette is enormous. If we assume that the lattice constant is $`a5`$ Å, the resulting magnetic field is $`B=h/(ea^2)10^4`$ T. This energy is much larger than the typical electronic energies per site, being of the order a few eV.

## IV Properties of fictitious fluxes

In Ref. \[\] we investigated the relations between the magnitude of the hopping $`\tau ^{\mathrm{𝐫𝐫}^{}}`$ and its complex phases for both ferromagnetic and antiferromagnetic spin-configurations. In the present section we shall briefly review and generalize these arguments.

We again consider a square lattice with a set of spins placed at each lattice site $`𝐫`$. The interior angles on the surface of the sphere are described by angles $`\alpha _i`$ as can be seen in Fig. 1. The fictitious flux through the plaquette is equal to half the solid angle covered by the plaquette in spin space, which by spherical geometry is given by the sum of the interior angles in excess of $`2\pi `$. It is obvious from Fig. 1 that if $`\theta _{\mathrm{𝐫𝐫}^{}}`$, the angle subtended by the arcs on the sphere is small, the area of the spherical parallelogram, and hence the fictitious flux, will be small as well. The following expressions give the size of the hopping and the fictitious flux through the plaquette:

$`|\tau ^{\mathrm{𝐫𝐫}^{}}|`$ $`=`$ $`t\mathrm{cos}{\displaystyle \frac{\theta }{2}}`$ (11)
$`\mathrm{\Phi }_{0123}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}\alpha _i\pi `$ (12)
For small values of opening angles on the spin sphere $`\theta _{\mathrm{𝐫𝐫}^{}}\theta `$ and $`\alpha _i\alpha _{2+i}`$ the flux is approximately given by $`\mathrm{\Phi }_{0123}\frac{\theta ^2}{2}\mathrm{cos}\frac{\alpha _1\alpha _2}{2}`$, showing that the flux is bounded by $`|\mathrm{\Phi }|\frac{\theta ^2}{2}`$.

If we instead turn our attention to antiferromagnetic configurations described by $`\theta \frac{\pi }{2}`$, the situation changes drastically. The path taken in spin space when going around a plaquette in an antiferromagnetic configuration is shown in Fig. 2. If we denote the antiferromagnetic (staggered) spin at site $`𝐫`$ by $`\widehat{\mathrm{\Omega }}_𝐫`$, the path taken is $`\widehat{\mathrm{\Omega }}_0\widehat{\mathrm{\Omega }}_1^{}\widehat{\mathrm{\Omega }}_2\widehat{\mathrm{\Omega }}_3^{}\widehat{\mathrm{\Omega }}_0`$, where $`\widehat{\mathrm{\Omega }}^{}=\widehat{\mathrm{\Omega }}`$ denotes the antipodal point on the sphere.

Redefining $`\theta `$ to denote the opening angle between two neighboring staggered spins, we find the following relationships:

$`|\tau ^{\mathrm{𝐫𝐫}^{}}|`$ $`=`$ $`t\mathrm{sin}{\displaystyle \frac{\theta }{2}}`$ (13)
$`\mathrm{\Phi }_{01^{}23^{}}`$ $`=`$ $`\pi {\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}(1)^i\alpha _i`$ (14)
We note that the size of the fictitious flux is now completely decoupled from the opening angle $`\theta `$, as long as $`\theta \frac{\pi }{2}`$. Furthermore, in the antiferromagnetic case the fictitious flux is staggered since the path on a neighboring plaquette will be traversed with the sublattices exchanged.

If we allow for the possibility of next-nearest-neighbor (diagonal) hopping, we have to take into account additional gauge invariant fluxes. There are four of these for each square plaquette, defined by the removal of one of the four vertices of the plaquette. In Ref. \[\] we show that there is a topological constraint relating these four fluxes. This constraint takes the form

$$\mathrm{\Phi }_{01^{}2}+\mathrm{\Phi }_{023^{}}\mathrm{\Phi }_{01^{}3^{}}\mathrm{\Phi }_{1^{}23^{}}=2\pi n,$$
(15)
where $`n=0`$ for a ferromagnet and $`n=\pm 1`$ for an antiferromagnet. In case of a ferromagnet a prime does not denote the antipodal point, but rather the point itself. This relation is easy to verify by looking at Fig. 2 noting that the sphere is exactly covered by the four regions in Eq. (15). The topological constraint, Eq. (15), does not rely on the assumption of a spherical parallelogram and does also hold in the presence of an external electromagnetic flux. Counting degrees of freedom, we know that there are two degrees of freedom per plaquette (or site) in choosing the spin-configuration and in addition we have one parameter coming from an external flux. All in all, there are three free parameters per plaquette and hence we expect that the four fluxes through the sub-triangles are related by a single constraint, given above.

## V Spinless free fermions with flux

Before turning to a more through investigation of the physics of the effective Hamiltonian, Eq. (9), we will discuss the effects of a flux through a system of free spinless fermions. Hasegawa et al.  investigated a system of free electrons on a square lattice with a uniform magnetic flux. They found that the energy is minimized when there is exactly one flux quantum per particle, i.e. the optimal flux per plaquette is related to the doping according to $`\mathrm{\Phi }=(1x)\mathrm{\Phi }_0`$, where $`\mathrm{\Phi }_0`$ is the flux quantum. To illustrate this effect we have diagonalized the Hamiltonian for such a system of free fermions in a uniform flux $`\mathrm{\Phi }`$ for different values of $`\mathrm{\Phi }`$, finding the single-particle energies $`ϵ_i(\mathrm{\Phi })`$. The total energy of the system is then found by summing up the single-particle energies according to

$$E(\mathrm{\Phi },n)=\underset{ϵ_i(\mathrm{\Phi })<ϵ_F(n,\mathrm{\Phi })}{}ϵ_i(\mathrm{\Phi }),$$
(16)
where $`ϵ_F(n,\mathrm{\Phi })`$ is the Fermi energy corresponding to filling $`n`$ and flux $`\mathrm{\Phi }`$.

In Fig. 3, we plot $`E(\mathrm{\Phi },n)E(0,n)`$, i.e. the energy per site for different fluxes compared to the flux free case. The figure clearly shows how the optimal flux changes with doping. We also note that the system is particle-hole symmetric and hence $`E(\mathrm{\Phi },n)=E(\mathrm{\Phi },1n)`$.

However, the flux that is generated by the antiferromagnetic skyrmions is staggered in which case the Hamiltonian can be exactly diagonalized. The spectrum of this system, assuming $`\theta _x=\theta _y=\theta `$, is

$`ϵ_𝐤(\mathrm{\Phi })`$ $`=\pm 2t\mathrm{sin}\left({\displaystyle \frac{\theta }{2}}\right)`$ (18)
$`\times \sqrt{\mathrm{cos}^2k_x+2\mathrm{cos}{\displaystyle \frac{\mathrm{\Phi }}{2}}\mathrm{cos}k_x\mathrm{cos}k_y+\mathrm{cos}^2k_y}.`$
In Fig. 4 we show a plot of the energy per site for different staggered fluxes. The figure shows that the optimal flux is either $`0`$ or $`\pm \pi `$ depending on doping, the only choices which are consistent with time-reversal invariance.

From this analysis we conclude that it is reasonable when searching for minimum energy spin textures to consider configurations supporting flux $`0`$ (coplanar configurations) and $`\pm \pi `$. In the next section we will construct a mean-field theory based on these observations.

### A Second-neighbor hopping

As we have seen in Section III, an antiferromagnetic spin-configuration strongly suppresses the effective nearest-neighbor hopping on the square lattice. However, a next-nearest-neighbor hopping is compatible with antiferromagnetic order. For the purpose of illustration let us consider the following purely kinetic Hamiltonian describing spinless fermions

$$H_{tt^{}}=t\underset{\mathrm{𝐫𝐫}^{}}{}\left[c_𝐫^{}c_𝐫^{}+\text{h.c.}\right]t^{}\underset{\mathrm{𝐫𝐫}^{\prime \prime }}{}\left[c_𝐫^{}c_{𝐫^{\prime \prime }}+\text{h.c.}\right],$$
(19)
where $`\mathrm{𝐫𝐫}^{\prime \prime }`$ denotes next-nearest-neighbor pairs. We remark that the sign of the nearest-neighbor hopping $`t`$ is irrelevant as it can be changed by the transformation $`c_𝐫(1)^𝐫c_𝐫`$. This transformation leaves the sign of $`t^{}`$ unchanged, and this sign is important. Without loss of generality, we assume $`t=1`$.

Another symmetry operation of interest is the particle-hole transformation $`c_𝐫(1)^𝐫c_𝐫^{}`$. Under this operation, the sign of the nearest-neighbor hopping is unchanged while the sign of $`t^{}`$ is changed, showing that next-nearest-neighbor hopping breaks the particle-hole symmetry. Furthermore, the number operator $`n_𝐫1n_𝐫`$ as particles are mapped into holes. This symmetry was seen in Figures 3 and 4.

Let us define $`E(n,t^{})`$ as the energy per site in the ground state of Eq. (19) with next-nearest-neighbor hopping $`t^{}`$. It is easy to show that as long as $`n<1/2`$, $`E(n,|t^{}|)<E(n,|t^{}|)`$ showing that for small fillings the energy is lower for the positive $`t^{}`$-case. If we instead consider the region $`n>1/2`$ the particle-hole transformation discussed above immediately tells us that $`E(n,|t^{}|)<E(n,|t^{}|)`$, showing that the negative $`t^{}`$-case is favorable. At precisely half-filling, the energy is independent of the sign of $`t^{}`$.

Assume that we consider the case $`n<1/2`$ and $`t^{}<0`$. Then, according to the discussion above, we would gain energy if we could somehow change the sign of $`t^{}`$. One way of accomplishing this would be to add a uniform flux through the system, corresponding to one flux quanta per square plaquette. This flux would not affect the nearest-neighbor hopping, but it would change the sign of $`t^{}`$ and therefore lower the energy of the system. Barford and Kim  generalized the result of Hasegawa et al. to include the $`tt^{}`$ model above and found that in the thermodynamic limit, the kinetic energy is minimized by a flux corresponding to one flux quantum per site plus or minus one flux quantum per particle.

## VI Hartree-Fock theory of uniform phases of $`H_{\mathrm{eff}}`$

It has been recognized for some time  that a plausible response of a Heisenberg antiferromagnet to doping, is to form a spiral spin-wave where the Néel order-parameter rotates uniformly around a fixed spin-axis as one moves along a symmetry axis in the lattice. This together with the fact that the doped electron gas favors a staggered flux close to half-filling, indicates that a state containing an antiferromagnetic spin-texture that generates both a staggered fictitious flux and a spiral-like order could lead to an energetically favorable state.

In the remaining sections of this paper, we will address two related questions. First, we investigate the effective Hamiltonian, Eq. (9), looking for the spin-textures which provide the energetically most favorable uniform states. In particular, we are interested in whether or not the system chooses to incorporate fictitious fluxes. Furthermore, it is known that the $`tJ`$ model has a tendency to phase-separate. Concerning the thermodynamic stability of the spiral states, Hu et al.  found that in a Hubbard model, for small dopings, the spiral phase is unstable against phase-separation (for $`U/t10`$) or domain wall formation (for $`U/t10`$). For larger dopings, there are regions in the phase-diagram, located around $`U/t10`$, where the spiral phases are thermodynamically stable. This indicates that, in the $`tJ`$ model, for small dopings the spiral state is not thermodynamically stable. An interesting question is; if the $`tJ`$ model prefers a flux-phase in some region of parameter space, can this thermodynamically stabilize the system, preventing it from phase-separation?

Inspired by the recent interest in striped phases, in Section VII we use our approach to model different domain walls between Néel-ordered regions. These domain walls have an appealing structure as they provide a smooth implementation of the antiphase boundary and at the same time provides the electrons in the doped channel with a fictitious flux.

### A MFT-formulation

We now look for different uniform phases of the effective Hamiltonian given in Eq. (9). The coefficients $`\tau ^{\mathrm{𝐫𝐫}^{}}`$ and $`K^{\mathrm{𝐫𝐫}^{}}`$ are now given by

$`\tau ^{\mathrm{𝐫𝐫}^{}}`$ $`=`$ $`t\mathrm{sin}{\displaystyle \frac{\theta _{\mathrm{𝐫𝐫}^{}}}{2}}\mathrm{exp}[i𝒜(\widehat{\mathrm{\Omega }}_𝐫,\widehat{\mathrm{\Omega }}_𝐫^{},\widehat{z})/2]`$ (20)
$`K^{\mathrm{𝐫𝐫}^{}}`$ $`=`$ $`{\displaystyle \frac{J}{4}}(1+\mathrm{cos}\theta _{\mathrm{𝐫𝐫}^{}}),`$ (21)
where $`\widehat{\mathrm{\Omega }}_𝐫`$ denotes the local staggered spin-orientation at site $`𝐫`$, and $`\mathrm{cos}\theta _{\mathrm{𝐫𝐫}^{}}=\widehat{\mathrm{\Omega }}_𝐫\widehat{\mathrm{\Omega }}_𝐫^{}`$. With this definition $`\mathrm{cos}\theta _{\mathrm{𝐫𝐫}^{}}=1`$ for a Néel state and $`\mathrm{cos}\theta _{\mathrm{𝐫𝐫}^{}}=1`$ for a ferromagnet, where $`𝐫`$ and $`𝐫^{}`$ are nearest neighbors. $`𝒜(\widehat{\mathrm{\Omega }}_𝐫,\widehat{\mathrm{\Omega }}_𝐫^{},\widehat{z})`$ is the solid angle of the spherical triangle spanned by the vectors $`\widehat{\mathrm{\Omega }}_𝐫`$, $`\widehat{\mathrm{\Omega }}_𝐫^{}`$, and $`\widehat{z}`$.

The approach we will use is a simple mean-field theory assuming a fixed spin-texture $`\{\widehat{\mathrm{\Omega }}_𝐫\}`$, defined as the direction of the quantization axis, $`\widehat{\mathrm{\Omega }}_𝐫`$. We will assume that $`\theta _{\mathrm{𝐫𝐫}^{}}=\theta _x`$ when $`𝐫`$ and $`𝐫^{}`$ are nearest horizontal neighbors, and $`\theta _{\mathrm{𝐫𝐫}^{}}=\theta _y`$ when they are nearest vertical neighbors. First of all we perform a standard mean-field decomposition of the Hamiltonian, allowing only for mean-fields carrying no charge and momenta zero or $`𝐐=(\pi ,\pi )`$. This results in the following Hamiltonian:

$`H_{\mathrm{MF}}`$ $`=`$ $`{\displaystyle \underset{𝐤\mathrm{BZ}^{}}{}}\mathrm{\Psi }_𝐤^{}\{2t[(\mathrm{sin}{\displaystyle \frac{\theta _x}{2}}\mathrm{cos}k_x+\mathrm{cos}{\displaystyle \frac{\mathrm{\Phi }}{2}}\mathrm{sin}{\displaystyle \frac{\theta _y}{2}}\mathrm{cos}k_y)\sigma _3+\mathrm{sin}{\displaystyle \frac{\mathrm{\Phi }}{2}}\mathrm{sin}{\displaystyle \frac{\theta _y}{2}}\mathrm{cos}k_y\sigma _2]`$ (24)
$`{\displaystyle \frac{J}{2N}}[(2+\mathrm{cos}\theta _x+\mathrm{cos}\theta _y)\mathrm{\Delta }^0𝟙(\mathrm{𝟙}+\mathrm{cos}\theta _𝕩)\mathrm{cos}𝕜_𝕩\mathbb{\Delta }_{𝕔𝕩}^\mathrm{𝟛}\sigma _\mathrm{𝟛}(\mathrm{𝟙}+\mathrm{cos}\theta _𝕪)\mathrm{cos}𝕜_𝕪(\mathbb{\Delta }_{𝕔𝕪}^\mathrm{𝟚}\sigma _\mathrm{𝟚}+\mathbb{\Delta }_{𝕔𝕪}^\mathrm{𝟛}\sigma _\mathrm{𝟛})]\}\mathbb{\Psi }_𝐤`$
$`+{\displaystyle \frac{J}{4N}}\left[(2+\mathrm{cos}\theta _x+\mathrm{cos}\theta _y)(\mathrm{\Delta }^0)^2(1+\mathrm{cos}\theta _x)(\mathrm{\Delta }_{cx}^3)^2(1+\mathrm{cos}\theta _y)\left((\mathrm{\Delta }_{cy}^2)^2+(\mathrm{\Delta }_{cy}^3)^2\right)\right],`$
where we have introduced a two-component vector $`\mathrm{\Psi }_𝐤=(c_𝐤,c_{𝐤+𝐐})^t`$, mixing momenta $`0`$ and $`𝐐`$. We have introduced the Pauli matrices as a basis for the $`2\times 2`$ matrices coupling the $`\mathrm{\Psi }_𝐤`$’s, although we want to emphasize that they have nothing to do with spin in this context. The sum over momenta is reduced to half the Brillouin zone, defined by $`\mathrm{BZ}^{}=\{|k_x|+|k_y|\pi :\pi k_x,k_y<\pi \}`$. Furthermore, we have only kept those four order-parameters  that turn out to be non-zero numerically. These four fields are defined through

$$\{\begin{array}{ccc}\hfill \mathrm{\Delta }^0& =& _{𝐤\mathrm{BZ}^{}}\mathrm{\Psi }_𝐤^{}𝟙\mathbb{\Psi }_𝐤\hfill \\ \hfill \mathrm{\Delta }_{cy}^2& =& _{𝐤\mathrm{BZ}^{}}\mathrm{cos}k_y\mathrm{\Psi }_𝐤^{}\sigma _2\mathrm{\Psi }_𝐤\hfill \\ \hfill \mathrm{\Delta }_{cx}^3& =& _{𝐤\mathrm{BZ}^{}}\mathrm{cos}k_x\mathrm{\Psi }_𝐤^{}\sigma _3\mathrm{\Psi }_𝐤\hfill \\ \hfill \mathrm{\Delta }_{cy}^3& =& _{𝐤\mathrm{BZ}^{}}\mathrm{cos}k_y\mathrm{\Psi }_𝐤^{}\sigma _3\mathrm{\Psi }_𝐤\hfill \end{array},$$
(25)
where the average $`\text{ }\text{ }`$ denotes a thermal expectation value with respect to the Fermi-distribution of quasiparticles of $`H_{\mathrm{MF}}`$. The order-parameter $`\mathrm{\Delta }^0`$ is simply the number of particles in the system, while the other three correspond to hopping induced through the term $`n_𝐫n_𝐫^{}`$ in the effective Hamiltonian, Eq. (9). In particular, we note that $`\mathrm{\Delta }^0`$ and $`\mathrm{\Delta }_{cx/y}^3`$ are diagonal and hence do not mix momenta $`𝐤`$ and $`𝐤+𝐐`$. On the contrary, $`\mathrm{\Delta }_{cx}^2`$ does mix the two, and therefore carries a momentum $`𝐐`$. As can be seen from Eq. (24), this term only exists in the kinetic term when there is a non-zero staggered flux which reduces the translational symmetry of the model.

### B Instability towards spiral-order at low dopings

Before turning to the numerical results, let us now discuss the electron gas in Eq. (9) at low dopings, confining our discussion to coplanar spin configurations and neglecting exchange effects in the Heisenberg term so that all order-parameters, except $`\mathrm{\Delta }^0`$, in Eq. (25) are zero. We choose a spin-structure $`\widehat{\mathrm{\Omega }}_𝐫\widehat{\mathrm{\Omega }}_{𝐫+\widehat{x}}=\mathrm{cos}\theta _x`$ and $`\widehat{\mathrm{\Omega }}_𝐫\widehat{\mathrm{\Omega }}_{𝐫+\widehat{y}}=\mathrm{cos}\theta _y`$, allowing for an asymmetry between the $`x`$\- and $`y`$-directions. This yields a trivial system which is exactly diagonalized. The total energy per site as a function of density $`n`$ and $`\theta _i`$ ($`i=x,y`$) is given by

$`E_{\theta _x,\theta _y}(n)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{ϵ_𝐤ϵ_F(n)}{}}ϵ_𝐤(\theta _x,\theta _y)`$ (27)
$`{\displaystyle \frac{1}{4}}Jn^2\left(\mathrm{cos}\theta _x+\mathrm{cos}\theta _y+2\right)`$
where $`ϵ_F(n)`$ is the Fermi energy corresponding to density $`n`$ and

$$ϵ_𝐤(\theta _x,\theta _y)=2t\left[\mathrm{sin}\frac{\theta _x}{2}\mathrm{cos}k_x+\mathrm{sin}\frac{\theta _y}{2}\mathrm{cos}k_y\right].$$
(28)
Our description of the spin-order in terms of the $`\theta `$-angles does not distinguish between spiral states and so called canted states. They both lack fictitious flux and both have the same relative angle between nearest-neighbor spins. Only second neighbor terms resolve this degeneracy. The difference between these two classes of states is illustrated in Fig. 5.

If we make a series expansion of the energy per site in terms of the doping $`x`$ we find the following expression to first order in $`x`$

$`E_{\theta _x,\theta _y}(x)`$ $`=`$ $`2t\left(\mathrm{sin}{\displaystyle \frac{\theta _x}{2}}+\mathrm{sin}{\displaystyle \frac{\theta _y}{2}}\right)x{\displaystyle \frac{1}{4}}J(12x)`$ (30)
$`\times (\mathrm{cos}\theta _x+\mathrm{cos}\theta _y+2)+O(x^2).`$
Remember that we have made the transformation $`\theta _i\pi \theta _i`$, expressing the order relative to the antiferromagnet instead of the ferromagnet. The energy is minimized by $`\theta _x=\theta _y=2\mathrm{arcsin}[2tx/J(12x)]4tx/J`$ for moderate dopings. The dependence on $`tx/J`$ clearly shows the competition between the kinetic energy, which drives the system towards ferromagnetism, and the Heisenberg term, which favors antiferromagnetism. These results are consistent with those of Schulz , who also found this instability. Shraiman and Siggia , using a more elaborate method, also found this instability, but their spiral state has its pitch-vector along the $`(10)`$-direction rather than the $`(11)`$-direction as is found here. In the Hubbard model, it is known from Hartree-Fock theory that the antiferromagnetic state is stabilized by the opening of a gap at the Fermi surface, and Schulz  has argued that a modulation of the spin-order along the $`(10)`$-direction opens up a gap more effectively than a spiral along the $`(11)`$-direction, i.e. it opens up a gap at a larger part of the Fermi surface. The first order theory described in this section does not take this fact into account. We note that the deviation from Néel-order is proportional to the doping and in the limit $`J0`$, the deviation from Néel-order becomes large for any finite doping. This is consistent with Nagaoka’s theorem , stating that in the limit $`J0`$, a single hole doped into a Heisenberg antiferromagnet drives the system into a ferromagnetic state. The energy of the spiral state is given by $`E=\frac{4t^2x^2}{J(12x)}J(12x)`$ and hence we have for the second derivative of $`E`$ with respect to the filling $`n`$,

$$\frac{^2E}{n^2}=\frac{8t^2}{J(12x)^3}.$$
(31)
We note that the energy versus filling is concave for dopings $`x<1/2`$, showing that the spiral state has an instability towards phase-separation for all $`J>0`$. Recall that our expansion is only valid for small dopings, and in this regime we expect the above conclusion to hold. Also, as mentioned earlier, we cannot distinguish between the spiral- and canting states within this approach. According to Kane et al.  quantum fluctuations seem to stabilize the spiral state compared to the canting state.

A similar analysis can be performed for fillings close to zero, i.e. $`n1`$. In this case, measuring $`\theta _i`$ with respect to the ferromagnetic configuration, the analogue of Eq. (30) becomes (now keeping terms up to second order in $`n`$)

$`E_{\theta _x,\theta _y}(n)`$ $`=`$ $`2t\left(\mathrm{cos}{\displaystyle \frac{\theta _x}{2}}+\mathrm{cos}{\displaystyle \frac{\theta _y}{2}}\right)n`$ (34)
$`+2\pi tn^2\sqrt{\mathrm{cos}{\displaystyle \frac{\theta _x}{2}}\mathrm{cos}{\displaystyle \frac{\theta _y}{2}}}`$
$`+{\displaystyle \frac{1}{4}}Jn^2(\mathrm{cos}\theta _x+\mathrm{cos}\theta _y2)+O(n^3).`$
In particular we note how the second term introduces a coupling between the spin-order in the $`x`$\- and $`y`$-directions, allowing for an asymmetry between $`\theta _x`$ and $`\theta _y`$. Minimizing the above energy with respect to the angles $`\theta _i`$, we find that the ferromagnetic state is stable up to a finite doping. At this point different things can happen depending on $`J/t`$, the system can pick a state where $`\theta _x=0`$ and $`\theta _y=\pi `$ (or vice versa), i.e. the system organizes itself ferromagnetically along the $`x`$-direction while being an antiferromagnet along the $`y`$-direction. Another possibility is that the system chooses a pitch-vector along $`(11)`$ with $`\theta _x=\theta _y=2\mathrm{arccos}\frac{(2n\pi )t}{Jn}`$.

The picture we have obtained is therefore that starting at zero filling, the system remains in a ferromagnetically ordered state up to some threshold value of the filling. This threshold increases with decreasing values of the ratio $`J/t`$. Above this threshold filling, the system can be a spiral spin-wave with pitch along $`(10)`$ or $`(11)`$. When the filling approaches $`n=1`$, a $`(11)`$ spiral state is optimal which continuously merges with the Néel-state as $`n1`$. But we find in all cases that the system is unstable against phase-separation for small dopings.

### C Numerics

Given the thermodynamic instability of the spiral spin waves, we use the full Hartree-Fock theory and take into account uniform phases that have a splay in the Néel order-parameter. The aim is to search the space of spin-textures, parametrized by $`(\theta _x,\theta _y,\mathrm{\Phi })`$, to determine the one which minimizes the free energy of the system. We will consider the following two types of spin-textures:

1. Coplanar states, described by $`\mathrm{\Phi }=0`$. (35)
2. $`\pi `$-flux state, described by $`\theta _x=\theta _y`$, $`\mathrm{\Phi }=\pi `$. (36)
The argument for only considering $`\mathrm{\Phi }=0`$ and $`\mathrm{\Phi }=\pi `$ states was given in Section V; Fig. 4 showed that the energy per site for a system of free fermions in a staggered flux is minimized by either of these two choices.

Our numerical algorithms work within the grand canonical ensemble, with a free energy $`G(T,\mu )=H_{\mathrm{MF}}TS\mu N`$, assuming a fixed chemical potential $`\mu `$. Diagonalizing the mean-field Hamiltonian, Eq. (24), we obtain a set of quasiparticle states specified by their momenta $`𝐤\mathrm{BZ}^{}`$, where $`\mathrm{BZ}^{}`$ is the reduced Brillouin zone corresponding to $`|k_x+k_y|\pi `$, and band-index. The band index refers to the two bands occurring because of the staggered flux. $`\alpha \{1,2\}`$. If we then minimize the free energy $`G(T,\mu )`$ with respect to the occupation numbers $`f_{𝐤\alpha }`$ of the quasiparticles, we find that they are distributed according to the Fermi-distribution. The entropy, $`S`$, introduced above is defined through,

$$S=k_B\underset{\alpha ,𝐤\mathrm{BZ}^{}}{}\left(f_{𝐤\alpha }\mathrm{ln}f_{𝐤\alpha }+(1f_{𝐤\alpha })\mathrm{ln}(1f_{𝐤\alpha })\right),$$
(37)
where $`\alpha \{1,2\}`$ labels the two bands.

In the analysis of the numerical data we would rather consider the free energy as a function of the number of particles $`N`$ than the chemical potential $`\mu `$. This can be achieved by forming the Helmholtz free energy through the following Legendre transformation; $`F(N,T)=G+\mu N`$. Having the Helmholtz energy, we can use the Maxwell construction to discuss the thermodynamic stability of the phases.

Analyzing the Hartree-Fock theory involves the following procedure: Given a set of coupling constants, a spin-configuration, temperature $`T`$, and chemical potential $`\mu `$; pick a set of initial values of the mean-fields. Then solve for the quasiparticles of Eq. (24) and calculate the new mean-fields using Eq. (25). The procedure is then iterated until the mean-fields have converged to a fixed point corresponding to a minimum in the free energy. Since we are interested in the spin-configuration minimizing the free energy $`F(T,N)`$, we will manually vary the spin-configuration parameters searching for a global minimum of the free energy.

In practice, rather than choosing a certain value of the chemical potential, we choose a fixed filling, successively adjusting the chemical potential during the iterations. Since the chemical potential was assumed to be fixed during the derivation of the self-consistency equations, we have to make sure that the algorithm converges to the correct fixed point. We have checked explicitly in several cases that the correct fixed point is found. An advantage with this method is that we can access all possible fillings, whether or not it is a thermodynamically stable region. This is not the case if we specify the chemical potential, since the function $`\mu (n)`$ is not invertible.

Let us start by considering a simple Hartree-approximation, see Figures 6 and 7. This corresponds to putting $`\mathrm{\Delta }_{cy}^2=\mathrm{\Delta }_{cx}^3=\mathrm{\Delta }_{cy}^3=0`$. Starting at zero temperature, we find that within the Hartree approximation there is a critical value $`J_c1`$ of the coupling $`J`$ below which the coplanar phase always dominates over the flux-generating configurations. When $`J>J_c`$, there appears a region around $`n0.5`$ where the flux-states are the energetically lowest states, see Fig. 7. When the flux-state minimizes the free energy, the $`\theta `$-angles are $`\pi /2`$, i.e. the maximally allowed values. This state corresponds to having the spins distributed along the equator with an angle $`\pi /2`$ between two successive spins. The solid angle spanned by this configuration covers half of the unit sphere, ensuring $`\mathrm{\Phi }=\pi `$. Thermodynamically, however, it seems to be favorable for the system to phase-separate into regions consisting of a hole-free antiferromagnet and a hole-rich coplanar structure, respectively. If we consider small dopings ($`x1`$), the $`\theta `$-angle of the optimal state (coplanar) is successively reduced to zero as $`x0`$. Since the maximum amplitude of the hopping for a staggered flux-state is $`|\tau |=\frac{t}{\sqrt{2}}`$ as determined by $`\theta \frac{\pi }{2}`$, it clearly has a disadvantage compared to the coplanar states, supporting $`|\tau |=t`$. The effect of this is that the coplanar state will always be favorable at fillings where the kinetic energy is dominant. However, as the filling is further increased, the Heisenberg energy becomes more important and it becomes favorable to reduce $`|\tau |`$ in order to improve the Heisenberg energy. At this point, if this occurs at a suitable filling, the flux-state can yield equally good Heisenberg and kinetic energies, while at the same time providing the fermions with a flux that can lower the energy even further. We know from the work of Hasegawa et al.  that the energy of a system of free electrons on a square lattice is minimized when there is exactly one flux quantum per particle. Since our flux-state carries a flux $`\pi `$ per plaquette, it will be most suitable close to $`n=1/2`$. The coplanar state on the other hand carries no flux, and will therefore be optimal when $`n=0`$ or $`n=1`$. This competition explains why the coplanar state becomes energetically favorable again when moving from a flux-phase towards higher fillings.

In our calculations, a slight technical point should be mentioned. When the system phase-separates into a coplanar part and an undoped antiferromagnet, the chemical potentials in the two subsystems are not equal. This can be seen, for instance in Fig. 6, where the Maxwell construction connects the antiferromagnetic $`n_1=1`$ point with a point $`n_00.55`$ in the coplanar phase. Since $`n_1`$ is an end-point of the free energy curve, the derivatives at $`n_0`$ and $`n_1`$ are unequal, and hence the chemical potential is different in the two phases.

### D Results from Hartree-Fock calculations

We now apply the Hartree-Fock scheme to the problem to understand whether or not exchange effects can resolve the near degeneracies found in the Hartree calculations. The spin-textures considered are those given by Eq. (36). Before starting the full calculations, we determined numerically which of the different order-parameters are the important ones. Initially there were ten  but numerically we find that only four are non-zero in the spin-configurations we have examined. To speed up the calculations we have set the other six to zero by hand. We find that besides the filling ($`\mathrm{\Delta }^0`$), $`\mathrm{\Delta }_{cy}^2`$, and $`\mathrm{\Delta }_{cx}^3`$ are important for the flux-phases, while $`\mathrm{\Delta }_{cx}^3`$ and $`\mathrm{\Delta }_{cy}^3`$ are important for the coplanar configurations.

Numerically we find that the $`\pi `$-flux-phase converges to a state where the non-zero order parameters are $`\mathrm{\Delta }_{cy}^2=\mathrm{\Delta }_{cx}^3=\mathrm{\Delta }`$, which leads to a quasiparticle spectrum of the form $`ϵ_𝐤=\stackrel{~}{\mu }\pm 2\stackrel{~}{t}\sqrt{\mathrm{cos}^2k_x+\mathrm{cos}^2k_y}`$, where $`\stackrel{~}{\mu }`$ and $`\stackrel{~}{t}`$ are renormalized values of the chemical potential and hopping amplitude, respectively. The momentum $`𝐤`$ belongs to the reduced Brillouin zone. At half-filling, this dispersion relation has four gapless Dirac-points where the energy vanishes linearly. These points are located at $`(k_x,k_y)=(\pm \pi /2,\pm \pi /2)`$. Similarly, for the coplanar states, the mean-fields renormalize the hopping amplitudes so that $`ϵ_𝐤=\stackrel{~}{\mu }\pm (\stackrel{~}{t_x}\mathrm{cos}k_x+\stackrel{~}{t_y}\mathrm{cos}k_y)`$.

If we look at a typical free energy plot, such as Fig. 6, we find, as in the Hartree case, that for low fillings the coplanar configuration is the optimal, where at low dopings it merges with the antiferromagnet. For small $`J`$’s the coplanar configuration clearly dominates over the flux-configuration for all fillings up to the point where they merge with the pure antiferromagnet. Increasing $`J`$ brings the flux-configuration energetically closer to the coplanar configuration. However, in contrast to the Hartree-case, it does not seem like the flux-state will become energetically favorable over the coplanar-states. Concerning phase-separation, the picture is very much the same as the one described above. For low temperatures and dopings smaller than roughly 0.5, the system favors a phase-separation into parts consisting of a hole-free antiferromagnet and a hole-rich coplanar state with doping $`x0.5`$.

Figures 6 and 7 show the free energies for the best coplanar- and flux-configurations plotted together with the free energy of the pure antiferromagnet in two different cases, $`J=0.75`$ and $`J=2.00`$. The energy-scale is fixed by $`t=1`$. The temperature is set to $`T=0.1`$ and the size of the system being considered is $`16\times 16`$ sites.

In Fig. 7, it is clearly seen how the flux-phase dominates over the coplanar-phase close to $`n0.5`$ at the Hartree-level, but not in the Hartree-Fock approximation.

We have generalized the $`tJ`$-model by including a nearest-neighbor Coulomb repulsion through a term $`V_{𝐫,𝐫^{}}n_𝐫n_𝐫^{}`$. In our effective model this corresponds to redefining $`K_{\mathrm{𝐫𝐫}^{}}K_{\mathrm{𝐫𝐫}^{}}+V`$. When including this term the order-parameter corresponding to a charge-density-wave, $`\mathrm{\Delta }^1=_{𝐤\mathrm{BZ}^{}}\mathrm{\Psi }_𝐤^{}\sigma _1\mathrm{\Psi }_𝐤`$, becomes important. As it turns out, a positive value of $`V`$ can favor the flux-phase as is seen in Fig. 8 where we have shown the free energy versus filling for a system described by $`J=1.25`$, $`K=4V/J=2.00`$, and temperature $`T=0.1`$. As can be seen from this figure, there appears a narrow region around $`n0.53`$ where the flux phase has the best energetics.

Finally, in Fig. 9, we show a phase-diagram as a function of filling and coupling $`K`$. In this data we have fixed the exchange coupling $`J=1.25`$ and temperature $`T=0.1`$. The figure shows how a narrow region of a staggered-flux phase occurs close to $`n0.5`$.

Summarizing our numerical results, we find that a staggered flux-phase can be energetically favorable compared to a coplanar spiral state. However, in the region of doping where this happens, the system generically seems unstable against phase-separation. Even if this indicates that the flux-phase does not occur as a thermodynamically stable phase in our effective model, it has good energetics and it would be interesting to investigate the effects of spin-fluctuations in this picture.

In the next section we discuss stripes within this framework. We know that stripes are antiphase boundaries between antiferromagnetically ordered regions. One way to model this is to have a spin-order twisting as the boundary is crossed. Since we also know that the stripes are doped, it is tempting to think that the twisting is such that a fictitious flux is generated. The picture is also appealing since it makes use of the instability of the doped antiferromagnet towards phase-separation. We will return to this topic in Section VII

### E Circulating currents

An important issue to address is whether or not the flux-states are accompanied by circulating currents in the system. To answer this question, we consider the current operator at site $`𝐫`$ in the $`\delta `$-direction, $`j_\delta (𝐫)`$, which can be identified from charge conservation together with the Heisenberg equation of motion,

$$\underset{\delta =\widehat{x},\widehat{y}}{}\left[j_\delta (𝐫)j_\delta (𝐫\delta )\right]=\frac{\rho (𝐫)}{t}=\frac{i}{\mathrm{}}[H_{\mathrm{eff}},\rho (𝐫)].$$
(38)
The result is a current operator taking the form

$$j_\delta (𝐫)=\frac{i}{\mathrm{}}(\tau _{𝐫,𝐫+\delta }c_𝐫^{}c_{𝐫+\delta }\mathrm{h}.\mathrm{c}.).$$
(39)
We decompose the current into uniform- and staggered parts as $`j_\delta (𝐫)=j_\delta ^u(𝐫)+(1)^𝐫j_\delta ^s(𝐫)`$. Using the mean-field decomposition we find that the expectation value of the uniform currents vanish. The uniform currents are proportional to some of the order-parameters that have been left out of the discussion. We know that our Hamiltonian is invariant under a lattice translation followed by a time-reversal operation, but this composite operation reverses the direction of the uniform currents which hence must vanish in a thermodynamic expectation value. The staggered currents take the form

$`j_x^s(𝐫)`$ $`=`$ $`{\displaystyle \frac{2t}{N\mathrm{}}}\mathrm{sin}\left({\displaystyle \frac{\theta _x}{2}}\right)\mathrm{\Delta }_{cx}^2`$ (40)
$`j_y^s(𝐫)`$ $`=`$ $`{\displaystyle \frac{2t}{N\mathrm{}}}\mathrm{sin}\left({\displaystyle \frac{\theta _y}{2}}\right)\left(\mathrm{\Delta }_{cy}^2\mathrm{cos}{\displaystyle \frac{\mathrm{\Phi }}{2}}\mathrm{\Delta }_{cy}^3\mathrm{sin}{\displaystyle \frac{\mathrm{\Phi }}{2}}\right).`$ (41)
The currents are gauge invariant, and the formal lack of symmetry between the currents in the $`x`$\- and $`y`$-directions is due to gauge choice. The order-parameters also are gauge dependent, and this restores the symmetry between currents in the two directions, see Fig. 10.

According to the previous subsection the non-zero order-parameters in the $`\pi `$-flux state are $`\mathrm{\Delta }_{cx}^3`$ and $`\mathrm{\Delta }_{cy}^2`$, and we conclude that there are no circulating currents in the $`\pi `$-flux state. Obviously, the flux-free case is also free of currents, as can be seen from the equations above, in this case it is $`\mathrm{\Delta }_{cy}^3`$ rather than $`\mathrm{\Delta }_{cy}^2`$ which is non-zero. Noting that in the case $`\mathrm{\Phi }=\pm \pi `$ the Hamiltonian is time-reversal invariant so the the one-particle states possibly carrying current must be degenerate in energy. In a thermodynamic ensemble, all states of the same energy are weighted equally, and the physical currents in a system cancel and there will be no staggered currents in the system. For any flux not being an integer multiple of $`\pi `$, this symmetry is lost and circulating currents appear.

If we calculate the above expectation values in a state having a flux $`\mathrm{\Phi }n\pi `$, we find circulating currents in the system.

A proper treatment of this problem should consist of a gauge invariant coupling to a real electromagnetic flux in addition to the spin-generated flux . Let $`\mathrm{\Phi }_{\mathrm{mag}}`$ denote the real electromagnetic flux through a plaquette, induced by circulating currents. The energy cost of creating the magnetic field should be added through a term

$$E_{\mathrm{mag}}=\frac{K}{2}\underset{\mathrm{plaquettes}}{}\mathrm{\Phi }_{\mathrm{mag}}^2.$$
(42)
The constant $`K`$ is given by $`K=\frac{dh^2}{\mu _0\mu _re^2a^2}`$, where $`a`$ is the two dimensional lattice constant, $`d`$ the distance between the copper oxide planes, and $`\mu _r`$ the relative permeability. As was discussed in Section III, this constant is huge compared to the typical electronic energies. The total energy as a function of the electromagnetic flux is then written as

$$E(\mathrm{\Phi }_{\mathrm{mag}})=E_{tJ}(\mathrm{\Phi }_{\mathrm{mag}})+E_{\mathrm{mag}}(\mathrm{\Phi }_{\mathrm{mag}}),$$
(43)
where $`E_{tJ}(\mathrm{\Phi }_{\mathrm{mag}})`$ denotes the energy of the $`tJ`$ model when there is an extra flux $`\mathrm{\Phi }_{\mathrm{mag}}`$ in addition to the spin-generated flux. Minimizing the energy with respect to the electromagnetic flux leads to an equation of the form $`g(\mathrm{\Phi }_{\mathrm{mag}})+K\mathrm{\Phi }_{\mathrm{mag}}=0`$, where $`g(\mathrm{\Phi }_{\mathrm{mag}})=E_{tJ}^{}(\mathrm{\Phi }_{\mathrm{mag}})`$. As we have pointed out, the magnetic energy-scale $`K`$ is much larger than the electronic energy which is of the order $`\mathrm{max}(t,J)`$. As a consequence of this, the magnetic flux will be suppressed and it is reasonable to put $`\mathrm{\Phi }_{\mathrm{mag}}=0`$ when we solve for the eigenstates of the system. For this to be a self-consistent solution there must be no currents in the system due to the spin-generated flux. This is true for the cases $`\mathrm{\Phi }=0`$ or $`\pi `$ which we have focused on.

## VII Striped structures

Striped structures forming antiphase domain walls between undoped antiferromagnetic regions have been experimentally observed in the doped high-temperature superconductors . We would like to understand if a striped phase can be explored using the effective $`tJ`$ Hamiltonian, Eq. (9). There are several facts that make this an appealing approach. We have already seen in the previous section that for low dopings, the uniform (spiral) states are unstable against phase-separation. Using our spin-polarized approach we can create a smooth antiphase boundary, successively changing the order-parameter from $`+\widehat{z}`$ to $`\widehat{z}`$. And, since all the holes of the stripe are located in the domain wall it could be favorable for the system to generate a flux in this region. This can be accomplished using a spin-texture, which simultaneously generates the antiphase boundary. Furthermore, the experimentally observed value of the doping of stripes in La<sub>2-x</sub>Sr<sub>x</sub>CuO<sub>4</sub> is $`0.5`$ holes per unit length of the stripe . This is close to the region where we have seen that a $`\pi `$-flux may be favorable. (See Section V.) Inspired by these nice properties we have undertaken an investigation of striped phases within our approach. Our main ambition has been to gain an understanding of what such a striped phase would look like, and, in particular, whether a fictitious flux is exploited or not in our model. The technique we use is a self-consistent Hartree calculation.

Let us consider a system with an antiphase boundary along the $`x`$-axis, and with the antiferromagnetic order being $`+\widehat{z}`$ at $`y=\mathrm{}`$ and $`\widehat{z}`$ at $`y=\mathrm{}`$. Two possible scenarios for an antiphase spin-order in a stripe comes to mind. First, the amplitude of the antiferromagnetic order-parameter may simply decrease, passing through zero and becoming negative as one passes through the stripe. This scenario preserves the rotational symmetry about the spin $`z`$-axis. The second possibility is that the Néel order-parameter starts to tilt as one approaches the stripe along the $`y`$-axis, lying within the spin $`x`$-$`y`$-plane at the center of the stripe and then rotating towards the positive $`z`$-axis as one goes to $`y=\mathrm{}`$. Introducing holes in the domain wall, the amplitude of the spin-order will decrease as it depends on the particle density through $`𝐒_𝐫=\frac{1}{2}n_𝐫\widehat{\mathrm{\Omega }}_𝐫`$, but it will not vanish. An illustration of these two scenarios is shown in Fig. 11. From an experimental point of view, Tranquada and coworkers argue  that there are indications speaking in favor of the first scenario, but these indications are not conclusive.

Four different stripe geometries with preserved uniaxial symmetry come to mind. First, stripes can go along either the $`(10)`$\- or $`(11)`$-directions and we can choose either site- or bond-centered stripes, all in all four possible combinations. The difference between site- and bond-centered stripes is shown in Fig. 12.

Our aim is to find the distribution of holes and spin-texture that minimizes the energy of an antiferromagnet, when we assume the holes are arranged in stripes. The configurations we consider have the structure shown in Fig. 13, consisting of a repeated structure of Néel-regions and antiphase domain walls. The order-parameter of the antiferromagnetic regions change sign every time a domain wall is passed.

We will assume that all the holes are located near the domain walls. Following the notation in Fig. 13, we will denote the width of a single antiphase boundary by $`\xi `$, and the width of an antiferromagnetic region by $`d`$. The average particle density in the domain wall will be denoted $`n_d`$, and we will also use the number of holes per unit length (along $`x`$) of the domain wall, $`\delta =(1n_d)\xi `$. Assuming an average filling $`n`$ of the system, we have $`(\xi +d)n=\xi n_d+d`$.

Before writing down the total energy of this configuration, we make the following observations. We note that the system in Fig. 13 is build up from units (antiferromagnets and domain walls) and we would like to write the total energy in terms of the energies of the individual units. To do this we define the energy per site of an antiferromagnetic unit, $`E_{AF}`$, as the energy per site of the antiferromagnet with periodic boundary conditions in the $`x`$\- and $`y`$-directions. The energy per site of the domain wall as a function of the number of holes per unit length of the stripe, $`E_d(\delta )`$, is similarly defined by putting periodic boundary conditions on the domain wall. In this case it is important that the edge of the domain wall is Néel-ordered and undoped. If this condition is not fulfilled, there will be surface energies associated with the gluing of a domain wall to an antiferromagnetic region. It is easy to see that minimizing the total energy of the system is equivalent to minimize the energy per introduced hole in the domain wall.

Our approach will therefore be to consider a single domain wall and minimize the energy per introduced hole with respect to the parameters describing the domain wall. These parameters contain spin-texture related parameters (which will be introduced shortly), the number of holes per unit length of the stripe ($`\delta `$), and finally we have the four options for the stripe geometry; site/bond-centered and direction (10)/(11).

Note that the definition of $`\xi `$ is somewhat arbitrary in the sense that there in practice may be a smooth crossover from the domain wall to the Néel-ordered region. For this reason, $`\delta `$ is a better measure of the stripe-doping than $`n_d`$ which depends on $`\xi `$. Knowing $`\delta `$ and $`\xi `$, the stripe periodicity is

$$l=2(\xi +d)=\frac{2\delta }{x},$$
(44)
where we have introduced the average doping of the system, $`x=1n`$. Thus we find that for low dopings, the separation between the stripes scales as $`x^1`$. This description is valid as long as the stripes remain separated so that we can neglect stripe-stripe interactions. This condition is fulfilled as long as $`x\frac{\delta }{\xi }`$. We note that as we change the overall doping of the system, the structure of the isolated stripes remains, at least as long as $`d1`$. Thus the wavevector describing the spin order is $`2\pi /l`$, and the wavevector of the charge order is twice that, i.e. $`4\pi /l`$.

### A Stripes in the $`(10)`$-direction

Let us start with a description of how a single $`(10)`$-stripe is modeled. Although we are interested in infinite domains, in the numerical simulations we are forced to work with finite stripe width, which we denote by the integer $`w`$. We will assume the local spin orientation of the stripe to described by a unit-vector field $`\widehat{\mathrm{\Omega }}_𝐫`$ as

$$\widehat{\mathrm{\Omega }}_𝐫=(1)^𝐫(\mathrm{sin}\theta _𝐫\mathrm{cos}\varphi _𝐫,\mathrm{sin}\theta _𝐫\mathrm{sin}\varphi _𝐫,\mathrm{cos}\theta _𝐫),$$
(45)
where the spherical angles $`\theta _𝐫`$ and $`\varphi _𝐫`$ are functions of position. The spatial dependence of these angles is assumed to take the form

$`\theta _𝐫`$ $`=`$ $`\mathrm{arccos}\left[\mathrm{tanh}{\displaystyle \frac{w+12y}{2\xi }}\right]`$ (46)
$`\varphi _𝐫`$ $`=`$ $`q_xx.`$ (47)
The construction can be thought of as a cylindrical projection (also known as Mercator projection) of the lattice into the spin sphere by associating latitude and longitude with the $`x`$\- and $`y`$-coordinates, respectively.

We are interested in the limit $`w\mathrm{}`$. This limit has two discrete values, depending on whether $`w`$ is even or odd. For odd values of $`w`$ the spin-order of Eq. (47) is such that the spins at $`y=(w+1)/2`$ are lying in the $`x`$-$`y`$-plane, and hence the stripe is site-centered. If we instead consider an even integer $`w`$, the center of the stripe is located in between two rows and the stripe is bond-centered.

In addition to the discrete choice of bond- or site-centered stripe, there are two free continuous parameters in Eq. (47); $`\xi `$, which determines the characteristic width of the spin-domain wall, and $`q_x`$, which is the pitch of the rotation of the spin about the $`z`$-axis along the length of the stripe. Fig. 14 illustrates this construction.

Making an apparently trivial point, which we return to in the next section, we note that there is no choice for $`\xi `$ and $`q_x`$ which unwraps the antiferromagnet into a ferromagnet; taking $`q_x=\pi `$ will create a ferromagnetic channel in the center of the domain wall, running along the stripe, but will not affect the $`z`$ component of the Néel order-parameter.

When writing down the Hamiltonian of this system we need to fix a gauge. If the spin-configuration around a square plaquette in terms of spherical angles is given by: $`(\theta _1,\varphi _1)`$, $`(\pi \theta _2,\varphi _1+\pi )`$, $`(\theta _2,\varphi _2)`$, and $`(\pi \theta _1,\varphi _2+\pi )`$, the product of the hopping elements around the plaquette is

$`\tau _{01}`$ $`\tau _{12}`$ $`\tau _{23}\tau _{30}=`$ (49)
$`\mathrm{sin}\theta _1\mathrm{sin}\theta _2\mathrm{sin}^2{\displaystyle \frac{\theta _1\theta _2}{2}}\mathrm{sin}^2{\displaystyle \frac{\varphi _1\varphi _2}{2}}0,`$
showing that the flux through a plaquette is exactly equal to $`\pi `$. The texture is therefore quite natural for the following reasons: It is periodic and has a uniform flux per square. The flux $`\pi `$ per plaquette favors heavily doped regions near the center of the stripe, while the effective hopping in the Néel regions, where the system is undoped, is vanishing. With flux $`\pm \pi `$ through the plaquettes there is no broken time-reversal symmetry, and there are no circulating currents or induced local magnetic fields.

In order to perform our calculations, we use a mixed representation using momentum space in the $`x`$-direction and real space in the $`y`$-direction due to the translational invariance of the system along the $`x`$-direction. The length of our system along the stripe will be denoted by $`L`$. Furthermore we assume that the number density at a site $`𝐫`$ only depends on the $`y`$-coordinate, i.e. $`n(𝐫)=n(y)`$. We will use the following Hartree-decomposition of the interaction term in Eq. (9)

$$n_𝐫n_𝐫^{}n(y)n_𝐫^{}+n_𝐫n(y^{})n(y)n(y^{}).$$
(50)
The effective Hamiltonian can then be written as

$`H_{\mathrm{bulk}}`$ $`=`$ $`2t{\displaystyle \underset{y=1}{\overset{w}{}}}(1)^y{\displaystyle \underset{k}{}}\mathrm{cos}(k)\tau _x(y)n_{k,y}t{\displaystyle \underset{y=1}{\overset{w1}{}}}\tau _y(y){\displaystyle \underset{k}{}}\left[c_{k,y+1}^{}c_{k,y}+\text{h.c.}\right]`$ (53)
$`+{\displaystyle \frac{J}{4}}{\displaystyle \underset{y=1}{\overset{w}{}}}(h_x(y)1){\displaystyle \underset{k}{}}\left[2n_{k,y}n(y)n(y)^2\right]`$
$`+{\displaystyle \frac{J}{4}}{\displaystyle \underset{y=1}{\overset{w1}{}}}(h_y(y)1){\displaystyle \underset{k}{}}[n_{k,y}n(y+1)+n_{k,y+1}n(y)n(y)n(y+1)],`$
where we have introduced the coordinate-dependent coupling constants

$`h_x(y)`$ $`=`$ $`\widehat{\mathrm{\Omega }}_𝐫\widehat{\mathrm{\Omega }}_{𝐫+\widehat{x}}`$ (54)
$`\tau _x(y)`$ $`=`$ $`\left({\displaystyle \frac{1+h_x(y)}{2}}\right)^{\frac{1}{2}},`$ (55)
and the analogous relations for $`h_y(y)`$ and $`\tau _y(y)`$. The factor $`(1)^y`$ in the hopping term of Eq. (53) comes from the flux $`\pm \pi `$ through each plaquette. In addition to the terms in Eq. (53), we have also added a local potential at the upper and lower edges of the system to simulate the effect of the adjoining antiferromagnetic region. Without this potential, which has the form

$$H_{\mathrm{b}.\mathrm{c}.}=\frac{J}{2}\underset{k}{}\left[n_{k,1}+n_{k,w}\right],$$
(56)
the system may gain energy by expelling the holes to the edges where they break fewer antiferromagnetic links. The total Hamiltonian is then given by

$$H_{10}=H_{\mathrm{bulk}}+H_{\mathrm{b}.\mathrm{c}.}.$$
(57)
The numerical calculation involves solving self-consistently for a charge profile described by $`n(y)`$, where

$$n(y)=L^1\underset{k}{}n_{k,y},$$
(58)
and the expectation value is with respect to the Fermi-distribution of quasiparticles of $`H_{10}`$. We use a chemical potential to control the overall number of particles in the system. The chemical potential is determined during each iteration of the self-consistency equation, Eq. (58). Our calculations are performed at temperatures close to zero, $`T=0.01t`$, and the quantity we focus on is the free energy as a function of the number of particles and the parameters describing the stripe, $`F(T,N)=H_{10}TS`$.

### B Stripes in the $`(11)`$-direction

The analysis of the $`(11)`$-stripes is similar, although the geometry of the stripe introduces some complications. In Fig. 15, we have shown the geometry used for these configurations. Due to the tilting of the lattice, we note that every second row of constant $`y`$ in the lattice is shifted to the right by half the lattice spacing in the $`x`$-direction. Moreover, if we denote the even- and odd sublattices by $`A`$ and $`B`$, respectively, we find that all points belonging to $`A`$ reside on points having odd values of $`y`$, while those belonging to $`B`$ are assigned even values of $`y`$. This is shown explicitly in Fig. 15.

The stripe is directed along the $`x`$-axis in the coordinate system defined in Fig. 15, and the spin-configuration is given by

$$\widehat{\mathrm{\Omega }}_𝐫=(1)^y(\mathrm{sin}\theta _𝐫\mathrm{cos}\varphi _𝐫,\mathrm{sin}\theta _𝐫\mathrm{sin}\varphi _𝐫,\mathrm{cos}\theta _𝐫).$$
(59)
We note that $`\theta _𝐫=0`$ and $`\theta _𝐫=\pi `$ correspond to the two Néel-states, and hence we can (as in the $`(10)`$-case) interpolate between the two by continuously changing $`\theta _𝐫`$ across the domain wall. To be explicit, we will use the following parametrization of the spherical angles:

$`\theta _𝐫`$ $`=`$ $`\mathrm{arccos}\left[\mathrm{tanh}{\displaystyle \frac{w+12y}{2\xi }}\right]`$ (60)
$`\varphi _𝐫`$ $`=`$ $`\{\begin{array}{cc}q_xx\hfill & y\text{ odd}\hfill \\ q_x(x+1/2)\hfill & y\text{ even}\hfill \end{array}.`$ (63)
Note that we have shifted $`x`$ by $`1/2`$ for even $`y`$-values to account for the shift of lattice points in the $`x`$-direction as discussed above. It is also important to stress that, contrary to the $`(10)`$-case, it is possible to recover a ferromagnetic configuration by a suitable choice of parameters, namely making $`\xi `$ large and taking $`q_x=2\pi `$. This corresponds to $`\theta _𝐫=\pi /2`$ and $`q_x=2\pi `$ for odd $`y`$’s and $`q_x=\pi `$ for even $`y^{}s`$, respectively. The effect of the rotation due to $`q_x`$ is therefore to rotate all spins belonging to sublattice $`B`$ by $`\pi `$ about the spin $`z`$-axis. The result is a ferromagnetic configuration, where all spins point along the positive $`x`$-axis. This difference between the $`(10)`$\- and $`(11)`$-stripes reflects the fact that the local field along either side of a $`(11)`$-stripe is ferromagnetic whereas in the $`(10)`$ case it is antiferromagnetic .

The next issue we will address is the properties of the fictitious fluxes generated by a certain spin-configuration in the domain wall. Fig. 16 shows the flux-pattern which is generated from the spin-configuration in Eq. (63).

We find that the flux is uniform along the $`x`$-direction, while it depends on the $`y`$-coordinate of the plaquette to which it belongs. Fig. 16 also fixes a gauge, defined by the phases $`\varphi _y`$. As we can see from the figure, there is one phase more than there are gauge invariant fluxes. Therefore, we can choose a gauge in which $`\varphi _1=0`$. Next, we consider the flux $`\mathrm{\Phi }_1`$ which is determined by $`\mathrm{\Phi }_1=\varphi _1+\varphi _2`$. More generally, we find $`\mathrm{\Phi }_y=(1)^{y+1}(\varphi _y+\varphi _{y+1})`$. In this way we can solve for the phases $`\varphi _y`$ in terms of the fluxes, finding

$$\varphi _{y+1}=(1)^y\underset{y^{}=1}{\overset{y}{}}\mathrm{\Phi }_y^{}.$$
(64)
Numerically, given the spin-configuration of the domain wall we calculate the fictitious fluxes, $`\mathrm{\Phi }_y`$, using Eq. (10) and then use Eq. (64) to find the appropriate phases that enter the Hamiltonian.

We define the amplitudes of the hopping and Heisenberg interactions as

$`h_y(y)`$ $`=`$ $`\widehat{\mathrm{\Omega }}_{x,y}\widehat{\mathrm{\Omega }}_{x,y+1}`$ (65)
$`h_d(y)`$ $`=`$ $`\widehat{\mathrm{\Omega }}_{x,y}\widehat{\mathrm{\Omega }}_{x+(1)^y,y+1}`$ (66)
$`\tau _y(y)`$ $`=`$ $`\left({\displaystyle \frac{1+h_y(y)}{2}}\right)^{\frac{1}{2}}`$ (67)
$`\tau _d(y)`$ $`=`$ $`\left({\displaystyle \frac{1+h_d(y)}{2}}\right)^{\frac{1}{2}},`$ (68)
which allows us to write the bulk part of the Hamiltonian as

$`H_{\mathrm{bulk}}`$ $`=`$ $`{\displaystyle \underset{y=1}{\overset{w1}{}}}{\displaystyle \underset{k}{}}\left[t\left(\tau _y(y)e^{i\varphi (y)}+\tau _d(y)e^{i(1)^yk}\right)c_{k,y+1}^{}c_{k,y}+\text{h.c.}\right]`$ (70)
$`+{\displaystyle \underset{y=1}{\overset{w1}{}}}{\displaystyle \frac{J}{4}}\left(h_y(y)+h_d(y)2\right){\displaystyle \underset{k}{}}\left[n_{k,y}n(y+1)+n_{k,y+1}n(y)n(y)n(y+1)\right].`$
For the same reasons as in the $`(10)`$-case, we will add a local potential to the vertical boundaries (note that each boundary site connects to two sites in the environment)

$$H_{\mathrm{b}.\mathrm{c}.}=J\underset{k}{}\left[n_{k,1}+n_{k,w}\right].$$
(71)
The total Hamiltonian, $`H_{11}`$, is the sum of the bulk- and boundary-contributions, i.e.

$$H_{11}=H_{\mathrm{bulk}}+H_{\mathrm{b}.\mathrm{c}.}.$$
(72)
The numerical procedures are completely equivalent to those used in the $`(10)`$-case.

### C The optimal stripe

According to our model, the physically relevant stripe configuration is that which minimizes the domain wall energy per introduced hole. We will use the undoped antiferromagnet as an energy reference state, this being the optimal state at zero doping. As we dope holes into the system, the total energy of the domain wall will be a discrete function of the geometry (site- or bond-centered and direction (10) or (11)) and a continuous function of $`N_h`$, $`q_x`$, and $`\xi `$. Note that since we work with a chemical potential, the number of holes in the domain wall, $`N_h`$, is not restricted to be an integer.

We label this energy $`E(N_h,q_x,\xi )`$. However, the physically interesting quantity is the number of holes per unit stripe length, $`\delta =N_h/L`$. We define the domain wall energy per hole according to

$$E_h(\delta ,q_x,\xi )=\frac{E(\delta L,q_x,\xi )E_{\mathrm{AF}}}{\delta L}.$$
(73)
As we argued in the beginning of this section, to find the optimal stripe-configuration we have to minimize this function with respect to $`\delta `$, $`q_x`$, and $`\xi `$ for bond- and site-centered $`(10)`$\- and $`(11)`$-stripes.

Turning to our numerical results, we have initially considered a system with a Heisenberg coupling $`J=0.40`$, where the energy scale is fixed by $`t=1`$. This value was chosen because it corresponds to a value of the exchange coupling constant which has been used by others to model the high-temperature superconductors. In Fig. 17 we show the optimal energy per hole as a function of $`\delta `$ for the four stripe geometries, i.e. we have plotted $`E_h(\delta )=\mathrm{min}_{q_x,\xi }E_h(\delta ,q_x,\xi )`$.

From Fig. 17 we can read off the optimal stripe configuration, which will make up the domain walls in the striped phase. As is indicated in this figure, the optimal domain wall is a site-centered stripe along the $`(10)`$-direction having $`\delta 0.46`$ holes per unit length of the stripe. This agrees with results from DMRG calculations by White and collaborators , who find stripes with $`\delta =0.5`$ for $`J=0.35`$. Furthermore, the experimental data indicate that $`\delta 0.5`$ . We also find that $`\xi `$ is very small for this optimal stripe, i.e. there is a sharp spin-domain wall with a single tilted row of spins. This row is ferromagnetically ordered as we find $`q_x=\pi `$, and the holes are tightly confined in the neighborhood of the domain wall. The spin- and charge-profiles are shown in Fig. 18.

Since the optimal domain wall is so narrow, the product of the effective hopping-amplitudes around any plaquette in the lattice is approximately zero. Thus we must conclude that the system does not take advantage of the fictitious $`\pi `$-flux through the plaquettes.

As Fig. 17 shows, the bond-centered $`(11)`$-stripe is energetically very close to the optimal $`(10)`$-stripe described above. An illustration of this domain wall configuration is shown in Fig. 19. This $`(11)`$-stripe is characterized by $`\delta 0.71`$, $`q_x=2\pi `$, and $`\xi 0.73`$. It is important to point out that $`q_x=2\pi `$ is not equivalent to $`q_x=0`$ since $`q_x=2\pi `$ performs a $`\pi `$-rotation about the $`z`$-axis of one of the sublattices, as we discussed in the previous subsection. Furthermore, we want to emphasize that this stripe-configuration does not induce any fictitious fluxes and consequently there are no currents that could energetically disfavor this configuration.

Since the optimal $`(10)`$-stripe is very close in energy to the $`(11)`$-stripes it is interesting to investigate what happens as we tune the strength of the Heisenberg interaction, $`J`$. Numerically we find that slightly increasing $`J`$ above $`J=0.40`$ favors the $`(10)`$-stripes compared to $`(11)`$-stripes, while decreasing $`J`$ favors the $`(11)`$-stripes. There will be a crossing point slightly below $`J=0.40`$, at approximately $`J0.36`$, where the $`(11)`$-stripes have lower domain wall energy than the $`(10)`$-stripes.

There is a technical point which is important when looking at the energy per hole as a function of $`\delta `$ for the $`(11)`$-configurations. As we mentioned in the previous subsection, it is possible to unwrap the $`(11)`$-domain walls into ferromagnets. If we follow the $`(11)`$-curves in Fig. 17 to larger values of $`\delta `$ we find that the energy successively decreases below what we called the optimal stripe configuration. The configurations that correspond to these low-energy states are close to ferromagnetic and with an almost uniform charge-distribution. Physically this corresponds to a global phase-separation into an undoped Néel region and a heavily doped ferromagnetic region, i.e. the phase-separation discussed in Section VI. Let us therefore compare the energetics of the optimal $`(10)`$-domain wall and the global phase-separation.

At the minimum of $`E_h(\delta )`$ in Fig. 17, we can read off the energy per hole of the optimal domain wall as $`E_{\mathrm{dw}}0.76`$. This is to be compared to the energy per hole for the totally phase-separated state $`E_{\mathrm{ps}}1.06`$ for the same value of $`J`$. This number was obtained from a Hartree calculation using Eq. (27) on a system of the same size as the stripe grid. This clearly shows that, within our approximation at least, phase-separation is energetically advantageous compared to domain wall formation. In a real system, the energy of the phase-separated state is raised due to the Coulomb interaction between the holes and this could make the domain wall thermodynamically stable.

The cuprates seem to favor the formation of rather narrow stripes, not phase-separation. We will take the point of view that there is some mechanism, not captured in our approach, such as a long-range Coulomb interaction, which prevents grouping all the holes together and instead favors the formation of stripes on some intermediate length scale. Therefore, we will only consider the stripe-configurations which are local in nature.

### D Including next-nearest-neighbor hopping

An unphysical feature of our simulation is that since hopping between the antiferromagnetic spins is forbidden, the effect of second-neighbor (diagonal) hopping becomes important in the Néel state. In our approximation of the $`tJ`$ model, hopping is frozen out for this state, so that electron transport will be dominated by any second-neighbor terms if they are non-zero. This term will permit electrons to diffuse into the Néel state and could therefore be expected to delocalize the holes from the stripe center. We will extend our model of the $`(10)`$-stripes by adding this hopping to the Hamiltonian through a term

$$H_{nnn}=t^{}\underset{\mathrm{𝐫𝐫}^{}}{}\left[\tau _{\mathrm{𝐫𝐫}^{}}^{}c_𝐫^{}^{}c_𝐫+\text{h.c.}\right],$$
(74)
where $`\mathrm{𝐫𝐫}^{}`$ denotes next-nearest-neighbor pairs. As before, the calculations reported below are performed with the Heisenberg coupling $`J=0.40`$ and at temperature $`T=0.01`$ in terms of energy units set by $`t=1`$. In this section, we will only consider the effect of second-neighbor hopping on the $`(10)`$-stripe, since this was found to be the optimal stripe-configuration for $`J=0.40`$.

It has been argued in the literature that the sign of $`t^{}`$ in the high-$`T_c`$ cuprates depends on whether the system is hole-doped or electron-doped . A hole-doped system corresponds to $`t^{}<0`$, while an electron-doped system has $`t^{}>0`$. In an antiferromagnet, the presence of a second-neighbor hopping is important since it allows for holes moving through the sublattices without disrupting the antiferromagnetic order. Typically, the value of $`t^{}`$ used in the literature for describing a hole-doped antiferromagnet is $`t^{}0.3`$.

If we consider the limit in which the nearest-neighbor hopping is completely frozen out, and there is only second-neighbor hopping, i.e.

$`H_{t^{}J}`$ $`=`$ $`t^{}{\displaystyle \underset{\mathrm{𝐫𝐫}^{}}{}}\left[\tau _{\mathrm{𝐫𝐫}^{}}^{}c_𝐫^{}^{}c_𝐫+\text{h.c.}\right]`$ (76)
$`+{\displaystyle \frac{J}{4}}{\displaystyle \underset{\mathrm{𝐫𝐫}^{}}{}}\left(\widehat{\mathrm{\Omega }}_𝐫\widehat{\mathrm{\Omega }}_𝐫^{}1\right)n_𝐫n_𝐫^{},`$
we note that the transformation $`c_𝐫(1)^xc_𝐫`$ leaves the Heisenberg term unchanged, while the second-neighbor term changes sign. This shows that in this limit the asymmetry between $`\pm t^{}`$ vanishes. Hence, the sign of $`t^{}`$ is irrelevant in the Néel-ordered regions, and it is only in the region where the spin-twist occurs that $`t`$ and $`t^{}`$ are simultaneously present and accordingly, the sign of $`t^{}`$ is important.

Introducing next-nearest-neighbor hoppings as in Section IV allows for new closed particle-orbits in the lattice, and hence also for new gauge invariant fluxes. There are four of these fluxes and they are defined in the right inset of Fig. 20.

Investigating the spin-structure of Eqs. (45) and (47), we find that the four fluxes $`\{\mathrm{\Phi }_k^y\}_{k=1}^4`$ associated with row $`y`$ are determined by a single parameter $`\mathrm{\Psi }_y`$ according to

$$\mathrm{\Phi }_1^y=\mathrm{\Psi }_y\text{ , }\mathrm{\Phi }_2^y=\pm \pi \mathrm{\Psi }_y\text{ , }\mathrm{\Phi }_3^y=\mathrm{\Psi }_y\text{ , }\mathrm{\Phi }_4^y=\pi +\mathrm{\Psi }_y.$$
(77)
Recalling our topological constraint, Eq. (15), we find that it is satisfied as

$$\mathrm{\Phi }_1^y+\mathrm{\Phi }_2^y\mathrm{\Phi }_3^y\mathrm{\Phi }_4^y=2\pi n,$$
(78)
with $`n=\pm 1`$, which is what we expect for an antiferromagnet. Which sign that applies to a certain plaquette depends on the sublattice associated with the plaquette as well as the spin-configuration. Fig. 20 defines a gauge by introducing the phases $`\alpha _y`$ and $`\beta _y`$. Using Eq. (77) it is easy to read off the parameters $`\alpha _y`$ and $`\beta _y`$ from Fig. 20. Doing this we find $`\alpha _y=\beta _y=\mathrm{\Psi }_y`$ for odd values of $`y`$, and $`\alpha _y=\beta _y=\pm \pi \mathrm{\Psi }_y`$ for even $`y`$’s.

As in the nearest-neighbor case, we will work in momentum space in the $`x`$-direction, while keeping the real space description in the $`y`$-direction. We note from Fig. 20 that the phases of the links form a staggered structure, doubling the size of the unit cell, and introducing a scattering between states of momenta $`k`$ and $`k+\pi `$. To deal with this we introduce the same two-component wavefunctions as was used in Section VI.

Using this model we have investigated how the optimal stripe evolves as the second-neighbor hopping amplitude $`t^{}`$ is changed from zero. For small values of $`t^{}`$, approximately, $`0.3<t^{}<0.15`$, the structure of the stripe is largely unchanged. It is still described by $`q_x=\pi `$ and vanishing $`\xi `$. The optimal number of holes per unit stripe length also remains ($`\delta 0.46`$). All that happens is basically that there are small redistributions of the holes within the stripe.

Concerning the tendency to global phase-separation we have considered the behavior of $`E_{\mathrm{dw}}(t^{})/E_{\mathrm{ps}}(t^{})`$, i.e. the ratio of the energy per hole of the domain wall and phase-separated states, respectively. We find that this ratio decreases as we increase $`t^{}`$ slightly from zero, and decreases when $`t^{}`$ becomes negative. This indicates that a negative $`t^{}`$ favors the domain wall configurations compared to the global phase-separation but, at least for small $`t^{}`$, the domain walls are still unstable against phase-separation.

When $`t^{}`$ becomes larger than $`0.15`$, the holes diffuse into the antiferromagnet and widens the domain wall. The optimal number of holes per unit stripe length increases. The spin-twist $`\xi `$ becomes non-zero making the antiphase boundary wider. Since we are working with a finite width $`w`$ of the domain wall, we get problems when the holes start diffusing away from the center of the domain wall. For our model to be valid, we must require that the domain wall is Néel-ordered and undoped at its vertical edges. This is to avoid surface energies when gluing together the domain wall with a Néel region.

In the case relevant for the hole-doped cuprate planes, i.e. $`t^{}<0`$, we find that decreasing the value of $`t^{}`$ below $`0.3`$ keeps the structure of the stripe rather intact in the sense that $`\xi `$ remains vanishingly small and that the holes are localized close to the antiphase boundary. The optimal number of holes per unit stripe length does however change, it is reduced as $`t^{}`$ is decreased, e.g. at $`t^{}=0.5`$ we find that the energy per hole is minimized by $`\delta 0.38`$. However, if we further decrease $`t^{}`$ the optimal doping $`\delta `$ will rise again as we reach $`t^{}0.7`$. For such large negative values of $`t^{}`$, the holes will spread into the antiferromagnet just as we found in the positive $`t^{}`$ case. Our numerical calculations indicate that $`\xi `$ remains small, i.e. the spin-twist still occurs over a very short distance.

## VIII Conclusion

We have considered an effective version of the two dimensional $`tJ`$-model where the electrons are considered as completely spin-polarized. The effect of such a spin texture is to generate a fictitious topological flux through the lattice. In this paper we have extended the discussion of a previous paper  concerning the properties of these fluxes.

Keeping in mind the result of Hasegawa et al. , where it is shown that the energy of a free electron gas on a square lattice is minimized when there is one flux quanta per particle, the possibility of the system prefering a flux-generating spin-configuration does not seem remote.

To check the above theory we have performed Hartree-Fock mean-field theory calculations for the system. In the Hartree approximation it seems like the system prefers to form a flux-phase for certain choices of doping and coupling constants. However when the exchange-terms are included this effect seems to vanish and the coplanar spiral phase is energetically more favorable than the flux phase. Introducing a nearest-neighbor Coulomb repulsion it is possible to make the flux-phase energetically most favorable also with the exchange terms present. However, the calculations indicate that for a wide range of dopings, these uniform phases are unstable against phase-separation into an undoped antiferromagnet and a highly doped coplanar spiral phase. Thus we have to conclude that from a thermodynamic point of view, we do not expect to find a flux-phase in the phase-diagram of the model we have considered.

The main part of the paper has been concerned with the generalization of this construction to describe stripes, directed along either the (10)- or (11)-direction of the lattice. These stripes are appealing as they provide a smooth realization of an antiphase domain wall, continuously merging two Néel-ordered regions with opposite signs of the order-parameter. The holes naturally reside within this domain-wall which generates a fictitious flux which further can reduce the energy of the holes in domain wall.

Using $`J=0.40`$, we find that the optimal antiphase domain wall is site-centered and runs along the $`(10)`$-direction of the lattice. The structure of this domain wall is such that there is a sharp spin-twist, basically only affecting a single row of spins, which aligns the spins ferromagnetically in a one-dimensional channel. The doping of the domain wall is approximately $`\delta =1/2`$ holes per unit length of the stripe. All holes are tightly bound to the domain wall, spreading over three rows of lattice sites. Due to the narrow spin-structure of this stripe, we found that the system does not exploit the $`\pi `$-flux generated through each plaquette by this domain wall configuration.

We also find that stripes directed along the $`(11)`$-direction are energetically very close the optimal $`(10)`$-stripe, and if we decrease $`J`$ they will become the most advantageous domain walls at $`J0.36`$, while increasing $`J`$ favors the $`(10)`$-stripe. Looking at the structure of the optimal $`(11)`$-stripes we found that they do not generate any fictitious flux.

Comparing the energetics of the domain walls with global phase-separation, we find that within our approximation the global phase-separation is favorable. Furthermore, we have incorporated second-neighbor hopping in the case of $`(10)`$-stripes. For small values of this hopping, the structure of the optimal domain wall remains, while at larger values the holes starts spreading out, widening the domain wall.

There are a number of questions that are left unanswered at this point, and which we believe are interesting to further investigate. Concerning the next-nearest-neighbor hopping it would be interesting to investigate more through how it affects the tendency towards global phase-separation, and if it can stabilize the stripes. Furthermore, it would be interesting to further examine the effect of this hopping on the structure of the stripes, and also investigate its effects on the $`(11)`$-stripes.