Effects of intercomponent couplings on the appearance of time-reversal symmetry breaking fermion quadrupling state in two-component London models

A detection of bosonic metallic state that breaks the $Z_2$ time-reversal symmetry has been recently reported in Ba$_{1-x}$K$_x$Fe$_2$As$_2$ with a doping level $x \approx 0.8$. This is a metallic state of fermionic quadruplets that breaks time-reversal symmetry. As such, it has no condensate of Cooper pairs but has a long-range order between fermionic quartets. In the present manuscript, we investigate the emergence of this phase in a two-component London model via Monte Carlo simulations as a function of various intercomponent couplings.


I. INTRODUCTION
Superconductivity arises as a consequence of the formation and condensation of electron pairs. The condensation of Cooper pairs in many materials is well described by the mean-field Bardeen-Cooper-Schrieffer (BCS) theory [1,2].
However, multicomponent systems may exhibit a different kind of ordering associated with the long-range order of electron quadruplets [3][4][5][6][7][8][9][10][11][12][13][14]. Related phases have also been discussed in strongly-correlated ultracold atoms [15][16][17] and other models [18,19]. This type of order arises when the energy cost of composite topological defects (i.e. bound states of defects in different fields) is significantly cheaper than the energy cost associated with elementary topological defects. In this case, thermal or quantum fluctuations can induce the proliferation of composite defects that break the phase coherence of Cooper pairs, without breaking the phase coherence of fermionic quartets. Thus leading to the appearance of a quartic fermionic state.
The recent experiment [20] reported the observation of a fermionic-quadrupling condensate that spontaneously breaks time-reversal symmetry in Ba 1−x K x Fe 2 As 2 .
The multiple broken symmetries required for the formation of this type of fermionic-quadrupling order can result from the presence of multiple bands crossing the Fermi surface that, in turn, can lead to a ground state where the inter-band phase difference is neither 0 nor π [21][22][23][24][25][26][27][28][29][30]. Such a ground state, has an additional Z 2 two-fold degeneracy corresponding to the spontaneous breaking of the time-reversal symmetry (BTRS). In this scenario, the superconducting (SC) ground state spontaneously breaks a U (1) × Z 2 total symmetry.
The muon-spin rotation experiments [27,31] have found evidence for such a multicomponent s + is superconducting state at low temperatures in Ba 1−x K x Fe 2 As 2 with a doping level x ≈ 0.8. The main evidence in favor of this state comes from the analysis of polarization of spontaneous magnetic fields [31][32][33], however the closely related s + id state [34] cannot be completely ruled out today.
At the level of mean-field BCS theory, the critical temperature associated with the spontaneous breakdown of the time-reversal symmetry is always smaller than or equal to the superconducting temperature: T c ≥ T Z2 c . [22][23][24]26]. The initial experimental results shown in [27,31] for a certain range of doping, were consistent with this picture. However, very recent experimental results [20] found in Ba 1−x K x Fe 2 As 2 a regime near the doping level x ≈ 0.8 where T Z2 c > T c , signaling the onset of a fermionic-quadrupling phase with spontaneously broken time-reversal symmetry. The conclusion was reached based on a number of experimental probes: muon-spin rotation experiments, conductivity, diamagnetic response, thermolelectric and ultrasound probes. These probes revealed a very interesting state of matter that required further exploration. For a recent theoretical work on the effective model, and some of the properties of this state see [35].
The previous theoretical studies beyond mean-field approximation [6,7,20] indicate that the formation and size of that state depend on the strength of direct intercomponent coupling, the magnetic field penetration length relative to other length scales, and on the strength of mixed-gradient terms. However, systematic investigation of the interplay of these parameters was not performed.
Here we present such an investigation.
The paper is organized as follows. In Sec. II, we introduce the continuum two-component London model. In Sec. III we discuss the Villain approximation and the details of the Monte Carlo (MC) simulations. In Sec. IV, we present the numerical results obtained in the different regimes investigated. Conclusions are given in Sec. V.

