Impurity induced magnetic ordering in Sr$_2$RuO$_4$

Ti substituting Ru in Sr$_2$RuO$_4$ in small concentrations induces incommensurate spin density wave order with a wave vector $\boldsymbol{Q} \simeq (2 \pi /3, 2 \pi /3)$ corresponding to the nesting vector of two out of three Fermi surface sheets. We consider a microscopic model for these two bands and analyze the correlation effects leading to magnetic order through non-magnetic Ti-doping. For this purpose we use a position dependent mean field approximation for the microscopic model and a phenomenological Ginzburg-Landau approach, which both deliver consistent results and allow us to examine the inhomogeneous magnetic order. Spin-orbit coupling additionally leads to spin currents around each impurity, which in combination with the magnetic polarization produce a charge current pattern. This is also discussed within a gauge field theory in both charge and spin channel. This spin-orbit coupling effect causes an interesting modification of the magnetic structure, if currents run through the system. Our findings allow a more detailed analysis of the experimental data for Sr$_{2}$Ru$_{1-x}$Ti$_{x}$O$_{4}$. In particular, we find that the available measurements are consistent with our theoretical predictions.


I. INTRODUCTION
Motivated by the discovery of unconventional superconductivity, Sr 2 RuO 4 has been studied extensively for well over two decades due to its many intriguing properties [1][2][3][4]. However, the nature of the superconducting phase as well as its underlying pairing mechanism are still under debate, in particular in view of recent experiments putting the chiral p-wave phase as the previously most prominent candidate for the pairing symmetry into question [5][6][7]. Due to its strongly correlated Fermi-liquid properties in the normal state Sr 2 RuO 4 with its singlelayer perovskite structure has been considered as a twodimensional analog of 3 He [8,9]. Moreover, it belongs to a Ruddelson-Popper series whose end member is the ferromagnetic compound SrRuO 3 . Consequently, the pairing was anticipated to be mediated by ferromagnetic fluctuations. Interestingly, ferromagnetism is not the dominant part of the spin correlations in Sr 2 RuO 4 . Instead, incommensurate (IC) correlations have been experimentally observed at the wave vector Q (2π/3, 2π/3) (lattice constant taken as unit of length) [10][11][12][13][14][15][16], very consistent with theoretical predictions [17][18][19][20][21].
The IC wave vector is attributed to the nesting of two out of three Fermi surface (FS) sheets, which are dominated by the 4d-t 2g orbitals of the Ru 4+ ions. The d xy orbital yields the two-dimensional γ band, while the d xz and d yz orbitals give rise to the α and β bands, which incorporate quasi-one-dimensional features [22,23]. Given their one-dimensional character favoring Fermi surface nesting, these latter two are responsible for the IC correlation [20,24]. Sr 2 RuO 4 falls barely short of establishing magnetic long-range order originating from the strong IC correlations. It was, however, found that replacing some Ru 4+ by non-magnetic Ti 4+ (3d 0 ) induces magnetic order already at rather low Ti concentrations, x > 0.025 in Sr 2 Ru 1−x Ti x O 4 [25]. Inelastic neutron scattering reveals a transition to a phase corresponding to a two- dimensional IC spin density wave state, whose moments are oriented along the c-axis and are modulated corresponding to wave vectors related to the Fermi surface nesting [26,27]. Naturally, Ti-doping suppresses the unconventional superconducting phase rather quickly even before magnetic order is established. Therefore there is no direct interference between the two ordered phases, although understanding magnetic correlations in Sr 2 Ru 1−x Ti x O 4 could shed light on the magnetic fluctuations of the pure compound. Interestingly, non-magnetic impurities have been found to cause magnetic order also in insulating spin systems with antiferromagnetic short-range correlations and a spin excitation gap. In particular, a series of cuprate systems show this feature upon replacing some Cu-ions by Zn impurities [28,29]. This type of systems has been investigated theoretically quite extensively in the context of the concept order by disorder and the modification of magnetic quantum phase transitions by impurities [30][31][32]. The nucleation of a pattern of magnetic order around a non-magnetic impurity observed in these model systems is rather similar to our metallic case. However, there are also obvious differences concerning the coupling of several impurities and the fact that our system is metallic. A recent, prominent case of a metallic system close to a magnetic quantum phase transition is FeSe, which shows antiferromagnetic order in films epitaxially grown on SrTiO 3 , while bulk FeSe does not exhibit magnetic order [33]. It has been suggested that non-magnetic defects and impurities could be responsible for this behavior, whereby in this case superconductivity could play a role [34][35][36].
The aim of our study is to explore the formation of magnetic order through non-magnetic impurities from a microscopic point of view. Since the magnetic state is clearly governed by the Fermi surface nesting of the α and β sheets, we will neglect the third (γ) sheet for simplicity in the following discussion. From the microscopic model of the α-β-bands including spin-orbit coupling we develop a mean field theory for a spatially dependent magnetization of the two involved 4d-orbitals and determine the pattern of magnetic order around a single impurity. We find that the dominant magnetic correlation in the pure system is connected with the four IC wave vectors Q (±2π/3, 0) and (0, ±2π/3), which, if strong enough, would lead to a spin density wave state based on these wave vectors. In this study, our focus lies on establishing the qualitative structure of the magnetization pattern induced by impurities, meaning we do not attempt to give a quantitative discussion.
A simple probe of the magnetic response in the presence of an impurity can be obtained by considering the spin susceptibility χ zz (r, r ) in real space on a square lattice using the electron structure of the α-β-bands without interactions. We observe that the dominant spin correlation yields a cross like pattern with the impurity at the center, as displayed in Fig. 1. Each beam of the cross can be attributed to one of the two orbitals, 4d yz and 4d zx , and the spin polarization is modulated by the corresponding nesting vectors Q (x) (2π/3, 0) and Q (y) (0, 2π/3) in x-and y-direction, respectively. Moreover, the sign is reversed on the two beams such that no net magnetic moment is associated with the impurity. Indeed, we will show below that this is the basic structure, which appears spontaneously due to the interaction around the impurities.
In the following, we will analyze the magnetic structure induced by non-magnetic impurities based on a microscopic model of the α-β-bands including the effects of spin-orbit coupling [Sec. II]. A first attempt using a theory based on averaged impurity scattering does not indicate any inclination of the system towards magnetic order. Thus, we turn towards a mean field approach with spatial resolution, which allows us to show that a nonvanishing spin polarization appears around the impurity due to the repulsive onsite potential. In order to further analyze our findings we consider also the multi impurity situation, which we examine using a phenomenological Ginzburg-Landau formulation [Sec. III]. The Ginzburg-Landau theory is subsequently supported by a field theory derived from the microscopic model. The different approaches fall very much in line and provide a clear picture of the structure of the impurity induced magnetic order. Extending the field theory into a gauge field theory of the charge and spin degrees of freedom reveals several more features of the system, such as (impurity) potential induced spin currents and the connection between spin polarization and charge currents [Sec. IV]. These phenomena may allow for interesting probes of the magnetic properties of such a system.

