Comment on `Oxygen vacancy-induced magnetic moment in edge-sharing CuO$_{2}$ chains of Li$_{2}$CuO$_{2}$'

In a recent work devoted to the magnetism of Li$_{2}$CuO$_{2}$, Shu et al. [New J. Phys. 19 (2017) 023026] have proposed a"simplified"unfrustrated microscopic model that differs considerably from the models refined through decades of prior work. We show that the proposed model is at odds with known experimental data, including the reported magnetic susceptibility $\chi(T)$ data up to 550~K. Using an 8$^{\rm th}$ order high-temperature expansion for $\chi(T)$, we show that the experimental data for Li$_{2}$CuO$_{2}$ are consistent with the prior model derived from inelastic neutron scattering (INS) studies. We also establish the $T$-range of validity for a Curie-Weiss law for the real frustrated magnetic system. We argue that the knowledge of the long-range ordered magnetic structure for $T<T_N$ and of $\chi(T)$ in a restricted $T$-range provides insufficient information to extract all of the relevant couplings in frustrated magnets; the saturation field and INS data must also be used to determine several exchange couplings, including the weak but decisive frustrating antiferromagnetic (AFM) interchain couplings.

Li 2 CuO 2 takes a special place among the still increasing family of frustrated chain compounds with edge-sharing CuO 4 plaquettes and a ferromagnetic (FM) nearest neighbor (NN) in-chain coupling J 1 [1]. This unique position is due to its ideal planar CuO 2 chain structure and its well-defined ordering characterized by a 3D Neél-type arrangement of adjacent chains whose magnetic moments are aligned ferromagnetically along the chains (b-axis). Li 2 CuO 2 is well studied in both experiment and theory (see e.g. [2][3][4][5][6][7][8][9][10][11]) and serves nowadays as a reference system for more complex and structurally less ideal systems. In particular, it is accepted in the quantum magnetism community that the leading FM coupling is the NN inchain coupling J 1 . (J 1 is also dominant but antiferromagnetic (AFM) in the special spin-Peierls case of CuGeO 3 [12].) There is always also a finite frustrating AFM next-nearest neighbor (NNN) coupling J 2 > 0, see figure 1, left. This inchain frustration is quantified by α = J 2 / |J 1 |. In the present case, and in that of the related Ca 2 Y 2 Cu 5 O 10 , there are only frustrating AFM interchain couplings (IC) with adjacent chains shifted by half a lattice constant b. In this lattice structure there is no room for unfrustrated perpendicular IC. This AFM IC with NN and NNN components plays a decisive role in the stabilization of the FM alignment of the magnetic moments along the chain direction. Although weak at first glance, with eight NN and NNN it is nevertheless significant enough (by a factor of two) to prevent a competing non-collinear spiral type ordering in Li 2 CuO 2 (with the frustration ratio the skew AFM frustrating IC's J 5 and J 6 between adjacent chains responsible for the FM inchain ordering but ignored in [13], similarly as the generic AFM NNN inchain coupling J 2 . There, the FM NN intrachain coupling J 1 is underestimated by a factor of four (see table 1). Green line: the weak NNN IC J 3 ≡ J a claimed to be dominant and FM by Shu et al. [13]). Right: Main: The inverse spin susceptibility for a magnetic field along the a axis of Li 2 CuO 2 (the sample that was used for INS studies, cf. figure 4 of [6]) fitted by the eight-order high-T expansion expression (10). Inset: Convergence of the two pseudo-CW parameters to their high-T asymptotic values: α > 1/4), as often observed for other members of this family with unshifted chains [1]. All these well-established features were practically excluded by Shu, Tian, Lin et al. (STL) [13], proposing instead (i) a very nonstandard unfrustrated model (dubbed hereafter as STL-model) with comparable couplings in all directions, and where the leading FM coupling is given by an unphysically large NNN FM IC J a (denoted as J 3 therein) perpendicular to the chains in the basal ab-plane (J a =−103 K for stoichiometric and −90 K in the presence of O vacations in Li 2 CuO 2−δ with δ = 0.16). (ii) The coupling between the NN chains, as derived from the inelastic neutron scattering (INS) data [6] and in qualitative accord with the results of LDA+U calculations [5], has been ignored and replaced ad hoc by an artificially large, "effective" non-frustrated AFM IC J (see the right of figure 9 in [13] with a 4-fold coordination) absent in real Li 2 CuO 2 .
In the present Comment, we show that this parametrization is a direct consequence of an incorrect analysis of their susceptibility χ(T ) data in addition to ignoring the highly dispersive magnon mode and its local softening observed by INS. We admit that the very question about an influence of O vacancies on the magnetic properties of Li 2 CuO 2 raised by Shu et al. [13] is interesting and should be studied; however, this must be done within a proper analysis based on a realistic phenomenological model reflecting the established large |J 1 | values exceeding 200 K [6] and excludes a simple Curie-Weiss (CW) law below 1000 K. In this context we mention similar mistakes made in the literature, where even an artificial AFM Θ CW < 0 has been found [2,[14][15][16] prior to 2009, when the large value of J 1 was not yet known.
Li 2 CuO 2 is a frustrated quasi-1D system that has been well studied during the last decades. The first two rows of table 1. (see table 6 of [13]) provide the J's suggested by Shu et al. from their qualitative simulation of the magnetic ordering and an analysis of the measured χ(T ) for two samples with different O content. The main striking difference between these sets from all previous ones is the absence of both magnetic frustration and of the quasi-1D regime with a dominant J 1 realized in all edge-sharing CuO 2 chain compounds. Moreover, the proposed sets are evidently at odds with the results of two INS studies [6,16]. The set derived from the INS in the bc plane [6] see the last row in table 1 (and table 6 of [13]) does explain the χ(T ) data for T > T N , especially when supplemented by a weak AFM IC to the a axis (in accord with the reported weakly dispersing magnon in that direction [16,17]).
We consider the Heisenberg spin-Hamiltonian where m enumerates the sites in the magnetic (Cu) lattice (see figures 1 and 9 in [13]), J r is the interaction of a pair of spinsŜ m andŜ m+r . Note the different notation [16,18,19], where the same interaction is denoted −2J r . The form of equation (1) implies positive (negative) signs for AFM (FM) couplings. Shu et al. [13] use the same notation in tables 5 and 6, but use the wrong signs in equations (8) and (10). A small anisotropy of the couplings seems to be unimportant for the analysis of χ(T ) and is ignored here.
The magnetic field H z is directed along the easy axis (i.e. the crystallographic a-axis in Li 2 CuO 2 ). Finally, µ B is the Bohr magneton and g a is the gyromagnetic ratio for this direction. Shu et al. use two approaches when analyzing χ(T ) of Li 2 CuO 2 . First, they fit it in the range 250 K< T < 550 K using a CW-law They then extract three effective couplings instead of six original ones by fitting χ(T ) with an RPA-like expression derived for quasi-1D systems. Note that Shu et al.
give an obviously erroneous form in their equation (7) with the factor [1 − 2 (z J + z J )], being a sum of the dimensionless value 1 and a value with the dimension of energy. No figure is shown for χ(T ) fitted by their curve, so it is impossible to evaluate their fit, nor to estimate its quality and validity range. However, in [18] (reference 33 of their paper), the correct expression reads: where J=J 1 is the inchain coupling and J =J 3 , J ∼ 2J 5 are IC's with z =2, z =4 the corresponding numbers of neighbors (see figures 9(b) in [13] for a simplified structure, which differs from the real one, see figure. 9(a) therein). In equation (3) we have accounted for the different notations of the exchange couplings between [13] (the same as ours) and [18]. Note that the quasi-1D regime assumed in equation (3) implies J (z J + z J ), which is obviously violated by the STL model. We recall that a CW-law exactly reproduces the high-T behavior of the spin susceptibility of any system described by the Heisenberg Hamiltonian (1) with the CW temperature Θ CW Equation (4) is the exact result of a high-T expansion (HTE) of the susceptibility [20] (see equation (27) of section IV.B in [21], see also [22], and references therein), which is valid for any Heisenberg system. As mentioned above, the expression for Θ CW given in equation (10) of [13] has a wrong sign. It gives positive (negative) contribution for AFM (FM) interactions in conflict with the physical meaning of Θ CW . Now we show that χ q1D also obeys a CW-law for large enough T . We recast equation (3) in the form where C denotes the Curie constant, z = 2 is the in-chain coordination number, and The Θ CW = +14 K for the STL-model given by equation (4) should coincide with the value obtained by the CW fit presented in tables 3 and 4 of [13], provided the CW fit is justified for the chosen range of T . We stress that a value of Θ CW > 0 doesn't cause a divergence of χ(T ) at T = Θ CW , in contrast to the erroneous claim of Shu et al. in speculations after their equation (10). In general, a divergence of the susceptibility χ(Q 0 , T → T 0 ) means the emergence of a long-range order in a magnetic system characterized by a spontaneous magnetization m(Q 0 ) for T < T 0 . A ferrimagnetic (FIM) and a FM ordering correspond to Q 0 at the center of Brillouin zone (BZ) and to the uniform component of χ(0, T ) ≡ χ(T ), respectively. Note that a FIM ordering and the divergence of χ(T → T C ) may occur for systems with purely AFM couplings and negative Θ CW due to the geometry of the spin arrangement (see e.g. figure 4 of [23]). A Q 0 = (0, 0, 0) corresponds to a helimagnetic or to an AFM ordering. Then, the uniform χ(T ) remains finite at T 0 . An AFM ordering corresponds to Q 0 located at the edge of the BZ. This is the case for the prior model of Li 2 CuO 2 , as well as for the STL-model. The range of validity of the CW-law for χ q1D [equation (3)] is As already noted, the approximate expression equation (3) is relevant for a quasi-1D system, where the inchain J's dominates. This is true also for the condition (7). Let us establish now a general condition for the applicability of a CW-law. For this aim it is convenient to use the inverted exact HTE [24] for spin-1/2 systems with equivalent sites (cf. equations 5(a) and (b) in [24]). Thus, a CW-law with Θ CW = −D 1 is valid in the range From this expression it is clear that for systems where both FM and AFM couplings are present, the CW behavior is valid at T |Θ CW |. For the unfrustrated STL-model (the first row in table 1), the condition (9) reads T 70 K, while for the prior INS based model (the last row in table 1) it is T 120 K. One should also take into account that the convergence of the HTE is slow. To show that the exchange values determined from the INS are compatible with the χ(T ) data, we reproduce in figure 1 the data measured on the same sample (i.e. with the same O content) used for the INS studies (see figure  4 of [6]). We have fitted the data in the range 340 < T < 380 K with the expression where χ 8 (T ) is the eighth-order HTE obtained by the method and programs published in [25] ‡, and N A is Avogadro's number. The HTE program of [25] we have used here can only treat systems with four different exchange couplings,only, so we decided to adopt one effective coupling in the a direction J a = J 3 + 2J 4 . Only one parameter g a was adjusted during the fit despite the small background susceptibility χ 0 which was set to zero for the sake of simplicity and its insensitivity to our fit. The [4,4]-Padé approximant for the HTE fits the data well down to T ∼ 20 K with the reasonable value g a = 2.34.
The inset shows C * (T ) and Θ * CW (T ), the two parameters of the pseudo-CW-law (cf. figure 4 of [6]), which is given by a tangent to the 1/χ(T ) curve at a given T . The pseudo-Curie "constant" C * (T ) exceeds its asymptotic value for all T . Hence, it can not be used to extract the number of spins in the system. The C from a more elaborate HTE is much better suited for this aim. For example, a clear crossover between a CW (without any divergence at Θ CW =+39 K) and a pseudo-CW regimes has been observed recently near 150 K for CuAs 2 O 4 , see figures 8 and 9 in [26].
We note that the difficulties encountered at intermediate T when avoiding the use of a HTE as demonstrated above can be circumvented by applying a strong magnetic field (up to saturation at H s where the AFM IC is suppressed). At very low T and in the isotropic limit, one has a very simple but useful relation The quantities on the r.h.s. can be deduced from the INS data [6,17]. Indeed, with J 1 = −230 K, α = 0.326, H s =55.4 T for H||b, g b = 2.047 [27], and J a ≈ 8 K, we arrive nearly at the same Θ CW = 54 K as in the HTE analysis of χ(T ). For collinear systems, H s is related to the inter-sublattice couplings (see equation (1) in [27], valid at T =0) where N IC =8 is the number of nearest IC neighbors, J IC =J 6 , and J IC =J 5 . Equation (12) gives (J IC + J IC )/k B =9.5 K, close to the value from our model. Conversely, equation  [7], which were explained within a pd-model consistent with the INS data. The detection of an intrachain Zhang-Rice singlet exciton at T > T N and significant increase of its weight at 300 K evidence the increase of AFM correlations in the chain and cannot be explained without a frustrating AFM intrachain coupling J 2 . Next, we remind the reader of some microscopic insights into the magnitude of the exchange couplings for Li 2 CuO 2 . Empirically, a dominant FM J 1 [6,17] and the minor role of J a [6,16,17] follow directly from the observed strong (weak) dispersion of the magnons along the b(a)-axis, respectively. Microscopically, the large value of J 1 as compared with -100 K suggested in [28] is a consequence of an enlarged direct FM exchange K pd between two holes on Cu and O sites within a CuO 4 plaquette and the non-negligible Hunds' coupling between two holes in different O 2p x and 2p y orbitals in edge-sharing CuO 2 chain compounds. The former is well-known from quantum chemistry studies (QCS) of corner-sharing 2D cuprates having similar CuO 4 plaquettes with K pd,π ≈ −180 meV for Φ = 180 • Cu-O-Cu bond angles [29]. Now we confirm similar values of the two K pd interactons in Li 2 CuO 2 using QCS [30], which yielded for Φ ≈ 90 o : K pxd ≈ K pxd,π/2 = K pyd,π/2 ≈ K py,d ∼ −100 meV within K pxd + K pyd < K pd,π < − K 2 pxd,π/2 + K 2 pyd,π/2 , (ignoring a weak crystal field splitting of the two O onsite energies). Thus, our two K pxd = K pyd values being naturally slightly different are more than twice as large as -50 meV for both direct FM couplings ad hoc adopted in [12,28] in the absence of QCS calculations for Li 2 CuO 2 and CuGeO 3 [12,31] 20 years ago. Concerning J a , QCS yield a small direct FM dd contribution ≈ -1 K, while our DFT-derived Wannier-functions point to a tiny AFM value of about 5 K. A more direct calculation within the LDA+U scheme for large supercells is in preparation. Anyhow, even with possible error bars of the DFT and QCS calculations, a value of J a exceeding ±10 K is not expected. Hence, a much larger J a ≈ −103 K as suggested by Shu et al. is either an artifact from the constructed "simplified" non-frustrated STL model and/or a consequence of the improperly analyzed χ(T )-data, as explained above.
Thus, we have shown that the "effective" unfrustrated STL-model proposed by Shu et al. together with its doubtful analysis of χ(T ) are at odds with our current understanding of real Li 2 CuO 2 and will only confuse readers. It may, however, be of some academic interest to illustrate the crucial effect of intra-and IC frustrations evidenced by a recent INS study. In case of finite NNN exchange J 2 with α > α c =1/4, as evidenced by the minimum of the magnon dispersion near the 1D propagation vector, there would be an inphase or antiphase ordering of non-collinear spirals for FM (AFM) perpendicular IC J [32]. Such a situation would also be realized for an STL-model corrected by α > α c , in sharp contrast to the observed collinear magnetic moments aligned ferromagneticaly along the chain direction with weak quantum effects. In this context, we note that we know of no example of a frustrated model that can be replaced comprehensively by a non-frustrated effective model as Shu et al. have tried to do. We admit that some selected special properties like the long-wave susceptibility at q = 0, if treated properly by a HTE, could be approximated this way. But quantities (e.g. χ(q, ω)) with the q-wave vector in the vicinity of the 1D 'propagation vector' defined by the maximum of the 1D magnetic structure factor S(q, 0)) certainly do not belong to them. This is clearly evidenced by the magnon minimum (soft spin excitation above about 2.5 meV, only, from the ground state) as seen in the INS [6] reflecting the vicinity to a critical point. Thus, the low-energy dynamics of both models are different since the microscopic mechanisms and the stability of the FM inchain ordering are completely distinct for these two approaches.
We hope that our detailed discussion of χ(T ) at intermediate T is of interest to a broader readership. We have drawn attention to the available HTE computer codes. Thereby we stress that the usual smooth behavior of χ(T ) doesn't allow one to deduce multiple exchange parameters, especially if the maximum and the submaximum regions at low-T are not involved in the fit. Other thermodynamic quantities such as the saturation field, magnetization, specific heat, and INS data must be included. What can and must be done when dealing with χ(T ), is to check the compatibility of already existing derived or proposed parameter sets with experimental χ(T ) data. The model proposed by Shu et al. has not passed our χ-check. Their argument against our FM Θ CW ≈ 50 K [6] due to a seemingly "incomplete selection of fitting parameters" is obsolete since a weak J a value has been included in our refined analysis resulting now in Θ CW ≈ 55 K. The stressed seemingly divergence of χ(T ) at T = Θ CW ≈ 50 K of our frustrated set doesn't reflect "an inconsistency" of our approach as claimed by Shu et al..
Rather, it provides just a proof and documents only that the authors did not consider the restricted regions of validity for any CW-law approaching even intermediate and of course lower T ! In this context, we mention that, if the authors would have inserted their own correct sign values for Θ CW > 0 derived from its definition by equation (4) (see table 1), they would be confronted themselves with the same pseudo-problem they have tried to accuse us! Their referring to the work [33] § as a 'first principle density functional theory (DFT)' explanation for why a spiral ordering is missing, is somewhat misleading. There the relative weak frustrating AFM IC's and the related and observed [6] dynamical spiral fluctuations were ignored but Xiang et al. [33] arrived nevertheless at a non-negligible finite J 2 value above α c = 1/4, in qualitative accord with the INS data [6]. To the best of our knowledge, there is no edge-sharing chain cuprate with a vanishing J 2 and a finite frustrating AFM value is considered to be generic for the whole family. There is no reason to prescribe a seemingly vanishing J 2 value to the presence of a significant level of O vacancies in high-quality samples. The INS [6,17] and RIXS data [10] provide clear evidence for a skew AFM IC with no need for a speculative and exotic 'order from disorder' scenario.
To summarize, the frustrated model for Li 2 CuO 2 , refined by the INS data [6,16,17] (s. the last row in table 1), our DFT and five-band pd-Hubbard model calculations [6], are consistent with the experimental data, including the χ(T ) [6], the field dependent magnetization, the saturation field value [27], and the T -dependent RIXS spectra [7]. This is in contrast to the SLT-model, which does not account for these data.  (7) is wrong. In the proposed corrected equation, they use the notations of Ref. [18] for the exchange couplings, i.e. the exchange interaction of a pair of spins is defined as −2J * , in contrast to our (and that of [13]) notations, where it reads +J. As we have already mentioned above, a positive sign of J corresponds to an AFM interaction.
(ii) Unfortunately, in the Reply [38], Shu et al. repeat their erroneous statement that a positive value of Θ CW would be allowed only for a system with an FM ground state (see the end of p. 5). Just for illustration, let us consider the case of a simple cubic magnetic lattice with a large NN FM interaction along the x-direction and a tiny NN AFM coupling along the y-and z-directions. It is clear that the system will have a Θ CW > 0 and an AFM ground state.
(iii) On pages 5,6 of Ref. [38], Shu et al. erroneously claim that HTE (they call it "HTSE") would be an approximation to the Curie-Weiss law (see the paragraph below Fig. 3 of Ref. [38]). Moreover, they state: "Finally, the polynomial HTSE fitting can always achieve high goodness-of-fit for nearly any curve without abrupt turns via higher order parameter adjustment within a limited temperature range". We may repeat here only that the high-temperature series expansion (HTE or HTSE) is an exact expansion and that the Curie-Weiss law is given by the first two terms of it (see equation (27) of section IV.B in [21]). It is not a result of the mean-field approximation. The terms of any order of a HTE are calculated from the Heisenberg model and are unambiguously expressed via the exchange parameters without any adjustment.
(iv) In section 4.5 of [38], Shu et al. erroneously claim that i<j J ij S i S j = 1 2 i,j J ij S i S j and doubt that we properly interpret our exchange values. We may mention that in the notations of Eq. (1), the total width of the dispersion of spin excitations is W ≈ 2|J 1 | [6]. The recently found full magnetic dispersion for Ca 2 Y 2 Cu 5 O 10 [36,39] amounts W ≈ 55 meV. So, the J 1 /k B > 200 K value (that give Θ CW ≈ +80 K) which is proved experimentally.