The origin of up-up-down-down magnetic order in Cu$_2$GeO$_4$

TL;DR

This study uses density-functional calculations to explain the unique up-up-down-down magnetic order in Cu$_2$GeO$_4$, highlighting the role of nearly canceled exchange interactions and structural buckling in its magnetic behavior.

Contribution

It reveals that weak symmetric anisotropy and structural buckling are key to the magnetic order, differing from previous assumptions about exchange interactions.

Findings

Negligibly small nearest-neighbor coupling J_1 due to exchange compensation.

Weak anisotropy induces the UUDD magnetic order.

Buckling of copper chains suppresses J_1, distinguishing Cu$_2$GeO$_4$ from similar magnets.

Abstract

We use density-functional band-structure calculations to explore the origin of the up-up-down-down (UUDD) magnetic order in Cu$_2$GeO$_4$ with the frustrated $J_1-J_2$ spin chains coupled into layers within the spinel-like crystal structure. In contrast to earlier studies, we find that the nearest-neighbor coupling $J_1$ may be negligibly small, owing to a nearly perfect compensation of the ferromagnetic direct exchange and antiferromagnetic superexchange. Under this condition, weak symmetric anisotropy of the exchange couplings gives rise to the UUDD order observed experimentally and also elucidates the non-trivial ordering pattern between the layers, whereas a small Dzyaloshinsky-Moriya interaction causes a spin canting that may generate local electric polarization. We argue that the buckling of the copper chains plays a crucial role in the suppression of $J_1$ in Cu$_2$GeO$_4$ and sets this compound apart from other $J_1-J_2$ chain magnets.

Figures

Tables

Table 1. Table 1: Isotropic exchange couplings (in meV) in Cu 2 GeO 4 calculated within DFT+ U 𝑈 U ( U d = 9.5 subscript 𝑈 𝑑 9.5 U_{d}=9.5  eV, J H = 1 subscript 𝐽 𝐻 1 J_{H}=1  eV, FLL) using three different band-structure codes: FPLO, ELK and VASP.

	FPLO	ELK	VASP
$J_1$	$-7.2$	$-3.3$	$-0.2$
$J_2$	6.0	5.0	5.6
$J$	9.0	7.8	8.5
$J_{ab}$	0.4	0.3	0.4
$J_c$	$-0.1$	0.5	$-0.1$

Table 2. Table 2: DFT and DFT+ U 𝑈 U ( U d = 9.5 subscript 𝑈 𝑑 9.5 U_{d}=9.5  eV, J d = 1 subscript 𝐽 𝑑 1 J_{d}=1  eV, FLL) relative energies E 𝐸 E of magnetic configurations (in meV) used for evaluating J 1 = ( E ↑ ↑ + E ↓ ↓ − E ↓ ↑ − E ↑ ↓ ) / 4 ​ S 2 subscript 𝐽 1 subscript 𝐸 ↑ absent ↑ subscript 𝐸 ↓ absent ↓ subscript 𝐸 ↓ absent ↑ subscript 𝐸 ↑ absent ↓ 4 superscript 𝑆 2 J_{1}=(E_{\uparrow\uparrow}+E_{\downarrow\downarrow}-E_{\downarrow\uparrow}-E_{\uparrow\downarrow})/4S^{2} ( S = 1 / 2 𝑆 1 2 S=1/2 ) within three different band-structure codes. The ↑ ↑ \uparrow and ↓ ↓ \downarrow symbols denote the magnetic moment alignment of two nearest-neighbour Cu 2+ ions coupled by the J 1 subscript 𝐽 1 J_{1} interaction, whereas magnetic moments of other Cu 2+ ions were fixed according to the procedure described in Ref.  Xiang  et al. ( 2011 ) . Magnetic moments of the two interacting copper sites, m 1 subscript 𝑚 1 m_{1} and m 2 subscript 𝑚 2 m_{2} (in μ B subscript 𝜇 𝐵 \mu_{B} ), are given for comparison. They are essentially similar between different spin configurations in DFT+ U 𝑈 U , but deviate from each other in spin-polarized DFT that underestimates local moments in correlated materials.

		DFT			DFT+$U$
		FPLO	ELK	VASP	FPLO	ELK	VASP
$\uparrow \uparrow$	$E$	166.2	165.7	156.1	21.2	20.2	28.2
	$m_1$	0.65	0.62	0.55	0.80	0.80	0.80
	$m_2$	0.65	0.62	0.55	0.80	0.80	0.80
$\downarrow \downarrow$	$E$	0	0	0	0	0	0
	$m_1$	-0.57	-0.56	-0.52	-0.79	-0.79	-0.79
	$m_2$	-0.57	-0.56	-0.52	-0.79	-0.79	-0.79
$\uparrow \downarrow$	$E$	61.2	62.5	56.8	14.1	11.6	14.1
	$m_1$	0.64	0.62	0.56	0.80	0.80	0.80
	$m_2$	-0.55	-0.54	-0.49	-0.79	-0.79	-0.79
$\downarrow \uparrow$	$E$	63.7	64.7	58.9	14.3	11.9	14.3
	$m_1$	-0.55	-0.54	-0.49	-0.79	-0.79	-0.79
	$m_2$	0.64	0.62	0.56	0.80	0.80	0.80
$J_1$		41.3	38.5	40.4	$-7.2$	$-3.3$	$-0.2$

Table 3. Table 3: The magnetic model parameters for Cu 2 GeO 4 , Li 2 CuO 2 , and CuGeO 3 in the basis of Wannier functions (in meV). The corresponding Coulomb parameters calculated by using random phase approximation (RPA). In particular, self-screening effects were extracted from Coulomb integrals (constrained RPA), whereas the direct exchange integrals were evaluated without excluding the self-screening effects. We use the notation U ~ i ​ j = U i ​ i − U i ​ j subscript ~ 𝑈 𝑖 𝑗 subscript 𝑈 𝑖 𝑖 subscript 𝑈 𝑖 𝑗 \widetilde{U}_{ij}=U_{ii}-U_{ij} for the effective Coulomb parameter.

	Cu2GeO4			Li2CuO2		CuGeO3
	$U_{ii}$ = 2 eV			$U_{ii}$ = 3.6 eV		$U_{ii}$ = 4.1 eV
	$J_1$	$J_2$	$J$	$J_1$	$J_2$	$J_1$	$J_2$
