Tomonaga-Luttinger liquid spin dynamics in the quasi-one dimensional Ising-like antiferromagnet BaCo$_2$V$_2$O$_8$

TL;DR

This study combines neutron scattering and simulations to explore the spin dynamics of BaCo$_2$V$_2$O$_8$, revealing a quantum phase transition and Luttinger liquid behavior in a quasi-one-dimensional Ising antiferromagnet.

Contribution

It provides experimental evidence of Tomonaga-Luttinger liquid behavior in a real material undergoing a quantum phase transition under magnetic field.

Findings

Observation of a quantum phase transition at 3.8 T

Reconfiguration of spin dynamics near the critical field

Emergence of longitudinal excitations characteristic of Luttinger liquids

Abstract

Combining inelastic neutron scattering and numerical simulations, we study the quasi-one dimensional Ising anisotropic quantum antiferromagnet \bacovo\ in a longitudinal magnetic field. This material shows a quantum phase transition from a N\'eel ordered phase at zero field to a longitudinal incommensurate spin density wave at a critical magnetic field of 3.8 T. Concomitantly the excitation gap almost closes and a fundamental reconfiguration of the spin dynamics occurs. These experimental results are well described by the universal Tomonaga-Luttinger liquid theory developed for interacting spinless fermions in one dimension. We especially observe the rise of mainly longitudinal excitations, a hallmark of the unconventional low-field regime in Ising-like quantum antiferromagnet chains.

Figures

Equations

$${\mathcal{H}}_{\mathrm{intra},\mu }=$$

$$-g_{zz}{\mu }_B{\mu }_{0}H\sum _n{S}_{n,\mu }^z,$$

$${\mathcal{H}}_{\mathrm{eff}}=J\sum _n(S_n^x{S}_{n+1}^x+S_n^y{S}_{n+1}^y+\Delta S_n^z{S}_{n+1}^z)-g_{zz}{\mu }_B{\mu }_{0}H\sum _n{S}_n^z+J^{\prime }\sum _n⟨S_n^z⟩S_n^z.$$

$$S_{n}(\bm{Q},\omega)\equiv\frac{d^{2}\sigma}{d\Omega dE}\propto\frac{|\bm{q}^{\prime}|}{|\bm{q}|}\sum_{\alpha,\beta=x,y,z}\Big{(}\delta_{\alpha\beta}-\frac{Q_{\alpha}Q_{\beta}}{|\bm{Q}|^{2}}\Big{)}|F(\textbf{Q})|^{2}S_{\alpha\beta}(\bm{Q},\omega),$$

$$S_{\alpha\beta}(\bm{Q},\omega)=\bigg{|}\mathrm{Im}\int dt\sum_{\bm{r}}e^{i(\omega t-\bm{Q}\cdot\bm{r})}C_{\rm ret}^{\alpha\beta}(\bm{r},t)\bigg{|}.$$

$$C_{\mathrm{ret}}^{\alpha \beta }(\mathbf{r},t)=-i{\vartheta }_{\mathrm{step}}(t)⟨[S^{\alpha }(\mathbf{r},t),S^{\beta }(0,0)]⟩,$$

$$S_{n}(\bm{Q},\omega)\propto\frac{|\bm{q}^{\prime}|}{|\bm{q}|}|F(\textbf{Q})|^{2}\Big{[}\Big{(}1+\frac{Q_{z}^{2}}{|\bm{Q}|^{2}}\Big{)}S_{xx}(\bm{Q},\omega)+\Big{(}1-\frac{Q_{z}^{2}}{|\bm{Q}|^{2}}\Big{)}S_{zz}(\bm{Q},\omega)\Big{]}.$$

$${\mathcal{H}}_{\mathrm{inter}}=J^{\prime }\sum _n\sum _{⟨\mu ,\nu ⟩}{\mathbf{S}}_{n,\mu }\cdot {\mathbf{S}}_{n,\nu }\simeq J^{\prime }\sum _{n,\mu }{\mathbf{S}}_{n,\mu }\cdot ⟨{\mathbf{S}}_n⟩.$$

