Unveiling a nematic quantum critical point in multi-orbital systems

Electronic nematicity, proposed to exist in a number of transition metal materials, can have different microscopic origins. In particular, the anisotropic resistivity and meta-magnetic jumps observed in Sr3Ru2O7 are consistent with an earlier proposal that the isotropic-nematic transition is generically first order and accompanied by meta-magnetism when tuned by a magnetic field. However, additional striking experimental features such as a non-Fermi liquid resistivity and critical thermodynamic behavior imply the presence of an unidentified quantum critical point (QCP). Here we show that orbital degrees of freedom play an essential role in revealing a nematic QCP, even though it is overshadowed by a nearby meta-nematic transition at low temperature. We further present a finite temperature phase diagram including the entropy landscape and discuss our findings in light of the phenomena observed in Sr3Ru2O7.


Introduction
A variety of transition metal materials such as the cuprates [1], Ru-oxides [2], and Fepnictides [3] has been proposed to harbour an electronic nematic phase [4]. Electronic nematic phases are broadly characterized by the presence of spontaneously broken rotational symmetry and viewed as the quantum counterpart of nematic classical liquid crystal phases. The theoretical proposal of nematic quantum liquid crystals became more concrete when experiments on ultra-pure bilayer ruthenate (Sr 3 Ru 2 0 7 ) samples subjected to a magnetic field along the c-axis revealed an unusual phase characterized by a pronounced residual resistivity in place of a putative meta-magnetic quantum critical point (QCP) [5]. Interestingly, Sr 3 Ru 2 O 7 was initially viewed as a prototype for the study of quantum phase transitions, exhibiting a striking non-Fermi liquid resistivity thought to originate from the putative magnetic field tuned QCP [6]. The unusual phase found in ultra-pure samples is delimited by two consecutive first order meta-magnetic transitions at low temperature and, remarkably, exhibits a significant in-plane magnetoresistive anisotropy when the external field is slightly tilted towards one of the in-plane crystal axes [2]. These observations strongly imply the formation of an anisotropic metallic, i.e. electronic nematic, phase in the bilayer ruthenate compound.
Based at first on the two consecutive metamagnetic transitions, an electronic nematic phase was proposed and generic features of nematic phase formation were theoretically explored early on [5,7,8]. It was found that the transition between the isotropic and nematic phase is generally first order, and that nematic order typically develops near a van Hove singularity (vHS) to avoid a Lifshitz transition. Varying the chemical potential, the nematic phase is bounded at low and high values by two isotropic phases, while the concomitant first order transitions lead to jumps in the electron density. When a magnetic field is applied (and the chemical potential is held fixed at, say, some low value), Zeeman coupling acts as a spin-dependent chemical potential term. Tuning the Zeeman field increases the volume of, e.g., the spin-up species which then enters the nematic phase near the vHS, resulting in a jump in the spin-up density at the isotropic-nematic transition. In contrast, the spin-down Fermi volume continuously decreases without passing a van Hove point. Such behavior gives rise to meta-magnetism, i.e. a sudden jump in the difference of up-to down-spin density at the isotropic-nematic transition [9]. Using more realistic band structures, microscopic models were suggested and the above conclusion was found to be robust [10,11].
Despite this mechanism of nematic phase formation, the origin of the critical signatures remains mysterious and recent thermodynamic data revealed additional complexity that requires new insight into the nematic theory [12]. For over a decade in temperature T , the resistivity behaves nearly perfectly linear [2], while the specific heat coefficient (C/T ) varies as log(T ) over the same range [12,13]. This behaviour resembles quantum criticality observed in a variety of other correlated materials, suggesting the existence of an unidentified QCP related to the nematic phase. How can one reconcile the first order nature of the isotropic-nematic transition and the critical behavior associated with a second order transition?
In this paper, we present a way to resolve this puzzle by introducing two nematic order parameters relevant for multi-orbital systems, and investigate as to how their interplay leads to interesting physics. While the current work is motivated by the bilayer ruthenate, our study on the isotropic-nematic transition is generally applicable 3.

