# Antiferromagnetism in doped anisotropic two-dimensional spin-Peierls systems

\[

## Abstract

We study the formation of antiferromagnetic correlations induced by impurity doping in anisotropic two-dimensional spin-Peierls systems. Using a mean-field approximation to deal with the inter-chain magnetic coupling, the intra-chain correlations are treated exactly by numerical techniques. The magnetic coupling between impurities is computed for both adiabatic and dynamical lattices and is shown to have an alternating sign as a function of the impurity-impurity distance, hence suppressing magnetic frustration. An effective model based on our numerical results supports the coexistence of antiferromagnetism and dimerization in this system.

PACS: 75.10 Jm, 75.40.Mg, 75.50.Ee, 64.70.Kb

\]

General interest for spin-Peierls (SP) systems was recently renewed by the discovery of CuGeO<sub>3</sub> , the first inorganic SP material. The SP transition is characterized by a freezing of the spin fluctuations below an energy scale given by the spin gap $`\mathrm{\Delta }_S`$ accompanied by a simultaneous lattice dimerization . Rich phase diagrams have been obtained experimentally upon doping this compound with non-magnetic impurities . In site-substituted systems such as (Cu<sub>1-x</sub>M<sub>x</sub>)GeO<sub>3</sub>, where M=Zn (Ref. ) or Mg (Ref. ), long range antiferromagnetic (AF) order is stabilized at low temperature while the dimerization still persists (D-AF phase). In Mg-doped compounds, for impurity concentrations larger than a critical value ($`x_c0.02`$), a first order transition occurs between the D-AF phase and a uniform AF (U-AF) phase where the dimerization disappears.  The coexistence of the two types of order in the D-AF phase is an intriguing phenomenon since lattice dimerization favors the formation of spin singlets on the bonds while low energy spin fluctuations exist in an AF phase.

Theoretically, the effect of impurity doping in SP systems was considered for fixed-dimerized , adiabatic  and quantum-dynamical  lattices. A single nonmagnetic impurity releases a soliton in the chain which can be viewed as a kink in the lattice distortion. In the absence of interchain couplings, such an excitation can freely propagate away from the impurity. On the other hand, the interchain elastic coupling $`K_{}`$ was shown to produce confinement within some distance from the impurity .

For a finite impurity concentration, the coexistence between SP and AF orders has been previously discussed either considering randomly distributed domain walls in a $`XX`$ chain  or assuming small fluctuations of the magnetic exchange constants . Despite their success to describe some experimental results, these models are rather limited since they do not take into account the microscopic origin of the soliton formation nor the interchain couplings. In this paper, a realistic microscopic model with interchain magnetic and elastic couplings is considered to describe the formation of a region with AF correlations in the vicinity of each impurity and which allows an estimation of the effective interaction between impurities in the two-dimensional (2D) system. Thus, we are able to construct and study an effective model in order to understand the effects of a finite impurity concentration. 

In a first step, the spin-phonon coupling is treated in the adiabatic approximation. The Hamiltonian $`=_{\mathrm{mag}}+_{\mathrm{el}}`$ is:

$`_{\mathrm{mag}}`$ $`=`$ $`J_{}{\displaystyle \underset{i,a}{}}(1+\delta _{i,a})𝐒_{i,a}𝐒_{i+1,a}+J_{}{\displaystyle \underset{i,a,b}{}}𝐒_{i,a}𝐒_{i,b},`$ (1)   
$`_{\mathrm{el}}`$ $`=`$ $`{\displaystyle \underset{i,a}{}}\{{\displaystyle \frac{1}{2}}K_{}\delta _{i,a}^2+K_{}\delta _{i,a}\delta _{i,a+1}\},`$ (2)   
where $`a`$ is a chain index and $`i`$ labels the sites along the chains. Atomic displacements are only considered along the chain direction, $`\delta _{i,a}`$ being here a classical variable related to the change of the bond length between sites $`(i,a)`$ and $`(i+1,a)`$. The magnetic part includes a magnetoelastic coupling $`J_{}`$ (hereafter set to unity) and an exchange interaction $`J_{}`$ connecting nearest neighbor (NN) chains. We eventually include in our model a next NN exchange interaction along the chain whose relevance for CuGeO<sub>3</sub> has been emphasized . $`_{\mathrm{el}}`$ is the elastic energy. The interchain elastic interaction ($`K_{}`$) is limited to NN chains. Stability of the lattice implies $`K_{}2|K_{}|`$. Typical values of the parameters for CuGeO<sub>3</sub> are $`J_{}0.1`$ (Ref. ) and $`K_{}/K_{}0.2`$ (Ref. ).

