How spin-orbital entanglement depends on the spin-orbit coupling in a Mott insulator

The concept of the entanglement between spin and orbital degrees of freedom plays a crucial role in our understanding of various phases and exotic ground states in a broad class of materials, including orbitally ordered materials and spin liquids. We investigate how the spin-orbital entanglement in a Mott insulator depends on the value of the spin-orbit coupling of the relativistic origin. To this end, we numerically diagonalize a one-dimensional spin-orbital model with Kugel-Khomskii exchange interactions between spins and orbitals on different sites supplemented by the on-site spin-orbit coupling. In the regime of small spin-orbit coupling with regard to the spin-orbital exchange, the ground state to a large extent resembles the one obtained in the limit of vanishing spin-orbit coupling. On the other hand, for large spin-orbit coupling the ground state can, depending on the model parameters, either still show negligible spin-orbital entanglement or evolve to a highly spin-orbitally-entangled phase with completely distinct properties that are described by an effective XXZ model. The presented results suggest that (i) the spin-orbital entanglement may be induced by large on-site spin-orbit coupling, as found in the 5 d transition metal oxides, such as the iridates; (ii) for Mott insulators with weak spin-orbit coupling of Ising type, such as, e.g., the alkali hyperoxides, the effects of the spin-orbit coupling on the ground state can, in the ﬁrst order of perturbation theory, be neglected.

ground states of LaTiO 3 , LaVO 3 , YVO 3 , and Ba 3 CuSb 2 O 9 [17][18][19][20][21][22], where spin and orbital degrees of freedom are intertwined and entangled. A similar situation occurs in MnP [23], the first Mn-based unconventional superconductor under pressure, and in some two-dimensional (2D) model systems [24][25][26]. In all these cases, the ground state can only be explained by invoking the joint spin-orbital fluctuations. Consequently, the mean field decoupling separating interactions into spin and orbital degrees of freedom fails and cannot be used.
In this paper, we study the spin-orbital entanglement which manifests itself when a quantum many-body system with interacting spin and orbital degrees of freedom is split into the subsystems with separated degrees of freedom; i.e., one attempts to write interacting spin and orbital wave functions separately [27][28][29]. The concept of entanglement has been first introduced in these systems to understand the violation of the so-called Goodenough-Kanamori rules [30,31] in the ground states of several transition metal oxides with partially filled 3d orbitals, strong intersite spin-orbital (super)exchange interactions but typically negligible value of the on-site spinorbit coupling. It was also realized that entanglement is important to understand the excited states where spin and orbital variables are intertwined. Good examples are the temperature evolution of the low-energy excitations in LaVO 3 [32,33] and the renormalization of spin waves by orbital tuning in spin-orbital systems due to the weak interactions with the lattice [34]. Furthermore, the spin-orbital entanglement is also crucial to understand the first unambiguous observations of the collective orbital excitation (orbiton) in Sr 2 CuO 3 and CaCu 2 O 3 [35,36] and their interpretation in terms of the spin and orbital separation in a 1D chain [37][38][39].
Crucially, the spin-orbital entanglement is expected as well in the oxides with strong on-site spin-orbit couplingprobably best exemplified by the partially filled 4d and 5d orbitals as found in the ruthenium and iridium oxides [40,41] or in the recently discovered 5d Ta chlorides [42]. For instance, the concept of spin-orbital entanglement was recently invoked to understand the inelastic x-ray spectrum of Sr 3 NiIrO 6 [43], the ground state of H 3 LiIr 2 O 6 [44], and, in a different physical setting, the photoemission spectra of Sr 2 RuO 4 [45,46]. Interestingly, the peculiarities of the interplay between the strong on-site spin-orbit coupling and the spin-orbital (super)exchange interactions allowed for the onset of several relatively exotic phenomena in this class of compounds-such as a condensed matter analog of a Higgs boson in Ca 2 RuO 4 [47,48] or the strongly directional, Kitaev-like, interactions between the low-energy degrees of freedom (pseudospins) in some of the iridates or ruthenates on a quasi-2D honeycomb lattice (Na 2 IrO 3 , Li 2 IrO 3 , α-RuCl 3 , H 3 LiIr 2 O 6 ) [44,49,50]. The latter might be described to some extent by the exactly solvable Kitaev model on the honeycomb lattice which, inter alia, supports the onset of a novel spin-liquid ground state with fractionalized Majorana excitations [51].