0.
1. to a t 2g -orbital system. We show that in the presence of moderate spin-orbit (SO) coupling the multi-orbital nature of the system is a key property for turning a nematic first order transition into a QCP. We also show that the nematic QCP is accompanied by a nearby jump in nematicity, dubbed here meta-nematic transition, which obscures the QCP at low temperature. This finding is in contrast to the current wisdom that a QCP is hidden under the nematic dome. Figure 1 (a) shows the finite temperature phase diagram of such a multi-orbital system, including the contour plot of the entropy difference between the nematic and the isotropic phase. The nematic QCP marked by × is found on the right low temperature phase boundary. Note also that the entropy inside the nematic phase near the meta-nematic transition indicated by the arrow is higher than in the isotropic phase. While the separation between the continuous and the meta-nematic transition depends on the underlying band structure and the SO coupling strength, thermodynamics and transport at finite temperature will be governed by the QCP aside from singular density of states (DOS) effects, which may be relevant for the phenomena observed in Sr 3 Ru 2 O 7 [2,5,6,12,14,15,16]. A more detailed discussion will be presented below.

Nematic QCP
Theoretical progress has recently been made in developing microscopic models for nematic phase formation in Sr 3 Ru 2 O 7 . One approach suggested that the nematic phase is a spontaneous Fermi surface (FS) distortion arising from a band with strong d xy character and a vHS near the Fermi level [11], while another proposed that it originates from orbital ordering, i.e., the density difference between d xz and d yz orbitals, ignoring the d xy orbital [17,18]. Both scenarios lead to broken x-y symmetry. To distinguish both types of ordering in the present paper, we call the former nematic and the latter orbital order. Although the underlying microscopic driving mechanisms are distinct, both approaches have noted the importance of SO coupling, which is in line with the results of angle resolved photoemission spectroscopy (ARPES) [19]. However, they failed to locate a QCP. Thus, given that all t 2g orbitals are coupled via SO interaction, we explore here a broader notion of x-y symmetry breaking in a multi-orbital system.
To understand the origin of nematicity and its consequences, we start from a tightbinding model, which reproduces the FS of Sr 3 Ru 2 O 7 [20,21,22] and was introduced in [11]. The model is based on the Ru 4d t 2g orbitals and includes moderate SO coupling H SO = 2λ i L i · S i . We also incorporate a staggered lattice potential to account for the effect of the staggered rotation of the RuO 6 octahedra, which doubles the real space unit cell [23,11,24]. In principle, the lattice potential can be momentum dependent (e.g. g(k) ∝ cosk x + cosk y as in [24]) but here we assume a constant valueḡ as the results do not depend on its momentum structure. The tight-binding band structure then has the form where consists of electron operators creating an electron with spin σ =↑, ↓ in one of the t 2g -derived orbitals α = yz, xz, xy. The orbital dispersions are given by ε yz/xz k = −2t 1 cosk y/x − 2t 2 cosk x/y + µ, ε xy k = −2t 3 (cosk x + cosk y ) − 4t 4 cosk x cosk y − 2t 5 (cos(2k x ) + cos(2k y )) + µ, while the hybridization between the quasi-1D orbitals is ε 1D k = −4t 6 sink x sink y . In the following the underlying band structure parameters are t 1 = 0.5, t 2 = 0.05, t 3 = 0.5, t 4 = 0.1, t 5 = −0.03, t 6 = 0.05,ḡ = 0.1 and λ = 0.14 (if not stated otherwise), while the chemical potential µ serves as a tuning parameter. All energies throughout the paper are expressed in units of the intra-orbital nearest neighbour hopping integral 2t 1 = 1.
Since the nematic phase reported in Sr 3 Ru 2 O 7 arises from a metallic state, it is reasonable to assume that long-range interactions are well screened, leaving moderately weak (compared with the bandwidth) on-site and nearest neighbour interactions. The microscopic interaction Hamiltonian is then with n α iσ = c α † iσ c α iσ the density operator for electrons in orbital α at site i and spin σ. Here, U ,Ũ , and V α represent repulsive intra-orbital on-site, inter-orbital on-site, and intra-orbital nearest neighbour interactions, respectively. We assume that V xz(yz) is finite along nearest neighbour x(y)-bonds, while V xy is finite along all four nearest neighbour bonds, and that V α ≡ V . Different instabilities compete within H 0 +H int . Most natural candidates are spin density wave and ferromagnetic instabilities. However, it was shown in [17] that orbital ordering between d yz and d xz orbitals dominates when the inter-orbital interaction is significant between the two quasi-1D orbitals in the presence of SO coupling. This important observation, however, did not take into account the 2D d xy orbital, which dominates the γ FS sheets of Sr 3 Ru 2 O 7 as observed by ARPES [20]. In particular, the γ 2 sheet possesses a vHS near the Fermi level implying that it is most susceptible to an instability. Reference [11] indeed suggests that nematic ordering in the xy dominated bands is the leading instability when the nearest neighbour interaction is taken into account. Since these theories imply that ferromagnetic and antiferromagnetic spin density wave instabilities are suppressed by SO coupling, we focus here on the interplay of nematic and orbital ordering.
Introducing nematic and orbital order parameters, ∆ σ N = N −1 k (cos k x − cos k y ) n xy kσ and ∆ O = N −1 k,σ n xz kσ −n yz kσ , respectively, one arrives at the following Hamiltonian  . Each panel represents the nematic order parameter for various SO coupling strengths λ (dotted lines visualize discontinuities). Note that the first order transition at the right phase boundary becomes continuous for 0.07 < λ < 0.2. The results in figure 1 and 2 are based on λ = 0.14. Very recent unpublished transport data [25] in fact bear a strong qualitative resemblance to this λ = 0.14 plot. The resistivity has a pronounced shoulder for fields higher than those of the two first-order transitions, and within this shoulder becomes anisotropic under the application of small in-plane magnetic fields.
where the effective interactions are given by Note that the on-site intra-orbital interaction U hinders orbital ordering (favouring magnetic ordering instead), but does not interfere with nematic ordering. On the other hand, the inter-orbital interactionŨ favours orbital ordering, while the nearest neighbour intra-orbital interaction V favours nematic over orbital order. The self-consistent solutions for nematic (∆ N = ∆ ↑ N = ∆ ↓ N ) and orbital order over the T −µ plane are shown in panels (I) and (II) of figure 2. Panel (III) shows the FS at selected points for zero temperature. As pointed out in figure 1, nematic and orbital order develop continuously from the isotropic phase at µ c3 = −0.47, in contrast to the generic first order transitions encountered when either nematic [7,8,11] or orbital [17] ordering is considered.
For the present discussion we have set V N = 0.8, while V O (≪ V N ) is chosen to be small or zero such that orbital order can only be induced by coupling to nematic order via SO interaction. However, orbital ordering can also be induced for larger V O and smaller V N values than the above choice (for example for V O = 0.2 and V N =0.6) but does not alter the main conclusion of the existence of the QCP and the nearby meta-nematic jump. The reason for preferring the nematic to be the leading instability (V N > V O ) is based on experimental constraints. We found that an independent orbital instability occurs for V O 0.55 (when V N = 0) but that it is also accompanied by a drastic change in the overall FS topology, exhibiting various strongly 1D-like open FS sheets inconsistent with the de Haas van Alphen experimental findings for Sr 3 Ru 2 O 7 . In contrast, the nematic instability gives rise to smooth changes of the FS before and after the nematic window as shown in figure 2 (III), in agreement with the de Haas van Alphen results [21,22]. Hence we consider the nematic instability as the relevant instability in this study. ‡ To understand the effect of SO coupling, we investigate how the phase transitions are modified by changes in the SO coupling strength λ. Figure 3 shows the nematic order parameter for different λ (orbital order is not shown, but is small and disappears for λ = 0). The nematic order parameter exhibits typical first order transitions for small λ. For intermediate λ, the first order transition at the higher critical chemical potential turns into a continuous transition preceded by a sudden meta-nematic jump. The change in nature from first order to second order is driven by the coupling between the two order parameters enabled by SO interaction and can be captured qualitatively by the following Ginzburg-Landau (GL) free energy analysis.

