Introduction

Organic metals and semiconductors1 possess a variety of features not shared by inorganic materials, e.g., light, flexible, and toxic-element-free. They have been rapidly developed over the past decades for use in consumer electronic devices, such as organic transistors, light-emitting diodes, and piezo actuators. These accomplishments, in combination with recent evolutions of inorganic spintronics based on spin current physics, have promoted a new field, i.e., organic spintronics. Now significant efforts are being made to elucidate spin transport phenomena in organic semiconductors2,3,4. However, organic spintronics devices are actually not purely organic but are hybrid with inorganic materials, because the generation of spin current basically requires an inorganic magnetic electrode. In fact, attempts for exploiting organic materials as the spin current generator are quite limited5,6.

Here, we theoretically propose a microscopic mechanism of spin current generation in organic materials utilizing an archetypal antiferromagnet. Figure 1a provides a schematic illustration of the present spin current generation in the antiferromagnetic (AFM) state, where the up and down spins aligned on the molecular checker plate play a role of a spin-rectifier converting a heat-current driven by a thermal gradient, or an electron current by an electric field, into the spin current. When we rotate the external field with respect to the crystal axes in the two-dimensional plane, the direction of the generated spin current rotates in the opposite direction as shown in Fig. 1b. The directional dependence is strikingly different from the conventional spin Nernst and spin Hall effects7,8,9,10,11,12, in which the spin current always flows perpendicular to the field direction. As a platform of this phenomenon, we focus on an organic antiferromagnet \(\kappa\)-(BEDT-TTF)\({}_{2}\)Cu[N(CN)\({}_{2}\)]Cl (abbreviated as \(\kappa\)-Cl).

Fig. 1
figure 1

Schematic illustrations of the spin current generation. a Flows of up- and down-spin magnons (electrons) and spin current driven by a thermal gradient (an electric field) in the AFM state. The red and blue ellipses represent the two kinds of molecular dimers, forming a checker-plate-type lattice. The arrows in the dimers represent the localized spin moments. b Field-angle dependence of the spin current generation (green arrows) driven by a thermal gradient or an electric field (gray arrows)

Results

Crystal structure and model

\(\kappa\)-Cl is a well-studied insulator, showing a variety of cooperative phenomena, e.g., AFM ordering, insulator-to-metal transition, and superconductivity, at low temperatures and/or under pressures13,14,15,16,17,18. The crystal structure is composed of an alternate stacking of two-dimensional conducting BEDT-TTF (abbreviated as ET) layers and insulating Cu[N(CN)\({}_{2}\)]Cl layers. Figure 2a shows the molecular arrangement (called \(\kappa\)-type) in the conducting layer, where four ET molecules in the unit cell form two kinds of dimers with different orientations, termed \(A\) and \(B\), connected by a glide operation (mirror and half translation).

Fig. 2
figure 2

Schematic illustration of the lattice structure of \(\kappa\)-Cl and the energy bands. a Molecular arrangement in the two-dimensional conducting layer. The red and blue circles represent the two kinds of ET dimers, termed \(A\) and \(B\), respectively, in the unit cell. The green line denotes the glide plane perpendicular to the \(xy\) plane. b Network of the dominant electron transfer bonds, \({\rm{a}}\) (orange bold line), \({\rm{b}}\) (dotted line), \({\rm{p}}\) (solid line), and \({\rm{q}}\) (broken line). The gray circles represent the ET molecules, and the red and blue ellipses show the \(A\) and \(B\) dimers, respectively. The arrows represent the local spin moments in the AFM phase. We note that another glide plane exists when considering the layer stacking, but it does not affect our discussions. c Energy band dispersion composed of the frontier orbitals of ET molecules in the PM metallic phase with the transfer integrals \(({t}_{{\rm{a}}},{t}_{{\rm{p}}},{t}_{{\rm{q}}},{t}_{{\rm{b}}})=(-0.207,-0.102,0.043,-0.067)\) eV (green solid line) and that of the single-band picture in the large dimerization limit (broken line). The average electron number in the unit cell is \(6\) and the Fermi energy \({\varepsilon }_{{\rm{f}}}\) is shown. d Energy band dispersion in the AFM insulating phase with the intra-molecular Coulomb interaction \(U=1\) eV, within the self-consistent mean-field theory. e Contour map of the spin splitting subtracting the down-spin energy from the up-spin energy of the top band in d in the first Brillouin zone. The trajectory shows the symmetric lines in c and d

This class of organic materials is known to show a simple electronic structure composed of frontier molecular orbitals19,20. In the \(\kappa\)-type materials, the frontier orbitals in each ET dimer become strongly hybridized by the intra-dimer transfer integral shown in Fig. 2b, and constitute bonding and antibonding orbitals. They result in four bands as there are two dimers in the unit cell: two lower(higher)-energy bands are from the (anti-)bonding orbitals, as shown in Fig. 2c. The system has three electrons per two dimers on average, and hence, the four bands are three-quarter filled.

