Quasiclassical Theory of C4-Symmetric Magnetic Order in Disordered Multiband Metals

Dzero, Maxim; Khodas, Maxim

doi:10.3389/fphy.2020.00356

ORIGINAL RESEARCH article

Front. Phys., 08 September 2020
Sec. Condensed Matter Physics
Volume 8 - 2020 | https://doi.org/10.3389/fphy.2020.00356

Quasiclassical Theory of C₄-Symmetric Magnetic Order in Disordered Multiband Metals

Maxim Dzero¹

Maxim Khodas²^*

¹Department of Physics, Kent State University, Kent, OH, United States
²The Racah Institute of Physics, The Hebrew University of Jerusalem, Jerusalem, Israel

Recent experimental studies performed in the normal state of iron-based superconductors have discovered the existence of the C₄-symmetric (tetragonal) itinerant magnetic state. This state can be described as a spin density wave with two distinct magnetic vectors ${\vec{Q}}_{1}$ and ${\vec{Q}}_{2}$ . Given an itinerant nature of magnetism in iron-pnictides, we develop a quasiclassical theory of tetragonal magnetic order in disordered three-band metal with anisotropic band structure. Within our model we find that the C₄-symmetric magnetism competes with the C₂-symmetric state with a single $\vec{Q}$ magnetic structure vector. Our main results is that disorder promotes tetragonal magnetic state which is in agreement with earlier theoretical studies.

1. Introduction

Quasiclassical approach to interacting many-body systems has proved to be a powerful tool in describing their transport and thermodynamic properties. Within this method, the quantum mechanical averages of an operator corresponding to a physical quantity are replaced with the averages of its classical counterpart over all classical trajectories. Alternatively, one can formulate the quasiclassical theory by using the quasiclassical functions which are obtained from the quantum mechanical single-particle propagators by integrating them over all single particle energies. Qualitatively, for a superconductor with pairing gap Δ and quasiparticles with Fermi momentum p_F and Fermi velocity v_F, this procedure corresponds to averaging over the short length scales of the problem $~ p_{F}^{- 1}$ and retaining the physics at long scales ~ v_F/Δ. Quasiclassical theory was particularly useful in the comparatively recent analysis of the problem of far-from-equilibrium order parameter dynamics in charge-neutral superfluids [1–5].

Most recently, several non-trivial phenomena have been observed in a family of iron-based superconductors and their alloys [6]. One example of such phenomena is an observation of the peak in the penetration depth in BaFe₂(As_1−xP_x)₂ as a function of phosphorus concentration [7–10], in Ba_1−xK_xFe₂As₂ as a function of potassium concentration [11] and, most recently in Ba(Fe_1−xCo_x)₂As₂ as a function of cobalt concentration [12]. Another example is the experimental observation of the spin-density-wave order which is characterized by two magnetic ordering vectors, ${\vec{Q}}_{1}$ and ${\vec{Q}}_{2}$ , in various iron-based superconducting alloys [13–20].

Due to the fact that in iron-based superconductors the superconductivity is often observed near magnetic instability, quasiclassical approaches initially developed for the purely superconducting states have been re-formulated to specifically include the effects of competition between superconducting and magnetic phases as well as the effects of disorder [21–24]. The experimental observations of the peak in the London penetration depth remains only partially understood [25, 26] which provides an additional motivation to look for possible explanations of this effect.

In turn, the experimental discovery of the double- $\vec{Q}$ magnetic state in iron-based superconductors has lead to an appearance of many theoretical works discussing the emergence of this state and its various properties as well as its relation with other magnetic states [27–35]. Most recently, the effects of disorder on the stability of the single- and double- $\vec{Q}$ states have been discussed [36]. In particular, it was found that disorder leads to suppression of the single- $\vec{Q}$ state in favor the double- $\vec{Q}$ one.

Inspired by the earlier work on this problem, in this paper we use a slightly simplified version of the model introduced in reference [36] to formulate a quasiclassical theory of the double- $\vec{Q}$ state in iron-based superconductors. Specifically, we consider the disordered model which incorporates both interband and intraband disorder. In agreement with the earlier results [36], we find that when the interband disorder can be ignored, the intraband disorder promotes the emergence of the double- $\vec{Q}$ state.

This paper is organized as follows. In the next section II introduce the model Hamiltonian. Section III is devoted to the formulation of the quasiclassical approach with the derivation of the quasiclassical equations. In section IV contains the results of the Landau expansion for the free energy using the quasiclassical equations. Section V contains the discussion of the results and comments related to the further development of the presented formalism in the context of the physics of iron-based superconductors. Sections with acknowledgments and Appendix with some technical details conclude the paper.

2. Model

In what follows we first introduce the model Hamiltonian, which consists of three terms:

