Spin-orbital mechanisms for negative thermal expansion in Ca2RuO$_4$

The phenomenon of negative thermal expansion (NTE) deals with the increase of the lattice parameters and the volume of the unit cell when the material is thermally cooled. The NTE is typically associated with thermal phonons and anomalous spin-lattice coupling at low temperatures. However, the underlying mechanisms in the presence of strong electron correlations in multi-orbital systems are not yet fully established. Here, we investigate the role of Coulomb interaction in the presence of lattice distortions in setting out the NTE effect, by focusing on the physical case of layered Ca$_2$RuO$_4$ with $d^4$ configuration at each Ru ion site. We employ the Slater-Koster parametrization to describe the electron-lattice coupling through the dependence of the $d-p$ hybridization on the Ru-O-Ru bond angle. The evaluation of the minimum of the free energy at finite temperature by fully solving the multi-orbital many-body problem on finite size cluster allows us to identify the regime for which the system is prone to exhibit NTE effects. The analysis shows that the nature of the spin-orbital correlations is relevant to drive the reduction of the bond angle by cooling, and in turn the tendency toward a NTE. This is confirmed by the fact that a changeover of the electronic and orbital configuration from $d^4$ to $d^3$ by transition metal substitution is shown to favor the occurrence of NTE in Ca$_2$RuO$_4$. This finding is in agreement with the experimental observations of a NTE effect which is significantly dependent on the transition metal substitution in the Ca$_2$RuO$_4$ compound.


