Unusual $H$-$T$ phase diagram of CeRh$_2$As$_2$ -- the role of staggered non-centrosymmetricity

Superconductivity in a crystalline lattice without inversion is subject to complex spin-orbit-coupling effects, which can lead to mixed-parity pairing and an unusual magnetic response. In this study, the properties of a layered superconductor with alternating Rashba spin-orbit coupling in the stacking of layers, hence (globally) possessing a center of inversion, is analyzed in an applied magnetic field, using a generalized Ginzburg-Landau model. The superconducting order parameter consists of an even- and an odd-parity pairing component which exchange their roles as dominant pairing channel upon increasing the magnetic field. This leads to an unusual kink feature in the upper critical field and a first-order phase transition within the mixed phase. We investigate various signatures of this internal phase transition. The physics we discuss here could explain the recently found $H$--$T$ phase diagram of the heavy Fermion superconductor CeRh$_2$As$_2$.


I. INTRODUCTION
The discovery of the heavy Fermion superconductor CePt 3 Si [1] triggered much interest in superconductivity in materials lacking an inversion center in their crystalline structure. Among the most intriguing features of such non-centrosymmetric superconductors we may mention the parity mixing of Cooper-pair states [2][3][4], the helical superconducting phase induced by a magnetic field [1,5,6], and the appearance of topological superconducting phases [7,8]. These properties emerge most spectacularly, if the superconductor is intrinsically prone to unconventional pairing as expected for heavy-Fermion compounds due to strong electron correlation. Parity mixing is generated by anti-symmetric spinorbit coupling present in non-centrosymmetric crystals, with Rashba spin-orbit coupling in systems which lack mirror symmetry the best-known example. This type of spin-orbit coupling is realized in a number of Ce-based heavy Fermion superconductors with tetragonal crystal symmetry, such as CeIrSi 3 [9] and CeRhSi 3 [10], as well as CePt 3 Si. The probably most striking properties seen in both, CeIrSi 3 and CeRhSi 3 is the enormous upper critical field for field directions perpendicular to the basal plane (lacking mirror symmetry), exceeding the paramagnetic limit by far [11,12]. This indicates also an extremely short coherence length which is not unusual for heavy-Fermion superconductors. The isostructural La-versions of these two compounds (replacing Ce by La), which superconduct with comparable critical temperature but without heavy-Fermion properties, have low upper critical fields dominated by orbital depairing [4]. The interplay of orbital and paramagnetic depairing has played an important role in other Ce-based superconductors, too, most notably in the centrosymmetric CeCoIn 5 , which besides displaying paramagnetic limiting behavior is also well-known for its low-temperature high-magnetic field phase, the so-called Q-phase, where superconductivity coexists with a spin density wave state [13,14]. This in-dicates that the Maki parameter α M = √ 2H c2 (0)/H p (0) is not small in this material, where H c2 and H p denote the paramagnetic and orbital upper critical field, respectively, [15].
Recently, it has been noticed that centrosymmetric materials may incorporate locally non-centrosymmetric units which can influence properties of a superconductor profoundly, in particular their behavior in a magnetic field [16,17]. Features of this origin are found, for example, in superlattices such as the regular stacks of superconducting CeCoIn 5 alternating with layers of YbCoIn 5 [18,19]. Here, local non-centrosymmetricity at the interfaces between Ce-and Yb-layers involves antisymmetric spin-orbit coupling which apparently protects Cooper pairs against paramagnetic depairing in this highα M system [20]. Moreover, it was suggested that a peculiar parity-mixing present in these superconductors [16] may give rise to a magnetic field phase transition changing the character of the pairing state [21].
The very recently discovered heavy-Fermion superconductor CeRh 2 As 2 belongs also to the class of locally noncentrosymmetric superconductors. While its tetragonal crystal structure has an inversion center, it consists of layers with alternating violation of inversion symmetry. The upper critical field directed along the c-axis (perpendicular to the staggered layers) extrapolates to ∼ 10 T at zero temperatures, which lies far beyond the paramagnetic limiting field H p ∼ 0.5 T for a critical temperature T c 0.4k [22]. The most intriguing aspect of this superconductor so far is its behavior in the c-axis magnetic field. The upper critical field shows a pronounced kink for a temperature T roughly half of T c . This anomaly strongly hints towards a switch in the order parameter symmetry upon increasing magnetic field.
Stimulated by this experimental finding, we investigate here a possible scenario explaining this behavior based on a model of a material with staggered Rashba spin-orbit coupling [16]. This symmetry property can be easily implemented into a generalized Ginzburg-Landau theory which is well suited to discuss the superconducting order parameter in the mixed phase. The basic idea of a field-induced phase transitions has been previously explored by Yanase and co-workers in the context of the Ce/YbCoIn 5 superlattices using a Bogolyubovde Gennes approach which ignored the presence of vortices [20,21,23]. In these studies, a few superconducting layers were considered whose two boundary layers constitute opposite non-centrosymmetric environments and, thus, Rashba spin-orbit coupling of opposite sign. The subsequent extension to the mixed phase was performed by Möckli et al. using a Ginzburg-Landau model [24]. In our investigation we base our analysis of the infinitelayer system with staggered Rashba spin-orbit coupling on this approach to analyze the situation of CeRh 2 As 2 .
In the following, we first introduce the Ginzburg-Landau model describing the effect of parity mixing by including an even-and odd-parity order parameter component. We show that this allows us to discuss two phases, which we refer to as A and B phase, which, in view of the crystalline inversion center, can be considered as even or odd, respectively. In a magnetic field along the crystalline c axis, we can reproduce qualitatively the kink feature of the upper critical field and also find an internal first-order phase transition within the mixed phase. The mixed phase is described here in a scheme using a cellular approximation of the vortex lattice unit cell [25]. Finally, we propose several ways to detect the internal phase transition in experimental.