\begin{array}{l} \hat{H} = {\hat{H}}_{0} + {\hat{H}}_{sdw} + {\hat{H}}_{dis} . & (1) \end{array}

The first term on the right hand side of this expression is a single-particle Hamiltonian which describes the band-structure consisting of three bands: hole-like band at the Γ point and two electron-like bands centered at ${\vec{Q}}_{X} = (π, 0)$ , ${\vec{Q}}_{Y} = (0, π)$ of the two-dimensional Brillouin zone. We use the compact notations to write down $\hat{H}$ ₀ using the six-component spinor ${\hat{Ψ}}_{k}^{†} = ({\hat{c}}_{k ↑}^{†}, {\hat{c}}_{k ↓}^{†}, {\hat{d}}_{k ↑}^{†}, {\hat{d}}_{k ↓}^{†}, {\hat{f}}_{k ↑}^{†}, {\hat{f}}_{k ↓}^{†})$ :

\begin{array}{l} {\hat{H}}_{0} = \sum_{k} {\hat{Ψ}}_{k}^{†} (\begin{matrix} ε_{Γ} (k) {\hat{σ}}_{0} & 0 & 0 \\ 0 & ε_{X} (k) {\hat{σ}}_{0} & 0 \\ 0 & 0 & ε_{Y} (k) {\hat{σ}}_{0} \end{matrix}) {\hat{Ψ}}_{k}, & (2) \end{array}

where ${\hat{σ}}_{0}$ is a unit 2 × 2 matrix and single particle energy spectra are given by ε_Γ(k) = ξ_k, $ξ_{k} = ϵ_{0} - k^{2} / 2$ , ε_X(k) = −ξ_k + δ₀ + δ₂ cos 2ϕ_k, ε_Y(k) = −ξ_k + δ₀ − δ₂ cos 2ϕ_k, ϵ₀ is the energy which amounts to the off-set between the bands and k = (k cos ϕ_k, k sin ϕ_k). Here δ₀ is a parameter which is defined relative to the chemical potential μ and describes the deviation from nesting: the bands are perfectly nested when δ₀ = 0. Lastly, δ₂ is an anisotropy parameter which accounts for the ellipticity of the corresponding Fermi pockets [36].

The second term, $\hat{H}$ _sdw, appearing in (1) accounts for the spin-density-wave order within the mean-field approximation:

\begin{array}{l} {\hat{H}}_{sdw} = - \sum_{k} {\hat{Ψ}}_{k}^{†} (\begin{matrix} 0 & {\vec{m}}_{X} \cdot \vec{σ} & {\vec{m}}_{Y} \cdot \vec{σ} \\ {\vec{m}}_{X} \cdot \vec{σ} & 0 & 0 \\ {\vec{m}}_{Y} \cdot \vec{σ} & 0 & 0 \end{matrix}) {\hat{Ψ}}_{k} . & (3) \end{array}

Here ${\vec{m}}_{X}$ , ${\vec{m}}_{Y}$ are the magnetizations corresponding to two structure vectors ${\vec{Q}}_{X}$ and ${\vec{Q}}_{Y}$ . In what follows, we will assume that magnetic state has Ising-like anisotropy, so we replace ${\vec{m}}_{X, Y} \cdot \vec{σ} \to m_{X, Y} {\hat{σ}}_{3}$ . Within the mean-field approach we have adopted here, the order parameters m_X,Y must be computed self-consistently.

Finally, the last term on the r.h.s. side of Equation (1) introduces the disorder potential in a system. In principle, the disorder should scatter quasiparticles within each band (intraband scattering) as well as between the bands (interband scattering). The disorder unavoidably leads to the suppression of itinerant magnetism. In this paper we will limit ourselves to the case of an intraband disorder only, for an interband disorder scattering only plays a crucial role in the problem of co-existence of magnetism and superconductivity [21–23], while for the problem at hand it will only lead the faster suppression of the magnetic order. Thus, we write for the last term in (1)

\begin{array}{l} {\hat{H}}_{dis} = u \int d^{2} r Ψ^{†} (r) Ψ (r) \sum_{i} δ (r - R_{i}) & (4) \end{array}

and the summation is performed over the impurity sites.

3. Quasiclassical Equations

In order to formulate the quasiclassical theory, we first introduce a single-particle correlation function

\begin{array}{l} G_{α β} (x, x^{'}) = - {〈 {\hat{T}}_{τ} ({\hat{Ψ}}_{α} (x) {\hat{Ψ}}_{β}^{†} (x^{'})) 〉}_{g.s.} & (5) \end{array}

in the Matsubara representation, ${\hat{Ψ}}_{α} (x) = {\hat{Ψ}}_{α} (r, τ)$ , and the averaging is performed over the ground state of the Hamiltonian (Equation 1). Next step consists in employing the equations of motion for the propagator (5):

