Low-temperature structure anomalies in CuNCN. Manifestations of RVB phase transitions?

Among the materials of the MNCN series (M  =Mn, Fe, Co, Ni), the phase CuNCN is the most peculiar example [1]. In variance with other members of the series exhibiting more or less standard antiferromagnetic behavior, CuNCN is a temperature independent (quasi-Pauli, see below) paramagnet at room temperature and, at lower temperatures, switches to a gapped (quasi-Arrhenius) temperature dependence. Nonetheless, it is not metallic in the temperature range where the Pauli paramagnetism occurs, hence no metal–insulator transition can be responsible for the quasi-Arrhenius behavior. It also does not manifest any magnetic neutron scattering [2] so that no long-range magnetically ordered (LRMO) state appears to which one could ascribe the susceptibility drop. These findings brought us [3–6] to the idea that CuNCN may form resonating valence bond (RVB) phases of the Cu²⁺ local spin 1/2 which are unequivocally observed in the Pauli paramagnetic phase with use of EPR [4]. This idea allowed us [3–6] to explain the magnetic and polarized neutron experiments by further assuming that the RVB phases are those of an anisotropic triangular antiferromagnetic Heisenberg model evolving in the ab crystallographic plane (figure 1). Can one expect structural manifestations of the conjectured transition between the RVB phases? Yes, because instead of the typical decrease of the unit cell dimensions one detects an anomalous temperature dependence of the lattice parameter a in synchrotron experiments, to be related with the 1D-RVB to 2D-RVB phase transition of the anisotropic triangular Heisenberg model which opens the gap in the quasiparticle spectrum along the a direction [5, 6].

Zoom In Zoom Out Reset image size

Puzzlingly, the experiment [6] also shows some irregularities for the lattice parameter c at about the same temperature (80–100 K) as for the lattice parameter a, and there are more irregularities in either the a- or the c-direction at about 30 K (see figure 2).

Zoom In Zoom Out Reset image size

The hypothesis of the triangular lattice is also not easy to reconcile with the intuitive picture [7] of the most important couplings extending along a and c. This brings us to the idea of considering more antiferromagnetic couplings and of formulating the Heisenberg model with the Hamiltonian

$$\begin{eqnarray} &&\mathop{\sum }\limits _{\mathbf{r}}\mathop{\sum }\limits _{\tau }{J}_{\boldsymbol{\tau}}{\mathbf{S}}_{\mathbf{r}}{\mathbf{S}}_{\mathbf{r}+\tau }, \end{eqnarray} \tag{ 1 }$$

where the translation vector τ takes four values τ_i,i = 1–4: τ₁ = (a,0);τ₂ = (0,c);τ₃ = (a,c);τ₄ = (a, − c), with an interaction of strength J_a along the lattice vector τ₁ (two neighbors), an interaction of strength J_c along the lattice vector τ₂ (two neighbors as well), and an interaction of strength J_ac along the lattice vectors τ₃ and τ₄ (two neighbors along each)—see figure 1. Either interactions along τ₁ and τ₂ or those along τ₃ and τ₄ taken separately must lead to an antiferromagnetic state. However, when considered simultaneously they interfere leading to a frustration not allowing the spins to arrange in any LRMO state. For similar systems a variety of RVB states have been proposed [8, 9]. Ground states of a similar, but spatially isotropic J₁J₂J₃ model have been treated recently by various methods, and it has been shown that spin-liquid (RVB) states are very probable [10]. This generalizes the known Nersesyan–Tsvelik model [11] which derives from ours by setting J_a = J_c, which is clearly not the case for CuNCN. In this paper, we consider in detail the RVB states of the above model in the mean-field approximation and apply this approximation to analyze the experimental data so far obtained for CuNCN.

2.1. Quasiparticle spectrum

Following the method [12] we use the fermion (spinon) representation of the spins,

$$\begin{eqnarray} &&{\mathbf{S}}_{i}=\frac{1}{2}{c}_{i \alpha }^{\dagger }{\boldsymbol{\sigma}}_{\alpha \beta }{c}_{i \beta }, \end{eqnarray} \tag{ 2 }$$

where ${c}_{i \sigma }^{\dagger }({c}_{i \sigma })$ are the fermion creation (annihilation) operators; σ_αβ are the elements of the Pauli matrices and summation over repeating indices is assumed. Applying the technique as described in appendix A (this method widely applies to various layered cuprates and is a standard treatment) we arrive at the quasiparticle spectrum,

$$\begin{eqnarray} &&\fl {E}_{\mathbf{k}}^{2}=9({J}_{a}^{2}{\zeta }_{a}^{2}{\cos }^{2}x+{J}_{c}^{2}{\zeta }_{c}^{2}{\cos }^{2}z+4{J}_{a c}^{2}{\zeta }_{a c}^{2}{\cos }^{2}x \hspace{0.167em} {\cos }^{2}z), \end{eqnarray} \tag{ 3 }$$

where we set x = k_x,z = k_z and introduce effective order parameters (OPs) ζ_a,ζ_c, and ζ_ac describing the states of the model.

The spectrum (3) is depicted in figure 3. It shows three principal regimes: (i) one with two pairs of lines of nodes (gapless 2D-RVB), (ii) one with a pair of lines of nodes (termed as 1D- and Q1D-RVB states), and (iii) one with two pseudogaps and four nodal points (pseudogapped 2D-RVB)—for the naming of the states see below.

Zoom In Zoom Out Reset image size

If either of the OPs ζ_a or ζ_c is the only nonvanishing OP, the quasiparticle spectrum acquires corresponding lines of nodes $z=\pm \frac{\pi }{2}$ (or $x=\pm \frac{\pi }{2}$) where the quasiparticles have zero energy. We call these states one-dimensional (1D-RVB) states since the dispersion of quasiparticles occurs in only one crystallographic direction (a or c—see, however, below). The quasiparticle density of states (qDoS) in the 1D-RVB states is constant at zero energy, but it diverges at the ceiling of the quasiparticle band due to the dispersionless ridge in the spectrum [3], similarly to the 1D-RVB state of the anisotropic triangular Heisenberg model found in [13]. If both OPs ζ_a,c vanish and the OP ζ_ac does not, two pairs of nodal lines exist along which the quasiparticles have zero energy. In this state the qDoS diverges logarithmically at zero energy.

If either of the nonvanishing OPs ζ_a,c is complemented by the nonvanishing OP ζ_ac, quasi-1D-RVB (Q1D-RVB) states appear. The difference from the true 1D-RVB states is that in the Q1D-RVB states one finds a nonvanishing dispersion transverse to the node lines, so that the quasiparticle spectra have local maxima and saddle points instead of a ridge. Thus, the qDoS develops a finite hop at the ceiling of the quasiparticle band and a logarithmic singularity at the pseudogap. By contrast, if the nonvanishing OPs are ζ_c and ζ_a, then (irrespective of the OP ζ_ac) there are no lines of nodes, but four nodal points ($\mathbf{k}=(\pm \frac{\pi }{2},\pm \frac{\pi }{2})$) of vanishing quasiparticle energies. The two possible states of this type are called 2D-RVB. The qDoS in 2D-RVB states vanishes at zero energy, being proportional to the energy well below the smaller pseudogap. Otherwise the quasiparticle dispersion law has saddle points at the two pseudogap energies and, thus, the qDoS of the 2D-RVB state develops two logarithmic singularities at the corresponding pseudogaps. In the Nersesyan–Tsvelik model (J_a = J_c) two logarithmic peaks coalesce and only one pseudogap is manifested. The exact forms of the qDoS of all RVB phases of the model are given in table C.1, appendix C. In the last column of table C.1 we show the analytical forms of the qDoS characteristic for the respective forms of the quasiparticle spectrum. It turned out quite unexpectedly that these qDoS can be found analytically for all phases of the c–a–ca-RVB model. Leaving the details of the derivation for a further publication we provide its sketch in appendix C.

2.2. Free energy and phase diagram of the model

The thermal behavior of the c–a–ca model derives from its free energy which can be written immediately [9] as

$$\begin{eqnarray} F&=&3{J}_{a}{\zeta }_{a}^{2}+3{J}_{c}{\zeta }_{c}^{2}+6{J}_{a c}{\zeta }_{a c}^{2}\nonumber\\ &&-~\frac{2 \theta }{4{\pi }^{2}}\int _{\mathrm{BZ}}\ln \left (2\cosh \left (\frac{{E}_{\mathbf{k}}}{2 \theta }\right )\right ){\mathrm{d}}^{2}\mathbf{k}, \end{eqnarray} \tag{ 4 }$$

where θ = k_BT and BZ stands for integration over the Brillouin zone. The minima of the free energy (4) with respect to the ζs correspond to possible phases of the system. The study of the ground state (zero-temperature limit) of the present model will be published elsewhere [14]. Here we focus on the results which can be obtained with use of the high-temperature expansion (technicalities are explained in appendix B). Minimizing the high-temperature expansion of (4) with respect to the OPs, we obtain the parameter phase diagram shown in figure 4 for the triple of the reduced exchange parameters subject to the condition ${J}_{a}^{\ast }+{J}_{c}^{\ast }+{J}_{a c}^{\ast }=1$ and a series of reduced temperatures θ* between 0.4 and 0.01 (the fractions of the above sum of the reduced exchange parameters are meant). Eight phases are possible. For the temperature above either of three critical ones,

$$\begin{eqnarray} &&{\theta }_{\tau }^{\mathrm{crit}}=\frac{3}{8}{J}_{\tau }, \end{eqnarray} \tag{ 5 }$$

