Supporting data and manuscript "Magnetic polarons due to spin-length fluctuations in d^4 spin-orbit Mott systems"

PHYSICAL REVIEW B 111, 195125 (2025) Magnetic polarons due to spin-length fluctuations in d4spin-orbit Mott systems Jan Revenda ,1Krzysztof Wohlfeld ,2,3and Jiˇ rí Chaloupka 1 1Department of Condensed Matter Physics, Faculty of Science, Masaryk University, Kotláˇrská 2, 61137 Brno, Czech Republic 2Institute of Theoretical Physics, Faculty of Physics, University of Warsaw, Pasteura 5, PL-02093 Warsaw, Poland 3Department of Physics and Quantum Centre of Excellence for Diamond and Emergent Materials (QuCenDiEM), Indian Institute of Technology Madras, Chennai 600036, India (Received 30 March 2025; revised 28 April 2025; accepted 29 April 2025; published 14 May 2025) Mott insulators based on 4dand 5dtransition-metal ions, where spin-orbit interaction plays a key role, can exhibit various forms of unusual magnetism. A particular example is the antiferromagnet Ca2RuO4containing d4Ru4+ions. Here the spin-orbit interaction stabilizes the nonmagnetic J=0 singlet ionic ground state, which gets dynamically mixed, via exchange interactions, with low-energy J=1 ionic excitations. Thanks to a sufficient strength of the exchange, these excitations condense and a long-range order emerges. The resulting ordered moments are soft and prone to fluctuations of their effective length. The corresponding amplitude mode appears as a prominent magnetic excitation and complements the conventional magnons involving rotations of the moments. Motivated by this peculiar kind of magnetic order and the specific spectrum of magnetic excitations, we study their influence on the propagation of doped carriers. To this end, we construct a microscopic model including both d4and d5degrees of freedom and address the propagation of an injected electron by employing self-consistent Born approximation. We find that the electron shows a combination of both free and a polaronic type of motion, where the mobile carrier strongly interacts with an accompanying cloud of magnetic excitations. Remarkably, in the latter case it is the exotic excitation, the amplitude mode, that is found to dominate over the contribution of magnons. Our soft-spin situation thus largely contrasts with spin polarons widely discussed in the context of doped Heisenberg-type magnets based on rigid spin moments. DOI: 10.1103/PhysRevB.111.195125 I. INTRODUCTION Mott insulators with active 3dorbitals in the valence shell have traditionally been described by the so-called KugelKhomskii spin-orbital superexchange models [1,2]. These Hamiltonians may host spin and orbitally ordered ground states as well as collective excitations carrying spin (magnon) [3–5] or orbital (orbiton) [6–8] quantum numbers. While such physics is rich and far from being completely understood in itself, another important element has been stressed by Jackeli and Khaliullin in 2009 [9], namely, the large spin-orbit coupling (SOC) present in ions with 4dor 5dvalence orbitals, which is capable to fundamentally change the low-energy description of the corresponding Mott insulators [10,11]. In their example, a Mott insulator with t5 2gions on the honeycomb lattice is, in the limit of dominant SOC, described by an effective J=1 2pseudospin model hosting frustrated interactions of Kitaev type. Such an unexpected possibility of a material realization of the highly interesting Kitaev model [12] featuring a spin-liquid ground state attracted a lot of attention and opened an entirely new field of “Kitaev materials” such as Na2IrO3or α-RuCl3[13–15]. Our paper is inspired by another peculiar spin-orbit Mott insulator, Ca2RuO4, featuring four valence electrons in t2g orbitals of 4d4ruthenium Ru4+ions forming square-lattice planes. Since its discovery in 1997 [16] and observation of the metal-insulator transition (MIT) at TS=357 K [17–19] as well as onset of the AF order at TN=110 K [20], this ruthenate is subject of ongoing extensive research. For instance, very recently several phenomena related to MIT due to surface doping [21], uniaxial strain [22], or dc current [23] were discussed. Here we choose Ca2RuO4because the effectively moderate value of its SOC leads to a rich physics, perhaps as fascinating as that of the above-mentioned Kitaev materials. In fact, the strongly competing spin-orbital superexchange and SOC in this compound have been intensely debated in the literature [24–35]: On the one hand, a sufficiently large SOC can rearrange the low-energy levels of t4 2g ions to form a nonmagnetic singlet ground state, separated in energy from triplet ionic excitations by the SOC value. On the other hand, superexchange interactions favor a dynamical mixing in this singlet-triplet level structure, promoting the condensation of triplet excitations to form long-range antiferromagnetic order. Consequently, if there is right balance between the singlet-triplet splitting and the strength of the superexchange interactions, then these two effects are competing and a soft-spin magnetism emerges [26,29]. This results in an unusual excitation spectrum, which features both magnonlike excitations and pronounced spin-length fluctuations in a form of amplitude oscillations of the condensate (an amplitude or Higgs mode) [29,36]. Interestingly, exactly this regime seems to be realized in 4d4Ca2RuO4, being facilitated by a substantial lowering of the triplet excitations in the presence of a tetragonal crystal field [29]. In this paper, we focus on the interplay of the above unusual soft-spin magnetic background with doped electronlike carriers. In general, doping Mott insulators remains one of the 2469-9950/2025/111(19)/195125(20) 195125-1 ©2025 American Physical Society REVENDA, WOHLFELD, AND CHALOUPKA PHYSICAL REVIEW B 111, 195125 (2025) major challenges in the study of strongly correlated electron systems [37–39]. Perhaps the only case that is relatively well understood is the case of a single hole (or electron) doped into a single-band Mott insulator. In this case, a doped hole disrupts the magnetic structure, leading to strong coupling between the hole and magnetic excitations [37]. Since the two-dimensional (2D) Mott insulator supports long-range antiferromagnetic (AF) order, the hole couples strongly to the collective excitations of the AF order (magnons), forming a well-defined magnetic polaron quasiparticle [40–47]. The single-band magnetic polaron concept has been rigorously tested and observed in the cold-atom simulations and in the ARPES spectra of copper oxides [48–53]. In contrast, the situation in Mott insulators with more than one orbital active in the valence band, both with and without SOC, seems more intricate: On the level of effective lattice models, earlier works on 3dtransition metal compounds, such as manganites [54], cobaltates [55], or vanadates [56], suggested onset of magnetic polarons, of spin or spin-orbital nature. Similarly, also in systems with non-negligible SOC, an onset of magnetic polarons was theoretically postulated [57]. Going to the more extreme case of Kitaev systems, the dynamic interplay between the doped carriers and the exotic excitations of the spin-liquid background lead to the suggestion of electron fractionalization [58–60]. However, DFT +Ucalculations often align well with ARPES results for Mott insulators with orbital degrees of freedom. This happens for the van der Waals magnet MnPS3[61], the spin-orbit-coupled 2-1-4 iridates [62–65], and, even more surprisingly, for the Kitaev candidate materials [66,67]. Thus, it seems as mean-field approaches are often sufficient to describe the photoemission spectra of these systems. Ca2RuO4clearly illustrates the difficulties mentioned above in identifying a spin polaron in an AF Mott insulator with active orbital degrees of freedom and non-negligible SOC. In fact, there is a notable inconsistency in the approaches used to describe the magnetic versus the electronic properties of Ca2RuO4: On one side, the magnetic phenomena are best described by a strongly correlated, localized electron picture. This is particularly evident in the observation of an amplitude (Higgs) mode in its inelastic neutron scattering data [29], which was successfully interpreted by employing a complex spin-orbital superexchange model [27,29]. In contrast, explanations for the key electronic properties of Ca2RuO4have largely relied on mean-field approaches. For example, the measured ARPES spectrum has been explained using a combination of DFT and single-site dynamical meanfield theory (DMFT) [68], where local correlation effects were fully incorporated and proved crucial for understanding the onset of Hund’s-induced splitting of the dxy band, as observed in ARPES. This result was recently confirmed by the GW +EDMFT approach [69]. One reason for the relatively good agreement of mean-field-based approaches with the ARPES data is that the spectra were measured at temperatures above TN, which largely suppresses intersite correlation effects. Although some weak signatures of spin polaron physics in the measured ARPES spectra of Ca2RuO4[68] were discussed in [70], the model employed there was a largely phenomenological S=1t-Jmodel and solely the response of the dxy/dyz orbitals was taken into account. The aim of this work is to examine the spin-polaronic features in a realistic spin-orbital t-J–like model for Ca2RuO4 derived on microscopic grounds and to understand the implications of the soft-spin nature of the magnetic background. We choose to study the inverse photoemission spectrum, i.e., the propagation of a single extra electron rather than the more widely discussed single-hole or photoemission case. The reason is that injecting a single electron into a d4background of Ru4+ions results in a nominal d5configuration, which, unlike the d3configuration relevant to ARPES, is largely affected by SOC and carries a pseudospin J =1 2. The paper starts with the derived effective model specified in Sec. II. Extending largely the earlier works [29,71], it captures the interactions among low-energy d4and d5degrees of freedom with complex internal structure reflecting the spin-orbital entanglement due to SOC. In Sec. III we review the essential physics of the d4background, discussing its phases and magnetic excitation spectra. The central part of the paper consists of Secs. IV and Vpresenting a study of the single d5carrier propagating through the d4background. Section IV provides the necessary technical basis in the form of self-consistent Born approximation (SCBA) equations. These equations are then numerically solved and the results interpreted in Sec. V. For a better understanding, we examine the behavior of the single carrier as a function of the relative strength of the tetragonal crystal field as compared to the SOC. While in Ca2RuO4the ratio of these two parameters is fixed, changing its value enables us to conveniently drive the system between its two phases: nonmagnetic (NM) phase built on ionic singlets and antiferromagnetically (AF) ordered phase associated with a condensate of ionic triplets. This way the main results are uncovered: In the NM phase, the carrier seems to primarily move as a free particle, however, a closer look reveals that a pseudogap is always present in the quasiparticle dispersion. This observation is further addressed in Sec. VI with the help of a simple toy model. In contrast, the carrier motion in the AF phase is predominantly of polaronic nature, albeit remnants of free-electron motion are still visible in the regime relevant to Ca2RuO4. A detailed analysis presented in Sec. Vreveals that the observed polaronic motion is far more heavily influenced by the coupling to the amplitude mode than to the magnon one. II. MODEL In this section we describe the microscopic model used to study the propagation of d5electronlike carriers in a soft-spin d4background. The model is obtained as an effective lowenergy limit of a three-orbital t2g-based Hubbard Hamiltonian on a square lattice. We consider the Hubbard model to be in the Mott regime and including SOC of sufficient strength so that the latter is decisive for the structure of the effective lowenergy degrees of freedom, as will be discussed in Sec. II A. The resulting exchange interactions within the d4background itself and its interaction with an injected d5carrier will be exposed in Secs. II B and II C, respectively. The values of the model parameters are chosen to describe the canonical spin-orbit-coupled d4material Ca2RuO4, though we keep the tetragonal crystal field as a free control parameter allowing us 195125-2 MAGNETIC POLARONS DUE TO SPIN-LENGTH … PHYSICAL REVIEW B 111, 195125 (2025) to examine how the physics changes once the system is driven through a quantum phase transition. A. Low-energy d4and d5ionic states We start by specifying the basis of our effective model. To this end, we determine the multiplet structure of d4and d5transition-metal ions adopting the LS coupling scheme when incorporating SOC and select the relevant low-energy states. Focusing on 4dor 5dtransition metals, all the valence electrons can be assumed to reside in t2gorbitals and form either spin S=1incaseoft4 2gconfiguration or spin S=1 2in case of t5 2gconfiguration. Let us first consider the case of crystal field of cubic symmetry implying degenerate t2gorbitals. In such situation, both t4 2gand t5 2gconfigurations carry effective orbital angular momentum L=1, which is combined with the spin by virtue of spin-orbit coupling into the total angular momentum J=L+Seigenstates having J=0,1,2 (i.e., singlet-triplet-quintuplet level structure) or J=1 2,3 2(doublet-quartet), respectively [11]. We now add a tetragonal crystal-field component that is generated, e.g., by an out-of-plane contraction or elongation of the metal-O6octahedra. When including this crystal field, the t2gorbitals split as indicated in Fig. 1(c) and consequently the above states get modified with their degeneracy being partially lifted. The resulting energy levels are presented in Figs. 1(a) and 1(b) as functions of measured by the single-electron spin-orbit coupling ζ. We consider only positive as realized in Ca2RuO4 (in this case /ζ was estimated to ≈1.5inRef.[29]). Using the effective model, we aim to address the lowenergy window and the regime ζ, we thus omit (i) the d4states derived from J=2, (ii) one of the original J=1 triplets that is pushed up in energy by positive , and (iii) the high-energy d5states. The remaining states included in our basis can be compactly expressed in terms of |LzSzstates as follows: First, we consider three d4eigenstates. The singlet d4 ground state takes the form |s=−sin ϑ0|0,0+cos ϑ0 1 √2(|+1,−1+|−1,+1) (1) with the parameter ϑ0determined by δ=/ζ via tan 2ϑ0= 2√2/(1 −2δ). The relevant two members of the original J= 1d4triplet that got lowered in energy by positive are based on Jz=±1 eigenstates |T±1=±sin ϑ1|0,±1∓cos ϑ1|±1,0(2) with ϑ1satisfying tan ϑ1=1/(√1+δ2−δ). For later convenience, we will utilize their linear combinations |Tx= i √2(|T+1−|T−1) and |Ty= 1 √2(|T+1+|T−1). Second, we consider the d5doublet |f↑=cos ϑ|+1,↓−sin ϑ|0,↑, |f↓=sin ϑ|0,↓−cos ϑ|−1,↑ (3) carrying pseudospin 1 2. The parameter ϑis determined by tan 2ϑ=2√2/(1 +2δ). We also give an explicit expression for the energy splitting between |sand |Tx,ystates, denoted here as ET=E(Tx,y) FIG. 1. d4and d5ionic states. Low-lying energy levels of t4 2g (a) and t5 2g(b) ions as functions of tetragonal crystal field measured by the SOC strength ζ. The energies are plotted relative to the ground-state energy (i.e., J=0or1 2levels, respectively). The vertical scale is common to both plots. (c) Splitting among t2gorbitals induced by a tetragonal crystal field >0. Within point-charge model, positive corresponds to a compression of metal-O6octahedra in the direction perpendicular to the xy plane defined by the square lattice as shown below. The particular AF order observed in Ca2RuO4is also indicated. (d) Auxiliary angles entering the factors in the ionic wave functions of Eqs. (1)–(3) and (e). (e) Schematic picture of selected relevant ionic states: t4 2gionic ground state |s, |T+1representing one of the triplet excitations, and |f↑from the t5 2g ground-state pair. We utilize hole representation taking t6 2gstate as a reference. The components of the twoor one-hole wave functions are depicted as polar plots of the density distribution with the spin polarization indicated by color. Orbital angular momentum is represented by the circular arrows. −E(s), as the most crucial parameter of our model arising from ionic physics: ET=1 4ζ[1 +(1 −2δ)2+8−21+δ2].(4) 195125-3 REVENDA, WOHLFELD, AND CHALOUPKA PHYSICAL REVIEW B 111, 195125 (2025) Note that the ETsplitting is largely reduced from its /ζ =0 value ET=ζ/2 already at moderate ≈ζ, which enables the exchange interactions to overcome the splitting and induce the condensation of Tas discussed later. The splitting vanishes completely at /ζ →∞. Pictorial representation of some of the above states shown in Fig. 1(e) emphasizes their complex spin-orbital entangled nature (as imposed by SOC), which gets reflected in the anisotropy of the exchange interactions. For example, pseudospin 1 2carried by |f↑,|f↓doublet is formally captured by spin 1 2, but due to their internal structure, these pseudospins may be subject to highly anisotropic interactions such as in the honeycomb Kitaev systems [9]. When increasing /ζ that enters via the three variables ϑ0,ϑ1,ϑ[see Fig. 1(d)], some components of the wave functions get gradually suppressed and the complexity stemming from the spin-orbital entanglement fades away. For instance, in the case of the three d4 states, the orbital angular momentum gets fully quenched in the /ζ →∞limit and we are left with pure spin-1 situation. B. Magnetic interactions in the d4background Having selected our local (ionic) basis, we continue by specifying the bond interactions constituting the effective model. We address first the exchange interactions within d4 background, that we have derived based on the square-lattice t2gHubbard model with SOC. The resulting d4part of the effective model is identical to the model used earlier in Ref. [29] to fit the neutron data. Our derivation complements this model with explicit expressions for the exchange interactions in terms of the parameters of the underlying Hubbard model with SOC. The effective d4model can be conveniently formulated as a hard-core boson model with the hard-core bosons being assigned to the three local-basis states s,Tx, and Ty. Particularly compact exchange expressions are obtained when introducing pseudospin-1 operating in the three-dimensional Hilbert space spanned by the selected d4states. Its components expressed using the hard-core boson operators read as  Sx=−i(s†Tx−T† xs),(5)  Sy=−i(s†Ty−T† ys),(6)  Sz=−i(T† xTy−T† yTx).(7) The auxiliary spin 1 has also a physical meaning: the in-plane components  Sx, Sycarry van Vleck magnetic moment while the out-of-plane component  Szcorresponds to the magnetic moment hosted by the Tstates. Details on this correspondence and the related g-factor values can be found in the Supplemental Material of Ref. [72]. The exchange model for the Mott limit is derived in the usual way by perturbatively (to second order) eliminating the nearest-neighbor hopping in the relevant t2gHubbard model on the square lattice. In this case, there are always two of the three t2gorbitals active for a given bond direction, for example, the hopping along the xbond ij(i.e., bond parallel to the xaxis) takes the form −t s=↑,↓ (d† xy,s,idxy,s,j+d† zx,s,idzx,s,j)+H.c.(8) Similarly, for the ybonds, the orbitals xy and yz are active. Together with the local terms explained below, the resulting Hamiltonian that describes the d4background can be written as Hd4= i∈sites ET Sz i2−1 2[( S i)2−( S⊥ i)2] + ij∈bonds Jx Sx i Sx j+Jy Sy i Sy j+Jz Sz i Sz j.(9) The first local term of Eq. (9) captures the splitting of the s and Tx,yd4ionic states discussed in the previous section [see Eq. (4)]. It can also be cast to the form ETnTwith nT= T† xTx+T† yTyand simply counts the number of Tparticles at a given site, penalizing them with the splitting energy ET.The second local term of Eq. (9) is included ad hoc to account for the in-plane anisotropy observed experimentally in Ca2RuO4. It fixes the ordered moment direction in the xy plane and also opens the magnon gap needed to properly fit the neutron data. The pseudospin component  Spoints along the ordered moment direction, in the case of Ca2RuO4along (x+y)/√2 [see Fig. 1(c)], i.e.,  S i=( Sx i+ Sy i)/√2. The component  S⊥is the perpendicular one:  S⊥ i=(− Sx i+ Sy i)/√2. The exchange parameters Jxand Jyof Eq. (9) are bond selective and take the values Jx=J±δJ,Jy=J∓δJwith the signs depending on the bond direction: upper sign for an xbond, bottom sign for a ybond. This selectivity is a consequence of spin-orbital entanglement combined with the orbital selection rules for the hopping processes. The derived rather lengthy expressions for J,δJ, and Jzin terms of the Hubbard model parameters are given in Appendix A. Let us note that in Ref. [29] the same model is formulated in a rotated coordinate frame, altering its form. Our parameter δJcorresponds to A of Ref. [29]. Utilizing the definitions (5)–(7), the interactions in Eq. (9) can be converted back to the language of hardcore bosons. The main exchange processes, associated with Jxand Jy, are sketched in Fig. 2(a). They include the effective hopping of Tparticles and their creation or annihilation in properly matched pairs. The bilinear exchange terms in Eq. (9) are accompanied by biquadratic ones, but their magnitude is found to be negligible. Figure 2(b) shows the evolution of d4model parameters with varying /ζ . It was constructed by fixing the ionic parameters at the values appropriate for Ru ions in Ca2RuO4: Hubbard repulsion U=3 eV, Hund’s coupling constant JH=0.5 eV, and SOC strength ζ=0.15 eV. The hopping amplitude was set to t=0.143 eV. These values lead to a good overall agreement with the estimated J≈5.8eV, ET≈25 meV obtained by fitting the INS data on Ca2RuO4 [29], considering the relevant /ζ interval of about 1.5–2.0. However, as mentioned earlier, we will not limit ourselves solely to the Ca2RuO4case and use /ζ in a wider range as a convenient handle to drive the model through different regimes. When applying such a virtual “straining” of the crystal, we will adopt the effective parameter values given in Fig. 2 195125-4 MAGNETIC POLARONS DUE TO SPIN-LENGTH … PHYSICAL REVIEW B 111, 195125 (2025) FIG. 2. Magnetic interactions in the d4background. (a) Cartoon representation of important exchange processes in the d4background. Within second-order perturbation theory, we consider pairs of virtual electron hoppings that preserve the d4valence at neighboring sites and collect them in Hd4of Eq. (9). The latter connects various combinations of bond states made of ionic ground-state singlets sand local Texcitations. The essential contributions are the effective hopping of T(top right) and a creation (and similarly an annihilation) of pairs of T(bottom right). (b) Parameters entering the d4magnetic model (9) as functions of /ζ for fixed U=3eV, JH=0.5eV,ζ=0.15 eV, and t=0.143 eV. The scale for ETis on the right. Dashed line marks the phase transition [quantum critical point (QCP)] between the NM and AFM phases obtained at the level of Eq. (14). The estimated /ζ ratio for Ca2RuO4is indicated by shading. (and later include those of Fig. 3). As observed in Fig. 2(b), the variations in /ζ mainly modify the value of ET, which gets reduced with increasing /ζ , enabling the exchange interactions to induce a condensation of Tparticles at a certain point. In the extreme limit /ζ →∞, the exchange reaches the isotropic Heisenberg form and the sand Tx,ylevels merge. In this case the orbital angular momentum gets fully quenched and we deal with a real spin-1 model. C. Propagation of doped d5electronlike carriers in the d4background By introducing extra electrons into the d4system, we get a mixture of d4and d5ions with the d5configuration acting as a mobile electronlike carrier. The latter essentially moves via the regular nearest-neighbor hopping tas given by Eq. (8). When deriving the low-energy model, this motion has to be expressed in the low-energy basis composed of s,Tx,y, and f↑,↓states introduced in Sec. II A. The resulting d5-d4 coupling Hamiltonian encompasses the motion of doped d5 carriers represented by fermionic fparticles and a simultaneous “counterflow” in the d4background involving the sand Tx,ybosons. Due to the nontrivial structure of the background, we get a number a contributions to the fpropagation: apart FIG. 3. Propagation of d5states in the d4background. (a) Motion of doped electronlike d5carriers in the d4background originates in the hopping processes between the d5and d4ions. These may lead, among other possibilities, to a simple motion without any pseudospin change (top right) or generate (or similarly absorb) an excitation in the background accompanied by a simultaneous change of pseudospin state of the doped carrier (bottom right). (b) Parameters entering the d5-d4part of the model (13) calculated for the same setup as in Fig. 2(b). from the usual hopping insensitive to the pseudospin 1 2carried by f, there are also pseudospin-selective hoppings as well as pseudospin-flip terms [see Fig. 3(a) for an example], all of these accompanied by a properly arranged d4counterflow. A detailed description of the individual processes and their role will be given in Secs. VBand VCwhen analyzing the actual numerical results on fpropagation. To efficiently write the d5-d4coupling, we introduce a bond analog of the pseudospin 1 2carried by the doped electron, defined as σν ij = s,s=↑,↓ f† si σν ssfsj,(10) where σνof the right-hand side stands for the Pauli matrices. The ν=0 component of the bond operator σij containing the identity matrix σ0corresponds to the regular hopping of the form σ0 ij =(f† ↑if↑j+f† ↓if↓j), the others to the pseudospin-selective hopping (σz ij) and to pseudospin flips (σx,y ij ). Similarly, the counterflow in the d4background is best expressed using bond analogs of the pseudospin-1 operators:  Sν ji =−i(s† jTνi−T† νjsi)(ν=x,y) (11) and  Sz ji =−i(T† xjTyi −T† yjTxi ),(12) which have the same structure as in Eqs. (5)–(7)butthe two bosonic operators are now assigned to the two sites of a nearest-neighbor bond. Altogether, the d5-d4coupling, which splits naturally into a spin-spin and density-density channel, 195125-5 REVENDA, WOHLFELD, AND CHALOUPKA PHYSICAL REVIEW B 111, 195125 (2025) reads as Hd5-d4=− ijAxσx ij Sx ji +Ayσy ij Sy ji −Bσz ij Sz ji +σ0 ij[C0s† jsi+CxT† xjTxi +CyT† yjTyi ]+H.c.. (13) As in the case of Eq. (9), the xand ycomponents of interaction parameters are bond selective, namely, Ax=A±δA, Ay=A∓δAand Cx=C∓δC,Cy=C±δCwith the upper and lower signs being applied in the case of xand ybonds, respectively. The signs in Eq. (13) are adjusted in such a way that all the parameters are positive in the whole /ζ > 0 range. The values of the interactions are shown in Fig. 3(b),assuming again the microscopic parameter setup appropriate for Ca2RuO4as in Fig. 2(b). The corresponding expressions are summarized in Appendix B(the interactions in a simplified form for the /ζ =0 case can be also found in the initial study [71]). Interestingly, despite the small bond selectivity found in the case of d4Hamiltonian, here it is found to dominate the pseudospin-flip Aterms. In the /ζ →∞limit, the states s,Tx,yare degenerate and participate equally in the d5-d4processes, leading to |Ax|=|Ay|=B=C0=C. Finally, let us comment on the trivial case of the nonmagnetic background composed exclusively of ssinglets with T particles completely suppressed. In such a case, obtained in large-ETlimit, the d5-d4coupling reduces to only a “simple” hopping contained in the ∝C0term with C0being the hopping amplitude. III. MAGNETISM OF d4BACKGROUND: PHASES AND EXCITATIONS In this paragraph we review the phases and unusual excitations hosted by the (soft-spin) d4background driven by Hamiltonian (9). The d4model will be solved on the level of linear flavor-wave theory employed also in earlier studies [27,29]. In general, this approach is reminiscent of the linear spin-wave theory, for it expresses the local states in terms of bosons and handles the bond interactions among them and constraints in a similar approximative way. It was initially developed for anisotropic spin-1 systems [73], later applied also to spin-1 nematics [74,75] as well as in other situations, e.g., to bilayer Heisenberg magnets [76] that can be captured by singlet-triplet models within the bond operator approach [77,78]. In Ref. [79], the method has been used to study phases and excitations of a 2D S=1 spin system with planar anisotropy, formally similar to our d4model (9), and successfully benchmarked against quantum Monte Carlo simulations. The approximation quite well describes the transition between the quantum paramagnet at large planar anisotropy and the ordered planar antiferromagnet for a smaller planar anisotropy. Moreover, it captures in a unified way both the excitonic excitations in the quantum paramagnet as well as the Goldstone magnons and amplitude modes of the planar antiferromagnet [80]. Note, however, that this approach has certain limitations at the spin-1 Heisenberg point with vanishing planar anisotropy [80](verylarge/ζ limit in our case). Applying the above method to our model of Eq. (9), we first approximate its ground state by the variationally optimized product state |=  R∈sites [1−ρs†+√ρ(d∗ xT† x+d∗ yT† y)]R|vac. (14) This trial wave function allows for a condensate of Tparticles, having a density ρand the internal structure captured by the site-dependent complex vectors d=(dx,dy). The simple form of (14) transparently highlights the essential feature of our model: the possibility of a coherent onsite mixing of sand T. By minimizing the average of Hd4in the variational ground state (14), taking ρand the set of das variational parameters, the model is found to support two phases. At large ET compared to the exchange interactions, a nonmagnetic phase with ρ=0 is realized, i.e., the d4background state is made of ionic ground states s. Once the exchange interactions reach the critical strength Jcrit =1 8(ET−1 2), the system switches, via the condensation of Tbosons, to an antiferromagnetic phase with ordered in-plane van Vleck moments  Sx,y, hosted by onsite superposition of sand Tx,y[see Eqs. (5) and (6)]. In the latter case d∗=ieiQ·R(cos φ,sin φ),(15) where Qis the AF ordering vector Q=(π,π) and φis the angle of the ordered moments in the xy plane that is fixed by term. Taking Ca2RuO4as an example, in the following we choose the corresponding φ=π/4. The condensate density obtained by minimization takes the value ρ=1 2[1 −(ET− 1 2)/8J] and grows with increasing exchange strength (or decreasing ET)uptoρ=1 2in the J/ET→∞limit. The phase diagram constructed by employing the model parameters as given in Fig. 2(b) is presented in Fig. 4(a). For simplicity, we ignore the minor parameter δJ. The system is tuned by varying /ζ, which influences mainly ETthat becomes significantly reduced and reaches the critical value at about /ζ ≈0.9. To determine the excitation spectra, Hd4is first transformed by rotating the bosonic operators (s,Tx,Ty)toanew set (a,b,c) following the implicit relations sR=cos θcR+ieiQ·Rsin θaR, TaR=ieiQ·Rsin θcR+cos θaR, TbR=bR(16) with the linear combinations Ta=cos φTx+sin φTyand Tb=−sin φTx+cos φTyadjusted to match the in-plane pseudospin-1 components  S(parallel to the ordered moment direction) and  S⊥(perpendicular to it) appearing in Eq. (9). The main rotation angle θis determined by the condensate density via sin θ=√ρ. The rotation is designed so that the new bosons play physically distinct roles. Boson cis associated with the condensate as seen by obtaining |= Rc† R|vacafter the rotation; on the other hand, bosons a and bare associated with two different kinds of elementary excitations in the system. The linear flavor-wave expansion is then performed by substituting c,c†→√1−na−nb≈ 1−1 2(na+nb) to account for the hard-core constraint in a dynamic manner, and collecting terms up to second order in aand b. The result is a sum of two quadratic Hamiltonians, 195125-6 MAGNETIC POLARONS DUE TO SPIN-LENGTH … PHYSICAL REVIEW B 111, 195125 (2025) FIG. 4. Ground and excited states of the d4background. (a) Phase diagram of the d4model obtained by a variational calculation at the level of Eq. (14) for the parameter values from Fig. 2(b). Above the critical value crit ≈0.9ζ, the initially nonmagnetic system develops an AF order characterized by van Vleck magnetic moment  Sand T-condensate density ρ. (b) Excitation spectra obtained via linear flavor-wave theory for several values of measured relative to its critical value crit. (c) Schematic view of the two excitations αand β. Upon excitation by an amplitude mode α,the moment length of a given site oscillates in time, keeping its direction intact. In contrast, the magnon mode βcorresponds to pure rotations of the moments. In the long-wavelength limit, the modes can also be understood as specific oscillations in the energy landscape shown on the right (the anisotropy term is not considered in this sketch). each of them separately involving either aor bbosons. At this point neglecting the δJterm is very convenient because such a term would lead to a mixing of aand bmodes. Both aand b contributions are diagonalized by successive Fourier and Bogoliubov transformations, giving rise to αand βeigenmodes described by the Hamiltonian q(ωαqα† qαq+ωβqβ† qβq). At the level of linear flavor-wave theory, they appear as sharp excitations with infinite lifetime, a damping would occur when going beyond the quadratic expansion [29]. The dispersions of the eigenmodes are given by ωαq=A2 aq−B2 aqand ωβq= A2 bq−B2 bq, with the subfactors Aaq=ET−1 2+4Jγq,(17) Abq=ET+1 2+4Jγq,(18) Baq=Bbq=4Jγq(19) valid for the NM phase or Aaq=4J(2 +γqcos22θ),(20) Baq=4Jγqcos22θ, (21) Abq=4J(2 +γq) cos2θ−4Jzγqsin2θ+,(22) Bbq=4(Jcos2θ+Jzsin2θ)γq(23) to be used in the case of the AF phase. Here γqdenotes the nearest-neighbor factor γq=1 2(cos qx+cos qy) for the square lattice. The evolution of the excitation spectra with increasing value of /ζ , and thus at different points of the phase diagram, is presented in Fig. 4(b). Starting deep in the NM phase, the two excitations are seen as nearly degenerate (the small splitting is due to the anisotropy). When the critical point is approached, the excitations soften at the AF momentum Q with the αmode touching zero level at the critical value of /ζ. The further evolution (past the critical point) of the two modes is very different. The αmode hardens and becomes a flat high-energy excitation at large /ζ .Theβmode remains dispersive, gradually softening around q=(0,0), and keeping a nearly constant gap at Qof about √8J.This distinct behavior can be interpreted based on the physical picture of the two modes provided by Fig. 4(c). The mode αrelated to boson acorresponds to oscillations of the balance between sand T, i.e., the amplitude of the condensate, which explains its increased stiffness as the condensate develops and becomes more robust at larger /ζ. The mode βoriginating from bbosons can be understood as a magnon gapped by the in-plane anisotropy . The softening of its dispersion around q=(0,0) reflects the progression towards Heisenberg-type situation due to vanishing single-ion anisotropy ETcontained in the Hamiltonian Hd4. IV. A d5CARRIER IN THE d4BACKGROUND: FORMAL MATTERS After briefly inspecting the features of the unusual magnetic background, we focus on the interactions of a doped electron f(corresponding locally to a d5ion) with the elementary magnetic excitations of the background. The resulting polaronic behavior of the doped electron is captured within the self-consistent Born approximation (SCBA) [43,57,70]. The required d5-d4interaction Hamiltonian is obtained by starting with Hd5-d4of Eq. (13) and applying the same steps as in the 195125-7 REVENDA, WOHLFELD, AND CHALOUPKA PHYSICAL REVIEW B 111, 195125 (2025) linear flavor-wave expansion, i.e., bosonic rotation (16) with a replacement of the condensed boson cby √1−na−nb, followed by an expansion in the aand bbosons and subsequent Fourier transformation. Finally, the magnetic bosons a and bare expressed in terms of the eigenmodes αand βas aq=uαqαq+vαqα† −qwith the Bogoliubov factors uαq=1 √2Aaq ωαq+1,vαq=1 √2Aaq ωαq−1sgnBaq (24) and similar for band β.TheAand Bfactors are listed in Eqs. (17)–(23). Two parts of the resulting coupling Hamiltonian are of interest here: (i) Zeroth-order terms in αand β corresponding to a “simple” electron hopping without any excitations of the background involved. However, this hopping is still affected by the very presence of the magnetic condensate blocking the fmotion. (ii) Terms linear in αand βproviding the basic coupling of the magnetic excitations to the doped electrons. In the following, we discuss the latter coupling in detail and describe how it gets incorporated into the SCBA scheme. We consider first the simple yet instructive case of the nonmagnetic phase before moving on to the antiferromagnetic one, which is of our main interest but leads to a rather complex formal structure. A. SCBA equations for the nonmagnetic phase In the nonmagnetic phase, the elementary excitations α,β are associated with the Tparticles themselves and the relevant zerothand first-order processes are therefore directly those illustrated in Fig. 3(a). Zeroth-order interaction term, which does not involve any T, corresponds to a “free” fhopping through the nonmagnetic d4background composed of s. According to Eq. (13), its amplitude equals to C0, which leads to the free motion given by Hamiltonian H(0) d5-d4= ks εkf† ksfks(25) with the bare dispersion relation of a doped electron: εk=−2C0(cos kx+cos ky).(26) First-order terms, on the other hand, involve a creation or annihilation of the Texcitations. Since T0has been removed by the tetragonal crystal field and the active Tx,yare based on T±1, the corresponding processes always include a flip of the electron pseudospin 1 2as depicted in Fig. 3(a). If, for the sake of brevity, we ignore the in-plane anisotropy in this paragraph, the αand βexcitations are degenerate, with a common dispersion ωαq=ωβq=ωqand identical Bogoliubov factors uαq=uβq=uq,vαq=vβq=vq. With these simplifications, the first-order coupling of doped electrons f to the magnetic excitations α,βis given by H(1) d5-d4= kqsMα kqsα† q+Mβ kqsβ† qf† k−q,−sfks+H.c.(27) with the matrix element Mα,β kqs∝δA(ηk−quq−ηkvq)±A(γk−quq−γkvq).(28) FIG. 5. (a) Diagrammatic representation of the coupling Hamiltonians (27)and(33). There is always a term flipping pseudospin of the doped electron upon emission or absorption of magnetic excitations α,β. In the AF phase, a new contribution appears that is pseudospin conserving but brings Q=(π,π) momentum shift. This contribution scales with condensate density as √ρ. (b) Dyson’s equations for the normal propagator of the doped electron and the anomalous one, involving a pseudospin flip and a simultaneous momentum shift Q. Correspondingly, two kinds of self-energies enter these Dyson’s equations. The input and output momenta and pseudospins are indicated at the top of the diagrams. The simple line corresponds to the only bare propagator G0=1/(E−εk). In the nonmagnetic phase, the anomalous propagator and self-energy are absent and the Dyson’s equations reduce to the first line with the last term omitted. (c) Self-energy contributions within SCBA for the AF phase. The wiggly line represents αand βexcitations (to be summed over) emitted and absorbed during electron motion. The internal electron lines corresponding to normal propagation are labeled by their momenta and pseudospins. In the nonmagnetic phase, the only nonvanishing term is the first contribution to the normal self-energy . Here, the ±sign applies to αand β, respectively, and the omitted proportionality factor equals −4ifor α,−4forβand s=↑, and +4forβand s=↓. The geometry of the hopping enters the matrix element via sand d-wave nearest-neighbor form factors for the square lattice: γk=1 2(cos kx+cos ky) and ηk=1 2(cos kx−cos ky). The renormalized electron propagator Gk(E)=[E−εk− k(E)]−1is obtained at the level of SCBA, summarized in Fig. 5for the more general case of the AF phase. In the nonmagnetic phase, the SCBA self-energy reduces to the first diagram in Fig. 5(c) that translates to a convolution of the renormalized electron propagator and propagators of the magnetic excitations. Note that even though the interaction in 195125-8 MAGNETIC POLARONS DUE TO SPIN-LENGTH … PHYSICAL REVIEW B 111, 195125 (2025) Eq. (27) flips the electron pseudospin, the balance between emitted and absorbed magnetic excitations during the electron motion restores the initial pseudospin in the electron propagator, which is then pseudospin conserving as well as pseudospin independent. For a convenient evaluation, we decompose the SCBA self-energy into real and imaginary parts =+i. The imaginary part  of the self-energy is calculated via  k(E)=−π q|Mkq|2Ak−q(E−ωq) (29) with the electron spectral function given by Ak(E)= −π−1Im Gk(E) and |Mkq|2denoting summed up αand β contributions |Mα kqs|2+|Mβ kqs|2: 32[(δA)2(ηk−quq−ηkvq)2+A2(γk−quq−γkvq)2].(30) The real part of the self-energy is subsequently obtained by Kramers-Kronig transformation  ks(E)=1 πP+∞ −∞  ks(ξ) E−ξdξ(31) and the steps are repeated until a self-consistent solution of the self-energy equation is found. B. SCBA equations for the antiferromagnetic phase The case of the antiferromagnetic phase is formally richer since the now nontrivial bosonic rotation (16) with θ= 0 and φ=π/4 brings up a number of new terms to the interaction. The contributions of zeroth order in α,βagain give rise to H(0) d5-d4of Eq. (25) with the bare electron dispersion εktaking this time the form εk=−2[C0−(C0+C)sin 2θ](cos kx+cos ky) −2δCcos 2φ(cos kx−cos ky).(32) The second term vanishes in our case of interest due to φ= π/4, the first one shows a reduction of the bare bandwidth compared to the nonmagnetic phase, the bare motion of the doped electron in the staggered magnetic structure of the AF condensate being increasingly inhibited with growing condensate density ρ=sin2θ. The linear coupling to α,β, represented diagrammatically in Fig. 5(a), has the structure H(1) d5-d4= kqsMα kqsα† q+Mβ kqsβ† qf† k−q,−sfks +¯ Mα kqsα† q+¯ Mβ kqsβ† qf† k−q+Q,sfks+H.c.(33) The matrix elements Mα,β kqsand ¯ Mα,β kqsare given in Appendix C. Compared to the nonmagnetic phase, the coupling is extended by pseudospin-conserving terms, which include a momentum shift by the AF wave vector Q=(π,π). The combination of pseudospin-flipping and pseudospin-conserving processes gives rise to a more complex structure of Dyson’s equations including additionally anomalous propagators of the type  Gk↑∼fk+Q↓f† k↑[see Fig. 5(b)]. They involve normal self-energy ks(E) as well as anomalous self-energy  ks(E), which is associated with both a pseudospin-flip and Qmomentum shift. By solving Dyson’s equations in Fig. 5(b),the normal electron propagator is obtained in the form Gk↑=E−εk−k↑− k↑ k+Q↓ E−εk+Q−k+Q↓−1 ,(34) where the energy arguments were omitted for compactness of the expression. The anomalous propagator takes the form  Gk↑=(E−εk−k↑)(E−εk+Q−k+Q↓)  k↑− k+Q↓−1 . (35) The self-energies in SCBA approximation include all combinations of the interaction vertices as depicted in Fig. 5(c).The  part of the normal self-energy reads as  ks(E) =−π qMα kqs2Ak−q,−s+¯ Mα kqs2Ak−q+Q,s +Mα kqs¯ Mα kqs∗ Ak−q,−s+¯ Mα kqsMα kqs∗ Ak−q+Q,sE−ωαq +βterms of the same structure ,(36) with Aand  Abeing the spectral functions corresponding to the normal and anomalous propagators Gand  G, respectively, all of them having the energy argument E−ωαq(αterms) or E−ωβq(βterms). Similarly, the   part of the anomalous self-energy can be expressed as   ks(E) =π qMα kqs¯ Mα kq,−s∗Ak−q,−s+¯ Mα kqsMα kq,−s∗Ak−q+Q,s +Mα kqsMα kq,−s∗ Ak−q,−s+¯ Mα kqs¯ Mα kq,−s∗ × Ak−q+Q,sE−ωαq+βterms .(37) The self-energy expressions (36) and (37) can be further simplified by using the symmetry relations Aks=Ak,−sfor the normal components and  Aks= A∗ k,−s, Aks= Ak+Q,sfor the anomalous ones. As in the case of the nonmagnetic phase, the system of Dyson’s equations and self-energy equations is solved by iterations until self-consistency is reached. V. A d5CARRIER IN THE d4BACKGROUND: NUMERICAL RESULTS Utilizing the SCBA approach described in the previous section, we are able to study the propagation of doped electrons in the unusual soft-spin background. The aim of this section is to present the numerical results and analyze the interplay of doped electrons with the excitations of the magnetic background, focusing in particular on the AF-ordered phase hosting both amplitude and magnon excitations. We will contrast the resulting polaronic behavior to spin polarons thoroughly discussed in context of regular spin-1 2Heisenberg antiferromagnets. We start by presenting in Fig. 6(a) the overall trends in the normal spectral function Ak(E) through the magnetic phase diagram. The parameter values used are those presented earlier in Figs. 2(b) and 3(b), with the exception of δJ.This 195125-9 REVENDA, WOHLFELD, AND CHALOUPKA PHYSICAL REVIEW B 111, 195125 (2025) contribution to the single-particle propagator decreases (since their “free-particle character” gets reduced). Small V : The hidden pseudogap limit. At V∼0.2the spectrum in Fig. 11(b) seems to contain merely a free dispersion. However, a closer look at the underlying eigenvalues reveals that the situation is more complex [see the left panel of Fig. 11(c)]: just as in the moderate Vregime there is a pseudogap in the free dispersion. The latter arises due to the fact that again this “free band” is actually formed by two eigenstates: the ground and the excited state that switch their characters depending on momentum. As a consequence, either the ground or the excited state shows up in the spectral function of cfor a given momentum. To sum up, the employed toy model shows that, as a result of the interplay between the free hopping and stringlike motion, a pseudogap, or a vanishing noninteracting quasiparticle spectral weight, may occur in the free dispersion. Interestingly, this happens not only for moderate ratios of the free dispersion with respect to the stringlike motion but even once the latter type of motion is vanishingly small. This phenomenon can also be verified by inspecting properties of the propagator (43): The quasiparticle spectral weight Z=1/(1 −∂Re/∂E|E=E) vanishes at E=Eudue to ∂Re/∂E|E=E=∞for any finite value of V. Altogether, the similarity between the behavior of the hedgehog toy model, which essentially is a noninteracting one, and the full polaronic model means that the pseudogaplike feature may be understood within a single-band picture as an effect of a scattering on a single-energy level. Formally, this correspondence stems from the rather crude, yet apparently viable, single-pole approximation to the self-energy of the full model. VII. CONCLUSIONS To conclude we have studied the motion of a single doped electron introduced into the 2D Mott insulator with d4valence shell that is subject to moderate to strong spin-orbit coupling. In this situation, the entanglement of spin and orbital degrees of freedom becomes essential for both the d5-doped carrier as well as for the d4-local moments in the magnetic background. We have formulated the corresponding exchange model for the undoped case by extending earlier works, which show that the interaction between the magnetic moments in the case of the prototypical d4spin-orbit Mott insulator Ca2RuO4can be described in terms of a 2D XYZ model for an effective pseudospin  S=1. The essential part of the model is strong planar single-ion anisotropy microscopically determined by the interplay of the spin-orbit coupling ζand tetragonal field splitting among the t2gorbitals [29]. Using /ζ as a control parameter, the model can be driven through a quantum critical point separating its two phases: the nonmagnetic (NM) phase at smaller /ζ and the antiferromagnetic (AF) phase for larger /ζ . The latter can be associated with a condensation of the excitoniclike magnetic excitations hosted by the NM phase and features a Goldstone magnon and an amplitude mode in its excitation spectrum. When a mobile d5carrier is introduced into such a magnetic background, its propagation proceeds via a rather complex set of processes, which include processes associated with a direct creation or annihilation of a magnetic excitation as well as those that apparently leave the background intact. Here we analyzed the competition of these two classes of processes in great detail and uncovered their respective role in the renormalization of the doped electron across the model phase diagram. Experimentally, the studied model at /ζ ∼1.5–2.0, i.e., moderately deep in the AF phase, describes the inverse photoemission performed on Ca2RuO4. The two main results of this work are as follows: First, both on the NM side of the quantum critical point as well as just past it, i.e., for the weakly ordered AF case, the doped carrier largely behaves as a free quasiparticle. In this regime the electron propagates via a combination of free hopping, which does not disturb the d4background built mostly of ionic singlets, and scattering processes originating in the spin-spin channel, i.e., from the coupling of pseudospin 1 2carried by the electron to the pseudospins  S=1 of the background. This scattering employs, without a particular preference, all kinds of magnetic excitations. It operates through the entire phase diagram, but in the above case and within our parameter regime, it is by itself not sufficient to induce a strong polaronic behavior of the doped carrier. Nevertheless, the coupling to the magnetic excitations is still visible in the spectrum. Namely, at particular values of momenta a pseudogap develops in the quasiparticle dispersion and this can be attributed to the onset of a stringlike polaronic motion. We have understood the interplay of the latter and the free dispersion in quite some detail using a (hedgehog) toy model. Second, deeper in the AF phase the character of the spectrum qualitatively changes and strong ladderlike features, that are signatures of the onset of the string potential [40,42,46], are visible. While this result describes the realistic case of Ca2RuO4, from the general point of view it is interesting in itself to see such a strong onset of the string potential in the model. Our detailed analysis shows that the strong string potential arises here due to the dominant coupling of the doped electron to the amplitude mode in the density-density channel. This type of scattering starts upon entry into the AF phase with the appearance of a magnetic condensate and steeply rises afterwards, following the growth of the condensate amplitude. In a broader context, the AF-phase results are somewhat reminiscent of the observed coupling of the doped carrier to the spin fluctuations deep in the antiferromagnetic phase of the bilayer Heisenberg magnets [81–85]. These are described by singlet-triplet models that share certain features with our model, such as the magnetic condensation and the structure of excitation spectra, though the details of the respective models are different. On the other hand, we note that our case is distinct from the ruthenate photoemission problem modeled in [70]. There, the coupling to the amplitude mode was neglected, for only the magnon mode of Ref. [88] was considered. Nevertheless, the ladder spectrum and the string potential was also observed, being traced back to the strong magnetic and hopping anisotropies of the t-J model used. 195125-16 MAGNETIC POLARONS DUE TO SPIN-LENGTH … PHYSICAL REVIEW B 111, 195125 (2025) ACKNOWLEDGMENTS We would like to thank G. Khaliullin for useful discussions. J.R. and J.C. acknowledge support by Czech Science Foundation (GA ˇ CR) under Project No. GA22-28797S and by the project Quantum Materials for Applications in Sustainable Technologies, Grant No. CZ.02.01.01/00/22_008/0004572. Computational resources were provided by the e-INFRA CZ project (Project ID No. 90254), supported by the Ministry of Education, Youth and Sports of the Czech Republic. K.W. thanks IIT Madras for an IoE Visiting Scientist position which enabled the completion of this work. DATA AVAILABILITY The data that support the findings of this article are openly available [89]. APPENDIX A: MICROSCOPIC FORMULAS FOR d4 EXCHANGE PARAMETERS The exchange parameters entering the effective pseudospin-1 model (9) that drives the d4magnetic background are given explicitly by J=7−c0+c1−7c0c1−2√2s0s1 48(1 −3η) +13−c0−5c1+29c0c1+13√2s0s1 96 +9−c0−c1+5c0c1+3√2s0s1 32(1 +2η),(A1) δJ=5−3c0+3c1−5c0c1+2√2s0s1 48(1 −3η) +−10−6c1+4c0c1−7√2s0s1 96 −2c0+2c0c1+√2s0s1 32(1 +2η),(A2) Jz=1−cos 4ϑ1 24(1 −3η)+17 −24c1+7 cos 4ϑ1 96 +7−8c1+cos 4ϑ1 32(1 +2η).(A3) The respective interaction constants are given in units of t2/Uand using a shorthand notation c0,1=cos 2ϑ0,1,s0,1= sin 2ϑ0,1, and η=JH/U.HereUdenotes the usual intraorbital Hubbard repulsion and JHHund’s coupling. APPENDIX B: MICROSCOPIC FORMULAS FOR d5-d4 INTERACTION PARAMETERS Interaction parameters of the coupling between d5electronlike carriers and d4background used in Eq. (13) may be compactly written as follows: A=[2c2c0−(cs0+sc0)s]1 √2c1,(B1) δA=(cs0+sc0)cs1+1 √2sc1,(B2) B=cs2 1−cc2 1+√2sc1s1c,(B3) C0=c2+c2 0+2csc0s0,(B4) C=1−s2s2 1+√2csc1s1,(B5) δC=cs1+1 √2sc1√2sc1.(B6) They are given in units of t/4 and expressed as polynomials in the wave-function factors c=cos ϑ,s=sin ϑ,c0=cos ϑ0, s0=sin ϑ0,c1=cos ϑ1, and s1=sin ϑ1. APPENDIX C: MATRIX ELEMENTS OF THE COUPLING BETWEEN fAND ELEMENTARY EXCITATIONS α,β The matrix elements entering the linear interaction of α, βelementary excitations with the doped electrons in Eq. (33) read as Mα kqs=4iWαs k−quαq−Wαs kvαq, Mβ kqs=4iWβs k−quβq−Wβs kvβq, ¯ Mα kqs=4i¯ Wα k−quαq+¯ Wα kvαq, ¯ Mβ kqs=4i¯ Wβs 1k−quβq+¯ Wβs 2kvβq.(C1) They combine the Bogoliubov factors uαq,vαq,uβq,vβq[see Eq. (24)] together with Wαs k=δAηk(−cos φ+is sin φ)−Aγk(cos φ+is sin φ), Wβs k=[δAηk(sin φ+is cos φ)+Aγk(sin φ−is cos φ)] ×cos θ, ¯ Wα k=[−(C0+C)γk+δCηkcos 2φ]sinθcos θ, ¯ Wβs 1k=(+isB γk−δCηksin 2φ)sinθ, ¯ Wβs 2k=(−isB γk−δCηksin 2φ)sinθ. (C2) In these formulas, shas to be understood as s=±1for pseudospin ↑,↓, respectively. The momentum dependence of the matrix elements is captured using the s-wave and d-wave nearest-neighbor factors for the square lattice: γk= 1 2(cos kx+cos ky) and ηk=1 2(cos kx−cos ky). Note that all ¯ Wand consequently the ¯ Mmatrix elements corresponding to pseudospin-conserving terms in Eq. (33) contain sin θmaking them to disappear in the nonmagnetic phase and to grow as √ρupon entering the AF phase. The simplified formulas for the nonmagnetic phase and =0 that are relevant to Eq. (27) are obtained by setting θ=0, φ=0. 195125-17 REVENDA, WOHLFELD, AND CHALOUPKA PHYSICAL REVIEW B 111, 195125 (2025) [1] K. I. Kugel and D. I. Khomskii, The Jahn-Teller effect and magnetism: transition metal compounds, Usp. Fiz. Nauk 136, 621 (1982) [Sov. Phys.–Usp. 25, 231 (1982)]. [2] Y. Tokura and N. Nagaosa, Orbital physics in transition-metal oxides, Science 288, 462 (2000). [3] F. Moussa, M. Hennion, J. Rodriguez-Carvajal, H. Moudden, L. Pinsard, and A. Revcolevschi, Spin waves in the antiferromagnet perovskite LaMnO3: A neutron-scattering study, Phys. Rev. B54, 15149 (1996). [4] A. M. Ole´ s, L. F. Feiner, and J. Zaanen, Quantum melting of magnetic long-range order near orbital degeneracy: Classical phases and Gaussian fluctuations, Phys.Rev.B61, 6257 (2000). [5] G. Khaliullin and S. Okamoto, Theory of orbital state and spin interactions in ferromagnetic titanates, Phys. Rev. B 68, 205109 (2003). [6] K. Wohlfeld, M. Daghofer, S. Nishimoto, G. Khaliullin, and J. van den Brink, Intrinsic coupling of orbital excitations to spin fluctuations in Mott insulators, Phys.Rev.Lett.107, 147201 (2011). [7] J. Schlappa, K. Wohlfeld, K. J. Zhou, M. Mourigal, M. W. Haverkort, V. N. Strocov, L. Hozoi, C. Monney, S. Nishimoto, S. Singh, A. Revcolevschi, J.-S. Caux, L. Patthey, H. M. Rønnow, J. van den Brink, and T. Schmitt, Spin–orbital separation in the quasi-one-dimensional Mott insulator Sr2CuO3, Nature (London) 485, 82 (2012). [8] L. Martinelli, K. Wohlfeld, J. Pelliciari, R. Arpaia, N. B. Brookes, D. Di Castro, M. G. Fernandez, M. Kang, Y. Krockenberger, K. Kummer, D. E. McNally, E. Paris, T. Schmitt, H. Yamamoto, A. Walters, K.-J. Zhou, L. Braicovich, R. Comin, M. M. Sala, T. P. Devereaux et al., Collective nature of orbital excitations in layered cuprates in the absence of apical oxygens, Phys. Rev. Lett. 132, 066004 (2024). [9] G. Jackeli and G. Khaliullin, Mott insulators in the strong spinorbit coupling limit: From Heisenberg to a quantum compass and Kitaev models, Phys.Rev.Lett.102, 017205 (2009). [10] W. Witczak-Krempa, G. Chen, Y. B. Kim, and L. Balents, Correlated quantum phenomena in the strong spin-orbit regime, Annu. Rev. Condens. Matter Phys. 5, 57 (2014). [11] T. Takayama, J. Chaloupka, A. Smerald, G. Khaliullin, and H. Takagi, Spin–orbit-entangled electronic phases in 4dand 5dtransition-metal compounds, J. Phys. Soc. Jpn. 90, 062001 (2021). [12] A. Kitaev, Anyons in an exactly solved model and beyond, Ann. Phys. 321, 2 (2006). [13] S. M. Winter, A. A. Tsirlin, M. Daghofer, J. van den Brink, Y. Singh, P. Gegenwart, and R. Valentí, Models and materials for generalized Kitaev magnetism, J. Phys.: Condens. Matter 29, 493002 (2017). [14] H. Takagi, T. Takayama, G. Jackeli, G. Khaliullin, and S. E. Nagler, Concept and realization of Kitaev quantum spin liquids, Nat. Rev. Phys. 1, 264 (2019). [15] Y. Matsuda, T. Shibauchi, and H.-Y. Kee, Kitaev quantum spin liquids, arXiv:2501.05608. [16] S. Nakatsuji, S.-i. Ikeda, and Y. Maeno, Ca2RuO4:NewMott insulators of layered ruthenate, J. Phys. Soc. Jpn. 66, 1868 (1997). [17] M. Braden, G. André, S. Nakatsuji, and Y. Maeno, Crystal and magnetic structure of Ca2RuO4: Magnetoelastic coupling and the metal-insulator transition, Phys.Rev.B58, 847 (1998). [18] C. S. Alexander, G. Cao, V. Dobrosavljevic, S. McCall, J. E. Crow, E. Lochner, and R. P. Guertin, Destruction of the Mott insulating ground state of Ca2RuO4by a structural transition, Phys. Rev. B 60, R8422 (1999). [19] E. Gorelov, M. Karolak, T. O. Wehling, F. Lechermann, A. I. Lichtenstein, and E. Pavarini, Nature of the Mott transition in Ca2RuO4,Phys.Rev.Lett.104, 226401 (2010). [20] G. Cao, S. McCall, M. Shepard, J. E. Crow, and R. P. Guertin, Magnetic and transport properties of single-crystal Ca2RuO4: Relationship to superconducting Sr2RuO4,Phys. Rev. B 56, R2916 (1997). [21] M. Horio, F. Forte, D. Sutter, M. Kim, C. G. Fatuzzo, C. E. Matt, S. Moser, T. Wada, V. Granata, R. Fittipaldi, Y. Sassa, G. Gatti, H. M. Rønnow, M. Hoesch, T. K. Kim, C. Jozwiak, A. Bostwick, E. Rotenberg, I. Matsuda, A. Georges et al., Orbital-selective metal skin induced by alkali-metaldosing Mott-insulating Ca2RuO4,Commun. Phys. 6, 323 (2023). [22] S. Riccò, M. Kim, A. Tamai, S. McKeown Walker, F. Y. Bruno, I. Cucchi, E. Cappelli, C. Besnard, T. K. Kim, P. Dudin, M. Hoesch, M. J. Gutmann, A. Georges, R. S. Perry, and F. Baumberger, In situ strain tuning of the metal-insulatortransition of Ca2RuO4in angle-resolved photoemission experiments, Nat. Commun. 9, 4535 (2018). [23] C. T. Suen, I. Markovi´ c, M. Zonno, N. Heinsdorf, S. Zhdanovich, N. H. Jo, M. Schmid, P. Hansmann, P. Puphal, K. Fürsich, S. Smit, C. Au-Yeung, V. Zimmermann, B. Zwartsenberg, M. Krautloher, I. S. Elfimov, R. Koch, S. Gorovikov, C. Jozwiak, A. Bostwick et al., Electronic response of a Mott insulator at a current-induced insulator-to-metal transition, Nat. Phys. 20, 1757 (2024). [24] T. Mizokawa, L. H. Tjeng, G. A. Sawatzky, G. Ghiringhelli, O. Tjernberg, N. B. Brookes, H. Fukazawa, S. Nakatsuji, and Y. Maeno, Spin-orbit coupling in the Mott insulator Ca2RuO4, Phys. Rev. Lett. 87, 077202 (2001). [25] M. Cuoco, F. Forte, and C. Noce, Probing spin-orbital-lattice correlations in 4d4systems, Phys.Rev.B73, 094428 (2006). [26] G. Khaliullin, Excitonic magnetism in Van Vleck-type d4Mott insulators, Phys.Rev.Lett.111, 197201 (2013). [27] A. Akbari and G. Khaliullin, Magnetic excitations in a spinorbit-coupled d4Mott insulator on the square lattice, Phys. Rev. B90, 035137 (2014). [28] S. Kunkemöller, D. Khomskii, P. Steffens, A. Piovano, A. A. Nugroho, and M. Braden, Highly anisotropic magnon dispersion in Ca2RuO4: Evidence for strong spin orbit coupling, Phys. Rev. Lett. 115, 247201 (2015). [29] A. Jain, M. Krautloher, J. Porras, G. H. Ryu, D. P. Chen, D. L. Abernathy, J. T. Park, A. Ivanov, J. Chaloupka, G. Khaliullin, B. Keimer, and B. J. Kim, Higgs mode and its decay in a twodimensional antiferromagnet, Nat. Phys. 13, 633 (2017). [30] S. Kunkemöller, E. Komleva, S. V. Streltsov, S. Hoffmann, D. I. Khomskii, P. Steffens, Y. Sidis, K. Schmalzl, and M. Braden, Magnon dispersion in Ca2Ru1−xTixO4: Impact of spinorbit coupling and oxygen moments, Phys.Rev.B95, 214408 (2017). [31] L. Das, F. Forte, R. Fittipaldi, C. G. Fatuzzo, V. Granata, O. Ivashko, M. Horio, F. Schindler, M. Dantz, Y. Tseng, D. E. McNally, H. M. Rønnow, W. Wan, N. B. Christensen, J. Pelliciari, P. Olalde-Velasco, N. Kikugawa, T. Neupert, A. Vecchione, T. Schmitt et al., Spin-orbital excitations in 195125-18 MAGNETIC POLARONS DUE TO SPIN-LENGTH … PHYSICAL REVIEW B 111, 195125 (2025) Ca2RuO4revealed by resonant inelastic x-ray scattering, Phys. Rev. X 8, 011048 (2018). [32] K. von Arx, P. Rothenbühler, Q. Wang, J. Choi, M. GarciaFernandez, S. Agrestini, K.-J. Zhou, A. Vecchione, R. Fittipaldi, Y. Sassa, M. Cuoco, F. Forte, and J. Chang, Resolving the orbital character of low-energy excitations in Mott insulator with intermediate spin-orbit coupling, arXiv:2312.04406. [33] G. Zhang and E. Pavarini, Higgs mode and stability of xy-orbital ordering in Ca2RuO4,Phys. Rev. B 101, 205128 (2020). [34] T. Feldmaier, P. Strobel, M. Schmid, P. Hansmann, and M. Daghofer, Excitonic magnetism at the intersection of spin-orbit coupling and crystal-field splitting, Phys.Rev.Res.2, 033201 (2020). [35] P. Strobel, F. Aust, and M. Daghofer, Magnetic phases for strongly correlated t4 2gelectrons on the square lattice: Impact of spin-orbit coupling and crystal field, Phys. Rev. B 104, 115148 (2021). [36] S.-M. Souliou, J. Chaloupka, G. Khaliullin, G. Ryu, A. Jain, B. J. Kim, M. Le Tacon, and B. Keimer, Raman scattering from Higgs mode oscillations in the two-dimensional antiferromagnet Ca2RuO4,Phys. Rev. Lett. 119, 067201 (2017). [37] P. A. Lee, N. Nagaosa, and X.-G. Wen, Doping a Mott insulator: Physics of high-temperature superconductivity, Rev. Mod. Phys. 78, 17 (2006). [38] M. Imada, A. Fujimori, and Y. Tokura, Metal-insulator transitions, Rev. Mod. Phys. 70, 1039 (1998). [39] D. I. Khomskii, Basic Aspects of the Quantum Theory of Solids: Order and Elementary Excitations (Cambridge University Press, Cambridge, 2010). [40] L. N. Bulaevskii, E. L. Nagaev, and D. I. Khomskii, A new type of auto-localized state of a conduction electron in an antiferromagnetic semiconductor, Zh. Eksp. Teor. Fiz. 54, 1562 (1968) [J. Exp. Theor. Phys. 27, 836 (1968)]. [41] S. Schmitt-Rink, C. M. Varma, and A. E. Ruckenstein, Spectral function of holes in a quantum antiferromagnet, Phys.Rev.Lett. 60, 2793 (1988). [42] C. L. Kane, P. A. Lee, and N. Read, Motion of a single hole in a quantum antiferromagnet, Phys. Rev. B 39, 6880 (1989). [43] G. Martínez and P. Horsch, Spin polarons in the t-Jmodel, Phys.Rev.B44, 317 (1991). [44] E. Manousakis, String excitations of a hole in a quantum antiferromagnet and photoelectron spectroscopy, Phys. Rev. B 75, 035106 (2007). [45] F. Grusdt, M. Kánasz-Nagy, A. Bohrdt, C. S. Chiu, G. Ji, M. Greiner, D. Greif, and E. Demler, Parton theory of magnetic polarons: Mesonic resonances and signatures in dynamics, Phys. Rev. X 8, 011046 (2018). [46]K.Bieniasz,P.Wrzosek,A.M.Ole ´ s,andK.Wohlfeld,From “weak” to “strong” hole confinement in a Mott insulator, SciPost Phys. 7, 066 (2019). [47] P. Wrzosek and K. Wohlfeld, Hole in the two-dimensional Ising antiferromagnet: Origin of the incoherent spectrum, Phys. Rev. B103, 035113 (2021). [48] C. S. Chiu, G. Ji, A. Bohrdt, M. Xu, M. Knap, E. Demler, F. Grusdt, M. Greiner, and D. Greif, String patterns in the doped Hubbard model, Science 365, 251 (2019). [49] J. Koepsell, J. Vijayan, P. Sompet, F. Grusdt, T. A. Hilker, E. Demler, G. Salomon, I. Bloch, and C. Gross, Imaging magnetic polarons in the doped Fermi-Hubbard model, Nature (London) 572, 358 (2019). [50] J. Koepsell, D. Bourgund, P. Sompet, S. Hirthe, A. Bohrdt, Y. Wang, F. Grusdt, E. Demler, G. Salomon, C. Gross, and I. Bloch, Microscopic evolution of doped Mott insulators from polaronic metal to Fermi liquid, Science 374, 82 (2021). [51] B. O. Wells, Z.-X. Shen, A. Matsuura, D. M. King, M. A. Kastner, M. Greven, and R. J. Birgeneau, Eversus krelations and many body effects in the model insulating copper oxide Sr2CuO2Cl2,Phys.Rev.Lett.74, 964 (1995). [52] F. Ronning, K. M. Shen, N. P. Armitage, A. Damascelli, D. H. Lu, Z.-X. Shen, L. L. Miller, and C. Kim, Anomalous highenergy dispersion in angle-resolved photoemission spectra from the insulating cuprate Ca2CuO2Cl2,Phys.Rev.B71, 094518 (2005). [53] B. Bacq-Labreuil, C. Fawaz, Y. Okazaki, Y. Obata, H. Cercellier, P. LeFevre, F. Bertran, D. Santos-Cottin, H. Yamamoto, I. Yamada, M. Azuma, K. Horiba, H. Kumigashira, M. d’Astuto, S. Biermann, and B. Lenz, Universal waterfall feature in cuprate superconductors: Evidence of a momentumdriven crossover, Phys.Rev.Lett.134, 016502 (2025). [54] R. Kilian and G. Khaliullin, Orbital polarons in the metalinsulator transition of manganites, Phys. Rev. B 60, 13458 (1999). [55] J. Chaloupka and G. Khaliullin, Spin polaron theory for the photoemission spectra of layered cobaltates, Phys. Rev. Lett. 99, 256406 (2007). [56] K. Wohlfeld, A. M. Ole´ s, and P. Horsch, Orbitally induced string formation in the spin-orbital polarons, Phys. Rev. B 79, 224433 (2009). [57] E. M. Pärschke, K. Wohlfeld, K. Foyevtsova, and J. van den Brink, Correlation induced electron-hole asymmetry in quasitwo-dimensional iridates, Nat. Commun. 8, 686 (2017). [58] F. Trousselet, M. Berciu, A. M. Ole´ s, and P. Horsch, Hidden quasiparticles and incoherent photoemission spectra in Na2IrO3,Phys.Rev.Lett.111, 037205 (2013). [59] G. B. Halász, J. T. Chalker, and R. Moessner, Doping a topological quantum spin liquid: Slow holes in the Kitaev honeycomb model, Phys. Rev. B 90, 035145 (2014). [60] B.-B. Wang, W. Wang, S.-L. Yu, and J.-X. Li, Dynamics of a single hole in the Heisenberg–Kitaev model: a self-consistent Born approximation study, J. Phys.: Condens. Matter 30, 385602 (2018). [61] J. Strasdas, B. Pestka, M. Rybak, A. K. Budniak, N. Leuth, H. Boban, V. Feyer, I. Cojocariu, D. Baranowski, J. Avila, P. Dudin, A. Bostwick, C. Jozwiak, E. Rotenberg, C. Autieri, Y. Amouyal, L. Plucinski, E. Lifshitz, M. Birowska, and M. Morgenstern, Electronic band structure changes across the antiferromagnetic phase transition of exfoliated MnPS3flakes probed by μ-ARPES, Nano Lett. 23, 10342 (2023). [62] B. J. Kim, H. Jin, S. J. Moon, J.-Y. Kim, B.-G. Park, C. S. Leem, J. Yu, T. W. Noh, C. Kim, S.-J. Oh, J.-H. Park, V. Durairaj, G. Cao, and E. Rotenberg, Novel Jeff =1/2 Mott state induced by relativistic spin-orbit coupling in Sr2IrO4,Phys. Rev. Lett. 101, 076402 (2008). [63] A. de la Torre, S. McKeown Walker, F. Y. Bruno, S. Riccó, Z. Wang, I. Gutierrez Lezama, G. Scheerer, G. Giriat, D. Jaccard, C. Berthod, T. K. Kim, M. Hoesch, E. C. Hunter, R. S. Perry, A. Tamai, and F. Baumberger, Collapse of the Mott gap and emergence of a nodal liquid in lightly doped Sr2IrO4,Phys. Rev. Lett. 115, 176402 (2015). 195125-19 REVENDA, WOHLFELD, AND CHALOUPKA PHYSICAL REVIEW B 111, 195125 (2025) [64] H. Zhang, K. Haule, and D. Vanderbilt, Effective J=1/2insulating state in Ruddlesden-Popper iridates: An LDA+DMFT study, Phys.Rev.Lett.111, 246402 (2013). [65] B. Lenz, C. Martins, and S. Biermann, Spectral functions of Sr2IrO4: Theory versus experiment, J. Phys.: Condens. Matter 31, 293001 (2019). [66] R. Comin, G. Levy, B. Ludbrook, Z.-H. Zhu, C. N. Veenstra, J. A. Rosen, Y. Singh, P. Gegenwart, D. Stricker, J. N. Hancock, D. van der Marel, I. S. Elfimov, and A. Damascelli, Na2IrO3as a novel relativistic Mott insulator with a 340-meV gap, Phys. Rev. Lett. 109, 266406 (2012). [67] X. Zhou, H. Li, J. A. Waugh, S. Parham, H.-S. Kim, J. A. Sears, A. Gomes, H.-Y. Kee, Y.-J. Kim, and D. S. Dessau, Angle-resolved photoemission study of the Kitaev candidate α-RuCl3,Phys.Rev.B94, 161106(R) (2016). [68] D. Sutter, C. G. Fatuzzo, S. Moser, M. Kim, R. Fittipaldi, A. Vecchione, V. Granata, Y. Sassa, F. Cossalter, G. Gatti, M. Grioni, H. M. Rønnow, N. C. Plumb, C. E. Matt, M. Shi, M. Hoesch, T. K. Kim, T.-R. Chang, H.-T. Jeng, C. Jozwiak et al., Hallmarks of Hunds coupling in the Mott insulator Ca2RuO4, Nat. Commun. 8, 15176 (2017). [69] F. Petocchi, V. Christiansson, and P. Werner, Fully ab initio electronic structure of Ca2RuO4,Phys.Rev.B104, 195146 (2021). [70] A. Kłosi´ nski, D. V. Efremov, J. van den Brink, and K. Wohlfeld, Photoemission spectrum of Ca2RuO4: Spin polaron physics in an S=1 antiferromagnet with anisotropies, Phys.Rev.B101, 035115 (2020). [71] J. Chaloupka and G. Khaliullin, Doping-induced ferromagnetism and possible triplet pairing in d4Mott insulators, Phys. Rev. Lett. 116, 017203 (2016). [72] J. Chaloupka, Emergent transverse-field Ising model in d4spinorbit Mott insulators, Phys.Rev.B109, L020403 (2024). [73] N. Papanicolaou, Pseudospin approach for planar ferromagnets, Nucl. Phys. B 240, 281 (1984). [74] N. Papanicolaou, Unusual phases in quantum spin-1 systems, Nucl. Phys. B 305, 367 (1988). [75] A. V. Chubukov, Fluctuations in spin nematics, J. Phys.: Condens. Matter 2, 1593 (1990). [76] T. Sommer, M. Vojta, and K. W. Becker, Magnetic properties and spin waves of bilayer magnets in a uniform field, Eur. Phys. J. B 23, 329 (2001). [77] S. Sachdev and R. N. Bhatt, Bond-operator representation of quantum spins: Mean-field theory of frustrated quantum Heisenberg antiferromagnets, Phys. Rev. B 41, 9323 (1990). [78] A. V. Chubukov and D. K. Morr, Phase transition, longitudinal spin fluctuations, and scaling in a two-layer antiferromagnet, Phys. Rev. B 52, 3521 (1995). [79] Z. Zhang, K. Wierschem, I. Yap, Y. Kato, C. D. Batista, and P. Sengupta, Phase diagram and magnetic excitations of anisotropic spin-one magnets, Phys. Rev. B 87, 174405 (2013). [80] D. G. Joshi and M. Vojta, Nonlinear bond-operator theory and 1/dexpansion for coupled-dimer magnets. II. Antiferromagnetic phase and quantum phase transition, Phys. Rev. B 91, 094405 (2015). [81] M. Vojta and K. W. Becker, Doped bilayer antiferromagnets: Hole dynamics on both sides of a magnetic ordering transition, Phys. Rev. B 60, 15201 (1999). [82] L. Rademaker, K. Wu, and J. Zaanen, Dynamics of a single exciton in strongly correlated bilayers, New J. Phys. 14, 083040 (2012). [83] L. Rademaker, K. Wu, H. Hilgenkamp, and J. Zaanen, The dynamical frustration of interlayer excitons delocalizing in bilayer quantum antiferromagnets, Europhys. Lett. 97, 27004 (2012). [84] J. H. Nyhegn, K. K. Nielsen, and G. M. Bruun, Equilibrium and nonequilibrium dynamics of a hole in a bilayer antiferromagnet, Phys. Rev. B 106, 155160 (2022). [85] J. H. Nyhegn, G. M. Bruun, and K. K. Nielsen, Wave function and spatial structure of polarons in an antiferromagnetic bilayer, Phys. Rev. B 108, 075141 (2023). [86] K. Wohlfeld, M. Daghofer, A. M. Ole´ s, and P. Horsch, Spectral properties of orbital polarons in Mott insulators, Phys. Rev. B 78, 214423 (2008). [87] C. Brünger and F. F. Assaad, Single-hole dynamics in the Kondo necklace and bilayer Heisenberg models on a square lattice, Phys. Rev. B 74, 205107 (2006). [88] J. Oitmaa and C. J. Hamer, Ground-state properties and oneparticle spectra for a spin-1 Heisenberg antiferromagnet from series expansions, Phys. Rev. B 77, 224435 (2008). [89] J. Revenda, K. Wohlfeld, and J. Chaloupka, Zenodo (2025), https://doi.org/10.5281/zenodo.15125675. 195125-20