In order to study numerically model (2), we treat exactly the single chain problem using exact diagonalization (ED) or Quantum Monte Carlo (QMC) methods, while the interchain magnetic coupling is treated in a self-consistent mean-field (MF) approximation. This is an standard procedure to include interchain couplings in the study of quasi-one-dimensional systems. Moreover, Inagaki and Fukuyama have used a similar MF approximation to treat the interchain coupling in the bosonized version of (2) within a self-consistent harmonic approximation. Thus, in our procedure, the interchain magnetic coupling is replaced by its MF form:

$`_{\mathrm{MF}}^{}=J_{}{\displaystyle \underset{i,a,b}{}}\{S_{i,a}^zS_{i,b}^z+S_{i,a}^zS_{i,b}^zS_{i,a}^zS_{i,b}^z\}.`$ (3)   
By extending a similar approach previously applied to one-dimensional (1D) chains  to the case of the 2D lattice, a sweep is performed in the transverse direction, i.e. $`aa+1`$. For each chain $`a`$, we compute the MF values $`S_{i,a}^z`$ and the classical variables $`\{\delta _{i,a}\}`$ by energy minimization, which is achieved by solving iteratively the equations

$`\delta _{i,a}=\{J_{}𝐒_{i,a}𝐒_{i+1,a}+K_{}(\delta _{i,a+1}+\delta _{i,a1})\}/K_{}.`$ (4)   
Then, these new values of the AF and SP order parameters enter as input for the chain $`a+1`$. This procedure is iterated until convergence is reached. In this way, we can study numerically finite clusters consisting of $`N`$ coupled chains with $`L`$ sites where, typically, $`N\times L=12\times 18`$ in ED and $`N\times L=6\times 40`$ in QMC, with toroidal boundary conditions.

A similar MF approach can be adapted to study a model equivalent to (2) but with quantum phonon degrees of freedom. In this case, phonon operators $`b_{i,a}^{}`$ and $`b_{i,a}`$ are introduced on each bond and the displacements $`\delta _{i,a}`$ become $`g(b_{i,a}^{}+b_{i,a})`$, where $`g`$ is the magnetoelastic constant. Then, the classical elastic term $`_{\mathrm{el}}`$ is replaced by its quantum version,

$`_{\mathrm{ph}}=\mathrm{\Omega }{\displaystyle \underset{i,a}{}}\{b_{i,a}^{}b_{i,a}^{}+\mathrm{\Gamma }(b_{i,a}^{}+b_{i,a}^{})(b_{i,a+1}^{}+b_{i,a+1}^{})\},`$ (5)   
where $`\mathrm{\Gamma }=K_{}/(2K_{})`$, and the phonon frequency $`\mathrm{\Omega }`$ is related to $`g`$ by $`\mathrm{\Omega }=2g^2K_{}`$. The adiabatic limit (2) is recovered when $`\mathrm{\Omega }0`$ (requiring $`g0`$ also). Similarly to the interchain magnetic term, the interchain elastic term of (5) is then treated in mean-field by introducing a lattice order parameter $`\delta _{i,a}^{\mathrm{MF}}=gb_{i,a}^{}+b_{i,a}`$. Then, the term (5) is replaced by:

$`_{\mathrm{ph},\mathrm{MF}}=\mathrm{\Omega }{\displaystyle \underset{i,a}{}}\{b_{i,a}^{}b_{i,a}^{}+{\displaystyle \frac{\mathrm{\Gamma }}{g}}(b_{i,a}^{}+b_{i,a}^{})\delta _{i,a+1}^{MF}\},`$ (6)   
Note that in this case it is not necessary to solve an equation similar to (4). To diagonalize the single chain Hamiltonian with $`L8`$, a Lanczos algorithm is used. The phononic degrees of freedom are treated within a variational formalism previously introduced . Note that inelastic neutron scattering experiments  on CuGeO<sub>3</sub> reveal a rather large phonon frequency $`\mathrm{\Omega }/J_{}2`$ suggesting large lattice quantum effects in this material.

