Magnetic orders in pnictide superconductors: the effect of biquadratic interaction

The undoped three-orbital spin-fermion model with the additional biquadratic interaction which arises from various factors, such as quantum or thermal fluctuation, is studied using the Monte Carlo method in order to investigate the spin orders in pnictide superconductors. The simulations show that the experimentally observed nematic state can be stabilized by the positive biquadratic interaction, suggesting that such interaction may be another origin for nematicity in pnictides in addition to the couplings to the lattice degrees of freedom. Furthermore, the so-called flux state at low temperatures is identified when a rather weak negative biquadratic interaction is introduced, which is the same as the phase predicted in earlier theoretical work.


Introduction
In the past few years, the magnetic orders of normal states of iron-based superconductors have drawn extensive attention since the electron pairing mechanism in these materials may be related to spin fluctuations [1][2][3][4]. Experimentally, it was reported that the long-range antiferromagnetic (AFM) order with a wave vector (π, 0) (inset of figure 1(a)) can be stabilized in the parent compounds of most pnictides below the Néel temperature T N . Furthermore, a structural transition from the tetragonal to the orthorhombic phase occurs at T S ⩾ T N [5][6][7].
In order to understand the close relationship and possible bifurcation between T N and T S , two different mechanisms involving the crucial component of spin [8,9] and orbital fluctuations [10,11] have been proposed, respectively. Specifically, it is suggested in the first scenario that the short range magnetic order (the nematic order) can be developed before the stabilization of the long range AFM order, and drives the structural transition [12]. This point of view is strongly supported by a recent work in which the evident scaling relation between the magnetic and lattice fluctuations has been revealed [13]. In addition, several theoretical attempts to understand the stabilization of the nematic order were reported [14]. For example, the classical spin model with a rather large biquadratic interaction has been explicitly studied using Monte Carlo (MC) simulations, and a small separation of the nematic and AFM transition temperatures was predicted [15,16]. This phenomenon has also been qualitatively reproduced in another report in which a quantum spin Heisenberg model for pnictides was investigated [12]. Furthermore, it has also been suggested in a recent theoretical work studying an itinerant electron model, that the orbital order in iron pnictides are induced by the nematic order [9].
In addition to the Heisenberg model, the spin-fermion (SF) model, which is somewhat similar to the double-exchange model for manganites, has also been proposed for studying the pnictide superconductors [17,18]. Both itinerant electrons and localized spins are included in the SF model, and the calculations have demonstrated the great success of this model for explaining complicated physics in these typical correlated electron systems. For example, not only the magnetic and spectral properties of the ground state [19], but also the in-plane resistivity anisotropy observed in BaFe 2 As 2 [20,21] can be qualitatively reproduced in MC simulations using the three-orbital SF model [22]. In addition, the lattice degrees of freedom have been introduced into this model to study the nematic state of pnictides [23]. It is suggested that both the spin-lattice and orbital-lattice couplings are indispensable in stabilizing the nematic order and lattice distortion observed in experiments.
However, the phase diagram of the isovalent-doping system BaFe 2 (As 1 − x P x ) 2 recently obtained by the magnetic torque measurement reveals that the nematicity (with transition point T*) can be developed well above both T N and T S [7]. To some extent, this experimental phenomenon suggests that these couplings to the lattice degrees of freedom may not be the only mechanism to stabilize the nematicity, and other possible factors should also be responsible for the nematicity. This viewpoint is partially supported by one of the earlier theories in which the structural and magnetic transitions were assumed to occur independently [24]. On the other hand, it is expected that the biquadratic interaction plays an essential role and can not be neglected in any model calculation describing these systems [25,26]. For pnictides, several possible origins for this term such as quantum fluctuation and thermal fluctuation have been proposed [9]. In fact, according to the electronic structure calculations, in iron pnictides, the biquadratic coupling of intrinsic electronic origin is much larger than that generated by thermal fluctuations or by magneto-structural coupling. As a matter of fact, this interaction has been included in several phenomenological Heisenberg models in earlier reports [12,16,27]. In these studies, only the localized spins are kept, while the itinerant character of electrons, which is essential in understanding the experimental phenomena such as the anisotropy of electrical transport, is completely ignored. In fact, the limitations of such a pure spin model was revealed in earlier reports, which indicated the importance of the itinerant character of electrons because most of these materials are metallic [16]. Overall, one may question whether the biquadratic interaction is really significantly important to the magnetic properties in real pnictides instead of in a toy model of pure spins. Therefore, here we carried out the MC study on the more sophisticated three-orbital SF model [28] with the additional biquadratic interaction. Our simulation finds that the experimentally reported T* > T N can be qualitatively reproduced when the positive biquadratic interaction is considered. Thus, our MC simulations suggest that such interaction may be another origin for the stabilization of the nematic order, in addition to those couplings to the lattice degrees of freedom which has been verified earlier [23]. Furthermore, it is demonstrated that instead of the well-known (π, 0) stripe AFM order, the exotic flux state [29,30] can be stabilized at low temperatures (T) when a rather weak negative biquadratic interaction is introduced to the SF model, in agreement with previous studies based on the classical spin model [16].