In the last few decades, extensive studies have been made for understanding the cooperative phenomena in this system21,22,23. Most of them, however, are based on the single-band picture, where the two fully occupied bands are disregarded (see the broken lines in Fig. 2c). This approach is justified in the large dimerization limit19, where the crystallographic distinction of the \(A\) and \(B\) dimers is lost. In other words, the glide symmetry in the molecular arrangement in the conducting layer was disregarded in the previous studies. In the following, we will discuss that the breaking of the glide symmetry by the AFM ordering plays an essential role in a peculiar spin current generation.

We investigate electronic structures and spin current transport properties of \(\kappa\)-Cl based on the Hubbard model, taking into account the distinct two types of dimers and the anisotropy in the transfer integrals between them19, \(({t}_{{\rm{a}}},{t}_{{\rm{p}}},{t}_{{\rm{q}}},{t}_{{\rm{b}}})=(-0.207,-0.102,0.043,-0.067)\) eV, evaluated by a first-principles calculation24 (see Fig. 2b). At three-quarter filling where the number of electrons in the unit cell is equal to \(6\), the ground state exhibits a metal-to-insulator transition from a paramagnetic (PM) phase to an AFM phase on increasing the intra-molecular Coulomb interaction \(U\)19,25.

Spin splitting

A crucial feature in the AFM state of \(\kappa\)-Cl is that up and down spins are situated on the dimers with the different orientations as shown in Fig. 2b, resulting in the glide symmetry breaking with respect to the \(yz\) plane. Here we consider the glide operation not acting on the spins. The molecular orientation makes the AFM state not invariant under the combination of time reversal and spatial translation operations, unlike simple N\(\acute{{\rm{e}}}\)el-type AFM state, e.g., on the square lattice. This situation gives rise to an energy band splitting depending on the spins, which has been overlooked previously. Figure 2d shows the band structure in the AFM state, calculated within the self-consistent mean-field theory (see Methods). The spin splitting appears in the whole Brillouin zone except on the \({k}_{x}\)-, \({k}_{y}\)-axes and the zone boundary as shown in Fig. 2e.

The origin of the spin splitting is understood from the real-space anisotropy induced by the AFM ordering as follows. Figure 3 shows the effective inter-dimer transfer integrals between the antibonding orbitals, calculated by the second-order perturbation with respect to the inter-orbital hybridizations (see Methods). In the PM phase, as shown in Fig. 3a, b, the \(A\) and \(B\) dimers show different real-space anisotropies owing to the molecular orientations, but the anisotropies are symmetric with respect to the glide operation and do not depend on the spin degree of freedom. In the AFM phase, in contrast, the transfer integrals for up-spin electrons on the \(A\) dimer (Fig. 3c) and down-spin electrons on the \(B\) dimer (Fig. 3f) are enhanced, whereas their counterparts (Fig. 3d, e) are reduced. This spin-dependent anisotropy leads to the spin splitting.

Fig. 3
figure 3

Spatial anisotropy of inter-dimer transfer integrals. a, b Effective transfer integrals between the central dimer (\(A\) in a and \(B\) in b) and the surrounding ones, obtained from the second-order perturbation processes with respect to the bonding-antibonding orbital hybridizations, in the PM phase. cf Effective transfer integrals calculated likewise for up- (c, d) and down-spin (e, f) electrons in the AFM phase (the local magnetic moment is about \(0.168\ \hslash\)). The areas of the red (blue) shaded circles represent the amplitudes of positive (negative) transfer integrals. The amplitudes are shown in the circles in unit of meV

The real-space anisotropies also show up in the effective spin exchange interactions in the Heisenberg model, derived from the above Hubbard model. Note that the system retains \(SU(2)\) symmetry because of the absence of the spin-orbit coupling. Figure 4a shows the spatial distributions of the nearest-neighbor (NN) exchange interactions \(J\) and \(J^{\prime}\), and the next-nearest-neighbor (NNN) interactions \(K\) and \(K^{\prime}\). \(K\) and \(K^{\prime}\) arise from fourth-order perturbation processes with respect to the NN transfer integrals (see Methods). As shown in Fig. 4b, \(K^{\prime}\) becomes much smaller than \(K\) for realistic parameters. Then, the AFM magnon dispersion of the Heisenberg model exhibits a spin splitting as shown in Fig. 4c, where we take \(K=2\) meV and \(K^{\prime} =0\) for simplicity, and \(J=80\) meV and \(J^{\prime} =20\) meV14 (see Methods). Similar spin splitting was reported in non-centrosymmetric systems with the spin-orbit coupling26,27, but the present mechanism requires neither non-centrosymmetry nor the spin-orbit coupling.

Fig. 4
figure 4

