Quasiclassical theory of $C_4$-symmetric magnetic order in disordered multiband metals

Recent experimental studies performed in the normal state of iron-based superconductors have discovered the existence of the $C_4$-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_4$-symmetric magnetism competes with the $C_2$-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.


I. 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 −1 F and retaining the physics at long scales ∼ v F /∆. Quasiclassical theory was particularly useful in the comparatively recent analysis of the problem of farfrom-equilibrium order parameter dynamics in chargeneutral superfluids. [1][2][3][4][5] Most recently, several nontrivial 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 2 (As 1−x P x ) 2 as a function of phosphorus concentration 7-10 , in Ba 1−x K x Fe 2 As 2 as a function of potassium concentration 11 and, most recently in Ba(Fe 1−x Co x ) 2 As 2 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, Q 1 and Q 2 , in various alloys iron-based superconducting alloys. [13][14][15][16][17][18][19][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 reformulated to specifically include the effects of competition between superconducting and magnetic phases as well as the effects of disorder. [21][22][23][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-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][28][29][30][31][32][33][34][35] Most recently, the effects of disorder on the stability of the singleand double-Q states have been discussed. 36 In particular, it was found that disorder leads to suppression of the single-Q state in favor the the double-Q one.
Inspired by the earlier work on this problem, in this paper we use a slightly simplified version of the model introduced in Ref. [36] to formulate a quasiclassical theory of the double-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-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 ironbased superconductors. Sections with acknowledgements and Appendix with some technical details conclude the paper.

II. MODEL
In what follows we first introduce the model Hamiltonian, which consists of three terms: 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: holelike band at the Γ point and two electron-like bands centered at Q X = (π, 0), Q Y = (0, π) of the twodimensional Brillouin zone. We use the compact notations to write downĤ 0 using the six-component spinor whereσ 0 is a unit 2 × 2 matrix and single particle energy spectra are given by 0 is the energy which amounts to the off-set between the bands and k = (k cos φ, k sin φ). Here δ 0 is an anisotropy parameter which is defined relative to the chemical potential µ, so that the bands are perfectly nested when δ 0 = 0. Lastly, δ 2 is an anisotropy parameter which accounts for the ellipticity of the corresponding Fermi pockets. 36 The second term,Ĥ sdw , appearing in (1) accounts for the spin-density-wave order within the mean-field approximation: Here m X , m Y are the magnetizations corresponding to two structure vectors Q X and Q Y . In what follows, we will assume that magnetic state has Ising-like anisotropy, so we replace m X,Y · σ → m X,Yσ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 Eq. (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) and the summation is performed over the impurity sites.

III. QUASICLASSICAL EQUATIONS
In order to formulate the quasiclassical theory, we first introduce a single-particle correlation function in the Matsubara representation,Ψ α (x) =Ψ α (r, τ ), and the averaging is performed over the ground state of the Hamiltonian, Eq. (1). Next step consists in employing the equations of motion for the propagator (5): HereĤ r acts on r, the self-energy partΣ is generated by the disorder potential and its action on the propagator iŝ The summation over the repeated indices is assumed. Next, we perform the Wigner transformation 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: Quasiclassical equations can now be derived after we multiply the first equation (6) from the left and the second equation from the right byM 3 . Subtracting the resulting first equation from the second one we find where we introduced the quasiclassical function, [f,ĝ] implies the usual commutation relation and 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):Σ Quasiclassical equation (10) is linear inĜ and therefore is not sufficient to findĜ 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 Equation for this matrix function can be easily derived from (10). It then follows that quasiclassical functions must satisfy the following normalization condition: In order to solve the quasiclassical equations (10) selfconsistently, we need to specify the matrix structure of the functionĜ.

A. 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Ĝ 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 toM 3 must appear in the expression forĜ. This conjecture also implies that there should also appear two other terms proportional tô M 1 andM 2 so we write the following ansatẑ The commutators which includeĤ sdw must lead to the appearance of the three more terms inĜ: 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 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: and Ω n = ω n − iδ 0 /2. Furthermore, given the expression (15) the constraint condition (13) reduces to the set of the following simple relations: Note, that by combining the first two relations with the last two ones one also finds q 2 x = 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 selfconsistency conditions for the magnetizations, which in terms of the quasiclassical functions have the following form: where f denotes averaging over φ k and g sdw is the coupling constant. The first two quasiclassical equations (16) can be rewritten in a compact form using relations (17). Indeed, by introducing the auxiliary variables the quasiclassical equations acquire the following form 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 a. Single-Q state. Since in this case p y = q x = g 2 = 0, we have In turn, function u 1 (iω n , φ) is determined by one of the two roots of the quadratic equation (first equation in (20) with u 2 = 0) which recovers the correct expression for the non-interacting propagator: where Z n (φ) = iΩ n +(δ 2 /2) cos(2φ) and γ is the prefactor which guarantees that in the limit when m 1 → 0, u 1 also vanishes.
b. Double-Q state. The solution of the equations (20) in this cases reduces to the solution of a single cubic equation Functions u 1 and u 2 can then be computed from 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 δ 2 → 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, Eq. (18), and minimum in free energy.
c. Results. We have used thes expressions to evaluate the dependence of the order parameters m 1 and m 2 on the anisotropy parameter δ 2 for a fixed value of δ 0 and fixed temperature. Naturally, we find that both m 1 and m 2 are the same for the same values of the model parameters. The results of the calculations for the temperature dependence of the magnetizations m 1 and m 2 are presented on Fig. 1(a). Perhaps it is not too surprising that we found the values of m 1 and m 2 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.

