Synthesis and Magnetic Properties of the Multiferroic [C(NH2)3]Cr(HCOO)3 Metal–Organic Framework: The Role of Spin–Orbit Coupling and Jahn–Teller Distortions

We report for the first time the synthesis of [C(NH2)3]Cr(HCOO)3 stabilizing Cr2+ in formate perovskite, which adopts a polar structure and orders magnetically below 8 K. We discuss in detail the magnetic properties and their coupling to the crystal structure based on first-principles calculations, symmetry, and model Hamiltonian analysis. We establish a general model for the orbital magnetic moment of [C(NH2)3]M(HCOO)3 (M = Cr, Cu) based on perturbation theory, revealing the key role of the Jahn–Teller distortions. We also analyze their spin and orbital textures in k-space, which show unique characteristics.


S1. Additional Sample Preparation Details
Solvents were purchased from Fisher, degassed, and stored under argon before use; Et 2 O was dried by reflux over potassium and stored a potassium mirror.CrCl 3 .6H 2 O and CrCl 2 were purchased from Aldrich, guanidinium carbonate and formic acid were purchased from Acros, sulfuric acid was purchased from Fisher and chromium metal was purchased from Hopkin and Williams Ltd.; all were used as obtained without further purification.
For the synthesis method use to produce single crystals of 1-Cr [C(NH 2 ) 3 ](HCOO) was made by suspending [C(NH 2 ) 3 ] 2 CO 3 (5.00g, 27.8 mmol) in isopropanol (20 cm 3 ).A solution of H 2 COO (2.5 cm 3 , 66 mmol) in isopropanol (10 cm 3 ) was added dropwise with stirring over the course of 20 minutes, resulting in significant effervescence and dissolution of the solid.The reaction mixture was stirred overnight, and the solvent removed in vacuo to give a colorless liquid.This was triturated with 30 cm 3 diethyl ether to give a colorless solid.The liquor was removed by canula filtration and the solid dried in vacuo to give crude guanidinium formate as a colorless solid (5.45 g, 52 mmol, 93% yield) which was used without further purification.
The solution of [C(NH 2 ) 3 ](HCOO) used for preparation of the bulk sample was made by suspending [C(NH 2 ) 3 ](CO 3 ) (1.81 g, 20 mmol) in H 2 O (10 cm 3 ) with H 2 COO (1.88 cm 3 , 50 mmol) then added dropwise.This induced rapid effervescence, which ceased after the reaction was stirred for a further hour, leading to the formation of a colorless solution.This solution was heated to boiling briefly to remove CO 2 and then allowed to cool to ambient temperature.

S2. Single Crystal Diffraction Measurements and Analysis
SCXRD studies utilized Cu-Kα (λ = 1.54184Å) radiation with single crystals mounted on MiTeGen microloops with temperature control via an Oxford Cryosystems Cryostream.
Unit cell determination, data reduction and absorption corrections were performed using CrysAlisPro 171.38.43.S1 The structures were solved with the SHELXT structure solution program S2 via intrinsic phasing controlled via the Olex2 GUI S3 and subsequently refined with the SHELXL refinement package using least squares minimisation S4 (See Table S1 for crystallographic details).Non-hydrogen atoms were refined anisotropically, and hydrogen atoms were included using a riding model with displacement parameters constrained to be 1.2 or 1.5 times that of the N or C to which they are attached.
The structure obtained at 100 K is compared with the theoretically predicted structure used for the DFT calculations.The lattice constants used for the DFT calculations deviate from the experimental lattice constants by −0.19 %, −0.79 %, and −1.56 % for a, b, and c axes, respectively.For the internal coordinates of atoms, the quantitative comparison was performed by using COMPSTRU of Bilbao Crystallographic Server.S5 In terms of the quantities defined therein, the differences are represented by the degree of lattice distortion S = 0.0060, the maximum distance d max.= 0.2705 Å, the arithmetic mean of the distance d av.= 0.1171 Å, and the measure of similarity ∆ = 0.038.From these values, we conclude that the two structures are close enough to justify our DFT description.Table S1: Crystallographic data for 1-Cr determined using single crystal X-ray diffraction between 100-400 K.The single crystal was treated as a racemic twin.
Table S2: Cation coordination bond angles in • at 300 K for 1-Cr.O2-Cr1-O4 88.9(4) O1-Cr1-O4 90.8(4) O3-Cr1-O4 178.7(4)The space group of the 1-Cr crystal is P na2 1 , comprising four crystalline symmetry operations: a twofold rotation 2 z about the z axis followed by a translation of half the c lattice constant, and two mirror operations m x , m y whose planes are perpendicular to the a and b axes respectively.One can define four magnetic order parameters per cell:

S3. Symmetry Analysis of Magnetic Ground State
where the positions of magnetic (Cr) ions are given by the 4a Wyckoff positions of the P na2 1 space group: Starting from the paramagnetic space group P na2 1 1 ′ , the transformation rules for the order parameters introduced before in Eq. ( S1) can be readily obtained and are given in Table S3.
Using these transformation rules, it is straightforward to identify allowed couplings between Table S3: Transformation rules of magnetic order parameters and corresponding magnetic irrep (isotropy subgroups also given in brackets).
different magnetic order parameters, as well as to tailor the possible allowed secondary order parameters generated by non-collinear arrangements of canted magnetic moments.In general, four magnetic space groups can entail the magnetic configuration of the P na2 1 structure, i.e., P na2 1 , P n ′ a ′ 2 1 , P na ′ 2 ′ 1 and P n ′ a2 ′ 1 .As consistent with the previous studies, S6,S7 our DFT(+U + J) calculations confirm that the A-type AFM always displays the lowest energy and thus ruling out the magnetic space group P n ′ a2 ′ 1 .Assuming that the primary order parameter is A, the following three situations are compatible with the symmetry properties: two of which are compatible with weak ferromagnetic components (parallel to a and c axes).
DFT calculations for 1-Cr show that A a corresponds to the most stable magnetic configuration.

S4. Spin Model
Magnetic properties are modeled via an isotropic 3D Heisenberg model for classical spins (with S = 2): Two different (bond-anisotropy) exchange interactions are considered: J c > 0, denoting the AFM exchange coupling along the c axis, and J ab < 0, denoting the FM interaction within ab planes.The isotropic spin model is supplemented by a site-dependent single-ion anisotropy (SIA) defined as S7 Here the local SIA tensor, accounting for the JT-induced orbital order, is given in the local reference frame of distorted CrX 6 octahedra by S6)   and R(θ t,i ) is a rotation matrix that accounts for the effect of tilting rotations on SIA: Assuming an A-AFM configuration aligned along the a-axis, we modeled the secondary ferromagnetic and AFM order parameters.The mean-field energy of the model can be evaluated assuming the following expressions for the local magnetic moments: where ϕ = π/4 to enforce the primary order parameter A a , while δ, ϵ account for the canting producing weak G b , M c components, respectively.The mean-field energy reads: The canting angles can be determined by minimizing Eq. ( S9) with respect to both δ and ϵ.Simple analytical formulas are obtained assuming δ = 0 or ϵ = 0, giving respectively (Note a corrected sign for ϵ to Ref. S7): We notice that the ferromagnetic canting angle scales roughly as 1/J c , while the G-AFM angle δ ∼ 1/J ab .(2.5,0.5) and set 2 (3,1), respectively.In a similar way, the optimal canting angle of the Model parameters can be estimated by total-energy mapping of DFT computations, as done in the Supplementary Information of Ref. S7 (Note that the method therein is formulated as 'per unit cell' values).In addition, constrained spin calculations were performed to obtain the optimal canting angle ϵ as shown in Fig. S6.Two sets of parameters have been used, as listed in Table 1 of the main text, corresponding to DFT+U results with fixed U − J = 2 eV but with (U, J) = (2.5, 0.5) or (U, J) = (3, 1).Since hybrid compounds containing Cr 2+ ions are rare, there is no available reference of the U parameter.Therefore, the choice of these parameters is arbitrary.Note that J c > |J ab | implies that the G-AFM canting dominates over the weak ferromagnetic component because ϵ ∼ 1/J c and δ ∼ 1/J ab .
We briefly comment on the possible role of the Dzyaloshinskii-Moriya interaction (DMI), S8,S9 that in the previous work S7 was shown to be incompatible with {A a ; M c } or {A c ; M a } spin orderings.When we consider two Cr ions neighboring along c-direction with either one of those two spin orderings, the octahedra tilting implies that the DM vector is parallel to a-axis (D ∥ â).On the other hand, the cross-product of two spins is parallel to the b-axis Since the DMI term in the Hamiltonian is usually described in the form of D • (S 1 × S 2 ) = 0, it can not affect the spin canting.However, this is not the case in the presence of the G b or C b orders, which are not considered in the previous work.S7 Let us first consider the {A a ; G b , M c } case (P n ′ a ′ 2 1 ).As a first-order expansion, we can write Hence, we can not strictly rule out the DMI as an origin of the spin canting.Nevertheless, we can ignore the DMI in practice for the 1-Cr where we adopted P n ′ a ′ 2 1 : {A a ; G b , M c } order because the DMI energy 2DS 2 δϵ is a second-order term in the small canting angles.
On the other hand, in the case of {A c ; C b , M a } (P na ′ 2 ′ 1 ), we can write S 1 = S(ϵ, δ, 1) and The DMI is linear to δ, which does not induce a net WFM moment.Therefore, the original argument of dismissing the DMI is still valid for P na ′ 2 ′ 1 , which corresponds to the 1-Cu.