Effective NNN exchange interactions, magnon dispersions, and heat-spin current conversion in the AFM insulating state. a Real-space distribution of the NN exchange interactions \(J\) (solid lines) and \(J^{\prime}\) (broken lines) between \(A\) and \(B\) dimers (left panel) and those of the NNN exchange interactions \(K\) (purple solid lines) and \(K^{\prime}\) (purple broken lines) between \(A\) dimers (middle panel) and \(B\) dimers (right panel). b, \({t}_{\rm{q}}\) dependences of \(K\) and \(K^{\prime}\) at \(U=1\) eV. The red arrow represents the value of \({t}_{\rm{q}}\) in \(\kappa\)-Cl. c Magnon dispersions at \((J,J^{\prime} ,K,K^{\prime} )=(80,20,2,0)\) meV within the linear spin-wave theory. The inset shows a contour map of the spin splitting between the up- and down-spin magnons in the first Brillouin zone. d \((T,K)\) dependences of the spin current conductivity under a thermal gradient, \({\chi }_{xy}^{{\rm{SQ}}}\). The other exchange interactions and the damping factor are \((J,J^{\prime} ,K^{\prime} ,\eta )\) = \((80,20,0,1)\) meV. e \(K\) dependences of the heat-spin current conversion rate \(\alpha (=|2J{\chi }_{xy}^{{\rm{SQ}}}/\hslash {\kappa }_{yy}|)\) at \({k}_{{\rm{B}}}T=0.5\) meV and \(1\) meV, where \({\kappa }_{yy}\) is the thermal conductivity along the \(y\)-axis. The red arrow represents the value of \(K\) in \(\kappa\)-Cl

Spin current by a thermal gradient

The spin-split magnon excitations lead to a spin current generation. Figure 4d shows the off-diagonal spin current conductivity, along the \(x\)-axis with respect to the thermal gradient along the \(y\)-axis, \({\chi }_{xy}^{{\rm{SQ}}}\), as a function of temperature \(T\) and the exchange interaction \(K\), calculated by the linear response theory (see Methods). The range of \(T\) is chosen well below the N\(\acute{{\rm{e}}}\)el temperature of \(\kappa\)-Cl, \(23\) K. The polarization of the spin current is parallel to the AFM moment, and the damping factor \(\eta\) is fixed at \(1\) meV. We obtain nonzero \({\chi }_{xy}^{{\rm{SQ}}}\) for \(T> 0\) and \(K> 0\), which monotonically increases in proportion to \({T}^{2}\) and \(K\). Remarkably, the conductivity tensor \({{\boldsymbol{\chi }}}^{{\rm{SQ}}}\) is symmetric, \({\chi }_{xy}^{{\rm{SQ}}}={\chi }_{yx}^{{\rm{SQ}}}\), with vanishing diagonal elements, \({\chi }_{xx}^{{\rm{SQ}}}={\chi }_{yy}^{{\rm{SQ}}}=0\). This leads to the peculiar field-angle dependence that we showed in Fig. 1b, which is distinct from the conventional spin Nernst effect where the spin current is always perpendicular to the thermal gradient.

This spin current generation is a direct consequence of the magnon dispersion in Fig. 4c which indicates that the up-spin magnon has high mobility along the \((1,1)\) and \((-1,-1)\) directions, while the down-spin magnon has along the \((1,-1)\) and \((-1,1)\) directions. When the temperature gradient is applied along the \(y\)-axis as shown in Fig. 1a, the up- and down-spin magnons are rectified toward the \((1,1)\) and \((1,-1)\) directions, respectively, in a symmetric way. Accordingly, a pure spin current, where the net magnon current is canceled out, is generated along the \(x\)-axis. This gives rise to the positive transverse \({{\boldsymbol{\chi }}}^{{\rm{SQ}}}\) in Fig. 4d. On the other hand, if the temperature gradient is parallel to the \((1,1)\) direction, while the transverse component disappears, a net up-spin magnon current is generated parallel to the field as a result of the incomplete cancellation (see Fig. 1b). This provides a finite longitudinal component of \({{\boldsymbol{\chi }}}^{{\rm{SQ}}}\) in the rotated coordinate, which is consistent with the symmetric form of the conductivity tensor.

We find that \({\chi }_{xy}^{{\rm{SQ}}}\) is inversely proportional to the damping factor \(\eta\) and diverges in the clean limit (\(\eta =0\)), in analogy with the diagonal thermal conductivities \({\kappa }_{xx}\) and \({\kappa }_{yy}\) (see Supplementary Fig. 1a). This indicates that the ratio \(\alpha \equiv |2J{\chi }_{\mu \nu }^{{\rm{SQ}}}/\hslash {\kappa }_{\nu \nu }|\), which is used in the literatures as the conversion rate from the heat-current to the spin current, does not depend on \(\eta\), however, depends on the field angle. Therefore, we choose \(\mu =x\) and \(\nu =y\) since \({\kappa }_{\nu \nu }\) is largest in this direction, considering its implication as the conversion rate. Figure 4e shows \(K\) dependences of \(\alpha\) at \({k}_{{\rm{B}}}T=0.5\) meV and \(1\) meV linearly increasing with \(K\), but almost independent of \(T\). The heat-spin current conversion efficiency reaches \(\sim 5\)\(\%\) for the case of \(\kappa\)-Cl, which is close to one-quarter of that in Pt due to the strong spin-orbit coupling28.

Spin current by an electric field