Full text

Quantum magnets offer an extremely rich variety of phases ranging from the conventional long-range ordered ones, dubbed spin “solids”, to various kinds of spin “liquids”. In the latter, the excitations have often an unconventional nature such as a topological character or fractional quantum numbers. Among such systems, one dimensional (1D) quantum magnets are especially interesting in that the topological excitations are the norm rather than the exception, and because the interplay between exchange coupling and extremely strong quantum fluctuations due to the reduced dimensionality gives rise to profuse physical phenomena giamarchi2004 .

On the experimental front, the recent realization of quantum magnets with relatively weak magnetic exchange has paved a new avenue to an efficient manipulation of systems with realizable magnetic fields, enabling novel phases and phenomena to be probed experimentally. Plentiful examples of such successful investigations exist, e.g. scaling properties of Bose-Einstein condensation giamarchi-ruegg-magnonBEC2008 ; batista-RMP-magnonBEC2014 , quantitative tests of Tomonaga-Luttinger liquid (TLL) theory klanjsek-PRL-LLinladder2008 ; bouillot-PRB-ladder2011 ; schmidiger-PRL-DIMPY2013 , scaling properties at quantum critical points zheludev-PRB-CriticalScale2017 ; zheludev-PRL-CriticalBPCB2018 , fractionalized excitations broholm-takagi-PRL-SrCuO2-spinon2004 ; thielemann-PRL-fractional-ladder2009 , topological phase transitions faure2018 , other exotic excitations zheludev-giamarchi-tsvelik-PRB-ladder2013 ; grenier2015 ; bera2017 ; wang2018 . The effect of an external magnetic field competing with the excitation gap associated to rung-singlets klanjsek-PRL-LLinladder2008 or to the Haldane state renard-PRL-MagneticFieldHaldaneChain1989 for instance is especially interesting. Quite remarkably, all these transitions fall into the same universality class, the so called Pokrovsky-Talapov commensurate-incommensurate (C-IC) phase transition Talapov1979 ; giamarchi2004 , which is also pertinent to the Mott transition in itinerant electronic systems. Hence there is a considerable interest in experimental analyses of such phenomena, and investigations have been conducted in systems such as bosons in a periodic lattice naegerl-Nature-pinning-transition2010 ; modugno-PRA-BosonMott2016 , spin-1 chains zvyagin-PRL-DTN-magfield2007 and spin-1/2 ladders klanjsek-PRL-LLinladder2008 ; bouillot-PRB-ladder2011 . However, in these realizations, magnetic excitations in the IC phase are dominated by spin-spin correlations transverse to the applied field, and a study of the opposite and more exotic case, where the longitudinal excitations are dominant, is still lacking.

In this paper, we focus on this particular case. We investigate the Ising-like compound BaCo2V2O8 under a magnetic field along the anisotropy axis by combining inelastic neutron scattering experiments and numerical simulations. We show that the quantum phase transition provoked by a longitudinal field of 3.8 T is indeed in the C-IC universality class through the analysis of spin-spin dynamical correlations. Furthermore, we demonstrate that most of the spectral weight in the IC phase consists in longitudinal excitations, which are a strong fingerprint of TLL dynamics with IC solitonic excitations.

BaCo2V2O8 consists of screw chains of Co2+ ions running along the fourfold $c$-axis of a body-centered tetragonal structure [Fig. 1(b)] wichmann1986 . Due to an anisotropic $g$ tensor kimura2006 , the Co2+ magnetic moments are described effectively by weakly coupled spin-1/2 $XXZ$ (Ising-like) chains abragam1951 . The Hamiltonian includes intrachain and interchain interactions ${\cal H}=\sum_{\mu}{\cal H}_{{\rm intra},\mu}+{\cal H}_{\rm inter}$, which write