II. MAGNETISM OF THE TWO-ORBITAL SYSTEM
After introducing the basic two-band model we analyze the magnetic properties first using a formalism based on configuration-averaged impurity scattering and then by focusing on a local mean field theory for a single impurity.

A. The model Hamiltonian
Our model Hamiltonian describes the α-and β-band of Sr 2 RuO 4 , which are derived from the 4d-t 2g -orbitals d xz and d yz , labelled as 'x' and 'y' in the following. It includes nearest-neighbor (NN) intra-orbital and next-nearest neighbor (NNN) inter-orbital hopping corresponding to the pattern depicted in Fig. 2, onsite spinorbit coupling and Coulomb interactions. We split the Hamiltonian into four parts, where H b describes the tight-binding bands of the two orbitals, H SO corresponds to the spin-orbit coupling, H U to the onsite Coulomb interaction between electrons and H imp is the impurity potential. These terms are given by where a † rαs (a rαs ) is the creation (annihilation) operator for electrons on site r, in orbital α = {x, y} and with spin s = {1 = ↑, −1 = ↓} and n r,α,s = a † r,α,s a r,α,s is the density operator. The vectorsx andŷ are primitive lattice vectors of the square lattice in x-and y-direction, respectively (see Fig. 2).We use the notationᾱ = {y, x} and restrict the lattice to 0 ≤ i, j < N , where N is the linear system size. The chemical potential and the spinorbit coupling constant are denoted by µ and λ, respectively. Impurities at the positions R j exhibit an onsite repulsive potential V > 0, which is assumed to be much larger than the band width. The interaction term includes intra-orbital U , inter-orbital K and Hund's rule coupling J. For the last one we neglected the pair hopping, for simplicity, without changing our results qualitatively. We also assume the standard relation between interaction strengths, U = K + 2J.
We can rewrite the Hamiltonian in k space using the spinor notation C k,s = (c kxs , c kys ) T , where Thus, the single-particle part combining H b , H SO and H imp yields with k ± = k ± q/2 and where ξ α,k = −2t cos k α − µ and γ k = 4t sin k x sin k y and V −q = Rj V e iq·Rj (s = ±1). The interaction part can Illustration of all possible square lattice hopping terms for the x and y orbitals. The intra-orbital (inter-orbital) hopping matrix element is denoted by t (t ). be written as and The interaction matrices in orbital space are given bŷ Ignoring impurities we obtain the two-dimensional band structure shown in Fig. 3 from H 0 , where the first Brillouin zone (BZ) is given by k x , k y ∈ [−π, π] with the lattice constant a = 1. In the following we take the NN hopping matrix element t as the unit of energy and choose t = 0.1t and λ = 0.1t. The chemical potential is fixed to keep the electron density at n = 8/3. Since the band structure is dominated by the NN intra-orbital hopping, the band structure retains a strong one-dimensional character leading to pronounced nesting wave vectors Q (2π/3, 0) and (0, 2π/3).