S5. Monte Carlo Results
Monte Carlo simulations have been performed on a L×L×L supercell, with L = 10, 20 or 24, comprising 4000, 32000, or 55296 spins, respectively.A standard Metropolis algorithm has been used, with 10 5 MC steps for equilibration and 10 5 MC steps for averaging.
The transition temperature can be estimated from the peak appearing in the specific heat per unit cell evaluated as where k B is the Boltzmann constant, β = 1/k B T , E the energy and N s the number of spins in the simulation cell, while ⟨...⟩ indicates statistical averages.One can similarly define the (generalized) magnetic susceptibilities as where M is the ferromagnetic order parameter L = A, C, G denotes the antiferromagnetic order parameter (all magnetic order parameters are defined per magnetic ion in MC calculations).A singularity in the generalized magnetic susceptibility allows identifying the primary magnetic order parameter.
Given the magnetic anisotropy, the full magnetic susceptibility tensor has also been evaluated, each component being defined accordingly as: The "powder average" of the magnetic susceptibility is simply defined as χ M av = ( α χ M αα )/3.Results of MC calculations for L = 24 are summarized in Figs.S7 and S8, showing specific heat, primary order parameter A and associated susceptibility as well as secondary order parameters.The estimated critical temperature is ∼ 14.5 K for set 1 and ∼ 9.5 K for set 2. In both parameter sets, the dominant secondary order parameter appears to be G b : the WFM component develops in both cases at the transition temperature, but it is one order of magnitude smaller than the primary AFM-A order parameter when using parameter set 1.In Fig. S9, the powder-averaged magnetic susceptibility is also shown for parameter set 2: despite the curve is still quite noisy, a qualitative trend can already be seen, showing a hump followed by a sizeable decrease as the temperature is lowered across the transition temperature T c .The square of the magnetic moment is µ 2 B S 2 in the classical approach, whereas µ 2 B S(S + 1) in the quantum mechanics.When the χ av is converted to the physical unit (emu mol −1 Oe −1 ), this results in a difference of 1.5 times between classical and quantum cases for S = 2.In Fig. 3 (a) of the main text, the quantum mechanically converted values are used.Note that the previous study S7 which ignored the U and J corrections predicted the critical temperature as about 40 K.

S6. DFT Parameter Dependence of Magnetic Moments
The calculated magnetic moments of the 1-Cr/Cu can depend on the DFT parameters, as shown for the spin magnetic moment of 1-Cr in the previous section.Table S4 shows the net spin and orbital weak ferromagnetic moments per unit cell (4 f.u.) and the magnitude of the orbital magnetic moment of a single Cr/Cu ion in the 1-Cr/Cu.One can see that the orbital magnetic moments are enhanced by the inclusion of the +U + J correction.