C. Main question(s) and organization of the paper
It may come as a surprise that the concept of the spinorbital entanglement has so far been rigorously investigated only for the systems where the spins and orbitals at neigh-boring sites interact, as a result of the spin, orbital, and spin-orbital (super)exchange processes in Mott insulators [27][28][29][52][53][54][55][56][57]. This case is physically relevant to all Mott insulators with negligible spin-orbit coupling of relativistic origin and with active orbital degrees of freedom [58]-e.g. to the above-mentioned case of transition metal oxides with partially filled 3d orbitals.
On the other hand, to the best of our knowledge, such analysis has not been done for the systems with strong on-site coupling between spins and orbitals [59]-as in the abovediscussed case of the transition metal oxides with partially filled 4d and 5d orbitals and strong spin-orbit coupling. We stress that in this case the spin and orbital degrees of freedom can get entangled as a result of both the nearest neighbor exchange interactions as well as on-site spin-orbit coupling. In fact, one typically implicitly assumes that the spin-orbital entanglement should be nonzero, since the spin S and orbital L operators couple at each site into a total angular momentum J = S + L [60]. The latter, "spin-orbital entangled" operators (also called pseudospins), then interact as a result of the exchange processes in the "relativistic" Mott insulators and are best described in terms of various effective pseudospin models, such as the Kitaev-like model discussed above [61]. Finally, very few studies discuss the problem of the evolution of a spin-orbital system between the limit of weak and strong spin-orbit coupling [62][63][64][65].
Here we intend to bridge the gap between the understanding of the spin-orbital physics in the above two limits. We ask the following questions: (i) What kind of evolution does the spin-orbital entanglement develop with increasing spin-orbit coupling? (ii) Can one always assume that in the limit of strong spin-orbit coupling the spin-orbital entanglement is indeed nonzero? (iii) How does the spinorbital entanglement arise in the limit of the strong spin-orbit coupling?
To answer the above questions, we formulate a minimal 1D model with S = 1 /2 spin and T = 1 /2 orbital (pseudospin) degrees of freedom. The model has the SU(2)⊗SU(2) intersite interactions between spins and orbitals which are supplemented by the on-site spin-orbit coupling of the Ising character-its detailed formulation as well as its relevance is discussed in Sec. II. Using exact diagonalization (ED), the method of choice described in Sec. III, we solve the 1D model and evaluate the various correlation functions used to study the entanglement. Next, we present the evolution of the ground-state properties as a function of the model parameters: in Sec. IV A for different values of the three model parameters, and in Sec. IV B for a specific choice of the relation between the two out of the three model parameters. We then show two distinct paths of ground-state evolution in Secs. IV B 1 and IV B 2. The evolution of the exact spectra of the periodic L = 4 chain is analyzed in Sec. IV B 3 for increasing λ. We discuss obtained numerical results utilizing mapping of the model onto an effective XXZ model in Sec. V, which is valid in the limit of the strong on-site spin-orbit coupling. We use the effective model to understand (i) how the spin-orbital entanglement sets in the model system and (ii) how it depends on the value of the on-site spin-orbit coupling constant λ. The paper ends with the conclusions presented in Sec. VI and is supplemented with an Appendix which discusses in detail the mapping onto the effective XXZ model in the limit of large spin-orbit coupling λ → ∞.

D. Practical note on the organization of the paper
We note that, whereas the main results of the paper are given in the extensive Sec. IV, some of the main results can be understood by using a mapping onto an effective XXZ model in Sec. V. Thus, we refer the interested audience looking for the more physical and intuitive understanding of (some of) the obtained numerical results to the latter section. Finally, we stress that all the important results of the paper are not only listed but also discussed in detail in Sec. VI A.