Model and method
In this work, we studied the three-orbital SF model with the biquadratic interaction. The Hamiltonian can be written as: In the above model Hamiltonian, the first term is the Fe-Fe hopping of itinerant electrons with fermion spin σ in the xz, yz and xy d orbitals denoted by α or β. The operator creates (annihilates) an electron at the orbital α of site i. Both nearest-neighbor (NN) and diagonal next-nearest-neighbor (NNN) hoppings are considered, as explicitly given in an earlier report [28]. The hopping amplitudes between different sites and orbitals are shown in table 1, and the detailed mathematical forms were described in the equations (1)-(3) of [28]. The second term is the AFM superexchange interaction between NN localized spins S with coupling J NN , while the third term is that between NNN localized spins with coupling J NNN . The fourth term is the Hund interaction between the localized spin and the itinerant-fermion one at site i. The last term is the biquadratic interaction J Bi between NN localized spins. Note that most of these parameters have been exactly calculated or reasonably estimated by the experimental data and first-principles calculations in earlier reports [31], and were consistent with the magnetic and spectral properties of the undoped compounds. In this work, these parameters are fixed to the same values as previously reported [22]. For example, J Hund is set to be 0.1 electron volts (eV), J NNN to be 0.01 eV, and the itinerant electron density to be 4/3. It has been demonstrated by linear response calculations and density-functional calculations that the NN exchange constant in iron pnictides is strongly anisotropic in the collinear state [32,33]. To mimic such anisotropy, the anisotropic NN couplings J NN x = 0.016 eV and J NN y = 0.014 eV along the x and y axis are adopted in the following simulations if not noted explicitly. The physical conclusion regarding the nematic state will not be affected by such a choice of anisotropic J NN , [34], which will be clarified in detail later. In contrast to previous studies [23], the lattice degrees of freedom are not included in this work, conveniently allowing us to focus on the effects of the biquadratic interaction.
The above modeled Hamiltonian was studied via a combination of exact diagonalization and MC simulation on a two-dimensional 8 × 8 lattice [35]. Typically, the first 10 4 MC steps are used for thermal equilibrium and another 2 × 10 4 MC steps are used for measurements at each T and for each set of parameters. To characterize the AFM spin orders, the spin structure factors at wave vector ⃗ q are calculated by performing Fourier transforms of the real-space correlation functions [36]: Table 1. Values of the Fe-Fe hopping amplitudes (t 1 -t 8 ) and the energy difference between the xy and the degenerate xz/yz orbitals Δ xy with the overall energy unit in electron volts.
where N is the number of sites. On the other hand, the spin nematic order can be characterized by short-range spin correlations, and the order parameter is defined as [37]: The magnetic and nematic transition temperatures (T N and T*) can be estimated from the positions of the peaks in the spin magnetic susceptibility χ Φ(π, 0) and the spin nematic susceptibility χ Ψ , respectively.