$t_{ij}$	119	86	121	70	81	205	76
$U_{ij}$	500	200	330	1510	680	1480	730
$4t_{ij}^2/{\tilde{U}}_{ij}$	37.8	16.4	35.1	9.4	9.0	64.2	6.9
$2J_{\mathrm{bare}}^F$	$-62.4$	$-15.0$	$-35.0$	$-88.6$	$-6.8$	$-72.4$	$-9.0$
$2J_{\mathrm{scr}}^F$	$-30.2$	$-2.8$	$-5.6$	$-43.0$	$-2.4$	$-34.4$	$-2.6$
$J_{ij}^{\mathrm{VASP}}$	$-0.2$	5.6	8.5	$-11.6$	5.5	18.7	3.7

Equations

$${\hat{\mathcal{H}}}^{\mathrm{Heis}}=\sum _{i>j}{J}_{ij}{\hat{\mathbf{S}}}_i{\hat{\mathbf{S}}}_j.$$

$${\hat{\mathcal{H}}}^{\mathrm{spin}}={\hat{\mathcal{H}}}^{\mathrm{Heis}}+\sum _{i>j}{\mathbf{D}}_{ij}[{\hat{\mathbf{S}}}_i\times {\hat{\mathbf{S}}}_j]+\sum _{i>j}{\hat{\mathbf{S}}}_i\text{tensor}{\Gamma }_{ij}{\hat{\mathbf{S}}}_j,$$

$${\mathcal{H}}_{LT}=\sum _{\mathbf{k}}{\mathbf{S}}_{\mathbf{k}}^{\ast }\text{tensor}J(\mathbf{k}){\mathbf{S}}_{\mathbf{k}},$$

$${\mathbf{S}}_i=\frac{1}{\sqrt{N}}\sum _{\mathbf{k}}{\mathbf{S}}_{\mathbf{k}}{e}^{i{\mathbf{kR}}_i}$$

$${\mathcal{H}}_{LT}=\sum _{\mathbf{k}\mu }{\omega }_{\mathbf{k}\mu }{S}_{\mathbf{k}\mu }^{\ast }{S}_{\mathbf{k}\mu }^{\ast },$$

$$\Gamma^{\mu\nu}_{J_{1}}=\left(\begin{array}[]{rrr}-0.03&0.00&0.00\\ 0.00&0.00&0.00\\ 0.00&0.00&0.03\end{array}\right).$$

$$\Gamma^{\mu\nu}_{J_{2}}=\left(\begin{array}[]{rrr}0.00&0.00&0.00\\ 0.00&0.00&0.02\\ 0.00&0.02&0.00\end{array}\right).$$

$$\Gamma^{\mu\nu}_{J}=\left(\begin{array}[]{rrr}-0.01&0.00&0.00\\ 0.00&0.00&0.03\\ 0.00&0.03&0.01\end{array}\right).$$

$$\Gamma^{\mu\nu}_{\sum_{J}}=\left(\begin{array}[]{rrr}-0.08&0.00&0.00\\ 0.00&0.00&0.10\\ 0.00&0.10&0.08\end{array}\right).$$

$$J_{ij}=\frac{4t_{ij}^2}{{\tilde{U}}_{ij}}+2J_{ij}^F,$$

$$F(\mathbf{q})=\int |W(\mathbf{r}){|}^2{e}^{-i\mathbf{qr}}d\mathbf{r}.$$

