Thank you for visiting nature.com. You are using a browser version with
limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off
compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site
without styles and JavaScript.

Subjects

Abstract

Toroidal quantum states are most promising for building quantum computing and information storage devices, as they are insensitive to homogeneous magnetic fields, but interact with charge and spin currents, allowing this moment to be manipulated purely by electrical means. Coupling molecular toroids into larger toroidal moments via ferrotoroidic interactions can be pivotal not only to enhance ground state toroidicity, but also to develop materials displaying ferrotoroidic ordered phases, which sustain linear magneto–electric coupling and multiferroic behavior. However, engineering ferrotoroidic coupling is known to be a challenging task. Here we have isolated a {CrIIIDyIII6} complex that exhibits the much sought-after ferrotoroidic ground state with an enhanced toroidal moment, solely arising from intramolecular dipolar interactions. Moreover, a theoretical analysis of the observed sub-Kelvin zero-field hysteretic spin dynamics of {CrIIIDyIII6} reveals the pivotal role played by ferrotoroidic states in slowing down the magnetic relaxation, in spite of large calculated single-ion quantum tunneling rates.

Introduction

The magnetic behavior of molecular coordination complexes continues to intrigue scientists around the world, revealing many interesting physical properties and offering many potential applications such as new storage and information processing technologies1,2,3. Fundamental research into, e.g., single-molecule magnets (SMMs)4,5,6, spin-crossover7, and magnetic systems with toroidal moments8,9,10,11,12,13,14,15,16,17,18,19 are recognized as important areas of molecular magnetism. SMMs exhibit slow relaxation of the magnetization, acting as nano-magnets below their blocking temperature2, 4,5,6. Molecular coordination complexes with a toroidal arrangement of local magnetic moments are rare, but are of great interest as they have several potential applications such as quantum computation, molecular spintronics devices, and the development of magneto–electric coupling for multiferroic materials20,21,22. Toroidal moments at a molecular level were first predicted12 and observed9, 17 in 2008, in strongly anisotropic metal complexes with a ring topology, having in-plane magnetic axes tangential to the ring, and weak nearest neighbor exchange coupling of the appropriate sign. In particular, the observation of a toroidal texture9, 17 in the ground state of the {DyIII3} triangular system8, generated great interest in this area, with state-of-the-art theoretical9, 14 and experimental10, 17, 18 techniques being employed to probe toroidal states of this prototype and subsequent related molecules. Other than triangular {Dy3} complexes, the toroidal moments have been reported as well for rhombus {Dy4}23,24,25 and {Dy6} wheel26 complexes.

If two molecular rings exhibiting toroidal states are connected, e.g., via a 3d ion, then coupling between the two toroidal moments leading to an enhancement of the collective toroidal moment may occur, which is a prerequisite to achieve a molecular ferrotoroidically ordered phase and the development of molecule-based multiferroics11, 13, 18, 27. In order to isolate materials with the above mentioned properties, we target coordination complexes containing anisotropic 4f ions. In contrast to the great deal of work on the synthesis of 3d–4f coordination complexes using 3d ions such as MnIII, FeIII, and CoII28,29,30,31,32, there have been few reports of studies, both structurally and magnetically, on mixed Cr(III)–Ln(III) systems. We have recently shown, however, that the combination of the isotropic CrIII ion and the anisotropic DyIII ion resulted in a family of SMMs with long relaxation times, relative to other lanthanide-based SMMs33,34,35,36. With this in mind, we have chosen to expand our studies, utilizing chromium(III) nitrate, with various lanthanide(III) ions, with carboxylic acid pro-ligands.

Herein, we describe the synthesis, structural characterization, and magnetic properties of a heterometallic complex of formula [CrIIIDyIII6(OH)8(ortho-tol)12(NO3)(MeOH)5]∙3MeOH (1), where ortho-tol = ortho-toluate. Complex 1 displays slow magnetic relaxation and SMM behavior at temperatures below 2 K. Furthermore, in {CrIIIDyIII6} we find, for the first time, a ferrotoroidically coupled ground state fully determined by dipolar coupling between the two con-rotating toroidal triangles. Our observations on 1 have been compared to earlier reported studies on coupled molecular {DyIII3} toroids10, 11, 18. The ferrotoroidically coupled ground state thus leads to an enhanced toroidal moment in the ground state for the {CrIIIDyIII6} complex, which is shown to play a central role in the observed magnetization dynamics featuring zero-field hysteresis.

Results

Synthesis and magnetic properties

Compound 1 was synthesized by the reaction of Cr(NO3)3∙9H2O and Dy(NO3)3·6H2O, with ortho-toluic acid in acetonitrile at ambient temperature. The solvent was then removed and re-dissolved in MeOH/iPrOH (Supplementary Methods). Suitable single crystals (pale purple color) for X-ray analysis were isolated after allowing the solvent to slowly evaporate.

Single-crystal X-ray analysis reveals that compound 1 crystallizes in the triclinic space group, P-1 (see Supplementary Table 1 for full crystallographic details). The asymmetric unit contains half the complex, (three DyIII ions and one half of the CrIII ion) which lies upon an inversion center. Compound 1 is a heterometallic heptanuclear complex consisting of a single CrIII ion and six DyIII ions (Fig. 1). The low CrIII to DyIII ratio of 1:6 in 1 is likely a consequence of the limited solubility of Cr(NO3)3∙9H2O in MeCN. The metallic core is based on two triangular DyIII units that lie above and below a single central CrIII ion, revealing two vertex sharing trigonal pyramids or tetrahedra. The metallic core is stabilized by 8 μ3-hydroxide, 12 ortho-toluate, with MeOH and [NO3]− ligands. Six of the μ3 hydroxide ligands bridge two DyIII ions to the central CrIII ion, while the remaining two bridge the three DyIII ions that make up each triangle. Six of the ortho-toluate ligands each bridge a DyIII–DyIII triangular edge, while six are found to chelate, each to a single DyIII ion. Terminal MeOH ligands coordinate to all six DyIII ions. It is found, however, that, at two of the DyIII sites, disordered MeOH and nitrate ions (50:50 occupancy) are present. The CrIII ion is six coordinate with an octahedral geometry, while the six DyIII ions are eight coordinate.