Results and discussion
3.1. Nematic phase above T N First, figure 1(a) shows the calculated Φ(π, 0) and Ψ for various T at J Bi = 0 as the original results without the biquadratic interaction. The two curves coincide with each other very well, indicating that the nematic and (π, 0) AFM orders are developed at identical T. This coincidence can be verified by the results of susceptibilities χ Φ(π, 0) and χ Ψ which are given in figure 1(b). It is clearly demonstrated from the positions of the peaks in χ Φ(π, 0) and χ Ψ that the nematic transition occurs simultaneously with the Néel transition at T N = T* = 112.6 K.
Then, the simulated results at J Bi = 0.003 eV are shown in figures 1(c) and (d). On the one hand, it is observed that both T N and T* are higher than the original ones at J Bi = 0, which can be straightforwardly understood from the energy landscape. Specifically, it is noted that the additional energy contribution due to the biquadratic interaction between NN spins can be written as H Bi = −J Bi cos 2 θ ij with θ ij as the angle between the NN spins. As a consequence, the NN spins tend to be parallel or antiparallel with each other to satisfy the positive biquadratic interaction, resulting in the stabilizations of the long range (π, 0) AFM order and the nematic order as J Bi increases. In other words, the effective magnetic coupling between NN spins are enhanced by this J Bi and stronger thermal fluctuations are needed to break these two orders, leading to the fact that both these transitions shift toward the high-T side, as revealed in our simulation.
On the other hand, an obvious separation between the nematic order parameter Ψ and the AFM order parameter Φ(π, 0) can be observed in the J Bi = 0.003 eV case. In the whole investigated T range, Ψ stays well above Φ(π, 0), suggesting that the nematic transition occurs at a temperature higher than T N . In addition, the difference in the positions of the two peaks of the susceptibility curves clearly shows the separation between the nematic and (π, 0) AFM phase transitions (T N = 118.4 K and T* = 123.4 K estimated from the peaks of χ Φ(π, 0) and χ Ψ ). As pointed out in earlier reports, the nematic phase has a broken Z 2 symmetry with short range NN spin-spin correlations [35]. In this report, it is demonstrated that the Z 2 symmetry can be broken by the additional positive biquadratic coupling. In other words, our simulation suggests that such an interaction may be another origin for the stabilization of the nematic order in pnictides [38], in addition to the couplings to the lattice degrees of freedom [23,39].
Other values of J Bi are also simulated. As a summary, the simulated phase diagram for positive biquadratic coupling (J Bi > 0) is presented in figure 2. For nonzero positive biquadratic J Bi , two phase transitions occur: the first transition from the paramagnetic (PM) phase to the nematic phase occurs at T*, and then a successive transition to the (π, 0) AFM order occurs at T N . Both transitions shift toward the high-T side with increasing J Bi , indicating that both the long range (π, 0) AFM order and the nematic order are favored by such an interaction, as discussed above. The interesting region in the phase diagram (with yellow background) is just the middle T range between T* and T N in which the nematicity is ordered but the global magnetism remains disordered. Thus, the experimentally reported nematic transition above T N for BaFe 2 As 2 can be qualitatively reproduced when the additional positive biquadratic coupling is introduced in the SF model.
Since an anisotropic J NN is adopted in the above simulation, which in principle breaks the tetragonal-type rotation symmetry, it is necessary to check its contribution to nematic order, which also breaks the symmetry between the x-axis and y-axis. Thus, the impact of the anisotropy of the superexchange interaction J NN has been systematically investigated in order to further verify this point. Figure 3 shows the simulated phase diagram in the (A J , T) parameter space for J Bi = 0.003 eV. Here, A J = J J NN x NN y −1 is defined to characterize the degree of anisotropy of the superexchange interaction. The separation between T* and T N persists even when A J is decreased to small values like 0.05. Therefore, the role of A J in our simulation is not to stabilize the nematic order but to provide a bias field to choose the preferred axis for  already-formed nematic order. A similar method has been used in earlier simulations to choose an axis in systems with degenerated axes [40]. Such a bias field is not mandatory in principle but practically helpful in simulations to avoid multi-axes of nematic domains. In addition, the simulation results shown in figures 1(a), (b) have definitely confirmed that the anisotropic J NN can not induce the nematic order above T N although both of them break the tetragonal-type rotation symmetry. In short, the formation of the nematic order is almost independent of the superexchange anisotropy.