GL free energy analysis
A generalized GL free energy takes the form where we assume that the higher order terms denoted by the ellipsis are negligible. The GL free energy analysis is a phenomenological approach, and it is important to understand how each term is allowed by symmetry, in particular the δ term. The nematic order parameter breaks only the x-y symmetry of the underlying lattice (or the rotational symmetry in the case of the Jellium model). Thus the GL free energy should be invariant under a sign change of the order parameter (e.g. ∆ N → −∆ N , corresponding to a π/2 rotation), which restricts the GL expansion to only even powers in the nematic order parameter. Using the Jellium model of the electron gas, the GL free energy with nematic order was analyzed in [26], where the isotropic-nematic transition is assumed to be continuous with α n changing sign at the transition, while β n stays positive. However, it was first pointed out in [8] using a single nematic order parameter that the free energy has a log-singularity due to a vHS (Lifshitz transition) leading to a first order isotropic-nematic transition in a simple square lattice model. An order parameter expansion shows that the β n coefficient becomes negative while α n and γ n remain positive at the transition point [8]. Since then, turning the first order transition into second order has become an interesting problem. In the present study of multi-orbital systems such as the t 2g -orbital Sr 3 Ru 2 O 7 compound, there are two x-y symmetry breaking order parameters as introduced in the previous section. Both orbital ordering (i.e. the density imbalance between d xz and d yz orbitals) and nematic ordering (defined on the d xy -orbital manifold) break the same x-y symmetry. In principle, the GL free energy contains a finite coupling term such as ∼ ∆ N · ∆ O , as long as the Hamiltonian contains mixing between d xy and quasi-1D (d xz and d yz ) orbitals. However, due to the nature of the orbitals and the inversion symmetry of the square lattice about the x-y plane, finite hopping integrals between d xy and the quasi-1D orbitals are absent. Therefore, the dominant mixing between d xy and quasi-1D orbitals is SO coupling. Under the symmetry considerations discussed above, the δ-term is allowed when SO interaction is present, i.e. δ ∝ λ + O(λ 3 ). Starting with positive coefficients only, no ordering is initially present. However, as was shown previously for a single band, ‡ While bilayer coupling may alter the role of orbital and nematic orders, the existence of a nematic QCP is not sensitive to the choice of leading and induced orders. the nematic instability is typically of first order due to β n < 0 and α n , γ n > 0 [8]. Therefore, without SO coupling (δ = 0) and α o , β o > 0 finite nematic order develops when β 2 n > 4α n γ n , while orbital ordering is absent (∆ O = 0). On the other hand, when δ = 0 nematic order starts to develop continuously for δ 2 > 4α n α o , inducing also orbital ordering (with ∆ O = − δ 2αo ∆ N up to second order GL terms in ∆ O ). A jump in the order parameters furthermore occurs when, e.g., the quartic order η-term becomes sufficiently negative.