II. MODEL
In this paper, we study a spin-orbital model H defined in the Hilbert space spanned by the eigenstates of the spin S = 1 /2 and orbital (pseudospin) T = 1 /2 operators at each lattice site of a 1D chain with periodic boundary conditions. The model Hamiltonian consists of two qualitatively distinct terms, The first term H SE describes the intersite (super)exchange interactions between spins and orbitals. The spin-orbital ("Kugel-Khomskii") exchange reads (2) where J > 0 is the exchange parameter and the constants α and β are responsible for the relative strengths of the individual spin and orbital exchange interactions. This 1D SU(2)⊗SU(2) symmetric spin-orbital Hamiltonian has been heavily studied in the literature-it is exactly solvable by the Bethe ansatz at the SU(4) point, i.e., when α = β = 1 /4 [66][67][68], has a doubly degenerate ground state at the so-called Kolezhuk-Mikeska point α = β = 3 /4 [69,70], and was studied using various analytical and numerical methods for several other relevant values of {α, β} parameters [39,[71][72][73][74][75]. In particular, the entanglement between spin and orbital degrees of freedom in such a class of Hamiltonians is extremely well understood [28,29,52]. Last but not least, it was suggested that this model may describe the low-energy physics found in NaV 2 O 5 and Na 2 Ti 2 Sb 2 O [72], CsCuCl 3 and BaCoO 3 [76], as well as in the artificial Mott insulators created in optical lattices [77,78], and the so-called Coulomb impurity lattices [79].
Altogether, this means that the spin-orbital exchange interaction has the simplest possible form [58] that can, nevertheless, describe a realistic situation found in the transition metal oxides. This, as already mentioned in the introduction, constitutes the main reason behind the choice of this form of spin-orbital intersite interaction. We note that the spin-orbital exchange would often has a more complex form. For instance, this would be the case if, e.g., three instead of two active orbitals were taken into account and the corrections from finite Hund's exchange were included (as relevant for the 5d iridates, whose spin-orbital exchange interactions are given by, e.g., Eq. (3.11) of Ref. [19]).
The second term in the studied Hamiltonian (1) describes the on-site interaction between the spin and orbital degrees of freedom and reads Here the parameter λ measures the strength of the on-site spinorbit coupling term (of relativistic origin). The above Ising form of the spin-orbital coupling was chosen as the simplest possible and yet nontrivial one. Moreover, exactly such a form of the spin-orbit coupling is typically realized in systems with two active orbitals. This is the case of, e.g., the active t 2g doublets in YVO 3 [80] and Sr 2 VO 4 [81], the molecular π orbitals of KO 2 [82], and optical lattices. In fact, such a highly anisotropic form of spin-orbit coupling is valid for any system with an active orbital doublet, either two p (p x and p y ) or two t 2g (xz and yz) orbitals.

III. METHODS AND CORRELATION FUNCTIONS
As we are interested in quantum entanglement, and moreover, the exchange Hamiltonian (2) itself bears a rather complex quantum many-body term ∝ (S i S j )(T i T j ), we opt for the exact diagonalization (ED) method, which preserves the quantum fluctuations in the numerically found ground state. More specifically, we choose the ED calculations, since (i) it allows us to investigate the system ground state in a numerically exact manner and in a completely unbiased way, which for the first study of its kind is usually selected as the method of choice; (ii) the analytically exact Bethe ansatz approach can only be applied to a few selected values of the model parameters; (iii) the ED calculations can be relatively easily repeated for a number of model parameters and can typically address the qualitative properties of the ground state rather well; and (iv) we are interested here in rather local correlations which follow from the local spin-orbit and nearest neighbor exchange interactions.
We calculate the properties of the ground state of model (1) on finite chains with periodic boundary conditions. We utilize chains of length L = 4n sites, where n is an integer number, in order to avoid a degenerate ground state appearing in the case of a (4n + 2)-site chain (see Table 1 of Ref. [67]). For chains L = 4, a standard full ED procedure is performed, while for L = 8, 12, 16, and 20 sites we restrict the ED calculations to the Lanczos method [83].
To capture the changes in the ground state of the spinorbital model at increasing spin-orbit coupling λ, we define and investigate the following correlation functions which will be used, besides the von Neumann spin-orbital entanglement entropy [84], to monitor the evolution of the ground state with changing model parameters: (i) The intersite spin-orbital correlation function C SO : 013353-3 which measures the intersite (nearest neighbor) correlation between the spin and orbital degrees of freedom and has already been used in the literature as a good qualitative estimate for the von Neumann spin-orbital entanglement entropy [55,56]. This correlator can also be used to monitor the failure of the mean field decoupling between the spins and orbital pseudospins once C SO = 0 [27].
(ii) The intersite spin S or orbital T correlation function, which measures the intersite (nearest neighbor) correlation between the spin (orbital) degrees of freedom and is therefore sensitive to the changes in the ground-state properties taking place solely in the spin (orbital) subspace. We emphasize that these two functions are defined on equal footing in the model with SU(2)⊗SU(2) spin-orbital superexchange.
(iii) The γ component S γ γ of spin scalar product where γ = x, y, z. This function measures the component γ of the scalar product and thus allows one to investigate possible anisotropy of the intersite (nearest-neighbor) correlations between the spin degrees of freedom. The orbital scalar product component T γ γ is defined analogously to Eq. (7). (iv) Crucial for the systems with finite spin-orbit coupling is the on-site spin-orbit correlation function O SO : which measures the correlations between the z components of the spin and orbital operators on the same site. The precise form of this correlator is dictated by the Ising form (3) of the spin-orbit coupling present in Hamiltonian (1). Conveniently, the function (8) is one of the generators of the SU(4) group [68], which proved to be quite useful for examining the range of the SU(4)-symmetric ground state.

A. Von Neumann entropy in a general case
The main goal of this paper is to determine how the spinorbital entanglement changes in the spin-orbital model (1) upon increasing the value of the spin-orbit coupling λ. To this end, we first define the entanglement entropy calculated for a system that is bipartitioned into two subsystems: A and B. Typically such a subdivision refers to two distinct parts of the real [10][11][12][13] or momentum [16] space. Here, however, it concerns spin (A) and orbital (B) degrees of freedom [27][28][29]53].
A standard measure of the entanglement entropy between subsystems A and B in the ground state |GS of a system of size L is due to von Neumann [84]. It is defined as S vN = −Tr A {ρ A ln ρ A }/L and is obtained by integrating the density matrix, ρ A = Tr B |GS GS| over subsystem B. Consequently, in this paper we use the following definition of the von Neumann spin-orbital entanglement entropy: where is the reduced spin-only (S) density matrix with the orbital (T ) degrees of freedom integrated out. The spin-orbital von Neumann entropy is calculated using ED on L-site chain for model (1) and is shown as function of the parameters {α, β} for three representative values of the spin-orbit coupling λ in Fig. 1. In perfect agreement with Refs. [29,52], the von Neumann entropy S vN is finite in a rather limited region of the {α, β} parameters for λ = 0, i.e., in the entangled spin-orbital phase near the origin α = β = 0. The nonzero entanglement in that case is well understood and attributed to the onset of the dominant antiferromagnetic (AF) and alternating orbital (AO) fluctuations in the ground state without broken symmetry [29,52]; see discussion below in Sec. IV B 4.  Interestingly, a finite but small spin-orbit coupling λ < λ CRIT (λ CRIT is discussed in more detail in Sec. IV B) does not substantially increase the region in the {α, β}parameter space for which the spin-orbital entropy is nonzero; cf. Figs. 1(a) and 1(b). A drastic change in the behavior of the spin-orbital von Neumann entropy only happens for the dominant spin-orbit coupling λ > λ CRIT . In this case, the region of nonzero spin-orbital entanglement is not only much larger but also takes place for different values of the {α, β}-parameter space. For instance, it is remarkable that the von Neumann entropy in the case of λ > λ CRIT almost does not depend on α along the lines of constant α + β. Moreover, finite entanglement is activated when α + β > −1/2-however, the value of the von Neumann entropy strongly decreases for α and β located above the stripe given by the inequalities −1/2 α + β 2 and showing the highest value of entropy. In fact, it will be shown later (see Sec. V) that the von Neumann entropy is expected to vanish in the limit of α + β → ∞.
Altogether, we observe that (i) in the limit of small λ < λ CRIT the spin-orbital entanglement entropy does not change substantially with regard to the case with vanishing spin-orbit coupling, (ii) in the limit of large λ > λ CRIT the spin-orbital entanglement can become finite even if it vanishes for λ = 0, though it can also happen that (iii) in the limit of large λ > λ CRIT the spin-orbital entanglement vanishes when α + β < −1/2.

