Evolution of Spin-Orbital Entanglement with Increasing Ising Spin-Orbit Coupling

Several realistic spin-orbital models for transition metal oxides go beyond the classical expectations and could be understood only by employing the quantum entanglement. Experiments on these materials confirm that spin-orbital entanglement has measurable consequences. Here, we capture the essential features of spin-orbital entanglement in complex quantum matter utilizing 1D spin-orbital model which accommodates SU(2)xSU(2) symmetric Kugel-Khomskii superexchange as well as the Ising on-site spin-orbit coupling. Building on the results obtained for full and effective models in the regime of strong spin-orbit coupling, we address the question whether the entanglement found on superexchange bonds always increases when the Ising spin-orbit coupling is added. We show that (i) quantum entanglement is amplified by strong spin-orbit coupling and, surprisingly, (ii) almost classical disentangled states are possible. We complete the latter case by analyzing how the entanglement existing for intermediate values of spin-orbit coupling can disappear for higher values of this coupling.


I. INTRODUCTION
Quantum complex matter is characterized by several intertwined degrees of freedom. Mott insulators with orbital degeneracy belong to this class of materials and are typically characterized by entangled spin-orbital states, either only in the excited states or in the ground state as well. Here we shall present the spin-orbital entanglement in a broader context, summarizing important results that have been accumulated over the years. With this background established, we will focus on the recent study [1] of the evolution of spin-orbital entanglement with increasing Ising spin-orbit coupling, enriching this story with further details.
By definition, quantum entanglement makes information available in one part of the quantum state dependable on the result of the measurement on the other part of the state [2]. A particular kind of quantum entanglement, involving quantum spin and orbital operators in the entire quantum system [3] rather than entanglement between two spatial regions, may appear in the strongly correlated regime of oxides with orbital degeneracy. The importance of orbitals, besides spins, in the effective description of such systems was pointed out long ago by Kugel and Khomskii [4]. After this pioneering work it has been realized that orbital operators are quantum in nature, enforcing equal treatment with electron spins and producing joint spin-orbital fluctuations [5], both in the ground [6,7] and in excited [8,9] states. Quantum fluctuations are further amplified by the spin-orbit coupling [10].
In theory, spin-orbital entanglement has been considered in dimerized phases, spin-orbital models with S = 1 2 spins, both in perovskites [11] and on triangular lattices [12][13][14], as well as exactly solvable models such as SU(2)⊗XY model [15], but is weakened in layered oxides suggested as a possible realization of superconductivity [16]. Several experimental observations in strongly correlated quantum matter cannot be understood however without spin-orbital entanglement. For instance, the low energy spectral weight in the optical spectroscopy for LaVO 3 can be explained only when full quantum coupling between spins and orbitals is included [17,18], in contrast to LaMnO 3 where orbitals disentangle from large S = 2 spins [19]. Entangled states are also found in disordered spin-orbital liquid states [20,21] and in Kitaev quantum liquids [22][23][24][25], where spatial long-range entanglement [26] is induced due to strong on-site spin-orbit coupling.
Entanglement in the vanadium perovskites has more manifestations. One of them are the magnon excitations in the C-type AF (C-AF) phase of YVO 3 , where the gap opens along the Γ-Z direction [27], indicating the dimerization on the bonds along the c axis which occurs jointly both in spin and orbital channel at finite temperature [28]. Another experimental confirmation of spin-orbital entanglement is the evolution of the Néel temperature in the C-AF phase with decreasing ionic radius of R ions in RVO 3 which is induced by the coupling to the orbital state [29]. Experimentally the entanglement is observed in the critical competition between the two spin-orbital ordered states [30] and was investigated recently in Y 0.70 La 0.30 VO 3 single crystals [31]. Finally, entangled states play a role in the excited states [32] and may help to identify quantum phase transitions in spin-orbital models [8], including the transitions at finite magnetic field [33].
Entangled states play also a role in doped spin-orbital systems. One example is colossal magnetoresistance in manganites, where the orbital polarons modify the ground state in a one-dimensional (1D) model [34]. Another example is vanadium perovskites with charged defects which influence occupied orbitals and lead to the collapse of orbital order [35]. The spin-orbital entangled states were also found among the surface states of topological insulators [36]. More examples of spin-orbital entanglement could be found in review articles [37,38].
As a matter of fact, both superexchange (i.e., intersite) [3] and spin-orbit coupling (on-site) terms can promote spin-orbital entanglement, where it may lead to an effective description in terms of the spin-orbital pseudospins with exotic interactions. The best example is the onset of the Kitaev-like physics in the iridium oxides [10] in the so-called relativistic Mott insulators [39]. Possible competition between these two mechanisms (superexchange bonds and spinorbit) or enhancement of entanglement by their joint action is a challenging problem. It adds another aspect to the complexity of quantum matter in high-T c superconductors [40][41][42][43][44], intimately connected with quantum magnetism in strongly correlated transition metal oxides [45]. Here we analyze the spin-orbital entanglement evolution on a minimal available model, consisting of SU(2)⊗SU(2) symmetric superexchange term [46] and Ising on-site spin-orbit coupling in one dimension. Below we recall the recently published thorough numerical and analytical study [1], supplementing it with the intermediate global stages and previously not considered cuts, thus completing the story about the disappearance of spin-orbital entanglement for certain parameter regime.

