Non-Fermi-liquid at (2+1)d ferromagnetic quantum critical point

We construct a two-dimensional lattice model of fermions coupled to Ising ferromagnetic critical fluctuations. Using extensive sign-problem-free quantum Monte Carlo simulations, we show that the model realizes a continuous itinerant quantum phase transition. In comparison with other similar itinerant quantum critical points (QCPs), our QCP shows much weaker superconductivity tendency with no superconducting state down to the lowest temperature investigated, hence making the system a good platform for the exploration of quantum critical fluctuations. Remarkably, clear signatures of non-Fermi-liquid behavior in the fermion propagators are observed at the QCP. The critical fluctuations at the QCP partially resemble Hertz-Millis-Moriya behavior. However, careful scaling analysis reveals that the QCP belongs to a different universality class, deviating from both (2+1)d Ising and Hertz-Millis-Moriya predictions.


I. INTRODUCTION
Understanding the behavior of gapless fermionic liquids in the vicinity of a quantum phase transition is at the heart of strongly correlated electron systems, dating back to the celebrated Hertz-Millis-Moriya framework [1][2][3]. In particular, the question of Fermi-liquid instabilities at a magnetic (quantum) phase transition [4][5][6][7][8][9] and its applications to heavyfermion materials [10,11] and transition-metal alloys (such as cuprates and pnictides [12,13]), is of vital importance and broad interest to the condensed matter and high energy physics communities.
On the other hand, to be able to obtain the understanding of quantum critical phenomena in itinerant electron systems is extremely challenging. Recently, extensive research efforts have been devoted to this question utilizing advanced renormalization group analysis, including the 2d Fermi surface coupled to a U(1) gauge field, to Ising-nematic or to spin-density-wave bosonic fluctuations [14][15][16][17][18][19][20][21][22]. Other approaches include dimensional regularization [23,24] and working in the limit of a large number of boson flavors [25]. Although important insights have been revealed from these studies, many fundamental questions still remain open. For example, will anomalous dimensions arise at such an itinerant quantum critical point? The Hertz-Millis-Moriya theory predicts meanfield scalings. However, utilizing the effective field theory derived in Ref. [17], it was shown that divergence at fourloop order can lead to anomalous dimensions deviating from mean-field exponents [21]. In addition, the stability of these quantum critical points is also an open questions, e.g., will non-analyticities in the momentum and frequency dependence of the theory leads to an instability towards a fluctuationinduced first-order transition or a fluctuation-induced secondorder transition [26,27]? Moreover, various studies have suggested that near such a quantum critical point, critical fluctuations may trigger some other instability, resulting in a new phase that covers the QCP and masks the quantum critical region [28]. Experimentally, a ferromagnetic QCP was reported in a heavy-fermion metal [11].
Since analytical approaches are facing difficulties to study itinerant QCPs in a controlled manner, here, we study such a QCP utilizing unbiased numerical calculations, i.e., the determinant quantum Monte Carlo (QMC) technique, which has been demonstrated as a very effective tool for such problems [29][30][31][32][33][34][35][36][37]. For the numerical studies on QCPs in itinerant fermionic systems, one major challenge lies in the fact that critical quantum fluctuations often trigger strong effective attractions between fermions, resulting in instabilities in the particle-particle channel. For example, in the recent studies on Ising nematic and charge-density-wave (CDW) QCPs, superconducting domes are observed covering the QCPs [29][30][31][32][33]35]. Although it is a very intriguing phenomenon that a QCP can induce unconventional superconductivity, for the study of quantum criticality itself, the induced superconducting dome makes it difficult to obtain direct information about the critical point for the following two reasons. (1) To measure accurately the critical exponents, it is important to examine the close vicinity of the QCP, which becomes challenging if the QCP is buried inside a new quantum phase. (2) It is known that quantum criticality is one source for non-Fermi liquid behaviors. In order to understand and characterize such a non-Fermi liquid induced by a QCP, it is important to suppress other ordering in the quantum critical region. As a result, obtaining a pristine QCP without other induced instabilities becomes crucial for studying these QCPs, which is one main objective of this paper.
In this paper, we construct a model of 2d fermions interacting with gapless ferromagnetic (Ising) fluctuations, and use the determinant QMC technique to solve this (2+1)d problem exactly [29][30][31][32][33][34][35][36][37]. Our QMC results are consistent with such a pristine continuous quantum phase transition in itinerant ferromagnet in (2+1)d, with no superconducting ordering at any coupling strengths and down to the lowest temperature that we can access. The absence of superconductivity allows us to study the close vicinity of the QCP, where we found clear signatures of non-Fermi-liquid behavior in the fermion propagators, induced by critical fluctuations. Furthermore, we find that due to the coupling between fermions and bosonic modes, the ferromagnetic QCP is different from an ordinary (2+1)d Ising transition [38], but also deviates from the predictions of the Hertz-Millis-Moriya theory [6]. Hence, our results support the existence of a ferromagnetic QCP with markedly non-Fermi liquid behavior. These results broaden the theoretical understanding of itinerant quantum criticality and make connections to the existing experiment phenomena.

II. MODEL AND METHOD
We consider a two-dimensional lattice model of itinerant fermions coupled to an Ising ferromagnet with a transverse field [ Fig. 1(a)]. The Hamiltonian is comprised of three parts, (1) The fermionic part,Ĥ f , describes spin-1/2 fermions on a square lattice with two independent orbitals per site (λ = 1, 2). It includes a nearest-neighbor hopping term that preserves the spin and orbital sub-indices, where i and j labels the sites, while σ and λ are the spin and orbital indices, respectively. We work in the grand canonical ensemble, where the fermion density is set by the chemical potential µ. All energy scales are measured in units of t. In addition, each site of the square lattice has an Ising spinŝ z i = ±1, whose quantum dynamics are governed by a ferromagnetic transverse field Ising model H s . The Ising spins and the fermions are coupled via an on-site Ising termĤ s f , whereσ z iλ = (n iλ↑ −n iλ↓ )/2 is the z component of the fermion spin at orbit λ on site i.
While generically, models describing ferromagnetic transitions in fermionic systems suffer from the sign-problem, here, the absence of the sign-problem is guaranteed by the introduction of the two orbitals λ = 1 and 2. The two-orbital model is invariant under the anti-unitary symmetry iτ y K, where τ y is a Pauli matrix in the orbital basis and K is the complex conjugation operator, and is thus sign-problem-free [39]. The details of the DQMC implementation are presented in Appendix A.
The Hamiltonian in Eq.(1) possesses a SU(2) × SU(2) × U(1) × U(1) × Z 2 symmetry, where the two SU(2) symmetries are independent rotations in the orbital basis for spin up and spin down, the two U(1) symmetries correspond to conservation of particle number with spin up and spin down, and the Z 2 symmetry interchanges spin up and spin down while flipping the Ising spins s z → −s z .
As h and T are reduced, the system undergoes a paramagnetic-ferromagnetic (PM-FM) phase transition, spontaneously breaking the Z 2 symmetry. In the absence of coupling between the Ising spins and the fermions, the transition belongs to the Ising universality class [40]. However, in the presence of the coupling ξ between the Ising spins and the fermions, the system becomes strongly correlated near the ferromagnetic QCP. Here, we focus on the properties of this exotic itinerant paramagnetic-ferromagnetic transition.

III. PHASE DIAGRAM
For ξ = 0, the Ising spins are decoupled from the fermions and the phase transition is governed by the transverse-field Ising model. The phase diagram is shown in Fig. 1(b), where the phase boundary is marked by grey data points. Our numerical studies confirm that the phase boundary ends at a QCP at T = 0 with (2+1)d Ising universality class [40]. Near the QCP, the transition temperature follows the scaling relation, (2) and νz = 0.63, consistent with the literature [34,38].
We now study the itinerant PM-FM transition by turning on the coupling between the fermions and the Ising spins. We begin by setting the coupling strength ξ = 1 and chemical potential µ = −0.5, which gives rise to a fermion density n iλ ≈ 0.8. As shown in Fig. 1(b), turning on the coupling shifts FM phase boundary (orange data points) to higher values of T and h. At this coupling strength, down to the lowest temperature that we have accessed, β = 100 (T = 0.01), we observe no signature of any additional ordered phases near the QCP. We identify the finite-temperature ferromagnetic transition by a finite-size scaling analysis of spin susceptibilities, as explained in Appendix C. Extrapolation towards zero temperature indicates that the itinerant PM-FM quantum phase transition occurs at h c = 3.270 (6), and is of second order, but the scaling behavior near the QCP deviates strongly from that of the (2+1)d Ising universality class. For example, as shown in Fig. 1 (b), the transition temperature T N scales as T N (h) ∼ |h − h c | c with h c = 3.270(6) and c = 0.77(4). Note that due to the itinerant nature of the QCP, the exponent c is no longer expected to obey the relation c = zν [6].
Due to the non-zero coupling ξ, the Fermi surface structure changes across the FM transition. The fermionic low-energy spectral weight shown in the inset of Fig. 1(b), extracted from the imaginary time Green's function G(τ = β/2) [30,41], reveals the location of Fermi surface. In the PM phase (h > h c ), fermions with up and down spins share the same Fermi surface, due to the spin degeneracy (Ising symmetry) of the Hamiltonian [right inset of Fig. 1(b)]. This degeneracy is lifted in the FM phase (h < h c ) due to spontaneous symmetry breaking, and thus the Fermi surface splits [left inset of Fig. 1(b)].
At the QCP, low-energy fermionic excitations near the Fermi surface become strongly coupled with the critical bosonic spin fluctuations, offering an ideal platform for investigating the itinerant FM-QCP. Below, we demonstrate that this QCP is stable down to low energy scales, and more interestingly, that it is characterized by a dramatic breakdown of Fermi liquid behavior.

IV. NON-FERMI LIQUID BEHAVIOR
One of the key theoretical prediction for an itinerant QCP is that the critical bonsonic fluctuations induce strong damping of the fermions, resulting in a non-Fermi liquid where the low-temperature fermionic quasi-particle weight vanishes at the Fermi surface [14, 16-18, 22-24, 42]. Following [43], we measure the quasi-particle weight from the Matsubarafrequency self-energy where Σ is obtained from our QMC simulations, and k F is the momentum at FS. In our finite temperature simulations, the Matsubara frequencies take discrete values, ω n = π(2n + 1)T.
As an estimator for Z k F , we use the first Matsubara frequency, The results are shown in Fig. 2(a) for h = h c (squares) and h > h c (circles). Here, we plot the quasi-particle weight at two different momentum points on the Fermi surface, i.e. k F along the k x direction (θ = 0) and along the k x = k y direction (θ = π 4 ). In the paramagnetic phase (h > h c ), Z(T) remains close to unity at low temperature, indicating welldefined quasi-particles on the Fermi surface, as expected in a Fermi liquid. In contrast, at the QCP (h = h c ), Z(T) is suppressed with decreasing temperature and extrapolates to zero at T → 0, which is the key signature of a non-Fermi liquid. This is one of the key findings of this study.
As shown in Fig. 2(a), at the QCP, as T → 0, Z(T) decreases faster for k F along the diagonal direction (θ = π 4 ) in comparison with other parts of the FS (e.g. θ = 0). This anisotropy is due to the anisotropy of the Fermi surface. The fermion density in our simulation is around 0.8, which is not far from perfect nesting (at half-filling). This near-nesting Fermisurface geometry results in soft fermion-bilinear modes around Q = (π, π), which we have directly observed by measuring the fermion spin susceptibility as shown in the Appendix H. These soft modes have a finite gap, and thus they are not the origin of the non-Fermi liquid behavior. However, scattering with these soft modes presumably introduces additional dampings for the fermions. Such scattering processes are stronger (weaker) for k F near the diagonal direction θ = π 4 (θ = 0), where 2k F is close to (far away from) Q, and thus can lead to the observed anisotropy in Z(T).
In Fig. 2 (b), we show the imaginary part of the selfenergy, Σ(ω n ) = G −1 0 − G −1 (where G(ω n ) and G 0 (ω n ) are the fermionic Green's function, and the Green's function of the non-interacting system, respectively) as a function of the Matsubara frequency ω n . In the paramagnetic phase (circle symbols), −Im(Σ) approaches zero linearly as ω n → 0, as expected for a Fermi liquid. Such a behavior is not seen at the QCP, however, where we observe an increase of −Im(Σ) upon decreasing ω n , indicating a strong damping of the fermions at low frequencies. One possible mechanism for such surprising frequency dependence is related to thermal fluctuations of the FM order parameter, as shown in appendix I.
The fermion Green's functions, ImG(k F , ω n ), are presented in Fig. 2 The fact that −Im(G) increases with decreasing ω n indicates that both for h = h c and h > h c and to the lowest temperature we have accessed, there is no visible suppression of the fermionic spectral weight at low frequency, i.e. no sign of an opening of a gap in the fermionic spectrum (at least down to frequencies of the order of ω ∼ πT). A sim-  2 ) at FS for T = 0.05 while the right inset is for T = 0.1. Although there is anisotropy in Z k F at different parts of the FS, the quasi-particle weight in k x and k x = k y directions both approach zero at FM-QCP, indicating a non-Fermi liquid behavior for the entire FS. The data in the PM phase shows the quasi-particle weight approaching a constant (very close to 1), indicating the system is a Fermi liquid. (b) −Im(Σ(k F , ω n )) at FM-QCP (h = h c , square symbol), it increases as ω n → 0 -signifying the system at QCP loses their quasi-particle weight with a power law -a non-fermi liquid behavior, while in PM phase (h = 3.60, circle symbol) the imaginary part of self energy approaches zero linearly as ω n → 0 -a fermi liquid behavior. (c) Imaginary part of the single-fermion Green's function at the FM-QCP (h = 3.27, square symbol) and in PM phase (h = 3.60, circle symbol). No signature of gap formation is observed. ilar conclusion can be reached by noting that −ImΣ is never much larger than ω 0 , i.e the self-energy never dominates over the bare frequency dependence of the Green's function. This observation is consistent with the fact that there is no signature of a nearby superconducting phase for any value of h, as shown in Sec. VI and Appendix G. Similarly, there is no signature of any other competing phase that emerges close to the QCP (see Appendix H).

V. QUANTUM CRITICAL SCALING ANALYSIS
Near the quantum critical point, the bosonic critical modes become strongly renormalized by the coupling to the gapless fermionic degrees of freedom. As a result, the universality class of the quantum critical point is different from that of an ordinary Ising transition in (2 + 1) dimensions. Our QMC results indicate that the behavior at the QCP resembles the be-havior predicted by Hertz-Mills theory, but also deviates from it in a significant way. Below, we first summarize the Hertz-Millis predictions and then present a modified Hertz-Millis scaling formula, which fits our QMC data for the Ising spin susceptibility at all the momenta and frequencies simulated.
The Hertz-Millis-Moriya theory is based on quantum dynamics obtained from the random phase approximation (RPA) [1][2][3]. Within this approximation, the Ising spin susceptibility, , takes the following form near the QCP, where c t , c h , c q , c ω are constants. Here, the c q q 2 +c ω ω 2 comes from the bare action of the Ising degrees of freedom, and the ∆(q, ω n ) term is the contribution of the fermionic fluctuations.
For an isotropic 2D Fermi fluid, and for q and ω n much smaller than the Fermi momentum and energy, respectively, where c HM is a constant and v f is the Fermi velocity. A key property of ∆(q, ω n ) is its singular behavior in the limit q → 0, ω → 0. Depending on whether one first takes the long-wavelength (q → 0) limit or the low-energy (ω → 0) limit, ∆(q, ω n ) converges to different values, lim q→0 lim ω n →0 ∆(q, ω n ) ∼ ω n /q → 0, while lim ω n →0 lim q→0 ∆(q, ω n ) = c HM . This singularity is of great importance for the quantum dynamics at the QCP. It leads to the following property: Beyond the RPA level, the scaling relation above can be modified by higher order terms, and the scaling analysis in Ref. 2 suggests that the exponent for T shifts from 2 to 1 up to logarithmic corrections. Although these two relations differ by a constant (c HM ), note that the q and ω n dependence has the same scaling exponent. The reason for this behavior is that at q = 0 and small ω n (or at ω n = 0 and small q), ∆(q, ω n ) becomes independent of q or ω n , respectively. Thus the ω n or q dependence in χ −1 is dominated by the bare action of the Ising degree of freedom. The fact that the exponents characterizing the ω and q dependence are the same reflects the emergent Lorentz symmetry of the bare Ising action.
Our QMC results share some characteristics with this predicted form, but with anomalous scaling dimensions.
As shown in Figs. 3(a) and (b), we indeed find that lim q→0 lim ω n →0 χ −1 (q, ω n ) differs from lim ω n →0 lim q→0 χ −1 (q, ω n ) by a constant c HM = 0.20(4) as predicted by the RPA. The ferromagnetic susceptibility is found to be well-described by the following formula: with an anomalous exponent a q = 1.85 (3) . This is different both from the exponent for an Ising transition in (2 + 1) dimensions, 1.96, and from the Hertz-Millis value of 2. The presence of an anomalous exponent is another key finding of this study. The functional form in Eq. (10) is analogous to the RPA prediction [Eq. (6)], but allows for non-mean-field exponents (a t , γ and a q ). This form is found to fit all the numerical data points [ Fig. 3(c)].
Note that Eq. (10) contains ∆(q, ω n ), which has the form of a free-fermion susceptibility. However, as will be shown below, within numerical resolution, as long as ∆ captures the singular behavior at q → 0 and ω → 0, the quality of the fit of χ(h, T, q, ω n ) is not sensitive to the detailed functional form of ∆(q, ω n ). More importantly, our key conclusions, e.g., the anomalous dimension η = 2−a q = 0.15, are fully independent of the particular choice of ∆(q, ω n ). Therefore, we use here the RPA form [Eq. (7)] where, for simplicity, we set v f to 2.
By fitting with χ(h, T) at q = 0 and ω n = 0, we find that c t = 0.13(1), a t = 1.48(4), c h = 0.7(1), γ = 1.18(4) , details can be found in Appendix D. As shown above, by fitting with χ at q = 0 or ω = 0 [ Figs. 3(a) and (b)], we find a q = 1.85(3), c q = 1.00(2), c ω = 0.10(2) and c HM = 0.20 (4) . With all the fitting parameters fixed, we can use Eq. 10 to collapse all χ data for all q, ω n , h and T. The results are present in Fig 3(c). Clearly, all the data collapse onto the same curve, especially at small q, ω n , low temperature T and h ∼ h c .
For ∆(q, ω n ), we find that as long as lim ω n →0 ∆(q, ω n ) = 0 and lim q→0 ∆(q, ω n ) = c HM = 0.20(4), modifying the functional form of ∆(q, ω n ) has little impact on χ within numerical error bars of the QMC data. This uncertainty in ∆(q, ω n ) prevents us from obtaining detailed information about the dynamics at the QCP (e.g., the precise value of the dynamical critical exponent z). However, since by definition ∆(q, ω n ) vanishes at ω n = 0, all the conclusions regarding the static limit (ω n = 0), including the anomalous dimensions η = 2 − a q = 0.15 obtained from the q-dependence in χ(q, ω n = 0) are independent of the details of ∆(q, ω n ).

VI. SUPERCONDUCTIVITY
We now turn to discuss the superconducting properties close to the FM-QCP. Unlike previous QMC studies of quantum criticality in metals, such as the Ising-nematic, SDW and other QCPs [29][30][31][32][33]35], the FM-QCP shows a substantial separation of scales between the onset of superconducting correlations and the phenomena described thus far, i.e non-Fermi liquid behavior and quantum-critical scaling, thus making this system an ideal platform for the investigation of quantum fluctuations close to criticality.  As described in Appendix E, the spin fluctuations induce an attractive interaction in the spin-triplet channel. The twoband structure of the model allows for a number of distinct superconducting order parameters, of which we find that the strongest pairing tendencies occur in the orbital-singlet, spintriplet channel, with the order parameter where σ =↑, ↓ is the spin index. Indeed, this channel is found to be the leading instability in a weak-coupling, mean-field analysis, see Appendix E. The two components ∆ ↑ and ∆ ↓ are related by the Z 2 (spin-flip) symmetry of the model, and are therefore of equal magnitude in the magnetically disordered phase. The order parameter transforms as a scalar under lattice rotations and reflections, and hence we expect single-fermion excitations to be fully gapped in the superconducting state.
Because fermions in our model only preserve a U(1) spin rotational symmetry, a finite-temperature triplet ordering is in principle allowed, in contrast to 2D systems with SU(2) spin rotational symmetry where triplet pairing may only occur at zero temperature. As shown in Appendix F, the two-component nature of the order parameter allows, in principle, for more exotic phases such as a charge 4e superconductor and a spin-nematic phase [44][45][46][47]. We have found no numerical evidence for the existence of these phases, and therefore focus on the triplet, charge 2e superconductor described above.
To probe for superconducting tendencies near the QCP, we explore three parameter sets, {ξ = 1.0, J = 1.0, µ = −0.5}, {ξ = 1.5, J = 0.5, µ = −2}, and {ξ = 3.0, J = 0.5, µ = −2}. First, the finite temperature FM to PM phase transitions are identified by finite-size scaling of the Ising spin susceptibilities, as discussed in Appendix D. Then, as depicted in Fig. 4, the FM-QCPs are located by extrapolating the finite temperature phase boundary towards T = 0. As expected, the larger the coupling ξ, the higher the critical field h c at the FM-QCP.
Next, we calculate the pairing correlations for all pairing channels, with on-site and nearest-neighbor form-factors (see Appendix G). Among all the available channels, only the order parameter defined in Eq. (11) has pairing correlations peaked at the QCP. We conclude that it is the only channel which is substantially enhanced by critical FM fluctuations. Fig. 5 shows the pairing structure factor C = At all couplings we find a peak of the pairing correlations at the QCP. At ξ = 1 and ξ = 1.5, the pairing correlations do not grow with the system size, indicating we are far from a superconducting transition. At the strongest coupling, ξ = 3, C grows, albeit very slowly, with L, indicating that the correlation length of the superconducting order parameter is moderately large.
Although pairing correlations in this channel are enhanced close to the QCP, long-range or quasi-long range superconducting order never develops down to T = 0.025 for any of the three parameter sets. This conclusion in corroborated by an analysis of the superfluid density, shown in Appendix G.

VII. DISCUSSION
In this work, we have constructed a lattice model that realizes a clean ferromagnetic QCP in an itinerant Fermi system. At the QCP, clear non-Fermi liquid behavior is observed. The static (ω n = 0) ferromagnetic susceptibility obeys scaling with an anomalous dimension η = 2 − a q = 0.15, which deviates from both the (2+1)d-Ising value (η = 0.036) and the mean-field value expected from Hertz-Millis theory (η = 0).
The ω n dependence of the bosonic susceptibility is highly nontrivial. In particular, the long-wave length limit (q → 0) and the low-frequency limit (ω n → 0) do not commute. To fully characterize the quantum dynamics at the QCP would requires to obtain detailed information about ∆(q, ω n ) in Eq. (10), especially in the low-energy limit. The presence of an anomalous dimensions in a q and in transition temperature T N suggest that ∆ probably also deviates from the Hertz-Millis form, which has a dynamic critical exponent z = 3 and thus predicts mean-field exponents. At the qualitative level, our conclusions and the observation of anomalous dimensions are in good agreement with the four-loop scaling analysis in Ref. [21]. However, although our studies reveal direct information on the qualitative features of ∆(q, ω n ) (i.e. the singular behavior in the limit q → 0, ω n → 0), the detailed functional form is beyond the resolution set by finite-size and finite-temperature effects of our QMC study.
It is interesting to contrast the results of the present study with those obtained for a related (but different) problem of an itinerant Ising-nematic QCP [30,35]. The two problems have essentially the same description within Hertz-Millis theory. Similarly to our results for an Ising ferromagnetic QCP, the Ising nematic QCP displays strong deviations from Fermi liquid behavior. A key difference, however, is that in the FM-QCP, Fermi-liquid behavior is lost everywhere along the Fermi surface, i.e there are no 'cold-spots' with long-lived quasi-particles. Furthermore, near the FM-QCP, strong deviations from Fermi-liquid behavior were found at temperatures much larger than any scale associated with superconductivity. In the Ising-nematic problem, the Ising order parameter correlations also show a singular behavior in the (q → 0, ω n → 0) limit. However, in the FM case presented here, the fermionic magnetic susceptibility is exactly conserved, unlike the Isingnematic case where an approximate conservation was found. An additional key difference is the anomalous exponents detected in the present work.
Perhaps the most striking difference between the two models lies in their superconducting properties. While the Isingnematic QCP was found to be strongly unstable towards s-wave superconductivity [35], the leading superconducting instability in the FM case is described by a two-component, spin triplet order parameter, with the possibility of exhibiting interesting, unconventional phases such as a charge 4e superconductor.
On the experimental side, our results may shed some new light on the study about 2d or quasi-2d metallic materials with Ising magnetic order, like Fe 1/4 TaS 2 [48] and CeCd 3 As 3 [49]. The determinantal quantum Monte Carlo (DQMC) formalism starts with the partition function of original Hamiltonian. To efficiently evaluate the trace in partition function, discretized imaginary time is used and β = M∆τ (∆τ = 0.05). As the studied model contains both fermion and Ising degrees of freedom, the trace will involve both a sum over Ising spin configurations and a determinant after tracing out the fermion degrees of freedom.
Let S = s z 1 · · · s z N denoting the Ising spins, then now we can trace out the fermion degrees of freedom, and obtain the configurational weight, with the Ising part where Λ 2 = sinh(∆τh) cosh(∆τh), γ = − 1 2 ln (tanh(∆τh)). For the fermion part, we have As an anti-unitary symmetry iτ y K (where τ y is a Pauli matrix in the orbital basis and K is the complex conjugation operator) make the Hamiltonian invariant, the fermion part ratio can be further rewritten as with K λσ the hopping matrix for orbital λ and spin σ. It turns out both the fermion weight and the Ising weight are always positive, thus there is no sign problem. To systematically improve the simulation, especially close to (quantum) critical point, we have implemented both local update in DQMC and space-time global update [34]. In the global update, we use Wolff algorithm [50] to propose space-time clusters of the Ising spins and then calculate the fermion weight to respect the detail balance as the acceptance rate of the update. Further attempts, with the recently developed self-learning determinantal quantum Monte Carlo scheme [51][52][53], which can greatly reduce the autocorrelation at (quantum) critical point and speedup the simulation with O(N) fold, be able to access larger L and lower T, are in progress.

Appendix B: z-direction flux
To reduce spurious finite size effects, we have used the techniques introduced in Ref. [54]. The basic idea is to introduce an effective z-direction flux by multiplying the hopping parameter by Peierls phase factors, i.e To make sure the model remains free of the sign problem, we take therefore the applied flux is not a true magnetic field, as it couples differently to fermions of different flavors. Translational invariance imposes restrictions on the magnitude of the effective magnetic field, namely, B λσ = n Φ 0 L 2 with n an integer and Φ 0 the flux quanta. We use the Landau gauge A λσ (r) = −B λσ yx in the bulk of the system, whereas special care is taken at the 'edges' While applying a flux in the z direction dramatically improves the convergence to the thermodynamic limit, it also breaks translation symmetry, making it impossible to extract the momentum dependence of fermionic correlations. Whenever such information is needed (such as for the momentum resolved, singlet-particle Green's function G k (ω n ), we apply a flux in the x or y direction, which is equivalent to twisting the boundary conditions. Just as in the case of the z directed flux, choosing (B1) ensures the absence of the sign-problem [55]. We have used six kinds of twisted boundary conditions, (0, 0), (π/2, 0), (π/2, π/2), (π, 0), (π, π/2) and (π, π), thus to get 16 times higher resolution of momentum.

Appendix C: Thermal phase transition
The thermal phase boundary in Fig. 1 (b) of the main text is controlled by the 2d Ising critical exponents γ = 7/4 and ν = 1, implying that the zero frequency and zero momentum Ising spin susceptibility around the finite temperature critical field h N satisfies Here, the Ising spin susceptibility is defined as   Fig. 6 (b) is the data collapse according to Eq. C1, from which we can obtain h N (T = 0.5) ≈ 3.06. Due to the coupling between the fermions and the Ising spins, the fermions go through the same finite temperature phase transition, as shown in Fig. 7 (a). The fermionic spin susceptibility χ F (h, T, 0, 0), where behaves much like χ(h, T, 0, 0), and a data collapse, shown in Fig. 7 (b) gives rise to the same h N (T = 0.5) ≈ 3.06. In the main text, we have discussed the dynamic Ising spin susceptibility, χ(h, T, q, ω n ), and performed the quantum critical scaling analysis. Here we reveal more details.
3. Triplet superconductor: both sectors have quasi-long range order. Here ∆ ↑ has power law correlations, The phase diagram of the same model (in different physical contexts) has been studied in Refs. [44][45][46]. To determine the phase diagram, we consider the criteria for stability against the appearance of a single vortex. In this model there are three kinds of vortices: 1. A vortex of one spin species and an antivortex in the other. Across the branch cut, θ c → θ c , θ s → θ s + 2π.

2.
A vortex of both the spin species. Across the branch cut, θ c → θ c + 2π, θ s → θ s .
The phase diagram can be derived from considering the free energy of a single unpaired vortex, given by F = E − T S. (As in the usual Berezinski-Kosterlitz-Thouless transition, such an analysis reproduces the phase diagram from a more rigorous renormalization group treatment.) The stability condition is F > 0. The energy of a vortex of type 3 is where L is the system size and a is some short range cutoff. Similarly E 2 = πK s log( L a ), and E 1 = πK c log( L a ). The entropy is the same in all cases, T S = log( L 2 a 2 ). The resulting phase diagram is given in Fig. 10. The stiffnesses K c and K s can be extracted from certain current-current correlation functions [56], Here, where and j x iλσ = ite iφ λσ i, i+x c † iλσ c i+x,λσ + H.c is the current density for fermions of orbital λ and spin σ.

Appendix G: Superfluid density and pairing correlations
In this appendix we provide further details on the numerical evidence for superconductivity.
To identify the leading pairing channel, we considered all possible on-site and nearest neighbor pairing order parameters and computed the pair structure factor for each channel. In all channels other than the orbital-singlet, spin-triplet channel defined in Eq. (11) of the main text, we find a weak response with no substantial system-size dependence or enhancement close to the FM QCP (not shown).
To test whether there are possible superconductivity instabilities close to the QCP, we measured the superfluid densities ρ c and ρ s which are related to the stiffnesses K c and K s defined in Eq. (F3). ρ c,s = K c,s /β. In Fig. 11 we show the temperature dependence of ρ c + ρ s . For both sets of parameters, ρ c + ρ s < 8 π T at the largest system size, and it decreases with the system size, implying the absence of quasi-long range superconducting order down to T = 0.025 (see Fig. 10). Note that upon increasing coupling strength and decreasing temperature, the finite-size estimates for the superfluid density grow. It is likely that there is a transition to a superconducting phase at higher coupling strengths or lower temperatures.  For completeness, in Fig. 12, 13 we show ρ c , ρ s , respectively.

Appendix H: Fermion spin susceptibility
In the Fig. 2(a) of the main text, we discussed the anisotropy of the quasiparticle fraction Z k F (T) along θ = 0 and θ = π 4 directions, and associate such anisotropy to the soft fermionbilinear mode around Q = (π, π) in the fermion spin susceptibility, besides the dominant Q = (0, 0) ferromagnetic fluctuations. Here we show the fermion spin susceptibility (Eq.(C3)) in Fig. 14, and reveal it is indeed the case. Fig. 14 shows the fermion spin susceptibility at h = h c , T = 0.05 from a L = 24 system. The strongest intensity is naturally at Q = (0, 0), however, near Q = (π, π) , there are soft modes (the signal is very weak and we have to fix the intensity range less than one and take a logarithm of original data).