S7. Matrix Representations of the Spin-orbit Coupling
To express the SOC Hamiltonian H SOC = ζS•L, let us denote the local coordinate unit vectors for spin operator as (x ′ , y ′ , z ′ ), and those for orbital angular momentum as (x, y, z).H SOC i.e., s x = sin θ cos ϕ, s y = sin θ sin ϕ, and s z = cos θ.We rotate the primed coordinates with respect to the unprimed coordinates.Then the primed coordinate unit vectors are x ′ = cos θ cos ϕx + cos θ sin ϕy − sin θz, y ′ = − sin ϕx + cos ϕy, and z ′ = sin θ cos ϕx + sin θ sin ϕy + cos θz.The SOC Hamiltonian is written as S10 The matrix representation for the orbital angular momentum operator is determined by the quantum mechanical relations for the angular momentum states L z |l, m⟩ = m |l, m⟩ and We take d-orbitals (l = 2) in the real spherical harmonics form as a basis: Having the order of basis in {d yz , d zx , d xy , d x 2 −y 2 , d z 2 }, matrix representations of L i 's are as follows in the atomic units.
The above orbital angular momentum operators can be redefined according to the JT effective Hamiltonian.The JT-rotated e g -orbitals in Eq. ( 2) of the main text define the unitary matrix where U is the unitary matrix for the whole d-orbitals.Consequently, the new orbital angular momentum operator matrices can be obtained by unitary rotation with this matrix, The SOC in the JT distorted system can be expressed by using these operators.
In the following sections, we omit the superscripts 'd' and 'new' of L i .

S8. Calculation of The Perturbed Orbital Angular Momentum
This section describes the detailed steps toward Eq.(3) in the main text.The unperturbed basis with the 'JT rotated' e g -orbitals are {|d 0 nσ ⟩} where n = {yz, zx, xy, −, +} and σ =↑ or ↓ spins.The first order corrected d-orbitals are where the α and β are combined indices of orbital species and spin.The orbital angular momentum expectation value for a perturbed d-orbital d n↑ is The t 2g orbital components are included in this expression.This is why the orbital angular momentum can appear.In the d 4 configuration, Matrix elements ⟨d 0 n |L i |d 0 m ⟩ are given by Eq. ( S20)-(S22).

S9. Details in The Model Analysis for 1-Cr/Cu
In order to apply Eq. ( 3) in the main text to the orbital magnetic moments of 1-Cr/Cu, we take one out of four Cr/Cu ions in a unit cell of 1-Cr/Cu, say Cr/Cu1, as a reference to describe the system (See Fig. S10 (a-c)).If we know the local MO 6 structure and the local moment of Cr/Cu ion at one site, those of other sites are determined by the space group and magnetic group symmetry.The transformation rules of the magnetic moment with respect to Cr/Cu1 by each symmetry operation belonging to P na ′ 2 ′ 1 and P n ′ a ′ 2 1 , where the identity operator 1 corresponds to Cr/Cu1 site, are listed in Table 2 in the main text.Next, rotate the octahedron by −π/4 around the c axis, then by tilting angle −θ t around the a axis consecutively.As a result, (x, y, z) can be written as follows.
The direction of the local spin magnetic moment of Cr/Cu1 is expressed in terms of the spherical coordinates θ spin and ϕ spin with respect to the local coordinate.If we ignore a small spin canting, the direction of spin is exactly c-direction in P na ′ 2 ′ 1 magnetic group.It correspond to θ spin = θ t and ϕ spin = − π 4 .For P n ′ a ′ 2 1 , spin direction is a and corresponding angles are θ spin = π 2 and ϕ spin = π 4 .For the P na ′ 2 ′ 1 symmetry, the total moment is four times the a-component of the moment of the reference Cr/Cu1.By using the above relations, Likewise, for the P n ′ a ′ 2 1 symmetry, the total moment is four times the c-component of the moment of the reference Cr/Cu1.
The inclusion of the small spin canting (ϵ < 1 • ) would not significantly alter the predicted orbital magnetic moment.For example, the canting of 1 • changes only about a factor of 0.0175 (1 • in radians) of the results.
Let us consider the detailed aspect of the parametrization of the θ JT by λ.In the 1-Cu, as the λ increases from 0 to 1, JT mode Q 2 changes linearly from 0 to 0.288 as shown in Fig. S10 (d).Meanwhile, Q 3 changes very little, so it can be considered a constant (Fig. S10 (e)).λ = 0 and λ = 1 correspond to θ JT = π and θ JT = 1.934 ≈ 0.616π, respectively (Fig. S10 (f)).This can be considered that tan(π − θ JT ) is proportional to Q 2 .This leads to the parametrization in the main text.The simplified value θ JT,λ=1 = 2π/3 corresponds to The parameters used to estimate A are as follows.We adopted SOC parameter ζ = 56.54 meV for Cu and ζ = 27.52 meV for Cr, S11 ∆E = E 0 + − E 0 t 2g = 2.5 eV for both of the 1-Cu and 1-Cr.∆E is roughly chosen from the projected density of states obtained from the DFT calculation.The tilting angles of MO 6 octahedron extracted from the structures are θ t = 31.61• for 1-Cu and θ t = 31.88• for 1-Cr.