Possible flux phase
It was predicted in earlier first-principles calculations that the biquadratic coupling might change sign in some other pnictides such as KFe 2 Se 2 [16]. Thus, the behavior of the SF model with the negative biquadratic interaction (J Bi < 0) is also investigated for integrity in this report. The simulated Φ(π, 0) and its susceptibility as a function of T for various negative values of J Bi are shown in figures 4(a), (b), respectively. For each J Bi , two peaks are observed in the χ Φ(π, 0) -T curve, indicating the successive two phase transitions with decreasing T. At the first transition point, spin structure factor Φ(π, 0) increases from the background baseline, while Φ(π, 0) remains small (figure 4(c)), indicating a standard transition to the (π, 0) AFM phase. When T falls down to the second transition point (T F ), Φ(0, π) steeply increases, while the value of Φ(π, 0) increases gradually. The MC snapshot (inset of figure 4(a)) at J Bi = −0.002 eV and T = 23 K shows the so-called flux state. It is noted that this state has four sites in one unit cell, which is clearly depicted in the snapshot. For J Bi < 0, NN spins tend to be perpendicular to each other to satisfy the negative biquadratic interaction, leading to the stabilization of the flux state. In fact, the flux transition point T F can also be estimated from the position of the peak in the susceptibility of Φ(0, π), as shown in figure 4(d).
For the ideal flux state, the value of Φ(0, π) should be identical to Φ(π, 0). However, due to the anisotropic J NN used in our simulation, they are not exactly equivalent. The anisotropic NN exchange interaction favors the (π, 0) AFM order more than the (0, π) AFM one. Thus, there is strong competition between the anisotropic exchange interaction and the negative biquadratic coupling, leading to the fact that Φ(π, 0) is larger than Φ(0, π), as revealed in our simulations. However, the real space pattern as shown in the inset of figure 4(a) is unambiguous and it is clear that the obtained state is the flux order.
As expected, the flux state can be strongly stabilized by the negative biquadratic interaction, and the transition point significantly shifts towards the high T side with the increasing absolute value of negative J Bi . As a result, the flux state occupies a considerable region in the parameter space for negative values of J Bi , as is clearly shown in the simulated phase diagram (figure 5).
In fact, the possible flux state has been proposed in some earlier theoretical works to explain the low value of magnetization in several pnictides [41][42][43]. In addition, the phase diagrams of the effective Heisenberg model obtained by mean-field approximations and MC simulations also show a rather large region with the flux state for negative biquadratic couplings, similar to that presented in this work [16]. However, it is observed that the anisotropy of the spin fluctuation spectrum increases with increasing T, which is inconsistent with the experimental report, demonstrating the limitation of such a model [44]. Thus, our work has verified the possibility of this flux magnetic ordering when a weak negative biquadratic interaction (J Bi = ∼J NN /10) is introduced into the SF model. In fact, this flux state has been predicted in the double-exchange models for manganites in earlier reports [28,45], and suggested to be associated with the anomalous Hall effect [29]. Thus, more interesting effects may be caused by the flux state in the SF model, which remains to be checked. However, this issue is beyond the scope of the present work.

Conclusion
In summary, the effect of the biquadratic interaction on the magnetic properties in iron-based superconductors has been investigated based on the three-orbital spin-fermion model. It is demonstrated that the experimentally reported nematic state can be developed when the additional positive biquadratic interaction is considered. Thus, our study suggests that such interaction may be another origin for the nematic order, in addition to the couplings to the lattice degrees of freedom that have been identified in earlier reports. In addition, the simulation shows that the so-called flux state can be stabilized by a rather weak negative biquadratic interaction, which goes beyond previous theoretical works based on the classical spin model.