the 'gray' phase persists in which all three OPs are equal to zero.

Zoom In Zoom Out Reset image size

Below these temperatures, gapless phases appear in respective corners of the parameter phase diagram. There are two 1D-RVB phases with OPs ζ_c or ζ_a ≠ 0 (red and green areas in figure 4). The gapless 2D-RVB phase with only one OP ζ_ac ≠ 0 occupies the blue area in figure 4.

Below the critical temperatures

$$\begin{eqnarray} &&{\theta }_{\tau \rightarrow \tau ,{\tau }^{\prime}}^{\mathrm{crit}}=\frac{{J}_{{\tau }^{\prime}}}{8}\left (1-\frac{2{J}_{{\tau }^{\prime}}}{3{J}_{\tau }}\right )^{-1}=\frac{3{J}_{\tau }{J}_{{\tau }^{\prime}}}{8\left (3{J}_{\tau }-2{J}_{{\tau }^{\prime}}\right )} \end{eqnarray} \tag{ 6 }$$

the corresponding phases with two nonvanishing OPs ζ_τ and ζ_τ' appear (the notation refers to a transition from the state with one nonvanishing OP, ζ_τ ≠ 0, to a state with two nonvanishing OPs, ζ_τ,ζ_τ' ≠ 0). The phases with two nonvanishing OPs are different as well. If one of the nonvanishing OPs is ζ_ac (magenta and cyan areas in figure 4) these are Q1D-RVB phases. If the two nonvanishing OPs are ζ_c and ζ_a, a 2D-RVB phase (orange) appears from the 1D-RVB phases (red and green).

The phase with three nonvanishing OPs (yellow) can only appear below the octal point ${J}_{a}^{\ast }={J}_{a c}^{\ast }={J}_{c}^{\ast }=\frac{1}{3};{\theta }^{\ast }=1/8$ where the phase with three vanishing OPs (gray) completely disappears and shows up from the above Q1D-RVB phases (magenta or cyan) at the critical temperatures

$$\begin{eqnarray} &&{\theta }_{\tau ,a c\rightarrow a,c,a c}^{\mathrm{crit}}=\frac{3{J}_{a c}{J}_{\bar {\tau }}}{8\left (3{J}_{a c}-2{J}_{\bar {\tau }}\right )}, \end{eqnarray} \tag{ 7 }$$

($\bar {\tau }=c,a$ for τ = a,c). It is transient and exists only above the critical temperature of

$$\begin{eqnarray} &&{\theta }_{a,c,a c}^{\mathrm{crit}}=\frac{3{J}_{a}{J}_{a c}{J}_{c}}{8\left (3{J}_{a}{J}_{a c}-5{J}_{a}{J}_{c}+3{J}_{a c}{J}_{c}\right )}, \end{eqnarray} \tag{ 8 }$$

where it switches to the 2D-RVB phase with two nonvanishing OPs (orange). The only difference between the dispersion laws in these two 2D-RVB phases is somewhat more pronounced dispersion along a 'ridge' in the case of the phase with three nonvanishing OPs. Otherwise both pseudogapped 2D-RVB phases have a qDoS with two van Hove singularities at the energies of their characteristic pseudogaps, and their physics has to be pretty similar. All phase transitions are of second order, that is to say that the OPs split from zero continuously at the corresponding transition temperatures. Obviously, in the high-temperature approximation only the classical value of the critical exponent $\beta =\frac{1}{2}$ can appear.

In general, one can expect that the high-temperature expansion will be valid down to temperatures of the order of J. This is, however, a too conservative estimate. It had been shown in [9] that the high-temperature estimate for the critical temperature of the transition between the Curie and Pauli paramagnetic phases of the anisotropic triangular Heisenberg model, which coincides with (5), precisely reproduces the numerical result, showing by this that the high-temperature expansion in a frustrated system remains valid for temperatures reasonably below 0.5J. Furthermore, the lower critical temperature given by (6) as well is in fair agreement with the numerical result of [12] for the anisotropic triangular Heisenberg model as shown in [3, 5]. This shows that the validity range of the high-temperature expansion in fact extends down to about 0.1J', where J' < J is a smaller (oblique) exchange parameter of the anisotropic triangular Heisenberg model. Thus at least in the most studied frustrated model the number and approximate positions of critical points are semiquantitatively reproduced by the high-temperature expansion which in that sense remains valid down to a pretty low temperature. This turns out to be true for a frustrated system—where the motions remain chaotic due to the frustration. That is why we cautiously rely on this method for our new frustrated Hamiltonian.

In order to make the not particularly intuitively clear minima of the free energy (4) more comprehensible we on the basis of the analysis of [9] represent pictorially the characteristic microstates contributing to the respective phases in figure 5. According to [9] the RVB wavefunction is a superposition (linear combination/resonance) of microstates written with respect to individual spins residing in the vertices of a lattice. In the case of the c–a–ca model these are the vertices of a rectangular lattice with the constants a (horizontal in figure 5) and c (vertical in figure 5). In the microstates all spins form singlet pairs (valence bonds) of the form $\frac{1}{\sqrt{2}}({c}_{\mathbf{r}\uparrow }^{\dagger }{c}_{{\mathbf{r}}^{\prime}\downarrow }^{\dagger }+{c}_{\mathbf{r}\downarrow }^{\dagger }{c}_{{\mathbf{r}}^{\prime}\uparrow }^{\dagger })$. The microstates are thus the products of valence bonds for different pairs of vertices rr' such that in a given microstate each spin enters in one and only one valence bond. Being summed with different amplitudes controlled in their turn by the values of the OPs symbolized by the ⊕ sign, the microstates produce specific RVB states. The states with a single nonvanishing OP are easiest to understand. These are the gapless 2D- and 1D-RVB states shown in two leftmost panels of figure 5. In these states the bonds appear only along the directions where the OPs are nonvanishing. Thus the 1D-RVB states differ, as one can imagine, by the directions (a or c) in which the bonds extend. Formally, being directed in, say, the c direction means that for all bonds the following holds: r = r' + nc, with an integer n, and analogously for the 1D-RVB state extended in the a direction. We use the color code of the respective phases in figure 4 to render the bonds extended in the corresponding directions. In the Q1D-RVB states the bonds extended in one of the lattice directions are complemented by those coupling the spins along the geometric diagonals: a ± c. In the (pseudo)gapped 2D-RVB states the bonds extend in either direction along the crystallographic axes. The 2D-RVB state with all three nonvanishing OPs (yellow) contains all possible bonds, whereas the 2D-RVB state with vanishing ζ_ac (orange) misses the bonds extended along the geometric diagonals, as can be seen in the two rightmost panels. The above considerations apply to the pure RVB states, i.e. at zero temperature. For the respective RVB phases at finite temperature the depicted ground states of each type are accompanied by excited states taken with respective statistical weights, in which some temperature dependent fraction part of the bonds is broken.

Zoom In Zoom Out Reset image size

As mentioned above and previously [3–5], the absence of magnetic scattering in CuNCN is perfectly explained by the hypothesis of the RVB character of its phases. The magnetic properties (qualitative behavior of the susceptibility as a function of temperature) as deriving from the characteristic features of the qDoS are designated in the first column of table C.1. It is not surprising that the temperature-independent paramagnetism (conditionally denoted as 'Pauli', although it goes about spinons rather than conductivity electrons) occurs in the (Q)1D-RVB phases with a constant qDoS at zero energy. The 'blue' gapless 2D-RVB phase with two pairs of intersecting node lines is manifested as a logarithmic singularity in the qDoS at zero energy. This singularity integrates and produces a logarithmic divergence of the susceptibility at zero temperature which qualitatively would be difficult to distinguish from the Curie paramagnetism. For two pseudogapped 2D-RVB phases (yellow and orange) with nodal points in the dispersion law and linear dependence of the qDoS in the low-energy range one has also to expect a linear temperature dependence of the susceptibility in the low-temperature region (well below the smaller pseudogap) superimposed with a quasi-Arrhenius behavior with characteristic energy of the pseudogap at higher temperature. This rich variety of possible phases in the c–a–ca-RVB model allows us to eventually explain the magnetic behavior of CuNCN. As previously, we assume that the (Q)1D-RVB phase sets in at some fairly high temperature which cannot be directly checked due to decomposition of the material [4]. Incidentally, Curie paramagnetism is not observed in CuNCN in the range of its thermal stability. The same applies to the gapless 2D-RVB phase with logarithmically divergent susceptibility. The Pauli paramagnetism, which changes to the quasi-Arrhenius gapped regime, is thus logically attributable to the (Q)1D-RVB phases, transforming to the pseudogapped 2D-RVB phases at lower temperatures. By this, the entire picture of the magnetic properties of CuNCN is qualitatively reproduced.

The observed magnetic behavior qualitatively reconciles with the temperature course of the lattice parameters a and c (figure 2). The sequence of the parameter phase diagrams shows which specific RVB phases emerge at the specified temperatures in the respective areas of the parameter space. The CuNCN species is characterized by a specific set of parameter values derived in appendix D (see also below), mapping to one point in the parameter space marked by the star. At different temperatures various phases occupy slightly different areas of the parameter space. In the lower row of figure 4 we can see how the point representing CuNCN while staying in its original place effectively moves from one phase to another: subsequently changing from the Q1D-RVB (magenta) to the transient 2D-RVB (yellow) and finally to the low-temperature 2D-RVB (orange).