\begin{array}{l} - \frac{\partial}{\partial τ} \hat{G} - {\hat{H}}_{r} \hat{G} - \hat{Σ} ○ \hat{G} = δ (x - x^{'}) 𝟙, \\ \frac{\partial}{\partial τ^{'}} \hat{G} - \hat{G} {\hat{H}}_{r^{'}} - \hat{G} ○ \hat{Σ} = δ (x - x^{'}) 𝟙 . & (6) \end{array}

Here $\hat{H}$ _r acts on r, the self-energy part $\hat{Σ}$ is generated by the disorder potential and its action on the propagator is

\begin{array}{l} \hat{Σ} ○ \hat{G} = \int_{0}^{1 / T} d τ^{″} \int d^{2} r^{″} Σ_{α γ} (x, x^{″}) G_{γ β} (x^{″}, x^{'}) . & (7) \end{array}

The summation over the repeated indices is assumed. Next, we perform the Wigner transformation

\begin{array}{l} \hat{G} (x, x^{'}) = \int \frac{d^{2} k}{{(2 π)}^{2}} e^{i k \cdot (r - r^{'})} \hat{G} (τ - τ^{'}; \frac{r + r^{'}}{2}, k) . & (8) \end{array}

In the presence of the quenched disorder, propagators will be dependent on R = (r + r′)/2. In what follows we assume that the disorder in uncorrelated and will average the propagator over the disorder distribution which corresponds to self-consistent Born approximation. Lastly, we introduce the following matrices:

\begin{array}{l} {\hat{M}}_{1} = (\begin{matrix} 0 & 0 & 0 \\ 0 & {\hat{σ}}_{0} & 0 \\ 0 & 0 & {\hat{σ}}_{0} \end{matrix}), {\hat{M}}_{2} = (\begin{matrix} 0 & 0 & 0 \\ 0 & {\hat{σ}}_{0} & 0 \\ 0 & 0 & - {\hat{σ}}_{0} \end{matrix}), \\ {\hat{M}}_{3} = (\begin{matrix} {\hat{σ}}_{0} & 0 & 0 \\ 0 & - {\hat{σ}}_{0} & 0 \\ 0 & 0 & - {\hat{σ}}_{0} \end{matrix}), {\hat{P}}_{X} = (\begin{matrix} 0 & {\hat{σ}}_{3} & 0 \\ - {\hat{σ}}_{3} & 0 & 0 \\ 0 & 0 & 0 \end{matrix}), \\ {\hat{P}}_{Y} = (\begin{matrix} 0 & 0 & {\hat{σ}}_{3} \\ 0 & 0 & 0 \\ - {\hat{σ}}_{3} & 0 & 0 \end{matrix}), {\hat{Q}}_{X} = (\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & {\hat{σ}}_{0} \\ 0 & {\hat{σ}}_{0} & 0 \end{matrix}) . & (9) \end{array}

Quasiclassical equations can now be derived after we multiply the first equation (6) from the left and the second equation from the right by ${\hat{M}}_{3}$ . Subtracting the resulting first equation from the second one we find

\begin{array}{l} [ω_{n} {\hat{M}}_{3}, \hat{G} (i ω_{n}, ϕ_{k})] + i δ_{0} [{\hat{M}}_{1}, \hat{G} (i ω_{n}, ϕ_{k})] \\ + i δ_{2} cos (2 ϕ_{k}) [{\hat{M}}_{2}, \hat{G} (i ω_{n}, ϕ_{k})] \\ + i [({\hat{H}}_{sdw} + {\hat{Σ}}_{dis} (i ω_{n})) {\hat{M}}_{3}, \hat{G} (i ω_{n}, ϕ_{k})] = 0, & (10) \end{array}

where we introduced the quasiclassical function, $[\hat{f}, \hat{g}]$ implies the usual commutation relation and

\begin{array}{l} \hat{G} (i ω_{n}, ϕ_{k}) = \frac{i}{π} \int {\hat{M}}_{3} \hat{G} (i ω_{n}, k) d ξ_{k}, & (11) \end{array}

where $\hat{G}$ (iω_n, k) is the Matsubara transform of the matrix function $\hat{G}$ (τ, k). The self-energy part is determined by the quasiclassical function and disorder scattering rate $Γ = π ν_{F} | u |^{2}$ (ν_F is the density of states at the Fermi level per valley per spin):

\begin{array}{l} {\hat{Σ}}_{dis} (i ω) = - i Γ \int_{0}^{2 π} \frac{d ϕ}{2 π} {\hat{M}}_{3} \hat{G} (i ω, i ϕ_{k}) . & (12) \end{array}