$$V_{ij}=\int \frac{W_i^{\ast }(\mathbf{r})W_i(\mathbf{r})W_j^{\ast }({\mathbf{r}}^{\prime })W_j({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d\mathbf{r}d{\mathbf{r}}^{\prime },$$

$$J_{ij}^F=\int \frac{W_i^{\ast }(\mathbf{r})W_j(\mathbf{r})W_j^{\ast }({\mathbf{r}}^{\prime })W_i({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d\mathbf{r}d{\mathbf{r}}^{\prime }.$$

$${\mathbf{P}}_{ij}\sim {\mathbf{ϵ}}_{ij}\times [{\mathbf{S}}_i\times {\mathbf{S}}_j],$$

$${\hat{\sigma }}^a{\hat{\sigma }}^b={\delta }_{ab}\hat{I}+i{ϵ}_{abc}{\hat{\sigma }}_c,$$

$$({\hat{\mathbf{S}}}_i{\hat{\mathbf{S}}}_j{)}^2=\frac{1}{16}({\hat{\mathbf{\sigma }}}_i{\hat{\mathbf{\sigma }}}_j{)}^2=\frac{1}{16}({\hat{\sigma }}_i^x{\hat{\sigma }}_j^x+{\hat{\sigma }}_i^y{\hat{\sigma }}_j^y+{\hat{\sigma }}_i^z{\hat{\sigma }}_j^z{)}^2=$$

$$=\frac{1}{16}({\hat{\sigma }}_i^x{\hat{\sigma }}_i^x{\hat{\sigma }}_j^x{\hat{\sigma }}_j^x+{\hat{\sigma }}_i^y{\hat{\sigma }}_i^y{\hat{\sigma }}_j^y{\hat{\sigma }}_j^y+{\hat{\sigma }}_i^z{\hat{\sigma }}_i^z{\hat{\sigma }}_j^z{\hat{\sigma }}_j^z+$$

$$+2{\hat{\sigma }}_i^x{\hat{\sigma }}_i^y{\hat{\sigma }}_j^x{\hat{\sigma }}_j^y+2{\hat{\sigma }}_i^x{\hat{\sigma }}_i^z{\hat{\sigma }}_j^x{\hat{\sigma }}_j^z+2{\hat{\sigma }}_i^y{\hat{\sigma }}_i^z{\hat{\sigma }}_j^y{\hat{\sigma }}_j^z)=$$

$$=\frac{1}{16}(3\hat{I}-2({\hat{\sigma }}_i^z{\hat{\sigma }}_j^z+{\hat{\sigma }}_i^y{\hat{\sigma }}_j^y+{\hat{\sigma }}_i^x{\hat{\sigma }}_j^x))=\frac{3}{16}\hat{I}-\frac{1}{2}{\hat{\mathbf{S}}}_i{\hat{\mathbf{S}}}_j.$$

Full text

††preprint: APS/123-QED

I Introduction

Copper oxides built by chains of edge-sharing CuO4 plaquettes serve as material prototypes of frustrated spin-$\frac{1}{2}$ chains with competing nearest-neighbor and next-nearest-neighbor interactions $J_{1}$ and $J_{2}$, respectively. This simple spin model received ample attention qua (2004) triggered by the prospects of chiral, multipolar, and spin-nematic phases that may occur therein Kolezhuk and Vekua (2005); Hikihara et al. (2008); Sudan et al. (2009); Heidrich-Meisner et al. (2009); Zhitomirsky and Tsunetsugu (2010); Sato et al. (2013); Balents and Starykh (2016). Whereas long-range order does not take place in one dimension (1D), interchain couplings in real materials will usually cause three-dimensional (3D) collinear or non-collinear order depending on the $J_{2}/J_{1}$ ratio. On the classical level, incommensurate spiral order appears for $J_{2}/|J_{1}|>\frac{1}{4}$, whereas at $J_{2}/|J_{1}|<\frac{1}{4}$ the second-neighbor coupling is not strong enough to tilt the spins, and the collinear ferromagnetic or up-down-up-down antiferromagnetic order form depending on the sign of $J_{1}$. Quantum effects preserve the spiral state in the case of ferromagnetic (FM) $J_{1}$ Zinke et al. (2009), but destroy the order and open a spin gap for antiferromagnetic (AFM) $J_{1}$ at $J_{2}/J_{1}>0.241$ White and Affleck (1996); Eggert (1996); Furukawa et al. (2010).

Real-world prototypes of the $J_{1}-J_{2}$ spin chains will typically follow one of these scenarios. The majority of quasi-1D copper oxides develop incommensurate spiral order Banks et al. (2009); Capogna et al. (2010); Zhao et al. (2012); Willenberg et al. (2012). Li2CuO2 Boehm et al. (1998); Lorenz et al. (2009), Ca2Y2Cu5O10, Fong et al. (1999); Kuzian et al. (2012), and CuAs2O4 Caslin et al. (2014) are notable exceptions, where $J_{1}$ is also FM, but spin alignment along the chains is purely ferromagnetic, owing to a smaller $J_{2}$. Spin-chain compounds with AFM $J_{1}$ are more rare, although tentative indications of the spin-gap formation at $J_{2}/J_{1}>0.241$ have been reported Lebernegg et al. (2017).

One puzzling case in this series is Cu2GeO4 Yamada et al. (2000) that reveals an unanticipated antiferromagnetic up-up-down-down (UUDD) order Zou et al. (2016) despite the prediction of FM $J_{1}$ and AFM $J_{2}$, both of the same magnitude Tsirlin et al. (2011a). This parameter regime would normally lead to the incommensurate spiral order, similar to LiCuVO4, CuCl2, and other $J_{1}-J_{2}$ cuprates. Here, we address this discrepancy and first analyze whether additional terms beyond $J_{1}$ and $J_{2}$ could destabilize the incommensurate order and give way to the UUDD state. This appears not to be the case, but instead $J_{1}$ is unusually weak in Cu2GeO4 and underlies the UUDD ground state of this compound.

The remainder of this paper is organized as follows. In Sec. II, we review the crystal structure of Cu2GeO4 and experimental information available for this material. Sec. III covers methodological aspects. In Sec. IV.1, we estimate both isotropic and anisotropic exchange interactions in Cu2GeO4, and in Sec. IV.2 analyze the ensuing magnetic ground state. Ferromagnetic direct exchange appears to be crucial and merits further analysis presented in Sec. IV.3 followed by the analysis of experimental magnetic susceptibility in Sec. IV.4 and a brief discussion and summary in Sec. V.

II Structure and properties of Cu2GeO4

Cu2GeO4 adopts a distorted spinel structure, where the Jahn-Teller effect inherent to Cu2+ transforms CuO6 octahedra into CuO4 plaquettes Hegenbart et al. (1981). The backbone of the structure is then formed by infinite chains of edge-shared plaquettes linked into a 3D network via the non-magnetic GeO4 tetrahedra (Fig. 1a).

Magnetic susceptibility measurements revealed a broad maximum around 80 K followed by an antiferromagnetic transition at $T_{N}\simeq 33$ K Yamada et al. (2000). This behavior is typical of low-dimensional and frustrated magnetism. In the case of Cu2GeO4, strong magnetic interactions are expected in the $ab$ plane, both along the chains of the Cu atoms ($J_{1},J_{2}$) and perpendicular to the chains ($J$), see Fig. 1b. The interactions $J_{c}$ between the planes are at least one order of magnitude weaker and form triangular loops together with $J_{1}$. This tentative magnetic model was confirmed by density-functional (DFT) band-structure calculations that yield $J_{1}\simeq-5.2$ meV (FM) and $J_{2}\simeq 6.9$ meV (AFM) as well as $J\simeq 11.2$ meV. Even if the leading coupling $J$ runs perpendicular to the copper chains, magnetic order along these chains is still determined by the competition between $J_{1}$ and $J_{2}$, similar to the 1D $J_{1}-J_{2}$ model. Detailed numerical analysis confirmed the stability of the spiral order along the copper chains as well as the collinear spin arrangement perpendicular to the chains, where no significant frustration occurs Tsirlin et al. (2011a).

Surprisingly, neutron diffraction data Zou et al. (2016) did not support this scenario and pinpointed the collinear UUDD order along the $J_{1}-J_{2}$ chains (Fig. 1b). This spin configuration is uncommon for cuprates and has never been seen in the $J_{1}-J_{2}$ compounds before. Biquadratic exchange was considered as the driving force of this unusual order Zou et al. (2016) and may explain it indeed Kaplan (2009), but appears irrelevant to Cu2GeO4, because biquadratic terms do not exist for spin-$\frac{1}{2}$ (they can be re-written as standard bilinear terms in the Hamiltonian Nagaev (1982); Mila and Zhang (2000), see Appendix A). Additionally, dielectric measurements revealed a clear anomaly in the permittivity at $T_{N}$, as well as a non-zero electric polarization that appears below $T_{N}$ in this formally centrosymmetric ($I4_{1}/amd$) crystal structure Zhao et al. (2018); Yanda et al. (2018). In the absence of spiral magnetic order that is typically associated with the electric polarization in chain cuprates Park et al. (2007); Naito et al. (2007); Seki et al. (2008); Mourigal et al. (2011), the origin of ferroelectricity in Cu2GeO4 remains controversial Zhao et al. (2018).

Here, we seek to throw some light on this problem from the ab initio perspective. The conclusion of Ref. Tsirlin et al., 2011a on the spiral order was based on the parametrization of an isotropic spin Hamiltonian, so it is natural to suspect, following Ref. Zou et al., 2016, that non-Heisenberg terms act against the spiral order and stabilize the UUDD one. We calculate such terms but find them to be small and affecting spin directions in the ordered state but not the nature of the ordered state itself. On the other hand, isotropic exchange couplings of Ref. Tsirlin et al., 2011a have to be revised, eventually giving a clue to the formation of the UUDD order in Cu2GeO4.

III Methods

In Ref. Tsirlin et al., 2011a, the magnetic behavior of Cu2GeO4 was analyzed on the level of the Heisenberg spin Hamiltonian,

[TABLE]

With $J_{1}\simeq-5.2$ meV and $J_{2}\simeq 6.9$ meV Tsirlin et al. (2011a), it leads to the spiral order along the copper chains at odds with the experiment. To account for the experimental UUDD order, additional terms may be invoked as follows,

[TABLE]

where $\mathbf{D}_{ij}$ are Dzyaloshinskii-Moriya (DM) vectors and $\Gamma_{ij}$ are symmetric anisotropy tensors. The latter favor collinear spins and, therefore, may stabilize the UUDD order over the spiral one, whereas the former do not stabilize collinear spin configurations per se, but may act against the spiral state. Specifically, in Cu2GeO4 the alternating directions of $\mathbf{D}_{1}$ (Fig. 1b) are incompatible with the continuous spin rotation in the spiral. Biquadratic and other higher-order corrections do not appear as independent terms in the spin-$\frac{1}{2}$ Hamiltonian Nagaev (1982); Mila and Zhang (2000), see also Appendix A.

Magnetic exchange parameters are obtained from DFT calculations performed within the generalized gradient approximation (GGA) Perdew et al. (1996) implemented in Vienna ab initio Simulation Package (VASP) Kresse and Furthmüller (1996a, b). Additionally, the full-potential FPLO Koepernik and Eschrig (1999) and ELK Dewhurst et al. codes were used. The crystal structure given in Ref. Hegenbart et al., 1981 was employed in all calculations, similar to Ref. Tsirlin et al., 2011a.

In the absence of electronic correlations, Cu2GeO4 features a metallic band structure with several bands crossing the Fermi level. The complex of four bands between $-0.6$ and 0.6 eV corresponds to four Cu atoms in the primitive cell and arises from $d_{x^{2}-y^{2}}$ orbitals that are half-filled in Cu2+. Electronic correlations split these bands and open a gap, in agreement with the insulating behavior of Cu2GeO4 expected from the green sample color Hegenbart et al. (1981); Zhao et al. (2018). The effect of correlations is modeled on the DFT+$U$+SO level, with all parameters of the spin Hamiltonian, Eq. (2), extracted from total energies of ordered spin configurations using the mapping procedure Xiang et al. (2011). Alternatively, we perform a model analysis based on hopping parameters of the uncorrelated band structure and additionally calculate ferromagnetic contribution to the exchange from the overlap of Wannier functions, as further explained in Sec. IV.3. Similar methodology has been used in the previous DFT study Tsirlin et al. (2011a), but several technical details were different and proved to be crucial, as we show below.

The DFT+$U$+SO method relies on the parametrization of the Coulomb and Hund’s exchange interactions $U_{d}$ and $J_{H}$, respectively.In Ref. Tsirlin et al., 2011a, $U_{d}=6.5$ eV and $J_{H}=1$ eV were chosen empirically along with the around-mean-field (AMF) double-counting correction scheme that is more suitable for correlated metals. Here, we use the double-counting correction in the fully localized limit appropriate for insulators and obtain $U_{d}-J_{H}\sim 8.5$ eV from the linear-response method Cococcioni and de Gironcoli (2005). Assuming $J_{H}=1$ eV, we find $U_{d}=9.5$ eV, which is similar to the parametrization that is typically used for copper oxides Janson et al. (2012); Lebernegg et al. (2013); Prishchenko et al. (2017) in conjunction with the FLL flavor of the double-counting correction.

Magnetic ground state of the spin Hamiltonian is obtained from the Luttinger-Tisza (LT) method considering spins as classical moments Luttinger and Tisza (1946); Lyons and Kaplan (1960),

[TABLE]

where $\mathbf{S}_{\mathbf{k}}$ is the Fourier transform of the spin:

[TABLE]

Diagonalization of Eq. (3) yields Liu et al. (2016); Sklan and Henley (2013)

[TABLE]

where $S_{\mathbf{k}\mu}=\mathbf{S}_{\mathbf{k}}\hat{\mathbf{e}}_{\mathbf{k}\mu}$, $\omega_{\mathbf{k}\mu}$, and $\hat{\mathbf{e}}_{\mathbf{k}\mu}$ are corresponding eigenvalues and eigenvectors of $\tensor{J}(\mathbf{k})$. The LT mode $S_{\mathbf{k}\mu}$ with the most negative eigenvalue $\omega_{\mathbf{k}\mu}$ is considered as an ”optimal” mode with the wave vector $\mathbf{Q}_{\rm LT}$. If the constructed spin state $\{\mathbf{S}_{i}\}$ is the linear combination of the optimal LT modes and complies with the ”strong constraint” of $|\mathbf{S}_{i}|^{2}=1$, it can be considered as a ground state Sklan and Henley (2013).

We also calculate magnetic susceptibility of Cu2GeO4 using the loop algorithm Todo and Kato (2001) of the ALPS simulation package Albuquerque et al. (2007). To this end, finite lattices with up to $16\times 16$ sites and periodic boundary conditions were used.

IV Results

IV.1 Microscopic magnetic model

Isotropic exchange couplings of the Heisenberg spin Hamiltonian, Eq. (1), are listed in Table 1. DFT calculations were performed in three different codes that delivered largely consistent results for $J_{2}$ but not for $J_{1}$ that varies between $-0.2$ meV in VASP and $-7.2$ meV in FPLO, with ELK returning an intermediate value. This large spread of $J_{1}$ leads to a highly ambiguous physical picture, because the competition between $J_{1}$ and $J_{2}$ may be either strong (FPLO) or nearly non-existent (VASP).

Other exchange couplings, including $J$ and $J_{2}$, are largely consistent between the different band-structure codes. This alone indicates that the ambiguity of $J_{1}$ does not stem from numerical inaccuracies, but reflects a complex nature of the coupling, which is short-range and combines dissimilar contributions of the direct exchange and superexchange, as opposed to the long-range couplings $J$ and $J_{2}$ dominated by the superexchange. To further exclude any technical issues related to the basis sets or energy convergence, we performed spin-polarized DFT calculations with $U_{d}=J_{H}=0$ and arrived at $J_{1}\sim$ 40 meV in all three codes 111The remaining spread of 3 meV, which is less than 10% of the absolute value, may be caused by the variation of magnetic moments between different spin configurations. This effect is unavoidable in uncorrelated calculations that underestimate local magnetic moment on Cu2+. (Table 2), thus confirming that the ambiguity of $J_{1}$ arises not from the different basis sets and not from the different treatment of the crystal potential, but from the way the DFT+$U$ correction is applied in each code. We also note that a variation of $U_{d}$ within the reasonable range of $1-2$ eV does not improve the consistency between VASP and FPLO.

We also considered whether the mapping procedure could accidentally fail for $J_{1}$ if, for example, different spin configurations produced largely different local moments. However, this was not the case (Table 2), and we are led to conclude that the problem with calculating $J_{1}$ is intrinsic. This coupling combines two contributions of different nature, the ferromagnetic direct exchange and antiferromagnetic superexchange, and different band-structure codes deliver very different estimates of these competing contributions.

We return to this $J_{1}$ problem in Sec. IV.3 but first consider anisotropic, non-Heisenberg terms that may also affect the ground state. These terms are obtained in VASP, because it delivers the most realistic estimate of $J_{1}$, as we show below.

Anisotropic exchange is driven by the spin-orbit (SO) coupling. The effect of SO can be seen from the weak band splitting near the Fermi level at some of the high-symmetry $k$-points (Fig. 3). The orbital moment of Cu2+ reaches its highest value of 0.18 $\mu_{B}$ for the direction perpendicular to the CuO4 plaquettes, similar to other cuprates Badrtdinov et al. (2016); Volkova et al. (2012).

DM components for $J$ and $J_{2}$ are forbidden by the inversion symmetry. Therefore, the only non-vanishing DM vector is $\mathbf{D}_{1}$ that should lie in the $ab$ plane and perpendicular to the copper chains by virtue of the two mirror planes, one of them containing both Cu atoms and the other one passing through the middle of the Cu–Cu bond. From DFT+$U$+SO we find ${\bf D}_{01}$ = (0.01, 0, 0) meV for the plane with the copper chains running along the $b$ direction. In the neighboring planes with the Cu chains along $a$, the $\mathbf{D}_{1}$ vector has the same length but points along $b$ instead of $a$.

Symmetric components of the anisotropy for the interacting copper pairs shown in Fig. 1b are similar in magnitude to the above DM vector (in meV),

[TABLE]

[TABLE]

[TABLE]

Here, the two-fold rotation axis along $c$ and the $bc$ mirror plane cancel all off-diagonal components of the $\tensor{\Gamma}_{J_{1}}$ tensor. In the case of $J_{2}$ and $J$, the rotation axis is missing, such that the nonzero $bc$ and $cb$ components become allowed. Taken together, these three tensors define the overall anisotropy matrix

[TABLE]

Its lowest eigenvalue defines $a$ as the easy direction for the layer with the copper chains running along $b$. Similar to $\mathbf{D}_{1}$, this easy direction changes to $b$ in the adjacent layers with the copper chains running along $a$.

Although small compared to the isotropic exchange couplings, these anisotropic terms play crucial role in choosing spin directions in the magnetically ordered state (Sec. IV.2). It is also worth noting that calculated anisotropic terms fulfill all symmetries of the system, and this fact lends credence to the DFT+$U$+SO results. Anisotropic interactions in cuprates are of superexchange nature Shekhtman et al. (1992, 1993) and thus easier to evaluate than $J_{1}$, which is strongly influenced by the direct exchange (Sec. IV.3).

IV.2 Model solution

We shall now use the LT method to determine the magnetic ground state. It is instructive to apply this method to the $J_{1}-J_{2}$ Heisenberg model first. Fig. 4 shows the LT wave vector depending on the $J_{1}/J_{2}$ ratio. The spiral order spans the region $-4<J_{1}/J_{2}<4$, as expected. However, at $J_{1}/J_{2}=0$ the wave vector $\mathbf{Q}_{\rm LT}=\pi/2$ corresponds to two possible solutions, the spin spiral with the pitch angle of $\pi/2$ (orthogonal spin configuration) and the collinear UUDD order that has been observed in Cu2GeO4 experimentally Zou et al. (2016). We may thus expect the UUDD order at $J_{1}=0$.

Fig. 5 shows the LT wavevectors obtained for the full set of in-plane exchange couplings ($J_{1}$, $J_{2}$, $J$, and $J_{ab}$) from FPLO and VASP. The FPLO results clearly lead to the incommensurate position of the minimum and stabilize the spiral order. Weak anisotropic terms presented above do not change the $\mathbf{q}$-vector significantly. On the other hand, the VASP results produce the minimum at $\mathbf{q}=(\pi/2,\pi/2)$ compatible with the UUDD order or with the spin spiral having the $\pi/2$ rotation along the copper chain. The former state is collinear and thus benefits from the symmetric anisotropy, the $xx$-term of $\Gamma^{\mu\nu}_{\sum_{J}}$. The spiral state will, on the other hand, gain less energy from $\Gamma^{\mu\nu}_{\sum_{J}}$, because different spin directions are present. This has been confirmed by a direct DFT+$U$+SO calculation for the UUDD state and for the spin spiral with the pitch angle of $\pi/2$. The UUDD state is lower in energy by 0.45 meV per Cu2+ ion.

We conclude that the UUDD order can be competitive with the spiral order, but around $J_{1}=0$ only. As soon as this condition is fulfilled, symmetric anisotropy present in Cu2GeO4 favors the UUDD order over the spiral one, because all spins can follow easy direction in the UUDD state but not in the spiral. This way, the negligibly small $J_{1}$ is the necessary condition for the formation of the UUDD order. We shall further justify this condition in Sec. IV.3 below, but first demonstrate that the combination of $J_{1}=0$ and weak symmetric anisotropy explains not only the UUDD order along the copper chains, but also all other features of the experimental magnetic structure.

In the absence of the interlayer coupling, the symmetric anisotropy $\Gamma^{\mu\nu}_{\sum_{J}}$ puts spins along $b$ in the layer where the copper chains run along $a$, and along $a$ in the layer where the copper chains run along $b$. This would lead to orthogonal spin directions in the neighboring layers (Fig. 6a) and becomes compatible with the scenario of frustrated interlayer couplings $J_{c}$. However, the UUDD order releases the frustration on those tetrahedra, where spins are parallel along the $J_{1}$ bonds, and such tetrahedra may gain energy from $J_{c}$. Then $\mathbf{a}+\mathbf{b}$ or $\mathbf{a}-\mathbf{b}$ are chosen as compromise spin directions between the two layers (Fig. 6b). In contrast, the tetrahedra with antiparallel spins along $J_{1}$ remain frustrated and enjoy the orthogonal spin arrangement (Fig. 6b). This leads to the peculiar magnetic order observed in Cu2GeO4. The spin direction alternates between $\mathbf{a}+\mathbf{b}$ and $\mathbf{a}-\mathbf{b}$ in every second layer in response to the frustration present on one half of the Cu4 tetrahedra and absent on the other half.

The DM interactions were not considered so far, because they neither stabilize nor destabilize the collinear UUDD state. They may, however, introduce weak spin canting as shown in Fig. 7. This canting is fully compensated within each chain and does not produce any net magnetic moment. From the $\mathbf{D}_{1}$ value obtained in Sec. IV.1 and from the direct relaxation of the magnetic structure within VASP, we estimate only a weak non-collinearity with the canted moment of about 0.005 $\mu_{B}$. Such a moment is clearly too small to be detected by powder neutron diffraction Zou et al. (2016), but is allowed by symmetry and may be relevant to the development of electric polarization, as we further explain in Sec. V.

IV.3 Direct exchange

Having established $J_{1}=0$ as the necessary condition for the UUDD order, we now discuss the microscopic origin of $J_{1}$ and the reasons for the full compensation of this coupling in Cu2GeO4. Magnetic couplings in insulators are generally composed of two contributions, the kinetic term due to the superexchange, $J^{\rm kin}_{ij}=4t^{2}_{ij}/\widetilde{U}_{ij}$, and the potential term due to the direct exchange interaction $J^{F}_{ij}$ arising from the direct overlap of the magnetic orbitals Anderson (1959). We then write an isotropic exchange coupling in the form

[TABLE]

where $t_{ij}$ is the hopping integral, $\widetilde{U}_{ij}=U_{ii}-U_{ij}$ is an effective screened Coulomb repulsion parameter Mazurenko et al. (2007); Badrtdinov et al. (2018), and $J^{F}_{ij}$ is the direct exchange.

The direct exchange depends on the overlap between the magnetic orbitals. This overlap is very sensitive to hybridization effects, because spin polarization spreads onto ligands, which contribute to the overlap and largely determine the $J^{F}_{ij}$ values in real materials. This hybridization effect can be captured using Wannier functions that serve as a realistic representation of the magnetic orbitals. Here, we use maximally localized Wannier functions for Cu2+ Marzari and Vanderbilt (1997) and illustrate the role of the hybridization by calculating three-dimensional magnetic form-factors $F(\mathbf{q})$ as Fourier transforms of the Wannier orbitals Mazurenko et al. (2015),

[TABLE]

In Fig. 8, we compare two scenarios: i) Wannier functions calculated for four Cu $d_{x^{2}-y^{2}}$ bands in the vicinity of the Fermi level; ii) Wannier functions calculated for all Cu $3d$ and O $2p$ states. The second case leads to the lower oxygen contribution and renders $F(\mathbf{q})$ more symmetric, similar to the purely ionic form-factor for Cu2+. In contrast, case i) captures the full Cu–O hybridization, makes $F(\mathbf{q})$ less symmetric, and causes a faster decay at higher $q$’s. Similar effects were reported for other Cu2+ oxides Mazurenko et al. (2015); Walters et al. (2009) and may be responsible for the reduced ordered moment of $0.89(5)$ $\mu_{B}$ determined by neutron diffraction using the ionic form-factor for Cu2+ Zou et al. (2016).

We are now in a position to calculate bare Coulomb (both on-site $V_{ii}$ and intersite $V_{ij}$) and non-local direct exchange $J^{F}_{ij}$ integrals in the basis of Wannier functions $W_{i}(\bm{r})$,

[TABLE]

and

[TABLE]

Screening effects were captured on the level of random phase approximation (RPA) Ehrenreich and Cohen (1959). However, energy bands at the Fermi level are also involved in the screening processes and cause a ”self-screening” that needs to be excluded in order to evaluate the partially screened, realistic Coulomb parameters Aryasetiawan et al. (2006). Therefore, we utilized constrained RPA and obtained the on-site $U_{ii}$ and intersite $U_{ij}$ Coulomb parameters listed in Table 3. As for the non-local direct exchange, its evaluation within constrained RPA requires a very accurate integration within the Brillouin zone and proved to be unfeasible. Therefore, we used standard RPA and calculated the fully-screened $J^{F}_{\rm scr}$ that gives the lower bound for the FM contribution and also allows for a comparison between different exchange pathways and different compounds.

Wannier representation of the band structure also gives access to the hopping integrals $t_{ij}$. This way, we obtain both FM and AFM contributions to the exchange couplings in Cu2GeO4, as listed in Table 3. In the case of $J_{2}$ and $J$, AFM superexchange clearly dominates over the fully screened direct exchange, and the overall AFM couplings ensue. On the other hand, $J^{F}_{\rm scr}$ for $J_{1}$ is only slightly smaller in magnitude than $J^{\rm kin}$, suggesting that $J_{1}$ may be close to zero, as VASP DFT+$U$ calculations predict.

It is also instructive to juxtapose Cu2GeO4 with other compounds containing copper chains. To this end, we choose Li2CuO2 with its large FM $J_{1}\simeq-19.6$ meV Lorenz et al. (2009) and CuGeO3 where $J_{1}$ was proposed to be AFM Castilla et al. (1995). The main structural difference between these compounds lies in the nearest-neighbor Cu–O–Cu angle that increases from $93.97^{\circ}$ in Li2CuO2 (at 1.5 K) Sapiña et al. (1990) to $99.24^{\circ}$ in CuGeO3 (at 20 K) Braden et al. (1996). The nearest-neighbor hopping is, consequently, enhanced and makes $4t_{1}^{2}/\tilde{U}_{ij}$ much larger than $2J^{F}_{\rm scr}$ (Table 3). This way, the crossover from FM $J_{1}$ in Li2CuO2 to AFM $J_{1}$ in CuGeO3 is caused by the increased Cu–O–Cu angle, while all other microscopic parameters of these compounds are similar.

Coming now to Cu2GeO4, we realize that its hopping integral $t_{1}$ and, thus, the AFM contribution to $J_{1}$ are intermediate between those of Li2CuO2 and CuGeO3. The FM contribution $2J^{F}_{\rm scr}$ is, on the other hand, reduced in magnitude. Both aspects create suitable conditions for the cancellation of the FM and AFM contributions, leading eventually to $J_{1}\simeq 0$. Such an unusual behavior may be rooted in peculiarities of the Cu2GeO4 crystal structure. The nearest-neighbor Cu–O–Cu angle of $91.57^{\circ}$ Zhao et al. (2018) is in fact lower than in Li2CuO2, so naively one would expect an even lower $t_{1}$, which is not the case. Weak buckling of the copper chains (Fig. 1a) appears to be crucial here, because it reduces the direct overlap responsible for $J^{F}_{\rm scr}$ and, simultaneously, shortens the nearest-neighbor Cu–Cu distance from 2.860 Å in Li2CuO2 Sapiña et al. (1990) to 2.796 Å in Cu2GeO4 Zhao et al. (2018), thus enhancing the direct $d-d$ hopping. A similar argument can be applied to linarite, PbCu(SO4)(OH)2, that also features buckled copper chains with an intermediate Cu–Cu distance of 2.823 Å. Indeed, its $J_{1}\simeq-8.6$ K Rule et al. (2017) or $-13.8$ K Cemal et al. (2018) is smaller in magnitude than in Li2CuO2, but still on the ferromagnetic side.

IV.4 Magnetic susceptibility

Further support for the $J_{1}\simeq 0$ scenario can be garnered by analyzing magnetic susceptibility of Cu2GeO4. According to Table 1, the minimal magnetic model for this compound should only include the coupling $J$ that forms spin chains perpendicular to the structural chains of the Cu atoms. Adding the coupling $J_{2}$ connects these spin chains into a rectangular spin lattice. We used both models to fit the experimental susceptibility data from Ref. Yamada et al., 2000.

The spin-chain model leads to a decent fit above 100 K with $J=12.1$ meV, $g=2.23$, and the temperature-independent term $\chi_{0}=-7.7\times 10^{-5}$ emu/mol, but at lower temperatures this model overestimates the experimental susceptibility (Fig. 9), suggesting that antiferromagnetic interchain couplings are at play. By including $J_{2}$, we obtain an excellent fit down to $T_{N}$ with $J=10.7$ meV, $J_{2}=5.3$ meV, $g=2.33$, and $\chi_{0}=-0.0001$ emu/mol (Fig. 9). The fitted values of $J$ and $J_{2}$ are in good agreement with the DFT estimates in Table 1. Moreover, we confirm that the experimental susceptibility of Cu2GeO4 is compatible with the $J_{1}\simeq 0$ scenario.

V Discussion and summary

We have shown that the UUDD magnetic structure of Cu2GeO4 can be obtained in the limit of the weak nearest-neighbor coupling $J_{1}\simeq 0$. Two additional ingredients, frustration on half of the Cu4 tetrahedra and orthogonal easy directions in the adjacent copper layers, explain all peculiarities of the experimental magnetic structure with its two spin directions, $\mathbf{a}\pm\mathbf{b}$, that change in every second layer Zou et al. (2016). We revise the previous ab initio results Tsirlin et al. (2011a) and establish the new magnetic model of Cu2GeO4 compatible with both experimental magnetic susceptibility and ground state, thus resolving the discrepancies regarding the magnetic behavior of this compound. Several remarks are in order, though.

First, $J_{1}$ must be small, but no symmetry argument requires a complete cancellation of this coupling. The possible range of the $J_{1}$ values is determined by the energy difference between the UUDD and spiral states, as compared to the energy gain from the symmetric anisotropy. Quantum effects neglected within our LT analysis may also play a role here Zinke et al. (2009). Detailed estimates go beyond the scope of our present manuscript but may be interesting if the symmetric anisotropy would be determined experimentally, e.g., by measuring magnon gap with electron spin resonance or THz spectroscopy.

Second, DFT proves incapable of estimating $J_{1}$ in a consistent manner (Table 1). Similar problems were encountered for other short-range couplings in copper and vanadium compounds Tsirlin et al. (2010, 2011b), although the Cu2GeO4 case appears to be most severe, because not only different flavors of DFT+$U$ but also different band-structure codes return largely different values of $J_{1}$. We attempted to vary the Coulomb repulsion $U_{d}$ and to change the double-counting correction, but were unable to reduce $|J_{1}|$ below 2 meV using FPLO as the full-potential code. This indicates that a lot of caution should be taken in analyzing the short-range couplings obtained from DFT+$U$. On the more positive side, Cu2GeO4 may be an excellent test case for ab initio methods, because the $J_{1}\simeq 0$ condition is very robust. Any significant deviation from it leads to the spiral order, which is not observed experimentally.

Third, spin canting caused by the DM interactions may give a clue to the formation of local electric polarization. The inverse DM mechanism triggers the polarization Katsura et al. (2005)

[TABLE]

where $\bm{\epsilon}_{ij}=(0,1,0)$ is a vector connecting the magnetic sites $i$ and $j$ along the copper chains, and in a given layer the spins are presented by $\mathbf{S}_{i}=(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},\delta)$, $\mathbf{S}_{j}=(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},-\delta)$ for the up-up pair or $\mathbf{S}_{i}=(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},-\delta)$, $\mathbf{S}_{j}=(-\frac{1}{\sqrt{2}},-\frac{1}{\sqrt{2}},-\delta)$ for the up-down pair, where $\delta$ is small canting. Both pairs produce electric polarization of the same sign directed along $c$. This way, each copper layer generates a finite electric polarization that, however, cancels out between the neighboring layers following the symmetry of the $I_{c}\bar{4}2d$ magnetic space group Zou et al. (2016). Nevertheless, it is conceivable that weak structural changes in the magnetically ordered state may reduce the symmetry, thus leading to a non-zero polarization. The step-like changes in the magnetic susceptibility and permittivity at $T_{N}$ Zhao et al. (2018), as well as the abrupt onset of the polarization in the magnetically ordered state Zhao et al. (2018); Yanda et al. (2018), may indicate a weak first-order nature of the magnetic transition, similar to $\alpha$-CaCr2O4, where the electric polarization has also been observed Singh et al. (2011) at odds with the symmetry of the magnetic structure Chapon et al. (2011). On the experimental side, further thermodynamic measurements probing the nature of the magnetic transition in Cu2GeO4, as well as dielectric measurements probing the direction of the electric polarization, can be useful.

