Nematic and Antiferromagnetic Quantum Criticality in a Multi-Orbital Hubbard Model for Iron Pnictides

TL;DR

This paper investigates quantum criticality in a multi-orbital Hubbard model for iron pnictides, revealing a single critical point where nematic and antiferromagnetic orders emerge simultaneously, using advanced variational Monte Carlo methods.

Contribution

It introduces a microscopic multi-orbital model analysis with non-perturbative electron correlations, identifying a quantum critical point relevant to iron-based superconductivity.

Findings

Evidence for a unique quantum critical point with coexisting nematic and antiferromagnetic orders.

The critical point occurs in the bad-metal regime of the phase diagram.

Implications for understanding iron-based superconductivity are discussed.

Abstract

The extent to which quantum criticality drives the physics of iron pnictides is a central question in the field. Earlier theoretical considerations were based on an effective field theory, and the proposed realization in P-doped iron arsenides has received extensive experimental evidence. To connect the quantum critical behavior with the underlying electronic physics, it is important to analyze it within microscopic models. Here, we do so for a multi-orbital model containing both Hubbard and Hund's interactions, by a variational Monte Carlo method based on Jastrow-Slater wave functions that allow for a non-perturbative treatment of electron correlations. We find strong evidence for the existence of a unique quantum critical point, where both nematic and $(\pi,0)$ antiferromagnetic orders develop together, in the bad-metal regime of the phase diagram. Implications of our results for the iron-based superconductivity are discussed.

Equations

$$\mathcal{H}$$

$$|{\Psi }_v⟩={\mathcal{J}}_s{\mathcal{J}}_n|{\Phi }_0⟩;$$

$${\mathcal{H}}_{\mathrm{aux}}$$

$${\mathcal{J}}_n$$

$${\mathcal{J}}_s$$

$$S(q)=\frac{1}{N}\sum _{i,j}⟨{\mathbf{S}}_i\cdot {\mathbf{S}}_j⟩e^{i\mathbf{q}\cdot ({\mathbf{R}}_i-{\mathbf{R}}_j)},$$

$${\sigma }_S=\frac{1}{N}\sum _i[⟨{\mathbf{S}}_i\cdot {\mathbf{S}}_{i+x}⟩-⟨{\mathbf{S}}_i\cdot {\mathbf{S}}_{i+y}⟩],$$

$$D=\frac{1}{N}\sum _{i,\alpha }{n}_{\alpha ,i,\uparrow }{n}_{\alpha ,i,\downarrow }.$$

Full text

Introduction – Iron-based superconductivity develops in correlated bad metals near an antiferromagnetic (AFM) order Kamihara et al. (2008); Dai et al. (2012); Hirschfeld (2016); Si et al. (2016). Furthermore, the critical temperature $T_{c}$ for the onset of superconductivity has a dome-shaped dependence as a function of doping or external pressure. These features share a considerable degree of similarity with those of the high-$T_{c}$ cuprates Lee et al. (2006); Dagotto (1994); Imada et al. (1998) and heavy-fermion systems Coleman and Schofield (2005); Si and Steglich (2010). Like in the latter systems, the role of quantum criticality is an important question in the physics of both the normal and superconducting states of the iron-based systems.

From the experimental point of view, the most extensive studies on quantum criticality in the iron-based materials have taken place in the so-called $122$ family, with the parent compound BaFe2As2. This family provides a unique opportunity to study the evolution of the electronic properties over a wide range of the phase diagram. In particular, the phosphorus-substituted system, BaFe2(As1-xPx)2, is rather clean since P-substitution does not induce appreciable scattering Shishido et al. (2010); Analytis et al. (2010); van der Beek et al. (2010). Extensive experimental observations (for earlier reviews, see Ref. Abrahams and Si (2011); Shibauchi et al. (2014)) have provided strong evidence of quantum criticality in this materials series Jiang et al. (2009); Kasahara et al. (2010); Hayes et al. (2016); Hu et al. (2018) and in related iron-based systems de la Cruz et al. (2010); Kuo et al. (2016).

From the theoretical perspective, quantum criticality in iron pnictides was proposed within an effective field theory that contains both AFM and nematic order parameters Dai et al. (2009); Wu et al. (2016). This theoretical analysis provided the basis for the early proposal of realizing quantum criticality in P-doped iron arsenides Dai et al. (2009). A quantum critical point (QCP) corresponds to a special class of second-order phase transitions that takes place at zero temperature, typically in a material in which the critical temperature has been driven to zero by non-thermal parameters Sachdev (2011); Coleman and Schofield (2005); Si and Steglich (2010). Whether a putative QCP exists in the phase diagram of the iron-pnictide materials or not remains debated. Nevertheless, it is widely believed that a QCP is hidden beneath the superconducting dome; deepening our understandings of the physics regarding a QCP may then shed light on the mechanism of unconventional superconducting pairing. In order to address the physics of these compounds, ab-initio calculations based on the standard local-density approximation have been performed Lebègue (2007); Singh and Du (2008); Cao et al. (2008). While this approach was useful to highlight the rather complex band-structure, with several $d$ orbitals participating in the low-energy properties, it does not properly treat the effects of strong correlations and electron pairing.

To address the physics of quantum phase transitions in microscopic studies requires a systematic treatment of correlated multi-orbital models with both Hubbard and Hund’s interactions. The bad-metal nature of the normal state implies that the interactions are of the same order as the kinetic energy and, thus, the need to treat correlated states beyond perturbative approaches. In particular, Hund’s interactions are important for the correlation effects in this regime. Additionally, the entwining of AFM and nematic correlations highlights the need to address spatial correlations.

In this letter, we aim at identifying the QCP and studying the physics around it within a multi-orbital model that contains both Hubbard (intra-orbital $U$ and inter-orbital $U^{\prime}$) and Hund’s ($J_{H}$) interactions. We perform variational Monte Carlo (VMC) simulations with Jastrow-Slater states. While VMC analyses have previously been carried out in multi-orbital Hubbard models Yang et al. (2011); Watanabe et al. (2012); Tocchio et al. (2016); De Franco et al. (2018), in our study we incorporate a spin Jastrow factor, in order to improve the description of correlation effects that originate from the Hund’s term. At large electron-electron interactions, we find a stable AFM order with pitch vector ${\bf Q}=(0,\pi)$ [or equivalently $(\pi,0)$] coexisting with a nematic order, i.e., a spontaneous breaking of $\pi/2$ rotations of the lattice. As the interaction decreases, a unique phase transition takes place, where both antiferromagnetism and nematicity disappear. This represents the most important result of the present work. Furthermore, by analysing the behavior of the double occupation, we identify the location of the metal-insulator (Mott) transition, which appears not too far above the AFM-nematic QCP.

Models and Methods – Theoretically, in the iron pnictides, there are six electrons occupying the $3d$ Fe orbitals. In these materials, it is fundamental to include this multiplicity of the involved orbitals Hirschfeld (2016); Si et al. (2016). The most generic Hamiltonian that describes $N_{\rm orb}$ orbitals with inter- and intra-orbital interactions is given by:

[TABLE]

where $c^{\phantom{{\dagger}}}_{\alpha,i,\sigma}$ is the electronic annihilation operator with spin $\sigma$ on orbital $\alpha$ and site $i$, and the density operator is $n_{\alpha,i,\sigma}=c^{\dagger}_{\alpha,i,\sigma}c^{\phantom{{\dagger}}}_{\alpha,i,\sigma}$, with $n_{\alpha,i}=\sum_{\sigma}n_{\alpha,i,\sigma}$; moreover, ${\bf S}_{\alpha,i}=\frac{1}{2}\sum_{\sigma,\sigma^{\prime}}c^{{\dagger}}_{\alpha,i,\sigma}\tau_{\sigma,\sigma^{\prime}}c^{\phantom{{\dagger}}}_{\alpha,i,\sigma^{\prime}}$, $\tau$ being the Pauli matrices, are the spin operators. Here, $t^{\alpha,\beta}_{i,j}$ are the intra- ($\alpha=\beta$) and inter-orbital ($\alpha\neq\beta$) hopping amplitudes, respectively; $\Delta^{\rm cf}_{\alpha}$ denotes the crystal field; $U$ and $U^{\prime}=U-2J_{H}$ are the intra- and inter-orbital on-site Hubbard interactions, while $J_{H}$ is the Hund’s coupling Georges et al. (2013). Hopping amplitudes and crystal fields can be estimated by density functional theory through a comparison with the experimental data of the Fermi surface. While a full description of the $3d$ shell of Fe orbitals would imply $N_{\rm orb}=5$ Hirschfeld (2016); Si et al. (2016), there are simplified versions restricting to $N_{\rm orb}=3$ or $2$ orbitals. While the former one includes all the $t_{2g}$ orbitals (with $xz$, $yz$, and $xy$ symmetry) Daghofer et al. (2010), the latter one restricts to the case in which only $xz$ and $yz$ orbitals are kept in the low-energy description of iron pnictides Raghu et al. (2008); Moreo et al. (2009). We mention that the 3-orbital model contains two Fe atoms per unit cell, but full translational invariance can be recovered after considering $c^{\dagger}_{\alpha,i,\sigma}\to(-1)^{{\rm R}_{i}}c^{\dagger}_{\alpha,i,\sigma}$ for $\alpha=xz$ and $yz$ (in the 2-orbital case, this transformation does not spoil translational invariance). For the 2-orbital model, there are no crystal fields and the actual values of the hopping amplitudes are $t^{1,1}_{i,i+x}=t^{2,2}_{i,i+y}=-1.3$, $t^{1,1}_{i,i+y}=t^{2,2}_{i,i+x}=1.0$, $t^{\alpha,\alpha}_{i,i+x+y}=t^{\alpha,\alpha}_{i,i+x-y}=-0.85$, and $t^{1,2}_{i,i+x+y}=-t^{1,2}_{i,i+x-y}=0.85$ (with symmetric amplitudes when exchanging orbital indices) Raghu et al. (2008); Moreo et al. (2009). Instead, for the 3-orbital model, we take the parameters given in Table I of Ref. Daghofer et al. (2010).

Within the VMC approach, different trial wave functions can be considered, in order to optimize the expectation value of the total energy. In this letter, we will choose the Jastrow-Slater wave functions, defined by Capello et al. (2005):

[TABLE]

here, the uncorrelated state $|\Phi_{0}\rangle$ is specified by the following auxiliary (quadratic) Hamiltonian Tocchio et al. (2016); De Franco et al. (2018):

[TABLE]

where ${\tilde{t}}^{\alpha,\beta}_{i,j}$, $\mu^{\alpha}$, and $\Delta^{\rm AFM}_{\alpha}$ are variational parameters. In particular, the presence of $\Delta^{\rm AFM}_{\alpha}\neq 0$ implies magnetic order on each orbital. By choosing ${\bf Q}=(\pi,\pi)$ or $(0,\pi)$ [or equivalently $(\pi,0)$], we can have either the Néel order or the collinear AFM order (CAFM). To allow for the treatment of the correlation effect associated with the Hund’s coupling, we go beyond previous studies of multi-orbital Hubbard models Yang et al. (2011); Watanabe et al. (2012); Tocchio et al. (2016); De Franco et al. (2018) by incorporating the spin Jastrow factor in addition to the usual density one:

[TABLE]

where $u^{\alpha,\beta}_{i,j}$ and $v^{\alpha,\beta}_{i,j}$ are translationally invariant variational parameters, which introduce quantum fluctuations into the wave function. In the context of the one-band Hubbard model, how the intersite spin Jastrow factor better captures the spin-spin correlations has been discussed before Becca et al. (2000). In Eq. (3), we have incorporated the magnetic order parameter through $\Delta_{\alpha}^{\rm AFM}$, which couples to the spins perpendicular to the $z$ direction. This improves the description of the ordered magnetic state, given that the spin fluctuations produced by the spin Jastrow factor are orthogonal to the direction of the ordered moment in the uncorrelated wavefunction Hu et al. (2015a).

We use the stochastic reconfiguration optimization method Sorella (2005) to optimize the variational parameters in both the auxiliary Hamiltonian and the Jastrow factors, in order to find the energetically favored state in the VMC scheme. In the following, we allow for states that break rotational symmetry in spin and lattice, but always keep the translational symmetry. We emphasize that the optimization includes the effective hopping amplitudes ${\tilde{t}}^{\alpha,\beta}_{i,j}$, which is also allowed to break the C4 rotational symmetry and this is crucial to realize the nematic order. In addition, the Jastrow factor is optimized up to the next-nearest neighbor in order to limit the number of variational parameters. In the following, we denote with $L$ the linear dimension of the system, corresponding to a $N=L\times L$ total number of sites.

Results – We have performed VMC simulations for both the 2- and the 3-orbital Hubbard model at half filling not with the Hund’s coupling $J_{H}=0.1U$, and considered both the Néel ordered phase and the CAFM one in the auxiliary Hamiltonian of Eq. (3). After optimization, we find that the CAFM order has lower energy than the Néel one, especially in the large $U$ region, for both the 2- and the 3-orbital Hubbard model. Therefore, in the following, we show numerical results only for the variational wave function with the CAFM order.

In order to identify the magnetic order, we compute the spin structure factor as

[TABLE]

where $\langle{\bf S}_{i}\cdot{\bf S}_{j}\rangle$ is the spin-spin correlation function. In addition, we introduce the nematic order parameter $\sigma_{S}$ to assess the rotational symmetry breaking (i.e., the bond nematicity):

[TABLE]

The results are reported in Figs. 1 and 2 for the 2- and the 3-orbital Hubbard model, respectively. For large electron-electron interactions, a sharp peak in the spin structure factor is clearly visible at ${\bf Q}=(0,\pi)$, typical of CAFM order. In order to assess the presence of a true long-range order, we define the magnetic order parameter as $m^{2}_{S}=S(0,\pi)/N$, which is also reported in Figs. 1 and 2, as a function of the interaction $U$ on different lattice sizes for both the 2- and the 3-orbital Hubbard model. It is clearly shown that a magnetic phase transition occurs between a CAFM ordered phase and a nonmagnetic one in both cases. According to the finite size scaling, we locate the magnetic transition points at $U_{c}\approx 7.5$ for the 2-orbital Hubbard model and at $U_{c}\approx 1.0$ for the 3-orbital one. Considering the half-bandwidth $W=6$ for the 2-orbital Hubbard model Raghu et al. (2008); Moreo et al. (2009) and $W=1.5$ for the 3-orbital one Daghofer et al. (2010), we observe that the magnetic critical points are located at $U_{c}/W\approx 1.25$ for the 2-orbital Hubbard model and $U_{c}/W\approx 0.666$ for the 3-orbital one.

While a CAFM order with pitch vector ${\bf Q}=(0,\pi)$ would naturally imply nematic order, since correlations along $x$ and $y$ directions are fundamentally different, a nematic order may be possible even without CAFM order. In order to clarify this important issue we report the behavior of the nematic order parameter $\sigma_{S}$ as a function of the interaction $U$, for both the 2- and the 3-orbital Hubbard models, see Figs. 1 and 2. Our numerical calculations indicate that the nematic phase transition is concurrent with the magnetic one, strongly suggesting that a pure nematic phase is hardly obtained at zero temperature.

Within the multi-orbital models that have been considered so far, the presence of a nonzero AFM order does not necessarily imply the presence of a charge gap. Therefore, it is important to assess the metallic/insulating behavior of the variational wave function. We consider the density of double occupation:

[TABLE]

The numerical calculations are shown in Fig. 3 for both the 2- and the 3-orbital Hubbard model. At $U=0$, the doublon density is $0.5$ for the 2-orbital model and $\approx 0.75$ for the 3-orbital one. When including the Hubbard $U$, $D$ monotonically decreases, showing a small kink connected to the development of AFM order and a bending that is directly associated to the opening of the charge gap Tocchio et al. (2011), namely to the metal-insulator transition. This can be located around $U=11$ and $U=2$ for the 2- and 3-orbital model, respectively. Since these values are not too much larger than $U_{c}$ for the AFM-nematic QCP, the system at $U_{c}$ falls in the bad-metal regime. In addition, our analysis highlights the possibility to have an intermediate metallic phase with both AFM and nematic orders, sandwiched between a paramagnetic metal at small values of $U$ and an AFM insulator with nematic order at strong coupling.

Finally, we report the effective filling of each orbital in both cases, see Fig. 4. In the 2-orbital model, both $xz$ and $yz$ are close to half filling for every value of the electron-electron interaction; instead, we find that in the 3-orbital model the filling number of the $xy$ orbital is different from the one of the $xz$ and $yz$ orbitals, when $U\lesssim 2$, i.e., in the metallic regime.

Discussions – Several points of considerations are in order. First, similar to the Hubbard-Heisenberg model of the one-band case Daul et al. (2000), we can introduce $J_{1}-J_{2}$ interactions into the multi-orbital Hamiltonian Eq. (1). We find the results to be qualitatively similar, including the concurrence of the nematic and $(\pi,0)$ AFM transitions.

Second, it is instructive to compare the result of our microscopic studies with that of an effective field theory for the AFM and nematic orders. There, the coupling between the two degrees of freedom originates from a quartic interaction of the underlying magnetic field. The relevant nature of this coupling turns the concurrent transition to be ultimately first order. However, because the coupling is only marginally relevant, it is a weakly first order transition: The jump of the order parameter is very small compared to the saturated moment, leaving an extended dynamical range for quantum criticality Dai et al. (2009). The same conclusion is derived by a saddle-point calculation in a large-$N$ limit of the field theory Wu et al. (2016), and is consistent with experimental observations Shibauchi et al. (2014); Hayes et al. (2016); Hu et al. (2018, 2015b). In our microscopic study, the numerical results point to a unique QCP for the transition of both orders; however, our data do not rule out the quantum phase transition to be weakly first order, with the order parameters experiencing a small jump.

Third, given the purely two-dimensional feature of the model, an infinitesimal temperature completely destroys magnetic order, while nematicity is expected to survive up to a critical temperature $T_{\rm n}$, since it is related to a discrete symmetry breaking. In a three-dimensional lattice with nonzero interlayer couplings, the AFM order is also stable against thermal fluctuations, and the critical temperature $T_{\rm af}$ will not be linked to $T_{\rm n}$, possibly allowing for a pure nematic order, with no long-range antiferromagnetism.

Conclusions – We have studied both the 2- and 3-orbital Hubbard models by means of the variational Monte Carlo method. Using optimized Jastrow-Slater wave functions, we have calculated both magnetic and nematic order parameters showing the existence of a unique quantum critical point separating paramagnetic and antiferromagnetic (with broken rotational symmetry) phases. Indeed, the size scaling analysis of both magnetic and nematic order parameters indicated the concomitant insurgence of these orders, limiting the possibility to have a genuine nematic ground state (with no magnetic order) to a tiny region, which is not detected in our calculations. Furthermore, the behavior of the doublon density suggests that the metal-insulator transition is located slightly inside the antiferromagnetic phase, indicating that quantum criticality occurs in the bad-metal regime and that a metallic phase with the CAFM order exists in both the 2- and the 3-orbital Hubbard model. Our results provide the theoretical basis to understand the quantum criticality observed in the experiments on the iron pnictide BaFe2As2 Jiang et al. (2009); Kasahara et al. (2010); Hayes et al. (2016); Hu et al. (2018). In addition, the concurrence of the antiferromagnetic and nematic transitions suggest that both type of quantum criticality will strongly influence the development of the dome of high temperature superconductivity in its vicinity.

Acknowledgement.— The work was supported in part by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-SC0018197 (W.-J.H., H.-H.L., and Q.S.), the Robert A. Welch Foundation Grant No. C-1411 (H.H., and Q.S.), a Smalley Postdoctoral Fellowship of the Rice Center for Quantum Materials (H.-H. L.), the National Science Foundation of China Grant number 11674392 and Ministry of Science and Technology of China, National Program on Key Research Project Grant number 2016YFA0300504 (R.Y.). The majority of the computational calculations have been performed on the Shared University Grid at Rice funded by NSF under Grant EIA-0216467, a partnership between Rice University, Sun Microsystems, and Sigma Solutions, Inc., the Big-Data Private-Cloud Research Cyberinfrastructure MRI-award funded by NSF under Grant No. CNS-1338099 and by Rice University, the Extreme Science and Engineering Discovery Environment (XSEDE) by NSF under Grant No. DMR160057. Q.S. acknowledges the hospitality and support by a Ulam Scholarship of the Center for Nonlinear Studies at Los Alamos National Laboratory and the hospitality of the Aspen Center for Physics (NSF grant No. PHY-1607611).

I The filling number of each orbital for iron pnictides

To estimate the electron filling in each orbital, we have performed calculations by using the U(1) slave-spin method Yu and Si (2012) on 5-orbital Hubbard models for iron pnictides LaOFeAs and BaFe2As2 with six electrons occupying the $5$ $3d$-orbitals of each Fe ion. The results for the Ising-only Hund’s coupling $J_{H}/U=0.25$ for LaOFeAs and BaFe2As2 are shown in Fig. 5(a,b). For large interactions in each model, one of the two $e_{g}$ orbitals is doubly occupied ($x^{2}-y^{2}$ orbital for LaOFeAs and $z^{2}-r^{2}$ orbital for BaFe2As2), while all the other orbitals are at half filling. Most importantly, for both models, at small and intermediate interactions, the filling numbers of the three $t_{2g}$ orbitals $xz$, $yz$, and $xy$ are close to half filling. These results are consistent with those in our VMC calculations shown in Fig. 4 and justifies setting the total filling number to half-filling in the 3-orbital Hubbard model in the VMC calculations. We have used a larger Ising-only Hund’s coupling to mimic the effect of SU(2)-invariant Hund’s coupling of $J_{H}/U=0.1$, but our conclusion is not sensitive to the precise value of $J_{H}/U$. We illustrate this point in Fig. 5(c), with a choice of the Ising-only $J_{H}/U=0.1$; here again, the filling numbers of the $xz$, $yz$, and $xy$ orbitals are close to half filling for small and intermediate interactions.

II Additional numerical results by variational Monte Carlo

In Fig. 6, we present the finite size scaling for the magnetic $m^{2}_{S}$ and the nematic $\sigma_{S}$ order parameters at different interactions $U$ in both the 2- and 3-orbital Hubbard model with the Hund’s coupling $J_{H}/U=0.1$. For the 2-orbital Hubbard model, at small interactions, i.e., $U<7.5$, the finite size scaling show vanishing values of $m^{2}_{S}$ and $\sigma_{S}$ in the thermodynamic limit; at $U\simeq 7.5$, small but nonzero $m^{2}_{S}$ and $\sigma_{S}$ suggest concurrent magnetic and nematic phase transitions. For the 3-orbital Hubbard model, nonzero $m^{2}_{S}$ and $\sigma_{S}$ after the finite size scaling are obtained around $U=1.0$, suggesting again concurrent magnetic and nematic phase transitions. We notice that the extrapolation to the thermodynamic limit is remarkably smooth, except in the metallic phase with CAFM order, where stronger fluctuations of the order parameters with the system size are observed. Nevertheless, a finite value of the order parameters in the thermodynamic limit, monotonically increasing with the interaction $U$, can be clearly obtained also in the metallic phase with CAFM order.

For the 3-orbital Hubbard model, we have calculated additional data for the filling numbers of each orbital, see Fig. 7. It clearly shows the difference between the behavior of the orbital $xy$ from orbitals $xz$ and $yz$, especially at small and intermediate interactions.