In the absence of disorder, quasiclassical equation (10) is linear in $\hat{G}$ and therefore is not sufficient to find $\hat{G}$ unambiguously. In order to define the problem completely, one has to complement (11) with a certain constraint. To derive this constraint, we introduce a new (matrix) function [37]

\begin{array}{l} \hat{B} (τ, τ^{'}; ϕ_{k}) = \int_{0}^{1 / T} \hat{G} (τ, τ^{″}; ϕ_{k}) \hat{G} (τ^{″}, τ^{'}; ϕ_{k}) d τ^{″} . \end{array}

Equation for this matrix function can be easily derived from (10). It then follows that quasiclassical functions must satisfy the following normalization condition:

\begin{array}{l} {\hat{G}}^{2} (i ω_{n}, ϕ_{k}) = 𝟙 . & (13) \end{array}

As we will demonstrate below, inclusion of disorder potential does not violate this condition. In order to solve the quasiclassical equations (10) self-consistently, we need to specify the matrix structure of the function $\hat{G}$ .

3.1. Clean System

We start by setting the disorder scattering rate to zero, Γ = 0, for it would allow us to keep the resulting expressions more compact. Most of the results derived in this section are easily generalized for the case when Γ ≠ 0 (see below).

In the absence of the magnetic order, the expression for the function $\hat{G}$ follows from (11) by comparing the solution of the quasiclassical equations with the expression found from the expression for the single-particle propagator, so that a term proportional to ${\hat{M}}_{3}$ must appear in the expression for $\hat{G}$ . This conjecture also implies that there should also appear two other terms proportional to ${\hat{M}}_{1}$ and ${\hat{M}}_{2}$ so we write the following ansatz

\begin{array}{l} {\hat{G}}_{0} = \frac{1}{2} (g_{1} + g_{2}) {\hat{M}}_{1} + \frac{1}{2} (g_{2} - g_{1}) {\hat{M}}_{2} + g_{3} {\hat{M}}_{3} . & (14) \end{array}

The commutators which include $\hat{H}$ _sdw must lead to the appearance of the three more terms in $\hat{G}$ : each one of the two of them being proportional to the corresponding magnetizations, while the third one being proportional to the product of m_X and m_Y. The calculation yields the following expression

\begin{array}{l} \hat{G} - {\hat{G}}_{0} = p_{x} {\hat{P}}_{X} + p_{y} {\hat{P}}_{Y} + q_{x} {\hat{Q}}_{X} . & (15) \end{array}

After plugging this ansatz into the quasiclassical equations and collecting the terms proportional to the same matrices (these matrices are different from those introduced above and will not be listed here), we derive the following set of quasiclassical equations:

\begin{array}{l} [2 i Ω_{n} + δ_{2} cos (2 ϕ)] p_{x} + m_{X} (g_{2} - 2 g_{3}) = - m_{Y} q_{x}, \\ [2 i Ω_{n} - δ_{2} cos (2 ϕ)] p_{y} + m_{Y} (g_{1} - 2 g_{3}) = - m_{X} q_{x}, \\ 2 δ_{2} cos (2 ϕ) q_{x} = m_{Y} p_{x} - m_{X} p_{y} & (16) \end{array}

and Ω_n = ω_n − iδ₀/2. Furthermore, given the expression (15) the constraint condition (13) reduces to the set of the following simple relations:

\begin{array}{l} q_{x} = - g_{1} \frac{p_{y}}{p_{x}} = - g_{2} \frac{p_{x}}{p_{y}}, {(g_{3} - g_{1} - g_{2})}^{2} = 1, \\ p_{x}^{2} = g_{1} (2 g_{3} - g_{1} - g_{2}), p_{y}^{2} = g_{2} (2 g_{3} - g_{1} - g_{2}) . & (17) \end{array}

Note, that by combining the first two relations with the last two ones one also finds $q_{x}^{2} = g_{1} g_{2}$ . With the help of relations (17) it is also straightforward to show that the third equation in (16) is redundant, so overall we have got the system of six non-linear equations with six unknowns. These equations must also be supplemented by the self-consistency conditions for the magnetizations, which in terms of the quasiclassical functions have the following form:

\begin{array}{l} m_{X, Y} = - 2 π ν_{F} g_{sdw} T Im \sum_{ω_{n} > 0} 〈 p_{x, y} (i ω_{n}, ϕ_{k}) 〉, & (18) \end{array}

where 〈f〉 denotes averaging over ϕ_k and g_sdw is the coupling constant.

The first two quasiclassical equations (16) can be re-written in a compact form using relations (17). Indeed, by introducing the auxiliary variables

\begin{array}{l} u_{1} = \frac{m_{X} \sqrt{g_{1}}}{\sqrt{2 + g_{1} + g_{2}}}, u_{2} = \frac{m_{Y} \sqrt{g_{2}}}{\sqrt{2 + g_{1} + g_{2}}}, & (19) \end{array}