Now we propose another way of a spin current generation, in carrier doped metallic regions. The carrier doping has recently been realized experimentally17,18. We here focus on the electron-doping case where the AFM metallic state is stable in our model. Figure 5a shows the off-diagonal spin current conductivity induced by the electric field, \({\chi }_{xy}^{{\rm{SC}}}(={\chi }_{yx}^{{\rm{SC}}})\), as a function of the Coulomb interaction \(U\) and the number of electrons in the unit cell \(n\) in the ground state (see Methods). \({\chi }_{xy}^{{\rm{SC}}}\) is zero in the PM metallic and the AFM insulating phases, while it turns finite in the AFM metallic phase where the Fermi energy lies in the top band in Fig. 2d, whose spin splitting was shown in Fig. 2e. We note that the sign of \({\chi }_{xy}^{{\rm{SC}}}\) changes around \(n=6.2\), associated with the change in the Fermi surface topology as shown in the insets of Fig. 5a. The conductivity tensor is also symmetric with zero diagonal components and inversely proportional to the damping factor (see Supplementary Fig. 1b). This means that the field-angle dependence is the same as that of \(\chi_{xy}^{\mathrm{{SQ}}}\) in the insulating case. It is comprehensible from the spin-dependent anisotropy of the electron transfers in Fig. 3c–e, reflected in the anisotropy in the spin-split band in Fig. 2e with the same character as the magnon band in Fig. 4c. We define the charge-spin current conversion rate by \(\beta \equiv |2e{\chi }_{yx}^{{\rm{SC}}}/\hslash {\sigma }_{xx}|\), in analogy with \(\alpha\) above (the electrical conductivity \({\sigma }_{\nu \nu }\) becomes largest in the \(\nu =x\) direction due to the quasi-one-dimensional Fermi surfaces). As shown in Fig. 5b, in the lightly doped region with small \({\sigma }_{xx}\), \(\beta\) becomes relatively large and approaches \(7\)\(\%\), comparable to the spin Hall effect in Pt29, while in the highly doped region, it decreases because of the suppression of the AFM ordering and the spin splitting.

Fig. 5
figure 5

Charge-spin current conversion in the electron-doped AFM metallic state. a \((U,n)\) dependences of the spin current conductivity to an electric field, \({\chi }_{xy}^{{\rm{SC}}}\). The broken line represents the phase boundary between the PM and AFM metallic phases, and the blue thick line shows the AFM insulating phase at three-quarter filling. The damping factor is \(\gamma =1\) meV. The insets show the Fermi surface structures of up-spin (red) and down-spin (blue) electrons at \((U,n)=(1\ {\rm{eV}},6.1)\) and \((1\ {\rm{eV}},6.4)\). The gray shaded areas denote the occupied states. b \(n\) dependences of the charge-spin current conversion rate \(\beta (=|2e{\chi }_{yx}^{{\rm{SC}}}/\hslash {\sigma }_{xx}|)\), \(|{\chi }_{yx}^{{\rm{SC}}}|\), and the diagonal electrical conductivity \({\sigma }_{xx}\) at \(U=1\) eV

The spin current generation in our mechanism is expected to be observed at sufficiently low temperature compared to the N\(\acute{{\rm{e}}}\)el temperature, which is not determined for doped \(\kappa\)-Cl. We anticipate it to be lower than the undoped case of \(23\) K, but to remain the same order in the lightly doped region30.

Discussion

The present spin current generation is strikingly different from the conventional spin Nernst and spin Hall effects. In the conventional mechanisms, a spin current is activated by the spin–orbit coupling in non-centrosymmetric lattice structures. The conductivity tensor is antisymmetric, namely, the generated spin current is always perpendicular to the applied field direction and the conversion rate is invariant under the rotation of the field. The transverse conductivity converges to a finite value in the clean limit because of the dominant inter-band contributions7,8. However, the strong spin–orbit coupling also disturbs the spin polarization of carriers via the spin-flipping process.

In stark contrast, the present mechanism requires neither the spin-orbit coupling nor spatial inversion symmetry breaking. The spin current conductivity is described by the symmetric tensor which results in the peculiar field-angle dependence shown in Fig. 1b, and diverges in the clean limit due to the intra-band contributions. In \(\kappa\)-type ET systems, the Dzyaloshinskii-Moriya interaction due to the spin-orbit coupling is estimated to be a few Kelvin15,31, which is much smaller than the NNN exchange interaction \(K\). Furthermore, this class of organic charge transfer salts is known to have relatively less impurities and lattice disorders compared to inorganic crystals and organic polymers. Indeed, the specific heat and thermal transport measurements32,33 suggest that the low temperature properties are well described by intrinsic contributions from electronic charge and spin degrees of freedom. These facts ensure a long spin lifetime in \(\kappa\)-Cl, which facilitates the experimental detection. Although the phenomenon has a similarity with the spin current generation in ferromagnetic metals in the sense that the time reversal symmetry is lost, the net magnetization is absent in our system; this enables us to generate a pure spin current in contrast to the spin-polarized current in ferromagnets and has the advantage of small field leakage as discussed in AFM spintronics. These considerations lead us to conclude that our proposal provides a new type of spin current generation essentially distinct from the other existing mechanisms.