We have examined the structural distortions at individual DyIII sites using SHAPE software37, 38. The geometry of each DyIII ion is best described by a triangular dodecahedron. The deviation of 2.7 for Dy1 and Dy1′, 1.2 for Dy2 and Dy2′, 1.5 for Dy3 and Dy3′ are observed with respect to the ideal triangular dodecahedron. Selected bond lengths and DyIII–O–DyIII and CrIII–O–DyIII bond angles are given in Supplementary Table 2. The DyIII–O bond lengths are in the range, 2.391–2.492 Å. In the {DyIII3} triangular unit, the bond distance between Dy1–Dy2, Dy1–Dy3, and Dy2–Dy3 is found to be 3.749, 3.767, and 3.780 Å, respectively and the Dy3–Dy1–Dy2, Dy1–Dy2–Dy3, and Dy2–Dy3–Dy1 bond angles are 55.99°, 60.52°, and 59.49°, respectively. The average Dy-(μ3-OH)-Dy angle is 106.0o, while the average Dy-(μ3-OH)-Dy angle also bridging to the CrIII ion is slightly smaller at 103.4o. The centroid to centroid distance between the two triangular units is found to be 5.38 Å. A shorter distance, compared to other linked {DyIII3} triangles11 (5.64 Å) is likely to yield stronger dipolar coupling between the two {DyIII3} units. Packing diagrams of 1, viewed along the a, b, and c axes and of a neighboring pair are shown in Supplementary Fig. 1. There are H-bonds linking adjacent {CrIIIDyIII6} complexes via MeOH···MeOH(solv)···O(carb) groups, combined with edge to face π–π interactions.

Dc (direct current) magnetic susceptibility data were collected for 1 and the variation of χMT with temperature is shown in Fig. 2a. The room temperature χMT product of 87.16 cm3 K mol−1 is in agreement with the value expected (86.9 cm3 K mol−1) for one CrIII (S = 3/2, g = 2, C = 1.875 cm3 K mol−1) and six DyIII (S = 5/2, L = 5,6H15/2, g = 4/3, C = 14.17 cm3 K mol−1) ions that are non-interacting. As the temperature is reduced the χMT product (H = 1 T) decreases gradually between room temperature and 50 K, before a more rapid decrease below this temperature, reaching a value of 14.37 cm3 K mol−1 at 1.8 K. The decrease in χMT at higher temperatures is attributed to the depopulation of the excited mJ Stark states of the DyIII ions, while the more rapid decrease at low temperatures is indicative of the presence of dominant antiferromagnetic exchange interactions. The low temperature χMT value of 14.37 cm3 K mol−1 at 2 K is higher than that expected for a single paramagnetic Cr(III) ion suggest that there are several close lying excited states including that of DyIII ion, which possess significant magnetic moment.

The isothermal M vs. H plots (Fig. 2b; Supplementary Fig. 2) at low fields reveal a non-linear, S-shaped curve at 2 and 3 K, which suggests the presence of a toroidal moment8,9,10,11,12,13,14 and/or possible blockage of the magnetization vector and therefore slow magnetic relaxation. Above 4 K, the plots display a rapid increase in magnetization up to ~2 T followed by a steady increase and almost saturating at 5 T (Fig. 2b). The simulation of the plots in Fig. 2 are discussed later.

To probe for SMM behavior, the magnetization dynamics were investigated via alternating current (ac) susceptibility measurements as a function of both temperature and frequency. A 3.5 Oe ac field was employed, utilizing both a 0 and a 2000 Oe static dc field. A non-zero out-of-phase magnetic susceptibility component (χ″) is observed, at Hdc = 0 Oe, however, no maxima are found upon reducing the temperature down to 1.8 K (Supplementary Fig. 3). This is also the case for when Hdc = 2000 Oe (Supplementary Fig. 3). This does not prove SMM behavior, but suggests the possibility of such, with a small energy barrier to magnetic reorientation and fast relaxation times, even at 1.8 K. As a consequence of the fast magnetic relaxation times, even at temperatures below 2 K, it is suggested that the low-field magnetization behavior points to the presence of a toroidal magnetic moment.

Single crystals of 1 were studied using the micro-SQUID technique at various temperatures and sweep rates for two different orientations of the molecule39. The curve displays a stepped shape of the magnetization (Fig. 3; Supplementary Fig. 4), similar to that observed by the archetypal triangular toroidal system8. The position of the step (Hs) depends on the orientation of the magnetic field with respect to the plane of the triangles, as expected on the basis of the in-plane anisotropy predicted by our model (vide infra). Thus, when the field is parallel to two inversion-related tangential Dy magnetic axes, and perpendicular to none, this leads to an in-plane easy axis, while a magnetic field perpendicular to that direction, thus perpendicular to two Dy easy axes and parallel to none, leads to an in-plane hard axis. In particular, when the field is applied along any Dy–Dy bond vector as in Fig. 3a, i.e., along the y-axis in Fig. 6, hysteresis is observed below 0.8 K, with the coercive field widening on cooling (Hc = 0.6 T at 0.03 K, with a sweep rate of 0.28 T s−1). This behavior is characteristic of an SMM, with slow zero-field relaxation. We then observe a large step in the magnetization at about Hs = 0.7 T. Application of the magnetic field perpendicular to the Dy–Dy bond vector results in a reduction of coercivity (Hc = 0.5 T at 0.03 K, with a sweep rate of 0.28 T s−1, Supplementary Fig. 4). Upon comparison of the hysteresis profile of complex 1 with the archetypal {Dy3} toroidal complex8, we find that profile for 1 appears to be superior comparing the coercitivity. This suggests that coupling between the two {DyIII3} triangles in 1 enhances the zero-field slow relaxation properties of the system. A theoretical model explaining this behavior is developed later in the paper.

Fig. 3

Single-crystal studies. Single-crystal magnetization (M) vs. applied field measurements (μ-SQUID) for complex 1 at a 0.03–0.8 K with the scan rate of 0.14 T s−1; b with different field sweep rates at 0.03 K