S10. The Second Order Energy Correction Term in The Perturbation Theory
In this section, we calculate the energy correction by the SOC in the perturbation approach, which gives rise to the MSIA.In the previous study, S7 only the fixed JT phase and the samespin contribution are considered for the MSIA.We improve the formulation by including the general JT phase and the opposite spin contribution.The first order energy correction ⟨d 0 nσ |H SOC |d 0 nσ ⟩ vanishes because the diagonal components of L d i are 0. The lowest order energy correction is the second order correction term, It is convenient to separate the H SOC into the same-spin block and opposite-spin block. and Similar to the orbital angular momentum case, the summation of the second-order correction to the energy with the non-degenerate assumption in the half-filling case with only up spins provides a shortcut for the calculation.
The first term of the last line is the summation of the same-spin contribution to energy correction which vanishes.On the other hand, the second term is the summation of the opposite-spin contribution and is non-vanishing in general.It makes the difference between the expression of the energy correction of the d 4 spin configuration and that of the d 9 .
The opposite-spin contributions from each orbital of the d 4 configuration are as follows The spin direction-dependent terms including θ and ϕ can be rewritten in physically intuitive expressions, where ŝ = s/|s| is the unit vector indicating the spin direction.Meanwhile, it implies the existence of the spin direction independent contribution to the energy correction.The samespin contribution to the correction is simply (∆E 2 d 4 ) ↑↑ = −(∆E 2 +↑ ) ↑↑ because the first term of the Eq.(S34) vanishes.
By summing up these terms, the spin direction-dependent part of the second order cor-rection to the energy in d 4 configuration is For the d 9 spin configuration, the second order correction to the energy can be obtained simply.
The spin direction dependent part is Terms in the spin direction-dependent part of the second order energy corrections are divided into three parts corresponding to the local coordinate components, (x, y, z).Each directional part is again divided according to the dependency on the JT phase.Orbital ordering following the JT distortion affects the energy correction in two ways.One is the JT transformation of the e g orbitals, which is explicitly expressed by the JT phase in Eq. (S41) and Eq.(S43).
The other is the changes in the SOC-unperturbed orbital energies E 0 i which is implicit in the expressions Eq. (S41) and Eq.(S43).Because the crystal field splitting is larger than the changes by the JT effect, we can consider that the factors in the trigonometric functions of the JT phase in the Eq.(S41) and Eq.(S43) are leading factors to determine the MSIA direction.We can ignore the JT dependency of E 0 i 's for simplicity.Note that the difference between the d 4 (Eq.(S41)) and d 9 (Eq.( S43)) cases rationalizes the difference in the favored axis of the spin between Cr-and 1-Cu, although a comparison by evaluating all those terms is too complicated.
In ∆E 2 d 4 (S), each directional part is composed of three sub-terms.Two of them are JT phase dependent and one is JT phase independent.Because we assumed E  Since we observed that the orbital magnetic moment is a dominant contribution to the magnetic moment of the 1-Cu and a local orbital magnetic moment is not parallel to the local spin moment in both 1-Cr/Cus, it is interesting to investigate the orbital texture.In our context, the orbital texture is the orbital angular momentum texture defined as l kn = ⟨ψ kn |L|ψ kn ⟩.For this purpose, we adopted another first-principles calculation package OpenMX, S12 which adopts the linear combination of the pseudo atomic orbital (LCPAO) method.Therefore, orbital angular momentum can be calculated at an atomic orbital angular momentum operator level.As a LCPAO basis, we adopted H6.0-s2p2, C6.0-s2p2, N6.0-s2p2, O6.0-s2p2d1, Cu6.0H-s2p2d2f1, and Cr6.0-s2p2d2f1, which indicate the number of the set of each s, p, d, and f orbital sets for each atomic species.Other computational settings are compatible with the setting for the VASP package.
As a preliminary step, we calculated the spin texture using OpenMX S13 to check if the two software packages give compatible results.As shown in Fig. S11, the main features of the spin textures are consistent: In the 1-Cu, the curly in-plane texture at k z = 0 and the alignment to z-direction at k z = π/c (Fig. S11 (a)); In the 1-Cr, the persistent type spin textures (Fig. S11 (b)).On the other hand, there are some differences in the case of 1-Cu.
The magnitude of the spins are larger in the OpenMX data at k z = 0, and the signs of the z components at k z = π/c show a different pattern.Nevertheless, we can expect that primary features will be shared in the OpenMX results.
The orbital textures calculated by OpenMX are shown in Fig. S12.The orbital angular momenta are calculated using an in-house post-processing program for OpenMX.Interestingly, the orbital textures are not parallel to the spin textures at all k-points.In the 1-Cu (Fig. S12 (a)), there are no out-of-plane components at k z = 0 plane, and their in-plane components show a different pattern from the spin texture.At k z = π/c, the out-of-plane orbital texture shows a similar texture with the spin besides their magnitude.However, there are also in-plane components in the orbital texture.In the 1-Cr (Fig. S12 (b)), interestingly, the in-plane orbital texture shows a Dresselhaus-type texture S14 at both k z = 0 and k z = π/c planes, i.e., the coexistence of the persistent spin texture and the Dresselhaus orbital texture.Also, there are small out-of-plane components in the orbital texture.