As a recent experimental progress relevant to our proposal, the three-dimensional AFM structures in several \(\kappa\)-type ET systems have been determined by the detailed analyses of the magnetization processes15. It was found that \(\kappa\)-Cl and \(\kappa\)-Br show the same intra-layer AFM structure as shown in Fig. 2b, but different inter-layer stackings; the “in-phase” stacking, where the inter-layer NN spins are ferromagnetically aligned, is realized in \(\kappa\)-Cl, while \(\kappa\)-Br shows the “anti-phase” stacking. This difference will give an effective way to verify the present spin current generation because in our mechanism the sign of generated spin current is reversed by the reversal of the AFM moment. This means that a net spin current is expected in \(\kappa\)-Cl while it will be canceled out in \(\kappa\)-Br. In addition, our mechanism also has the inverse effect similar to the inverse spin Nernst or spin Hall effect, i.e., the generation of a thermal gradient or an electrical voltage by spin current injection parallel to the AFM ordered ET layers, which will give another experimental approach.

The present spin current generation arises from AFM ordering in spatially-oriented molecular orbitals. The molecular orbital degrees of freedom are fundamental and ubiquitous in organic materials. Meanwhile, similar orbital degrees of freedom are also found in inorganic materials, such as transition metal and rare-earth compounds. Thus, our new mechanism can be applied to a wide range of AFM materials. In this perspective, therefore, our finding strikes out a new direction of materials exploration for spintronics without relying on the spin-orbit coupling.

Methods

Mean-field approximation

The Hamiltonian of the Hubbard model on the \(\kappa\)-type lattice is given by

$$\begin{array}{ccc}{{\mathcal{H}}}_{{\rm{Hubb}}}&=&U{\sum }_{i\mu }{n}_{i\mu \uparrow }{n}_{i\mu \downarrow }+{t}_{{\rm{a}}}{\sum }_{i\sigma }({c}_{ia\sigma }^{\dagger }{c}_{ib\sigma }+{\rm{H.c.}})\\ &&+{\sum }_{\langle ij\rangle \mu \mu ^{\prime} \sigma }{t}_{ij}^{\mu \mu ^{\prime} }({c}_{i\mu \sigma }^{\dagger }{c}_{j\mu ^{\prime} \sigma }+{\rm{H.c.}}),\end{array}$$
(1)

where \({c}_{i\mu \sigma }\) and \({n}_{i\mu \sigma }(={c}_{i\mu \sigma }^{\dagger }{c}_{i\mu \sigma })\) are the annihilation operator and the number operator of an electron with a spin \(\sigma (=\uparrow ,\downarrow )\), on the frontier orbital of molecular site \(\mu (=a,b)\) in the \(i\)th dimer, respectively, \(U\) is the intra-molecular Coulomb interaction, and \({t}_{{\rm{a}}}\) and \({t}_{ij}^{\mu \mu ^{\prime} }\) are the inter-molecular transfer integrals shown in Fig. 2b. We treat the Coulomb interaction term within the mean-field approximation as \({n}_{i\mu \uparrow }{n}_{i\mu \downarrow }\simeq {n}_{i\mu \uparrow }\langle {n}_{i\mu \downarrow }\rangle +\langle {n}_{i\mu \uparrow }\rangle {n}_{i\mu \downarrow }-\langle {n}_{i\mu \uparrow }\rangle \langle {n}_{i\mu \downarrow }\rangle\), and determine the expectation values self-consistently so as to minimize the total energy of the system.

Effective electron transfer integrals

We divide the mean-field Hamiltonian in the AFM phase into three terms as \({{\mathcal{H}}}_{{\rm{MF}}}={{\mathcal{H}}}_{{\rm{intra}}}+{{\mathcal{H}}}_{{\rm{inter}}}+{{\mathcal{H}}}_{{\rm{AFM}}}\), where the first and second terms represent the intra-orbital and inter-orbital transfer integrals, respectively, and the third term is the local AFM field. By taking the linear combinations of the original electron operators, we define the annihilation operator of an electron in the antibonding (bonding) orbital on the \(i\)th dimer as \({\tilde{c}}_{i\alpha (\beta )\sigma }=({c}_{ia\sigma }-(+){c}_{ib\sigma })/\sqrt{2}\), and the three terms are given by

$${{\mathcal{H}}}_{{\rm{intra}}}={t}_{{\rm{a}}}{\sum }_{i\sigma }({\tilde{n}}_{i\beta \sigma }-{\tilde{n}}_{i\alpha \sigma })+{\sum }_{\langle ij\rangle \nu \sigma }({\tau }_{ij}^{\nu \nu }{\tilde{c}}_{i\nu \sigma }^{\dagger }{\tilde{c}}_{j\nu \sigma }+{\rm{H.c.}}),$$
(2)
$${{\mathcal{H}}}_{{\rm{inter}}}={\sum }_{\langle ij\rangle \nu \sigma }({\tau }_{ij}^{\nu \bar{\nu }}{\tilde{c}}_{i\nu \sigma }^{\dagger }{\tilde{c}}_{j\bar{\nu }\sigma }+{\rm{H.c.}}),$$
(3)
$${{\mathcal{H}}}_{{\rm{AFM}}}=\frac{U\delta }{4}({\sum }_{i(\in B)\nu \sigma }-{\sum }_{i(\in A)\nu \sigma })\sigma {\tilde{n}}_{i\nu \sigma },$$
(4)