the quasiclassical equations acquire the following form

\begin{array}{l} u_{1} [2 i Ω_{n} + δ_{2} cos (2 ϕ) - u_{1} - u_{2}] = m_{X}^{2}, \\ u_{2} [2 i Ω_{n} - δ_{2} cos (2 ϕ) - u_{1} - u_{2}] = m_{Y}^{2} . & (20) \end{array}

Perhaps, for the clarity of our subsequent discussion it would be useful to mention that in the case when magnetizations are vanishingly small, m_X,Y ≪ πT, functions $g_{1, 2} \propto m_{X, Y}^{2}$ , q_x ∝ m_Xm_Y, while p_x,y ∝ m_X,Y.

We have to analyze the solution of the equations (20) in two special cases only: (i) single- $\vec{Q}$ state for which we set m_Y = 0 and m_X = m₁ and (ii) double- $\vec{Q}$ state in which $m_{X} = m_{Y} = m_{2} / \sqrt{2}$ .

3.1.1. Single- $\vec{Q}$ State

Since in this case p_y = q_x = g₂ = 0, we have

\begin{array}{l} p_{x} (i ω_{n}, ϕ) = (\frac{m_{1}^{2} + u_{1}^{2}}{m_{1}^{2} - u_{1}^{2}}) \frac{2 m_{1}}{2 i Ω_{n} + δ_{2} cos (2 ϕ)} . & (21) \end{array}

In turn, function u₁(iω_n, ϕ) is determined by one of the two roots of the quadratic equation [first equation in (20) with u₂ = 0] which recovers the correct expression for the non-interacting propagator:

\begin{array}{l} u_{1} = Z_{n} (ϕ) - γ \sqrt{Z_{n}^{2} (ϕ) - m_{1}^{2}}, & (22) \end{array}

where Z_n(ϕ) = iΩ_n + (δ₂/2) cos(2ϕ) and γ is the pre-factor which guarantees that in the limit when m₁ → 0, u₁ also vanishes.

3.1.2. Double- $\vec{Q}$ State

The solution of the equations (20) in this cases reduces to the solution of a single cubic equation

\begin{array}{l} (x + 2 i Ω_{n}) [x^{2} - δ_{2}^{2} {cos}^{2} (2 ϕ) + m_{2}^{2}] = 2 i Ω_{n} m_{2}^{2} . & (23) \end{array}

Functions u₁ and u₂ can then be computed from

\begin{array}{l} u_{1, 2} = \frac{1}{2} (1 + \frac{2 i Ω_{n}}{x_{a}}) [x_{a} \pm δ_{2} cos (2 ϕ)], & (24) \end{array}

where x_a is one of the roots of Equation (23).

It is a priori not clear which one of the three roots must be chosen. An additional difficulty in choosing the correct root consists in the fact that after finding an analytic expressions for the roots (23) it turns out that depending on the limiting case (m_X,Y → 0 or δ₂ → 0, for example) different roots recover the correct expressions for the quasiclassical functions. The procedure we have adopted consisted in analyzing all three complex roots of (23) and picking up the one for which all the equations (16,17) are satisfied and in addition Im[p_x,y] < 0. The latter condition guarantees the positive contribution to magnetization (Equation 18), and minimum in free energy.

3.1.3. Results

We have used these expressions to evaluate the dependence of the order parameters m₁ and m₂ on the anisotropy parameter δ₂ for a fixed value of δ₀ and fixed temperature. Naturally, we find that both m₁ and m₂ are the same for the same values of the model parameters. The results of the calculations for the temperature dependence of the magnetizations m₁ and m₂ are presented on Figure 1A. Perhaps it is not too surprising that we found the values of m₁ and m₂ equal to each other within the error bars of the numerical calculations. Therefore, self-consistency equations cannot be used to determine which of the two states would be more favorable and we will have to compute the free energy for each state.

FIGURE 1

Figure 1. Dependence of magnetizations m₁ and m₂ on temperature and anisotropy parameter δ₂ obtained by the numerical analysis of the self-consistency equations (18). Within the numerical accuracy, we found that m₁ = m₂. (A) Magnetization as a function of temperature are plotted for various values of scattering rate Γ and δ₂ = 0.5δ₀ with δ₀ = 2πT_s where T_s is a Curie temperature in isotropic system without impurities. (B) Magnetization as a function of the Fermi surface anisotropy parameter δ₂ and T = 0.1T_s.

3.2. Disordered System

Quasiclassical equations for the disordered system naturally have similar form as equations (16) for the fact that the matrix structure of the quasiclassical function does not change as soon as Γ becomes non-zero. The calculation of the commutation relations (10) yields