As a preliminary study, we apply the MF procedure to the case of an homogeneous system without impurities. In models (2) and (5), $`J_{}`$ is expected to stabilize the AF state while a small $`K_{}`$ (or large magnetoelastic coupling $`g`$) tends to favour SP order. For each value of the couplings $`J_{}`$ and $`K`$ (where $`K=K_{}2|K_{}|`$ is the relevant parameter in the SP phase) we obtain the ground state without imposing any restriction on the MF parameters. We found only two different phases, the SP phase where $`S_z=0`$ and $`\delta _{i,a}0`$ and the antiferromagnetic (AF) phase with $`S_z0`$ and $`\delta _{i,a}=0`$. Then, the phase diagram in the $`KJ_{}`$ plane, can be obtained in a more efficient way by a direct comparison of the energies of the Néel antiferromagnetic phase, and of the uniformly dimerized phase. The phase diagram shown in Fig. 1 exhibits a transition line between AF and SP phases. In the adiabatic case, this line could be fitted by a law $`J_{}=\frac{A}{K}+B`$ with $`A=0.3656`$ and $`B=0.06`$ (this artificial small negative value may be a consequence of small finite size effects, see e.g., Ref. ). A phase boundary with a form $`J_{}=\frac{A}{K}`$ has been predicted by Inagaki and Fukuyama . However their bosonized approach do not fix unambiguously the value of $`A`$.

In the case of the adiabatic calculation finite size effects were shown to be small for $`K<2`$. On the other hand, for the dynamical lattice, the calculation is reliable only for larger lattice couplings (i.e. smaller values of $`K`$) due to stronger finite size effects. As expected, for very small $`K`$, lattice quantum fluctuations are less effective in dimerizing the chain than the adiabatic lattice. Then, the phase boundary obtained with quantum phonons is located below the adiabatic one. This tendency becomes clear as $`\mathrm{\Omega }`$ is increased, as shown in the figure. On the other hand, it has been suggested that dynamical phonons induce an effective magnetic frustration. This frustration which becomes relatively important for larger $`K`$ destabilizes the AF phase thus moving the phase boundary upwards as seen in the figure. Consistently with this behaviour, if a next NN exchange term is included in the Hamiltonian, in the adiabatic approximation, the SP phase is more stable and the phase boundary is located above the corresponding curve for $`\alpha =0`$. This behavior is shown in Fig. 1 for the realistic value of $`\alpha =0.36`$ obtained for CuGeO<sub>3</sub> (Ref. ), where $`\alpha `$ is the value of the next NN exchange coupling constant in units of $`J_{}`$. We have checked that the set of realistic parameters $`K20`$ and $`K_{}`$, $`J_{}`$ mentioned above corresponds to a point in the SP phase.

To start our analysis of impurity doping, we consider a single impurity in order to investigate the appearance of AF correlations in the SP phase. As mentioned above, the impurity releases in the chain a topological spin-1/2 solitonic excitation  characterized by a change of parity of the dimerization order which occurs in a finite region of longitudinal size $`\xi _{}`$ given by the soliton width. The local magnetization on each chain $`a`$ can be decomposed into uniform and staggered components , $`S_{i,a}^z=M_{i,a}^{\mathrm{unif}}+(1)^{i+a}M_{i,a}^{\mathrm{stag}}`$. In fact, the excess uniform component $`S_{\mathrm{sol}}^z=\pm \frac{1}{2}`$ and the soliton, characterized by a broad maximum of $`M_{i,a}^{\mathrm{stag}}`$, remain confined in the chain with the impurity. However, as seen in Fig. 2(a),