Below the octal point θ* = 1/8 (the lower row of ternary diagrams in figure 4), the thermal evolution consists en gros in squeezing out all the phases by the 2D-RVB phase with vanishing ζ_ac (orange). The 2D-RVB phase with three nonvanishing OPs (yellow) is transient, so its parameter area is never large and it is subject to deformations and displacements under the 'pressure' of the orange phase. The parameter estimate of appendix D positions CuNCN on the parameter phase diagram (figure 4) at the point with the reduced parameters (barycentric coordinates) ${J}_{a}^{\ast }=0.1 7 2,\hspace{0.167em} {J}_{c}^{\ast }=0.4 0 0$, i.e., within the magenta phase, but close to and somewhat below the quadruple point of the red, magenta, orange, and yellow phases in the leftmost graph in the lower row. In this row one can observe the following sequence of transitions between the RVB phases:

$$\begin{eqnarray} &&\begin{array}{@{}c@{}} \displaystyle \mathrm{Q1D}\text{-}\mathrm{RVB}\rightarrow \mathrm{2D}\text{-}\mathrm{RVB}\rightarrow \mathrm{2D}\text{-}\mathrm{RVB}\\ \displaystyle {\zeta }_{c},{\zeta }_{a c}\rightarrow {\zeta }_{c},{\zeta }_{a c},{\zeta }_{a}\rightarrow {\zeta }_{c},{\zeta }_{a}\\ \displaystyle \mathrm{magenta}\rightarrow \mathrm{yellow}\rightarrow \mathrm{orange} \end{array} \end{eqnarray} \tag{ 9 }$$

upon decreasing temperature. The thermal dependence of the OPs within these phases is described by the expressions given in table 1. One should not expect that these expressions derived from the high-temperature expansion for the free energy will be exactly valid at low temperature. Specifically, the prefactors θ should not be taken seriously since it has been shown [14] that the exact ζ_a,c OPs flow to some finite values as the temperature flows to zero. Thus, the OPs split from zero values at the critical temperatures given next to them. Since the exchange parameters satisfy the condition J_ac⪆J_c, the OP ζ_ac is always small and decreases with decreasing temperature, thus we do not consider it explicitly further. The OP ζ_c increases either in the magenta or in the yellow phase according to the same law (the reference/critical temperatures for this OP in these two phases are equal), but we assume that in the interesting temperature range $\theta \ll {\theta }_{c,a c}^{\mathrm{crit}}={\theta }_{a,a c\rightarrow a,c,a c}^{\mathrm{crit}}=4 3 5~\mathrm{K}$ it has almost reached its zero-temperature limit and does not significantly change any more¹. The OP ζ_a splits from zero at the critical temperature of ${\theta }_{c,a c\rightarrow a,c,a c}^{\mathrm{crit}}\approx 1 0 0~\mathrm{K}$ (transition to the yellow phase). Finally, at the lowest accessible critical temperature ${\theta }_{a,c,a c}^{\mathrm{crit}}$ (8) a transition to the orange phase takes place². The transition to the orange phase (evanescence of ζ_ac) at ${\theta }_{a,c,a c}^{\mathrm{crit}}$ whatever it is affects neither the bandwidth, since this OP is never large, nor the character of the temperature dependence of ζ_a, although some variation of its slope can be expected. At this phase transition, however, a remarkable change can be expected in the character of the temperature dependence of ζ_c. Namely, it switches from an increase to a decrease when the temperature decreases. This happens because of an instantaneous change of its reference temperature shown in table 1 from positive ${\theta }_{a,a c\rightarrow a,c,a c}^{\mathrm{crit}}$ in the yellow phase to ${\theta }_{a\rightarrow a,c}^{\mathrm{crit}}$ in the orange phase, which is negative in the relevant area (J_ac⪆J_c > J_a) of the exchange parameter space. Changing the sign of the reference temperature causes the change of the temperature course of the OP ζ_c which, in the orange phase, starts to decrease with temperature decrease instead of being constant (or very slowly growing) in the yellow phase.

Table 1.  Temperature dependences of the OPs for possible phases of CuNCN in the high-temperature approximation of the c–a–ca-RVB model. The areas of existence of the corresponding phases are given by the condition of positiveness of the expressions under the square roots.

Color code	OPs versus θ
Blue	$\left \{ \begin{array}{@{}ll@{}} \displaystyle {\zeta }_{c}={\zeta }_{a}=0 \\ \displaystyle {\zeta }_{a c}=\frac{4 \theta }{3{J}_{a c}}\sqrt{1-\theta /{\theta }_{a c}^{\mathrm{crit}}}&\displaystyle {\theta }_{a c}^{\mathrm{crit}}=3{J}_{a c}/8 \end{array}\right.$
Magenta	$\left \{ \begin{array}{@{}ll@{}} \displaystyle {\zeta }_{a}=0\\ \displaystyle {\zeta }_{c}=\frac{4 \theta }{3{J}_{c}}\sqrt{1-\theta /{\theta }_{c,a c}^{\mathrm{crit}}}&\displaystyle {\theta }_{c,a c}^{\mathrm{crit}}=\frac{3{J}_{a c}{J}_{c}}{8\left (3{J}_{a c}-2{J}_{c}\right )} \\ \displaystyle {\zeta }_{a c}=\frac{8}{3}\sqrt{\frac{2}{3}}\frac{\theta }{{J}_{a c}}\sqrt{\frac{\theta \left ({J}_{a c}-{J}_{c}\right )}{{J}_{a c}{J}_{c}}} \end{array}\right.$
Yellow	$\left \{ \begin{array}{@{}ll@{}} \displaystyle {\zeta }_{a}=\frac{4 \theta }{3{J}_{a}}\sqrt{1-\theta /{\theta }_{c,a c\rightarrow a,c,a c}^{\mathrm{crit}}}&\displaystyle {\theta }_{c,a c\rightarrow a,c,a c}^{\mathrm{crit}}=\frac{3{J}_{a c}{J}_{a}}{8\left (3{J}_{a c}-2{J}_{a}\right )} \\ \displaystyle {\zeta }_{c}=\frac{4 \theta }{3{J}_{c}}\sqrt{1-\theta /{\theta }_{a,a c\rightarrow a,c,a c}^{\mathrm{crit}}}&\displaystyle {\theta }_{a,a c\rightarrow a,c,a c}^{\mathrm{crit}}=\frac{3{J}_{a c}{J}_{c}}{8\left (3{J}_{a c}-2{J}_{c}\right )} \\ \displaystyle {\zeta }_{a c}=\frac{4}{3}\sqrt{\frac{1}{3}}\frac{\theta }{{J}_{a c}}\sqrt{\theta /{\theta }_{a,c,a c}^{\mathrm{crit}}-1}&\displaystyle {\theta }_{a,c,a c}^{\mathrm{crit}}\text{ given by }(8) \end{array}\right.$
Orange	$\left \{ \begin{array}{@{}ll@{}} \displaystyle {\zeta }_{a}=\frac{4 \theta }{\sqrt{1 5}{J}_{a}}\sqrt{1-\theta /{\theta }_{c\rightarrow a,c}^{\mathrm{crit}}}&\displaystyle {\theta }_{c\rightarrow a,c}^{\mathrm{crit}}=\frac{3{J}_{a}{J}_{c}}{8\left (3{J}_{c}-2{J}_{a}\right )} \\ \displaystyle {\zeta }_{c}=\frac{4 \theta }{\sqrt{1 5}{J}_{c}}\sqrt{1-\theta /{\theta }_{a\rightarrow a,c}^{\mathrm{crit}}}&\displaystyle {\theta }_{a\rightarrow a,c}^{\mathrm{crit}}=\frac{3{J}_{a}{J}_{c}}{8\left (3{J}_{a}-2{J}_{c}\right )} \\ \displaystyle {\zeta }_{a c}=0 \end{array}\right.$

Previously [5, 6], we were able to explain the temperature anomaly of the a lattice parameter by attributing it to opening a gap in the quasiparticle spectrum in the crystallographic direction a and formation of a gapped 2D-RVB state in the triangular anisotropic Heisenberg model. This is accompanied by the split of an RVB OP from zero, which is coupled to the lattice constant through 'magnetostriction' terms. Although the previous model (triangular anisotropic Heisenberg lattice) could not be fully substantiated, similar moves absolutely apply in the present model with three exchange parameters and respective OPs. We assume as previously a linear relationship between the exchange parameters J_τ and geometry parameters ρ_λ,

$$\begin{eqnarray} &&{J}_{\tau }={J}_{\tau 0}+\mathop{\sum }\limits _{\lambda }J_{\tau ,\lambda }^{\prime}{\rho }_{\lambda }. \end{eqnarray} \tag{ 10 }$$

Following [15] we assume that zero values of ρs correspond to a hypothetical structure, which the CuNCN crystal would have provided the exchange interactions J_τ were turned off. The observed geometry of the crystal where the exchange interactions result in the formation of one of the RVB phases corresponds to the minimum with respect to ρ_λ of the free energy (4) to which the elastic energy (E.1) is added, and the exchange parameters are replaced according to the recipe (10). Using the special form of the RVB free energy (4) we arrive at the explicit expression (for details see appendix E)

$$\begin{eqnarray} &&{\rho }_{\mu }=\mathop{\sum }\limits _{\lambda }(K)_{\mu \lambda }^{-1}\left (\mathop{\sum }\limits _{\tau }{A}_{\tau }\zeta _{\tau }^{2}{J}_{\tau ,\lambda }^{\prime}\right ), \end{eqnarray} \tag{ 11 }$$

the sought relation between the RVB OPs and their structure manifestations (A_τ are the numerical coefficients: 3 for τ = a,c; 6 for τ = ac) which further generalizes the famous bond-length–bond-order relation to the RVB states.

The temperature dependence of the structure described by the parameters ρ_λ is thus that of the relevant combination of the squares of the corresponding OPs. In the magenta and yellow phases the OP ζ_c is (almost) temperature independent, thus ${\zeta }_{c}^{2}$ does not contribute to the visible variation of the lattice constants; therefore we focus first on the contribution of the OP ζ_a (our previous consideration allows us to neglect the OP ζ_ac).

Before switching to rationalizing the low-temperature behavior of the lattice parameters shown in more detail in figure 6, we stress once again that our purpose is to provide a unified view of the data obtained in polarized neutron scattering [2], susceptibility [1, 4], and structure (synchrotron [6] and neutron [16]) experiments. Thus, in our reasoning we follow the above unifying concept in the hope that minor inconsistencies (possibly coming from the uncertainty in the measurements or sample irregularities) can be resolved in the future. In order to proceed we introduce geometry variables ρ_a = δa = (a − a₀), 2ρ_c = δc = (c − c₀), and assume that the exchange constant J_a is independent of the deformation in the c-direction (see below). Thus, ${J}_{a,c}^{\prime}=0$, and we apply the general relation (11). In this setting the variation of two lattice parameters is

$$\begin{eqnarray} &&\left (\begin{array}{@{}c@{}} \displaystyle \delta a\\ \displaystyle \delta c \end{array}\right )=\frac{3{J}_{a,a}^{\prime}{\zeta }_{a}^{2}}{\vert K\vert }\left (\begin{array}{@{}c@{}} \displaystyle {K}_{c c}\\ \displaystyle -{K}_{a c} \end{array}\right ), \end{eqnarray} \tag{ 12 }$$

where |K| is the determinant of the 2 × 2 matrix of the force constants. Clearly, the off-diagonal force constant K_ac makes the lattice parameter c sensitive to variations of the OP ζ_a. Moreover, the effect on c is predicted to have the opposite sign to that on a (a increases, c decreases, and vice versa). When the pseudogap opens in the a direction (about 100 K), two lattice parameters have to start changing in opposite directions, as observed. The same behavior can be traced in the neutron scattering data [16] which are going to be analyzed in detail in a forthcoming publication.

Zoom In Zoom Out Reset image size

The crucial point is the sign of the effect. At first glance the situation seems to be counterintuitive since in order to conform with the experimental increase of the a-parameter in the 2D-RVB phase with two pseudogaps (orange or yellow) as compared to the (Q)1D-RVB phase (red or magenta) the exchange parameter J_a has to increase with an increasing Cu–Cu interatomic separation. However, the effective value of the (antiferromagnetic) exchange parameter J_a is a sum of contributions of different signs,

$$\begin{eqnarray} &&{J}_{a}={J}_{a}(\mathrm{antiferro})-{J}_{a}(\mathrm{ferro})\gt 0, \end{eqnarray} \tag{ 13 }$$

where both J_a(antiferro) and J_a(ferro) are positive. The antiferromagnetic contribution is accumulated on numerous superexchange paths and is weakly affected by the Cu–Cu interatomic separation since no direct matrix elements between the states of the two Cu atoms affect J_a(antiferro) ³. By contrast, the ferromagnetic contribution strongly depends on the angle $\widehat{\mathrm{CuNCu}}$ (see [7, 17]), and in the range of $\widehat{\mathrm{CuNCu}}\gt 9 0^{\circ}$, which is the case for CuNCN, decreases, while the $\widehat{\mathrm{CuNCu}}$ angle and thus the Cu–Cu separation increases. Likewise, the overall effective antiferromagnetic exchange parameter J_a increases while the counterpoising ferromagnetic contribution decreases in absolute value. This explains the overall positive sign of ${J}_{a,a}^{\prime}$ and, thus, the experimental fact of increasing lattice parameter a in the phase with the pseudogap developing along the a-direction. Due the above sign relation between the effect of the a-pseudogap opening on the a- and c-lattice parameters, the experimentally observed extra decrease of the c lattice parameter at 100 K is explained as well.

The available data are insufficient to solve the subtle question of the critical exponents of the OPs. However, assuming the classical value of the critical exponent $\beta =\frac{1}{2}$ prescribed by the high-temperature expansion as a plausible interpolation in a way confirmed by fitting the susceptibility run in the 2D-RVB phases (see appendix D), the lattice parameters must obey the linear temperature evolution, i.e., one with the exponent 2β. The latter is approximately observed for the lattice parameter a in a wide range below the (Q)1D- to 2D-RVB transition at about 100 K which, however, changes to a somewhat more chaotic behavior at about 30 K. It incidentally corroborates with the temperature course of the c lattice parameter measured in the synchrotron experiment [6] (figure 6(b)) which has not received due attention so far. As one can see, the lattice parameter c generally decreases with temperature. This behavior is in line with the usual thermal expansion and thus does not require any special explanation. It shows noticeable irregularities (not changing the sense of the course) at about 80–100 K, where the magnetic susceptibility changes from the quasi-Pauli to the quasi-Arrhenius regime, which has been explained above through its interplay with the opening of the pseudogap in the a-direction right in this temperature range. The sense of the temperature course of c changes at 30 K where the lattice parameter a starts to show irregularities in its turn. This, however, nicely fits in the above picture of the RVB phase transitions. Since the available amount of experimental data does not suffice to reliably estimate quite a few magnetostriction and other parameters required for the quantitative description, we try to derive a qualitative description in terms of the slopes of the lattice parameters versus temperature at critical points. In the minimal setting (i.e. neglecting the off-diagonal force constants and magnetostriction terms) we immediately obtain for the slope of c versus θ

$$\begin{eqnarray} &&{\frac{\mathrm{d}\delta c}{\mathrm{d}\theta }\biggr \vert }_{\theta \rightarrow {\theta }^{\mathrm{crit}}}=\frac{6{J}_{c,c}^{\prime}{K}_{a a}}{\vert K\vert }\left [{\zeta }_{c}\left (\theta \rightarrow 0\right )\right ]^{2}\left (-\frac{1}{{\theta }^{\mathrm{crit}}}\right ). \end{eqnarray} \tag{ 14 }$$

At the transition from the yellow to the orange phase this slope has to change its sign since that of θ^crit instantaneously changes from positive in the yellow to negative in the orange phase, as explained above, and that is precisely what happens. We also conclude that the magnetostriction parameter ${J}_{c,c}^{\prime}$ has an intuitively understandable sign: it is negative since this exchange parameter is dominated by the electron transfers through the π-systems of the NCN²⁻ units which in their turn strongly depend on the Cu–N interatomic separations affected by the deformation in the c-direction (see footnote 3). The temperature slope must therefore be negative in the orange (low-temperature) phase and positive in the yellow (transient) one, and that is again what happens.

In a more general setting the contribution of the OP ζ_c to the slopes of the lattice parameters versus temperature is

$$\begin{eqnarray} &&\frac{\mathrm{d}}{\mathrm{d}\theta }\left .\left (\begin{array}{@{}c@{}} \displaystyle \delta a\\ \displaystyle \delta c \end{array}\right )\right \vert _{\theta \rightarrow {\theta }^{\mathrm{crit}}}&=&\frac{3\left [{\zeta }_{c}\left (\theta \rightarrow 0\right )\right ]^{2}}{\vert K\vert }\left (-\frac{1}{{\theta }^{\mathrm{crit}}}\right )\nonumber\\ &&\times ~\left (\begin{array}{@{}c@{}} \displaystyle \left ({K}_{c c}{J}_{c,a}^{\prime}-2{K}_{a c}{J}_{c,c}^{\prime}\right )\\ \displaystyle \left (-{K}_{a c}{J}_{c,a}^{\prime}+2{K}_{a a}{J}_{c,c}^{\prime}\right ), \end{array}\right ) \end{eqnarray} \tag{ 15 }$$

so that it is clear that the contribution of the variation of the OP ζ_c to the slope of a versus θ changes its sign when going from the yellow to the orange phase. The visible decrease of the negative slope of a in absolute value as seen below 30 K (in the orange phase as compared to the yellow phase) requires that the addition to it coming from ζ_c is positive. Semiquantitatively, this requires that $\left ({K}_{c c}{J}_{c,a}^{\prime}-2{K}_{a c}{J}_{c,c}^{\prime}\right )\gt 0$ since in the orange phase θ^crit in the above expression equals the negative ${\theta }_{a\rightarrow a,c}^{\mathrm{crit}}$. Incidentally, the term $-2{K}_{a c}{J}_{c,c}^{\prime}(\gt 0)$ has the required sign as estimated above. Thus, ${J}_{c,a}^{\prime}$ has to be either positive, which is in principle conceivable as the above example of ${J}_{a,a}^{\prime}$ shows, or at least not too strongly negative. This tentatively explains the anomalies observed in figure 6.

A new form of the frustrated spatially anisotropic antiferromagnetic Heisenberg Hamiltonian close to the popular J₁J₂J₃ model with exchange parameters J_c, J_a, and J_ac extending along the c,a, and a ± c directions of a two-dimensional rectangular lattice is proposed. When applied to model the fascinating physics of copper carbodiimide (CuNCN) it explains the absence of magnetic order in CuNCN down to 4 K by assuming a resonating valence bond (RVB) character of its emerging phases. The quasiparticle spectrum of the RVB model of the proposed Hamiltonian shows three principal regimes: (i) a state with two pairs of lines of nodes, (ii) states with a pair of lines of nodes (termed as 1D- and Q1D-RVB states), and (iii) states with two pseudogaps and four nodal points (2D-RVB states). The extraordinarily rich parameters–temperature phase diagram of the model contains eight different phases whose magnetic behavior includes Curie and quasi-Pauli paramagnetism (1D- and Q1D-RVB phases), and (pseudo)gapped (quasi-Arrhenius) paramagnetism (2D-RVB phases). Both the quasi-Pauli and quasi-Arrhenius regimes are reliably observed in experiment. By adding magnetostriction and elastic terms to the free energy of the model, the temperature dependence of the CuNCN crystal structure is explained through a sequence of transitions between different RVB phases when the temperature decreases. Confronting the model with magnetic susceptibility and structure data recorded as functions of temperature between about 10 and 200 K, one finds remarkably good agreement between the theoretical predictions and the experiments.

Although structural manifestations are well known to accompany various magnetic transitions, e.g. [19], the specific result of this work is the unified description of the somewhat contradictory combination of absence of magnetic order, nongapped (Pauli) and gapped (Arrhenius) temperature courses of the magnetic susceptibility in the respective temperature ranges and structural effects coupled to the latter. The RVB hypothesis has a crucial importance not for the structural effects themselves, although we concentrate on them in this paper, but as an 'integrating factor' allowing us to explain from a unified viewpoint the facts coming from different experimental methods.

The spin-Peierls state in the c-direction could be a good ground state candidate within the picture of a strong exchange extended in this direction mediated by the NCN²⁻ units. Had the spin-Peierls be there it would explain at least two key features of CuNCN: the activation dependence of the susceptibility due to opening of the spin-Peierls gap and the absence of magnetic long-range order visible in the neutron experiments due to the formation of alternating Cu–Cu spin pairings along the c-direction. Unfortunately, this intellectually attractive picture has to be abandoned for the time being since no signature of superstructure in either direction could be traced in the available structural data.

This work has been performed with the support of the DFG. Dr Andrej Zorko of Jožef Stefan Institute (Ljubljana, Slovenia) is acknowledged for sending numerical data on the ESR susceptibility [4]; Dr Wouter van Beek of Swiss–Norwegian Beamlines at ESRF and Professor Dr Uwe Ruschewitz of the Institute of Inorganic Chemistry, University of Cologne are acknowledged for sending the structure data of [6] in numerical form. The referees are acknowledged for instructive comments and suggestions.

The equations of motion are based on the Heisenberg representation in which each operator obeys

$$\begin{eqnarray} &&\mathrm{i}\hbar \dot {A}=\left [A,H\right ], \end{eqnarray} \tag{ A.1 }$$

where [...,...] stands for the commutator and the '$\dot {~}$' for the time derivative. Applying this to the creation and annihilation operators ${c}_{\mathbf{r}\sigma }^{\dagger }({c}_{\mathbf{r}\sigma })$ and performing commutation, mean-field decoupling and Fourier transformation as done previously [3] results in their mean-field equations of motion,

$$\begin{eqnarray} \mathrm{i}\hbar {\dot {c}}_{\mathbf{k}\sigma }&=&-\frac{3}{2}\mathop{\sum }\limits _{\tau }{J}_{\boldsymbol{\tau}}{\xi }_{\boldsymbol{\tau}}\cos (\mathbf{k}\tau ){c}_{\mathbf{k}\sigma }\nonumber\\ &&-~\frac{3}{2}\mathop{\sum }\limits _{\boldsymbol{\tau}}{J}_{\boldsymbol{\tau}}{\Delta }_{\tau }\cos (\mathbf{k}\tau )c_{-\mathbf{k}-\sigma }^{\dagger },\nonumber\\ \mathrm{i}\hbar \dot {c}_{\mathbf{ k}\sigma }^{\dagger }&=&\frac{3}{2}\mathop{\sum }\limits _{\tau }{J}_{\boldsymbol{\tau}}{\xi }_{\boldsymbol{\tau}}\cos (\mathbf{k}\tau )c_{\mathbf{k}\sigma }^{\dagger }\nonumber\\ &&-~\frac{3}{2}\mathop{\sum }\limits _{\boldsymbol{\tau}}{J}_{\boldsymbol{\tau}}\Delta _{\tau }^{\ast }\cos (\mathbf{k}\tau ){c}_{-\mathbf{k}-\sigma }. \end{eqnarray} \tag{ A.2 }$$

The latter reduce to the set of 2 × 2 eigenvalue problems for each wavevector k,

$$\begin{eqnarray*} \left (\begin{array}{@{}cc@{}} \displaystyle {\xi }_{\mathbf{k}}&\displaystyle {\Delta }_{\mathbf{k}}\\ \displaystyle {\Delta }_{\mathbf{k}}^{\ast }&\displaystyle -{\xi }_{\mathbf{ k}} \end{array}\right )\left (\begin{array}{@{}c@{}} \displaystyle {u}_{\mathbf{k}}\\ \displaystyle {v}_{\mathbf{k}} \end{array}\right )={E}_{\mathbf{k}}\left (\begin{array}{@{}c@{}} \displaystyle {u}_{\mathbf{k}}\\ \displaystyle {v}_{\mathbf{k}} \end{array}\right ), \end{eqnarray*}$$

with

$$\begin{eqnarray*} &&{\xi }_{\mathbf{k}}=-3\mathop{\sum }\limits _{\tau }{J}_{\boldsymbol{\tau}}{\xi }_{\tau }\cos (\boldsymbol{k\tau}),\nonumber\\ &&{\Delta }_{\mathbf{k}}=3\mathop{\sum }\limits _{\tau }{J}_{\boldsymbol{\tau}}{\Delta }_{\tau }\cos (\boldsymbol{k\tau}). \end{eqnarray*}$$

The summation over τ extends to τ_i;i = 1–4, which results in the eigenvalues

$$\begin{eqnarray*} &&{E}_{\mathbf{k}}=\pm \sqrt{{\xi }_{\mathbf{k}}^{2}+\left \vert {\Delta }_{\mathbf{k}}\right \vert ^{2}}, \end{eqnarray*}$$

whose eigenvectors are combinations of the annihilation and creation operators with the coefficients u_k,v_k. This set of equations is closed by the self-consistency conditions

$$\begin{eqnarray} &&{\xi }_{\boldsymbol{\tau}}=-\frac{1}{2 N}\mathop{\sum }\limits _{\mathbf{k}}\exp (\mathrm{i}\mathbf{k}\tau )\frac{{\xi }_{\mathbf{k}}}{{E}_{\mathbf{k}}}\tanh \left (\frac{{E}_{\mathbf{k}}}{2 \theta }\right ),\nonumber\\ &&{\Delta }_{\boldsymbol{\tau}}=\frac{1}{2 N}\mathop{\sum }\limits _{\mathbf{k}}\exp (-\mathrm{i}\mathbf{k}\tau )\frac{{\Delta }_{\mathbf{k}}}{{E}_{\mathbf{k}}}\tanh \left (\frac{{E}_{\mathbf{k}}}{2 \theta }\right ) \end{eqnarray} \tag{ A.3 }$$

for the order parameters (OPs). The lattice-symmetry considerations allow us to restrict ourselves to the OPs: ξ_a,ξ_c,ξ_ac;Δ_a,Δ_c,Δ_ac. For the complex OPs we introduce a polar representation,

$$\begin{eqnarray} &&{\Delta }_{\boldsymbol{\tau}}={\eta }_{\boldsymbol{\tau}}{\mathrm{e}}^{\mathrm{i}{\varphi }_{\tau }}. \end{eqnarray} \tag{ A.4 }$$

The standard moves foreseen for the SU(2) symmetric solutions are used to exclude the cross terms in the OPs from ${E}_{\mathbf{k}}^{2}$, which leads to the system of three equations

$$\begin{eqnarray} &&{\xi }_{\boldsymbol{\tau}}{\xi }_{{\boldsymbol{\tau}}^{\prime}}=-{\eta }_{\boldsymbol{\tau}}{\eta }_{{\boldsymbol{\tau}}^{\prime}}\cos ({\varphi }_{\tau }-{\varphi }_{{\tau }^{\prime}}). \end{eqnarray} \tag{ A.5 }$$

Introducing the relative phases as θ_a = φ_a − φ_ac, θ_c = φ_c − φ_ac, we arrive at equations of the form⁴

$$\begin{eqnarray*} &&{\xi }_{a}{\xi }_{c}+{\eta }_{a}{\eta }_{c}\cos ({\theta }_{a}-{\theta }_{c})=0,\nonumber\\ &&{\xi }_{a}{\xi }_{a c}+{\eta }_{a}{\eta }_{a c}\cos {\theta }_{a}=0,\nonumber\\ &&{\xi }_{c}{\xi }_{a c}+{\eta }_{c}{\eta }_{a c}\cos {\theta }_{c}=0, \end{eqnarray*}$$

similar to those derived in [8] which can be satisfied, e.g., by setting

$$\begin{eqnarray} &&\begin{array}{@{}c@{}} \displaystyle {\xi }_{a}{\xi }_{c}=-{\eta }_{a}{\eta }_{c}\not = 0,\\ \displaystyle {\theta }_{a}={\theta }_{c}=\frac{\pi }{2},\\ \displaystyle {\eta }_{a c}\not = 0,\tqs {\xi }_{a c}=0. \end{array} \end{eqnarray} \tag{ A.6 }$$

Of course, the OPs can also be vanishing.

Under the above conditions the spectrum of quasiparticles acquires the form

$$\begin{eqnarray} {E}_{\mathbf{k}}^{2}&=&9({J}_{a}^{2}({\xi }_{a}^{2}+{\eta }_{a}^{2})\hspace{0.167em} {\cos }^{2}x+{J}_{c}^{2}({\xi }_{c}^{2}+{\eta }_{c}^{2})\hspace{0.167em} {\cos }^{2}z\nonumber\\ &&+~4{J}_{a c}^{2}({\xi }_{a c}^{2}+{\eta }_{a c}^{2}){\cos }^{2}x~{\cos }^{2}z), \end{eqnarray} \tag{ A.7 }$$

where we set x = k_x;z = k_z. Obviously, the solutions of the problem depend only on the effective OPs ${\zeta }_{a}=\sqrt{{\xi }_{a}^{2}+{\eta }_{a}^{2}}$, ${\zeta }_{c}=\sqrt{{\xi }_{c}^{2}+{\eta }_{c}^{2}}$, and ${\zeta }_{a c}=\sqrt{{\xi }_{a c}^{2}+{\eta }_{a c}^{2}}$.

At high temperature we can use the expansion

$$\begin{eqnarray} &&\ln \left (2\cosh \left (\frac{{E}_{\mathbf{k}}}{2 \theta }\right )\right )\approx \ln 2+\frac{1}{2}\left (\frac{{E}_{\mathbf{k}}}{2 \theta }\right )^{2}-\frac{1}{1 2}\left (\frac{{E}_{\mathbf{k}}}{2 \theta }\right )^{4}\nonumber\\ \end{eqnarray} \tag{ B.1 }$$

which, when inserted in the formula (4), integrates explicitly. For determining the critical temperatures to the first approximation it suffices to restrict ourselves to the second-power terms. This results in an expression quadratic in ζ_τ. Combining the thus obtained 'kinetic' energy

$$\begin{eqnarray} &&-2 \theta \cdot \frac{9}{8{\theta }^{2}}\left (\frac{1}{2}{J}_{a}^{2}{\zeta }_{a}^{2}+\frac{1}{2}{J}_{c}^{2}{\zeta }_{c}^{2}+{J}_{a c}^{2}{\zeta }_{a c}^{2}\right ) \end{eqnarray} \tag{ B.2 }$$

with the potential energy terms from equation (4), we arrive at

$$\begin{eqnarray} &&{F}_{\mathrm{HT}}=\mathop{\sum }\limits _{\tau }\left [-\frac{9}{8 \theta }J_{\tau }^{2}+3{J}_{\tau }\right ]\zeta _{\tau }^{2}. \end{eqnarray} \tag{ B.3 }$$

This result can be improved with use of the Ginzburg–Landau approximate free energy F_GL(ζ_τ,θ), which involves the terms up to the fourth power in ζ_τs. They appear from the integration of the fourth power of the spectrum which is performed explicitly and yields the 'kinetic' energy of the form

$$\begin{eqnarray*} &&\fl \frac{1}{7 6 8{\theta }^{3}}[3 6{\zeta }_{a}^{2}{J}_{a}^{2}(2 7{\zeta }_{a c}^{2}{J}_{a c}^{2}+9{\zeta }_{c}^{2}{J}_{c}^{2}-2 4{\theta }^{2})\nonumber\\ &&+~2 4 3{\zeta }_{a}^{4}{J}_{a}^{4}+3(3 6{\zeta }_{a c}^{2}{J}_{a c}^{2}(9{\zeta }_{c}^{2}{J}_{c}^{2}-1 6{\theta }^{2}))\nonumber\\ &&+~4 8 6{\zeta }_{a c}^{4}{J}_{a c}^{4}-2 8 8{\theta }^{2}{\zeta }_{c}^{2}{J}_{c}^{2}+8 1{\zeta }_{c}^{4}{J}_{c}^{4}], \end{eqnarray*}$$

which together with the potential energy yields the free energy F_GL(ζ_τ,θ) used for further analysis.

In section 2.2 we gave an impression of the complexity of the phase diagram of the RVB model with three exchange parameters. We also gave a brief description of the most characteristic features of the qDoS in various RVB phases. Here we provide a brief sketch of the derivation of the qDoS given in table C.1.

Table C.1.  The qDoS in the possible phases of the c–a–ca-RVB model in the high-temperature approximation. Their characteristic graphs are given in figure 3. The pseudogap/bandwidth parameters in this table are A = 3J_aζ_a, C = 3J_cζ_c, B = 3J_acζ_ac. In the last column K stands for the complete elliptic integral of the first kind. Details of the derivation will be communicated elsewhere.

Character	Color code	g(ε)
Curie	Gray	δ(ε)
Pauli	Red	$\frac{2}{\pi }\frac{1}{\sqrt{{C}^{2}-{\varepsilon }^{2}}}$
Pauli	Green	$\frac{2}{\pi }\frac{1}{\sqrt{{A}^{2}-{\varepsilon }^{2}}}$
$\frac{2{\mu }_{\mathrm{B}}^{2}}{{\pi }^{2}B}\ln \frac{3 2 B}{\pi {\mathrm{e}}^{\gamma }\theta }$	Blue	$\frac{2}{{\pi }^{2}B}\mathsf{K}\left (\frac{\sqrt{4{B}^{2}-{\varepsilon }^{2}}}{2 B}\right )$
Pauli	Magenta	$\frac{4}{{\pi }^{2}}\frac{\mathsf{K}\left (\frac{2 B}{\sqrt{4{B}^{2}+{C}^{2}-{\varepsilon }^{2}}}\right )}{\sqrt{4{B}^{2}+{C}^{2}-{\varepsilon }^{2}}}$
Pauli	Cyan	$\frac{4}{{\pi }^{2}}\frac{\mathsf{K}\left (\frac{2 B}{\sqrt{4{B}^{2}+{A}^{2}-{\varepsilon }^{2}}}\right )}{\sqrt{4{B}^{2}+{A}^{2}-{\varepsilon }^{2}}}$
Arrhenius	Orange	$\left \{ \begin{array}{@{}l@{}} \displaystyle \frac{4 \varepsilon }{{\pi }^{2}A C}\mathsf{K}\left (\frac{\varepsilon \sqrt{{A}^{2}+{C}^{2}-{\varepsilon }^{2}}}{A C}\right ),\tqs \varepsilon \lt \min (A,C) \\ \displaystyle \frac{4}{{\pi }^{2}}\frac{\mathsf{K}\left (\frac{A C}{\varepsilon \sqrt{{A}^{2}+{C}^{2}-{\varepsilon }^{2}}}\right )}{\sqrt{{A}^{2}+{C}^{2}-{\varepsilon }^{2}}},\tqs \min \left (A,C\right )\lt \varepsilon \lt \max \left (A,C\right ) \\ \displaystyle \frac{4 \varepsilon }{{\pi }^{2}A C}\mathsf{K}\left (\frac{\varepsilon \sqrt{{A}^{2}+{C}^{2}-{\varepsilon }^{2}}}{A C}\right ),\tqs \varepsilon \gt \max \left (A,C\right ) \end{array}\right.$
Arrhenius	Yellow	$\left \{ \begin{array}{@{}l@{}} \displaystyle \frac{4}{{\pi }^{2}}\frac{\varepsilon \mathsf{K}\left (\sqrt{\frac{{\varepsilon }^{2}\left ({A}^{2}+4{B}^{2}+{C}^{2}-{\varepsilon }^{2}\right )}{{A}^{2}{C}^{2}+4{B}^{2}{\varepsilon }^{2}}}\right )}{\sqrt{{A}^{2}{C}^{2}+4{B}^{2}{\varepsilon }^{2}}},\tqs \varepsilon \lt \min \left (A,C\right ) \\ \displaystyle \frac{4}{{\pi }^{2}}\frac{\mathsf{K}\left (\sqrt{\frac{{A}^{2}{C}^{2}+4{B}^{2}{\varepsilon }^{2}}{{\varepsilon }^{2}\left ({A}^{2}+4{B}^{2}+{C}^{2}-{\varepsilon }^{2}\right )}}\right )}{\sqrt{{A}^{2}+4{B}^{2}+{C}^{2}-{\varepsilon }^{2}}},\tqs \min (A,C)\lt \varepsilon \lt \max (A,C) \\ \displaystyle \frac{4}{{\pi }^{2}}\frac{\varepsilon \mathsf{K}\left (\sqrt{\frac{{\varepsilon }^{2}\left ({A}^{2}+4{B}^{2}+{C}^{2}-{\varepsilon }^{2}\right )}{{A}^{2}{C}^{2}+4{B}^{2}{\varepsilon }^{2}}}\right )}{\sqrt{{A}^{2}{C}^{2}+4{B}^{2}{\varepsilon }^{2}}},\tqs \varepsilon \gt \max \left (A,C\right ) \end{array}\right.$

The definition of the qDoS reads

$$\begin{eqnarray} &&g(\varepsilon )=\frac{1}{4{\pi }^{2}}\int _{\mathrm{BZ}}\delta \left (\varepsilon -{E}_{\mathbf{k}}\right ){\mathrm{d}}^{2}\mathbf{k}. \end{eqnarray} \tag{ C.1 }$$

Following [20] we insert the integral representation for the Dirac δ-function,

$$\begin{eqnarray*} g(\varepsilon )&=&\frac{1}{2 \pi }\frac{1}{4{\pi }^{2}}\int _{\mathrm{BZ}}\int \nolimits _{-\infty }^{\infty }\mathrm{d}t\exp \left (\mathrm{i}t\left (\varepsilon -{E}_{\mathbf{ k}}\right )\right ){\mathrm{d}}^{2}\mathbf{k}\nonumber\\ &=&\frac{1}{2 \pi }\int \nolimits _{-\infty }^{\infty }\mathrm{d}t\left [\exp \left (\mathrm{i}t \varepsilon \right )\frac{1}{4{\pi }^{2}}\int _{\mathrm{BZ}}\exp \left (-\mathrm{i}t{E}_{\mathbf{k}}\right ){\mathrm{d}}^{2}\mathbf{k}\right ]. \end{eqnarray*}$$

For all phases having lines of nodes the integration over one of the components of the wavevector k in the BZ is performed and yields an intermediate result in terms of the Bessel and Struve functions of arguments dependent on the Fourier transformation variable t and the remaining component of the wavevector k. The Fourier transforms with respect to t can be performed for the intermediate answers of that form. The integration yields integrands of the elliptic integrals over the remaining component of the wavevector k. This solves the problem of calculating the qDoS for the dispersion laws with the lines of nodes in the BZ. The results are given in the respective cells of table C.1.

In order to cope with the remaining two phases whose spectra contain only nodal points one more trick solves the problem, namely performing the previous moves for the squared spectrum and thus obtaining the distribution of the states as a function of their squared energy. The distribution of squares of the quasiparticle energies is

$$\begin{eqnarray} &&\varrho ({\varepsilon }^{2})=\frac{1}{4{\pi }^{2}}\int _{\mathrm{BZ}}\delta \left ({\varepsilon }^{2}-{E}_{\mathbf{ k}}^{2}\right ){\mathrm{d}}^{2}\mathbf{k}. \end{eqnarray} \tag{ C.2 }$$

Then the sought qDoS is given by [21]

$$\begin{eqnarray} &&g\left (\varepsilon \right )=2 \varepsilon \varrho ({\varepsilon }^{2}). \end{eqnarray} \tag{ C.3 }$$

To obtain ϱ(ε²) we again use the integral representation of the Dirac δ-function,

$$\begin{eqnarray*} \varrho ({\varepsilon }^{2})&=&\frac{1}{2 \pi }\frac{1}{4{\pi }^{2}}\int _{\mathrm{BZ}}\int \nolimits _{-\infty }^{\infty }\mathrm{d}t\exp \left (\mathrm{i}t\left ({\varepsilon }^{2}-{E}_{\mathbf{ k}}^{2}\right )\right ){\mathrm{d}}^{2}\mathbf{k}\nonumber\\ &=&\frac{1}{2 \pi }\int \nolimits _{-\infty }^{\infty }\mathrm{d}t\left [\exp (\mathrm{i}t{\varepsilon }^{2})\frac{1}{4{\pi }^{2}}\int _{\mathrm{BZ}}\exp \left (-\mathrm{i}t{E}_{\mathbf{k}}^{2}\right ){\mathrm{d}}^{2}\mathbf{k}\right ]. \end{eqnarray*}$$

Remarkably enough, sequential integrations of the squared spectrum over one of the components of the wavevector k and Fourier transformation with respect to t yield expressions of the same form as integration of the spectrum itself: the Bessel function and an integrand of the elliptic integral. Thus, the final integration over the remaining component of the wavevector k yields some elliptic integrals given in the respective cells of table C.1.

For the semiquantitative treatment we use the standard definition of the magnetic susceptibility per spin,

$$\begin{eqnarray} &&\chi =-2{\mu }_{\mathrm{B}}^{2}\int g(\varepsilon )\frac{\partial f\left (\frac{\varepsilon }{\theta }\right )}{\partial \varepsilon }\mathrm{d}\varepsilon , \end{eqnarray} \tag{ D.1 }$$

where f(...) is the Fermi distribution function. For a pseudogap to open in a state with, say, nonvanishing ζ_c (red or magenta phase) one needs that at an observable critical temperature the OP ζ_a splits from zero (orange or yellow phase) at some observable temperature. We model the temperature dependence of the pseudogap

$$\begin{eqnarray} &&3{J}_{a}{\zeta }_{a}\left (\theta \rightarrow 0\right )\left (1-\frac{\theta }{{\theta }^{\mathrm{crit}}}\right )^{\beta } \end{eqnarray} \tag{ D.2 }$$

so that it appears below the critical temperature θ^crit (where by θ^crit one of the temperatures ${\theta }_{c\rightarrow c,a}^{\mathrm{crit}}$ or ${\theta }_{c,a c\rightarrow c,a,a c}^{\mathrm{crit}}$ is meant) and perform the numerical integration of the expression (D.1) for the susceptibility with the qDoS for the RVB phase with two pseudogaps where we also set ζ_ac = 0.⁵

The values of θ^crit (${\theta }_{c\rightarrow c,a}^{\mathrm{crit}}$ or ${\theta }_{c,a c\rightarrow c,a,a c}^{\mathrm{crit}}$), 3J_aζ_a(θ → 0), 3J_cζ_c(θ → 0), and β can be adjusted to reproduce the experimental course of the susceptibility [4] shown in figure D.1. The scale of 3J_cζ_c(θ → 0) is given by the magnitude of the Pauli paramagnetic susceptibility at a higher temperature and yields the parameters of the original model (1). The fact that the classical exponent $\frac{1}{2}$ coming from our simple high-temperature treatment allows for an acceptable fit of the experimental susceptibility is pretty remarkable. However, one has to realize that this is not the true critical exponent (for which as we noticed above we do not have enough data), but rather a value permitting us to interpolate the susceptibility and thus the OP in a wider range of temperatures. For the sake of simplicity we assume that the phase with the finite qDoS at zero energy is the 1D-RVB one with nonvanishing ζ_c. Then, using the zero temperature limiting value of this OP, ${\zeta }_{c}\left (\theta \rightarrow 0\right )=\frac{1}{\pi }$ [3], we get J_c = 1310 K. This nicely agrees with the original estimate of [1] (about 1000 K). It is also what one can intuitively expect relying on the Goodenough–Kanamori rules [7]⁶. Using it in the high-temperature estimate for the critical temperature (5) we obtain for the highest critical temperature (that of the transition from the Curie paramagnetic to the 1D-RVB state) the value ${\theta }_{c}^{\mathrm{crit}}=4 9 0~\mathrm{K}$, which lies a fair amount above the decomposition temperature [1]. In the case where the Q1D-RVB phase is assumed to be responsible for the Pauli paramagnetism we take J_ac = 1400 K as a plausible estimate. By doing so the critical temperatures ${\theta }_{a c}^{\mathrm{crit}}=5 3 5~\mathrm{K}$ and ${\theta }_{c,a c}^{\mathrm{crit}}=4 3 5~\mathrm{K}$ derive to be considerably higher than the decomposition temperatures. In that sense these estimates are consistent (neither of phases with divergent paramagnetic susceptibility is observed). For the exchange parameter J_a we notice that its values are pretty stable. Considering the fitted critical value θ^crit as the critical temperature ${\theta }_{c,a c\rightarrow c,a,a c}^{\mathrm{crit}}$ given by (7) yields J_a = 560 K. Combining this with the zero-temperature limit of the pseudogap, 3J_aζ_a(θ → 0) = 140 K, we get ζ_a(θ → 0) = 0.08.

Zoom In Zoom Out Reset image size

The semiquantitative agreement of the data shown in figure 6(a) is achieved by confronting the slope of a versus θ in the yellow phase, −7 × 10⁻⁶ Å K⁻¹ (green line), with the model (12), which yields for the slope

$$\begin{eqnarray} &&\left .\left (\frac{\mathrm{d}}{\mathrm{d}\theta }\delta a\right )\right \vert _{\theta \rightarrow {\theta }^{\mathrm{crit}}}=\frac{3{J}_{a,a}^{\prime}}{{K}_{a a}}\left [{\zeta }_{c}\left (\theta \rightarrow 0\right )\right ]^{2}\left (-\frac{1}{{\theta }^{\mathrm{crit}}}\right ) \end{eqnarray} \tag{ D.3 }$$

(θ^crit = 100 K). Then the ratio of parameters describing the spin–phonon coupling, ${J}_{a,a}^{\prime}/{K}_{a a}=0.0 3 4~\mathring{A}$, seems to fairly fit the expectations.

The observed geometry of the crystal with RVB phases corresponds to the minimum with respect to ρ_λ of the free energy (4) to which the elastic energy necessary to deform the hypothetical structure with no exchange interaction

$$\begin{eqnarray} &&\frac{1}{2}\mathop{\sum }\limits _{\mu \lambda }{K}_{\mu \lambda }{\rho }_{\mu }{\rho }_{\lambda } \end{eqnarray} \tag{ E.1 }$$

for each nearest neighbor Cu–Cu pair is added and the exchange parameters are replaced according to the recipe (10). This yields the mechanic equilibrium condition of the crystal in an RVB phase,

$$\begin{eqnarray} &&\frac{\partial F}{\partial {\rho }_{\lambda }}+\mathop{\sum }\limits _{\mu \lambda }{K}_{\mu \lambda }{\rho }_{\mu }=0. \end{eqnarray} \tag{ E.2 }$$

We notice that the free energy F given by expression (4) has a special form. The first term ('kinetic' energy) is an integral of a function of the dispersion law E_k whose argument has the form

$$\begin{eqnarray} &&\mathop{\sum }\limits _{\tau }A_{\tau }^{2}{J}_{\tau }^{2}{\zeta }_{\tau }^{2}{f}_{\tau }^{2}\left (\mathbf{k}\right ), \end{eqnarray} \tag{ E.3 }$$

where A_τ are numerical coefficients (3 for τ = a,c; 6 for τ = ac); f_τ(k) are trigonometrical expressions (cosk_τ for τ = a,c; cosk_acosk_c for τ = ac). The 'potential' energy contribution to the free energy (4) is

$$\begin{eqnarray} &&\mathop{\sum }\limits _{\tau }{A}_{\tau }{J}_{\tau }\zeta _{\tau }^{2}. \end{eqnarray} \tag{ E.4 }$$

Due to the above special form of the 'kinetic' and 'potential' energies, the self-consistency equations for ζ_τ each have the form

$$\begin{eqnarray*} &&\int _{\mathrm{BZ}}\tanh \left (\frac{{E}_{\mathbf{k}}}{2 \theta }\right ){E}_{\mathbf{k}}^{\prime}{A}_{\tau }^{2}{J}_{\tau }^{2}{\zeta }_{\tau }{f}_{\tau }^{2}\left (\mathbf{k}\right ){\mathrm{d}}^{2}\mathbf{k}={A}_{\tau }{J}_{\tau }{\zeta }_{\tau }, \end{eqnarray*}$$

to be solved simultaneously for all τ. Apparently, whatever (sub)set of ζ_τ = 0 satisfies the equations, for the nonvanishing OPs the self-consistency equations acquire the form

$$\begin{eqnarray} &&\int _{\mathrm{BZ}}\tanh \left (\frac{{E}_{\mathbf{k}}}{2 \theta }\right ){E}_{\mathbf{k}}^{\prime}{f}_{\tau }^{2}\left (\mathbf{k}\right ){\mathrm{d}}^{2}\mathbf{k}=\frac{1}{{A}_{\tau }{J}_{\tau }}. \end{eqnarray} \tag{ E.5 }$$

Now we can turn to the structure determination. The derivative of the kinetic energy with respect to the geometry parameters ρ_λ reads

$$\begin{eqnarray} &&-2\mathop{\sum }\limits _{\tau }A_{\tau }^{2}\left ({J}_{\tau }{\zeta }_{\tau }^{2}{J}_{\tau ,\lambda }^{\prime}+{J}_{\tau }^{2}{\zeta }_{\tau }\frac{\partial {\zeta }_{\tau }}{\partial {\rho }_{\lambda }}\right )\nonumber\\ &&\times ~\int _{\mathrm{BZ}}\tanh \left (\frac{{E}_{\mathbf{k}}}{2 \theta }\right ){E}_{\mathbf{k}}^{\prime}{f}_{\tau }^{2}\left (\mathbf{k}\right ){\mathrm{d}}^{2}\mathbf{k}. \end{eqnarray} \tag{ E.6 }$$

The value of the above integral for the equilibrium values of the OPs is given by the self-consistency conditions (E.5). Thus, the derivative of the kinetic energy is rewritten as

$$\begin{eqnarray} &&-2\mathop{\sum }\limits _{\tau }{A}_{\tau }\left (\zeta _{\tau }^{2}{J}_{\tau ,\lambda }^{\prime}+{J}_{\tau }{\zeta }_{\tau }\frac{\partial {\zeta }_{\tau }}{\partial {\rho }_{\lambda }}\right ). \end{eqnarray} \tag{ E.7 }$$

Combining this with the derivative of the potential energy and that of the elastic energy, we get

$$\begin{eqnarray} &&\fl -2\mathop{\sum }\limits _{\tau }{A}_{\tau }\left (\zeta _{\tau }^{2}{J}_{\tau ,\lambda }^{\prime}+{J}_{\tau }{\zeta }_{\tau }\frac{\partial {\zeta }_{\tau }}{\partial {\rho }_{\lambda }}\right )\nonumber\\ &&+~\mathop{\sum }\limits _{\tau }{A}_{\tau }\left (\zeta _{\tau }^{2}J_{\tau ,\lambda }^{\prime}+2{J}_{\tau }{\zeta }_{\tau }\frac{\partial {\zeta }_{\tau }}{\partial {\rho }_{\lambda }}\right )+\mathop{\sum }\limits _{\mu \lambda }{K}_{\mu \lambda }{\rho }_{\mu }=0.\nonumber\\ \end{eqnarray} \tag{ E.8 }$$

The terms including the derivative $\frac{\partial {\zeta }_{\tau }}{\partial {\rho }_{\lambda }}$ stemming from the kinetic and potential energies cancel each other, a remarkable consequence of the Hellmann–Feynman and virial theorems, which immediately results in

$$\begin{eqnarray} &&-\mathop{\sum }\limits _{\tau }{A}_{\tau }\left (\zeta _{\tau }^{2}{J}_{\tau ,\lambda }^{\prime}\right )+\mathop{\sum }\limits _{\mu \lambda }{K}_{\mu \lambda }{\rho }_{\mu }=0, \end{eqnarray} \tag{ E.9 }$$

which after some trivial algebra results in

$$\begin{eqnarray} &&{\rho }_{\mu }=\mathop{\sum }\limits _{\lambda }(K)_{\mu \lambda }^{-1}\left (\mathop{\sum }\limits _{\tau }{A}_{\tau }\zeta _{\tau }^{2}{J}_{\tau ,\lambda }^{\prime}\right ). \end{eqnarray} \tag{ E.10 }$$

For the lattice parameters a and c, we introduce the geometry variables ρ_a = δa = (a − a₀), 2ρ_c = δc = (c − c₀) and apply the general relation (11). Then the variation of the two lattice parameters is

$$\begin{eqnarray} &&\left (\begin{array}{@{}c@{}} \displaystyle \delta a\\ \displaystyle \delta c \end{array}\right )&=&\frac{3}{\vert K\vert }\left [\left (\begin{array}{@{}c@{}} \displaystyle \left ({K}_{c c}{J}_{a,a}^{\prime}-2{K}_{a c}{J}_{a,c}^{\prime}\right )\\ \displaystyle \left (-{K}_{a c}{J}_{a,a}^{\prime}+2{K}_{a a}{J}_{a,c}^{\prime}\right ) \end{array}\right ){\zeta }_{a}^{2}\right .\nonumber\\ &&\left .+~\left (\begin{array}{@{}c@{}} \displaystyle \left ({K}_{c c}{J}_{c,a}^{\prime}-2{K}_{a c}{J}_{c,c}^{\prime}\right )\\ \displaystyle \left (-{K}_{a c}{J}_{c,a}^{\prime}+2{K}_{a a}{J}_{c,c}^{\prime}\right ) \end{array}\right ){\zeta }_{c}^{2}\right ], \end{eqnarray} \tag{ E.11 }$$

where we denote by |K| the determinant of the 2 × 2 matrix of the force constants and make use of our previous consideration allowing us to neglect the OP ζ_ac. It seems, however, to be reasonable that the exchange constant J_a is independent of the deformation in the c-direction and, thus, ${J}_{a,c}^{\prime}=0$. We may therefore write

$$\begin{eqnarray} &&\left (\begin{array}{@{}c@{}} \displaystyle \delta a\\ \displaystyle \delta c \end{array}\right )&=&\frac{3}{\vert K\vert }\left [{J}_{a,a}^{\prime}{\zeta }_{a}^{2}\left (\begin{array}{@{}c@{}} \displaystyle {K}_{c c}\\ \displaystyle -{K}_{a c} \end{array}\right )\right .\nonumber\\ &&\left .+~{\zeta }_{c}^{2}\left (\begin{array}{@{}c@{}} \displaystyle \left ({K}_{c c}{J}_{c,a}^{\prime}-2{K}_{a c}{J}_{c,c}^{\prime}\right )\\ \displaystyle \left (-{K}_{a c}{J}_{c,a}^{\prime}+2{K}_{a a}{J}_{c,c}^{\prime}\right ) \end{array}\right )\right ]. \end{eqnarray} \tag{ E.12 }$$

Turning to the temperature dependence of the lattice parameters as given by their slopes versus temperature at the respective critical points, we notice that assuming a power law for the OPs like in the interpolation formula (D.2),

$$\begin{eqnarray} &&{\zeta }_{\tau }\left (\theta \right )={\zeta }_{\tau }^{\ast }\left (\theta \right )\left (1-\frac{\theta }{{\theta }_{\tau }^{\mathrm{crit}}}\right )^{\beta }, \end{eqnarray} \tag{ E.13 }$$

where the ${\zeta }_{\tau }^{\ast }\left (\theta \right )$s are smooth in ${\theta }_{\tau }^{\mathrm{crit}}$ and ${\theta }_{\tau }^{\mathrm{crit}}$ are the reference temperatures given in table 1 for the corresponding OPs in the relevant phases, immediately yields

$$\begin{eqnarray} \left .\left (\frac{\mathrm{d}}{\mathrm{d}\theta }\left [{\zeta }_{\tau }\left (\theta \right )\right ]^{2}\right )\right \vert _{\theta \rightarrow {\theta }^{\mathrm{crit}}}&=&-2 \beta \left (1-\frac{\theta }{{\theta }_{\tau }^{\mathrm{crit}}}\right )^{2 \beta -1}\nonumber\\ &&\times ~\left (\frac{1}{{\theta }_{\tau }^{\mathrm{crit}}}\right )\left [{\zeta }_{\tau }^{\ast }\left ({\theta }_{\tau }^{\mathrm{crit}}\right )\right ]^{2},\nonumber\\ \end{eqnarray} \tag{ E.14 }$$

which for the interpolated value of the exponent β = 0.5 provides the results used in our analysis.