B. Von Neumann entropy for β = −α: Two distinct evolutions for increasing λ
In order to better understand the physics behind the observations (i) and (ii) discussed at the end of the previous subsection, here we study in great detail the onset of the spinorbital entanglement once β = −α. As shown in Fig. 1, for these values of the model parameters the region of the nonzero spin-orbital entanglement increases dramatically with the increasing value of the spin-orbit coupling λ.
We present in Fig. 2 the von Neumann spin-orbital entanglement entropy S vN (9) and the three spin-orbital correlation functions in the ground state of Hamiltonian (1) with α = −β, calculated using ED on an L = 12-site chain. We begin the analysis by comparing the values of the three spin-orbital correlation functions (4), (8), and (5) against the von Neumann spin-orbital entanglement entropy; see Figs. 2(b), 2(c), and 2(f). We observe that only the intersite spin-orbital correlation function C SO can be used as a qualitative measure for the von Neumann entropy, consistent with previous studies [55,56].
In particular, the on-site spin-orbit correlation function O SO cannot be used to monitor the entanglement entropy, for it measures the correlations between spins and orbitals locally and on the Ising level only. Nevertheless, both O SO as well as S vN can be used to identify various quantum phases obtained in the {α, λ} parameter space of the Hamiltonian, as suggested before for system with negligible spin-orbit coupling [56].
Next, we study the evolution of the von Neumann spinorbital entanglement entropy with increasing spin-orbit coupling λ for various values of the parameter α; see Fig. 2(a). We observe that in the representative α ∈ [−1, 1] interval there exist three distinct regimes of the value of the von Neumann entropy: (i) two compact areas in the {α, λ} parameter space for which the von Neumann entropy is vanishingly small, which exist in the large parameter range of ; and (iii) the compact area in the {α, λ} parameter space for which the von Neumann entropy is neither negligible nor takes maximal value, which exists in the relatively small parameter range between cases (i) and (ii). In order to understand the onset of these three distinct regimes, we study below two qualitatively different cases of the von Neumann entropy evolution with the increasing spin-orbit coupling: case A with |α| = 0.5 and case B with α = 0 [shown with dashed lines in Fig. 2(a)].

From a product state to highly entangled state
When |α| = 0.5 (case A), the evolution of the von Neumann spin-orbital entanglement entropy S vN (9) with increasing spin-orbit coupling can be well approximated by a logistic function; see Fig. 2(d). The von Neumann entropy has infinitesimally small values for λ/J 10 −1 , experiences a rapid growth for λ/J ∈ (10 −1 , 10 1 ), and saturates at approximately 0.5 for λ/J 10 1 . Comparable behavior is observed for the intersite spin-orbital correlation (C SO ), which, as already discussed, is a good and computationally not expensive qualitative measure for the von Neumann entropy; see Figs. 3(a), 3(c), and 3(e). Crucially, the latter calculations are obtained for the spin-orbital chains of different length and [as well visible in Figs. 3(a), 3(c), and 3(e)] show relatively small finite-size effects. This means that indeed the von Neumann entropy S vN depends here mainly on short-range processes and can remain negligibly small for a finite value of the spin-orbit coupling even in the thermodynamic limit. Finally, Figs. 2 and 3 allow us to define the critical value λ CRIT for case A as being located in an interval of rapid growth of the spin-orbital entanglement: λ CRIT /J ∈ (10 −1 , 10 1 ).
While the nature of the quantum phase for large spin-orbit coupling λ > λ CRIT is discussed in detail in Sec. V B, here we merely mention that in this case the value of the spin-orbital entanglement entropy saturates at about 0.5 (0.504 for L = 12 site chain) per site. Hence, we call this quantum phase a highly entangled state. Besides, in this case also the absolute value of the on-site spin-orbit correlation function O SO takes its maximal value, while the spin (and orbital) correlation function S (T ) is weakly AF (AO).
Next, we focus on the properties of the ground state obtained for small λ < λ CRIT . To this end, we investigate the evolution of the two other correlation functions, the onsite spin-orbit correlation function O SO and the spin correlation function S, for λ < λ CRIT and |α| = 0.5; see Figs. 2(c) and 2(f). We observe that whereas the on-site spin-orbit correlation function shows vanishingly small values in this limit, the spin correlation function S 0.25 (S −0.45) for α = 0.5 (α = −0.5), thus behaving similarly to the 1D FM (AF) chain, respectively. We note that the (unshown) analogous nearest neighbor orbital correlation function T calculated for α = ±0.5 takes complementary values to the spin correlation function for α = ∓0.5, i.e., T = −S. Such behavior is again observed for chains of various lengths, with S(T ) better approximating the expected AF Bethe ansatz value for larger chains; see Fig. 4. Altogether, this shows that the quantum phase that is observed for λ < λ CRIT qualitatively resembles the phases obtained in the limit of λ = 0: The FM ⊗ AO (AF ⊗ FO) for α = 0.5 (α = −0.5), respectively.
The above discussion contains just one caveat. Let us look at the evolution of the anisotropic spin (and orbital) correlation function S γ γ (and T γ γ ) with the increasing spinorbit coupling λ; see Figs. 4(a), 4(c) and 4(e). We notice that whenever λ/J > 0 for α = 0.5 there exist an anisotropy between the zz (solid red lines) and the planar (xx, yy, solid green lines) correlation functions-which is absent for λ = 0.
However, for λ/J 3 × 10 −1 the anisotropy is only partial, being absent in the strongly AF T γ γ correlations, in contrast to the S γ γ correlations. In fact, S δδ (where δ = x, y) stays positive as in λ = 0 case while S zz becomes negative. In this way, the energy coming from the finite spin-orbit coupling is minimized in the ground state without qualitatively changing the nature of the FM ⊗ AO and AF ⊗ FO ground states, allowing, however, for a very small value of the spin-orbital entanglement. This is the reason why, in what follows, this quantum phase is called a perturbed FM ⊗ AO product state.