DOS, entropy and finite temperature phase diagram
The finite temperature phase diagram is displayed in figure 1 (a). The phase boundary between the isotropic ground state and the nematic phase is continuous except near µ c1 = −0.56, where a first order transition appears, extending up to T ≈ 0.015 in temperature. Furthermore, the meta-nematic transition appearing at low temperature near µ c2 = −0.497 disappears at a slightly higher temperature. Remarkably, the slope of the phase boundary near the QCP at µ c3 = −0.47 is positive, indicating that the nematic phase has a higher entropy than the isotropic phase.
To check this, we calculated the entropy where f (E) is the Fermi function, which is also shown in figure 1. Over most of the region the entropy associated with the bare band structure is larger than the one for the nematic phase and peaks (or, rather, logarithmically diverges) at low temperature near µ = −0.54 due to the vHS in the DOS associated with the saddle points in the γ 2 band as shown in figure 1 (b). This is expected since previous work [8,10], as discussed in the introduction, showed that the formation of the nematic phase helps to avoid a large DOS near the Fermi level, and hence a large entropy, by splitting a logarithmic DOS peak into two peaks located further away from the Fermi energy. One of the peaks is indeed shifted to lower µ values in the nematic phase, which is illustrated in figure 4 (see e.g. the γ 2 peak(s) in the panels for µ = −0.56 and −0.558 near the first order transition at µ c1 ). Note also that the DOS peak deriving from γ 1 band near the Fermi level clearly reacts to the formation of nematic order by splitting as well.
There is another peak in the entropy landscape at the meta-nematic transition µ c2 inside the nematic phase where the nematic entropy exceeds the entropy of the isotropic phase. This singular behaviour is unexpected and occurs since finite nematicity persists beyond µ c2 due to the coupling of nematic and orbital order as discussed in GL free energy analysis. This unusual feature is also responsible for the concave curvature of the phase boundary in the finite temperature phase diagram in figure 1 (a). In the DOS the meta-nematic transition corresponds to a sudden changeover from one side of the Fermi level to the other of the lower of the two split γ 2 van Hove peaks (cf. the sequence from µ = −0.498 to −0.492 in figure 4), while only at the QCP (µ c3 = −0.47) do the γ vHS merge into one (respectively for γ 1 and γ 2 ).