II. THE MODEL AND BASIC QUESTIONS
The full 1D spin-orbital Hamiltonian considered here has the form [1], The superexchange H SE is a 1D model with two spin (S = 1 2 ) and two orbital (pseudospin) (T = 1 2 ) degrees of freedom per site in a Mott-insulating regime described by SU(2)⊗SU(2) symmetric expression [47], where J ≡ 1 is the energy unit and the last term αβ cancels a constant. Spin operators in the scalar product are, {S i } ≡ {S x i , S y i , S z i }, and similar for pseudospins, The superexchange depends on two parameters, α and β, and the spin/pseudospin exchange could be of either sign. Hence, also spin [pseudospin] order could be of either sign, ferromagnetic (FM) or antiferromagnetic (AF) [ferro-orbital (FO) or alternating-orbital (AO)]. The competition between different channels of exchange gives a rather rich phase diagram of Eq. (2) [47], with three product phases, FM⊗AO, FM⊗FO, and AF⊗FO, see Fig. 1. Yet, the entanglement plays a crucial role and the phase diagram includes the fourth spin-orbitally entangled phase near the exactly solvable SU(4) point [47], where the entanglement entropy of the von Neumann type is finite [48]. Also finite (but much smaller) spin-orbital entanglement characterizes the parameter range where α and β are large positive, promoting alternating arrangement of spins (pseudospins) with accompanying quantum fluctuations [48].
The above model (2) can be derived from the microscopic degenerate Hubbard model with the diagonal-only hopping t and the Coulomb repulsion U -then J = 4t 2 /U and α = β = 1 4 [49]. The highly symmetric form of spin and pseudospin interactions follows from the degeneracy of singlet and triplet excited states. These conditions are not fulfilled in widely studied geometries with corner-sharing MO 6 octahedra (M is a transition metal atom) as in perovskites, but the model is a realistic Kugel-Khomskii model for 1D chains of face-sharing MO 6 octahedra [49]. Such materials include for example hexagonal crystals such as BaCoO 3 , BaVS 3 , or CsCuCl 3 . Recently it was shown that the specific heat and the susceptibility show anomalies at finite temperature which can be attributed to phase transitions even in the regime of spin-orbital liquid [50]. Here we treat it as a generic model which stands for spin-orbital superexchange and entanglement.
As the competing interaction in Eq. (1) we take Ising on-site spin-orbit coupling [1], It has Z 2 -symmetry and favors the opposite orientation of spin S z i and orbital T z i moments at site i. This reduced type of spin-orbit coupling occurs, for instance, when orbital order concerns two t 2g orbitals and the third one is inactive as in LaVO 3 or YVO 3 [51]. Then, for an idealized case of untilted oxygen octahedra the spin-orbit coupling takes the form 2λS x i T y i . Surprisingly, similar 2λS z i T y i form can also be derived for KO 2 , where peculiar form of superexchange includes p orbitals of molecular oxygen [52]. Taking the symmetry of Hamiltonian (1) these forms are equivalent since they would modify the quantum states in the same way as the formula in Eq. (3), up to the rotation in the spin/pseudospin sectors. However, in case of finite orbital-lattice interactions this should be carefully studied.
For completeness, we add that when the full spin-orbit coupling is considered [10], effective entangled pseudospiñ J operators replace the sum of spin S and pseudospin T-theJ multiplets then define the ground state and excited states of the system. For instance, spin-orbital entanglement gives the breakdown of singlets and triplets at the Fermi surface in Sr 2 RuO 4 [53][54][55] and their linear combinations decide about the orbital character of the pairing components [56]. We believe that future work will also address the role played by entanglement for the d 4 ions in more detail where the pseudospin could vanish and the magnetic moments are hidden in the ground state but would become active in excited states [57][58][59][60][61].