II. THE MODEL
As discussed in detail in [36], the free-energy density of a clean three-band superconductor with a bilinear intraband Josephson interaction, can, under certain conditions, be approximated by a model with only two components and a biquadratic Josephson coupling. Similar model arises for s + is state generated by impurities in dirty two-band superconductors [37]. Finally, except for certain differences in structure of gradient terms, similar models also result from other systems such as s + id [34], p + ip [12,38,39], s + ig [28] and systems that break time-reversal symmetry.
In this work, we will study a two-component London model in three spatial dimensions, as a model for U (1) × Z 2 superconductors in the type-II limit.
The free-energy density of the system reads: Here φ i ∈ [0, 2π) are phase, coupled by the gauge field A. The two condensates have the same electric charge: e 1 = e 2 = e. The model (1) accounts for three different intercomponent interactions: the coupling of the two charged condensates via the fluctuating gauge field A, the interaction via the coupling constant ν, and the second-order biquadratic Josephson interaction with coupling constant η 2 . The coefficient ν sets the intercomponent current-current coupling. Such dissipationless-drag interaction (i.e. mixed-gradient) terms are generically present in multicomponent systems and originate for example from Fermi-liquid corrections or strong correlations as shown in various physical contexts [16,18,[40][41][42][43][44][45]. From Eq. (3), one can easily derive the stability condition of the system, i.e. the condition ensuring the free energy to be bounded from below, being: ν < √ ρ 1 ρ 2 .
Without loss of generality, in what follows we will fix ρ 1 = 1, tuning the disparity of the components via the parameter α = ρ 2 /ρ 1 . This corresponds to a reduction of the free parameters of the model by a proper rescaling of the coupling constants, the gauge field, the electric charge and the free energy.
In the absence of the Josephson interaction, the model (1) has a U (1)×U (1) symmetry. Phase transitions in various regimes in models with U (1) × U (1) symmetry were considered as a function of magnetic-field penetration length and strength of mixed-gradient terms in [3-5, 8-10, 15-18, 20, 46-48]. When such symmetry is present, the phase transition is driven by the proliferation of vortex loops in the two-component condensate (unless the system is strongly type-I). These topological excitations can be denoted by a pair of two integers corresponding to the phase winding in each condensate: Together with single-component vortices, which have phase windings in only one of the two condensates being (1, 0) or (0, 1), the system can also proliferate composite vortices (1, 1) or (1, −1) leading to phases with composite order.
Various phases supported by the model can be seen more clearly by rewriting Eq.(1) in terms of charged and neutral modes [10,35,46,49,50]: where: Consider first the case η 2 = 0. The most detailed Monte-Carlo calculations of the resulting phase diagram in three dimensions was presented in [10]. As one can see from Eq. (3), the gauge field is coupled with the phase-sum mode. As a consequence, for e = 0 the energetic cost per unit length of a (1, 1) composite vortex is finite and can be made arbitrarily small by increasing the value of the electric charge. When (1, 1) vortices proliferate the system retains order in the phase difference. That state is not a superconductor and has nontrivial magnetic properties, described beyond the London limit by an effective model related to the Skyrme model [35]. The energy cost of (1, −1) in turn depends on the parameter ν. When (1, −1) vortices proliferate, but (1, 0) or (0, 1) do not proliferate, the system retains order in the phase sum, representing a charge-4e superconductor.
Let us now consider the case where η 2 = 0 and ν > 0. The presence of an intercomponent biquadratic Josephson interaction changes the phase diagram of the model by explicitly breaking the U (1) × U (1) symmetry down to a U (1) × Z 2 symmetry, where the U (1) symmetry is associated with the charged phase-sum mode.
For η 2 < 0 the phase difference φ 1,2 = φ 1 − φ 2 in the ground state can be either 0 or π; while for η 2 > 0, either π/2 or −π/2. In this work we will focus on the latter case, where the spontaneous breakdown of the Z 2 symmetry is associated with the spontaneous breaking of a time-reversal symmetry, since the complex conjugation of the order parameter leads to a different ground state.
The presence of a Josephson coupling also brings in a different type of topological excitations. These appear as domain walls separating in space two energetically equivalent states, namely a state where φ 1,2 = π/2 from another where φ 1,2 = −π/2. So, just as the proliferation of composite-vortices of the kind (1, 1) will restore the U (1) continuum symmetry, the proliferation of domain walls will restore the Z 2 time-reversal symmetry. The phase diagram of the system will thus depend on the relative energetic cost of these two topological defects. In particular, a scenario where T Z2 c > T c can arise if the energetic cost of domain-walls nucleation is sufficiently high with respect that of (1, 1) vortices. The parameters of the model, i.e. the gauge field coupling that parametrizes the magnetic field penetration length, the strength of the Josephson coupling, and the mixedgradient coupling, each affect the relative cost of vortex and domain-wall excitations. However, the peculiarity of the model is such that simple energy arguments cannot be used to map out its phase diagram. Firstly, entropic factors are important and nontrivial. Secondly, domain walls are strongly and nontrivially interacting with vortices and have, under certain conditions, a tendency to form composite objects: Skyrmions [33,51,52]. Hence we use Monte-Carlo approach to study the phase diagram of the system.

III. DETAILS OF THE MONTE CARLO SIMULATIONS
A. Villain lattice model To perform Monte Carlo simulations, we need to provide a discrete lattice representation of the continuum model Eq.(1). We consider a three-dimensional cubic lattice of size L 3 and lattice spacing h = 1. The phases are defined on the lattice vertices φ r,j with j = 1, 2 labeling the two components. On the other hand, both the phase gradient, which is defined as the phase difference between two neighbouring sites ∂ µ φ(r) j → ∆ µ φ r,j , and the gauge field A r,µ are associated with the link connecting the vertex r with its neighbour r + µ, being µ =x,ŷ. The lattice curl of the gauge field, defined around a unitary plaquette, reads (∇ × A(r)) → ( κ,η µκη ∆ κ A r,η ), being µκη the Levi-Civita symbol.
We have performed Monte Carlo (MC) simulations of the Villain Hamiltonian Eq. (5), locally updating the two phase fields φ 1 , φ 2 ∈ [0, 2π) as well as the gauge field A by means of the Metropolis-Hastings algorithm. A single MC step here consists of the Metropolis sweeps of the whole lattice fields. To speed-up the thermalization, we also implemented a parallel tempering algorithm, allowing swap of field configurations between neighbouring temperatures. Typically, we propose one set of swap after 32 MC steps. For most of the numerical simulations, we performed a total of 2×10 5 Monte Carlo steps, discarding the transient time occurring within the first 50000 steps. For the simulation performed in the limit of large Josephson coupling we implemented a cluster update to prevent the system from getting stuck in metastable states and we extended the total MC time up to 4 × 10 5 steps discarding the first 150000 steps. We have considered different values of the linear size L, as is needed to properly assess the critical points of the model. In the limiting case e = 0, the U (1) transition is associated with the onset of a superfluid phase, captured by the helicity modulus of the phase sum [10,17]. The helicity modulus Υ measures the energetic cost associated with an infinitesimal twist of the order parameter phase across the system. In a multicomponent system, one can define several helicity moduli corresponding to different linear combinations of individual phase twist. In the twocomponent case, for each choice of the coefficients {a i } in one can define a corresponding helicity modulus where and . . . stays for the thermal average over the MC steps. Since the U (1) superfluid phase transition is associated with the phase-sum mode, the relevant observable is the phase-sum helicity modulus Υ + defined by the choice a 1 = a 2 = 1. We computed Υ + for different values of the linear sizes L and determined the critical temperature T c by means of the finite-size crossings, extrapolated to the thermodynamic limit, of the quantity LΥ + . In the case where e = 0, we locates the U (1) superconducting transition by computing the dual stiffness [20,[54][55][56], that accounts for the onset of the SC state by measuring the Meissner effect: We compute the dual stiffness in the z direction at the smallest relevant wave vector in the x direction, i.e. q x min = (2π/L, 0, 0). In what follows, we denote ρ z (q x min ) simply as ρ. This observable accounts for the longrange fluctuations of the magnetic field, suppressed in the superconducting phase and finite in the normal phase. Thus, contrarily to the superfluid stiffness, it is expected to be zero in the superconducting phase and non-zero in the normal one. Finally, we use finite-size crossings of Lρ, extrapolated to the thermodynamic limit, in order to locate superconducting transitions.

The Z2 phase transition
We define a Z 2 Ising order parameter m associated with the two possible degeneracy of the ground state. In particular, we set m to be equal to +1 or −1 according with the sign of the phase difference φ 1,2 .
In order to locate the Z 2 critical temperature, we compute the Binder cumulant U [57,58] for the order parameter m: which is expected to be a universal quantity at the critical point.
The critical temperature T Z2 c , associated with the spontaneous breaking of the time-reversal symmetry, is thus determined by means of finite-size crossings of the Binder cumulant U extrapolated to the thermodynamic limit.
The error bars of all the observables are estimated via a bootstrap resampling method. In the figures shown, when not visible, the estimated error bars are smaller than the symbol sizes.

IV. RESULTS
Here we present the Monte Carlo numerical results obtained for different regimes of the two-component London model (1).
Let us start discussing the limiting case of infinite penetration length λ → ∞, i.e. e = 0. In the isotropic case where α = 1, the free energy (1) reads: where we fixed η 2 = 0.1. The increase of the drag coupling ν reduces the energetic cost of nucleating (1, 1) composite vortices, whose proliferation restores the superconducting U (1) symmetry, with respect to that domain walls, associated with the the Z 2 time-reversal symmetry.
The two inverse critical temperatures, respectively β c (U (1)) and β c (Z 2 ) are reported in Fig.1(a) as function of the coupling ν. For low values of ν, the simulations indicate a single phase transition dividing a low-temperature superconducting phase, where both the U (1) and the Z 2 time-reversal symmetry are spontaneously broken, from a high-temperature phase where the system is in its normal metal phase. On the other hand, for sufficiently large ν, the two phase transitions split apart with β c (U (1)) > β c (Z 2 ). In this regime, the system shows a fermionic-quandrupling phase arising between the low-temperature superconducting phase and the high-temperature normal phase. Such a state is disordered in the phase-sum mode, i.e. it is not superconducting, but still ordered in the phase-difference mode due to the spontaneous breaking of the Z 2 timereversal symmetry. The order parameter characterizing this phase is thus the two-component phase difference, fourth order in terms of fermionic fields being: The region in the phase diagram Fig.1(a) where this BTRS quartic metal phase arises becomes wider for higher values of ν.
In the present case the free energy reads: with ζ = 1 + α − 2ν. Since here we are studying the model for fixed values of ν as function of α it is convenient to rewrite the stability condition of Eq. (14) as: α > ν, which automatically fulfill α > ν 2 , being ν < 1. The Monte Carlo numerical results for the model (14) are shown in Fig.2(a)-(b) for ν = 0.4 and for ν = 0.6 respectively. In both cases, decreasing α, the two transitions move to lower critical temperatures. However, for ν = 0.4 ( Fig.2(a)), we find that the inverse critical temperature associated with the time-reversal symmetry breaking β c (Z 2 ) grows faster than the one associated to the U (1) superconducting symmetry β c (U (1)), with a resulting split of the two transitions for α < 1. In this case, the split occurs with β c (Z 2 ) > β c (U (1)), implying a relative increase of the time-reversal invariant superconducting state. On the other hand, for the case ν = 0.6 ( Fig.2(b)), we find the opposite scenario: by decreasing α the inverse critical temperature β c (U (1)) increases faster than β c (Z 2 ). After presenting the numerical results obtained for the extreme type-II limit, let us now consider the phase diagram as a function of magnetic field penetration length relative to other length scales. We consider two equivalent superconducting components ρ 1 = ρ 2 , i.e. α = 1.
The free energy of the system reads: We fixed the value of the Josephson coupling η 2 = 0.1 and the dissipationless drag interaction to ν = 0.2, studying Eq. (15) for different values of the electric charge e. The resulting phase diagram is shown in Fig.3 (a). For small values of e, the phase diagram for ν = 0.2 does not qualitatively change with respect to the neutral case: we can resolve only one phase transition associated with the spontaneous breaking of the total U (1) × Z 2 symmetry. However, by further increasing e, the energy of (1, 1) vortex-loops decreases. This reduces the critical temperature associated with the U (1) symmetry breaking, without equally affecting the Z 2 order. As a result, the two phase transitions separate for large enough e, with the appearance of the BTRS quartic metallic phase.
Consistently with the results obtained for ν = 0 in a different discretization scheme of the London model [7], we obtain that for strong enough coupling with the gauge field, the splitting of the two phase transitions appears in this model even for ν = 0 Fig.4. We finally consider the limiting case of very large Josephson coupling η 2 at fixed values of the electric charge e and the dissipationless drag interaction ν. As η 2 increases, the energy cost of nucleating a domain wall increases, with a resulting increase of the critical temperature associated with the Z 2 time-reversal symmetry breaking. However, the Josephson coupling is a local coupling, hence for very large values of η 2 the energy cost of a domain wall will saturate since the domain-wall width cannot be larger than the lattice spacing. Thus one would expect the saturation of T c (Z 2 ) in the lattice London model. As reported in Fig.5, for the model (1) with ν = 0 and e = 0 in our numerical results we cannot resolve the splitting for the Z 2 and U (1) phase transitions, at least for the system sizes that we simulated.
Nevertheless, large values of the Josephson coupling η 2 widen the region of the phase diagram Fig.1(a) where the BTRS quartic metal phase appears. For ν = 0.4, as reported in Fig.6, the two phase transitions separate significantly for η 2 = 1.

V. CONCLUSIONS
In conclusion, the pairwise interaction of electrons usually leads to the formation of pair condensates. However, in systems that break multiple symmetries, there could be fluctuations-induced states corresponding to fermionic-quadrupling condensates. These states appear from the proliferation of bound states of topological defects that partially restore symmetry and lead to an effective fermionic-quadrupling interaction.  1)). In such temperature range between β c (U (1)) and β c (Z 2 ) there is a BTRS quartic metallic phase.
We studied the appearance of a fermionic-quadrupling condensate that breaks time-reversal symmetry, as reported in a recent experimental work [20] in Ba 1−x K x Fe 2 As 2 with doping level x ≈ 0.8. The phase shows order only in the phase differences between the components of complex fields. While the precise microscopic model for this compound is not established yet, we studied the appearance of such a state in a threedimensional two-component London model. We presented the phase diagram of the model and appearance of the quartic state as a function of the magnetic field penetration length (parameterized by the gauge field coupling constant), the strength of biquadratic Josephson, and mixed-gradient couplings.