Discussion and summary
Quantitatively, the present model is based off the bilayer ruthenate compound, where the formation of a nematic phase is driven by a moderate magnetic field applied along the c-axis [2]. Although one of the tuning parameters of the present phase diagram is the chemical potential, let us attempt to relate our results to the experimental, magnetic field tuned phase diagram [12]. Within the weak coupling approach, the leading effect of a magnetic field along c-axis is to shift the location (in momentum and energy) of the relevant van Hove points in the γ bands due to the Zeeman term [11]. In the absence of SO coupling, the Zeeman term largely acts as a uniform (in momentum space), spin-dependent chemical potential as discussed in section 2. In the presence of SO coupling, which mixes up-and down-spins, the effective chemical potential shift becomes momentum dependent. The shift is not uniform but determined by the spin composition of each band at a given momentum k. As discussed, the relevant γ 2 and γ 1 vHS arise from a small region near (±π, 0) and (0, ±π) in momentum space with a large d xy admixture. The effective shift of the γ is band (where s = ± denotes the pseudospin and i = 1, 2) therefore is determined by the spin nature of the contributing d xy orbital (i.e. either c xy k↑ or c xy k↓ ). Hence we do not expect qualitative changes to our main conclusion, i.e. the existence of a QCP and a meta-nematic transition, when the tuning parameter is an external magnetic field oriented along the c-axis.
However, one may expect that for a tilted field the interplay of SO interaction, Zeeman term, and additional orbital couplings may affect the underlying band structure in a complicated manner. The effect of an arbitrary magnetic field on nematic, orbital and magnetic instabilities therefore is complex and requires a selfconsistent theory with a number of different order parameters including independent spin-up and -down nematic and orbital orderings and orbital dependent magnetization channels, so that the present result has limited applicability. Moreover, the phase boundary at low µ and low T bends towards the nematic phase, opposite to the experimental observations. However, additional bilayer effects possibly change the underlying band structure such that the curvature becomes convex due to DOS effects and may even give rise to an additional nematic QCP. Further studies of the interplay of nematicity, orbital ordering, bilayer coupling, magnetization, and SO coupling in the presence of a magnetic field (accounting for both orbital and spin effects) are highly desirable and will be presented elsewhere [27].
Naturally, critical fluctuations accompanying the nematic QCP give rise to a quantum critical fan at finite temperature. Considering DOS and band structure effects only, one can generally expect that thermodynamic quantities diverge at most logarithmically. In contrast, the experimental specific heat and entropy curves for Sr 3 Ru 2 O 7 follow a power law ∝ H c /(H − H c ) as a function of magnetic field H when approaching the nematic phase, indicating that critical fluctuations are important [12,28] and substantiating our finding. In addition to entropy, critical non-Fermi liquid response was observed early on in resistivity measurements when approaching the nematic regime from finite temperature [6], exhibiting a ∼ T α dependence with α ≈ 1.0 − 1.5. Theoretically, nematic quantum critical fluctuations in continuum and lattice systems have been investigated as well and shown to lead to non-Fermi liquid behaviour [29,30,31,26]. However, a critical fluctuation theory in the presence of a singular DOS and in particular the magneto-thermal critical behaviour in close proximity to a meta-nematic transition at low temperature is still missing. It is also worthwhile to note that magnetic susceptibility measurements will be most sensitive to the divergences at the first order transition and the meta-nematic transition inside the nematic phase, which obscures the QCP. Such a study lies beyond the scope of the present paper and will be addressed in the future. While the details of the current results are based on the band dispersion of the bilayer ruthenate, our GL analysis with multiple order parameters, describing both continuous and disguised first order (jump in the order parameters) transitions, is in principle applicable to a range of quantum phase transitions in itinerant multi-orbital materials.
In summary, we have shown that SO coupling in a multi-orbital t 2g system can have a considerable impact on the formation of a nematic phase. Although the nematic instability arises from the d xy band, it also induces a density imbalance between d yz and d xz orbitals in the presence of SO interaction, which in turn changes the nematic phase transition from first to second order, such that the FS topology remains qualitatively unchanged across the transition. Under what circumstances a continuous or a first order phase boundary is more favourable in the presence of SO coupling likely depends on the detailed band structure and energetics as discussed in our analysis of the GL energy. Furthermore, since the continuous transition does not preempt a vHS, which in the present case happens via a nearby meta-nematic transition, the entropy diverges inside the nematic phase, causing a concave curvature of the phase boundary. The meta-nematic transition also hinders the observation of the QCP at low temperature, making further experimental analysis of the thermodynamics near the nematic phase boundaries in Sr 3 Ru 2 O 7 necessary to confirm the present theoretical proposal.