Figure S1 :
Figure S1: Asymmetric unit of 1-Cr at 300 K with atoms shown as thermal ellipsoids with 40 % probability.Atoms are colored as in Fig. 1.

Figure S2 :
Figure S2: Plot of the Cr-O bond distances in 1-Cr as a function of temperature.

Figure S3 :
Figure S3: Plot of lattice parameters versus temperature for 1-Cr.

Figure S4 :
Figure S4: Depiction of the evolution of the framework angles along the a and b axis for 1-Cr.Insert indicates the orientation of these framework angles in the structure, atoms are shown in the same color in Fig. 1.Hydrogen bonds are shown between the [C(NH 2 ) 3 ] + cation and the anionic [Cr(HCOO) 3 ] − framework as dashed gold lines; each hydrogen of the [C(NH 2 ) 3 ] + cation forms a single hydrogen bond with an oxygen atom from a formate ligand with D...A distances of 2.933(19) to 3.019(13) Å.

Figure S5 :
Figure S5: The spin configurations corresponding to the ferromagnetic (M ), G-type AFM (G), A-type AFM (A), and C-type AFM (C) orders, respectively.Upward red and downward blue arrows mean the spin-up and spin-down, respectively.The dotted lines represent the neighboring spin pairs and the solid line box is a half of the unit cell of 1-Cr/Cu by {a, b, c/2}.The figure shows the c-component case of each order, such as A c .The b(c)component can be obtained by setting the b(c)-axis as the spin-axis.
and θ t,i the structural tilting angle of CrX 6 octahedra.Note that we are using the different notation θ t from Ref.S7 (α) for the tilting angle for consistency with the main text.The global cartesian reference frame has the spin quantization axis z parallel to the crystal axis c, while the x, y axes are rotated by 45 • with respect to the a, b axes.S7