From SU(4) singlet to a highly entangled state
We now investigate how the von Neumann spin-orbital entanglement entropy S vN evolves with the spin-orbit coupling once α = 0 (case B), i.e., from its finite value for the SU(4)singlet ground state at λ = 0 [29,52] to an even higher value obtained in the limit of large λ/J in the highly entangled state (i.e., the state already encountered in case A). To this end, we first note that the von Neumann entropy S vN at α = 0 changes with the spin-orbit coupling in a qualitatively different manner than in the case of |α| = 0.5; see Fig. 2(e). While we again encounter a monotonically growing function in λ, which saturates at about 0.5 for λ/J 0.2, this function seems to be discontinuous at three particular values of λ and three kinks (for L = 12 sites) that can be easily identified in Fig. 2(e). A similar behavior is encountered in the qualitative measure for the von Neumann entropy-the spin-orbital correlation function C SO ; see Fig. 2

(b) and Figs. 3(b), 3(d) and 3(f).
As a side note, let us mention that once λ = 0 and α = β = 0 the model (1)  where x ≡ 1/L. As a result, the position of the first kink k 1 converges to λ/J = 0 when L → ∞, and it is indeed reasonable to expect that infinitesimal λ modifies weakly spin-orbital correlations in the thermodynamic limit. In contrast, the final kink k f would then shift to λ/J 0.201. Therefore, for the case B we define λ CRIT as a single number: λ CRIT /J 0.2. Altogether, this means that the quantum phase encountered for 0 < λ < λ CRIT does not disappear in the thermodynamic limit and that its spin-orbital entanglement grows with the increasing spin-orbit coupling in a continuous way. To contrast this intermediate phase with the one showing the maximal value of entanglement at λ > λ CRIT , we call it an intermediate entangled state.
To better understand the properties of this phase, we also consider the spin correlation function S, the anisotropic spin S γ γ , and the orbital T γ γ correlation functions; see Figs. 2(c) and 2(f); Figs. 3

(b), 3(d) and 3(f); and Figs. 4(b), 4(d) and 4(f).
Similarly to the von Neumann entropy, also O SO or C SO correlation functions show kinks due to finite-size effects which are expected to disappear in the thermodynamic limit. Noticeably, the behavior of S, S γ γ , and T γ γ is quite distinct with regard to the one observed both for the highly entangled phase and seemingly for the perturbed FM ⊗ AO or AF ⊗ FO phases. This shows that the intermediate entangled phase observed at α = 0 and for 0 < λ < λ CRIT is indeed qualitatively different and constitutes a "genuine" quantum phase.

Exact spectra for L = 4 at increasing λ
We also note that the phase transition to the highly entangled phase with increasing λ is detected by level crossing in Fig. 6(b) and by the discontinuity in the derivative ( ∂E ∂λ ), which appears as the only kink for L = 4 -site chain; cf. Fig. 5(a). Other phase transitions occur by varying α-here at λ/J = 0.1 a phase transition is found from the FM (FO) phase with S tot = 2 and T tot = 0 (S tot = 0 and T tot = 2) to an entangled SU(4) phase (with all 15 generators being equal to 0) at |α| 0.08; see Fig. 6(c). At this latter phase transition, one finds also a discontinuous change of the von Neumann entropy [56]. For convenience, we introduce here the energỹ E which does not decrease with increasing λ as in Fig. 6(b). For a chain of length L = 4, it is defined as follows:  (12), with their degeneracies d, obtained for the periodic L = 4 chain at α = β = 0 describing the spin-orbital model Eq. (1) at λ = 0 and for two representative values of λ (λ/J = 0.1 and 10), standing for weak and strong spin-orbit coupling. To get more insight into the evolution of the spectra with increasing spin-orbit coupling ∝ λ, we consider in more detail the exact spectra of the L = 4 periodic chain; see Table I. At λ = 0, the ground state is the SU(4) singlet with energy E = −0.75. The degeneracies of the excited states follow from the S and T quantum numbers. Indeed, several states with higher values of S and T exhibit huge degeneracies. Weak spin-orbit coupling λ/J = 0.1 perturbation of the superexchange introduces the splittings of degenerate excited states and in fact the spectrum is quite dense; see Fig. 6(c) and the data in Table I. However, the SU(4) singlet ground state is still robust as shown by the correlation functions C SO and O SO which do not change from their λ = 0 values; see the dotted line in Fig. 5(a). It indicates that the spin-orbit term does not align here spin S z and orbital T z components. The energy of the ground state (12) is just moved by 2λ from −0.75 toẼ = −0.55 (Table I), and the spin-orbit coupling does not modify the ground state.
At λ CRIT = 0.14219J, the energies of the ground state and of the lowest energy excited state cross and a completely different situation arises-then the von Neumann entropy changes in a discontinuous way to the value corresponding to the strongly entangled state (for L = 4), and the spin-orbit correlation O SO drops to −0.25; see Fig. 5(a). Since λ CRIT is defined based on the changes in the ground state, for λ > λ CRIT the full energy spectrum does change further and consists of several bands of states, separated by gaps of the order of λ; see Fig. 6(d). The energies and their degeneracies within the lowest band of states are shown in the last two columns of Table I at large λ/J = 10 for the chain of length L = 4.
We have verified that the 16 low-energy states displayed in Table I for λ/J = 10 collapse to the spectrum of the XY model in the limit of λ → ∞, with the degeneracies 1, 2, 10, 2, and 1, as expected and discussed in more detail in Sec. V A. Thus, in general, the spectrum consists of energy bands with the energies increasing in steps of λ, depending on the number of sites at which the spin-orbit coupling aligns the expectation values of spin and orbital operators, S z i T z i , at each site i. In this regime, the spectra are dominated by the spin-orbit coupling. Note that the highly entangled phase can indeed be regarded as a qualitatively unique phase, irrespectively of the value of α-provided that β = −α and that λ > λ CRIT .