A. Entanglement entropy and correlation functions
In Ref. [1] we have verified that for (perturbed) product phases of the SU(2)⊗SU(2) model (see Fig. 1), the spinorbital entanglement evolution with increasing λ experiences only marginal finite-size effects. Here we want to show the qualitative effects caused by large spin-orbit coupling and for that the presented studies of small systems suffice. Therefore, we restrict ourselves to the exact diagonalization (ED) of a 4-site chain (L = 4) with periodic boundary conditions (PBC). We examine the spin-orbital entanglement via entanglement entropy of the von Neumann type per site, where is the reduced density matrix with the orbital degrees of freedom traced out, estimating the entanglement by verifying how mixed the reduced state is. Independently, we quantify spin-orbital entanglement via correlation on a bond i, i + 1 between nearest neighbors: At finite λ, we monitor the response of the system to the Ising spin-orbit coupling via suitable on-site correlation, Finally, the spin and orbital states are investigated via spin-only and orbital-only correlation, They allow to identify the product phases, FM⊗AO, FM⊗FO, and AF⊗FO in a straightforward way. The entanglement entropy was investigated previously in the SU(2)⊗SU(2) model at λ = 0, and was found to be finite beyond the product phases where at least one (spin or orbital) component has suppressed quantum fluctuations in a ferro-state [48]. We have investigated the level of spin-orbital entanglement in the ground state in the presence of finite spin-orbit coupling and identified small and large λ regimes along the special α + β = 0 line in the model (1) [1]. While small spin-orbit terms ∝ λ can be treated in the perturbation theory, below we focus on the more interesting regime of large λ.

B. Effective Hamiltonian in the large λ limit
In Ref. [1], following the general scheme set in Ref. [10], we have described the low-energy physics in the λ → ∞ limit by the model utilizing effective pseudospins {J i }. Namely, the Hamiltonian (2) projected onto the low-energy part of H SOC eigenbasis takes the XXZ form, It is remarkable that the Ising part of the above Hamiltonian, 2(α + β)J For this state one finds that the mean field approximation for the spin-spin correlations is exact and C SO = 0 (6), i.e., spin-orbital entanglement vanishes.
In Fig. 2(c) the highly entangled state exists in a close neighborhood of the β + α = 0 line. More precisely, highly entangled 1D spin-liquid class of states, called also 1D quantum antiferromagnet in the literature, could be found in a "quantum stripe" with one sharp boundary at α + β = − 1 2 and system size-dependent soft boundary seemingly approaching AF Heisenberg line α + β = 1 2 for larger chains, thus roughly fitting within the two outer white dashed lines. We observe that the entanglement is maximal between the α + β = − 1 2 and α + β = 0 lines, and decays at a rate depending on the system size. This behavior is verified for finite systems, i.e., for rings of L = 4 sites (shown here) and of L = 8, 12 sites (studied before in [1]). From α + β = 1 2 upwards we find interpolation between the 1D Heisenberg antiferromagnet and the classical Néel state in the limit of α + β → ∞ where the quantum fluctuations are quenched.