I. INTRODUCTION
Nowadays it is widely accepted that strong electron correlations play a prominent role in achieving and controlling a large variety of quantum phases, physical phenomena and effects related to cooperative behaviors that cannot be described from the properties of individual electrons. Emblematic cases are the unconventional pairing in high T c superconductors [1], as well as metalinsulator transitions [2] or emergent excitations such as monopoles [3] and skyrmions [4]. Transition metal oxides (TMOs) are paradigmatic quantum materials [5] in this context, as they are marked by strong correlations that result in several types of electronic phases, including charge-and spin-ordered states. Prominent examples are found among the 3d TMOs, where unconventional transport properties, ordering phenomena, and unusual spectroscopic properties are observed. It was argued that the comparably weak spatial extension of 3d orbitals may lead to large electronic Coulomb interactions, competing with kinetic contributions. Depending on crystal-field, hybridization, Hund's exchange, and band filling, this interplay can lead to renormalized metallic behavior or induce Mott insulating behavior. In 4d and 5d oxides, spin-orbit coupling acts on an energy scale comparable to the other energy scales of the system, and the observed electronic state is the result of a complex interplay of Coulomb interactions, spin-orbit splitting, and crystal field effects.
Interestingly, when the electron correlations are somewhat weaker, as in 4d or 5d TMOs, other challenging effects may arise due to the subtle competition between electronic and lattice degrees of freedom. One of most extraordinary one is the negative thermal expansion (NTE) in Ca 2 RuO 4 -derived layered materials. The phenomenon of NTE [6][7][8] deals with the observation of a shrinking of the lattice parameters when the material is heated or vice versa with a lattice expansion by thermal cooling.
Although the phenomenon is not broadly observed, it has a remarkable impact for several applications in electronics, optics, and for the design of thermal engines or medical products [9,10]. The origin of NTE is quite intricate as, apart from the structural modifications, it can involve the coupling of electrons, spin, orbital degrees of freedom. In materials such as ZrW 2 O 8 and ZnCr 2 Se 4 , the NTE is mainly due to oxygen vibrational modes [11][12][13][14][15] and anharmonicity related to spin canting and spin-lattice coupling. Hence, while the manifestation of the effect is directly provided by the modification of the lattice parameters, a fundamental question arises on whether the electronic correlations, and the resulting spin-orbital couplings, would cooperate or compete with the tendency of the lattice to exhibit a NTE. This issue is particularly relevant in materials with Mott physics and nontrivial magnetic ordering, with electronic degrees of freedom being dominated by local Coulomb interaction, especially in the presence of multi-orbital configurations having significant spin-orbital entanglement.
Looking more inside at the Ca 2 RuO 4 family, it has been reported that a tiny percentage of Cr substituting Ru produces a 1% volume reduction that has been ascribed to the collapse of the orbital and magnetic ordering [16]. Remarkably, substituting Ru with M ion (M =Cr, Mn, Fe, or Cu) generally yields NTE effects in Ca 2 Ru 1−x M x O 4 [17][18][19], with the onset following the behavior of the metal-insulator and magnetic ordering transition temperatures [17]. Therefore, the doping by transition metal (TM) ions plays an important role in unveiling a generic tendency toward the phenomenon of NTE and altering the strength of the spin-orbital interactions while keeping the system in the insulating regime [20,22,23].
It is thus compelling to attribute the observed NTE to a mechanism where electronic correlations play a critical role. Insulating Ca 2 RuO 4 is a representative case of a multi-orbital correlated TMO, where both spins and orbitals develop their own dynamics and are coupled to each other through lattice distortions. In this frame, the mechanism of NTE can thus contain a high degree of complexity. In TMOs, the distance of the TM ions is related to the TM−O−TM bond angle which in turn is determined by the rotation of the TM−O 6 octahedra around a given crystallographic axis. Then, a reduction of the lattice parameters can be obtained by a deviation from 180 • of the TM−O−TM bond angle (see Fig. 1). While the reduction of the TM−O−TM angle is expected to reduce the kinetic energy in an insulator described by a single orbital, the energy balance for multi-orbital electronic systems is not obvious, especially in the presence of Coulomb interaction.
Therefore, various cooperative effects may materialize. The electron localization reduces the kinetic energy of the electrons and, simultaneously, the lattice expansion may cost energy due to electron-lattice interaction. Hence, the lattice expansion may influence the kinetic energy variation even further. Moreover, the spin-lattice coupling also plays an important role, particularly when strong competition between ferromagnetic (FM) and antiferromagnetic (AF) interaction occurs since this interplay often leads to a strong bond frustration or lattice anomaly. Interestingly, the frustration may be released in a phase transition by which a certain type of magnetic order is stabilized.
Such competing mechanisms depend on the TM−O−TM bond angle, which directly enters in setting out the most favorable structural configurations. In this context, the short-range orbital or magnetic or spin-orbital correlations play a crucial role and decide about the energy gain or loss when a variation of the lattice parameters occurs.
Taking into account these multiple components entering into the phenomenology of the NTE, we aim to investigate such physical scenario by assessing the role of Coulomb interaction and spin-orbital correlations in setting out the thermal evolution of the TM−O−TM lattice configuration. In order to determine the impact of the TM−O−TM configuration in the electron correlated configuration, we employ the Slater-Koster parametrization to describe the electron-lattice coupling through the dependence of the d − p hybridization on the TM−O−TM bond angle. The evaluation of the minimum of the free energy at finite temperature is achieved by a complete solution of the multi-orbital many-body problems on finite-size clusters. This approach allows us to identify the electronic regime where the system can exhibit NTE effects. The analysis shows that the nature of spin-orbital correlations is relevant to drive the reduction of the Ru−O−Ru bond angle by thermal cooling, and in turn the tendency to develop a NTE. Indeed, the changeover of the electronic and orbital configuration from d 4 to d 3 by TM substitution clearly supports the role of electron correlations to account for the occurrence of NTE in Ca 2 RuO 4 . This finding agrees with the experimental observation of a NTE effect, which exhibits significant dependence on the TM substitution in the Ca 2 RuO 4 compound.
The paper is organized as follows. First, in Sec. II we consider a single TM−O−TM bond with two d 4 ions and investigate whether this bond may generate a configuration compatible with the occurrence of a NTE effect. Second, we show that the TM−O−TM bond with one d 3 ion substituting the d 4 ion has a tendency to favor bond distortions that in general support the NTE in a planar perovskite. The search for an optimal bond angle configuration in the presence of Coulomb interaction at the TM site is further analyzed with a 2 × 2 plaquette. We derive an effective low-energy model that properly includes all the d − p charge transfer processes, the bond angles, and the local Coulomb interaction. In Sec. III we investigate the case with d 4 charge configurations while in Sec. IV we consider the case of a single d 3 impurity replacing a d 4 site. The study of the planar plaquette allows us to address the role of orbital correlations in setting out electronic states which are compatible with the NTE effects. In particular, orbital configurations or electronic mechanisms that break the rotation by 90 • in the square lattice are particularly relevant for favoring lattice distortions that support the NTE effects in Ca 2 RuO 4 . Our main conclusions are summarized in Sec. V. The superexchange Hamiltonians are explicitly derived and reported together with the complete set of their coefficients in the Appendices A and B. In the Appendix C we provide details about the evolution of the optimal angle as a function of the Hund's coupling and of the spin-orbit interaction. In the Appendix D we present the results of the effective spin-orbital model for a single bond with d 4 − d 4 and d 3 − d 4 configurations. Finally, in the Appendix E we introduce a phenomenological electron-lattice potential to describe the lattice feedback on the electronic system that is able to limit the bond angle within a physical range. is realized when Coulomb interaction at the TM element is fully taken into account. To this end, we consider a multi-orbital p − d electronic system suitable for oxide materials. This includes all the atomic terms arising from the local Coulomb interaction and the orbital-dependent connectivity between a TM element and the oxygen ligands that are allowed in the planar square lattice geometry, typical for tetragonal perovskites. First, we consider the case of a single bond made of two TM elements, i.e., TM 1 and TM 2 , linked via an oxygen (O) ion between them, as schematically depicted in Fig. 1.
Note that the deviation of the bond angle away from 180 • is directly related to the angle θ, formed by the TM−O distance and the crystal axis. Hence, in the following, we will refer to θ as a measure of the bond angle, with θ=0 corresponding to an undistorted 180 • bond, and the maximal θ = 45 • corresponding to a 90 • TM−O−TM bond. The fundamental physical motivation for this study is to assess a role that can be played by the electron-electron correlations in setting the bond angle between TM elements and including their effects on the thermal evolution of the bond angle.
Within this framework, we aim to evaluate the role of spin, orbital and charge-transfer processes in setting out the energetically most favorable TM−O−TM angle. In particular, we consider how the Coulomb interaction at the TM site leads to the optimal TM−O−TM bond angle and whether the resulting configuration is related to the spin-orbital correlations or other electronic parameters associated with the tetragonal distortions and spin-orbit coupling. The analysis is performed by solving exactly a model Hamiltonian for the TM 1 −O−TM 2 cluster, with two TM ions at its ends.
To describe the TM 1 −O−TM 2 cluster, we employ a  [24][25][26] includes the complete Coulomb interaction projected on the t 2g electrons, the spin-orbit coupling, and the tetragonal crystal field (CF) potential. Therefore, the H loc at site i is, Here, the on-site Coulomb, spin-orbit, and CF terms are expressed at site i by

(2.4)
In these expressions, i labels the site and {α, β} are indices running over the three orbitals in the t 2g sector, i.e., α, β ∈ {d xy , d xz , d yz }, and d † iασ is the creation operator of an electron with spin σ at the site i in the orbital α. The interaction is parametrized by the Kanamori parameters: intra-orbital Coulomb interaction U and Hund's exchange J H .
The strength of the tetragonal distortions is expressed by the amplitude δ, with δ = (ε xy −ε z ). We also consider the possibility of having an orthorhombic splitting, δ ort of the {xz, yz} orbitals by assuming that ε yz = ε z + δ ort and ε zx = ε z − δ ort . U is the intra-orbital Coulomb interaction, J H Hund's exchange, and (U − 1 2 J H ) sets the strength of the inter-orbital electron-electron interaction, J is the pair hopping term. We assume a rotational invariant condition for the Coulomb amplitudes, so that U = U + 2J H , and J = J H . The operator − → l is the projection of the angular momentum operators to the t 2g subspace, (l k ) αβ = i kαβ such as The matrices for the orbital operators in the t 2g manifold are: Concerning the local Hamiltonian at the O site, it only includes the on-site energy term which is introduced to take into account the energy difference between the occupied orbitals of O and TM.
H O el (j) = ε x n j,px + ε y n j,py + ε z n j,pz (2.5) Furthermore, we consider the TM-oxygen hopping, which includes all the allowed symmetry terms according to the Slater-Koster rules [27,28] for a given bond connecting a TM to an oxygen atom along a given symmetry direction, e.g. the x-axis. Here, we allow for rotation of the octahedra around the c-axis assuming that the TM−O−TM bond can form an angle θ as depicted in Fig. 1. The case with θ = 0 corresponds to the tetragonal undistorted bond, while a non-vanishing value of θ arises when the TM−O 6 octahedra are rotated at the corresponding angle around the c-axis.
For the TM−O hopping term, we assume that for a generic bond connecting the TM to the O atoms along the x-direction, the d − p hybridization includes all the allowed terms, and thus where the hopping t dα,p β is a function of the bond angle θ. In particular, according to the Slater-Koster rules, one has: with n x = cos θ and n y = sin θ. In a similar way, one can also express the p − d hybridization amplitude for the O-TM 2 bond. By symmetry correspondence, one can write down the other hybridization terms along the other TM-O-TM bonds for the y-direction.
Since we are interested in the total free energy at finite temperature as a function of the bond tilting angle, it is useful to introduce the expression which is generally given by where β = 1/k B T is the inverse temperature with k B being the Boltzmann constant, while E i are the eigenvalues of the Hamiltonian evaluated by exact diagonalization.
To proceed further with the analysis, we determine the whole energy spectrum and the corresponding eigenstates for the TM−O−TM cluster. Having in mind the Ca 2 RuO 4 system, we start our analysis by fixing the atomic electronic parameters in a regime matching with the AF ground state of the monolayer compound [29][30][31][32]. There, the severe flattening of the RuO 6 octahedron occurring below the structural transition, corresponds to negative values for δ, while, according to first principle calculations or estimates employed to reproduce the resonant inelastic x-ray and the magnetic properties, its amplitude is in the range |δ| ∈ [200, 300] meV [30,33].
Furthermore, it is useful to point out that materialspecific values such as λ = 0.075 eV, U ∈ [1.5, 2.5] eV, and J H ∈ [0.35, 0.5] eV is taken as a reference for the current study. Similar values for δ, U and J H have been used for calculations of electronic spectra in Ca 2 RuO 4 and the ratio g = δ/(2λ) is typically considered to lie in the range ∈ [1.5, 2] for modeling the spin excitations observed by a neutron, Raman, and resonant inelastic x-ray scattering [30,[34][35][36]. For the hopping amplitudes, according to first principle calculations, one can assume the p − d values to be in the range of [1.0,2.0] eV [38]. These values are in the range of those employed to describe the electronic and magnetic properties of ruthenate oxides [37][38][39]. For V pdπ about larger than V pdσ , one finds a tendency to favor a small tilt of the bond angle with an amplitude that is close to zero. Otherwise, in the remaining portion of the phase diagram, the trend is to stabilize a large TM−O−TM bond angle, of π/4 (then all bonds are 90 • ). The transition between the two regimes is rapid and generally occurs in proximity to the line V pdσ ∼ V pdπ .
We point out that a similar trend is also obtained for the TM(d 3 )−O−TM(d 4 ) bond configuration and that a variation of the electronic parameters (e.g. U , J H , etc.) within the physical regions previously defined do not alter significantly the structure of the phase diagram. This trend can be qualitatively captured by considering the dependence of the d − p hybridization on the TM-O-TM bond angle. For large V pdπ the d − p charge-transfer and the resulting kinetic energy is optimized by θ 0. In the opposite regime with large V pdσ a significant deviation from θ 0 increases the kinetic energy due to the d − p hopping processes.
Let us then consider the thermal evolution of the optimal bond angle. To this aim, at each temperature, we has a nonzero value in the ground state, thus implying that the bond is distorted at low temperatures. Then, we track the evolution of the angle θ opt that minimizes the free energy as a function of temperature, by varying the Coulomb interaction U and the strength of Hund's exchange.
In Figs. 3(a) and 3(b) we show a representative behavior of the optimal bond angle for the TM(d 4 )-O-TM(d 4 ) configuration as a function of the temperature. As one can see by inspection of Fig. 3, the intra-orbital Coulomb interaction U is able to drive a transition from a regime where the bond angle is insensitive to the temperature change and it optimizes the total electronic energy by keeping the bond angle undistorted, to another a regime of strong coupling where the high-temperature bond angle is large and exhibits a tendency to decrease by reducing the temperature. On the other hand, we find that the strength of Hund's exchange also affects the critical value of the Coulomb interaction where the transition to an undistorted bond angle occurs, such transition being favored for higher values of J H . Moreover, it affects the crossover the temperature associated with a variation of the bond angle. Such dependence on Hund's coupling clearly indicates that the spin and orbital correlations are linked to the thermal changeover of the bond angle configuration.
In order to assess the microscopic mechanisms behind the setting of the energetically most favorable bond angle for the TM−O−TM configuration, we focus on the spin correlations between the spin moments at the TM sites by looking at both the scalar and vector product between spin operators. This study aims to search for a link between the spin configuration and the optimal bond distortion.
In Fig. 4 we report few representative cases for TM(d 4 )−O−TM(d 4 ) configuration to illustrate the overall trend. The bond-angle evolution of the thermal average of the scalar product between the spin moments at the TM sites (i.e., the Heisenberg interaction) indicates that the amplitude is almost temperature independent when x or z spin-spin correlation is considered, whereas y-y spin correlation tends to become of AF-type with decreasing the temperature [see Figs. 4(a) and 4(c)], if the bond becomes undistorted, i.e., when θ = 0. Moreover, looking at the vector spin-spin correlation, we note that its out-of-plane component gets larger than the planar one, whose value is almost vanishing. We also notice that the amplitude of the bond angle that minimizes the free energy does not correspond to an angle configuration that maximizes the spin correlations for both the parallel and perpendicular spin configuration. pation has been changed at one TM site, by replacing d 4 with a d 3 impurity. The temperature dependence of the optimal bond angle is completely different if compared to the TM(d 4 )-O-TM(d 4 ) bond. Indeed, we find that the general trend is to observe a decrease in the bond angle toward a more undistorted configuration for all the considered values of the Coulomb interaction. Although the minimization of the free energy yields large values of the bond angle at high temperatures it is relevant to point out the trend of the electronic driving force in distorting the TM-O-TM bond. Interestingly, an opposite behavior is observed when U increases and J H decreases; in this case indeed the bond angle increase when the temperature is lowered.
Focusing on the spin-spin scalar and vector product correlations between the TM sites, we find that the evolution of the thermal average of the scalar product upon the bond angle indicates that the amplitude is maximal when the bond is undistorted, i.e., θ = 0 [see Figs. 5(c) and 5(d)]. On the other hand, we note that the the amplitude of the bond angle that minimizes the free energy at low temperatures do correspond to an angle configuration which favors in-plane parallel spin configurations, with a prevalent component along the bond direction.
While this is expected in a single-band picture where the Heisenberg exchange depends on the hopping amplitude, i.e. J ∼ t 2 /U , and the charge transfer strength t is typically reduced by tilting the bond angle, in a multiorbital scenario as that one upon examination the compe-tition between the various orbital channels does not allow to have a direct prediction. On the other hand, the vector product spin correlations (i.e., the Dzyaloshinskii-Moriya (DM) interaction) are activated by distorting the bond as expected due to the inversion symmetry breaking arising from the bond tilting, see Figs. 5(e) and 5(f).
The general outcome is to have a maximal value of the vector product spin correlations for a bond angle that is nonvanishing. Moreover, as we explicitly demonstrate in Fig. 5(e)&5(f), the maximal amplitude of the DM exchange is also temperature dependent. These results indicate that the setting of the optimal angle in a TM(d 3 )-O-TM(d 4 ) bond is a consequence of the subtle competition between the Heisenberg and DM interactions, which are in turn strongly tied to the orbital degrees of freedom and to spin-orbital coupling strength. A second relevant outcome of the analysis is that the optimal bond angle for the ground state in terms of the hybridization d − p parameters generally show that there are two electronic regimes corresponding to a region of the phase space above (below) the diagonal in the {V pdσ , V pdπ } with small (large) bond angles, respectively. This is a general trend and the boundary can slightly move when varying the electronic parameters as the local Coulomb interaction and other atomic terms as the crystal field potential.
Finally, the evolution of the magnetic state as a function of the bond angle indicates that the configuration with minimum energy does not correspond to a state where there are maximal spin correlations for both the scalar and vector product between the spins at the TM sites. This implies that the variation of the bond angle results from a competition between the tendency to have parallel spin moments (bond angle about zero) or perpendicular ones when the inversion symmetry breaking is more pronounced (bond angle different from zero). In this context turns out that spin-orbit coupling is a relevant term, as it couples the spins to the orbital momenta and thus to the spatial orientation of the bond angle.
In the TM(d 3 )−O−TM(d 4 ) configuration, the spin-orbit coupling is inactive at the d 3 site, so we argue that its suppression can account for the tendency to observe an optimal bond angle which at low temperatures is always close to zero.
We point out that according to the relation U = U − 2J H , if J H /U > 1/3 then J H > U , which in turn sets a boundary for the crossover from weak to strong Hund's coupling. The solution of the multi-orbital Hubbard model for the single bond has been explored to investigate the changeover from weak to strong Hund's coupling regime. Indeed, the amplitude of Hund's interaction, as demonstrated in Fig. 3 (a,b), plays a relevant role in setting out the thermal negative expansion effects. This is evident, for instance, in the case of d 4 − O − d 4 bond where one finds a modification of the optimal angle with a positive thermal expansion effect (i.e. increase of the bond angle at low temperature), that turns into a negative thermal one when the amplitude of Hund's coupling grows. Along this line, an enhancement of the NTE effects, with a larger negative variation of the optimal angle and a shift of its onset towards higher temperatures are also observed for the doped d 3 − O − d 4 bond configuration ( Fig. 3 (c,d)).
The presented analysis is applicable to all oxide materials. However, since there is a clearcut experimental evidence of anomalous negative volume effects in the insulating phase of the Ca 2 RuO 4 compound and derived doped materials belong to the same family, we focus on an electronic regime that is relevant for ruthenium oxides. There, the Ru atom is in a d 4 configuration and one allows for atomic substitutions that are isovalent and can lead to a change into d 3 or d 2 electronic state at the TM site. This variation can be, for instance, induced by replacing the Ru with Mn or Cr ions for achieving a d 3 or d 2 configuration, respectively.

III. PLAQUETTE CONSISTING OF d 4 IONS: LOW-ENERGY SPIN-ORBITAL DESCRIPTION
In order to account for the electronic behavior in two dimensions, a minimal requirement is to simulate bonds in two inequivalent directions. To this aim, we choose to study a cluster based on a 2 × 2 plaquette. This choice is motivated by the fact that it is highly computationally demanding to deal with the full d-p multiband Hubbard model as the size of the Hilbert space rapidly grows. Therefore, we resort to an approximate treatment valid in the insulating regime with frozen charge degrees of freedom due to a large U coupling leading to an effective spin-orbital exchange. The spin-orbital Hamiltonian for the d 4 configuration with three t 2g orbitals has been already derived [40,41] for a basic symmetric configuration.
Here, instead, we need two new microscopic ingredients: (i) a non-trivial TM−O−TM bond angle, and (ii) the oxygen degrees of freedom in the virtual excited states. Note that (i) can be implemented without (ii) if we neglect the oxygen degrees of freedom so the effect of bond distortion can, to some extent, be captured by a simplified model. Nevertheless, we use the full model with oxygens because of the lack of the oxygen degrees of freedom cannot properly describe a lattice configuration with a variable bond angle. Indeed, we have found that the simplified model with projecting out oxygen degrees of freedom never favors a non-zero bond angle in the ground state which does not match the results of the d − p multi-band Hubbard model.
The detailed derivation of the TM−O−TM spinorbital exchange up to the fourth order is given in Appendix A. The resulting Hamiltonian couples two spins S = 1 and orbitals L = 1 along a bond. The bond Hamiltonian depends on the direction of the bond in the lattice and the bond angle. We adopt an idealized 2 × 2 plaquette with fixed TM−O distances and the distances between TM ions being reduced by cooperative distortions. The schematic view of the plaquette and the convention chosen for the site indices and for the angles indicating the rotations of RuO 6 octahedra are given in Fig. 6. Therefore, the Hamiltonian for the 2 × 2 cluster of four equivalent d 4 ions can be written as, where a † i , b † i and c † i are the hard-core boson operators that create a double occupation of the yz, zx and xy orbital, respectively.
We study the ground state and finite-temperature properties of the Hamiltonian (3.1), namely the behavior of the bond angle θ that minimizes the free energy of the system as given by where the sum runs over the eigenvalues {E n (θ)} of H. We also investigate the behavior of the bond angle in presence of a magnetic order with ferromagnetic alignment along one direction and AF in the other, thus realizing a so-called stripe C-AF pattern. This configuration is stabilized by adding an effective Zeeman field that breaks the symmetry and aligns the spin S and orbital L angular moments along the chosen orientation z: with J i = S i + L i . In Fig. 7(a) we show a phase diagram of a plaquette as a function of the hybridization amplitudes π and σ, per analogy to the one obtained for a single bond shown in Fig. 2. Here, at first glance, we see only two phases: the one with zero bond angle (straight bonds) when π hybridization is dominating, and the one with π/4 bond angle when σ hybridization dominates. However, looking closer one notices a very narrow intermediate phase along the phase boundary when the bond angle takes intermediate values. This is indicated by the behavior of a representative free energy curve found in this phase, see Fig. 7(b), having a minimum for an intermediate value of the optimal bond angle θ opt . One can then immediately track the behavior when going from one phase to the other. We point out that this peculiar behavior of the optimal angle, with only a small window of parameters for the d−p hybridization amplitude showing a rapid variation from small to large angles, can be regularized by including the effects of the lattice potential. This effect is demonstrated in the Appendix E. It is shown indeed that, depending on the electron-lattice coupling amplitude, one can enlarge the range of electronic parameters for the d − p hopping amplitudes where a smooth variation from large to small angle takes place. In Fig. 7(c) we show the curve of the optimal bond versus pd − π hybridization when crossing the phase boundary. We observe that it interpolates smoothly between π/4 and 0, leading to a second-order type transition between the two phases. Interestingly, when searching for an observable that would follow the behavior of the optimal angle we find that it is optimally reproduced by the cross-product of the spins along a given nearest-neighbor bond. From the symmetry property of the system, it follows that only the z component of such a cross-product takes non-vanishing values. Hence, in Fig. 7(d) we show the average of ( S 1 × S 2 ) z in the ground state (note that here all the bonds of the plaquette are equivalent). We see that the profile of the curve resembles very closely the shape of θ opt . The non-vanishing average of the crossproduct of spins means that the system tends to stabilize a non-collinear magnetic order. Therefore, we can conclude that in the examined plaquette the change of the optimal bond angle is closely related to the tendency of the system to develop a non-collinear spin texture.
The presence of an intermediate phase suggests a frustrated configuration of the plaquette in a well-defined parameter area. Thus, we study the thermal behavior of the system for the intermediate phase with the goal to achieve a configuration with a NTE character. Here, the manifestation of this effect can be tracked by identifying the configurations with a positive derivative of the  ) typical curves of the optimal bond angle and the ground state averages of S1 × S2 z and S1 × S3 z versus δ(V pdπ ) = V pdπ − V0 for V pdσ = 1.5 eV in the intermediate phase between θ = 0 and θ = π/4 (see Fig. 6); (c,d) thermal dependencies of the optimal bond angle and average S1 × S2 z and S1 × S3 z in the intermediate phase for optimal angle with respect to temperature, i.e., the system should become more distorted (bond angle different from zero) when the temperature grows. In Fig. 7(e) we show the behavior of the optimal angle as a function of temperature in the region of the parameter that is in close proximity to the intermediate phase. We see that initially θ opt stays constant when the temperature grows up to about T = 1.5 meV but then rapidly drops to zero at around T = 4.0 meV and stays zero until reaching a high-temperature regime. This means that no NTE effect can be obtained in this regime of parameters. Moreover, the behavior of θ opt is non-analytical at the two transition temperatures and this behavior is again well reproduced by the chiral spin correlations shown in Fig.  7(f). This outcome indicates that the thermal behavior of the system is also related to the non-collinearity of the spin texture. Finally, to trigger the NTE effect we have considered the differences between the plaquette and the former single bond analysis where the NTE effect was indeed obtained. We argue that the main difference between the two systems is the underlying symmetry. The plaquette obeys a fourfold in-plane rotation symmetry with which both main phases are compatible. This suggests that the NTE effect could be obtained by lowering the symmetry of the plaquette. To achieve this we have employed the striped C-AF magnetic order introduced at the beginning of this section.
The results at zero and finite temperatures are shown in Fig. 8-note that we do not show the phase diagram again because it changes only slightly. The behavior of the optimal angle when passing through the intermediate phase is shown in Fig. 8(a). Note that the angle still interpolates smoothly between 0 and π/4 but the drop of the function is much more rapid than before. Again, this behavior is well reproduced by the cross-product of the spins, see Fig. 8(b). Looking at the thermal behavior of θ opt in the intermediate phase, see Fig. 8(c), we see a sharp increase of the angle starting roughly at 20 meV and at around 34 meV. Therefore we get the NTE effect within this temperature range. Interestingly the θ opt seems to be non-analytical only at the larger temperature. Finally, we observe again in Fig. 8(d) that the spin cross-product follows the optimal angle, including the non-analytical point.
The qualitative outcome is not affected by the increase of the J H /U ratio, as inferred from the results reported in Appendix C where we have analyzed how the thermal profile of the bond angle is affected by a change of the Hund's coupling and spin-orbit coupling. Moreover, as shown in Appendix D where we plotted the results for the low energy spin-orbital model for the case of a single bond, we find that the qualitative trend is similar to that obtained by solving the full multi-orbital model in the regime of strong U coupling. Similarly to what has been done in the previous section, we have derived a superexchange Hamiltonian of the plaquette with a single d 3 dopant to be able to solve the problem in the limit of large U . A single dopant in the d 4 plaquette results in two d 4 − d 4 bonds between the host sites and two-hybrid d 3 − d 4 bonds between impurity and host atoms. Therefore, the hybrid bonds need to be derived including the effect of the oxygen degrees of freedom and bond distortions, as in the pure d 4 case. The full derivation of the effective model is reported in Appendix B. Note that the case with no distortion and no oxygens was already addressed in Ref. [20].
It is worth mentioning here that while for pure host bonds, having large U alone is enough to stabilize the low-energy spin-orbital model, in the case of the hybrid bonds it is necessary to add a mismatch between the atomic levels of the d 3 atom with respect to the d 4 one in order to secure different valences of these atoms. As a result, one gets an exchange model between spin S = 1 and orbital L = 1 at one site and the impurity spin S = 3/2 T �meV� with no orbital degree of freedom at the impurity. Taking the same angle convention as that shown in Fig. 6, the Hamiltonian with a d 3 impurity at site i = 4 reads as follows [cf. Eq. (3.1)], Here, we follow the same logic as in the previous section without d 3 impurity. The results are shown in Fig.  9. As in the previous case, the phase diagram of Fig.  9(a) exhibits two phases, one with zero bond angle and the other with an angle π/4. Their positions are similar  as in the pure d 4 case, however, in contrast to that configuration, we find no phase with an intermediate angle phase between these two. This follows from the behavior of the free energy, see Fig. 9(b), which exhibits two local minima at θ = 0 and θ = π/4. When one passes through the phase boundary the positions of the minima do not change but the values of the minima do. As a result, we get a discontinuous, first-order transition between the two phases.
In this context, the impact of having a finite temperature is quite peculiar. If we place ourselves close to the phase boundary we find that the system jumps between two optimal angles θ = 0 and θ = π/4 as the temperature grows, see Fig. 9(c). We argue that in the infinite system there will be phase separation at the temperatures when two phases are degenerate which will make the jumps smoother, leading to the temperature intervals hosting the NTE effect with NTE.
The situation becomes more evident in the C-AF phase. The transition at zero temperature remains first order but the optimal angle first drops slightly from the π/4 value before jumping to zero when crossing the phase boundary. If we now place ourselves in the region where the angle is decreased and look at the thermal behavior, see Fig. 9(d), we get a NTE effect in the temperature range roughly between 1.5 and 3.2 meV. The θ opt curve shows a non-analytic behavior at the latter point. Unlike in the pure d 4 case, looking at the spin-spin cross products, see Figs. 9(e-f), we do not observe that the curves follow θ opt very closely, however in the NTE interval they exhibit high gradients. Note that due to the presence of both the impurity atom and the C-AF order, all the bonds are now inequivalent. The low and highangle phases shown in Fig. 9 differ in a rather subtle way.
Looking at the spin-spin correlations in Fig. 10 as a function of the π hybridization, we can clearly see the discontinuities across the transition. The S α i ·S α j correlations along the host bonds exhibit only small jumps and do not change much whereas at the impurity bonds the changes are more notable, see Figs. 10(a-d). Specifically, the behavior of the S y i S y j correlator around the impurity is interesting-it changes signs across the transition. This change in sign is also visible in the non-collinear S i × S j z correlator which grows in absolute values as we increase V pdπ towards the transition just to drop to zero and stay at zero in the entire zero-angle phase.
As discussed in Sec. IV, we point out that for the doped configuration the thermal profile of the bond angle is not qualitatively affected by a change of the Hund's coupling and spin-orbit coupling (see Appendix C). Moreover, the outcomes of the single effective d 3 − d 4 bond within the spin-orbital model are qualitatively consistent with those of the full multi-orbital Hubbard model (as reported in Appendix D).