In summary, we have shown that the collinear UUDD magnetic order in Cu2GeO4 is only possible in the $J_{1}\simeq 0$ limit and should be traced back to the nearly perfect compensation of the FM and AFM contributions to this exchange coupling. The UUDD order along the copper chains removes the frustration on half of the Cu4 tetrahedra and, together with the weak symmetric anisotropy, leads to the peculiar magnetic structure with two different spin directions, as observed experimentally.

Acknowledgements.

We would like to thank Yaroslav O. Kvashnin (Uppsala University) for his help with the ELK code. The work of D.I.B. was funded by RFBR according to the research project No 18-32-00018. The work or V.V.M. was supported by the Russian Science Foundation, Grant No 18-12-00185. A.A.T. acknowledges financial support by the Federal Ministry for Education and Research through the Sofja Kovalevskaya Award of Alexander von Humboldt Foundation.

Appendix A Biquadratic term in case of $S$ = 1/2

Here we show that in the case of spin-$\frac{1}{2}$ the biquadratic term $(\hat{\mathbf{S}}_{i}\hat{\mathbf{S}}_{j})^{2}$ can be re-written in the bilinear form $\hat{\mathbf{S}}_{i}\hat{\mathbf{S}}_{j}$. To this end, we use the property of the Pauli matrices,

[TABLE]

where $\delta_{ab}$, $\hat{I}$, and $\epsilon_{abc}$ are the Kronecker delta, identity matrix, and Levi-Civita symbol, respectively. Using the above relation and the commutation rule for two different sites, $[\hat{\sigma}_{i}^{a},\hat{\sigma}_{j}^{b}]=0$, it is straightforward to show that

[TABLE]