II. GINZBURG-LANDAU MODEL FREE ENERGY
A. Free energy functional and pairing symmetries Even though the crystal lattice of our system possesses global inversion symmetry and a superconducting order parameter can thus be classified as even or odd under inversion, the unit cell of the superlattice includes two locally non-centrosymmetric subsystems stacked alternatingly on top of each other. The lack of local inversion symmetry gives rise to a staggered Rashba spin-orbit coupling, leading to the formation of Cooper pairs of mixed parity within each layer [3,16,17]. Within the Ginzburg-Landau formalism, this is accounted for by introducing even-and odd-parity order-parameter components on each layer j, Ψ e,j (r, ϕ e,j ) = ψ e,j (r) exp(iϕ e,j ) and Ψ o,j (r, ϕ o,j ) = ψ o,j (r) exp(iϕ o,j ), respectively, with r the in-plane coordinate.
We write the free energy functional as with A denoting the vector potential. The free energy density of the jth layer takes the form, with B = ∇×A the magnetic induction perpendicular to the plane and the summation running over both the even (l = e) and odd (l = o) pairing channels. The parameters a l = a 0,l (T c,l − T ), with T c,l representing the bare critical temperature of the respective pairing channel, and b l are phenomenological constants. We choose both critical temperature T c,l ≥ 0, indicating that the corresponding Cooper pairing channels have comparable attractive interactions, but with T c,e > T c,o . The parameter m l represents the in-plane Cooper-pair mass of the respective order-parameter component and the covariant derivative D = (−i ∇ + 2eA(r)/c) is restricted to the two inplane coordinates. The last term in the sum incorporates the coupling of the order parameters between neighboring layers with coupling strengthJ l , assuming a quasi-twodimensional electronic structure. Furthermore, we included the paramagnetic depairing effect, which directly affects the even-parity component and whose strength is determined by the parameterQ [4,21,26]. Note, however, that due to the parity mixing f (j) eo discussed below, paramagnetic depairing is effectively detrimental to both order-parameter components. We assume the field to be equal for all layers.
To lowest order, the spin-orbit-coupling-induced parity mixing between even and odd components in the jth layer takes the form [27] with˜ j the coupling strength on the jth layer, which, due to the staggered nature of the spin-orbit coupling, possesses alternating sign on neighboring layers, i.e.˜ j = (−1) j˜ with˜ > 0. The free energy is constructed in a way to be invariant under both time-reversal symmetry as well as local U(1) gauge transformations Φ At the onset of superconductivity, both order parameters become non-zero, and the resulting state is invariant under time-reversal operation. Furthermore, the relative phase of the two order parameter components is gaugeinvariant and anti-symmetric under time-reversal operation T , such that [27] (ϕ e,j − ϕ o,j ) mod 2π = 0, π.
To minimizes the interlayer coupling, the dominant order-parameter component maintains its sign across all layers, whereas the subdominant alternates in sign on neighboring layers [20]. Thus, the allowed orderparameter structures are either globally even under inversion, Fig. 1). The global parity is established by an inversion center lying between two layers and the corresponding parity operation exchanges two layers as well as changes the sign of the odd-parity order parameter component. To avoid confusion, we refer to the (globally) even order-parameter structure, which is favored at low fields, as phase A and the odd, high-field phase as phase B. Finally, the phase difference ϕ e,j − ϕ o,j adapts on each layer such that the overall free energy is minimal for the given order parameter structure, in particular, In order to limit the number of free parameters, we assume several coefficient in the free energy to be the same for both order-parameter components, in with the Ginzburg-Landau coherence length ξ given as this is equivalent to re-scaling the vector potential by Φ 0 /[2πξ(T )], where Φ 0 = hc/2e represents the flux quantum. In this formulation, the superconducting flux quantum is given by φ 0 = 2π/κ 0 and the dimensionless upper critical field B c2 (orbital depairing field) is equal to κ 0 , with κ 0 = λ/ξ the temperature-independent Ginzburg-Landau parameter. Lastly, we normalize the orderparameter components by a e (T )/b and denote them in dimensionless units by Ξ l,j (ρ, ϕ l,j ) = η l,j (ρ) exp(iϕ l,j ). In summary, we use The dimensionless free energy functional F is related to the free energy F through We rewrite the free energy density of Eq. (2) in dimensionless form, and split it up into the following five parts: where the dimensionless basic free energy is given as The Ginzburg-Landau parameter κ 0 is equal for both even and odd components under the assumptions mentioned above. The term quadratic in the triplet component retains an explicit temperature dependence as the order parameter is normalized by the singlet quantity a e = a 0,e (T c,e − T ). Next, we write the free energy density of the magnetic orbital part as where B(ρ) = ∇×A(ρ) is the dimensionless magnetic induction and χ n represents the normal-state susceptibility. Note that in the last line we exploited the constant phase difference, Eq. (5), leading to ∇ϕ e,j = ∇ϕ o,j and performed the gauge transformation A → A + ∇ϕ/κ 0 [24]. Finally, the interlayer coupling energy f with J, and Q corresponding to the coupling constants J,˜ ,Q in the dimensionless formulation. The interlayer coupling parameter has the unit of an energy (per area) and is therefore re-scaled in the following way: and thus J =J/(a e λ 2 ) = 2J/( κ 0 ) 2 . For the same reasoning, the spin-orbit coupling strength transforms as = 2˜ /( κ 0 ) 2 . Lastly, the phenomenological constant Q is related toQ through Q(T ) =Qa e /b, resulting in an explicit temperature dependence of f (j) p in the dimensionless formulation. For simplicity we denote Q(0) = Q.
In our formulation, the Maki parameter can be expressed through the system parameters, the normal state susceptibility χ n and Ginzburg-Landau parameter Imposing the order-parameter structure for the A and B phases, one can directly evaluate the sum over the layer index j. Note that is the only term in the free energy that includes the order parameter phase difference between the layers. In the remainder of the paper we choose the parameter values κ 0 = 100, α M = 20, = J = 1, a 0,e = 1, a 0,o = 0.4, b = 1, T c,o /T c,e = 0.6, and Q = 0.1 · 4πa e /b unless otherwise stated.