Figure S6 :
FigureS6: DFT total energies of 1-Cr depending on the canting angle ϵ in ac-plane.Without the U and J corrections (no U), the total energy is minimum at ϵ = 0.63 • .With the DFT+U + J correction, the minima are at ϵ = 1.19 • and ϵ = 1.86 • for the parameter set 1 (2.5,0.5) and set 2 (3,1), respectively.In a similar way, the optimal canting angle of the G b order is δ = 2.01 • (no U).

Figure S7 :Figure S8 :
Figure S7: Monte Carlo simulations (set 1).Left panel: specific heat and A-type AFM susceptibility per unit cell as a function of temperature, both displaying a peak at temperatures within 14 and 15 K. Central panel: evolution with temperature of the Atype AFM order parameter alongside secondary order parameters.Note that the unit of susceptibility here is β = 1/k B T in meV, as a result of dimensionless spin in the simulation.Right panel: zoomed view of secondary order parameters, highlighting the small G-type AFM and weak FM components developing at the critical temperature.Results shown for L = 24.

Figure S9 :
Figure S9: Monte Carlo simulations (set 2).Powder average of the magnetic susceptibility converted to a physical unit, showing a hump around the transition temperature.Results shown for L = 24.The labels mean in which way the unit is converted.
only the second term remains up to the first order in ζ.For the summation over the corresponding terms to Eq. (S24) of occupied orbitals, the following property can be exploited.If we sum up all the orbital angular momentum expectation values of the perturbed d-orbitals with up-spin, n ⟨d n↑ |L i |d n↑ ⟩ = L i d 0 m↑ + c.c. = 0 (S25) because each of the terms is canceled with the term whose n and m are exchanged.The same holds for down spin.It makes it easy to calculate the orbital angular momentum for d 4 and d 9 .For the spin-up high spin configuration of the Cr 2+ ion (d 4 ), occupied d-orbitals are {d yz↑ , d zx↑ , d xy↑ , d −↑ }.Therefore, ⟨L i ⟩ d 4 = − ⟨d +↑ |L i |d +↑ ⟩ for d 4 .For the Cu 2+ ion (d 9 ), only d +↓ is unoccupied.For d 9 , by using d 0 m↑ H SOC d 0 n↑ = − d 0 m↓ H SOC d 0 n↓ from the opposite signs of the diagonal components of the spin matrix in Eq. (S15) and d 0 n↑

Figure S10 :
Figure S10: (a-c) P na2 1 structure 1-Cu (λ = 1).The reference Cu1 is labeled.In (c), elongated directions are drawn by double-arrows.JT mode amplitudes (d) Q 2 and (e) Q 3 of the reference Cu1 and Cr1 and their (f) JT phase are shown with respect to λ.
we can determine whether each of the JT phase dependent sub-part is the energy-lowering or -raising term according to the sign of these energy-related factors.They are classified in Table.S5.JT phase independent terms are always energy-raising terms.Due to the phase difference by π of the JT phase dependent factors, the energy-lowering and -raising terms result in the same tendencies of the favored spin direction at a given θ JT .On the other hand, ∆E 2 d 9 (S) has only the energy-lowering

Figure S11 :
Figure S11: The spin textures calculated by using OpenMX.(a) 1-Cu and (b) 1-Cr, at the k z = 0 (upper panels) and k z = π/c (lower panels) planes.The in-plane x and y components are represented by arrows of which the length is the in-plane magnitude with respect to the reference spin 1/2 above the figure.The z components are represented by colormap of dots.The inner boundary of each figure coincides with the Brillouin zone boundary.

Figure S12 :
Figure S12: The orbital angular momentum textures calculated by using OpenMX.(a) 1-Cu and (b) 1-Cr, at the k z = 0 (upper panels) and k z = π/c (lower panels) planes.The in-plane x and y components are represented by arrows of which the length is the in-plane magnitude with respect to the reference orbital angular momentum 1 above the figure.The z components are represented by colormap of dots.The inner boundary of each figure coincides with the Brillouin zone boundary.

Table S4 :
DFT parameter dependence of the calculated magnetic moments.Units are Bohr magneton µ B .|M orb | is the magnitude of orbital magnetic moment of single Cr/Cu ion.