IV. RESULTS AND DISCUSSION
A. Gradual evolution of spin-orbital entanglement with increasing Ising spin-orbit coupling Here we first complement the analysis of Ref. [1] with the data for the entanglement entropy in the intermediate regime of λ/J, building a global picture of spin-orbital entanglement under increasing λ. Namely, we add 6 entanglement entropy scans in Figs. 2(d)-2(i), fitting in between the panels (b) and (c). Results presented here for L = 4 include small finite size effects and thus one finds finite spin-orbital entanglement above the stripe for large λ. Nevertheless, the evolution of perturbed product phases is qualitatively the same as before for all 4n chains [1], showing that spin-orbital entanglement is pretty local.
Finite spin-orbital entanglement in the central region of λ = 0 slice [ Fig. 2(a)] increases and spreads out with growing λ, see panels (b) and (d)-(g). In particular, infinitesimal values of entanglement entropy appear on the boundaries of perturbed product phases containing either AF or AO component, see panels (b), (d), and clearly visible on panels (e)-(f )). At the same time, FM⊗FO phase does not show any entanglement, which cannot be induced in this classical state by also classical spin-orbit Ising coupling (3). Interestingly, the border between completely classical region, originally hosting FM⊗FO ground state, and the spin-orbitally entangled region gradually changes its shape. The change starts for infinitesimal λ > 0 but can be resolved starting from panel (d). A sharp boundary, at first restricting the classical state to the third quarter of panel (b), moves so that a much wider classical region is visible in panel (i). This result is quite close to the classical triangle seen in panel (c).
There is a huge difference in size between the maximal and minimal spin-orbital entanglement for varying λ, displayed in Figs. 2(a)-2(b), and 2(c), correspondingly. At the λ = 0 case the spin-orbital term, (S i S i+1 )(T i T i+1 ), can produce substantial entanglement in the ground state only in a small approximately triangular region, where β and α scales of purely spin or purely orbital interactions are small enough. This corresponds to the critical spinorbitally entangled region in the phase diagram of the SU(2)⊗SU(2) model [47]. In λ → ∞ limit, the "quantum stripe" characterized by the highly entangled state can be also related to the (S i S i+1 )(T i T i+1 ) joint interaction. Namely, when projecting the Hamiltonian (2) on the lower part of H SOC eigenbasis term by term, one observes that (S i S i+1 )(T i T i+1 ) reduces to the XY term, independent of α and β, while β S i S i+1 and α T i T i+1 transform into ∝ βJ z iJ z i+1 and ∝ αJ z iJ z i+1 , respectively. As a result, on the α + β = 0 line the Ising terms cancel out, leaving the XY model with entangled ground state, and on the entire "quantum stripe" the ratio Ising/XY supports the highly entangled state. The left boundary of the stripe is well defined by the α + β = − 1 2 line where FM Heisenberg ground state appears. The degenerate energy splits immediately below the above line, leaving Ising 2-fold degeneracy.
B. Entanglement evolution for the quantum phases along the α + β = 0 line Next, we recall the evolution of spin-orbital entanglement on the α + β = 0 line, marked by the middle diagonal white dashed line in different panels of Fig. 2 and investigated originally in Ref. [1] (see Fig. 2 there). This line goes through the two product phases, FM⊗AO and AF⊗FO, and through the entangled phase which separates them (similar to the white line in Fig. 2 at α + β = 1 2 ). Two distinct types of evolution emerged on the vertical lines in Fig. 3 (under increasing λ): (i) evolution of type A starting from AF⊗FO or FM⊗AO state with no spin-orbital entanglement at λ = 0, and (ii) evolution of type B starting from the spin-orbitally entangled critical SU(4) singlet ground state [48]. Altogether, the behavior of the entanglement entropy and the correlation functions displayed in Fig. 3 for the L = 4 system is very similar to that shown in Fig. 2 of Ref. [1] for the L = 12 system which confirms that especially for the evolution of type A the spin-orbital entanglement is a local phenomenon and the finite size effects are small in the quantities analyzed here.
In case (i), the ground state has large degeneracy at λ = 0 due to FM or FO sector. From the energetic point of view the level crossing is at λ = 0 in the degenerate state. Then λ = 0 lifts the degeneracy and leads to non-degenerate perturbed product states. From that on, all changes in the ground state have a crossover character which is visible in the smooth increase of the entanglement entropy S vN along the line A, see Fig. 3(d). Finite entanglement is then transferred to the product phases at increasing λ. But the most important is that even when the energy scales fully separate at λ/J 1, finite spin-orbital entanglement is found on superexchange bonds, see Fig. 3(b). The finite value of C SO (7) is induced by the (S i · S i+1 + α) (T i · T i+1 + β) term to be reduced to the XY model and freed from α and β-dependent Ising contributions in λ → ∞ limit. Figure 3 shows that this limit is approximately valid as soon as the non-zero entropy and C SO values spread on α + β = 0 line. This can be also verified by the maximal response in on-site correlations O SO (7), see Fig. 3(c).
Entanglement entropy S vN in case (ii) starts from a finite value at λ → 0 in spin-orbitally entangled critical singlet phase [48]. Such singlet correlations are robust along the path B and thus the spin-orbital entanglement (described by the entropy and the nearest neighbor correlations C SO ) can change only in a discontinuous way, by the kink shown in Fig. 3(e). We have seen before that the evolution of spin-orbital entanglement for 4n-site chains occurs along the line of type B through a series of n kinks [1]. Entanglement gradually gets more robust up to λ crit /J 0.2, above which the highly-entangled state is reached. With this, we completed a short summary of the main results of Ref. [1]. C. Entanglement evolution for the quantum phases along the α + β = −1 line Figure 2 revealed the gradual evolution of the spin-orbital entanglement in λ. Moreover, we have seen how classical region in both spin and orbital sectors, initially restricted to FM⊗AO phase in the third quarter of the diagram, gradually spreads out, reaching extended triangular shape within the boundaries shown in Fig. 2(c). Within this triangular region, the evolution of FM⊗FO phase with increasing λ is rather trivial as long as we have the Ising form of the spin-orbit coupling (3). However, it is not obvious how the lower parts of AF⊗FO or FM⊗AO phases, see Fig.  1, approach the Ising FM ground state below the α + β = − 1 2 line (the lowest white dashed line in Fig. 2). Here, we concentrate on the evolution in λ along the representative α + β = −1 line (red dashed line in Fig. 2), including both of the above cases, completing the review of the λ evolution for the Hamiltonian Eq. (1).
We have selected this line because it crosses all three product phases AF⊗FO, FM⊗FO, and FM⊗AO at λ = 0 ( Fig. 1) and lies below the "quantum stripe" for large λ. The line runs close enough to the center of the {α, β} maps to bear quite insightful information. The entanglement evolution in λ along the line is shown in Fig. 4. Panels (a) and (b) show, respectively, the entanglement entropy and the C SO maps. One can observe that AF⊗FO and FM⊗AO, represented again by cross section C, at first gain gradually spin-orbital entanglement with growing λ, just like in the evolution along the line A in Fig. 4. S and T correlations confirm that those are indeed perturbed AF⊗FO and FM⊗AO phases. [S correlation is shown in Fig. 4(f ).] However, for some particular value of λ the growth of entanglement is sharply cut off. For parallel lines lying deeper within FM⊗FO phase, the effects discussed here will appear as well but for higher values of λ (see the trends in Fig. 2).
The origin of this evolution can be resolved by inspecting the energy spectrum of the L = 4-site ring, plotted in Fig.  4(d) at α = −2.5 and β = 1.5 and varying λ. Here we note that for λ > 0 the ground state degeneracy is completely lifted. It is quite difficult to spot the level crossing suggested by sharp changes in entanglement. Therefore, in Fig.  4(g) we analyzed the energy differences between the excited energies and the ground state energy, ∆E = E i − E 0 , plotted in gray. For convenience the ground state energy is renormalized as ∆E = E i − E 0 = 0 and plotted in blue. One observes that at λ ≈ 8.175 ± 0.025, 2-fold degenerate excited level crosses with the non-degenerate ground state, leaving 2-fold degenerate ground state. In this way the perturbed product state becomes an excited state and gradually gains entanglement, even crossing another excited 2-fold degenerate state for further increasing λ, while the Ising ferromagnet becomes the ground state.
For completeness, we describe also the evolution of the FM⊗FO state, choosing cross section B at α = β = − 1 2 , see Fig. 4(h). The degenerate 25-fold FM⊗FO level splits immediately when λ > 0, leaving only 2-fold degeneracy. Consistently, panels (a) and (b) confirm no spin-orbital entanglement, while panels (c) and (h) show that the response of the system to spin-orbit coupling is instantaneous.