\begin{array}{l} [2 i {\tilde{Ω}}_{n} + Δ_{n} (ϕ)] p_{x} + (m_{X} - i Γ 〈 p_{x} 〉) (g_{2} - 2 g_{3}) \\ = i Γ 〈 q_{x} 〉 p_{y} - (m_{Y} - i Γ 〈 p_{y} 〉) q_{x}, \\ [2 i {\tilde{Ω}}_{n} - Δ_{n} (ϕ)] p_{y} + (m_{Y} - i Γ 〈 p_{y} 〉) (g_{1} - 2 g_{3}) \\ = i Γ 〈 q_{x} 〉 p_{x} - (m_{X} - i Γ 〈 p_{x} 〉) q_{x} . & (25) \end{array}

In these equations ${\tilde{Ω}}_{n} = Ω_{n} + (Γ / 4) (4 〈 g_{3} 〉 - 〈 g_{1} 〉 - 〈 g_{2} 〉)$ and Δ_n(ϕ) = δ₂ cos(2ϕ) + i(Γ/2)(〈g₁〉 − 〈g₂〉). Just like in the case Γ = 0 the third equation is redundant and therefore is not listed here.

Equations (25) show that disorder renormalization plays out differently for single- $\vec{Q}$ and double- $\vec{Q}$ states. Given these disorder renormalizations, in order to solve the self-consistency equation (18), the angular averages (〈g₃〉 and 〈p_x〉 in a single- $\vec{Q}$ state, for example) had to be computed by iterations. We found that the values of the corresponding magnetizations still remain essentially identical for non-zero Γ (Figure 1B). We also found, that qualitative behavior of both m₁(δ₂) and m₂(δ₂) does not change with an inclusion of disorder.

Lastly, we would like to mention that the inclusion of the interband disorder with scattering rate Γ_π would not change the dependence of the magnetization on the anisotropy parameters, but only leads to a faster suppression of the magnetization with an increase in Γ_π.

4. Free Energy

To derive an expression for the free energy in terms of the quasiclassical functions, we can employ an expression for the effective action corresponding to the model Hamiltonian (1). Omitting the disorder potential for now, we have [36] $F (m_{X}, m_{Y}) = (m_{X}^{2} + m_{Y}^{2}) / g_{sdw} - S (λ = 1)$ with

\begin{array}{l} S (λ) = T \sum_{i ω_{n}} \int_{k} Tr log (\hat{𝟙} + λ {\hat{G}}_{0} (i ω_{n}, k) \hat{W}) . & (26) \end{array}

Here $\hat{G}$ ₀(iω_n, k) is the single-particle propagator for the non-interacting system, $\hat{W} = - m_{X} {\hat{S}}_{X} - m_{Y} {\hat{S}}_{Y}$ and

\begin{array}{l} {\hat{S}}_{X} = (\begin{matrix} 0 & {\hat{σ}}_{3} & 0 \\ {\hat{σ}}_{3} & 0 & 0 \\ 0 & 0 & 0 \end{matrix}), {\hat{S}}_{Y} = (\begin{matrix} 0 & 0 & {\hat{σ}}_{3} \\ 0 & 0 & 0 \\ {\hat{σ}}_{3} & 0 & 0 \end{matrix}) . & (27) \end{array}

The expression for the free energy in terms of the quasiclassical functions can be derived by following the steps in the calculation of reference [38]. First, we note

\begin{array}{l} \frac{\partial S}{\partial λ} = - i π ν_{F} \int_{0}^{π} \frac{d ϕ_{k}}{π} Tr [{\hat{M}}_{3} {\hat{G}}_{λ} (i ω_{n}, ϕ_{k}) \hat{W}], & (28) \end{array}

where ${\hat{G}}_{λ} (i ω_{n}, ϕ_{k})$ is found from solving the quasiclassical equations (10) in which order parameters have been rescaled by parameter λ, m_X,Y → λm_X,Y. The resulting expression for the free energy reads

\begin{array}{l} F (m_{X}, m_{Y}) = \frac{m_{X}^{2} + m_{Y}^{2}}{g_{sdw}} \\ - 2 i π ν_{F} T \sum_{i ω_{n}} \int_{0}^{1} d λ (m_{X} 〈 p_{λ x} 〉 + m_{Y} 〈 p_{λ y} 〉) . & (29) \end{array}

This expression can also be employed for the case of non-zero disorder by using the solution of equations (25) with the rescaled magnetizations.

It is a hopeless task to evaluate the free energy (29) analytically, but it is amenable to the numerical analysis. However, our numerical computation of the free energy for the single-Q and double-Q states ran into an unexpected problem: the difference between the free energies of the corresponding states fall within the numerical error of the calculation. Thus, in order to determine which one of the two magnetic states will be energetically favorable, below we derive the Landau expansion.

4.1. Free Energy Expansion in Powers of the Magnetization