B. Self-consistent T-matrix approximation
We first discuss the effect of Ti-impurities by considering configuration-averaged impurity scattering based on the self-consistent T-matrix approximation. This scheme uses the Green's function formalism and allows us to calculate straightforwardly the spin susceptibility, where S z α,q = (1/N ) k,s,s c † k+qαs σ z ss c kαs with c † kαs and c kαs denoting the creation and annihilation operators, respectively, for electrons of momentum k in orbital α with spin s. We restrict ourselves to the z-axis spin polarization because spin-orbit coupling yields an anisotropy, which boosts the z axis component compared to the inplane polarizability for q ≈ Q, which is most relevant for the magnetic correlations in our discussion [20]. The static spin susceptibility can be expressed by means of Green's functions, G αα with the Fermionic (Bosonic) Matsubara frequencies ω n = (2n + 1)π/β (ν l = 2lπ/β). We define the susceptibility of the orbital α as χ α zz (q) = α χ αα zz (q). For the impurity-free system we use the bare Green's function given bŷ which is derived straightforwardly from the equation, withĤ k,s of Eq. (8). Because the bare Green's function is diagonal in the spin index, we use the shorter notation G 0,s (k, iω n ) with a single spin index.
In line with the impurity potential introduced in the Hamiltonian [Eq. (4)] we write the intra-orbital singleimpurity scattering term as with V k,k = V , which corresponds to s-wave scattering. The impurity scattering leads to the renormalized Green's functionĜ ss (k, iω n ) defined through the Dyson equation, (19) whereΣ(k, iω n ) denotes the self-energy. Note that the Green's function is diagonal in the spin indices.
To account for a large impurity potential we employ the T-matrix approach, which includes multiple scatterings at the impurity [37][38][39]. Due to the s-wave character of the scattering the momentum dependence drops out of the self-energy, which is directly connected with the Tmatrix,T (iω n ), througĥ This relation is only valid for small impurity concentrations n imp , where interference effects between scattering events at several impurities can be neglected.
Inserting it in Eq. (19) we determineT (iω n ) selfconsistently, In Fig. 4 the static spin susceptibility χ (x+y) zz (q), summed over both orbitals, is depicted for the clean (a) and the impurity-doped (b) system. In both cases we observe pronounced peaks at the four wave vectors Q = (±1, ±1)Q 0 with Q 0 ≈ 2π/3, which is not shifted noticeably by disorder. This is obviously a feature originating from the Fermi surface nesting as mentioned above. Impurity scattering introduces, however, a decrease in quasiparticle life time, which broadens the peaks while simultaneously reducing their height. This behavior is systematic as can be seen in Fig. 4 (c). Within this approach we cannot explain the observed magnetic ordering due to Ti-doping. On the contrary, impurity scattering yields apparently a weakening of magnetic correlations. This is a shortcoming of the averaging procedure we use here, which leaves the system formally translationally invariant. The magnetic order we will describe below, however, is not at all homogeneous, but strongly modulated in the system. Thus, a different approach is necessary to analyze the magnetic instability. for (a) no impurities and (b) a finite concentration nimp = 0.05 of impurities in the non-interacting system (U = K = J = 0). The susceptibilities are normalized with the Pauli susceptibility, χP . In (c) the peak maximum of the susceptibility is plotted against the impurity concentration.