V. SUMMARY AND CONCLUSIONS
In summary, both in Ref. [1] and here, we have accumulated evidence for two distinct regimes of Ising spin-orbit coupling: (i) when it is weak, i.e., λ < λ crit , one is in the perturbative regime and the superexchange is almost undisturbed, while (ii) strong spin-orbit coupling λ > λ crit reduces entanglement due to superexchange but may induce it in the product phases where quantum fluctuations are finite at λ = 0, in either spin or orbital channel. The fluctuations enhance strongly the entanglement of already spin-orbital entangled central region. Note that here we define λ crit differently for the evolution of type A and B in Fig. 3 and type C and D in Fig. 4.
For the evolution of type A on the α+β = 0 line, λ crit stands for the dynamical region of rapid entanglement growth (in logarithmic scale), see Fig. 3. It connects the perturbed FM⊗AO and AF⊗FO phases to the highly entangled state at λ > λ crit via crossover [1]. The evolution of type C on the α + β = −1 line is more subtle. The perturbed FM⊗AO (AF⊗FO) rapidly gains entanglement in a smooth way but the crossover process is interrupted by a level crossing with the Ising ferromagnet. These transitions seem to be immune to finite-size effects.
For the evolution of type B on the α + β = 0 line, analyzing the spin-orbital entangled critical SU(4) singlet state through n kinks for L = 4n-site system, we can set λ crit (L → ∞) ≈ 0.2J [1]. In other words, finite-size scaling analysis in Ref. [1] suggests that the entanglement growth from the singlet state will be continuous in the thermodynamic limit, with a phase transition at λ crit ≈ 0.2J. Formally λ crit could be set to 0, since the large FM (FO) ground state degeneracy at λ = 0 is immediately reduced to 2-fold degeneracy at λ = 0, matching the Ising state. However, the simplest possible state representing the ground state energy is disentangled in the whole λ range. The case of FM⊗FO is rather special as the quantum fluctuations are suppressed in both spin and orbital channel at λ = 0. They are also absent in the Ising spin-orbit coupling (3), so irrespectively of the actual value of λ, there is no driving force to generate entangled states. As a result, this phase remains disentangled and has just maximal on-site spin-orbital correlations O SO = − 1 4 which are unable to induce any intersite entanglement. We emphasize that this will change when spin-orbit coupling is quantum and many-body states fully develop. For instance, such states are found in two-dimensional (2D) iridates Ba 2 IrO4 and Sr 2 IrO 4 , where it was demonstrated that electron-and hole-doped iridates are fundamentally different [62]. A similar situation is found in CuAl 2 O 4 spinel, where spin-orbital entangled Kramers doublets form and again the electron addition spectra are well described by the J ef f = 1 2 hole state, whereas electron-removal spectra have a rich multiplet structure [63]. Although the low-energy physics depends on the iridate geometry, it would be challenging to investigate in a similar way the interplay of spin-orbital superexchange with quantum spin-orbit coupling.
We also reported an interesting "emergence" of the spin-orbital entanglement along the α + β = −1 line in the Ising limit. In this case, both for the negligibly small λ as well as for the large λ > λ crit the spin-orbital entanglement vanishes. Surprisingly, however, for a moderate value of λ (and a finite value of α) the spin-orbital entanglement becomes finite. This shows the intricate complexity of the studied system. Crucially, the analysis at λ > λ crit is consistent with the phase diagram of the effective XXZ Hamiltonian (10). In particular, the ground state with vanishing spin-orbital entanglement appears when the quantum fluctuations vanish in the ground state of an effective model. This is the case of an Ising ferromagnet obtained at α + β ≤ − 1 2 . Finally, we argue that the results presented here may play an important role (after extension) in the understanding of the strongly correlated systems with non-negligible spin-orbit coupling-such as, e.g. the 5d iridates, 4d ruthenates, 3d vanadates, the 2p alkali hyperoxides, and other to-be-synthesized materials. More discussion of the entanglement in these materials was given before [1]. As discussed also in Ref. [1], the results presented for the 1D model can be (on a qualitative level) extended to the higher-dimensional hypercubic lattice. This is due to the fact that the mapping to the effective model is valid independently of the dimension. Therefore, the spin-orbital correlation function (and consequently the spin-orbital entanglement) is nonzero for the similar range of model parameters in 2D or threedimensional systems as for the 1D model reported here.