Having found an expression for the free energy, we consider the temperatures slightly below the critical temperature, so that both magnetizations are sufficiently small compared to πT. Then, we can formally obtain the solution of the quasiclassical equations (25) by expanding functions p_x and p_y in powers of m_X and m_Y.

4.1.1. Clean Case

In the case of the clean system the expression up to the fourth order in powers of magnetization reads

\begin{array}{l} F (m_{X}, m_{Y}) = a_{2} (m_{X}^{2} + m_{Y}^{2}) + b_{4} {(m_{X}^{2} + m_{Y}^{2})}^{2} \\ - g_{4} {(m_{X}^{2} - m_{Y}^{2})}^{2} + O (m^{6}), & (30) \end{array}

where the corresponding coefficients are given by b₄ = (a₄ + a_XY)/2, g₄ = (a_XY − a₄)/2 with

\begin{array}{l} a_{2} = \frac{1}{g_{sdw}} - 8 ν_{F} T \sum_{ω_{n} > 0} \int_{0}^{π} \frac{ω_{n} d ϕ_{k}}{4 ω_{n}^{2} + {(δ_{0} \pm δ_{2} cos (2 ϕ_{k}))}^{2}}, \\ a_{4} = 4 ν_{F} Im {T \sum_{ω_{n} > 0} \int_{0}^{π} \frac{d ϕ_{k}}{{[2 i Ω_{n} \pm δ_{2} cos (2 ϕ_{k})]}^{3}}}, \\ a_{X Y} = - 8 ν_{F} Re {T \sum_{ω_{n} > 0} \int_{0}^{π} \frac{Ω_{n} d ϕ_{k}}{{[4 Ω_{n}^{2} + δ_{2}^{2} {cos}^{2} (2 ϕ_{k})]}^{2}}} . \end{array}

The sign of the coefficient g₄ is crucial for it determines which one of the two states becomes energetically more favorable. Indeed, let us assume that we choose the model parameters such that both m₁ and m₂ are much smaller than πT. For a fixed value of m₁ = m₂ it follows that when g₄ > 0 the single- $\vec{Q}$ will have the lower energy compared to the double- $\vec{Q}$ one. However, one needs to keep in mind that this line of arguments holds only when the coefficients in the free energy expansion are all of the order O(1) and coefficient b₄ remains positive for a given set of values of parameters δ₀/2πT and δ₂/2πT.

4.1.2. Disordered Case

The question arises as to how non-zero disorder will affect the stability of the single- $\vec{Q}$ state [36]. The calculation of the quasiclassical functions is similar to the one in the clean case, with the only exception that the averages over the angle ϕ need to be computed self-consistently. For example, the first order corrections to functions p_x and p_y are

\begin{array}{l} p_{x}^{(1)} = \frac{2 (m_{X} - i Γ 〈 p_{x}^{(1)} 〉)}{2 i (| ω_{n} | + Γ) sign (ω_{n}) + δ_{0} + δ_{2} cos (2 ϕ)}, \\ p_{y}^{(1)} = \frac{2 (m_{Y} - i Γ 〈 p_{y}^{(1)} 〉)}{2 i (| ω_{n} | + Γ) sign (ω_{n}) + δ_{0} - δ_{2} cos (2 ϕ)} . & (31) \end{array}

After integrating both parts of these expressions over ϕ, we can easily solve for $〈 p_{x}^{(1)} 〉$ and $〈 p_{y}^{(1)} 〉$ .

The calculation of the expressions for the coefficients of the Landau expansion in this case gives

\begin{array}{l} A_{4} = 4 ν_{F} Im {T \sum_{ω_{n} > 0} {(\frac{η (i ω_{n}) - i Γ}{η (i ω_{n}) + i Γ})}^{3} \\ \times \int_{0}^{π} \frac{d ϕ_{k}}{{[2 i (Ω_{n} + Γ) \pm δ_{2} cos (2 ϕ_{k})]}^{3}}}, \\ A_{XY} = - 8 ν_{F} Re {T \sum_{ω_{n} > 0} {(\frac{η (i ω_{n}) - i Γ}{η (i ω_{n}) + i Γ})}^{3} \\ \times \int_{0}^{π} \frac{(Ω_{n} + Γ) d ϕ_{k}}{{[4 {(Ω_{n} + Γ)}^{2} + δ_{2}^{2} {cos}^{2} (2 ϕ_{k})]}^{2}}} \\ - 2 π ν_{F} Γ Re {T \sum_{ω_{n} > 0} {(\frac{η (i ω_{n}) - i Γ}{η (i ω_{n}) + i Γ})}^{3} z^{2} (i ω_{n})} . & (32) \end{array}