C. Local mean field approximation
In contrast to the spatial average of the previous section, we analyze now the magnetic properties around a single impurity using a local mean field approach. For this purpose we decouple the interaction term H U in our Hamiltonian [Eq. (5)] by means of a Hartree-Fock mean field approximation [40], where n rαs denotes the expectation value for the ground state. The resulting single-particle Hamiltonian is then solved self-consistently for the spatially resolved expectation values n rαs , which yield the local spin polarization, As mentioned above, the preference for z-axis orientation of the magnetic moments is caused by spin-orbit coupling.
For our numerical calculation we choose the finite size system of a N × N -lattice with N = 21 and a single im-purity in the center. The use of periodic boundary conditions yields effectively a regular lattice of impurities, which are N lattice constants apart. The fact that N can be divided by 3 makes the cells commensurate with the nesting vector Q 0 = 2π/3. In order to reduce other finite-size effects we average over twisted boundary conditions in the evaluation of the self-consistent equation, introducing phase factors to the hopping matrix elements through the boundaries.
Initially we use the following model parameters for the band structure: t = λ = 0.1t, and for the Coulomb interaction U = 1.5t with J = 0 (keeping U = K + 2J). The impurity potential V is several times the band width. The basic result displayed in Fig. 5 (a) shows a crossshaped pattern of the mean field spin polarization with the impurity position at the center. The narrow beam along the x-axis (y-axis) is strongly dominated by the x-(y-)orbital. The pattern has a basic oscillatory modulation with the wave vector Q 0 in both directions and decays rather rapidly. The spin polarization has opposite signs for the two beam directions. Hence, its symmetry follows the irreducible representation B 1 of the point group C 4v of the square lattice with the basis function Φ B1 (r) = x 2 − y 2 . This means that there is no net magnetic moment associated with the impurity. This basic finding is entirely in line with the structure obtained from the analysis of the local spin susceptibility around an impurity, discussed in the introduction [ Fig. 1].
To probe the influence of the different ingredients of our model on the structure, we vary some of the model parameters. The inter-orbital hopping t and the spinorbit coupling λ are essential for the opposite sign of the magnetization in x-and y-direction (B 1 symmetry). Neglecting the two, t = λ = 0, leads to the degeneracy with the structure where both beams have the same sign (belonging to the irreducible presentation A 1 ). Increasing the Coulomb repulsion U strengthens the spin polarization keeping the same modulation [ Fig. 5 (b)] (note the color scales given for each figure). Turning on the Hund's coupling raises the spin polarization as well [ Fig. 5 (c)]. On the other hand, increasing the spinorbit coupling reduces the spin polarization [ Fig. 5 (d)], whereby the beams gain in width, a feature which is not easy to observe in our figures, but has been tested in our numerical results. The mechanism lies in the transfer of the spin polarization to the orbital polarization through spin-orbit coupling. The angular momentum operator, L z ∝ i a † rxs a rys − a † rys a rxs , involves both orbitals. Thus, on the beam along the x-axis the spin polarization of the x-orbital is shifted to the y-orbital, which has preference to extend along the y-axis. This shift yields a reduction of the spin polarization directly on the beam axis.

III. PHENOMENOLOGICAL APPROACH FOR A SYSTEM WITH MANY IMPURITIES
The mean field theory of the previous section considered the spin polarization induced by a single impurity. The fact that we used periodic (twisted) boundary conditions for a finite system led effectively to a regular lattice of impurities, which are rather far apart from each other and, thus, weakly coupled. In this section we would like to examine the interplay of impurities at nearby locations and the effect on the magnetization pattern. While this can, in principle, also be dealt with in our previous mean field setup, a phenomenological approach is here more efficient. We formulate a Ginzburg-Landau theory, which we compare with an effective field theory derived from our microscopic model in the subsequent section.

A. Ginzburg-Landau theory
As in the mean field approach we consider the z-axis spin polarization in each orbital (x and y) as a position dependent order parameter, m x (r) and m y (r). These two order parameter components can be combined to belong to two irreducible representations of the square lattice point group C 4h , where m x and m y can be considered to transform like the basis functions x 2 and y 2 , respectively. Following Landau's scheme one could construct the Ginzburg-Landau theory based on the two order parameters m A and m B , describing the spontaneous symmetry breaking in one of the two representations. However, it is actually more advantageous to formulate the free energy expansion with m x and m y . Eventually, multiple impurities will reduce the point group symmetry and lead to a coupling of the order parameters of the two representations. The scalar free energy functional is given by where a(T ) = a (T − T 0 ) and c, b 1,2,3 and K 1,2,3 > 0. Note that we used the combinations m A and m B for the fourth order terms for compactness. The spatial modulation of the order parameter with the wave vector Q 0 is introduced by including higher order derivatives into the gradient terms where κ = 1/2Q 2 with Q ≈ Q 0 = 2π/3. The last term represents the effect of the impurities where the coupling γ < 0 facilitates the local nucleation of the order parameter. We ignore any internal structure of the impurity sites and approximate it by a δ-function for simplicity.
In order to analyze the nucleation of the impurity induced magnetic pattern we restrict ourselves to the linearized Ginzburg-Landau equation obtained through the variation of F[m x , m y ] with respect to the two order parameter components and by neglecting non-linear terms in m x,y , which is written in q-space using the Fourier transformation The symmetric matrix A αα in q-space reads and the right hand side of Eq. (28) uses The linearized GL equation (28) includes also the bulk instability. Setting γ = 0, the bulk instability to an incommensurate spin density wave state occurs at the temperature where the first of the two eigenvalues of the 2 × 2-matrixÂ(q, T ) changes sign for a certain q ( ≈ Q 0 ) from a positive to a negative value at a temperature T = T N < T 0 . In Sr 2 RuO 4 the bulk instability condition would lead to T N < 0. We restrict ourselves to a temperature range T > T N , where for all q the ma-trixÂ(q, T ) is strictly positive definite. We introducê η(q, T ) =Â(q, T ) −1 and solve Eq. (28) formally, which with Eq. (32) leads to an eigenvalue equation, where we define the coupling From App. A we see that η αα (q) is related to the spin susceptibility such that Λ αα jj has the features of an RKKY-type of coupling between impurity sites, although we deal with non-magnetic impurities.
We solve Eq. (34) as an eigenvalue equation, which gives the criterion for the nucleation of magnetic order around the impurities. In particular, it allows us to determine the critical temperature T N (or the critical γ for fixed temperature) together with the structure of the spin polarization pattern through the eigenvector Π αj and use of Eqs. (29) and (33).
Before the discussion of the multi impurity problem in the next subsection, we briefly analyze the single impurity part of the magnetic instability, which leads to the reduced equation where Λ αα jj (T ) = d 2 q η αα (q, T ).
By symmetry we find Λ(T ) = Λ xx jj (T ) = Λ yy jj (T ) and Λ (T ) = Λ xy jj (T ) = Λ yx jj (T ). Taking the single impurity as the symmetry center the point group C 4v is preserved and we can separate the solutions according to the two representations into which yields two decoupled equations, The comparison with the microscopic theory shows that Λ xx jj > 0, Λ xy jj < 0 and Λ xx jj − Λ xy jj > 0. As a result, the solution of the B 1 representation would be the first to nucleate in accordance with the finding of the previous chapter. Returning to the real space magnetization using Eqs. (29) and (33), we find a cross-like pattern, as depicted in Fig. 6, which is qualitatively similar to the one found in Fig. 5.