V. SUMMARY AND CONCLUSIONS
The phenomenon of NTE effect deals with the increase of lattice distances and, in turn, of the volume of the unit cell, if the material is thermally cooled. Here, we studied the conditions for obtaining negative thermal effects in materials characterized by strong electron-electron correlations involving spin degrees of freedom and multiple orbitals. In this context, when considering for instance perovskite oxides with transition metal elements inside the oxygen octahedral cages, the change of the unit cell size is linked to the octahedral rotations and thus to the amplitude of the TM-O-TM bond angle. An undistorted (distorted) TM-O-TM bond with zero (non-vanishing) bond angle corresponds to large (small) volume of the unit cell, respectively. Hence, the TM-O-TM bond angle is a key parameter in the phenomenon of NTE. The TM-O-TM bond angle, on the other hand, sets out the electrons connectivity within the crystal lattice through the hybridization of the orbitals at the transition metal and oxygen atoms. Hence, the amplitude of the electronic kinetic energy is determined by the TM-O-TM bond angle. For the case of an electron correlated insulating phase with multi-orbital configurations, the Coulomb interaction prevents the charge motion and thus the kinetic energy is substantially given by virtual charge transfer pro-cesses and by the resulting spin-orbital exchanges. Then, while for a single orbital model one would expect that the kinetic energy is optimal by having an undistorted bond, the resulting configuration in the case of multi-orbital electronic state cannot be a priori predicted and requires the knowledge of the orbital dependent magnetic correlations. This aspect underlines the complexity of the fundamental problem we have been facing in this paper: to establish how the spin-orbital correlations determine the TM-O-TM distortions and single out the conditions for achieving NTE effects. On this basis, it is evident that one needs to evaluate the dependence of the spin-orbital exchange on the bond angle. For this reason, we employed a theoretical approach that explicitly includes the bond angle in the amplitude of the d − p hybridization processes. In this way, we have accounted for the modification of the free energy as a function of the TM-O-TM bond configuration and assess the role of the electron correlations in a multi-orbital electron system. From our study we have been able to identify few mechanisms which are crucial for achieving NTE effects. First, the anisotropy of the orbital correlations, as pointed out above, can play an important role. In particular, the orbital anisotropy or preferential orbital occupation, which is driven by the electron correlations, leads to a ground state that is prone to exhibit NTE effects. The occurrence of these anisotropic orbital correlations drives the bond direction to get less distorted bond. This is what one can find at low temperatures where the spin-orbital correlations are stronger. On the other hand, thermal fluctuations tend to weaken the spin-orbital correlations, a condition which favors more distorted bonds. The thermal disorder also manifests in the spin channel with the formation of non-colinear magnetic correlations that energetically support a distorted bond. If such conditions can be realized, then the optimal bond angle will decrease by increasing temperature and thus we will have resulting NTE effects. A second key element that emerges from our analysis is related to the presence of low point group symmetry conditions. In the regime of strong Coulomb interaction, we find that the breaking of the C 4 rotation symmetry can turn a positive thermal expansion effect into a NTE behavior when tracking the TM-O-TM bond angle. In our study, this effect has been studied by considering a magnetic pattern that breaks the C 4 rotation symmetry. However, we expect that similar outcomes for the NTE can be also obtained by considering orbital or electronic patterns that reduce the point group symmetry of the system. In this respect, a significant modification of the NTE effects can be foreseen in thin films under strains or by applying electric field which would naturally provide microscopic conditions with low point group symmetry. Another relevant aspect concerns the role played by the strength of Hund's coupling to induce the occurrence of NTE effects. The solution of the multi-orbital Hubbard model for the single bond allows us to explore a regime of strong Hund's interaction, for either weak or strong intra-orbital Hubbard U , since it takes into account at the same level the charge, spin and orbital degrees of freedom. Indeed, the amplitude of Hund's coupling can play an important role in setting out the thermal negative expansion effects. This is evident, for instance, in the case of d 4 −O −d 4 bond where one finds a changeover of the optimal angle with a positive thermal effect profile that turns into a negative thermal one by increasing the amplitude of Hund's coupling. Additionally, we also find an enhancement of the NTE effects, with a larger negative variation of the optimal angle, and a shift of its onset towards higher temperatures, for the doped d 3 − O − d 4 bond configuration. Since this type of changeover from positive to NTE effects is not observed when considering the large U limit, as evaluated within the low-energy model employed for the simulation of the square cluster, we argue that the Hund's coupling in the intermediate regime of U interaction is relevant for the occurrence of negative thermal effects in multi-obital materials. Although our model analysis has not been expanded to investigate a Hund's exchange driven correlated metal [42], on the basis of the results of the single bond for the multiorbital Hubbard model, we argue that the effect of large Hund's coupling can favor the occurrence of NTE. In this context, we are confident that our analysis may further stimulate the search for NTE effects in doped Hund's metal exhibiting strong correlations, as for instance in iron based materials [43,44], or in the proximity to a metal-to-insulator orbital selective transition. Concerning real materials where our study could be applicable, we recall that the analysis has been substantially motivated by the NTE observed in the Ca 2 RuO 4 compound. In this context, our analysis is able to account for the occurrence of NTE effects in the Ca 2 RuO 4 family as due to the transition metal substitution. The analysis specifically refers to the substitution of Ru (d 4 ) with Mn (d 3 ). We find that both in the intermediate and in the strong coupling regimes we can achieve NTE effects close to the impurities. As mentioned above, since the impurity can act as an orbital polarizer it can sustain less bond distortions at low temperatures. Finally, we mention that the inclusion of a phenomenological electron-lattice potential, as discussed in Appendix E, allows to regularize the maximal bond distortions that are obtained within a purely electronic description and to make smoother the transitions from large to small bond angles across the phase space spanned by the p − d hybridization hoppings. This is physically plausible and indicates how one can simulate the conditions which might be more realistic for material applications. The analysis of correlated models that microscopically include the electron-phonon interaction will be the subject of future studies. The perturbative expansion involving direct p − d exchange, as well as Anderson and Goodenough d − p − d superexchange, yields six different virtual excitation energies given by −q −1 i expressed as: where U and J H are Hubbard and Hund's couplings at the metal ions and p is on-site energy of the oxygen orbital. The hoppings between metal and oxygen can be described by the matrices where T M1−O dp describes hopping between the first TM ion and oxygen and T M2−O dp between second transition metal ion and oxygen in presence of in-plane octahedral rotation. The rows of these matrices refer to the {d yz , d zx , d xy } orbital states of the TM ion and columns to {p x , p y , p z } orbital states of the oxygen. The hopping amplitudes are given by Slater-Koster rules, i.e., where θ is the angle of the bond between metal and oxygen with respect to 100 direction in the (001) plane.
To derive exchange we use a second-order perturbative expansion. In the atomic ground state we have two spins S = 1 and orbitals L = 1 at every metal ion. Oxygen is fully occupied and yields no spin or orbital degree of freedom. In the second order perturbation expansion with respect to the hopping, we get the d 4 1(2) p 6 d 5 1(2) p 5 exchange is given by the Hamiltonian, Here the coefficients in Eq. (A4) take the following form, β In the fourth order perturbation expansion with respect to the hopping we have, inter alia, the Anderson processes where the electron from the second metal goes through oxygen to the first one and back again, i.e.,    the coupling constant is proportional to the optimal angle that is obtained by minimizing the free energy of the purely electronic system. Such an assumption is physically plausible because the lattice acts as a restoring force and it is reasonable to expect that its strength, through the electron-phonon coupling, can be dependent on the amplitude of the optimal angle that is set by the electronic degrees of freedom. The proposed phenomenological form of the total free energy can be hence expressed as: where k is an elastic constant. Therefore, we assume that the lattice deformation potential is harmonic with the stiffness being proportional to θ el , the angle that minimizes the electronic part F el (θ).
We have implemented the phenomenological model in Eq. (E1) for the case of the d 4 plaquette. In Fig. 14(a) ��Vpd�� �meV� Results for the undoped d 4 plaquette with C-AF magnetic order (see the inset) and with JH and λ increased by a factor of χ = 3.2 compared to Sec. III: (a,b) typical curves of the optimal bond angle and the ground state averages of S1× S2 z and S1× S3 z versus δ(V pdπ ) = V pdπ −V0 for V pdσ = 1.5 eV in the intermediate phase between θ = 0 and θ = π/4; (c,d) thermal dependencies of the optimal bond angle and average S1 × S2 z and S1 × S3 z in the intermediate phase for we show the typical curves of the electronic free energy at zero temperature as a function of the bond angle. We note that the curves have a single minimum that moves from 0 to π/4 in a small window of energy variation for the V pdπ parameter (of the order of µ eV). In Fig. 14(b) we present the optimal angle versus V pdπ and versus temperature (inset) for the electron-phonon model of Eq. Results for a pure d 4 plaquette as obtained within the effective low-energy spin-orbital exchange approach: (a) electronic free energy versus bond angle for a given V pdσ = 1.5 eV and different values of V pdπ = V0 +δ(V pdπ ) chosen in the intermediate-angle phase, (b) optimal angle versus V pdπ for a purely electronic free energy (red curve) and for the phenomenological eletronic-phononic free energy [Eq. (E1)] assuming a representative value for the coupling constant, i.e., k = 6 meV. The inset of (b) shows the optimal angle versus temperature dependence for V pdπ = 1.2902 eV. The other electronic parameters are: U = 8.0, p = −4.5, δ = 0.35, δort = 0.09, and λ = 0.075, all values are in eV.
(E1). We notice that by a proper choice of the constant coupling k, which should be of the same order as the electronic free energy change between 0 and π/4 away from the transition region, one can stabilize an intermediateangle phase at a scale of up to tens of eV in V pdπ . Such a phase is demonstrated to exhibit a NTE effect at finite temperature, as it can be seen from an inspection of the inset in Fig. 14(b)).