B. 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 nonzero. The calculation of the commutation relations (10) yields In these equations Ω n = Ω n + (Γ/4)(4 g 3 − g 1 − g 2 ) and ∆ n (φ) = δ 2 cos(2φ) + i(Γ/2)( g 1 − g 2 ). 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-Q and double-Q states. Given these disorder renormalizations, in order to solve the self-consistency equation (18), the angular averages ( g 3 and p x in a single-Q state, for example) had to be computed by iterations. We found that the values of the corresponding magnetizations still remain essentially identical for nonzero Γ, Fig. 1(b). We also found, that qualitative behavior of both m 1 (δ 2 ) and m 2 (δ 2 ) 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 Γ π .

IV. 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 HereĜ 0 (iω n , k) is the single-particle propagator for the non-interacting system,Ŵ = −m XŜX − m YŜY and The expression for the free energy in terms of the quasiclassical functions can be derived by following the steps in the calculation of Ref. 38. First, we note whereĜ λ (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 This expression can also be employed for the case of nonzero 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.

A. 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 .
a. Clean case. In the case of the clean system the expression up to the fourth order in powers of magnetization reads where the corresponding coefficients are given by b 4 = (a 4 + a XY )/2, g 4 = (a XY − a 4 )/2 with The sign of the coefficient c 4 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 1 and m 2 are much smaller than πT . For a fixed value of m 1 = m 2 it follows that when g 4 > 0 the single-Q will have the lower energy compared to the double-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 4 remains positive for a given set of values of parameters δ 0 /2πT and δ 2 /2πT . b. Disordered case. The question arises as to how nonzero disorder will affect the stability of the single-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 2i (|ω n | + Γ) sign(ω n ) + δ 0 + δ 2 cos(2φ) , After integrating both parts of these expressions over φ, we can easily solve for p of the Landau expansion in this case gives Functions η(iω n ) and z(iω n ) appear as a result of disorder renormalization and are listed in Appendix. The coefficient g 4 in free energy is now given by g 4 = (A XY − A 4 )/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 4 on disorder can be easily analyzed numerically. The results of the numerical computations are shown in Fig. 3.

B. Phase diagram
To determine the phase diagram in the space of anisotropy parameters δ 0 and δ 2 , we need to find a point where the free energies of both states become degenerate, g 4 (δ 0c , δ 2c ) = 0. In Fig. 2 we show the phase diagram for the clean system. It agrees qualitatively with the one obtained previously: 36 for small values of δ 2 /2πT 1, single-Q state becomes energetically favorable when the value of electron-hole asymmetry δ 0 is above a critical value δ 0c /2πT ∼ 0.3.
With an addition of disorder, phase diagram is modified and the results are presented on Fig. 4 For small disorder the critical line separating two phases slightly moves to higher values of δ 2 . Perhaps unexpectedly, a small region of single-Q state appears at large (compared to δ 2 ) values of δ 0 . Upon further increase in the values of the disorder scattering rate, the phase boundary separating two states moves to higher values of δ 2 and also extends to higher values of δ 0 . Overall, we may conclude that disorder promotes double-Q state over the single-Q state.

V. 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-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-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-Q for moderately high values of the Fermi surface anisotropy parameter, δ 2 .

VI. ACKNOWLEDGMENTS
We would like to thank Alex Levchenko for bringing this problem to ou rattention and many fruitful conversations. Useful discussions with R. M. Fernandes and E. König are gratefully acknowledged. This work was financially supported by the U.S. Department of Energy, Basic Energy Sciences, grant de-sc0016481 (MD) and by the Israel Science Foundation, Grant No. 1287/15 (MK).
Appendix A: Coefficients in the free energy expansion I this Section we provide the details of the calculation for the Landau free energy expansion. Both p λx and p λy can be determined approximately for small values of m X and m Y from the quasiclassical equations. We start with the derivation for the clean case, Γ = 0.

Third order corrections
The second order correction to p λj is zero. To determine the third order correction, we first need to compute the second order corrections to g j 's. To do that, we first use equations (20) (and presume for simplicity that ω n > 0): so that λ3 = g (2) λ1 + g (2) λ2 . (A3) In addition, for the function q x we find q (2) λx = − 2λ 2 m X m Y (2iω n + δ 0 ) 2 − δ 2 2 cos 2 (2φ) .

(A4)
The choice of sign follows from considering the trivial case of δ 2 = 0.
Given all these expressions, we go back to equations (16) to obtain the following expression: (A6) After plugging these expressions into Eq. (29) and grouping the similar terms, we arrive to Eq. (30).