Spin response and topology of a staggered Rashba superconductor

Inversion symmetry is a key symmetry in unconventional superconductors and even its local breaking can have profound implications. For inversion-symmetric systems, there is a competition on a microscopic level between the spin-orbit coupling associated with the local lack of inversion and hybridizing terms that `restore' inversion. Investigating a layered system with alternating mirror-symmetry breaking, we study this competition considering the spin response of different superconducting order parameters for the case of strong spin-orbit coupling. We find that signatures of the local non-centrosymmetry, such as an increased spin susceptibility in spin-singlet superconductors for $T\rightarrow 0$, persist even into the quasi-three-dimensional regime. This leads to a direction dependent spin response which allows to distinguish different superconducting order parameters. Furthermore, we identify several regimes with possible topological superconducting phases within a symmetry-indicator analysis. Our results may have direct relevance for the recently reported Ce-based superconductor CeRh$_2$As$_2$ and beyond.


I. INTRODUCTION
In superconductors with inversion symmetry, evenand odd-parity gap functions are distinct by symmetry and, thus, correspond by the Pauli principle also to spin-singlet and spin-triplet pairing states, respectively. As a consequence, magnetic response can be used to distinguish the two cases. When the system lacks inversion, however, spin-singlet and spin-triplet can mix [1]. In addition, the spin-orbit coupling associated with the broken symmetry-an example is the spinorbit coupling of Rashba type for broken in-plane mirror symmetry-strongly restricts the possible spin-triplet order-parameter components, in other words it fixes the direction of the d vector [2]. Moreover, even in the case of dominant spin-singlet or spin-triplet order parameters, the magnetic response, such as the spin susceptibility or critical fields, is not a feasible distinguishing probe anymore, as for spin-singlet superconductors, a finite spin susceptibility for T → 0 and unusually high critical fields can be expected. Beyond that, it was shown that topological properties of the phases, such as vortex bound states or surface flat bands, can serve as fingerprints of the respective phases [3,4].
Signatures of non-centrosymmetry can survive even in inversion-symmetric systems [5,6]. In particular, a crystal comprising weakly-coupled sublattices whose subunits locally lack inversion, such as in the hexagonal SrP-tAs [7] or even in some high-temperature cuprates [8], can exhibit an unconventional magnetic response or intriguing spin textures. In addition to specific crystal structures, this local non-centrosymmetricity can arise in artificial superlattices such as the regular stacks of superconducting CeCoIn 5 alternating with layers of YbCoIn 5 [9][10][11]. Note that the effect of local noncentrosymmetricity depends on the relative strength of inversion-breaking-induced spin-orbit coupling and intersublattice hybridization. The focus of most studies has, thus, been on quasi-two-dimensional systems with weak c-axis dispersion.
The recently discovered heavy-Fermion superconductor CeRh 2 As 2 with its tetragonal crystal structure belongs also to the class of locally non-centrosymmetric superconductors. In particular, it consists of layers with alternating inversion-symmetry breaking. The upper critical field directed along the c axis (perpendicular to the staggered layers) extrapolates to ∼ 14 T at zero temperatures, which lies far beyond the paramagnetic limiting field H p ∼ 0.5 T for a critical temperature T c ≈ 0.26K [12]. Furthermore, the upper critical field shows a pronounced kink for a field H ≈ 4 T. This anomaly strongly suggests a change in the order-parameter symmetry upon increasing magnetic field as also found in recent Ginzburg-Landau studies [13,14]. Note that the critical field for in-plane directions, on the other hand, extrapolates to only approximately 2T.
Unlike most staggered systems studied so far, CeRh 2 As 2 is expected to have a rather strong c-axis dispersion, in other words it is a three-dimensional (3D) system. Being a Ce-based superconductor, we expect also a sizable spin-orbit coupling. Motivated by these observations, we revisit the physics of locally noncentrosymmetric superconductors in situations where both, the inter-sublattice hopping and the spin-orbit coupling strength are comparable to each other and the overall band width.
In particular, we investigate in detail a microscopic model of a layered system, where mirror symmetry is broken in a staggered fashion. By design, our model allows us to investigate how the system evolves from the 2D limit to the truly 3D case. For this purpose, we first study the spin susceptibility in the normal state, which leads us to identify four distinct regions, going from quasi-twodimensional (q2D) all the way to truly 3D. Then we discuss possible order parameters and analyze their spin response, showing how both in-plane and out-of-plane fields are necessary to distinguish them. Eventually, we address the topological phases in the fully gapped case. Note c FIG. 1. Schematic of the system with the triangles indicating the mirror symmetry breaking. Note that, for simplicity, we choose the layers to be spaced equidistantly in the z direction with unit-cell size c. that since the system retains inversion, we can use the recently developed concept of symmetry indicators [15][16][17][18] to show that the system can realize both first-and second-order topological superconducting phases.

A. Microscopic model
We consider in the following a system of stacked layers, where each layer lacks a mirror symmetry in such a way that the resulting Rashba spin-orbit coupling alternates in sign, see Fig. 1. Such a system has centers of inversion, which lie in between two neighboring layers. However, electrons moving within an individual layer are subject to a spin-orbit coupling of Rashba type. For concreteness, we choose these layers to consist of a square lattice with point group symmetry C 4v , while the full structure involves the tetragonal point group D 4h . On a single layer, the electrons are governed by the Hamilontian where c † ks (c ks ) creates (annihilates) an electron with momentum k and spin s. Setting the in-plane lattice constant a = 1, describes the dispersion due to nearest-neighbor hopping, enters the expression for the Rashba spin-orbit coupling, and σ denote the Pauli matrices. The three-dimensional system with a staggered stacking of such layers, as shown in Fig. 1, is then described by the Hamiltonian with the 4 × 4 matrix where the inter-layer hopping is given by and (ξ − k ) 3 = 0. Note that we have set the z-axis lattice constant c = 1. Furthermore, we have introduced the Pauli matrices τ i , i = 1, 2, 3, acting on the sublattice space of layers for the operators ψ † ks = (c † eks , c † oks ) denoted by even (e) and odd (o) sublattice index.
The eigenenergies of this Hamiltonian are doubly degenerate and are given by which we denote in the following as ξ α with α = ± neglecting the momentum index for shorter notation. Further, we introduce the chemical potential µ. For concreteness, we use in the following λ = 0.5t with t the energy unit, and we choose µ so as to fix the density of electrons to n tot = 0.15.

B. Magnetic Response
It is instructive to first consider the normal-state magnetic response of this system. For this purpose, we introduce the normal state Green's function, which is defined through with ω n = (2n + 1)πT the fermionic Matsubara frequencies with n ∈ Z and T the temperature. This 4 × 4 Green's-function matrix can explicitly be calculated by inverting Eq. (10), where we introduced andf In order to investigate the system's response to a magnetic field, we calculate the (normal-state) uniform, static spin susceptibility (q = 0, ω = 0), which reads where the trace runs over spin and layer indices. Performing the trace first, we find that due to spin-orbit coupling the susceptibility has a generic form with two contributions: the first is with n F (ξ) the Fermi distribution. This term corresponds to the Pauli-like susceptibility. In other words, it is proportional to the spectral density at the Fermi level, S k (µ), and after the k-integration, the density of states of the two bands at the Fermi level, N (µ). The second is and originates from inter-band processes due to the spinorbit coupling and thus, we refer to it as a van Vleck term. The susceptibility is a combination of these Pauli and van Vleck contributions, in particular for fields in the z direction, i = j = z. For fields along the x direction, we find and analogously for fields in y direction. This result generalizes the pure Rashba case [1,19], which is recovered when setting t z = δt z = 0. In the other limit, t z , δt z λ, the weight of the van Vleck term is strongly reduced and the Pauli susceptibility dominates. Consequently, we can use the ratio of van Vleck susceptibility to the total susceptibility for fields along the z axis as a measure of the inversion-symmetry breaking in the system. Figure 2 shows the ratio of van Vleck susceptibility to total susceptibility of the normal state as a function of the interlayer hopping t z . The susceptibility shows a distinct behavior depending on the Fermi surface topology as indicated by the insets: In the quasi-two-dimensional (q2D) regime, where both bands are partially filled and with open Fermi surfaces, the van Vleck contribution drops rapidly with increasing t z . Then, in the first quasi-3D regime (q3DI) with one closed Fermi surface, χ 0 vV /χ 0 decreases essentially linearly, before dropping again more rapidly once the closed Fermi surface disappears in the second quasi-3D regime (q3DII). Finally, the remaining Fermi surface closes, at which point the full system is 3D and behaves completely centrosymmetric. As we will discuss in Sec. IV, the different Fermi-surface topologies are also connected to different possible topological phases for a triplet superconductor.
Before we move to calculating the spin susceptibility in the superconducting state, we can use the above finding to discuss what we expect for a spin-singlet gap function. As the Pauli term depends on the density of states at the Fermi level, which vanishes for a full (spin-singlet) gap, the corresponding susceptibility is expected to vanish as well. The van Vleck term, however, should stay constant, at least for | f k | 2 + | ξ − k | 2 ∆, with ∆ the superconducting gap. Comparing Eqs. (18) and (19), we further expect the susceptibility for fields in-plane to be half the size of the susceptibility for fields along z for T → 0, independent of the microscopic details, as long as the system has C 4 symmetry.

A. Superconducting order parameters
Before calculating the spin response of possible superconducting phases, we discuss the possible pairing states of interest here. While mixing of spin-singlet and spintriplet order parameters is allowed within the plane and microscopically supported by the Rashba spin-orbit coupling, the order parameter can still be classified as even or odd due to the global inversion symmetry [5]. Similar to the non-centrosymmetric situation [1,2,19], most spintriplet order parameters are suppressed by the spin-orbit coupling in the staggered case, too [5]. The energetically most stable gap structures correspond to intra-band pairing, in other words to gap functions that are diagonal in the bands with energies given by Eq. (9) [20,21]. As these order parameters can also be treated analytically, we will focus exclusively on them. As a consequence, the order parameters we study have no in-plane mixing of spin-singlet and spin-triplet channels. Note that we assume in the following the superconducting gap to be the smallest energy scale in the problem.
For this intra-band pairing, we focus on the most symmetric gap functions leading to three different gap structures: a spin-singlet that transforms like A 1g , a spintriplet order parameter with its d-vector in plane that transforms like A 2u , and a spin-triplet order parameter with the corresponding d-vector along the z axis, that transforms like A 1u in D 4h [5]. We briefly discuss in the following all three order parameters.
For the spin-singlet channel, the intra-band order parameter has the form This order parameter describes fully-gapped s-wave spin- For the spin-triplet gap, Cooper pairing can be both within the plane and between neighboring planes: The first, where d k = (d x k , d y k , 0) T f k , is similar to the case of non-centrosymmetric superconductors [1,19] and transforms with A 2u . This gap vanishes for k x = k y = 0, such that for closed Fermi surfaces as found in the q3DI and 3D regime, the gap has point nodes. Notice that in the non-centrosymmetric case with point group C 4v , the gap functions of A 1g and A 2u symmetry mix, resulting instead in line nodes for closed Fermi surfaces and dominant triplet contribution.
In addition the A 2u order parameter with out-of-plane spin-triplet pairing is allowed due to the staggered nature of the mirror-symmetry breaking with a nearest-neighbor gap function While we can in general write a gap function of this form that is entirely intra-band, we will for simplicity consider in the following δt z = 0, such that ∆ ⊥ (k) = (d z k σ z )iσ y τ 1 with d z k ∝ sin k z /2 as the lowest-order basis function [22]. This corresponds to an inter-layer pairing between nearest layers. This gap function has a line node for k z = 0, which can be removed by mixing in the other spin-triplet component of A 1u symmetry, namely d k = (sin k x σ x + sin k y σ y )iσ y τ 0 . This order parameter thus in general allows for a full gap.
While knowledge of the pairing interaction is needed to determine the leading instability in the system, we would expect the dominant channel to be in-plane in the quasi-two-dimensional limit. Consequently, the most probable order parameters are the spin-singlet and the spin-triplet of A 2u symmetry. Only when approaching the three-dimensional limit, in other words t z ≈ t, an inter-layer pairing becomes feasible. As mentioned, this last pairing channel can lead to a full gap even in the three-dimensional limit, which in turn allows for topological superconductivity as we discuss in Sec. IV.

B. General spin response
The spin response of a superconductor can be calculated similarly to the one in the normal state, Eq. (15), using [23] where G(k, ω n ) and F (k, ω n ) are the normal and anomalous Green's functions, respectively. Given the intraband order parameters of Eqs. (20), (21), or (22), these can be calculated explicitly using the Gor'kov equations, see App. A. The resulting structure of both Green's functions is very similar to the normal state Green's function of the previous section, Eq. (11). In particular, the normal Green's function reads with and E α = ξ 2 α +|∆| 2 . The normal Green's function does not depend on the specifics of the intra-band gap function, other than its (momentum-dependent) magnitude, in other words |∆| 2 = |ψ| 2 for the spin-singlet, |∆| 2 ≡ |∆ k | 2 for the intra-layer spin-triplet, and |∆| 2 ≡ |∆ ⊥ k | 2 for the inter-layer spin-triplet case.
For the anomalous Green's function, on the other hand, the gap structure enters explicitly, To proceed, we separate the susceptibility in Eq. (23) into a 'normal' and an 'anomalous' part, χ ij = χ n ij + χ a ij . The trace of the normal part is the same for all three order parameters and, as in the discussion in the previous section, yields and similarly (29) As the gap structure enters the anomalous Green's function, we discuss in the following first the susceptibility for the case of the spin-singlet and then the case of the spin-triplet order parameters.

Spin-singlet order parameter
Using the s-wave singlet gap function in the anomalous Green's functions, Eq. (26), we find for the anomalous part of the susceptibility and similarly (31) We can again separate the susceptibility into two contributions, namely and The total susceptibility after Matsubara summation is given by for fields along z and, for fields along the x or y direction, the susceptibility has the form The Matsubara sums are evaluated in App. B and yield for the Pauli term and, for ∆ | ξ − k | 2 + | f k | 2 , we can approximate for the van Vleck susceptibility.

Intra-layer spin-triplet pairing
For spin-triplet pairing, we start with the intra-layerpairing case with order parameter ∆ (k) = ( d k · σ)iσ y τ 0 , where d k f k . For out-of-plane fields, i = j = z, the anomalous part of the susceptibility yields For in-plane fields, i = j = x, we find (39) Combining Eqs. (28) and (29) with (38) and (39), we thus find for the total susceptibility two additional contributions compared to the spin-singlet case. First, which, after performing the Matsubara sum (App. B) yields The second, yields for | ξ − k | 2 + | f k | 2 ∆ the same van Vleck contribution as Eq. (37), χ − vV (k) ≈ χ + vV (k) ≡ χ vV (k). For fields out of plane, we thus find the susceptibility and for the susceptibility for in-plane fields

Inter-layer pairing
For the inter-layer pairing with ∆ ⊥ (k) = (d z k σ z )iσ y τ 1 , the traces of the anomalous part of the susceptibility in Eq. (23) yield for out-of-plane fields and for in-plane fields, i = j = x, we find For fields out of plane, we thus find the susceptibility and for in-plane fields the susceptibility Because the Pauli susceptibility in Eq. (36) vanishes for T → 0 for the spin-singlet state, only the temperatureindependent van Vleck contribution of the normal state remains. Moreover, comparing Eqs. (34) and (35), we find for the spin-singlet pairing χ s z = 2χ s x . In order to discuss the behavior of the susceptibility for the triplet states in the T → 0 limit, we rewrite the term in Eq. (41) as which, for T ∆ t, reduces to a derivative of a step function around µ and, thus, approximately to the normal-state Pauli contribution to the susceptibility. Hence, we find that the spin susceptibility for outof-plane fields does not decrease for intra-layer pairing, while it follows the trend of the spin-singlet case for interlayer pairing. For in-plane fields, however, the intra-layer pairing state reduces the spin susceptibility, while the inter-layer pairing state is not paramagnetically limited. These findings are, thus, in accordance with the general expectation that for a d-vector perpendicular to the magnetic field the susceptibility stays constant, while the contribution parallel to the field reduces the susceptibility. Note that here, the van Vleck term is generally not affected by superconductivity, since it arises form interband contributions, and the band splitting is governed Figure 3 shows the susceptibilities for all three order parameters as a function of the z-axis hopping t z , where, as mentioned, we use δt z = 0 for simplicity. As discussed above, the zero-temperature in-plane susceptibility is always reduced to half of the out-of-plane susceptibility for (fully-gapped) spin-singlet pairing, while this is not the case for spin-triplet order parameters. Note that the inter-layer state has a larger susceptibility in these calculations due to the line nodes at k z = 0. These nodes can in principle be gapped out, such that the out-of-plane susceptibility of the spin-singlet and the inter-layer spintriplet would be equal. Finally, note that the in-plane (out-of-plane) pairing state is limited by orbital depairing for out-of-plane (in-plane) fields.

IV. TOPOLOGICAL CONSIDERATIONS
In this section, we investigate possible non-trivial topology of the superconducting phases in the different regimes identified in Sec. II. In systems invariant under inversion with an odd-parity superconducting order parameter, the topology of the superconducting phase can be identified using so-called symmetry indicators (SI) [15][16][17]. These indicators are defined in terms of the inversion eigenvalues of the occupied normal-state bands at time-reversal-symmetric momenta (TRIMs), where −k ≡ k up to reciprocal lattice translations. In particular, the 3D inversion-symmetry indicator κ 3D is a Z 8 -valued quantity, such that an odd SI indicates a strong topological superconductor, while κ 3D = ±2 and κ 3D = 4 indicate second-and third-order topological phases, respectively.
The inversion-symmetry indicator for a time-reversal symmetric odd-parity superconductor is defined as where n ± k,N is the number of occupied bands at the TRIMs with inversion eigenvalue ±1. To calculate the indicator, it is convenient to rewrite the dispersion in the z direction, Eqs. (7) and (8) as This corresponds to a gauge transformation such that the Fourier transform does not resolve the unit-cell structure anymore. The advantage, however, is that now IH k I −1 = H −k with I = σ 0 τ 1 and we can straightforwardly read off the inversion eigenvalues of the normalstate bands. We are interested in inversion eigenvalues at two TRIMs Γ = (0, 0, 0) and Z = (0, 0, π). At these two momenta, the two bands with energies given by Eq. (9), namely ξ Γ,± = −4t ∓ 2t z and ξ Z,± = −4t ∓ 2δt z , have inversion eigenvalues ±1. Depending on the Fermi-surface topology and hence, the number of bands occupied at the −π 0 π −π 0 π t z = 0.05 4. Projection of the Fermi surface onto kx-kz plane for different values of interlayer hopping parameters tz and δtz. Two colors correspond to two different momentum sectors with inversion eigenvalues respectively +1 and −1. TRIMs, the system can realize a first-or second-order topological superconductor, see Fig. 4. A third-order topological superconductor, however, is not possible at least with this simplified band structure, since for both bands occupied, the inversion eigenvalues are opposite and cancel in Eq. (50). Table I summarizes the symmetry indicators for the two spin-triplet (odd-parity) order parameters considered here. In particular, the A 2u state, which pairs electrons within the layers and is, thus, more probable in the q2D limit, allows for a second-order phase when only one open Fermi surface exists. In the other potentially non-trivial cases, this gap structure has point nodes for k x = k y = 0. The A 1u state, which can be fully gapped for any Fermisurface topology, thus allows in principle for both firstand second-order phases. As the A 1u state corresponds to an anisotropic Balian-Werthammer state known from the B phase of 3 He, a non-trivial topology might not be surprising [24]. Note, however, that this gap structure is rather unlikely in the quasi-2D limit.

V. CONCLUSION
In some globally centrosymmetric systems, the lack of inversion symmetry in subunits can influence the physical properties significantly. Such remnants of noncentrosymmetricity are particular interesting in the context of superconductivity, where inversion symmetry yields a strict distinction between spin-singlet and spintriplet Cooper pairing. In locally non-centrosymmetric systems, however, the two may mix in a very characteristic way.
In our work, we have investigated signatures of local noncentrosymmetricity for the case of a layered system with staggered layer-intrinsic Rashba spin-orbit coupling. Hereby, we have focused on different regimes going from a quasi-2D to fully 3D band structure upon tuning the interlayer hybridization. First looking at the normal state spin susceptibility, we identified four different regimes, characterized by their Fermi surface topology.
Unlike the case of globally non-centrosymmetric superconductors, where only one kind of spin-triplet state is feasible, we identified two such spin-triplet states in the staggered case. We analyzed the resulting three order parameters, namely a generic spin-singlet and those two spin-triplet phases and show that the comparison between their behavior under in-plane and out-of-plane magnetic fields allows for the order parameters' distinction. In particular, their low-temperature spin susceptibility displays different behavior that influence the paramagnetic limiting effects for different field orientations.
For the spin-triplet order parameters, we have finally explored possible topological phases within the symmetry indicator framework, and identified both first-and second-order topological phases. Third-order topology is, on the other hand, not possible in the minimal model considered here.
Finally, we comment on the relevance of our results concerning CeRh 2 As 2 [12]. For c-axis fields, the lowfield state seems to be paramagnetically limited, which could have either spin-singlet or inter-layer spin-triplet character. However, as the in-plane fields show a smaller critical field for CeRh 2 As 2 , it is rather likely that the lowfield phase constitutes spin-singlet pairing. Note that the spin-triplet phase that would in this scenario emerge at high fields, is certainly topologically trivial, as there are no first-order TSC in 3D. A generic possibility would, however, be a Weyl superconductor. As the field-induced phase transition indicates that the spin-triplet state is close by in parameter space, it might be possible to stabilize it also with other means, thus inducing possibly topological superconductivity.

Pauli-like terms
To evaluate the Matsubara sum of the first two contributions, Eqs. (36) and (41), we need to calculate the residues of (z + ξ λ ) 2 ± |∆ k | 2 (−z 2 + ξ 2 λ + |∆ k | 2 ) 2 n F (z), at the singularities not stemming from the Fermi distribution function. We find two second-order poles, namely z = ±E = ± ξ 2 + |∆| 2 (for simplicity, we omit the index α here). The residue is then Using that is an even function of z, we find for the sum of the two residues Res E + Res −E = 1 4T cosh 2 (E/2T ) .
Similarly, replacing |∆| 2 with −|∆| 2 in Eq. (B2) for χ − P , we find (B5) Finally, we write for the two Pauli contributions (B7) The terms with the cosh vanish for T → 0, such that only the first term of the latter equation survives,

van Vleck-like terms
The second set of contributions, Eqs. (37) and (42), has four poles at ±E ± and they are all first order. At T = 0, the poles with positive energy thus lead to a vanishing Fermi function while the Fermi functions for the negative energies yield 1. We thus find The energy scale for the van Vleck susceptibility is given by | f k | 2 + | ξ − k | 2 and since we assume | f k | 2 + | ξ − k | 2 ∆, we proceed by setting the gap to zero. Then, we find χ ± vV (k) ≈ χ vV (k) = 2µ 2 B (|ξ + | + ξ + )(|ξ + | + ξ − ) This expression yields exactly zero if both, ξ + and ξ − , are either positive or negative. However, when only one of them is negative, while the other is positive, i.e. n F (ξ + )− n F (ξ − ) = 0, we find the van Vleck susceptibility.