the interchain magnetic coupling $`J_{}`$ generates a large staggered component with the same parity, i.e. $`M_{i,a}^{\mathrm{stag}}`$ keeping the same sign, in the neighboring chains. Simultaneously, the amplitude of the SP dimerization is significantly suppressed compared to the bulk value i.e. far away from the impurity. Large AF correlations can be seen up to more than four chains away from the impurity chain for magnetic couplings as small as $`J_{}=0.1`$, in particular in the vicinity of the SP $``$ AF transition line of Fig. 1. The transverse range of the AF ‘polarization cloud’ around the impurity increases strongly with the transverse coupling $`J_{}`$.

A crucial feature of the polarization surrounding the impurity-soliton area is that the sign of $`M_{i,a}^{\mathrm{stag}}`$ is unambiguously fixed by the orientation ($`S_{\mathrm{sol}}^z=\pm \frac{1}{2}`$) of the soliton and by the position $`(i_0,a_0)`$ of the impurity in such a way that $`\mathrm{sign}\{\mathrm{M}_{\mathrm{i},\mathrm{a}}^{\mathrm{stag}}\}=\mathrm{sign}\{\mathrm{S}_{\mathrm{sol}}^\mathrm{z}\}(1)^{\mathrm{a}_0+\mathrm{i}_0+1}`$. This fact can be simply understood in the strong dimerization limit ($`\delta 1`$) where the introduction of the impurity on a given site releases a spin-1/2 on one of its neighboring sites by breaking a singlet bond. For smaller lattice coupling, the excess spin can effectively hop from site to site on the same sublattice (due to the underlying dimerization), hence producing AF correlations with the parity defined above.

Let us now consider two impurities introduced simultaneously on two sites $`(i_1,a_1)`$ and $`(i_2,a_2)`$ of different chains ($`a_1a_2`$). When the two polarization clouds associated to each soliton-impurity ‘pair’ start to overlap, one expects their interaction to depend on the relative orientation of the two solitons. As seen in Figs. 2(c,d), quite different patterns correspond to the singlet and triplet arrangements of the two spin-1/2 solitons. As confirmed by our calculations, the lowest energy is always obtained for a spin state which leads to the same parity of the staggered magnetization associated to each impurity, i.e. which avoids completely magnetic frustration. The simple argument developed above for a single impurity then suggests that a triplet $`S=1`$ (singlet $`S=0`$) configuration is favored when the two impurities are located on the same sublattice (opposite sublattices). It is then appropriate to define an effective magnetic coupling between the AF clouds associated to each impurity by $`J_{\mathrm{eff}}=E_{\mathrm{S}=1}E_{\mathrm{S}=0}`$. For a wide range of parameters leading to a SP state in the bulk, we have numerically found that $`J_{\mathrm{eff}}`$ is ferromagnetic (F) if the two impurities belong to the same sublattice and antiferromagnetic in the opposite case. This implies that the coupling between the two local moments associated to the impurities is either F or AF in such a way that no frustration occurs. The magnetic coupling, for physical values of the parameters, can be fairly extended in space as seen in Fig. 3(a). Its range is directly controlled by the overlap of the polarization clouds. It follows roughly a behavior like

$`J_{\mathrm{eff}}J_0(1)^{\mathrm{\Delta }a+\mathrm{\Delta }i+1}\mathrm{exp}(C_{}{\displaystyle \frac{\mathrm{\Delta }a}{\xi _{}}})\mathrm{exp}(C_{}{\displaystyle \frac{\mathrm{\Delta }i}{\xi _{}}}),`$ (7)   
where $`J_0J_{}`$, $`C_\alpha `$ are of order unity and $`\xi _\alpha J_\alpha /\mathrm{\Delta }_S`$.

To get an insight on how finite size effects might affect our results we have increased our cluster size ($`12\times 18`$) in both directions. The change of the transversal dimension has very little effect because the polarization clouds are almost independent when the impurities are separated by more than four chains (see the very small values of $`J_{\mathrm{eff}}`$ for $`\mathrm{\Delta }a=4`$ in Fig. 3(a)). On the other hand, by increasing the length of each chain, some changes occur but only when the impurities are located at the largest distances. However, the exponential decay of the effective interactions with distance (see below) is not qualitative changed as it can be seen in 3(a)). Therefore, we expect that the numerical values of the fitting parameters would not be much affected by finite size effects leaving the overall behaviour essentially unchanged. In particular, we believe that the presence of long range AF order in the effective model (see discussion below) is a robust feature not affected by finite size effects.