where the number operator is given by \({\tilde{n}}_{i\nu \sigma }={\tilde{c}}_{i\nu \sigma }^{\dagger }{\tilde{c}}_{i\nu \sigma }\), and \(\bar{\nu }=\beta\) (\(\alpha\)) for \(\nu =\alpha\) (\(\beta\)). The transfer integral between the neighboring dimers is given by \({\tau }_{ij}={\cal{U}}{t}_{ij}{\cal{U}}^{{\rm{T}}}\), by using the two-by-two unitary matrix \(\cal U\) satisfying \({({\tilde{c}}_{\alpha \sigma },{\tilde{c}}_{\beta \sigma })}^{{\rm{T}}}={\cal{U}}{({c}_{a\sigma },{c}_{b\sigma })}^{{\rm{T}}}\). The amplitude of the local AFM field is given by \(\delta =\langle {\tilde{n}}_{i\in A\uparrow }\rangle -\langle {\tilde{n}}_{i\in A\downarrow }\rangle =\langle {\tilde{n}}_{i\in B\downarrow }\rangle -\langle {\tilde{n}}_{i\in B\uparrow }\rangle\), where \({\tilde{n}}_{i\sigma }={\sum }_{\nu }{\tilde{n}}_{i\nu \sigma }\).

We treat \({{\mathcal{H}}}_{{\rm{inter}}}\) as the perturbation term and calculate the effective transfer integrals for the bonding and antibonding orbitals up to \({\mathcal{O}}({{\mathcal{H}}}_{{\rm{inter}}}^{2})\). In the \({\bf{k}}\) space, the mean-field Hamiltonian is described by the matrix form as \({{\mathcal{H}}}_{{\rm{MF}}}={\sum }_{{\bf{k}}\sigma }{{\bf{d}}}_{{\bf{k}}\sigma }^{\dagger }({H}_{{\bf{k}}\sigma }^{(0)}+{V}_{{\bf{k}}\sigma }){{\bf{d}}}_{{\bf{k}}\sigma }\), where \({H}_{{\bf{k}}\sigma }^{(0)}\) and \({V}_{{\bf{k}}\sigma }\) are the unperturbed and perturbed terms, respectively, given by \(4\times 4\) matrices. \({{\bf{d}}}_{{\bf{k}}\sigma }\) is the vector of the annihilation operators of the Bloch states, which is chosen so as to diagonalize the unperturbed term as \({\hat{H}}_{{\bf{k}}\sigma }^{(0)}|{\bf{k}}{\nu }_{\sigma }^{\xi }\rangle ={\varepsilon }_{{\bf{k}}\nu {\sigma} }^{\xi }|{\bf{k}}{\nu }_{\sigma }^{\xi }\rangle\), where \(\xi\)(\(=1,2\)) indicates the two bands in the bonding and antibonding bands each originating from the two dimers in the unit cell. The second-order perturbation term \({H}_{{\bf{k}}\sigma }^{(2)}\) is decomposed into two \(2\times 2\) matrices for the antibonding (\(\alpha\)) and bonding bands (\(\beta\)) as \({H}_{{\bf{k}}\sigma }^{(2)}={h}_{{\bf{k}}\alpha \sigma }^{(2)}\oplus {h}_{{\bf{k}}\beta \sigma }^{(2)}\). The matrix element of \({h}_{{\bf{k}}\nu \sigma }^{(2)}\) is given by

$${h}_{\nu ;\xi \xi ^{\prime} }^{(2)}={\sum }_{\eta =1,2}\frac{\langle {\nu }^{\xi }|\hat{V}|{\bar{\nu }}^{\eta }\rangle \langle {\bar{\nu }}^{\eta }|\hat{V}|{\nu }^{\xi ^{\prime} }\rangle }{{\varepsilon }_{\nu }^{\xi ^{\prime} }-{\varepsilon }_{\bar{\nu }}^{\eta }},$$
(5)

where the indices \({\bf{k}}\) and \(\sigma\) are omitted for simplicity. By the Fourier transformation of \({H}_{{\bf{k}}\sigma }^{(2)}\), we obtain the effective transfer integrals shown in Fig. 3.

Next-nearest-neighbor exchange interactions

From the Hubbard model in Eq. (1), we derive the effective NNN exchange interaction in the restricted space where each antibonding orbital is occupied by one hole due to the strong Coulomb interaction \(U\). The NNN exchange interaction is derived from the fourth-order perturbation process with respect to the inter-dimer transfer integrals, which is given by

$${{\mathcal{H}}}^{(4)}={\mathcal{P}}{\mathcal{V}}{\left(\frac{1}{{E}_{I}-{{\mathcal{H}}}^{(0)}}{\mathcal{Q}}{\mathcal{V}}\right)}^{3}|I\rangle\langle I|,$$
(6)

where \({\mathcal{P}}\) and \({\mathcal{Q}}\) are the projection operators onto inside and outside of the restricted space, respectively, and satisfy \({\mathcal{P}}+{\mathcal{Q}}=1\), \({\mathcal{V}}\) is the perturbation given by the third term in Eq. (1), \({{\mathcal{H}}}^{(0)}\) is the unperturbed Hamiltonian given by the first and second terms in Eq. (1), and \({E}_{I}\) is the energy of the initial eigenstate \(|I\rangle\) of \({{\mathcal{H}}}^{(0)}\). The resultant exchange interaction on the NNN bond between the \(A\) dimers denoted by \(K\) in the middle panel of Fig. 4a is given by

$${{\mathcal{H}}}_{ijk}=\tilde{J}\left({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}+{{\bf{S}}}_{j}\cdot {{\bf{S}}}_{k}-\frac{1}{2}\right)+K\left({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{k}-\frac{1}{4}\right),$$
(7)

where the indices \(ijk\) denote the three neighboring dimers, \({{\bf{S}}}_{i}\) is the spin operator of the \(i\)th dimer, and \(K\) is the NNN exchange constant. \(\tilde{J}\) is the NN exchange constant arising from the fourth-order perturbation process, which does not contribute to the magnon splitting. The details of the fourth-order process and the explicit form of \(K\) (and \(K^{\prime}\)) are given in Supplementary Note 2.

Linear spin-wave approximation

The effective Heisenberg model involving the NNN exchange interaction is given by

$$\begin{array}{ccc}&&{{\mathcal{H}}}_{{\rm{Heis}}}=J{\sum }_{\langle ij\rangle }{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}+J^{\prime} {\sum }_{\langle ij\rangle ^{\prime} }{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}\\ &&\,+K{\sum }_{\langle \langle ij\rangle \rangle }{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}+K^{\prime} {\sum }_{\langle \langle ij\rangle \rangle ^{\prime} }{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j},\end{array}$$
(8)

where \(\langle ij\rangle\) and \(\langle ij\rangle ^{\prime}\) stand for the diagonal and horizontal NN bonds on the equilateral triangular lattice, \(\langle \langle ij\rangle \rangle\) and \(\langle \langle ij\rangle \rangle ^{\prime}\) are the NNN bonds shown in Fig. 4a. By using the Holstein-Primakoff transformation, we obtain the linear spin-wave Hamiltonian given by

$$\begin{array}{ccc}{{\mathcal{H}}}_{{\rm{Heis}}}\simeq {{\mathcal{H}}}_{{\rm{LSW}}}&=&\frac{1}{2}{\sum }_{{\bf{k}}}\left[{A}_{{\bf{k}}}{a}_{{\bf{k}}}^{\dagger }{a}_{{\bf{k}}}+{B}_{{\bf{k}}}{b}_{-{\bf{k}}}^{\dagger }{b}_{-{\bf{k}}}\right.\\ &&\,\left.+{C}_{{\bf{k}}}({a}_{{\bf{k}}}^{\dagger }{b}_{-{\bf{k}}}^{\dagger }+{a}_{{\bf{k}}}{b}_{-{\bf{k}}})\right],\end{array}$$
(9)

where \({a}_{{\bf{k}}}\) and \({b}_{{\bf{k}}}\) are the Fourier transforms of the annihilation operators of magnons on the \(A\) and \(B\) dimers, respectively. The coefficients are given by

$$\begin{array}{ccc}{A}_{{\bf{k}}}&=&4J+2J^{\prime} [\cos ({\bf{k}}\cdot {{\bf{a}}}_{x})-1]\\ &+&2K[\cos ({\bf{k}}\cdot ({{\bf{a}}}_{x}+{{\bf{a}}}_{y}))-1]\\ &+&2K^{\prime} [\cos ({\bf{k}}\cdot ({{\bf{a}}}_{x}-{{\bf{a}}}_{y}))-1],\end{array}$$
(10)
$$\begin{array}{ccc}{B}_{{\bf{k}}}&=&4J+2J^{\prime} [\cos ({\bf{k}}\cdot {{\bf{a}}}_{x})-1]\\ &+&2K[\cos ({\bf{k}}\cdot ({{\bf{a}}}_{x}-{{\bf{a}}}_{y}))-1]\\ &+&2K^{\prime} [\cos ({\bf{k}}\cdot ({{\bf{a}}}_{x}+{{\bf{a}}}_{y}))-1],\end{array}$$
(11)

and

$$\begin{array}{ccc}{C}_{{\bf{k}}}&=&2J\left[\cos ({\bf{k}}\cdot ({{\bf{a}}}_{x}+{{\bf{a}}}_{y})/2)\right.\\ &+&\left.\cos ({\bf{k}}\cdot ({{\bf{a}}}_{x}-{{\bf{a}}}_{y})/2)\right],\end{array}$$
(12)

where \({{\bf{a}}}_{x}\) and \({{\bf{a}}}_{y}\) are the primitive translational vectors. \({{\mathcal{H}}}_{{\rm{LSW}}}\) is easily diagonalized by the standard Bogoliubov transformation, and the magnon energy dispersion shown in Fig. 4c is obtained.

Spin current conductivity to a thermal gradient

The spin current and energy current operators in the magnon system34 are given by

$${{\bf{J}}}_{{S}^{z}}=\frac{1}{i\hslash }[{{\bf{P}}}_{{S}^{z}},{{\mathcal{H}}}_{{\rm{LSW}}}]$$
(13)

and

$${{\bf{J}}}_{E}=\frac{1}{i\hslash }[{{\bf{P}}}_{E},{{\mathcal{H}}}_{{\rm{LSW}}}],$$
(14)

respectively. \({{\bf{P}}}_{{S}^{z}}\) and \({{\bf{P}}}_{E}\) are the spin polarization and the energy polarization operators defined by \({{\bf{P}}}_{{S}^{z}}=\hslash {\sum }_{i}{S}_{i}^{z}{{\bf{R}}}_{i}\) and \({{\bf{P}}}_{E}={\sum }_{i}{{\mathcal{H}}}_{i}{{\bf{R}}}_{i}\), respectively, where \({{\bf{R}}}_{i}\) is the position vector of the center of the \(i\)th dimer and \({{\mathcal{H}}}_{i}\) is the local energy density defined by \({{\mathcal{H}}}_{{\rm{LSW}}}={\sum }_{i}{{\mathcal{H}}}_{i}\), by the Fourier transformation of Eq. (9). In the magnon system where the chemical potential is zero, the heat-current operator \({{\bf{J}}}_{Q}\) is identical to the energy current operator \({{\bf{J}}}_{E}\). We note that the spin is a conserved quantity and the spin current is well defined here since our model does not include the spin-orbit coupling. In the linear response theory, the spin current conductivity to a static thermal gradient is given by

$$T{\chi }_{\mu \nu }^{{\rm{SQ}}}=\mathop{\mathrm{lim}}\limits_{\omega \to 0}\frac{{Q}_{\mu \nu }^{{\rm{SQ}}}(\omega )-{Q}_{\mu \nu }^{{\rm{SQ}}}(0)}{i\omega },$$
(15)

where \(\mu\) and \(\nu\) represent the spatial axes \(x\) and \(y\). The spin-current-heat-current response function \({Q}_{\mu \nu }^{{\rm{SQ}}}(\omega )\) is given by the Kubo formula

$${Q}_{\mu \nu }^{{\rm{SQ}}}(\omega )=\frac{i}{\hslash V}\int _{0}^{\infty }dt{e}^{it(\omega +i\eta )}{\langle [{J}_{{S}^{z}}^{\mu }(t),{J}_{Q}^{\nu }]\rangle }_{{\rm{eq}}},$$
(16)

where \({{\bf{J}}}_{{S}^{z}}(t)\) is the Heisenberg representation of the spin current operator, \(\eta\) is the damping factor, \(V\) is the volume of the system, and \({\langle \cdots \rangle }_{{\rm{eq}}}\) represents the thermal average under the temperature \(T\).

Spin current conductivity to an electric field

The spin current and charge current operators are defined by

$${{\bf{J}}}_{{s}^{z}}=\frac{1}{i\hslash }[{{\bf{P}}}_{{s}^{z}},{{\mathcal{H}}}_{{\rm{MF}}}]$$
(17)

and

$${\bf{J}}=\frac{1}{i\hslash }[{\bf{P}},{{\mathcal{H}}}_{{\rm{MF}}}],$$
(18)

respectively. \({{\bf{P}}}_{{s}^{z}}\) and \({\bf{P}}\) are the spin \({s}^{z}\) polarization and the electric polarization operators defined by \({{\bf{P}}}_{{s}^{z}}=\hslash {\sum }_{i}{s}_{i}^{z}{{\bf{r}}}_{i}\) and \({\bf{P}}=-e{\sum }_{i}{{\bf{r}}}_{i}\), respectively, where \({s}_{i}^{z}=\frac{{n}_{i\uparrow }-{n}_{i\downarrow }}{2}\) is the spin operator of the \(i\)th molecule at the position vector \({{\bf{r}}}_{i}\). The spin current conductivity to an static electric field is given by

$${\chi }_{\mu \nu }^{{\rm{SC}}}=\mathop{\mathrm{lim}}\limits_{\omega \to 0}\frac{{Q}_{\mu \nu }^{{\rm{SC}}}(\omega )-{Q}_{\mu \nu }^{{\rm{SC}}}(0)}{i\omega }.$$
(19)

The spin-current-charge-current response function \({Q}_{\mu \nu }^{{\rm{SC}}}(\omega )\) is given by the Kubo formula

$${Q}_{\mu \nu }^{{\rm{SC}}}(\omega )=\frac{i}{\hslash V}\int _{0}^{\infty }dt{e}^{it(\omega +i\gamma )}{\langle [{J}_{{s}^{z}}^{\mu }(t),{J}^{\nu }]\rangle }_{0},$$
(20)

where \({{\bf{J}}}_{{S}^{z}}(t)\) is the Heisenberg representation of the spin current operator, \(\gamma\) is the damping factor, and \({\langle \cdots \rangle }_{0}\) represents the average with respect to the ground state.