B. Magnetic order with several impurities
In order to illustrate the interplay between different impurities we consider the case of two impurities in multiple spatial configurations. In view of the solutions of the local mean field calculation in the last chapter, we use parameters yielding a correlation length of a few lattice constants at the nucleation. This means also that the inter-impurity interaction Λ αα jj decays on this length scale. Thus, we consider impurity arrangements in a comparable range of distances as depicted in Fig. 7 (af). In order to be closer to the microscopic description, we use for these simulations not the matrixÂ(q) of the Ginzburg-Landau description, butÂ(q), which we will derive through an effective field theory in the next chapter. This has the advantage that the q-dependence is treated more adequately.
We analyze here the nucleation of magnetic order by solving the corresponding eigenvalue equation, Eq. (34), for the critical temperature, which is slightly different for each configuration. From Π α,j we obtain the real space magnetization pattern shown in Fig. 7.
It is important to note that the symmetries leading to Eqs. (38) and (39) are lost in the presence of more than one impurity. There may be still some residual symmetries depending on the relative location of the impurities. This leads to a coupling of the two representations and a discussion in terms of Π x,j and Π y,j is more accessible.
In Fig. 7 (a-d) we display configurations where the impurities are placed along a main axis (x-axis for y j = 0) at different distances. For longer distances the cross structure of each impurity is like that of a single impurity. As the impurities get closer the overlap of the magnetization pattern increases such that the beam along the x-axis grows in strength and the y-axis beam becomes comparatively weaker (note the color scales of each panel). A further feature is the change of the relative sign of the two impurity patterns, which is a result of the oscillatory RKKY-like interaction structure and connected with the wave vector Q 0 . Note that we do not include the effect of electron depletion due to the impurity. Consequently, the polarization is generally not suppressed at the impurity sites. Only if the polarizations of both orbitals cancel each other, we find m x + m y = 0 at the impurity sites. Figs. 7 (e, f) show configurations where the two impurities are placed along the different directions. This then results in a more complex pattern, where the crossshape however is still visible. The increased polarization between the two impurities goes in line with a slightly enhanced nucleation temperature. The critical temperature T N for a single impurities is, therefore, on a mean field level a lower bound for the nucleation temperature of the many-impurity system, where the transition may be rather inhomogeneous.

IV. EFFECTIVE FIELD THEORY
In this chapter, we introduce a field theory based on the microscopic model to validate our Ginzburg-Landau approach. It will also allow us to consider further aspects of the system that have been omitted so far by including U (1) and SU (2) gauge fields, which give access to charge and spin currents in the system.

A. Formulation and basic results
Our starting point is the partition function Z expressed as a path integral [41], where β = 1/T and the Lagrangian L is derived from the Hamiltonian in Eq. (1), in particular Eqs. (7) and (9), In a first step we perform a standard Hubbard-Stratonovich transformation decoupling the Coulomb interaction with the auxiliary field φ q = (φ x q , φ y q ) T , which are conjugate fields to (m qx , m qy ) [41][42][43]. This leads to with The Hamiltonian matrixĤ k , the impurity potential V q and the interaction matrixÛ are defined in Sect. II A and φ q = diag(φ x q , φ y q ). For simplicity, we omit the charge density fields, ρ q [Eq. (10)], which can be absorbed in the chemical potential, and, thus, also neglect the effect of the repulsive inter-orbital coupling K. We can extend the Lagrangian by introducing the electromagnetic U (1) gauge field A and the spin SU (2) gauge field a, which couple to charge and spin currents, as will be discussed later. The complete Lagrangian is given in App. A in Eq. (A2).
After integrating out the Fermionic degrees of freedom we end up with an effective action (for technical details see App. A). Neglecting the gauge fields for now, we find up to second order in the auxiliary field φ α q , with the bare susceptibility and Note that we use here the short notation q = (q, iν l ).
Minimizing this action we obtain the mean field solution for φ α q , or for m α q , using the relation, Obviously, the renormalized couplingÂ(q) incorporates the random phase approximation for the susceptibility [see also Eq. (B5)]. Note that Eq. (50) corresponds directly to the variational equations obtained within the Ginzburg-Landau theory [Eq. (28)]. However, we should be aware thatÂ(q) is only a qualitative approximation ofÂ(q), since it is a small-q expansion, but attempts to mimic the behavior around Q 0 . After a comparison we find for the coefficients of the gradient terms of the Ginzburg-Landau free energy [Eq. (27)] that K 1 ∼ K 3 K 2 . The essentially vanishing K 2 ensures the clear crossstructure.