[TABLE]

and ${\cal H}_{\rm inter}=J^{\prime}\sum_{n}\sum_{\langle\mu,\nu\rangle}{\bf S}_{n,\mu}\cdot{\bf S}_{n,\nu}$. Here ${\bf S}_{n,\mu}$ is a spin-1/2 operator, $n$ the site index along the chain, $\mu,\nu$ label different chains, $J(>0)$ is the antiferromagnetic (AF) intrachain interaction, and $\Delta$ the Ising anisotropy. $g_{zz}\mu_{B}\mu_{0}H\sum_{n}S^{z}_{n,\mu}$ is the Zeeman term from the longitudinal field, with $g_{zz}$ the Landé factor and $\mu_{B}$ the Bohr magneton. The $a,b,c$ crystallographic axes coincide with the spin $x,y,z$ axes, respectively. The interchain coupling is treated by mean field theory supmat1 . At $H=0$ and $T\leq T_{N}$ ($T_{N}=5.4$ K), BaCo2V2O8 is in a gapped AF phase and the magnetic moments point along the Ising $c$-axis [Fig. 1(b)]. The elementary excitations are spinons, which are confined by the interchain coupling to form spinon bound states. They give rise to two series of discretized energy levels dispersing along the $c$-axis (and only weakly in the perpendicular directions), which have longitudinal ($\Delta S^{z}=0$) and transverse ($\Delta S^{z}=\pm 1$) character with respect to the anisotropy axis grenier2015 . The ground state phase diagram of a single spin-1/2 $XXZ$ chain under the application of a longitudinal magnetic field is shown in Fig. 1(a). In the Ising-like case ($\Delta>1$), $H>H_{c}$ is required to enter the TLL phase and close the excitation gap. The TLL phase is characterized by spatial spin-spin correlations transverse $C^{xx}(r)\equiv\langle S^{x}_{r}S^{x}_{0}\rangle\propto r^{-1/(2K)}$ and longitudinal $C^{zz}(r)\equiv\langle S^{z}_{r}S^{z}_{0}\rangle-m_{z}^{2}\propto r^{-2K}$ to the field direction, where $m_{z}$ is the field-induced uniform magnetization per site. The decay of $C^{zz}(r)$ and $C^{xx}(r)$ are dictated by the TLL parameter $K$. The field dependence of $K$ causes a crossover at $H^{*}$ from a low-field regime [red-shaded area in Fig. 1(a)] where the physics is dominated by $C^{zz}(r)$ to a high-field regime [grey-shaded area in Fig. 1(a)] dominated by $C^{xx}(r)$. The dispersion of low-energy excitations is expected to become gapless at both C and IC wave vectors $q=\pi,2\pi m_{z}$ for transverse excitations (captured by the space-time correlation $\langle S^{x}_{r}(t)S^{x}_{0}(0)\rangle$), and at $q=0,\pi(1-2m_{z})$ for longitudinal excitations (captured by $\langle S^{z}_{r}(t)S^{z}_{0}(0)\rangle$) muller1981 ; chitra1995 ; chitra1997 ; fath2003 . For BaCo2V2O8 ($\Delta=1.9$), the quantum phase transition occurs at $\mu_{0}H_{c}=3.8$ T from the Néel phase to the longitudinal spin density wave (LSDW) with an IC wave vector, both ordered phases stabilized by weak interchain couplings. In the latter phase, the magnetic moments are parallel to the field (and Ising) direction while their amplitude is spatially modulated [Fig. 1(c)]  kimura2008a ; kimura2008b ; canevet2013 ; supmat1 . When the external field is further increased, the LSDW phase is replaced by a canted AF order with staggered moments perpendicular to the $c$-axis above $\mu_{0}H^{*}\approx 9$ T grenierPRB2015 ; klanjsek2015 , which corresponds to the crossover from the TLL longitudinal to transverse-dominant correlations, before the magnetization saturates at higher field ($H_{sat}$).