Functions η(iω_n) and z(iω_n) appear as a result of disorder renormalization and are listed in Appendix. The coefficient g₄ in free energy is now given by g₄ = (A_XY − A₄)/2. Compared with the clean case, we see that expression for the coefficient A_XY contains an extra term proportional to Γ. The dependence of g₄ on disorder can be easily analyzed numerically. The results of the numerical computations are shown in Figure 4.

4.2. Phase Diagram

To determine the phase diagram in the space of anisotropy parameters δ₀ and δ₂, we need to find a point where the free energies of both states become degenerate, g₄(δ_0c, δ_2c) = 0. In Figure 2 we show the phase diagram for the clean system. It agrees qualitatively with the one obtained previously [36] for small values of δ₂/2πT ≪ 1, single- $\vec{Q}$ state becomes energetically favorable when the value of electron-hole asymmetry δ₀ is above a critical value δ_0c/2πT ~ 0.3.

FIGURE 2

Figure 2. Results of the numerical analysis of the coefficient g₄ in the free energy expansion for the clean system, Γ = 0. The solid line marks the first order transition line along which the coefficient g₄ is zero and the energies of the single- $\vec{Q}$ and double- $\vec{Q}$ states are degenerate.

FIGURE 3

Figure 3. Plot of the nematic coupling constant g₄ which appears in free energy as a function of the disorder scattering rate for various values of the parameters δ₀ and δ₂.

With an addition of disorder, phase diagram is modified and the results are presented on Figure 3 For small disorder the critical line separating two phases slightly moves to higher values of δ₂. Perhaps unexpectedly, a small region of single- $\vec{Q}$ state appears at large (compared to δ₂) values of δ₀. Upon further increase in the values of the disorder scattering rate, the phase boundary separating two states moves to higher values of δ₂ and also extends to higher values of δ₀. Overall, we may conclude that disorder promotes double- $\vec{Q}$ state over the single- $\vec{Q}$ state.

FIGURE 4

Figure 4. Results of the numerical analysis of the coefficient g₄ in the free energy expansion for the disordered system. With the increase in the value of the disorder scattering rate, the single- $\vec{Q}$ state is a ground state for higher and higher values of the anisotropy parameter δ₂ which accounts for the ellipticity of the electron-like pockets.

5. Discussion

As we have already pointed out in the Introduction, our main goal was to demonstrate how the quasi-classical method can be applied to analyse the competition between magnetic states in multiband metals in the presence of disorder. Having accomplished that goal, we can now generalize it to investigate the problem of an interplay between superconductivity and magnetism. It is already well-established that by including the interband disorder scattering Anderson-Abrikosov-Gor'kov theorem makes it possible for superconductivity and magnetism to co-exist in a certain region of the phase diagram, which size is determined by the ratio of the intra- and inter-band scattering rates [22, 23]. The question is then would be to check if superconducting order may provide an additional contribution in determining which of the two competing magnetic states would be energetically favorable. These results may be employed to provide a qualitative understanding as to why nematicity has been observed in stoichiometric iron selenide in contrast to electron-doped iron selenide.

Lastly, we would like to mention that the inclusion of the interband disorder scattering would not affect our results in any substantial way. Indeed, compared to the case of intraband disorder, the inclusion of the interband scattering leads primarily to the faster suppression of the critical temperature, without affecting the ground state energies of the single- and double- $\vec{Q}$ states significantly.

To summarize, in this paper we have formulated the quasi-classical approach to analyze the relative stability of the single- and double- $\vec{Q}$ spin-density-wave states with respect to band and effective mass anisotropy as well as disorder scattering. Generally, we find that with an increase in intraband disorder scattering rate, the system favors the single- $\vec{Q}$ for moderately high values of the Fermi surface anisotropy parameter, δ₂.

Data Availability Statement

All datasets generated for this study are included in the article/supplementary material.

Author Contributions

MD has performed numerical and analytical calculations. MK has performed the analytical calculations. All authors participated in the discussions of the ideas and the ways to solve the problem, and participated in writing the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was financially supported by the National Science Foundation grant NSF-DMR-BSF-2002795, and in part, by the U.S. Department of Energy, Basic Energy Sciences, grant DE-SC0016481 (MD). MK was supported by the Israel Science Foundation 1287/15.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Alex Levchenko for bringing this problem to our attention and many fruitful conversations. Useful discussions with R. M. Fernandes and E. König were gratefully acknowledged.

References

1. Barankov RA, Levitov LS, Spivak BZ. Collective Rabi oscillations and solitons in a time-dependent BCS pairing problem. Phys Rev Lett. (2004) 93:160401. doi: 10.1103/PhysRevLett.93.160401

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Yuzbashyan EA, Altshuler BL, Kuznetsov VB, Enolskii VZ. Solution for the dynamics of the BCS and central spin problems. J Phys A Math Gen. (2005) 38:7831. doi: 10.1088/0305-4470/38/36/003