B. Correlation functions connecting gauge fields and potentials
The Hamiltonian H b + H SO in Eqs. (2) and (3) incorporates a peculiar feature. The topology of the tightbinding matrix contains a hopping by one lattice constant along the x(y)-direction for the x(y)-orbital. Together with the diagonal inter-orbital hopping in the square plaquette a triangle path is formed. It can be closed in the plaquette corner by spin-orbit coupling connecting the two orbitals. Since spin-orbit coupling introduces a phase sπ/2 depending on the spin s, the triangles host a circulating spin current (see Fig. 8). Combining all four triangles of a plaquette leads to a net circular current on the edges, while it cancels on the diagonals. In the homogeneous system the circular currents of neighboring plaquettes compensate each other. Interestingly, an impurity potential on a site interrupting four closed triangular paths leads to a net circular spin current around the impurity, as illustrated in Fig. 8.
This effect, which is due to the spin-orbit coupling and the specific form of the inter-orbital hybridization, can be incorporated into our effective field theory by means of gauge fields. In Eq. (A6) we have included the U (1) gauge field A and the SU (2) gauge field a, where for the latter we restricted ourselves to the z-axis component of the spin. The expansion to lowest order leads to new terms, which couple the gauge fields to the potentials for the charge and spin density. The extended effective action takes the form with where we introduce the correlation functions andC µµ (q) = 1 β k,s,iωn α,α ,α,α We have introduced here the current matrix in k-space, defined through the relation Setting = 1 we find for the x-and y-component cos k x sin k y and W y k = 4t cos k y sin k x . These current matrices involve both the nearest and next-nearest neighbor contributions and do not depend on the spin s.
Note that in Eq. (54) we have already omitted those terms which vanish entirely after summing over (k, n, s). The finite terms are summarized in Table I. SC   TABLE I. The lowest order coupling terms of S g eff , which appear by including the U (1) gauge field A and the SU (2) gauge field a, and their respective correlation functions (CF). The origin of each term is given in the last column (CD: charge density, SD: spin density, CC: charge current, SC: spin current).
We will work in the following with the density-current part and list the current-current contribution just for completeness. The correlation function C µ α (q) has the form and while C x x (q) = C y y (q) = 0. It is worth noting that these correlation functions include in the denominator the combination of hoppings along the main axes (ξ x/y,k → t) and the plaquette diagonals (γ k → t ), as well as the spinorbit coupling (λ), which constitute exactly the ingredients used above for the circular spin currents in Fig. 8.
A straightforward symmetry analysis of the expressions for C x y (q) and C y x (q) shows that where the function Υ(q x , q y ) is even in q and Υ(−q x , q y ) = Υ(q x , q y ), Υ(0, 0) = 0 and generally Υ(q x , q y ) = Υ(q y , q x ). We will use this property to consider features such as the spin currents around an impurity as anticipated above, but also the connection between spin density polarization and charge currents, which are suggested by our gauge field terms.

Impurity induced spin currents
We first consider the spin current pattern due to a point like impurity potential. In this case the potential is structureless, V q = V , and the spin current is simply obtained as It is actually easier to consider the vorticity of this current rather than the current pattern itself, i.e. Ω s (r) = [∇ × j s (r)] z , It is obvious that q 2 x Υ(q y , q x ) + q 2 y Υ(q x , q y ) has the full symmetry of the point group C 4v . The vorticity Ω s (r) as a function of r has thus the same symmetry, which belongs to the irreducible representation A 1 , as Fig. 9 displays. This means that we find a circular spin current around the impurity site analogous to the schematic picture shown in Fig. 8. However, the correlation function has the property that it peaks for q x ≈ ±Q 0 and q y ≈ ±Q 0 , which yields an oscillation with distance such that the spin current reverses sign several times for increasing r until it fades away. It is also obvious from Eq. (63) that the net (integrated) vorticity vanishes.
Furthermore, it is possible to compare these findings for our square lattice model numerically, neglecting interactions. For this purposed we define the current operators for nearest neighbor bonds as and for next-nearest neighbor bonds analogously as This allows us to obtain charge (spin) currents by adding (subtracting) the bond currents of different spins. Using this formulation we determine the bond spin currents around a single impurity and obtain the pattern shown in Fig. 10. Close to the impurity the current pattern agrees qualitatively well with the picture of Fig. 8, as shown in the enlarged figure (b). However, further away the bond current pattern becomes more difficult to analyze. A better view is obtained by considering the spin current density passing through the sites which is obtained by simply averaging over all bonds adjacent to a site. In this way the circular character of the current is well visible, including its oscillations in size and orientation. Within our numerical calculations we have tested that the sign change of either t (inter-orbital hopping) or the spin orbit coupling λ lead to a reversal of the current pattern.