To probe the transition from the Néel to LSDW phase and their spin dynamics, we performed inelastic neutron scattering experiments at the cold-neutron triple axis spectrometer TASP (PSI, Switzerland). We used a horizontal cryomagnet, applying magnetic fields up to 6.8 T. Two BaCo2V2O8 single crystals, grown by floating zone, were co-aligned with an accuracy better than $1^{\circ}$. The magnetic field was applied along the $c^{*}-$axis of the $(a^{*},c^{*})$ scattering plane, hence along the magnetic moment direction. The data were measured at the base temperature of 150 mK with various fixed final wave vectors ranging from 1.06 to 1.3 Å*-1* (yielding an energy resolution from 70 to 150 $\mu$eV). In BaCo2V2O8, the crystallographic zone centers (ZC) are at ${\bf Q}=(h,k,l)$ positions with $h+k+l={\rm even}$. The magnetic Bragg peaks of the Néel phase appear at the AF points ${\bf Q}=(h+1,k,l)$ corresponding to the ${\bf k}_{AF}=(1,0,0)$ propagation vector canevet2013 . The presence of four screw-chains per unit cell folds the excitation branches and replication from the ZC positions is added to the usual contribution from AF points.

Energy scans with constant $Q$ have first been recorded for different magnetic fields at the AF position ${\bf Q}=(2,0,1)$. At $H=0$, the measured lowest energy peak corresponds to the doubly degenerate transverse excitation grenier2015 . The field produces a Zeeman splitting that lifts this degeneracy kimura2007 ; faure2018 , leading to the linear decrease of the lowest transverse mode up to the transition at $H_{c}$, as observed in Fig. 2 (red open circles). The same feature is seen at ZC wave vectors (red closed circles) due to the folding.

In the LSDW phase, the propagation vector becomes ${\bf k}_{LSDW}=(1,0,\delta)$. The field dependence of the IC modulation $\delta$ has been determined from $Q_{L}$-scans. In agreement with the TLL theory and previous report canevet2013 , we have found that it increases with the field as $\delta=2\pi m_{z}$, i.e., the period for the spatial modulation of the magnetic moments becomes shorter  supmat1 . The transition at $H_{c}$ into the LSDW phase also manifests as a change of magnetic excitation spectrum [from circles to triangles in Fig. 2]. To obtain the overall behavior of spin dynamics in this LSDW phase, constant-$Q_{L}$ energy scans have been collected along the $c^{*}$ direction across the AF point ${\bf Q}=(2,0,1)$ at $\mu_{0}H=4.2$ and $6$ T. Figures 3(a) and 3(c) show the corresponding maps as a function of energy transfer and $Q_{L}$. At 4.2 T, a strong excitation is observed, forming an arch bridging the IC positions $(2,0,1\pm\delta)$ over the AF center $(2,0,1)$. The dispersion has minima at the IC positions of the LSDW phase, which is a key signature of this field-induced TLL phase. The data show that the arch-like dispersion expands from 4.2 T to 6 T, while $\delta$ becomes about twice larger: The energy minimum at $(2,0,1\pm\delta)$ remains equal to $\approx$0.1 meV while the energy at the AF point increases.

From the above-mentioned folding, replications are observed around the ZC position ${\bf Q}=(3,0,1)$. This is illustrated in Figs. 4(a)-4(d) showing individual constant-$Q_{L}$ energy scans through ${\bf Q}=(3,0,Q_{L})$ positions with $Q_{L}$ ranging between $1$ and $1.16$ and energy up to 5 meV at 4.2 T. The map gathering such scans is displayed in Fig. 4(e) with a zoom in Fig. 4(f). These results show that most of the intensity is concentrated in an arch-like excitation with minimum energy of the dispersion $\simeq 0.65$ meV at $(3,0,1+\delta)$, the satellite position of the LSDW phase. At 6 T, the intense arch feature expands similarly to the result around the AF position supmat1 . Weaker excitations are also visible around 0.4, 0.8, and 1.5 meV. Further away from the ZC position along $Q_{L}$, only a broad feature remains, possibly corresponding to a continuum of excitations [see Fig. 4(d) for $Q_{L}=1.16\simeq 1+2\delta$]. Although the excitations in the AF and ZC regions show strong similarities, the energy gap at the IC wave vector is significantly smaller at the AF satellite than at the ZC one. This is also visible in Fig. 2 displaying the energy of the intense modes at the two IC positions $(2,0,1+\delta)$ and $(3,0,1+\delta)$. This is ascribed to the finite dispersion perpendicular to the chain direction caused by the interchain coupling, and also observed in zero field grenier2015 .

Aiming at a deeper understanding of the spin dynamics in the LSDW phase, we performed numerical simulations of the $XXZ$ model with a longitudinal magnetic field [Eq. (1)]. We obtained the ground state of the system by density matrix renormalization group white1992 and calculated the retarded correlation function by time-evolving block decimation vidal2003 . The inelastic neutron scattering cross section $S_{n}$ was derived as the Fourier transform of this correlation function faure2018 ; takayoshi2018 . The calculations were performed by considering the full magnetic structure factor of BaCo2V2O8 with the values $J=3.05$ meV and $\Delta=1.9$ J obtained from our previous investigation faure2018 . The agreement is best for interchain coupling $J^{\prime}=0$ and deteriorates with increasing it, especially near the C-IC transition point supmat1 , in contract with our previous estimation of $J^{\prime}=0.17$ meV grenier2015 . This may be due to a mean-field overestimation of its effect particularly in the critical region or to its possible dependence on the longitudinal field since it is an effective coupling derived from a complex set of interactions klanjsek2015 . All the numerical calculations presented here were therefore performed with $J^{\prime}=0$.

The calculated field dependence of $\delta(H)$ globally agrees with the experiment except near the transition supmat1 . We present the comparison of the measured vs calculated excitation spectra in Figs. 3(a) vs 3(b) and 3(c) vs 3(d) around the AF position at 4.2 and 6 T respectively, as well as in Figs. 4(e)-4(f) vs 4(g)-4(h) around the ZC position at 4.2 T. Note that the calculated peaks are broadened (0.3 meV resolution) compared to the experimental ones due to the finite time effect, i.e. the limitations of the calculations within the finite real time domain $0\leq t\leq T$. The main features, i.e. the dispersion of the low energy excitation bridging the two neighboring IC wave vectors and its spectral weight, are well reproduced. The relative intensity of the weaker branches at 0.8 and 1.5 meV at $(3,0,1+\delta)$ is less accurately reproduced, maybe due to the omission of interchain interaction in the calculations.

The nature of the fluctuations can be further analyzed by the numerically calculated transverse and longitudinal parts of the dynamical structure factor, $S_{xx}$ and $S_{zz}$, which are shown in Figs. 3(e)-3(f) around $(2,0,Q_{L})$ at 6 T and in Figs. 4(i)-4(j) around $(3,0,Q_{L})$ at 4.2 T. The most striking result is that the arch-like excitation has longitudinal character around both AF and ZC positions. For $(3,0,Q_{L})$ at 4.2 T, the weaker transverse excitations $S_{xx}$ give rise to two branches going softer toward the C position with minimum energies close to zero and 1.5 meV [Fig. 4(i)], both of which are seen in the experimental data of Fig. 4(f). This result proves that the spin dynamics is dominated by longitudinal fluctuations strongly excited at the IC wave vectors near the AF positions, which replicate around the ZC ones.

A recent THz spectroscopy investigation of the spin dynamics was performed under a longitudinal magnetic field in the gapless regime of SrCo2V2O8, the sister compound of BaCo2V2O8 yang2017 ; wang2018 . In this experiment, only transverse excitations ($S_{xx}$) at C positions could be probed, such as string and (anti)psinon-(anti)psinon, dressing the field-polarized ground state of 1D quantum antiferromagnets described by the Bethe Ansatz karbach2002 ; kohno2009 . Our neutron spectroscopy study opens up new avenues. We could first follow the dispersion in reciprocal space of the psinon-psinon and 2-string excitations corresponding to the weak transverse modes visible near zero and at 1.5 meV for $Q_{L}=1$ in Fig. 4(f). Moreover, both transverse and longitudinal fluctuations could be probed and we have proven that most of the intensity actually comes from longitudinal excitations missed by THz spectroscopy. This finding is essential to understand a growing number of experiments performed on similar systems with other probes. Our results finally pave the way to further investigations of unexplored regimes of the TLL physics in spin systems, such as the crossover from longitudinal to transverse dominant spin-spin correlations at higher magnetic field or the influence of interchain interactions.

In summary, our combined neutron scattering and numerical investigations of the LSDW phase in BaCo2V2O8 show that the quantum phase transition from the Néel to LSDW phase is described by the $XXZ$ model. Clear Tomonaga-Luttinger liquid signatures are observed such as the field-dependent incommensurability of the low energy excitations and the arch-like dispersion. The most striking result concerns the longitudinal nature of the excitations in the LSDW phase, which is a remarkable quantum signature of the field-induced TLL in Ising-like spin 1/2 1D antiferromagnets.

I Field-dependence of the incommensurability

Figure S1 shows the measured and calculated field-dependence of the incommensurability wave vector in the longitudinal spin density wave (LSDW) phase. We have calculated the ground state numerically by density matrix renormalization group (DMRG) for a finite size system consisting in 200 sites and obtained the local magnetization $M_{r}^{z}\equiv\langle S_{r}^{z}\rangle$. $\delta$ is determined from the peak that occurs in its Fourier transform $M^{z}(q)\equiv|\sum_{r}e^{iqr}M_{r}^{z}|$.

II Nature of the LSDW phase

A simple interpretation for the nature of the LSDW phase can be obtained by considering the softening of the lowest mode with increasing the external magnetic field. At zero magnetic field, the lowest energy excitation is the doubly degenerate transverse mode ($\Delta S^{z}=\pm 1$), which splits due to the Zeeman effect. The energy of an excitation with $\Delta S^{z}=+1$ decreases while that with $\Delta S^{z}=-1$ increases. When the excitation gap is closed, the domain wall excitations with $\Delta S^{z}=+1$ condensate and the quantum phase transition happens. The number of these domain walls proliferates with increasing the magnetic field so as to minimize the Zeeman energy. This increase of the number of domain walls is related with the decrease of their average distance of separation $1/\delta$. The domain walls actually have some intrinsic width as shown by nuclear magnetic resonance (NMR) measurements MK , so that this array of walls coincides with the spin density wave deduced from neutron diffraction. W e looked by neutron diffraction measurements for third order harmonics in the LSDW phase. Their presence would demonstrate a squaring of the sinusoidal amplitude modulation of the magnetic structure. However, such Bragg reflections could not be observed. The analysis of the error bars then shows that they must be at least 30 times smaller than the first order harmonics. This remains consistent with the expectations from the NMR line profile.

III Measured spin-dynamics at 6 T along (3, 0, $Q_{L}$)

Constant-$Q$ energy scans have been recorded at $\mu_{0}H=6$ T along the $c^{*}$ direction across points $\bm{Q}=(3,0,Q_{L})$ with $Q_{L}$ ranging between 1 [zone center (ZC) point] and 1.25. The resulting intensity map as a function of energy transfer and $Q_{L}$ is shown Fig. S2. It features an intense arch-like excitation, which takes the minimum at the incommensurate positions near the ZC position $(3,0,1)$.

IV Numerical simulations

In this section, we explain the method used to perform the numerical simulations. The principle of the calculations are the same as in Refs. faure2018 ; takayoshi2018 . Note that BaCo2V2O8 consists of the stacking of Co chains, each of which can be considered as a spin-1/2 Heisenberg model with Ising (easy-axis) anisotropy. Taking into account the interchain coupling by the mean field theory, we obtain an effective one-dimensional Hamiltonian

[TABLE]

Here, $x,y$ and $z$ coincide with the $a,b,c$ crystallographic axes. The local magnetization $\langle S^{z}_{n}\rangle$ is determined self-consistently. The parameters $J=3.05$ meV, $\Delta=1.9$ and $g$ factor along the $z$ axis $g_{zz}=6.07$ were determined so that they reproduce the neutron cross-section in zero-field at the scattering vector $\bm{Q}=(2,0,0)$ in the unit of $(2\pi/a,2\pi/b,2\pi/c)$, where $a,b,c$ are the lattice constants faure2018 . Note that a different convention was used in Ref. faure2018 , explaining the different numerical values of $J$: in the present paper, $J$ replaces $\epsilon J$ and $J\Delta$ replaces $J$ with $\Delta=1/\epsilon=1.9$. The differential neutron scattering cross section is represented as

[TABLE]

where $F(\bm{Q})$ is the magnetic form factor and $\bm{q},\bm{q}^{\prime}$ are the initial and final wave vectors, respectively ($\bm{Q}=\bm{q}-\bm{q}^{\prime}$). The dynamical structure factor $S_{\alpha\beta}(\bm{Q},\omega)$ is given as

[TABLE]

Here $C_{\rm ret}^{\alpha\beta}(\bm{r},t)$ is the retarded correlation function

[TABLE]

where $\vartheta_{\mathrm{step}}(t)$ is the step function. When the system has a rotational symmetry around the $z$ axis, as is the case of the $XXZ$ chain under a longitudinal field (S1), Eq. (S2) is recast into

[TABLE]

We first obtain the ground state using DMRG white1992 , then perform the time-evolution with time-dependent block decimation (TEBD) vidal2003 and calculate space-time correlation functions for the Hamiltonian (S1). In the calculations, the system size is $N=200$ and time interval is taken to be $0\leq t\leq 60J^{-1}$ with the discretization $dt=0.05J^{-1}$. The truncation dimension (i.e., the bond dimension of matrix product states) is $M=60$. For the Fourier transform in Eq. (S3), the summation is taken over the actual positions $\bm{r}$ of Co2+ ions.

V Effects of the interchain interaction in the numerical calculations

In this section, we examine the effects of the interchain interaction. Although it is likely that the interchain coupling consists of a complex set of interactions including further than the nearest neighbor klanjsek2015 , we consider for simplicity the interchain coupling only between the nearest neighbor sites and treated it in a mean-field theory as stated in the previous section,

[TABLE]

In Fig. S3, we show the results of numerical calculations around the AF position $\bm{Q}=(2,0,1)$ for $\mu_{0}H=4.2$ T [Fig. S3(a)] and 6 T [Fig. S3(b)] with varying $J^{\prime}$ from 0 to 0.17 meV. The agreement between the numerics and the experimental data, presented in the main article and in Fig. S2 is best for $J^{\prime}=0$ and worsens with increasing $J^{\prime}$. The deviation becomes larger as $H$ approaches the critical field as far as the spectral weight distribution is concerned. With increasing $J^{\prime}$, we can see that the Néel order (the weight at $Q_{L}=1$) is reinforced to the detriment of the incommensurate spin density wave. This is because an effective staggered field is induced by the Néel order through the interchain coupling. This results in an overestimation of the Néel order, and this effect becomes stronger as $H$ is closer to the phase transition point.