B. Superconducting instability -linearized GL equations
Before accounting for the spatial modulations of the order parameter due to the flux-line lattice, we solve the linearized GL equations to obtain the onset of superconductivity in a magnetic field. For phase A, these equations read and for phase B For B = (0, 0, B), we choose the gauge A = (0, Bx, 0) and use the ansatz Ξ = e ikyy u Ξ (x), such that we can solve giving and leading to the Landau levels E n = B(2n + 1)/κ 0 . The lowest Landau level is E 0 = B/κ 0 , and thus, the critical field for the A phase is determined by solving the equation, An analogous equation is obtained for the B phase. The resulting B-T phase diagram is shown in Fig. 2 and agrees well with the upper critical field of the model developed in Sec. III.

C. Internal phase transition
As paramagnetic limiting only directly affects the singlet component, the low-field phase A is paramagnetically limited. On the other hand, the high-field phase B is dominated by the triplet component, which is relatively immune to paramagnetic pair breaking. As we consider systems with a large Maki parameter, we expect that orbital depairing only plays a subordinate role when discussing the transition from phase A to B.
To gain an understanding of this internal phase transition, we thus use a simplification of Eq. (2) that ignores the in-plane spatial modulations of the order parameter as well as the orbital depairing. For this analysis the terms second order in the order parameters are sufficient. The free energy for both phases is then written as Expressing the quadratic form through a symmetric matrix and imposing its determinant to be zero, we obtain the critical field for both phases the larger of which determines the field where superconductivity disappears (see Fig. 3).
To determine the internal phase transition from phase A to phase B, their corresponding free energies are minimized and compared. As shown in App. A, this leads to the expression for the critical field of the internal phase transition B c,AB which for finite is larger than T c,e and converges to T c,e for → 0. Similarly, by equating (29) with (30) we obtain the temperature at which the two critical fields meet, Equation (33) indicates that the spin-orbit coupling also increases the temperature range in which the B phase can be observed. Furthermore, it shows the importance of the interplay between spin-orbit coupling and the interlayer coupling for the formation of the high-field B phase.

A. Circular cell method
To analyze the situation for a regular flux-line lattice, we employ the vortex cell method introduced by Brandt et al. [25] and previously used for a trilayer system [24]. For this purpose, we decompose the system in plane into a triangular lattice of hexagonal cells each containing one vortex in the center. Following the Wigner-Seitz idea, we approximate the cells by a circle of the same area A cell as the hexagon. The free energy restricted to a single cell with the corresponding boundary conditions then reads with N the number of layers. This expression can be converted into a dimensionless formulation as done before. The cell approach is based on the trial function for the vortex structure [28], where l = e, o and η e,∞ , η o,∞ , and ξ c are treated as variational parameters for the even and odd bulk magnitude as well as the vortex core size. The trial functions, Eq. (35), reduce the system to a single Ginzburg-Landau equation, obtained by the variational derivative with respect to A, which reads where we defined By taking the curl of Eq. (36), one obtains the differential equation [29] 1 1 + 4πχ n 1 ρ d dρ whose general solution is given in Ref. [25], Here, I n and K n are the nth-order modified Bessel functions of first and second kind, respectively, and we define P 2 = ρ 2 B + ξ 2 ∆ with ρ B being the radius of the vortex cell, i.e., A cell = πρ 2 B = 2π/(κ 0B ) withB the mean magnetic induction in a vortex core. The coefficients c 1 and c 2 are determined through the boundary conditions of vanishing current density at the cell boundary as well as the flux threading the unit cell equalling a flux quantum φ 0 (for χ n = 0). The latter is formalized by integrating the left side of Eq. (39), which together with the first boundary condition leads to Using χ n 1, the coefficients c 1 and c 2 reduce to With this analytic expression for the magnetic induction, we can use the result by Hao and Clem for the magnetic part, F m =BB(0) [30]. Next, we integrate the free energy density within the cell which leads to the analytical expressions for the different parts of the free energy, and we introduced (50) In the following, we minimize the the free energy numerically using the Nelder-Mead algorithm with respect to the variational parameters η e,∞ , η o,∞ , and ξ c .

B. Results
We choose a large Ginzburg-Landau parameter of κ 0 = 100, ensuring that the lower critical field H c1 is very small compared to the upper critical field. The numerical results for the order-parameter magnitudes, as well as the vortex-core size as a function of the field (at zero temperature), are displayed in Fig. 4. The numerical results show that at low fields, the singlet-dominated phase A is indeed favored. This changes at sufficiently high fields, where the odd component becomes competitive. Figure 5 shows a typical B-T phase diagram, obtained by reinstating the explicit temperature dependence. The two outer phase boundaries are of second order, whereas the internal phase boundary is of first order. The inset shows the comparison of this critical field B c,AB with the one in Eq. (31) derived from the simplified free energy. The excellent agreement suggests that orbital depairing and the flux line lattice only play a minor role in determining the critical field for the internal phase transition, which is dominated by the paramagnetic limiting effect.
Although the expression of the critical field for the internal phase transition, Eq. (31), does not explicitly contain the spin-orbit or Josephson interlayer coupling constants and J, they play a crucial role in the formation of this high-field phase. This is illustrated in Fig. 6, where the minimal Maki parameter required for phase B to form decreases with increasing spin-orbit coupling strength. The opposite applies to J, as a strong interlayer coupling suppresses the admixed subdominant odd order parameter component in phase A, thereby diminishing the effect of the spin-orbit-induced parity mixing. A possible way to experimentally observe the A-B transition is by measuring the magnetization curve. To this end, we calculate the dimensionless magnetic field H using The analytic expression for H is included in App. B. The magnetization and magnetic susceptibility can then be extracted from 4πM = B − H and χ = M/H. A typical magnetization and susceptibility curve is shown in Fig. 7 showing a jump at the A-B transition and thus providing an experimental signature of the internal first-order phase transition.
We can approximately express the magnetization jump by (see App. B1) (52) Unlike the critical field for the internal phase transition, Eq. (52) explicitly contains the variational parameters ψ ∞ , η ∞ , and ξ c and can, thus, not be directly evaluated. We observe that the magnetization jump increases with growing J and . Consequently, increasing the interlayer coupling strength J and/or the spin-orbit coupling leads to a larger magnetization jump as the relative importance of the orbital depairing term shrinks. This behavior can be seen in Fig. 7, where we compared the analytical approximation of Eq. (52) with the full numerical solution.
Apart from the phase diagram and magnetization curve, a further experimentally accessible quantity is the latent heat. However, it must be noted that observing the A-B transition by changing the temperature at constant magnetic field may prove difficult for typical parameter values, as the critical fieldB c only shows a weak temperature dependence. To obtain the thermodynamic quantities, we must reintroduce the explicit temperature dependence. As expected for a first-order transition, the entropy S = −∂ T F is discontinuous at the internal phase transition, see Fig. 8.

IV. DISCUSSION
A previous study by Yoshida et al. has investigated non-centrosymmetric bi-and trilayer systems within the framework of Bogoliubov-de-Gennes [21]. They similarly find that at low fields, the sign of the singlet order parameter remains constant due to interlayer coupling, whereas at high fields, a phase forms in which the order parameter adapts to the sign of the staggered Rashba spin-orbit coupling. Furthermore, while the results for the trilayer system show certain additional peculiar features due to the centrosymmetric middle layer, the resulting phase diagram for the bilayer system is qualitatively similar to Fig. 5. However, the state they investigated neglected the in-plane spatial modulations of the order parameter in the mixed phase. Thus, important aspects of the internal phase transition, including the magnetization jump and the role of the Maki parameter, could not be thoroughly investigated.
The two anomalies at the internal phase boundary line, the jump in the magnetization and the entropy (latent heat), are not the only features which may be used to observe this phase transition. Additionally, we expect that anomalies are observable in the ac-susceptibility as well as in ultrasound velocity (a change in the elastic constants) for certain modes. A further intriguing feature could result from the fact that in the B phase, the odd-parity exceeds the even-parity component. It has been proposed that in this case helical subgap edge states would appear at the surfaces with in-plane normal vectors [31,32]. Such helical states would be visible in in-plane quasiparticle tunneling spectroscopy as zero-bias anomalies. In the A phase with dominant evenparity component, such helical edge states would be absent. Thus, we expect a clear qualitative change of the quasiparticle I-V -characteristics when passing through the transition. Finally, the experimental analysis of the magnetocaloric effect may allow to detect the internal phase transition [33].

V. CONCLUSION
Staggered non-centrosymmetric superconductors with Rashba spin-orbit coupling alternating in sign from layer to layer involve an order parameter of mixed parity which can be influenced by external magnetic fields. While an even-parity spin-singlet pairing state may dominate over an odd-parity spin-triplet state at low magnetic fields, paramagnetic limiting effects may lead to a suppression of the former in favor of the latter in a growing magnetic field. This may trigger a phase transition to a phase with dominant odd-parity component provided that the Maki parameter is sufficiently large, i.e. that orbital depairing is weak due to a very short coherence length.
We suggest that such a situation is realized in the recently found heavy-Fermion superconductor CeRh 2 As 2 , which shows an anomaly in the upper critical field which exceeds the expected paramagnetic limiting field drastically. We show that besides the kink-like anomaly of the upper critical fields to the normal state, there is also an internal phase transition within the mixed phase between the low-field phase A with dominant even-parity and the high-field phase B with dominant odd-parity pairing. This internal phase boundary may be observed by various means as detailed in the previous section. In particular, this phase boundary should connect to the kink of the upper critical field. An intriguing feature of the B phase is the fact that it may have subgap edge states observable in quasi-particle tunneling spectroscopy.
We have not considered here features expected for inplane magnetic fields. As discussed in Ref. [23], there may exist a high-field superconducting phase, which has a spatially modulated order parameter similar to the helical phase [1,5,6]. Also in this case, the comparative strength of spin-orbit and interlayer coupling is essential, while the Maki parameter does not play an essential role.
2(3 +Bκ 0 ξ 2 c ) (2 +Bκ 0 ξ 2 c ) 2 + ln 1 − 2 2 +Bκ 0 ξ 2 c · (η 4 e,∞ + η 4 o,∞ ) + Qχ nB η 2 e,∞ 2C k +B where the derivative of C k is given as (B2) We can calculate the jump in magnetization at the internal transition whereB =B c in the case J e = J o =: J. Then, as shown in App. A the order parameters swap values at the first-order transition. Thus, only the terms originating from the magnetic orbital part, ∆ ∞ , ξ ∆ , and P , change under the permutation η e,∞ ↔ η o,∞ . Therefore, it is useful to split magnetization jump in to two parts where M m is the magnetization jump arising from the orbital term and ∆M 0 the remaining parts. The latter can be easily calculated by exploiting the fact that the variational parameters η e,∞ and η o,∞ swap values at the A-B transition For the orbital contribution ∆M m we note that the term (K 1 (∆ ∞ ξ ∆ )I 1 (∆ ∞ P ) − K 1 (∆ ∞ P )I 1 (∆ ∞ ξ ∆ )) −2 (1 + 4πχ n )Bκ 2 0 ξ 2 c P 2 (B6) has a strong dependence on the arguments of the modified Bessel functions as K 1 (x) diverges like 1/x near zero. The other originating from the magnetic orbital part, written on the second line of Eq. (B1), is almost unaffected by the permutation of the order parameters and thus neglected. Upon inspecting the definitions of P 2 , ∆ 2 ∞ , and ξ 2 ∆ we observe that, as Qχ n , ξ 2 c << 1, the change in ξ ∆ ≈ ξ c is very small when η e,∞ and η o,∞ swap values. Therefore, we ignore the effect of the permutation η e,∞ ↔ η o,∞ on P and ξ ∆ . Furthermore, we only account for the change in ∆ ∞ in the most singular term K 1 (x). Replacing the modified Bessel function of the second kind with its asymptotic behaviour near zero, K 1 (x) ≈ 1/x, the contribution to the magnetization jump can be written as where we omitted the pre-factor 1/2(1 + 4πχ n ) for compactness and with∆ ∞ the same as ∆ ∞ apart from ψ e,∞ and η o,∞ being permuted. Using I 1 (x) ≈ x/2 the term in the square bracket of Eq. (B7) simplifies to Reinserting the definitions of P , ∆ ∞ and ξ ∆ and settingB =B c we obtain the total magnetization jump at the internal phase transition