When two impurities are introduced on the same chain ($`a_1=a_2`$) two cases have to be distinguished. If the impurities are located on the same sublattice a similar behavior is observed as described above (compare Fig. 2(b) and Fig. 2(d)), i.e the effective interaction is ferromagnetic. However, the magnitude of $`|J_{eff}|`$ is $`0.4`$ (very slowly decaying as $`\mathrm{\Delta }i`$ increases) for the parameters of Fig. 3(a), i.e. much larger than the values corresponding to impurities in different chains. If the impurities belong to different sublattices then a chain with even number of sites is cut into two segments with even number of sites each. In the lowest energy configuration ($`S=0`$) no soliton-antisoliton pair was observed for separations $`\mathrm{\Delta }i`$ up to $`20`$ in agreement with previous work  and the triplet excitation energy remains large ($`\mathrm{\Delta }_S`$). Then, one can expect that for larger chains, when the formation of soliton becomes favourable, the effective interaction between them will be AF and their magnitude will be of the order of $`\mathrm{\Delta }_S`$. In summary, one can assume that the effective interaction between impurities on the same chain has the same form of Eq. (7) with $`\mathrm{\Delta }a=0`$ and $`J_0`$ being now $`\mathrm{\Delta }_S`$. This form is similar to the one adopted in Ref.  for impurities in a single chain except that these authors do not consider the sublattice sign alternation. Nevertheless, it should be noticed that this AF interaction should also decay for very large $`\mathrm{\Delta }i`$.

When lattice quantum fluctuations are introduced ($`\mathrm{\Omega }>0`$), the qualitative properties of the effective interaction $`J_{\mathrm{eff}}`$ are preserved. Consistently with the relatively larger stability of the AF phase in the small $`K`$ region shown in Fig. 1 with respect to the adiabatic case, we have found that lattice dynamics lead to an increase of the size of the AF cloud associated to each soliton. Therefore, the magnitude of the magnetic coupling $`J_{\mathrm{eff}}`$ increases with the phonon frequency $`\mathrm{\Omega }`$ as shown in Fig. 3(b).

The final part of our study is the analysis of a simple effective two-dimensional spin-1/2 Heisenberg model between impurities with a long range interaction given by Eq. (7). The ‘bare’ parameters are the same as in Fig. 3 and the parameters of the expression (7) have been obtained by fitting the curves shown in that figure and similar data for the case of $`\mathrm{\Delta }a=0`$. In the direction perpendicular to the chains we have neglected the effective interactions beyond a distance $`\mathrm{\Delta }a=5`$. We have also assumed that even segments are associated to a soliton-antisoliton pair.

A given number of spin-1/2 impurities $`4N_{\mathrm{imp}}16`$ is thrown at random on systems of coupled chains of sizes up to $`40\times 40`$. Then, the staggered magnetization $`M_{\mathrm{stag}}=(1/N_s^2)(_{i,a}(1)^{i+a}S_{i,a}^z)^2`$, where $`N_s=N\times L`$, is computed and averaged over, typically, 12,000-16,000 random realizations. The square root of this quantity is shown in Fig. 4. By extrapolating to the bulk limit for a fixed impurity doping using a polynomial in (1/$`\sqrt{N_s}`$) we found that $`M_{\mathrm{stag}}`$ is finite, implying long range AF order, and slowly decreasing as $`x`$ goes to zero. This behaviour is consistent with experimental results suggesting that $`M_{\mathrm{stag}}`$ decays exponentially to zero as $`x0`$.

In conclusion, spin-1/2 solitons released in 2D anisotropic SP systems by the introduction of impurities were shown to experience spatially extended F or AF exchange interactions depending on their relative positions. These exchange interactions coexisting with the SP order are calculated from realistic microscopic models and used to construct a simple effective model which in turn enables us to show the establishment of long range AF order and to compute the AF order parameter as a function of the impurity doping.

This work was supported in part by the ECOS-SECyT A97E05 program. We thank IDRIS (Orsay, France) and Florida State University for using their supercomputer facilities.