Summary
We have discussed in detail the evolution of the spin-orbital entanglement and its impact on the quantum phases, with the increasing value of the spin-orbit coupling λ for two representative values of the parameter α. We can now extend the above reasoning to the other values of α, keeping β = −α. However, in order to obtain a quantum phase diagram of the model, we still need to investigate whether the transitions between the obtained ground states could be regarded as phase transitions or are rather just of the crossover type. Dependence of the ground-state energy on the model parameters (see Fig. 6) as well as the analytic characteristics of the von Neumann entropy [see Figs. 2(d)-2(e); cf. Refs [87,88]] suggests that the transitions along cuts A and B [ Fig. 2(a)] are of distinct character. Whereas in case A the energy (as well as the von Neumann entropy) shows an analytic behavior across the transition [ Fig. 6(a)], in case B such behavior (both in energy as well as in von Neumann entropy) is clearly nonanalytic [ Fig. 6(b)]. This points to a crossover (phase) transition in case A (B), respectively.
Altogether, this allows us to draw, on a qualitative level, a quantum phase diagram in the {α, λ} parameter space (with β = −α); see Fig. 7 (colorful vertical plane). As already discussed in Sec. IV, there are four distinct ground states (first two shown in Fig. 7 in blue, and the other two in green and red, respectively): (i) the perturbed FM ⊗ AO state for α 0.08 and λ < λ CRIT , (ii) the perturbed AF ⊗ FO state for α −0.08 and λ < λ CRIT , (iii) the intermediate entangled state for |α| 0.08 and 0 < λ < λ CRIT , and (iv) the highly entangled state for λ > λ CRIT and for all values of α. The latter state is discussed in more detailed in Sec. V. The four clearly distinct states are supplemented by two crossover regimes (shown in yellow in Fig. 7), which separate phases (i) and (ii) from phase (iv)-see also discussion above.
It is instructive to place the above phase diagram in the context of the one already known from the literature and obtained for Hamiltonian (1) in the limit of the vanishing spin-orbit coupling λ but varying values of both α and β [29,39,52,[66][67][68][69][70][71][72][73][74][75]. As can be seen on the horizontal plane of Fig. 7, the λ = 0 phase diagram consists of three simple product phases (AF ⊗ FO, FM ⊗ AO, and FM ⊗ FO) as well as two spin-orbital entangled phases (cf. Fig. 1): a phase with previously mentioned global SU(4)-symmetric singlet ground state and gapless excitations [68] and a phase with the ground state breaking the Z 2 symmetry and opening a finite gap by forming the two nonequivalent patterns of the spin and orbital dimers [29,52]. We would like to emphasize at this point that the finite-size effects for the spin-orbital model (at λ = 0) calculated on chains of length L = 16 (the maximal size studied in Ref. [52]) and L = 20 (the maximal size studied here) are already relatively small [52]. This may suggest that the schematic phase diagram of Fig. 7 is qualitatively correct also in the thermodynamic limit.  (1) in the discussed regime of the parameters. The limit of α = −β is depicted by the colorful vertical plane and is based on the results from Sec. IV B, whereas the four distinct phases are depicted with their names and separated by solid lines (IES stands for the intermediate entangled state) and the two crossover regimes are denoted by yellow color and separated by the dashed lines. The limit λ = 0 is depicted by the horizontal plane and is adopted from Fig. 1 of Ref. [52]-see text for further details. We note that the shape of the phase boundaries depends on the logarithmic scale of λ, chosen here for convenience. The schematic phase diagram is based on the ED results on small clusters; see text for the validity of these results in the thermodynamic limit.

A. Effective XXZ model
To better understand the numerical results obtained in Sec. IV for the Hamiltonian (1) in the limit of the large spin-orbit coupling, λ > λ CRIT , we derive an effective lowenergy description of the system. In fact, as already discussed in the introduction, such an approach has become extremely popular in describing the physics of the iridium oxides [61], for it has led to the description of the latter in terms of effective Heisenberg or Kitaev-like models. To obtain such an effective description for the case of large spin-orbit coupling, λ > λ CRIT , we first obtain the eigenstates of the spin-orbit coupling Hamiltonian (3): These are two doublets, separated by the gap E = λ. Next, we restrict the Hilbert space to the lowest doublet {|↑− , |↓+ }, where |↓ (|− ) denotes the state with S z = − 1 /2 (T z = − 1 /2) quantum number. Lastly, we project the intersite Hamiltonian (2) onto the lowest doublet (see the Appendix for details) and obtain the following effective model: whereJ z i = − 1 2 (n i,|↑− + n i,|↓+ ) is an effectiveJ z = 1 /2 pseudospin operator.
Interestingly, it turns out that this effective Hamiltonian describes exactly a spin 1 /2 XXZ chain. Moreover, in the limit of α = −β the Ising interaction in Eq. (13) disappears and we obtain an AF XY model. Thus, resembling the iridate case [61], the effective model in the limit of large spin-orbit coupling has a surprisingly simple form.

B. Validity of the effective XXZ model:
Benchmarking α = −β case First, let us show that the effective XXZ model indeed gives the correct description of the ground state of the full spin-orbital model (1) in the limit of λ > λ CRIT . To this end, we compare the spin-orbital correlation functions calculated using the effective and the full models.
We first express the spin-orbital correlation function C SO , the on-site spin-orbit correlation function O SO , and the anisotropic spin (orbital) correlation functions S γ γ (T γ γ ) in the basis spanned by the two lowest doublets per site {|↑− , |↓+ }-see the Appendix for the explicit formula. Next, we compare the values of the correlation functions in the two special α = −β cases, already discussed above: (i) case A with |α| = 0.5 and (ii) case B with α = 0. As can be seen in Figs. 3 and 4, the correlation functions calculated using the two distinct models agree extremely well once λ/J 1-100 (λ/J 0.2) in case A (B), respectively. We note that calculations performed for other values of the {α, β} parameters (not shown) also show that the effective model describes the ground-state properties in the limit of λ > λ CRIT well. Moreover, once λ/J 10 6 , the ground and lowest lying excited states are quantitatively the same in the full and the effective models.

C. Why the spin-orbital entanglement can vanish
Having derived the effective model-and having shown its validity-we now discuss how it can help us with understanding one of the crucial results of the paper: How can the spinorbital entanglement vanish in the limit of large spin-orbit coupling λ > λ CRIT ?
We start by expressing the measure for the spin-orbital entanglement for nearest neighbors, the spin-orbital correlation C SO , in the basis of the effective model (see the Appendix for details): where the averages are calculated in the ground state. To evaluate Eq. (14), we calculate expectation values of the effective pseudospin operators using ED, which we show in Fig. 8. (We note in passing that the presented ED results for an XXZ L = 10 site chain agree well with those which were published earlier, cf. Ref. [89].) The obtained ground state of the effective Hamiltonian (13)  α + β < −1/2. In conclusion, the spin-orbital entanglement for α + β < −1/2 vanishes because not only the on-site interaction between spins and orbitals but also the intersite interactions in the ground state are of purely Ising type, and effectively the ground state is just a product state with no spin-orbital entanglement.

D. Why the spin-orbital entanglement can be finite
The effective model (13) can also be used to explain the presence of finite spin-orbital entanglement in the limit of large spin-orbit coupling λ > λ CRIT while it vanishes in the λ = 0 limit. Let us first look at the already discussed in detail β = −α case: In this case and in the λ = 0 limit, the term in (2) which is explicitly responsible for the spin-orbital entanglement, ∝ (S i S j )(T i T j ), can become relatively small for large α or β due to the presence of the αT i T j and βS i S j terms. Consequently, the region of significant spin-orbital entanglement is quite small without spin-orbit coupling along the β = −α line; see Fig. 1(a). This situation, however, drastically changes in the limit of large λ > λ CRIT , as discussed below.
Specifically, downfolding the exchange Hamiltonian (2) term by term onto the effective Hamiltonian (13) should reveal the origin of the spin-orbital entanglement in the large spin-orbit coupling limit. First, the αT i T j and βS i S j terms of Eq. (2) upon projecting onto spin-orbit coupled basis produce αJ z iJ z j and βJ z iJ z j , resulting in the Ising terms in the effective model (13). Note that in the case that β = −α, these Ising terms disappear. Second, the term responsible for the spinorbital entanglement, i.e., (S i S j )(T i T j ) (cf. above), reduces exactly to the XY terms in the effective model. These terms do not vanish once β = −α. In fact, in this special limit the whole effective Hamiltonian is obtained from the term that is fully responsible for the spin-orbital entanglement in the original Hamiltonian. Finally, as the ground state of the XY Hamiltonian carries spatial entanglement in pseudospins J, we expect the spin-orbital entanglement to be finite in the limit of large spin-orbit coupling λ > λ CRIT and once β = −α.
The above reasoning is confirmed by calculating the two contributions to the intersite spin-orbital correlation func-tionC SO in the effective model once β = −α. This can be done analytically for the XY model: J z iJ z i+1 = −1/π 2 and J x iJ x i+1 = J y iJ y i+1 = −1/(2π ). (These results agree with the correlations calculated using ED and presented in Fig. 8.) The above discussion can now be extended to the case that β = −α and α + β > −1/2, for which finite, though increasingly small for large and positive α + β, spin-orbital entanglement can be observed; see Fig. 1(c). Such result can be understood by using the effective model and by noting that the intersite spin-orbital correlationC SO is always finite provided that α + β is finite and α + β > −1/2. This is be- It is then only in the limit α + β → ∞ that the spin-orbital entanglement can vanish, for the ground state of the XXZ model is pure Ising antiferromagnet. (A completely different situation occurs once α + β < −1/2, i.e., for the FM ground state of the effective XXZ model, as already discussed in the previous subsection-that explains why the spin-orbital entanglement can sometimes vanish even in the limit of large spin-orbit coupling, λ > λ CRIT .)

A. Entanglement induced by spin-orbit coupling
In conclusion, in this paper we studied the spin-orbital entanglement in a Mott insulator with spin and orbital degrees of freedom. We investigated how the spin-orbital entanglement gradually changes with the increasing value of the on-site spin-orbit coupling. The results, obtained by exactly diagonalizing a 1D model with the intersite SU(2)⊗SU(2) spin-orbital superexchange ∝ J and the on-site Ising-type spin-orbit coupling ∝ λ, reveal the following: (1) For small λ < λ CRIT [90], (a) in general, the spin-orbital entanglement in the ground state is not much more robust than in the λ = 0 case; (b) if the ground state had finite spin-orbital entanglement for λ = 0, it is driven into a novel spin-orbital strongly entangled phase upon increasing λ; and (c) if the ground state did not show spin-orbital entanglement for λ = 0, it still shows none or negligible spin-orbital entanglement upon increasing λ.
(2) In the limit of large λ > λ CRIT , (a) in general, the spin-orbital entanglement in the ground state is far more robust than in the λ = 0 case; (b) the ground state may be driven into a novel spinorbitally entangled phase even if it does not show spinorbital entanglement for λ = 0; and (c) the ground state may still show vanishing spinorbital entanglement, but only if the quantum fluctuations vanish in the ground state of an effective model (as is the case of an Ising ferromagnet). The statements mentioned under point 2. above, concerning λ > λ CRIT , constitute, from the purely theoretical perspective, the main results of this paper. In particular, they mean that (i) the spin-orbital entanglement between spins and orbitals on different sites can be triggered by a joint action of the on-site spin-orbit coupling (of relativistic origin) and the spin-orbital exchange (of the Kugel-Khomskii type); and yet, (ii) the onset of the spin-orbital entanglement in such a model does not have to be taken for granted, for it can vanish even in the large spin-orbit coupling limit.
Crucially, we have verified that the spin-orbital entanglement can be induced by the spin-orbit coupling, for the latter interaction may enhance the role played by the spinorbitally entangled (S i S j )(T i T j ) term by quenching the bare spin (S i S j ) and orbital (T i T j ) exchange terms in an effective low-energy Hamiltonian valid in this limit. Interestingly, such mechanism can be valid even if the spin-orbit coupling has a purely classical Ising form (as, for example, in the case discussed in this paper). For a more intuitive explanation of these results, in Sec. V we presented a detailed analysis of the effective low-energy pseudospin XXZ model.

B. Consequences for correlated materials
The results presented here may play an important role in the understanding of the correlated systems with nonnegligible spin-orbit coupling-such as, e.g., the 5d iridates, 4d ruthenates, 3d vanadates, the 2p alkali hyperoxides, and other to-be-synthesized materials. To this end, we argue that, even though obtained for a specific 1D model, some of the results presented here are to a large extent valid also for these 2D or 3D systems: First, this is partially the case for the results obtained in the limit of large λ > λ CRIT . In particular, the mapping to the effective XXZ model is also valid in 2D and 3D cases. Moreover, one can easily verify that the spin-orbital correlation function [C SO , Eq. (14)], which measures spin-orbital entanglement, never vanishes also in the 2D and 3D cases, unless the quantum fluctuations completely disappear (as is the case of the 2D or 3D Ising ferromagnet or antiferromag-net). Therefore, the main conclusions from Secs. IV C and IV D are also valid in 2D and 3D cases and consequently also point 2 of Sec. VI A holds. This means that, for example, the results obtained here would apply to any Mott insulator with two active t 2g orbitals with small Hund's coupling and with λ > λ CRIT (such as, e.g., Sr 2 VO 4 [81]).
Naturally, the question remains as to what extent one could use the reasoning discussed here to the understanding of the spin-orbital ground state of the probably most famous Mott insulators with active orbital degrees of freedom and large spin-orbit coupling-the 5d iridates (such as, e.g., Sr 2 IrO 4 [41], Na 2 IrO 3 , Li 2 IrO 3 , etc. [50]). Here we suggest that, while the situation in the iridates might be quite different in detail and requires solving a distinct spin-orbital model with three active t 2g orbitals and an SU(2)-symmetric spin-orbit coupling (which is beyond the scope of this work), we expect point 2(b) of the concluding Sec. VI A to hold also in this case: In fact, the quantum nature of the Heisenberg spin-orbit coupling of the iridates (in contrast to the classical Ising spin-orbit coupling studied in this paper), should only facilitate the onset of the spin-orbital entanglement. Thus, we suggest that in principle also for the iridates the ground state may be driven into a novel spin-orbitally-entangled phase even if it does not show spin-orbital entanglement for λ = 0.
Second, we suggest that also the fact that the spin-orbit coupling does not induce additional spin-orbital entanglement in the limit of small λ < λ CRIT will carry on to higher dimensions and to spin-orbital models of lower symmetry-for a priori there is no reason why the tendency observed in a 1D (and highly symmetric) model, toward a more classical behavior should fail in dimensions higher than one (and for more anisotropic models). Thus, in general the spin-orbital entanglement of the systems with weak spin-orbit coupling λ < λ CRIT and Ising-like spin-orbit coupling [80], such as, e.g., the alkali hyperoxides with two active molecular 2p orbitals (e.g., KO 2 [82]), should not qualitatively depend on the value of spin-orbit coupling. This means that to simplify the studies, one may, in the first order of approximation, neglect the spin-orbit coupling in the effective models for these materials.

APPENDIX: EFFECTIVE XXZ MODEL
Let us consider the Hamiltonian (1) of the main text: where the intersite interaction H SE and on-site spin-orbit coupling are described by The characteristic scales for H SE and H SOC are intersite exchange parameter J and on-site SOC λ, respectively. In the strong spin-orbit coupling limit, λ > λ CRIT , H SE can be considered as a perturbation to H SOC . The eigenstates of the full Hamiltonian (A1) in zeroth order are then obtained by the diagonalization of the on-site spin-orbit part H SOC . In our simple case, H SOC is already diagonal with two doubly degenerate energies ±λ/2.
We then project the Hamiltonian (A2) onto spin-orbit coupled basis {J,J }: H SOC SE = U † H SE U . As we are interested in the low-energy physics, we truncate Hilbert space to the lowest doubletJ and obtain effective Hamiltonian (13) from the main text: To analyze the effective model (A5) and obtain important correlation functions, we first need to establish a link between operators describing correlation functions in original {|T z , S z } basis and spin-orbit coupled {J,J } basis. To this end, we project each of the spin-orbital operators, O r = {S γ r , T γ r }, γ = {x, y, z}, r = {i, i + 1} entering the original correlation functions (4) To express the on-site spin-orbit correlation function O SO , which does not include intersite terms, in the same basis, we multiply it by a 2 × 2 identity matrix representing the neighboring site: