A Replica Stoner Theory for Dirty Ferromagnets

TL;DR

This paper develops a theoretical framework combining replica and Stoner mean-field theories to analyze how disorder influences ferromagnetism, revealing disorder-enhanced ferromagnetism and the potential emergence of spin-glass phases.

Contribution

It introduces a novel combined replica-Stoner mean-field approach to study disorder effects on ferromagnetic metals, including phase diagram construction and disorder-induced ferromagnetism.

Findings

Disorder can enhance ferromagnetism in dirty ferromagnets.

A spin-glass phase may appear at high disorder levels.

Disorder can induce ferromagnetism in systems non-ferromagnetic in the clean limit.

Abstract

This paper investigates the effect of disorder on a ferromagnetic metal with repulsive interactions. We assume that, in the clean limit, the ferromagnetic state can be described by Stoner mean-field theory and study how disorder affects the the system by using a combined replica + Stoner mean-field approach. At zero temperature, we find that a replica-symmetric ferromagnetic mean-field solution exists in the presence of disorder with a modified Stoner criteria where the ferromagnetism is enhanced by disorder. At finite temperature, a Landau theory is employed to construct the phase diagram, revealing that beyond a critical disorder strength, a spin-glass phase may exist between the high-temperature paramagnetic phase and the low-temperature ferromagnetic phase. For weak (repulsive) interaction where the system is non-ferromagnetic in the clean limit, the possibility of a disordered-induced ferromagnetic ground state is observed both at zero temperature and finite temperature. The potential applicability of this framework to realistic materials is briefly discussed.

Equations

$$H_0=\sum _{ij,\sigma }{t}_{ij}{c}_{i\sigma }^{†}{c}_{j\sigma }+\sum _{i,\sigma }{W}_i{n}_{i\sigma },$$

$$H_{\mathrm{int}}=\sum _{i}U{c}_{i\uparrow }^{†}{c}_{i\uparrow }{c}_{i\downarrow }^{†}{c}_{i\downarrow }$$

$$\sum _{i}U{n}_{i\uparrow }{n}_{i\downarrow }=\sum _i-\frac{1}{3}U{\mathbf{S}}_i\cdot {\mathbf{S}}_i+\sum _i\frac{U}{4}{n}_i{n}_i,$$

$$\sum _{i}U{n}_{i\uparrow }{n}_{i\downarrow }\to \sum _i-\frac{1}{3}U{\mathbf{S}}_i\cdot {\mathbf{S}}_i.$$

$$H=\sum _{k\sigma }{\xi }_k{c}_{k\sigma }^{†}{c}_{k\sigma }-\sum _{kq,lp}{U}_{kq,lp}{\mathbf{S}}_{kq}\cdot {\mathbf{S}}_{lp},$$

$$U_{kq,lp}=\sum _{i}U{\varphi }_k^{\ast }({\mathbf{r}}_i){\varphi }_q({\mathbf{r}}_i){\varphi }_l^{\ast }({\mathbf{r}}_i){\varphi }_p({\mathbf{r}}_i).$$

$${\mathbf{S}}_{kq}\cdot {\mathbf{S}}_{lp}\to {\mathbf{S}}_{kq}\cdot ⟨{\mathbf{S}}_{lp}⟩+⟨{\mathbf{S}}_{kq}⟩\cdot {\mathbf{S}}_{lp}-⟨{\mathbf{S}}_{kq}⟩\cdot ⟨{\mathbf{S}}_{lp}⟩$$

$$\sim {\delta }_{lp}{\delta }_{kq}({\mathbf{S}}_k\cdot \text{expectationvalue}{\mathbf{S}}_p+\text{expectationvalue}{\mathbf{S}}_k\cdot {\mathbf{S}}_p-\text{expectationvalue}{\mathbf{S}}_k\cdot \text{expectationvalue}{\mathbf{S}}_p)$$

$$H_{MF}=\sum _{k\sigma }{\xi }_k{c}_{k\sigma }^{†}{c}_{k\sigma }-\sum _k\left(2{\mathbf{h}}_k\cdot {\mathbf{S}}_k-{\mathbf{h}}_k\cdot ⟨{\mathbf{S}}_k⟩\right),$$

$${\mathbf{S}}_i\cdot {\mathbf{S}}_i\sim 2⟨{\mathbf{S}}_i⟩\cdot {\mathbf{S}}_i-⟨{\mathbf{S}}_i⟩\cdot ⟨{\mathbf{S}}_i⟩.$$

$${\xi }_{k\alpha }^{\mathrm{HF}}{\varphi }_{k\alpha }^{\mathrm{HF}}({\mathbf{r}}_i)=\sum _{\beta ,j}[t_{ij}{\delta }_{\alpha \beta }+{\delta }_{ij}({\overline{W}}_i{\delta }_{\alpha \beta }-U\text{expectationvalue}{\mathbf{S}}_i\cdot {\mathbf{\sigma }}_{\alpha \beta }\left)\right]{\varphi }_{k\beta }^{\mathrm{HF}}({\mathbf{r}}_j).$$

$$\sum _i{\varphi }_{k\alpha }^{\ast }({\mathbf{r}}_i)U\text{expectationvalue}{\mathbf{S}}_i\cdot {\mathbf{\sigma }}_{\alpha \beta }{\varphi }_{k\beta }({\mathbf{r}}_i)\sim {\mathbf{h}}_a\cdot ⟨{\mathbf{\sigma }}_{\alpha \beta }{⟩}_a$$

$$\mathcal{F}[\mathbf{h},\mathbf{S}]=\sum _k{g}_k(h_k)-\sum _{k\ne p}{U}_{kp}{\mathbf{S}}_k\cdot {\mathbf{S}}_p+\sum _{k}2{\mathbf{h}}_k\cdot {\mathbf{S}}_k$$

$$g_k=-\frac{1}{\beta }\log$$

$$+e^{-\beta (2{\xi }_k+U_{kk})}).$$

$$\frac{\partial g_k}{\partial {\mathbf{h}}_k}=-2{\mathbf{S}}_k\text{ and }{\mathbf{h}}_k=\sum _{p(\ne k)}{U}_{kp}{\mathbf{S}}_p,$$

$$P(U_{kp})\sim \frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{(U_{kp}-\overline{U})}{2{\sigma }^2}}$$

$$\text{expectationvalue}{U_{kp}}_d\sim \frac{U}{L^d}\times P_{kp}=\frac{U}{V}$$

$$\text{expectationvalue}{U_{kp}^2}_d-\text{expectationvalue}{U_{kp}}_d^2\sim \frac{U^2}{L^{2d}}\times P_{kp}=\frac{U^2}{L^{d}V}\equiv \frac{h^2}{V}$$

$$\text{expectationvalue}{U_{kk}}_d\sim \frac{U}{L^d}$$

$$-\beta \text{expectationvalue}{\mathcal{F}}_d=\int \prod _{k\ne p}(\mathrm{d}{U}_{kp}P(U_{kp}))\log \mathcal{Z}[U_{kp}]=\text{expectationvalue}{\log \mathcal{Z}}_d$$

$$\text{expectationvalue}{\log \mathcal{Z}}_d=\underset{n\to 0}{\lim }\frac{\text{expectationvalue}{\mathcal{Z}}^n-1}{n}.$$

$$\text{expectationvalue}{f}_d=\frac{\text{expectationvalue}{\mathcal{F}}_d}{V}=3{\phi }_{+}{\phi }_{-}+\frac{{\mathbf{m}}^2}{4U}-\frac{1}{\beta V}\sum _k\int \frac{\mathrm{d}{\mathbf{y}}_k}{(2\pi {)}^{3/2}}{e}^{-\frac{{\mathbf{y}}_k^2}{2}}\log {\mathcal{Z}}_k[{\mathbf{y}}_k]$$

$${\mathcal{Z}}_k[{\mathbf{y}}_k]=\int \mathrm{d}{\mathbf{h}}_k\mathrm{d}{\mathbf{S}}_k{e}^{-\beta F_k[{\mathbf{y}}_k,{\mathbf{h}}_k,{\mathbf{S}}_k]}$$

$$F_k[{\mathbf{y}}_k,{\mathbf{h}}_k,{\mathbf{S}}_k]=g_k(h_k)-h{\phi }_{-}{\mathbf{S}}_k\cdot {\mathbf{S}}_k-(\mathbf{m}+\sqrt{2h\bar{\phi }}{\mathbf{y}}_k-2{\mathbf{h}}_k)\cdot {\mathbf{S}}_k.$$

$$M=\frac{m}{2U}=\frac{1}{V}\sum _k\text{expectationvalue}{\text{expectationvalue}{S}_k^z}_R$$

$${\phi }_{-}=\frac{h}{3\sqrt{2h\bar{\phi }}}\frac{1}{V}\sum _k\text{expectationvalue}{\text{expectationvalue}{\mathbf{y}}_k\cdot {\mathbf{S}}_k}_R$$

$${\phi }_{+}=\frac{h}{3}\frac{1}{V}\sum _k\text{expectationvalue}{\text{expectationvalue}{\mathbf{S}}_k\cdot {\mathbf{S}}_k}_R-\frac{h}{6\beta \sqrt{2h\bar{\phi }}}\frac{1}{V}\sum _k\text{expectationvalue}{\text{expectationvalue}{\mathbf{y}}_k\cdot {\mathbf{S}}_k}_R.$$

$$\text{expectationvalue}{\text{expectationvalue}\mathcal{O}[{\mathbf{S}}_k]}_R=\int \frac{\mathrm{d}{\mathbf{y}}_k}{(2\pi {)}^{3/2}}{e}^{-\frac{{\mathbf{y}}_k^2}{2}}(\frac{1}{{\mathcal{Z}}_k}\int \mathrm{d}{\mathbf{h}}_k\mathrm{d}{\mathbf{S}}_k\mathcal{O}[{\mathbf{S}}_k]e^{-\beta F_k})$$

$$N(\xi )\sim N(0)+\xi N^{\prime }(0)+\frac{1}{2}{\xi }^2{N}^{\prime \prime }(0).$$

$$M=M_0+\frac{h^2}{6U^2}(6a+5-2UN(0))$$

$$M$$

$$+(\frac{N^{\prime \prime }(0)}{6}+\frac{N^{\prime \prime }(0)}{2\sqrt{\frac{h^2}{6}N(0)U_{kk}}}\sqrt{\frac{2}{\pi }}{U}_{kk})U^2{M}^3.$$

$${\mathcal{F}}_{\mathrm{L}}$$

$$+\alpha (h{\phi }_{-})+\gamma (h{\phi }_{-}{)}^2+\eta (M^{2}h{\phi }_{-})+...$$

Full text

††preprint: APS/123-QED

I Introduction

The influence of disorder on magnetic systems has garnered significant interest within both experimental and theoretical physics communities [1, 2, 3, 4, 5, 6, 7]. In particular, the study of itinerant electron magnetism in the presence of disorder has been extensively explored through various approaches, including diagrammatic perturbation techniques [8, 9], mean-field theories [10], Monte-Carlo simulations [11], and self-consistent first-principle numerical calculations [12, 13, 14]. These studies have provided us understandings on the stability, phase transitions, and transport phenomena in disordered ferromagnets. Experimental investigations have also revealed intriguing properties of ferromagnetic materials subjected to different types of disorder [15, 16, 17, 18]. Notably, some experiments suggest that disorder can enhance ferromagnetism across a range of temperatures, while others indicate that strong disorder suppresses magnetic order [15].

In this work, we propose a semi-quantitative theory of ferromagnetism in a dirty metal. Specifically, we focus on systems where ferromagnetism can be described by Stoner mean-field theory [19] in the clean limit. For dirty metals, we combine the Stoner mean-field theory and replica-trick to study the effect of disorder on Stoner ferromagnetism. Our microscopic model considers electrons on a lattice with a random on-site potential and on-site repulsive interactions. Disorder effects are encoded in two order parameters within the mean-field equations derived from the replica-Stoner approach. We shall first consider the zero-temperature mean-field solutions in our theory and then construct the corresponding finite temperature Landau free energy to study the finite temperature magnetic phase transitions. Our results indicate that disorder can enhance ferromagnetism in ferromagnetic materials and even induce magnetism in systems which are non-magnetic in the clean limit. We also find that beyond a critical disorder strength, a spin-glass phase may exist between the high-temperature paramagnetic phase and the low-temperature ferromagnetic phase. These findings qualitatively align with recent experimental studies [16, 17, 18] and theoretical studies [8, 10, 14]. The framework presented here, combining mean-field theory and replica-trick, can be extended to other systems such as dirty superconductors [20] and spin-liquid states.

II Formulation

We begin with a Hamiltonian $\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{\rm{int}}$ describing interacting electrons on a lattice with random lattice potential and a Hubbard-type interaction:

[TABLE]

where $c^{\dagger}_{i\sigma}(c_{i\sigma})$ is the electron creation (annihilation) operator at site $i$. $W_{i}$ is a random on-site potential and $n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$. We focus on the case that the on-site interaction is repulsive ($U>0$) and the ferromagnetism is described by a Hartree mean-field theory which is essentially Stoner theory.

The mean-field theory is most conveniently performed by writing the interaction parts of the Hamiltonian $\mathcal{H}_{\rm{int}}$ as the sum of spin-spin interactions and the density-density interactions:

[TABLE]

in which $n_{i}=n_{i\uparrow}+n_{i\downarrow}$ and $\mathbf{S}_{i}=\frac{1}{2}c^{\dagger}_{i\alpha}\bm{\sigma}_{\alpha\beta}c_{i\beta}$ are the density and spin operator on site $i$, respectively, with $\bm{\sigma}_{\alpha\beta}=(\sigma^{x}_{\alpha\beta},\sigma^{y}_{\alpha\beta},\sigma^{z}_{\alpha\beta})$ representing the Pauli matrices. As we are interested only in magnetism, we assume that the density-density interaction gives rise only to a modified random potential $W_{i}\rightarrow\overline{W}_{i}=W_{i}+U\expectationvalue{n_{i}}/4$ which will be neglected in the following, i.e., we replace

[TABLE]

We shall absorb the factor $\frac{1}{3}$ into definition of $U$ in the following.

In the eigenstate basis given by $[c_{k\sigma},H_{0}-\mu N]=\xi_{k}c_{k\sigma}$ where $N$ is the total particle number operator and $c_{k\sigma}=\sum_{i}\phi_{k}(\mathbf{r}_{i})c_{i\sigma}$ where $\phi_{ki}$ are the eigenstates of $H_{0}$, (we set $\hbar=1$), the Hamiltonian can be written as (see also ref. [21])

[TABLE]

where $\xi_{k}$ is the single-particle energy of state $k$ measured from the chemical potential $\mu$, $\mathbf{S}_{kq}=\frac{1}{2}c^{\dagger}_{k\alpha}\bm{\sigma}_{\alpha\beta}c_{q\beta}$ and

[TABLE]

We note that $k$ represents eigenstates of $H_{0}$ which are not plane-wave states. We shall assume that $\xi_{k}$ forms a continuous spectrum across the Fermi surface.

II.1 Mean-Field Theory

We introduce the mean-field decoupling

[TABLE]

where $\mathbf{S}_{p}=\mathbf{S}_{pp}$ and $\langle...\rangle$ denotes thermal average. Denoting $U_{kk,pp}\equiv U_{kp}$, the corresponding mean-field Hamiltonian is

[TABLE]

where $\mathbf{h}_{k}=\sum_{p(\neq k)}U_{kp}\langle\mathbf{S}_{p}\rangle$. The corresponding eigen-energies are $\xi_{k\sigma}=\epsilon_{k}-\sigma|{\mathbf{h}_{k}}|$, where $\sigma=\pm 1$, denoting eigenstates $\phi_{k\sigma}(\mathbf{r}_{i})$ with spin parallel and anti-parallel to effective magnetic fields $\mathbf{h}_{k}$.

To understand our mean-field theory, we compare it with the usual Hartree theory based on the mean-field decoupling

[TABLE]

The corresponding eigenstates $\phi^{\mathrm{HF}}_{a}(\mathbf{r}_{i})$ with index $a=(k,\alpha=\pm 1)$ satisfies the mean-field equation

[TABLE]

To arrive at our mean-field theory we assume that disorder is strong enough so that all states $\phi^{\mathrm{HF}}_{a}(\mathbf{r}_{i})$ are localized. We also assume that $U\expectationvalue{\mathbf{S}_{i}}$ remains more-or-less constant within the region where the state $\phi^{\mathrm{HF}}_{a}(\mathbf{r}_{i})$ exists, and we approximate: $U\expectationvalue{\mathbf{S}_{i}}\sim{\mathbf{h}_{a}}$ where $\mathbf{h}_{a}$ represents mean value of $U\expectationvalue{\mathbf{S}_{i}}$ in this region. Notice we do not assume that ${\mathbf{h}_{a}}$ is the same for all states $\phi^{\mathrm{HF}}_{a}(\mathbf{r}_{i})$. Correspondingly

[TABLE]

where $\langle\bm{\sigma}_{\alpha\beta}\rangle_{a}=\sum_{i}\phi^{*}_{k\alpha}(\mathbf{r}_{i})\bm{\sigma}_{\alpha\beta}\phi_{k\beta}(\mathbf{r}_{i})$. With this approximation, we recover our mean-field theory given in Eq. (6) and Eq. (7) after some simple algebra.

The above approximation yields a clean ferromagnet when ${\mathbf{h}_{k}}\sim{\mathbf{h}}$ is constant as in disorder-free systems. For localized electronic states, we expect that fluctuation in ${\mathbf{h}_{k}}$ exists because of randomized coupling between localized states that are far apart. For states $k$ close to the Fermi surface with $|{\mathbf{h}_{k}}|>\xi_{k}$, the state $(k{\downarrow})$ becomes unoccupied resulting in formation of a local magnetic momentum polarized along the direction of ${\mathbf{h}_{k}}$.

It is convenient to introduce the free energy

[TABLE]

where

[TABLE]

We note that the mean-field decoupling excludes the self-interaction between the same $k$-eigenstate. The $U_{kk}$ term is treated exactly in evaluating the free energy of the doubly occupied $k$-state in $g_{k}$. The different treatment of $k\neq p$ and $k=p$ arises because they have very different statistical properties which we shall discuss in the next section.

Minimizing the free energy with respect to $\mathbf{h}_{k}$ and $\mathbf{S}_{k}$, we obtain

[TABLE]

corresponding to our mean-field theory (7). We shall see in the following that in the absence of disorder the interaction become constant $U_{kp}\rightarrow U_{0}$ and the mean-field equations reduce to the Stoner theory of ferromagnetism where a uniform solution $\mathbf{S}_{k}\equiv\mathbf{S}_{0}$ is obtained when $U_{0}$ satisfies the Stoner Criteria $U_{0}N(0)\geq 1$ [22], $N(0)$ is the density of states at the Fermi Surface. In the presence of disorder, $U_{kp}$ is no longer a constant and the solution $\mathbf{S}_{k}$’s are in general different for different $k$’s. For strong enough disorder, the average of $\mathbf{S}_{k}$ might vanish, resulting in a spin glass state as in classical spin models with random interaction between spins [23, 24].

To study these possibilities, we have to compute the disorder-averaged free energy $\expectationvalue{\mathcal{F}}_{d}$ as well as other physical observable $\expectationvalue{\mathcal{O}}_{d}$, where $\expectationvalue{...}_{d}$ denotes the average over disorder configurations. We shall employ the replica trick [24] to compute the disorder averages in this paper.

II.2 The Replica Approach

We discuss how the replica trick is implemented in this section. We note that with our approximate free energy Eq. (11) the disorder effects appears only through the interaction matrix elements $U_{kp}$ and $U_{kk}$.

To perform the replica trick, we assume that $U_{kp},U_{k^{\prime}p^{\prime}}$ are mutually uncorrelated for $k(p)\neq k^{\prime}(p^{\prime})$ and follows a Gaussian distribution:

[TABLE]

where $\overline{U}$ is the mean and $\sigma^{2}$ the variance. We caution that this is a very crude approximation as all the effects of disorder is absorbed in only two parameters in our approximation. In particular, the distribution is $(k,p)$-independent and spatial information is lost, corresponding to infinite range interaction in $k$-space.

The parameters $U$ and $\sigma^{2}$ have been estimated for localized single particle wavefunction in $d$-dimension space where $|\phi_{k}(\mathbf{r}_{i})|^{2}\sim\frac{1}{L^{d}}e^{-\frac{|\mathbf{r}_{k}-\mathbf{r}_{i}|}{L}}$, where $L$ is the localization length. The details of the calculation can be found in ref. [25]. It was found that the mean and variance of $U_{kp}(k\neq p)$ satisfy:

[TABLE]

where $V$ is the system volume, and $P_{kp}\sim\frac{L^{d}}{V}$ signifies the probability of finding two localized states $(k,p)$ within a distance $\sim L$ to each other. We define $h^{2}=\sigma^{2}V=U^{2}/L^{d}$ for later convenience. Please note that both the mean and variance of the distribution are of order $\mathcal{O}(1/V)$, which is necessary for obtaining finite result in thermodynamic limit [26]. In particular, we note that for extended states where $L^{d}\rightarrow V$, the effect of disorder vanishes in our approach, suggesting that a more refined approach beyond our crude approximation for $P(U_{kp})$ is needed in this case.

For $k=p$, we obtain

[TABLE]

as $P_{kk}=1$. It’s important to realize that $\expectationvalue{U_{kk}}_{d}\sim\mathcal{O}(1)$ whereas $\expectationvalue{U_{k\neq p}}_{d}\sim\mathcal{O}(1/V)$. This important distinction suggests that $U_{kp}$ and $U_{kk}$ should be treated differently.We have treated $U_{kk}$ exactly in writing down $g_{k}$. We now employ the replica trick to treat the fluctuations in $U_{kp}(k\neq p)$ in the following.

The disorder averaged free energy is given by

[TABLE]

where $\mathcal{Z}[U_{kp}]$ is the partition function for a given configuration of $\{U_{kp}\}$. The replica trick computes $\expectationvalue{\log\mathcal{Z}}_{d}$ by using the identity

[TABLE]

For integer $n$, the integral over $dU_{kp})$ can be evaluated easily with the free energy (11) and Eq. (14) for $P(U_{kp})$. As the distribution $P(U_{kp})$ is $(k,p)$-independent, the interaction term $\sum_{kp}U_{kp}\mathbf{S}_{k}\cdot\mathbf{S}_{p}$ is effectively infinite-range and a replica-mean-field theory can be employed in calculating $\expectationvalue{\mathcal{Z}^{n}}_{d}$ [24]. We consider only the replica-symmetric solutions in this paper. The limit $n\rightarrow 0$ is taken afterward through an analytical continuation. The resulting mean-field equation is presented in the following with the mathematical details given in Appendix A. At zero temperature, we find that the replica symmetric (RS) solution is stable against replica symmetry breaking (RSB) terms. The detail of the stability analysis is presented in Appendix D.

After some lengthy algebra, we obtain the disorder-averaged free energy

[TABLE]

where $\bar{\varphi}=\varphi_{+}-\frac{1}{2\beta}\varphi_{-}$. The free energy is characterized by three mean-field parameters $\mathbf{m}$, $\varphi_{+}$ and $\varphi_{-}$ which are determined by the mean-field equations

[TABLE]

Here $M$ is the average spin polarization characterizes average magnetization ($g\mu_{B}M$) and we have chosen a particular direction for magnetization $\mathbf{M}=M\hat{z}$ in writing down the equation for $M$. $\varphi_{-}\sim$ Edward-Anderson order parameter $q_{EA}$ [23] which characterize spin-glass order while the first term in $\varphi_{+}$ characterizes the average spin amplitude. The bracket $\expectationvalue{\expectationvalue{\mathcal{O}[\mathbf{S}_{k}]}}_{R}$ is a shorthand of:

[TABLE]

meaning that the configuration average over the $k$-state $\expectationvalue{\mathcal{O}[\mathbf{S}_{k}]}$ is performed first, followed by averaging over $\mathbf{y}_{k}$.

$\mathcal{Z}_{k}$ and $F_{k}$ can be regarded as effective partition function and free energy for a given $k$-state. The replica-symmetric mean-field theory leads to the replacement of coupling between different $k$-states in Eq. (11) by a random effective magnetic field $\mathbf{B}_{k}=\mathbf{m}+\sqrt{2h\bar{\varphi}}\mathbf{y}_{k}-2\mathbf{h}_{k}$ that couples to individual spin $\mathbf{S}_{k}$ and the disorder average is replaced by averaging over the Gaussian random variable $\mathbf{y}_{k}=(y_{k}^{x},y_{k}^{y},y_{k}^{z})$ that appears in $\mathbf{B}_{k}$.

At zero-temperature $\beta\rightarrow\infty$, the configuration average on state $k$ can be calculated by replacing $\expectationvalue{\mathcal{O}[\mathbf{S}_{k}]}$ with its saddle-point value $\mathcal{O}[\mathbf{S}_{k}^{(0)}]$ in which $\mathbf{S}_{k}^{(0)}$ minimizes the effective free energy $F_{k}$. This procedure recovers our mean-field result discussed above. A Landau theory is also developed from the free energy Eq. (20) to study the finite temperature phase transition. The mathematical details are shown in Appendices  A and  C. The three self-consistent equations at zero-temperature are solved numerically and the result is presented in sec. III. The corresponding finite temperature behaviors are presented in Sec. IV.

III Mean-Field results at zero temperature

The numerical solutions of mean-field parameters at zero-temperature are presented in this section. We have considered the mean-field phase diagram in two dimensions $(d=2)$ for various values of $L,U$. For weak interaction, only states at the vicinity of Fermi surface are important. We assume a continues, structureless Fermi surface where the only important information enters our analysis is the electronic density of states and have expanded the density of states around the Fermi surface as

[TABLE]

According to Eq. (15)-(17), we have set $U_{kk}\sim a\frac{h^{2}}{U}$ in solving the mean-field equations in which $a\sim\mathcal{O}(1)$ is a positive factor and we have set $a=1$ in the following numerical calculation, without loss of generality. We use a typical data of density of states at Fermi surface $N(0)=1.56\ \rm{eV}^{-1}$ per electron spin from a monolayer Stoner ferromagnets [27] and use a normalized Gaussian function to approximate its density of states and obtained the approximate data $N^{\prime}(0)=-2.84\ \rm{eV}^{-2}$, $N^{\prime\prime}(0)=-23.23\ \rm{eV}^{-3}$ per electron spin. We note that a negative $N^{\prime\prime}(0)$ is required for Stoner magnets [28] and using other reasonable choices of density of states data affects the numerical results but does not affect our qualitative conclusions.

We present in Fig. 1 the mean-field parameters $M,\varphi_{-}$ and $\varphi_{+}$ as a function of localization length $L$ (in unit of lattice spacing) for $UN(0)=1.2$. Denote $M_{0}$ as the magnetization in clean limit, we show $m,h\varphi_{-}$ and $\sqrt{h\varphi_{+}}$ in unit of $m_{0}=2UM_{0}$ in this plot, keeping them in a same dimension.

We observe that the spin glass order parameter and spin amplitude parameter are both zero while magnetization are nonzero from weak localization regime and increase as the disorder strength increases. The coexistence of ferromagnetism and spin glass order (all three parameters are nonzero) is similar to classical spin glasses [24, 29, 26] and represents a disordered ferromagnet. The enhancement in magnetization when the disorder strength increases is in agreement with previous experimental and theoretical results [18, 8].

In the limit $L^{-1}\rightarrow 0$, the self-consistent equations Eq. (21)-(23) can be studied analytically by expansion in powers of $h=U/L$ (2D). We found that in the leading order at $UN(0)>1$:

[TABLE]

in which $M_{0}=\sqrt{\frac{6(UN(0)-1)}{-U^{3}N^{\prime\prime}(0)}}$ is the magnetization in the clean limit $h=0$. It shows that the magnetization is enhanced by disorder for moderate strength of interaction $UN(0)$. We caution, however that our result is not reliable in the regime $UN(0)\ll 1$, as the precise electronic band structure of the system will become important in determining $N(\xi)$, and our expansion around the Fermi surface (Eq. (25)) would not be valid.

The mean-field equations in the regime $UN(0)<1$ are also studied where several solutions of magnetization (exact values of $M$) around $UN(0)\sim 1$ are plotted in the inset of Fig. 1. We find that magnetization can be induced by disorder even for $UN(0)<1$ as $L^{-1}$ increases.

The regime $UN(0)\ll 1$ is interesting. Numerical solutions of the mean-field equations are shown in Fig. 2 in which the exact values of magnetization and spin-glass order parameter are plotted as a function of the $UN(0)$ for different values of localization length.

We see that ferromagnetism develops at $UN(0)\ll 1$ and increases with increasing disorder. From Fig. 2(a), we see that magnetization is developed at $UN(0)\geq 0.28$ while from Fig. 2(b), we see that the spin-glass order persists throughout the whole range of $UN(0)$ and a cusp appears in the solution for $\varphi_{-}$ at the phase transition point $UN(0)\sim 0.28$.

The behavior of magnetism at small $UN(0)\ll 1$ can be studied analytically in the weak localization limit $h^{2}\rightarrow 0$. Using $U_{kk}\sim a\frac{U}{L^{d}}\sim a\frac{h^{2}}{U}$, we see that $U_{kk}\gg h^{2}N(0)$ in this limit suggesting we may solve the mean-field equation first in the limit $U_{kp}=0\ (k\neq p)$ while keeping $U_{kk}$ finite and then look at the effect of $U_{kp}$ perturbatively. At the starting limit, $U_{kp}=0$, it is easy to show that all states $k$ satisfying $\xi_{k}<U_{kk}/2$ are fully polarized but the directions of polarization are random when $U_{kp}=0$, i.e., a random magnetic moment state [25] is obtained in this limit. To leading order in $U_{kp}$, the equation for magnetization is (see Appendix B)

[TABLE]

We have assumed $N^{\prime\prime}(0)<0$, which is the case for usual ferromagnetic materials. We note that $\sqrt{\frac{16aUN(0)}{3\pi}}\ll 1$ for small $UN(0)\ll 1$ and magnetization $M=0$ in this limit. The system remains in a spin-glass state and magnetization develops only at $aUN(0)\geq 0.59$, in qualitative agreement with what we observed in our numerical solutions.

IV Finite temperature Landau Theory

Phase transitions between paramagnetic and ferromagnetic and/or spin glass states are expected to take place at non-zero temperatures. In this section we derive the Landau theory from the disorder averaged free energy at temperature near the phase transitions and study various plausible phase transition between the three states as a function of disorder and interaction. We see from Eq. (20a)-(20b) that $\expectationvalue{\mathcal{F}}_{d}$ is a functional of $\mathbf{h}_{k}$, $\mathbf{S}_{k}$ as well as $M$, $h\varphi_{-}$ and $h\varphi_{+}$. To derive the Landau theory, We have to eliminate $\mathbf{h}_{k}$, $\mathbf{S}_{k}$ and $h\varphi_{+}$ using the mean-field equations, and expand $\expectationvalue{\mathcal{F}}_{d}$ in powers of order parameters $M$ and $h\varphi_{-}$, noting that these parameters are small near phase transitions. We note that from Eq. (23) that $h\varphi_{+}$ is composed of $h\varphi_{-}$ and spin magnitude parameter $\sum_{k}\expectationvalue{\expectationvalue{\mathbf{S}_{k}\cdot\mathbf{S}_{k}}}_{R}$, which is not an order parameter representing a separate phase. The mathematical details of this straightforward but tedious calculation is presented in Appendix. C. The resulting Landau free energy is of the form [30]

[TABLE]

in which $h_{ext}$ is the external magnetic field. Factors $a$, $b$, $c$ and $\alpha$, $\gamma$, $\eta$ are functions of temperature and disorder strength $h$, expressed as:

[TABLE]

where the coefficients $a_{k}$ and $b_{k}$ are coming from the expansion (see Appendix. C for detailed expressions)

[TABLE]

Minimizing the Landau free energy Eq. (28) with respect to the order parameters $M$ and $h\varphi_{-}$, we obtain

[TABLE]

We observe that at the limit $h\rightarrow 0$ the magnetic susceptibility recovers the behaviors of pure Stoner magnets (see also Appendix. C). The qualitative behavior of the system in the presence of $h\varphi_{-}\neq 0$ can be understood by assuming that the coefficients $b(>0)$, $c(>0)$, $\gamma(>0)$, and $\eta(<0)$ are temperature independent, with $a(T)=a_{0}(T-T_{c})$ and $\alpha(T)=\alpha_{0}(T-T_{f})$, where $T_{c}$ and $T_{f}$ are the ferromagnet and spin-glass transition temperatures, respectively, assuming that $T_{f}$ and $T_{c}$ are close to each other.

There are two possible scenarios with this approximation, $T_{f}>(<)T_{c}$. In the first case $T_{f}>T_{c}$, a spin-glass transition is observed at $T_{f}$, with both $h\varphi_{-}=0$ and $M=0$ at $T>T_{f}$, and $h\varphi_{-}\neq 0$ and $M=0$ at $T_{f}>T>\overline{T}_{c}$, where $\overline{T}_{c}$ is determined by $a+\eta(h\varphi_{-})=0$, or

[TABLE]

as $\eta<0$, i.e., the ferromagnetic transition temperature is * enhanced* by he spin-glass order. An analogous conclusion can also be reached for $T_{c}>T_{f}$. Notice that the magnetic susceptibility exhibits a cusp at temperature $T_{f}(>T_{c})$ according to Eq. (31c) because of the appearance of nonzero $h\varphi_{-}$ at $T<T_{f}$.

To confirm the above results, we have computed the coefficients $a,b,c,\alpha,\gamma,\eta$ numerically to determine the phase boundaries at finite temperature for various values of $h^{2}$ and $UN(0)$. We observe that in agreement with the results at zero temperature, at finite temperature the ferromagnetism can also appear at $UN(0)\ll 1$.

The phase diagram at $UN(0)=0.32$ is shown in Fig. 3(a). The solid red line and dashed blue line represent the spin-glass and ferromagnet transitions, respectively. Both the spin-glass and ferromagnet transition temperatures are found to increase with disorder, with the ferromagnet transition temperature $T_{c}$ higher than the spin-glass transition temperature $T_{f}$ at weaker disorder (longer localization length) and with the order of transition reversed at stronger disorder (shorter localization length), in agreement with the general expectation that stronger disorder favors spin-glass state over ferromagnet. We note that the splitting between the spin-glass and ferromagnet transition temperatures is rather small ($\sim\mathcal{O}(0.1K$)) at the regime $L\gg 1$ where our approximation for the density of states Eq. (25) is valid. The magnetic susceptibly $\chi(T)$ at a strong disorder $L=7.0$ is shown in Fig. 3(b), the divergence at $T\rightarrow T_{c}$ and the cusp at $T\rightarrow T_{f}$ is clear.

Further numerical calculations suggests that for small localization length (large $L^{-1}$ for strong disorder strength), a non-magnetic material ($UN(0)\ll 1$) can exhibit net ferromagnetism at low temperature, and for large localization length (small disorder strength) the ferromagnetic phase boundary is close to the clean limit results. In both cases of non-magnetic materials and ferromagnetic materials the Curie temperature can be strongly enhanced as the disorder strength increases.

V Discussion

In this paper, we construct a replica-symmetric Stoner (RS-Stoner) mean-field theory to study how magnetization evolves in the presence of disorder in Stoner ferromagnet. To achieve this goal two crucial approximations are made: (i) We replace the usual on-site Hartree mean-field theory by a simplified mean-field theory where we do not have to determine the electronic eigenstates in the presence of disorder self-consistently. With this approximation the effects of disorder appear only in the effective interaction matrix element $U_{kp}$’s. (ii) To implement the replica trick, we make a very simple and crude approximation for the distribution of matrix elements $P(U_{kp})$ where we assume that $P(U_{kp})$ obeys a Gaussian distribution that depends only on two parameters. We also neglect the detail band structure of the material and assume that it is sufficient to expand the density of states around the Fermi surface Eq. (25) in our mean-field theory analysis. With these approximations a simple replica-symmetric Stoner mean-field theory can be constructed which can be analyzed rather easily. The aim of our paper is to see whether our calculation can qualitatively match the results from recent theoretical studies [8, 9, 10, 11, 12, 13, 14] and to see whether new insights can be obtained.

We find that at zero-temperature the net magnetization $M$ is enhanced by disorder for moderate interaction strength $U$, including the situation $UN(0)<1$ which is non-magnetic in clean case, in agreement with recent theoretical work [8]. For very weak interaction $UN(0)\ll 1$ We predict a spin glass state can be found. However, we caution that interaction between far away spins is neglected in our simple approximation of $P(U_{kp})$, and a ferromagnetic state may still form at exponentially small temperature [25].

Our constructed Landau theory at finite temperature predicts that the Curie temperature is enhanced by disorder, consistent with results from other recent studies [10, 14]. It also predicts that ferromagnetism can exist for $UN(0)\ll 1$, in agreement with our zero temperature calculation and multiple experimental results [16, 17, 18]. in addition, it suggests that a narrow region of spin-glass state can exists between the paramagnetic high temperature state and ferromagnetic low temperature when disorder is strong enough. The existence of spin-glass state shows up as a discontinuity in the derivative of magnetic susceptibility against temperature at temperature above the ferromagnetic transition.

Given the many approximations made in our theory, our results can at best be trusted only qualitatively. In particular, We expect that our Landau theory constructed from replica-symmetric solution is quantitatively not correct in the spin-glass state at finite temperature where replica-breaking terms become important, although it is expected that it can gives the phase boundaries correctly [24, 29]. We note that the replica-symmetric solution considered in our theory is stable at zero-temperature, suggesting that the finite temperature replica-symmetric mean-field theory may still be qualitatively correct.

Acknowledgements.

The project is supported by the School of Science, the Hong Kong University of Science and Technology.

Appendix A Details of the Replica mean-field Theory at Zero Temperature

We show here mathematical details of the replica mean-field theory at zero temperature. We begin with calculating $\expectationvalue{\mathcal{Z}^{n}}_{d}$,

[TABLE]

where $i$ is the replica index, $\mathop{}\!\mathcal{D}\mathbf{S}_{i}\mathop{}\!\mathcal{D}\bm{h}_{i}=\prod_{k}\mathop{}\!\mathrm{d}\mathbf{S}_{ik}\mathop{}\!\mathrm{d}\mathbf{h}_{ik}$ and $P(U_{kp})=\sqrt{\frac{V}{2\pi h^{2}}}e^{-\frac{V(U_{kp}-U/V)}{2h^{2}}}$ where $V$ is the volume of the system. Integrate over $U_{kp}$ which are Gaussian integrals, we obtain:

[TABLE]

where $\omega_{kp}=\sum_{i}\Big(S^{x}_{ik}S^{x}_{ip}+S^{y}_{ik}S^{y}_{ip}+S^{z}_{ik}S^{z}_{ip}\Big)$. We write

[TABLE]

and use the Hubbard-Stratonovich (HS) transformation to decouple $(\gamma^{\alpha\beta}_{ij})^{2}$, $\kappa^{\alpha}_{i}\kappa^{\alpha}_{i}$ and $\kappa^{\alpha\beta}_{k}\kappa^{\alpha\beta}_{k}$. Introducing the HS fields $q_{ij}^{\alpha\beta}$, $\chi^{\alpha}_{i}$ and $\chi^{\alpha\beta}_{i}$, we obtain:

[TABLE]

We shall employ a saddle point approximation for the parameters $q_{ij}^{\alpha\beta}$, $\chi^{\alpha}_{i}$ and $\chi^{\alpha\beta}_{i}$ in the following as the system has infinite-range interactions in $k$-space. Assuming unbroken spin rotational symmetry, i.e., $q^{\alpha\beta}_{ij}\rightarrow q_{ij}\delta_{\alpha\beta}$, $\chi^{\alpha}_{i}\rightarrow\chi_{i}$ and $\chi^{\alpha\beta}_{ij}(\alpha\neq\beta)\rightarrow 0$. In this case, the n-replica partition function is

[TABLE]

Next we consider the term $\frac{\beta U}{V}\sum_{kp}\mathbf{S}_{ik}\cdot\mathbf{S}_{ip}\equiv\frac{\beta U}{V}(\sum_{k}\mathbf{S}_{ik})^{2}$ where we perform another HS transformation:

[TABLE]

where $\mathop{}\!\mathrm{d}\mathbf{m}=\mathop{}\!\mathrm{d}m_{x}\mathop{}\!\mathrm{d}m_{y}\mathop{}\!\mathrm{d}m_{z}$ and $\mathbf{m}_{i}$ plays a role as an effective magnetic field. Treating $\mathbf{m}_{i}$ by another saddle-point approximation, which corresponds to a Hartree mean-field theory. We obtain:

[TABLE]

Assuming replica-symmetric solutions: $q_{ij}=q$, $\chi_{i}=\chi_{0}$, $\mathbf{m}_{i}=\mathbf{m}_{0}$, the non-diagonal term $\sum_{i\neq j}q_{ij}\sum_{k}\mathbf{S}_{ik}\cdot\mathbf{S}_{jk}$ can eliminated by another HS transformation

[TABLE]

Gathering all these HS fields, we obtain the n-replica symmetric partition function:

[TABLE]

where we define

[TABLE]

From Eq. (41), we can expand $\expectationvalue{\mathcal{Z}^{n}}_{d}$ w.r.t. $n$ as:

[TABLE]

and the disorder averaged free energy density is obtained by applying the identity of replica trick:

[TABLE]

Defining $\chi-q=\varphi_{-}$, $\chi+q=2\beta\varphi_{+}$, $\bar{\varphi}=\varphi_{+}-\frac{1}{2\beta}\varphi_{-}$ and $\mathbf{m}=\frac{\mathbf{m}_{0}}{\beta}$ we obtain

[TABLE]

The self-consistent equations for $\varphi_{+}$, $\varphi_{-}$ and $\mathbf{m}$ are obtained by minimizing the free energy:

[TABLE]

in which the average $\expectationvalue{\expectationvalue{...}}_{R}$ is defined as

[TABLE]

We next derive the mean-field equations at the zero temperature. For simplicity we have assume that $U_{kk}=constant$ from now. We see that the interaction term $U_{kk}$ appears only in doubly occupied states in $g_{k}(h_{k})$ and it is convenient to shift the single-particle energy $\xi_{k}\rightarrow\xi_{k}-U_{kk}/2$, obtaining at zero temperature:

[TABLE]

Moreover, $\bar{\varphi}\rightarrow\varphi_{+}$ at zero temperature and

[TABLE]

At zero temperature we can calculate the average $\expectationvalue{\expectationvalue{\mathcal{O}[\mathbf{S}_{k}]}}_{R}$ by replacing $\expectationvalue{\mathcal{O}[\mathbf{S}_{k}]}$ with the saddle-point value $\mathcal{O}[\mathbf{S}_{k}^{(s)}]$ where $\mathbf{S}_{k}^{(s)}$ is given by minimizing the effective free energy $F_{k}$. Minimizing it w.r.t. $\mathbf{h}_{k}$ and $\mathbf{S}_{k}$ we obtain

[TABLE]

where $\vec{e}_{k}$ is the direction vector of the vector $\mathbf{Y}_{k}=(\mathbf{m}+\sqrt{2h\varphi_{+}}\mathbf{y}_{k})$. We note that $\mathbf{S}^{(s)}_{k}\neq 0$ in the region $|\xi_{k}|<U_{kk}/2$ even when $\bm{h}^{(s)}_{k}=0$, i.e., spin are polarized in this region even in the absence of magnetic field, which is a result first pointed out in ref. [25].

As long as $\mathbf{h}_{k}^{(s)}\neq 0$, the saddle-point solutions show that $\mathbf{S}_{k}^{(s)}$ and $\mathbf{h}_{k}^{(s)}$ are pointing at the same direction for all states $k$ ($\mathbf{Y}_{k}\cdot\mathbf{S}_{k}=Y_{k}S_{k}\cos\theta$ is minimized when $\cos\theta=1$) and $S^{(s)}_{k}$ takes value $\frac{1}{2}$ for $-|\bm{h}^{(s)}_{k}|-U_{kk}/2<\xi_{k}<|\mathbf{h}^{(s)}_{k}|+U_{kk}/2$. Choose $\hat{z}$ as the direction of magnetization $\mathbf{m}=m\hat{z}$, we obtain for $m$:

[TABLE]

Using Eq. (25)) for the density of states, we obtain

[TABLE]

In terms of $M=\frac{m}{2U}$ which is the physical magnetization, this equation reduces to the equation from Stoner theory of ferromagnetism at $h=0$, $M=N(0)UM+\frac{N^{\prime\prime}(0)}{6}(UM)^{3}$. The Equations for $\varphi_{+}$ and $\varphi_{-}$ can be considered similarly, we obtain:

[TABLE]

[TABLE]

The above equations are solved numerically to obtain the mean-field parameters $M$ and $\varphi_{\pm}$.

Appendix B Solution of mean-field equation at $UN(0)\ll 1$ at zero temperature

In this section we analyze the mean-field equations (51)-(53) at the limit $UN(0)\ll 1$. Since $U_{kk}\sim a\frac{h^{2}}{U}$, this limit is equivalent to $U_{kk}\gg h^{2}N(0)$. Thus we may first solve the mean-field equation at $U_{k\neq p}=0$ with $U_{k\neq p}$ being treated as perturbation. At zeroth order, $U_{k\neq p}=0$ and $\mathbf{m}=0$ by definition. We obtain the directional vector of saddle-point solution $\mathbf{S}_{k}^{(s)}$ as

[TABLE]

In this situation, spins are fully polarized in the region $[-U_{kk}/2,U_{kk}/2]$ with fully random directions, i.e., a random magnetic moment state. Using the mean-field equation for $\varphi_{+}$, we obtain

[TABLE]

which is proportional to $h^{2}N(0)$. To first order in $U_{k\neq p}$, magnetization and spin glass order becomes nonzero. Denoting $\sqrt{2h\varphi_{+}}\sim\sqrt{\frac{h^{2}}{6}N(0)U_{kk}}\equiv C$, the directional vector of $\mathbf{S}^{(s)}_{k}$ can be expressed as:

[TABLE]

to first order of correction due to $\mathbf{m}$ (higher order corrections are not necessary because mean field equation for $m$ is naturally cubic). Similarly, we have

[TABLE]

Using (51) and (53) for $m$ and $\varphi_{-}$, and replacing the directional vector $\vec{e}_{k}$ by (56) and $|\mathbf{m}+\sqrt{2h\varphi_{+}}|$ by (57), we obtain after performing the Gaussian integral over $\mathbf{y}_{k}$ the spin glass order parameter to first order in $h^{2}N(0)$:

[TABLE]

and the magnetization $m$,

[TABLE]

where we have kept only lowest order term in $UN(0)$. Take $m=2UM$, we have the equation for physical magnetization $M$ at $UN(0)\ll 1$:

[TABLE]

Notice that $\sqrt{\frac{16aUN(0)}{3\pi}}-1<0$ in the limit $UN(0)\ll 1$ (for $a\sim\mathcal{O}(1)$), $M=0$ in this limit.

Appendix C Details of Obtaining the Landau Free Energy

We begin with the disorder averaged free energy, see Eq. (20a), we first expand $g_{k}(h_{k})$ in powers of $h_{k}$, assuming that $h_{k}$ is small near the phase transitions:

[TABLE]

where $g_{0}=-\frac{1}{\beta}\log(1+3e^{-\beta\xi_{k}}+e^{-\beta(2\xi_{k}+U_{kk})})$ and

[TABLE]

We then apply the mean-field equation $\frac{\mathop{}\!\partial g_{k}}{\mathop{}\!\partial\mathbf{h}_{k}}=-2\mathbf{S}_{k}$ to express $\mathbf{h}_{k}$ in terms of $\mathbf{S}_{k}$:

[TABLE]

With the expansion, we can eliminate $\mathbf{h}_{k}$ in the free energy, obtaining

[TABLE]

We then apply the saddle-point approximation:

[TABLE]

where $\mathop{}\!\overline{\mathbf{S}}_{k}$ represents the saddle-point solution. The corresponding saddle-point free energy is

[TABLE]

To obtain the Landau free energy near phase transitions, we assume that $h\varphi_{-}$ is small near the phase transitions, and expand:

[TABLE]

On the other hand, using mean-field equation for $h\varphi_{+}$ Eq. (23) we also express it as powers of $h\varphi_{-}$ and $\mathbf{m}$:

[TABLE]

Substituting these expansions into the saddle-point solution $\mathop{}\!\overline{\mathbf{S}}_{k}$ Eq. (66) and performing the integral over $\mathbf{y}_{k}$, we obtain the Landau free energy:

[TABLE]

where we have substitute $M=\frac{m}{2U}$. Please note that this Landau free energy is only valid in the vicinity of phase transitions and contributions of factors are approximately linear in temperature ($\sim T-T_{f}$ and $\sim T-T_{c}$). Other contributions such as $\mathcal{O}(h)^{4}$ in $h\varphi_{-}$ are neglected in this expansion. The factors of $\sum_{k}a_{k}^{n}$ can be evaluated numerically by energy integral

[TABLE]

If an infinitesimal uniform magnetic field term $\sum_{i}\mathbf{h}_{ext}\cdot c^{\dagger}_{i\alpha}\bm{\sigma}_{\alpha\beta}c_{i\beta}$ is added to microscopic Hamiltonian Eq. (1), the magnetic field will appear in $\mathbf{h}_{k}\rightarrow\mathbf{h}_{k}+\mathbf{h}_{ext}$ resulting in a linear term $-\sum_{k}Ua_{k}Mh_{ext}^{z}$ in the Landau free energy. We can see the longitudinal magnetic susceptibility can be expressed as:

[TABLE]

In the limit of $h=0$, the Landau free energy and magnetic susceptibility are reduced to the original form for Stoner ferromagnets.

Appendix D Stability of Replica-symmetric solution at Zero Temperature

We now show that at zero temperature the replica-symmetric solution is stable according to a replica-symmetry breaking (RSB) term. We shall write $q_{ij}=q+\delta q_{ij}$ in which $\delta q_{ij}$ is the small RSB term and see the related terms in the free energy. We begin with the n-replica partition function, Eq. (41), write as:

[TABLE]

where Trace denotes summing over all configurations of $\{\mathbf{S}_{ik}\}$ and $\{\mathbf{h}_{ik}\}$ as well as the disorder average over $\mathbf{y}_{k}$. Write $q+\delta q_{ij}$, we can express $H_{R}$ as:

[TABLE]

Because at zero temperature $\beta\rightarrow\infty$ the n-replica partition function is dominated by the configuration of $\{\mathbf{S}_{ik}\}$ minimizing $H_{R}$, we need to minimize

[TABLE]

and expand the solution of $\mathbf{S}_{ik}$ as powers of $\delta q_{ij}$. To achieve this, we write $\mathbf{S}_{ik}=\mathbf{S}_{k}+\mathbf{x}_{ik}$ in which $\mathbf{S}_{k}$ is given by the replica-symmetric mean-field solution and $\mathbf{x}_{ik}$ represents the contribution of the RSB term. Up to second order of $\mathbf{x}_{ik}$, we obtain:

[TABLE]

After minimizing $H_{R}$ w.r.t. $\mathbf{x}_{ik}$, we have the solution:

[TABLE]

Substitute the solution into $H_{\rm{eff}}$, we find that:

[TABLE]

We can conclude that the at $\beta\rightarrow\infty$ the quadratic term of $H_{R}$ is

[TABLE]

in which $z\geq 0$ obviously. For the quadratic term $H_{\rm{eff}}^{(2)}$ being at least positive semi-definite, the replica-symmetric solution can be stable. We have checked the eigenvalues of the Hessian matrix at at zero temperature and they are given by $-z(n-2)$ and $-z(2n-2)$. At the limit $n\rightarrow 0$ they are all non-negative. Our results are consistent with the well-known results of replica-symmetry breaking [31].