Charge currents induced by magnetization pattern
Apart from the impurity induced spin currents, we also see that the spin polarization induces charge currents, which are derived from the variation of S g eff with respect to the vector potential A, It is again illustrative to look at the vorticity, We consider now the pattern around a single impurity. For φ α q , we have then the relation where both are even functions of q. In this way, it is obvious that the expression in the integrand has the structure of a basis function of B 1 in C 4v , i.e. like x 2 − y 2 . This is reflected by the vorticity pattern illustrated in Fig. 11, which clearly exhibits the same B 1 symmetry. The charge currents generate the same magnetic field pattern as the spin magnetization.

Spin polarization induced by external charge currents
An interesting effect of the interplay between current and spin polarization is found when we consider an external current running through the system, which induces a spin polarization above the nucleation temperature of the magnetization. To illustrate this feature we use again the expression for the charge currents given in Eq. (66). Instead of looking at the reaction of the system to a charge current, we determine the current density for a specific spin polarization, which extends from the impurity (at the origin) along the y-axis and has opposite signs for positive and negative y. Thus, we assume that which involves strong localization along the x-direction.
Inserting φ y −q and taking into account that the current is source and drain free, i.e. ∇ · j c (r) = 0, requires that also φ x −q is finite. This leads to the current densities The symmetry of the current density corresponds to the pattern of a flow passing an obstacle (impurity) with We may also invert the point of view. This means that a uniform current parallel to a main axis in our system would induce a magnetization pattern around an impurity, which is perpendicular to the current and antisymmetric under reflection at the current axis.
To test this we go back to our mean field discussion of the microscopic model and consider a square shaped system with one impurity at the center. We induce a current by means of a uniform vector potential along the x-direction (Peierl's phase φ on nearest-neighbor hopping matrix elements), which then adapts to the situation around the impurity site. Note that for technical reasons we average over the twisted boundary conditions in order to reduce finite-size effects.
In accordance with our previous discussion, a spin polarization pattern emerges around the impurity, which extends as a narrow beam along the y-direction starting at the impurity and has opposite signs on the two sides of the impurity, as depicted in Fig. 12. Evidently, the impurity acts again as a nucleation center facilitating the appearance of a finite spin polarization. The results shown are for two different currents, which correspond to the Peierl's phases φ = 0.06π and 0.31π. Analyzing the model parameter dependence, we observe that the magnitude of the pattern depends not only on the strength of the current, but also strongly on t and λ such that the pattern vanishes when one or both parameters are zero. This affects also the distribution of the spin polarization on the two orbitals.

V. CONCLUSION
The goal of our study has been to find the origin and structure of the magnetic order induced by non-magnetic impurities in Sr 2 Ru 1−x Ti x O 4 . Surprisingly, already a few percent of Ti-ions replacing Ru-ions lead to the magnetic instability [25][26][27], which can be described well within a relatively simple model of two bands derived from the 4d-t 2g orbitals, d yz and d zx , of the Ru-ions. The pronounced nesting feature of the two bands explains the spatial modulation of the magnetization with wave vectors |Q 0 | ≈ 2π/3. Moreover, spin-orbit coupling is an additional ingredient, which favors the spin polarization to be along the z-axis, as discussed in Ref. [20]. Both features agree very well with experimental findings [26].
Because a simple impurity averaged approach is not capable of capturing the impurity induced magnetic order, we concluded that a more appropriate technique is necessary. Restricting ourselves to the two relevant bands, we could show based on a position-dependent mean field theory that a repulsive impurity potential acts as a nucleation center for the spin polarization. Around a single Ti-impurity the magnetization pattern takes a cross-like structure with beams along the crystalline main axis, each of which is dominated by one of the two involved orbitals. The orientation of the magnetization of the two beam directions is opposite such that no net magnetic moment arises. For these properties the inter-orbital hybridization as well as the spin-orbit coupling are essential. The mean field calculation also suggests that the nucleation of magnetic order occurs at a temperature higher than the transition to the incipient, uniform, incommensurate magnetic order. This would be consistent with the notion of impurity induced order.
The interplay between different impurities is a complex problem. In order to analyze it we introduced a phenomenological approach based on a Ginzburg-Landau free energy functional with a magnetic, two component order parameter, m x and m y , for the z-axis oriented magnetization of the two orbitals, d xz and d yz , respectively. It is rather straightforward to formulate a free energy functional reproducing the single impurity result of the space-dependent mean field theory. Considering two impurities, we find that their relative position and orientation are highly important for the interaction. While the beams along the crystal axis and their overlap play a crucial role, also the matching of the nesting wave vector with the separation of the impurities influences the emerging pattern. The analysis of the interaction incorporates features reminiscent of the RKKY-type coupling among localized magnetic moments in a metal, despite the fact that Ti enters as a nominally non-magnetic impurity. We show that one can map the Ginzburg-Landau model to a model, where each impurity has two magnetic moments along the z-axis, one for each orbital.
The Ginzburg-Landau functional can be derived through a field theory based on the microscopic model, which yields a more adequate spatial structure of the magnetization. It confirms the basic structure and adds details about the phenomenological parameters. In particular, it is obvious that the repulsive potential of the Ti-impurity yields an "attractive" nucleation center for the magnetic moments.
The field theory allows us to uncover a few further aspects of the system, which were not noticed in the mean field theory and the Ginzburg-Landau treatment. After adding the U (1) and SU (2) gauge fields to the field theory we examined the effects of charge and spin currents. The particular topology of orbital hybridization and spinorbit coupling leads to circular spin currents around the impurities, as any potential changing the electron density leads to spin currents. This kind of spin currents has been previously noticed at surfaces [44]. At the same time, we find that a non-vanishing magnetization yields charge currents. We have investigated these for a single impurity, where the magnetization yields a clover-like current pattern. In addition to that, we also find that a uniform charge current could induce a spin polarization around an impurity even above the spontaneous nucleation of the magnetization pattern. While this feature may be too weak to be easily observed, it is interesting to notice that currents could also influence and possibly deform the magnetically ordered pattern.
With our theory we have illustrated that Ti-impurities naturally lead to inhomogeneous magnetic order. The complex structure of this order depends on the random positions of the Ti-ions. Since we are interested here on the possibility of nucleating magnetic order by Ti-doping, we leave the discussion of the random multi impurity system to future studies. Our study also does not intend to give a quantitative account of the phenomenon.
In this section we elaborate on the details of the derivation of the effective gauge field theory used in Sec. IV A, starting with the functional integral form of the partition function of the microscopic model.
We use the auxiliary fields φ ij = (φ x ij , φ y ij ) to decouple the Coulomb interaction by means of a Hubbard-Stratonovich transformation, which results in the partition function including the integral over φ ij , whereĤ k,s [Eq. (8)] denotes the Hamiltonian andÎ k [Eq. (58)] the current matrix. We use here the shorthand notation k ± = k ± q/2 and the inverse interaction matrixÛ as well asφ −q = diag(φ x −q , φ y −q ). The potential of the impurities is given by V −q = j V e iq·Rj , where R j are the impurity positions and V > 0.
After integrating over the Fermionic degrees of freedom we end up with the effective action, which allows us to express the partition function as with M (k−,iωn,s),(k+,iωm,s ) = (−iω n +Ĥ k,s )δ ωn,ωm δ q,0 + A −q ·Î k + sa −q ·Î k The bare Green's function denoted asĜ 0,s is defined through with Fermionic (Bosonic) Matsubara frequencies, ω n = (2n + 1)π/β (ν l = 2lπ/β). We may now expand the ln(M ) term using Restricting ourselves to terms up to second order in the field φ α q , we obtain for the effective action the expression given in Eq. (46), which can now be used to analyze the magnetic stability of the system through variation with respect to φ α q in the static case, 0 = ∂S This corresponds to a linearized version of a Ginzburg-Landau theory. So far we have written the effective action in terms of the auxiliary fields φ α q . Alternatively, we can also express S We may now use these expressions to estimate the parameters of the Ginzburg-Landau free energy functional given in Sec. III A. The symmetric matrixÂ(q) can be compared withÂ(q). The impurity term involves some spatial structure not present in the Ginzburg-Landau theory, where the impurity term is simply approximated by a delta-function, Most notably, we find that for positive V the average of the left side is negative, which implies that γ < 0. The analysis of the instability by means of the field theory yields the result displayed in Fig. 6. We observe qualitatively the same cross like structure as in the local mean field theory of Sect. II C.
which included in the effective action [Eq. (46)] leads to Ignoring the impurities, the saddle point equation in the paramagnetic phase for the static field φ α q is and the magnetic moment m α −q can be defined through Now combining Eqs. (B3) and (B4) we obtain for the magnetization which defines the renormalized susceptibility on the level of the random phase approximation. On the other hand, we can derive from the two Eqs. (B3) and (B4) the relation between m α q and φ α q , which we used in the previous section.