We have also recorded EPR spectra at the X-band frequency at 5, 10, and 20 K (Fig. 4). The EPR spectrum recorded at 5 K reveals distinct features at very large g-values (g ~ 14.2). When we increased the temperature we found that the intensity of this signal decreases. There are also weak features at g ~ 2.05, g ~ 1.2, and g ~ 1.03 and the intensities of these features also decrease upon increasing the temperature. To gain an understanding on the nature of these EPR spectra, we have simulated the EPR spectrum by means of the XSOPHE simulation suite40, 41, using a {CrIIIDyIII3} model employing a pseudo S = 1/2 state for each DyIII ion and a S = 3/2 state for the CrIII ion. Ab initio computed g-anisotropies, directions, and J values are given as inputs (see below for details) along with the Dy···Dy and Dy···Cr distances from the X-ray structure. A small perturbation to the Euler angles without altering any other parameters yield reasonable fit to the experimental spectrum recorded at 5 K (Fig. 4; The further details of simulation are given in Supplementary Note 1), offering confidence on the estimated parameters. However, the lines appearing at g ~ 1.23 and g ~ 1.03 are much broader than they appear in the simulation and this may be attributed to a fact that only {CrIIIDyIII3} has been employed in the simulation and not the full {CrIIIDyIII6} Hamiltonian. Multi-frequency EPR including HF-EPR spectra are required, in future, to independently obtain the spin Hamiltonian parameters42.

Fig. 4

EPR spectroscopy. Powder EPR spectra of 1 at X-band frequency at 5, 10, and 20 K and the simulated curve

Theoretical analysis and characterization of a ferrotoroidic ground state

To explain the experimental observations, we performed ab initio calculations of electronic structure and magnetic properties, using MOLCAS 7.843, on individual DyIII and CrIII centers. The computed orientation of the anisotropy axes is shown in Fig. 5. In particular, we employed the ab initio MJ decomposition of the single-ion thermally isolated ground Kramers doublet (KD) wavefunctions along the ab initio g-tensor principal axis, to set up a model Hamiltonian for intramolecular magnetic coupling including dipolar coupling between all pairs of ions, which is parameters-free, intra-ring DyIII–DyIII superexchange interactions parameterized by a single coupling constant J2, and DyIII–CrIII superexchange interactions parameterized by a single coupling constant J1. The coupling parameters J1 and J2 were evaluated via DFT calculations. The well-known dipolar Hamiltonian reads:

where Sq (SCr) are the true spin moments of the Dy (Cr) ions, with the primed and unprimed subscripts labeling Dy ions belonging to different triangles. We note here that when the simple isotropic exchange Hamiltonian is projected on the thermally isolated ground KDs, it becomes a strongly anisotropic non-collinear Ising Hamiltonian, a widely employed protocol for {3d–4f} systems.

Fig. 5

Orientation of the magnetic anisotropy axes. The directions of the local anisotropy axes in the ground Kramers doublet on each Dy site (dotted lines) in 1. Blue arrows are the local magnetic moment in the ground exchange doublet. Black arrows show the con rotation of the toroidal magnetic moment and the yellow arrow is the S6 symmetry axis

To estimate the low-energy wavefunctions and magnetic anisotropy for each of the seven ions in 1, we have undertaken CASSCF+RASSI-SO calculations on the individual DyIII centers43. The calculations yielded the following g-tensor principal values: (Dy1; gx = 0.0523, gy = 0.0927, and gz = 19.5707); (Dy2; gx = 0.0737, gy = 0.0979, and gz = 19.4723); (Dy3; gx = 0.0233, gy = 0.0361, and gz = 19.6059) and (Cr; gx = gy = gz = 2.002) (Supplementary Tables 3, 6). The symmetry-related Dy1′, Dy2′, and Dy3′ ions possess essentially the same g-tensor. Although the two triangles are equivalent, Dy1 and Dy1′ slightly differs due to coordination of methanol in Dy1 and nitrate in Dy1′. These data, along with the J value are found to yield a good fit to the experimental susceptibility data (Fig. 2, vide supra). A qualitative mechanism developed based on single-ion Dy(III) is discussed in detail in the ESI (Supplementary Fig. 6; Supplementary Note 2).

The core structure of the molecule has a pseudo S6 axis passing through the CrIII ion and the center of both of the {DyIII3} triangular units (Fig. 5). The local principal anisotropy axes are found from the calculations to lie in the {DyIII3} plane with an out-of-plane angle of 0.29, 4.5, and 4.7° for Dy1, Dy2, and Dy3, respectively. The DyIII magnetic axes are also found to be almost perfectly aligned with the tangents to an ideal circumference enclosing the triangles (the angle of the anisotropy axis with these tangential directions are in the range of 1.1–7.9°). The computed energies of the eight low-lying KDs reflect that there are three types of DyIII ions in the complex (Supplementary Tables 4, 5), although the ground states of all DyIII ions consist of almost pure atomic | ± 15/2 > KDs. The energy gap between the ground and the first excited state KDs are found to be 142.8, 121.9, and 152.7 cm−1 for Dy1, Dy2, and Dy3, respectively. For the CrIII ion calculations yield isotropic g-tensors (Supplementary Table 3) and an axial zero-field splitting of ~1.0 × 10−4 cm−1.

Due to the large value of the ground KDs, angular momenta projections in \(\left| {{\bf m},{M_{{\rm{Cr}}}}} \right\rangle \), Hamiltonians Eqs. (1) and (2) are both diagonal on such basis, and \(\left| {{\bf m},{M_{{\rm{Cr}}}}} \right\rangle \) represent the low-energy exchange-coupled states of 1. The corresponding energies can be written as a sum of a superexchange contribution, an intra-ring dipolar contribution, an inter-ring dipolar contribution, and DyIII–CrIII dipolar contribution.

For the exchange energies we get (sum over q is understood modulus 3, so that m3+1 = m1):

The energy spectrum resulting from magnetic coupling can now be evaluated summing up Eqs. (3)–(7), and is reported in Fig. 7.

Fig. 7

Energy spectrum describes magnetic relaxation and toroidal states. The energy spectrum of {Dy6Cr} calculated using the energy expressions Eqs. (3)–(7). Energies are in cm−1. The number of dashes for each degenerate energy level indicates the number of degenerate states associated to that level. Besides a few energy levels, it provides a graphical representation of one of the corresponding non-collinear Ising quantum states \(\left| {{\bf m},{M_{{\rm{Cr}}}}} \right\rangle \), where the red (blue) thick arrows at the DyIII sites indicate a +15/2 or −15/2 quantum number along the tangent direction in the lower (upper) triangular ring, a yellow thick arrow at the center indicates the spin projection of the CrIII ion relative to that of the DyIII ions, and a resulting magnetic moment in the upper (lower) ring is indicated with thin blue (red) arrow. Gap between the ground FT state \(\left| { \pm {\tau _1}, \pm {\tau _2},{M_{{\rm{Cr}}}}} \right\rangle \), and the first (AFT) \(\left| { \pm {\tau _1}, \mp {\tau _2},{M_{{\rm{Cr}}}}} \right\rangle \) excited state are indicated

The dipolar coupling energies, Eqs. (5)–(7), contain no free parameter, as they only depend on the specific set of quantum numbers \(\left| {{\bf m},{M_{{\rm{Cr}}}}} \right\rangle \), on the experimental geometrical parameters r and h (Fig. 6), and on the average ground state ab initio magnetic moment given by 10ημB, with η = 0.975. Interestingly, from Eq. (5) it follows that intra-ring dipolar interactions always favor a toroidal texture on each triangle, penalizing the formation of a magnetic moment by an energy gap of ~4 cm−1.

While inter-ring dipolar coupling is smaller than intra-ring coupling, due to larger separation between DyIII ions, from Eq. (6) we learn that in fact this interaction splits FT and AFT states. In particular, the first term in Eq. (6) stabilizes a FT (AFT) ground state if \(h < r\sqrt {7/2} \) (\(h >r\sqrt {7/2} \)), while the second term (describing interactions between inversion-related centers) always favors FT coupling, and is stronger than the first term. In {CrIIIDyIII6}, \(h >r\sqrt {7/2} \), hence the two terms are in competition, the first (second) favoring AFT (FT) coupling. Since the second term is larger, we find here that dipolar interactions stabilize a FT ground state in {CrIIIDyIII6} (Fig. 5).

Moreover, from Eq. (6), we note that a structural design aimed at reducing the distance h between the two triangles will enhance FT coupling. We estimate the dipolar-induced FT/AFT splitting to be ~0.28 cm−1. Note that, in our symmetric model, Dy–Cr dipolar interactions within FT and AFT states (Eq. (7)) are exactly zero (since \({{\cal M}_{\bf m}} = 0\) in FT and AFT states), leading to eight-fold degenerate FT and AFT manifolds. Thus in FT and AFT states, the spin of Cr(III) is freely fluctuating, and FT/AFT toroidal excitations, involving the collective flipping of three DyIII spins, are fully determined by dipolar coupling. The energy gap of 3.3 cm−1 reported in Fig. 7 corresponds instead to the lowest magnetic excitation, obtained upon flipping a single DyIII spin. In the excited states, the Cr(III) spin is blocked in the direction of the in-plane Dy6 magnetic moment.

Considering now superexchange interactions in Eq. (3), we note that the dipolar-induced FT ground state will survive provided that the intra-ring Dy–Dy coupling is antiferromagnetic (i.e., J2 < 0), or ferromagnetic but smaller than dipolar coupling. To estimate the two exchange coupling constants J1 and J2 appearing in Eq. (3), we have employed DFT calculations, replacing the DyIII with GdIII ions in the X-ray structure. The computed coupling constants for the CrIII–GdIII pairs were then rescaled by the ratio between the spin of DyIII (S = 5/2) and that of GdIII (S = 7/2), while GdIII–GdIII were rescaled by the ratio between the square of S = 5/2 and the square of S = 7/233, 34 We obtained J1 = –0.08 cm−1 (CrIII–DyIII coupling) and J2 = –0.043 cm−1 (intra-ring DyIII–DyIII coupling), indicating antiferromagnetic coupling. The estimated antiferromagnetic interaction within the {DyIII3} triangles reinforces the effect of the intra-ring dipolar coupling, leading to a toroidal moment on each isolated triangle.

While a degenerate FT quantum ground state is compatible with the inversion symmetry of the molecule, such symmetry is not compatible with a ferrotoroidically ordered phase. Thus upon FT ordering, a concomitant structural phase transition should occur to rid the crystal of the inversion center, which, in turn, would allow the appearance of linear magneto–electric coupling.

We now use the developed theoretical model to simulate the experimental molar susceptibility, and magnetization. We report the results in Figs. 2 and 8. While the magnetization at 2 K and for fields up to 5 T is expected to be dominated by the low-energy states described by our model (Eqs. (1) and (2)), the molar susceptibility will be dominated by these states only at low temperatures, while at temperatures much higher than the coupling strength it will be dominated by the single-ion response, including the population of excited KDs. To reproduce the correct temperature dependence of χMT within our non-empirical model, we, therefore, evaluated the molar susceptibility as χT = (χT)LT + (χT)HT − (χ0T), where (χT)LT is computed from our model Eqs. (1) and (2), (χT)HT is computed as a sum of the ab initio susceptibilities of the six DyIII ions and the central CrIII ion, while (χ0T) = 0.125 × (6M(Dy)2 + 4 × 3/2 (3/2 + 1)) is the uncoupled contribution to the Curie molar susceptibility (cm3 K mol−1) arising from the six DyIII ground KD’s using the ab initio magnetic moment M(Dy) = 9.75μB and the ground quartet on CrIII. The results of the calculations of χT compared with the experimental data are reported in Fig. 2a (orange line), which shows an excellent agreement between theory and experiment.

Fig. 8

Interpretation of single-crystal magnetization experiments. a Single-crystal magnetization experiment (blue curve) measured at T = 0.03 K and a sweep rate of 0.1 T s−1 for a magnetic field oriented along the in-plane easy axis (y-axis in Fig. 6), superimposed on the theoretical equilibrium magnetization (orange curve) computed for the same temperature and field direction; b Energy levels as function of magnetic field (Zeeman spectrum) computed using the model presented in the text. The magenta (green) arrows connecting a and b associate the steep (smooth) raise of the theoretical (experimental) magnetization with slow/quasi-forbidden (fast/allowed) transition mechanisms occurring at level crossing between the two ground states (involving excited states) encircled in magenta (green). c Magnification of the low-field region of the Zeeman spectrum, highlighting degeneracy points (level crossings) of states between which faster 1-flip and 2-flip tunneling transitions are allowed (green circles and hexagons), and also highlighting 3-flip or higher-flip processes that are essentially forbidden (magenta circles). Some of the allowed phonon emission processes are also indicated with dashed green arrows; d High-field region of the Zeeman spectrum, where the four onion states antiferromagnetically coupled to Cr spin states become isolated from excited states, leading the system’s magnetization to saturate

Note also that despite the absence of fitting parameters in our model, the theoretically calculated magnetization curve reported in Fig. 2b and Supplementary Fig. 7 (orange curve) reproduces very well the experimentally measured M vs. H curve at 2 K and higher temperatures (blue data points). In particular, the predicted FT and AFT low-energy states display a strongly reduced magnetic response, as is evident from the S-shape behavior of the M vs. H curve. A uniform magnetic field does not interact with a toroidal moment, so that only the CrIII ion responds to the field for low H. As the Zeeman energy increases, at H ~ 0.4–0.5 T, a level crossing occurs (vide infra), corresponding to a magnetic excited state becoming the ground state, explaining the S-shape of the M vs. H curve. Finally we present a detailed comparison in the Supplementary Note 3 and 4 where various possible models are discussed and our findings are compared with previous studies of coupled molecular {DyIII3} toroids (Supplementary Fig. 9)10, 11, 18.

Theoretical analysis of the zero-field hysteretic spin dynamics

As a first attempt to interpret the single-crystal magnetization reported in Fig. 8a (blue curve), measured at T = 0.03 K and a field sweep rate of 0.1 T s−1, we calculated the equilibrium magnetization at the same temperature and magnetic field orientation (along the easy axis, y-direction in Fig. 6), also reported in Fig. 8a (orange curve). Experimental and theoretical curves share a few common features, both displaying a rise of the magnetization for intermediate fields, separating two magnetization plateaus corresponding to low- and high-field regions. Theoretical and experimental plateaus are seen to coincide, with the theoretical low-field plateau corresponding to saturation of the free fluctuating Cr spin (MSAT ~ 3µB) within the FT ground state. According to the Zeeman spectrum in Fig. 8b, c, the steep magnetization step observed in the theoretical magnetization at the level crossing field BLC ~ 0.4 T can be identified with the switching of the ground state from the low-field weakly magnetic FT Zeeman state \(\left| { \pm \tau , \pm \tau } \right\rangle \equiv \left| { \pm 1 \pm 1 \pm 1 \pm 1 \pm 1 \pm 1, - 3/2} \right\rangle \) (Supplementary Fig. 8a), to the high-field onion magnetic state \(\left| { + m, + m} \right\rangle \equiv \left| { + 1 + 1 - 1 - 1 + 1 + 1, + 3/2} \right\rangle \), in which half of the Dy spins circulate clockwise, the other half anticlockwise (Supplementary Fig. 8d), adding up to a magnetic moment M ~ 40µB polarized along the field.

As expected, the calculated equilibrium magnetization cannot reproduce inherently dynamical features observed in the experiments, as in the calculations transitions between energy levels are instantaneous on the timescale of the experiment. There are, in particular, two important differences between theory and experiment displayed in Fig. 8a: in the experiment, a hysteresis loop opens up around the zero-field region between –BLC and BLC, and the experimental magnetization step separating low- and high-field plateaus is not abrupt, but occurs within ~0.5 T field range starting at BLC.

Thus, to further analyze the measured magnetization dynamics, we explicitly consider finite-rate transitions between Zeeman states. We are not interested in the microscopic details of such transitions, other than these are driven by terms in the Hamiltonian that have so far been neglected (either because small, or because describing coupling with the surroundings), and other than broadly separating such processes in tunneling transitions characterized by coupling constants Γi, active between degenerate energy levels, and phonon-induced transitions characterized by coupling constants Γi, from higher to lower-energy states, with a probability that grows as the cube of the energy gap. The coupling constants γi and Γi are related to the matrix elements of the relevant Hamiltonian between initial and final states.

Crucially, we propose that a hierarchy should exist in the magnitude of γi and Γi, based on the number of Dy3+ spins that need to be flipped to change initial to final Zeeman states. We only consider in our model the 56 (out of 256) most relevant low-energy states involved in the relaxation dynamics, which are represented in Supplementary Fig. 8. Particularly relevant are the low-field ground FT states\(\left| { \pm \tau , \pm \tau } \right\rangle \), the high-field ground onion state \(\left| { + m, + m} \right\rangle \), and the excited intermediate magnetic states \(\left| { \pm \tau , + m} \right\rangle \) and \(\left| { + m, \pm \tau } \right\rangle \) (Supplementary Fig. 8c). Let us label the coupling constants γi (Γi) as: γCr (ΓCr) for transitions that only flip the Cr3+spin; γ1 (Γ1) for transitions flipping only one Dy3+spin; γ2 (Γ2) for transitions involving two simultaneous spin-flipping events; γ3 (Γ3) for transitions involving three or more spin-flipping events. Thus the hierarchy we invoke here reads: γCr≫γ1 > γ2 > γ3 ~ 0, and ΓCr≫Γ1 > Γ2 > Γ3 ~ 0.

Moreover, to quantitatively investigate the origin of the observed zero-field hysteresis loop, we implement these ideas in a dynamical model based on generalized Pauli master equations describing the dissipative dynamics of the non-equilibrium thermal populations of the CrDy6 states, coupled to an equilibrium reservoir of acoustic phonons at T = 0.03 K, and a source of random stray fields inducing incoherent tunneling between resonant energy levels. The states are obtained from our model Hamiltonians Eqs. (1) and (2), also including a time-dependent Zeeman term describing the interaction with the sweeping magnetic field-oriented along the easy axis, at sweeping rate 0.1 T s−1, varying between B = ± 1.4 T (a triangle wave signal is used). As outlined in the methods section, the relevant dissipative equations of motion for diagonal elements σii (populations) of the reduced density matrix σ, in the adiabatic approximation, are:

We numerically solved Eq. (8) for the time-dependent populations σ(t) of the 56 CrDy6 states considered here (Supplementary Figs. 8, 10). The magnetization curve \(M\left( t \right) = {\rm{Tr}}\left[ {\sigma (t)M} \right]\) is then computed (M is the magnetic moment operator along the field), parametrically plotted vs. the sweeping field, and reported in Fig. 9 with the experimental magnetization.

Fig. 9

Simulation of single-crystal magnetization. Single-crystal experimental magnetization (blue curve) measured at T = 0.03 K and a sweep rate of 0.1 T s−1 for a magnetic field oriented parallel to the triangle’s planes and along the easy axis (y-axis in Fig. 6), superimposed on the simulated dynamical magnetization at the same temperature, sweep rate, and field orientation, by solving Eq. (8) on the basis of 56 out of the 256 states obtained from our model and reported in Supplementary Fig. 8, for the following numerical values of the transition rates appearing in the equations: ΓCr = 105 Hz/(cm−1)3≫Γ1 = 10−7 × ΓCr > Γ2 = 10−3 × Γ1, and γCr = 1016 Hz2≫γ1 = 1012 Hz2 > γ2 = 10−3 × γ1

Quite remarkably, the opening of the hysteresis in the field region between +BLC and –BLC is captured by our model, together with the narrowing and closing of the hysteresis loop at fields ~BLC. Also, the simulated dynamical magnetization now displays a smooth increase between the low- and high-field “plateaus”, despite the low temperature, on account of the cascade of indirect transitions between the low-field (FT) and high-field (onion) ground states mediated by intermediate excited magnetic states.

There are, of course, shortcomings in our simulation. For instance, the zero-field hysteresis loop, after closing down as in the experiment, opens up again at larger fields, a feature that is only marginally present in the experiment. Also, the theoretical magnetization step covers a smaller field range than the experimental one. While we cannot exclude that a more thorough exploration of the parameter space might improve the fitting, we expect that these shortcomings will be partially overcome if all 256 states are included in the simulation, as the additional excited states would participate to the multi-step magnetization relaxation, further widening the field range covered by the magnetization step. Further discussion of our dynamical model, including a partitioning of the contributions to M(t) from individual states, is reported in the Supplementary Note 5.

In summary, a heptanuclear {CrIIIDyIII6} complex has been synthesized and structurally characterized. Experimental evidence in conjunction with theoretical calculations reveal and explain the presence of both SMM and single-molecule toroidic behavior. The toroidal states in the individual {DyIII3} triangles are found to be ferrotoroidically coupled. We note here that the CrIII ion does not play any fundamental role in the predicted FT coupling, and the coupling between the two toroidal wheels can in fact be fully explained solely in terms of dipolar interactions, which depend solely on the structural parameters of the complex. The fundamental structural elements influencing the strength of dipolar FT coupling are the staggered arrangement of the two triangles with respect to each other, and the distance h between the two wheels, which we found not to be optimal for maximizing FT coupling (i.e., \(h >r\sqrt {7/2} \)). The design and synthesis of similar {MIIIDyIII6} complexes, maintaining the staggered arrangement of the two triangular units, but featuring a smaller, even diamagnetic ion, is expected to lead to a smaller inter-ring distance h, or even to \(h < r\sqrt {7/2} \), thus according to Eq. (6), to a stronger FT coupling. This route is currently being explored in our labs. Importantly, our findings indicate, for the first time, how coupling between toroidal moments can be manipulated by structural design. Finally, for the first time, the experimental single-crystal magnetization dynamics of a polynuclear Dy complex, displaying zero-field opening of a hysteresis loop, is simulated via a theoretical dynamical model, showing that the FT ground state plays a pivotal role in hindering the flipping of magnetic onion states in a sweeping field, thus slowing down the zero-field magnetic relaxation.

Methods

Synthesis of [DyIII6CrIII(OH)8(ortho-tol)12(MeOH)5(NO3)]∙3MeOH (1)

Cr(NO3)3·9H2O (0.4 g, 1 mmol) and Dy(NO3)3·6H2O (0.22 g, 0.5 mmol) were dissolved in MeCN (20 ml), followed by the addition ortho-toluic acid (0.14 g, 1.0 mmol) and triethylamine (0.55 ml, 4.0 mmol), which resulted in a pale purple solution. This solution was stirred for 4 h after which time the solvent was removed to give a purple oil. The oil was re-dissolved in MeOH/iPrOH (1:1) and the solution allowed to slowly evaporate. Within 10–15 days, pale purple crystals of 1 had appeared, in approximate yield of 25% (crystalline product). Microanalysis for CrDy6C104H124NO43: expected (found); C 40.25 (39.86), H 4.02 (3.86), N 0.45 (0.62).

The synthesis reaction was carried out under aerobic conditions. Chemicals and solvents were obtained from commercial sources and used without further purification.

X-ray crystallography

X-ray measurements for 1 were performed at 100(2) K at the Australian synchrotron MX1 beam line45. The data collection and integration were performed within Blu-Ice46 and XDS47 software programs. Compound 1 was solved by direct methods (SHELXS-97)48, and refined (SHELXL-97)49 by full least matrix least squares on all F2 data within X-Seed50 and OLEX-2 GUIs51. Crystallographic data and refinement parameters are summarized in Supplementary Table 1.

Magnetic measurements

The magnetic susceptibility measurements were carried out on a Quantum Design SQUID magnetometer MPMS-XL 7 operating between 1.8 and 300 K for DC-applied fields ranging from 0 to 5 T. Microcrystalline samples were dispersed in vaseline in order to avoid torquing of the crystallites. The sample mulls were contained in a calibrated gelatine capsule held at the center of a drinking straw that was fixed at the end of the sample rod. Ac susceptibilities were carried out under an oscillating ac field of 3.5 Oe and frequencies ranging from 0.1 to 1500 Hz.

EPR instrumentation

The X-band measurements were made on Bruker spectrometers, at IIT Bombay, with a helium gas-flow cryostat. The X-band measurements were carried out at 5, 10, and 20 K.

Computational details

Even though there are only three crystallographic nonequivalent DyIII centers, we performed the calculations on all six DyIII ions to determine the direction of the local anisotropy axis. Using MOLCAS 7.843, ab initio calculations were performed on the six DyIII ions using the crystal structure of 1. The structure of the modeled Dy fragment employed for calculation is shown in Supplementary Fig. 5, where the neighboring DyIII ions are replaced with diamagnetic LuIII ions and the CrIII ion replaced with ScIII ion. We have employed this methodology to study a number of DyIII/ErIII SMMs52,53,54,55. The relativistic effects are taken into account based on the Douglas−Kroll Hamiltonian56. The spin-free eigenstates are achieved by the Complete Active Space Self-Consistent Field (CASSCF) method57. We have employed the [ANO-RCC…8s7p5d3f2g1h.] basis set58 for the Dy atoms, the [ANO-RCC…3s2p.] basis set for the C atoms, the [ANO-RCC…2 s.] basis set for H atoms, the [ANO-RCC…3s2p1d.] basis set for the N atoms, the [ANO-RCC…4s3p1d.] basis set for the Sc atom, the [ANO-RCC…5s4p2d1f.] basis set for the La atom, and the [ANO-RCC…3s2p1d.] basis set for the O atoms. First we performed the CASSCF calculation including nine electrons across seven 4f orbitals of the Dy3+ ion. With this active space, we have computed 21 roots in the configuration interaction procedure. After computing these excited states, we have mixed all roots using RASSI-SO59 and spin–orbit coupling is considered within the space of calculated spin-free eigenstates. Moreover, we have considered these computed SO states into the SINGLE_ANISO60 program to compute the g-tensors. The DyIII ion has eight low-lying KDs for which the anisotropic g-tensors have been computed. The Cholesky decomposition for two electron integrals is employed throughout our calculations. We have extracted the crystal field parameters using the SINGLE_ANISO code as implemented in MOLCAS 7.8.

DFT calculations were performed using the B3LYP functional61 with the Gaussian 09 suite of programs62. To estimate the exchange constant between CrIII–DyIII and DyIII–DyIII ions, the dysprosium ions were replaced with the spin-only GdIII ions in order to investigate the exchange interaction between the DyIII ions, which was then rescaled to the spin of dysprosium ions. We have used the LanL2DZ ECP basis set for Cr63, 64, the double-zeta quality basis set employing Cundari–Stevens (CS) relativistic effective core potential on Gd atom65, 6–31G* basis set for the rest of the atoms. The DFT calculations combined with the Broken Symmetry (BS) approach66 have been employed to compute the magnetic exchange J value. The actual energy spectrum and wavefunction and magnetic texture calculations, together with simulation of the magnetic measurements, was carried out by inserting the ab initio data into a model intramolecular magnetic coupling Hamiltonian, the development of which is described in the theoretical analysis section (see above).

Method of simulation of the magnetization dynamics

As is well known67, in the derivation of the Pauli equations for the diagonal reduced density matrix σ from the Liouville–Von Neumann equations for the total density matrix describing the coupled system–phonon reservoir, it is expedient to switch to the interaction picture, so that the unperturbed quantum system and reservoir Hamiltonians do not explicitly appear in the transformed equations of motion. The equations of motion in the interaction picture are thus propagated only by the spin–phonon coupling Hamiltonian, which in second order and in the Born–Markov limit generates the dissipative relaxation dynamics of the reduced density matrix. In our particular case, where time dependence arises both from the random dissipative relaxation fields, and from the periodic time dependence of the Zeeman Hamiltonian, a useful interaction picture can still be achieved by including the time-dependent Zeeman Hamiltonian in the zeroth order Hamiltonian \({H_0}(t) \equiv {H_0}\) together with Hamiltonian Eqs. (1) and (2). To avoid the complications related to the resulting time dependence of H0, which would imply a zeroth order time evolution operator only defined exactly via time-ordered products, we assume valid adiabatic approximation, so that at any time the system is assumed to be described by the eigenfunctions \(\left| {i(t)} \right\rangle.\) and eigenvalues Ei(t) of the time-dependent Hamiltonian H0, so that \({H_0}\left( t \right)\left| {i\left( t \right)} \right\rangle = {E_i}\left( t \right)\left| {i\left( t \right)} \right\rangle \), and the time evolution operator can be simplified as an exponential operator that remains diagonal in the basis of the 256 system’s eigenfunctions \(\left| {i(t)} \right\rangle \) at any time. Using well-known approximations to avoid at each given time, the integration of eigenvalues at all previous times, as described for instance in the paper68 and within the secular approximation67, two decoupled sets of equations of motion are obtained for the elements of the reduced density matrix \({\sigma _{ij}} = \left\langle {i{\rm{|}}\sigma {\rm{|}}j} \right\rangle \) between the adiabatic eigenstates of H0, one for the diagonal (populations, i = j) and one for the off-diagonal (coherences, \(i \ne j\)) elements of the CrDy6 reduced density matrix σ, which can be written in the usual Schrodinger’s picture as:

where \({\dot \sigma _{ij}} \equiv {\rm{d}}{\sigma _{ij}}/{\rm{d}}t\), \(\left[ {{H_0},\sigma } \right]\) is the commutator between the operators H0 and σ, δij is a Kronecker delta, \({W_{l \to m}}(t)\) is the transition rate from eigenstate \(\left| {l(t)} \right\rangle \) to eigenstate \(\left| {m(t)} \right\rangle \), which become themselves time dependent because of the time-dependent Zeeman field, both via the Zeeman eigenvalues Ei(t), but in principle also via the eigenstates \(\left| {l(t)} \right\rangle \), as Cr spin states are now coupled by the magnetic field. Finally, we have \({\mu _{ij}}(t) = 1/2\mathop {\sum }\nolimits_k \left( {{W_{i \to k}}(t) + {W_{j \to k}}(t)} \right).\)

The details of the spin–phonon coupling Hamiltonian contributing to the transition rates for this polynuclear system coupled to its crystal phonons represents in principle a formidable problem, whose detailed analysis goes beyond the scope of the current paper. Hence we describe here the coupling of the CrDy6 states to an idealized equilibrium acoustic phonon reservoir in terms of well-known phenomenological transition rates \({W_{l \to m}}(t)\) obtained within the Debye model2 and given by:

where Ei(t) are now the time-dependent Zeeman energies of the CrDy6 quantum states, while \({{{\Gamma }}_{ij}}\) are numerical parameters measuring the transition rate ΓCr, or Γ1, or Γ2, according to how many spin flips connect states \(\left| {i(t)} \right\rangle \) and \(\left| {j(t)} \right\rangle \) (we set \({{{\Gamma }}_{ij}} = 0\) if the two states are connected by three or more spin flips).

Finally, in order to account for the relaxation dynamics associated to quantum tunneling processes induced by, e.g., random stray magnetic fields produced by the fluctuations of nuclear dipole moments or neighboring molecular magnetic moments, we follow the approach proposed by Leuenberger and Loss69 (see in particular the derivation of their Eqs. (35) and (36)), and further correct the zeroth order Hamiltonian H0 in Eq. (9) with an additional term V describing tunneling between the Dy3+ and Cr3+ magnetic states via, e.g., coupling to fluctuating magnetic fields arising from neighboring spins. We note in fact that in Eq. (9), the diagonal matrix elements of the commutator between the eigenstates of H0 are exactly zero in the adiabatic approximation, and thus the differential equations for the populations are there decoupled from those for the coherences. If on the other hand we now introduce the tunneling operator V, so that in Eq. (8) \({H_0}\left( t \right) \to {H_0}\left( t \right) + V\), since V by definition has non-zero off-diagonal matrix elements between the eigenstates of H0, this correction will re-introduce coupling between populations and coherences in Eq. (9). Assuming a steady-state approximation for the coherences, so that \({\dot \sigma _{ij}}(t) \approx 0\) on the timescale over which the populations \({\sigma _{kk}}\left( t \right)\) display appreciable changes, Eq. (9) can be transformed in a system of differential equations for the populations only, a set of generalized Pauli equations that reads (this is Eq. (8), reported here for convenience):

where \({\omega _{ij}}\left( t \right) = \left( {{E_i}\left( t \right) - {E_j}\left( t \right)} \right)/\hbar \). We note that Eq. (12) is in fact the equation of a Lorentzian lineshape with a maximum at the resonance \({\omega _{ij}}\left( t \right) = 0\), thus describing incoherent tunneling processes between Zeeman states in proximity of level crossing. While in principle, this derivation predicts via Eq. (12) that the broadening of the Lorentzian lineshape should also be treated as time-dependent, we choose here to fix the broadening as a constant parameter λ to simplify our dynamical model.

The task of solving Eq. (11) is computationally not trivial, and to make the calculations faster and more stable, instead of including all the 256 states of the CrDy6 system arising from our low-energy model, we decided to include only the lowest energy states over the magnetic field range explored, which beside the 16 FT and AFT states, they correspond to Dy-based magnetic states with anisotropy axes aligned along the sweeping field. These states correspond in fact to the configurations shown in Supplementary Fig. 8. We report in Supplementary Fig. 10 the 56 Zeeman levels entering our dynamical model (Eq. (11)), as function of the magnetic field, which can be compared with the full plot of the 256 Zeeman levels in Fig. 8b–d in the main text.

Note that in principle we have at least seven free parameters entering Eq. (11): the three spin–phonon relaxation rates ΓCr, Γ1, and Γ2, with ΓCr≫Γ1 > Γ2, the three squares of tunneling relaxation rates γCr, γ1, and γ2 with γCr≫γ1 > γ2, and the Lorentzian broadening λ. The actual maximal relaxation rate at the Lorentzian maximum (i.e., at exact level crossing) is in fact given by γk/λ, with γk corresponding to the square of the tunneling splitting at level crossing expressed in Hz. Some of these parameters can in fact be fixed within reasonable ranges. For instance, from ref. 70., we learn that the typical range of values for ΓCr, i.e., the spin–phonon relaxation rate for simple paramagnetic ions, is in the range 3 × 103–3 × 105 Hz/(cm−1)3. We further assume that fluctuating stray fields are of the order of 1 mT, given that the magnitude of the non-fluctuating dipolar field induced at any Cr3+ or Dy3+ site within one molecule varies between ~8 mT (inter-wheel interactions) and ~200 mT (intra-wheel interactions), and that due to the strongly anisotropic and non-collinear character of the local magnetic moments, orientational effects will greatly reduce these fields, which can only induce tunneling when oriented perpendicular to the local anisotropy axes. Considering that the transition magnetic moment between Cr3+ spin states whose MS quantum number differ by one unit is ~2µB, we get the following rough estimation for \({\gamma _{{\rm{Cr}}}}\sim {\left( {2{\mu _B} \times 1{\rm{mT}}} \right)^2}/{\hbar ^2}\sim {10^{16}}{\rm{H}}{{\rm{z}}^2}\) (corresponding to a maximal tunneling frequency in the absence of broadening of ~0.1 GHz). Also, from our single-ion ab initio calculations for the Dy3+ ions in CrDy6, we find that the average value of the transition matrix element between the two components MJ~ ± 15/2 of the ground KD on each ion is of the order of \(\sim {10^{ - 2}}{\mu _B}\). Given that the 1-flip tunneling transitions only involve one Dy3+ ion at the time, we obtain \({\gamma _1}\sim {\left( {{{10}^{ - 2}}{\mu _B} \times 1{\rm{mT}}} \right)^2}/{\hbar ^2}\sim {10^{11}}- {10^{12}}{\rm{H}}{{\rm{z}}^2}\) (i.e., a maximal tunneling frequency of 1 MHz). While these rough considerations reduce the possible values for three of the seven free parameters, we still need to arbitrarily fix four of them: Γ1, Γ2, γ2, and λ.

Given the approximate nature of the model, it is not our aim to attempt a fully satisfying fitting of the remaining four parameters. Thus guided by the relations ΓCr = 105 Hz/(cm−1)3≫Γ1 > Γ2, and γCr = 1016 Hz2≫γ1 = 1012 Hz2 > γ2, we found that a reasonably good agreement between the calculated magnetization \(M\left( t \right) = {\rm{Tr}}\left[ {\sigma (t)M} \right]\) and experimental magnetization can be obtained for Γ1 = 10−7 × ΓCr, Γ2 = 10−3 × Γ1, γ2 = 10−3 × γ1, and λ = 1010 Hz, which are the parameters used to obtain both Fig. 9 and Supplementary Fig. 11.

Data availability

The X-ray crystallographic coordinates for structure reported in this study have been deposited at the Cambridge Crystallographic Data Centre (CCDC), under deposition number 1435033. These data can be obtained free of charge from the Cambridge Crystallographic Data Centre via www.ccdc.cam.ac.uk/data_request/cif.

Acknowledgements

G.R. would like to acknowledge the financial support from DST-India, and IIT Bombay for the high performance computing facility. K.S.M. and G.R. thank the Australia-India AISRF program for support. K.R.V. is thankful to the IITB-Monash Research Academy for a PhD studentship. G.R. acknowledges funding from SERB-DST (EMR/2014/000247) for financial support. A.S. acknowledges support from the Australian Research Council, Discovery Grant ID: DP15010325. W.W. acknowledges the Alexander von Humboldt foundation.

Author information

Affiliations

IITB-Monash Research Academy, IIT Bombay, Mumbai, 400076, India

Kuduva R. Vignesh

School of Chemistry, University of Melbourne, Melbourne, VIC, 3010, Australia

Alessandro Soncini

School of Science and the Environment, Division of Chemistry, Manchester Metropolitan University, Manchester, M15 6BH, UK

Search for Keith S. Murray in:

Search for Gopalan Rajaraman in:

Contributions

A.S., K.S.M., S.K.L., and G.R. visualized and designed the project. K.R.V. and S.K.L. carried out the syntheses and characterized the materials. S.K.L. and K.R.V. performed the synchrotron X-ray scattering measurements, analyzed the data, and solved the crystal structure. W.W. carried out micro-SQUID measurements. K.R.V. and A.S. carried out the ab initio calculations. A.S. developed the theoretical models for magnetic coupling, and for the hysteretic dynamics of the magnetization, and used them for the simulation and interpretation of the magnetic data. A.S., K.S.M., and G.R. wrote the manuscript. All the authors discussed the results and contributed to the manuscript.

Competing interests

The authors declare no competing financial interests.

Corresponding authors

Electronic supplementary material

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.