Controlling Majorana modes by $p$-wave pairing in two-dimensional $p+id$ topological superconductors

We show that corner Majorana zero modes in a two-dimensional $p+id$ topological superconductor can be controlled by the manipulation of the parent $p-$wave superconducting order. Assuming that the $p$-wave superconducting order is in either a chiral or helical phase, we find that when a $d_{x^2-y^2}$ wave superconducting order is induced, the system exhibits quite different behavior depending on the nature of the parent $p$-wave phase. In particular, we find that while in the helical phase, a localized Majorana mode appears at each of the four corners, in the chiral phase, it is localized only along two of the four edges. We furthermore demonstrate that the Majoranas can be directly controlled by the form of the edges, as we explicitly show in case of the circular edges. We argue that the application of strain may provide additional means of fine-tuning the Majorana zero modes in the system, in particular, it can partially gap them out. Our findings may be relevant for probing the topology in two-dimensional mixed-pairing superconductors.

New platforms for the realization of the MZMs have recently appeared with the advent of higher-order topological states [24][25][26][27][28][29][30][31][32][33][34] which generalize the standard bulkboundary correspondence. Within this class of states, higher-order topological SCs can localize MZMs at interfaces of codimension m > 1. In particular, two-(three-)dimensional second-(third-)order topological SC may host corner MZMs [35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68] with codimension m = d in d spatial dimensions. In a two-dimensional (2D) topological state with an insulating or superconducting bulk gap, the corner zero modes can be obtained by gapping out the first-order edge states with a mass term that features a domain wall in momentum space, realizing a special case of the hierarchy of higher-order topological states [31,69]. In this respect, several concrete ways for the realization of corner MZMs in two dimensions have been proposed so far, as, for instance, by inducing a superconducting gap for the edge states of a topological insulator, either intrinsically [55], or via the proximity effect [35,38,40,56,[58][59][60]. It has also been shown that they may emerge when pairing of an appropriate symmetry is combined with a spin-dependent field, such as spinorbit coupling [41,61,62]. Furthermore, several works have recently proposed means by which the order of a topological superconductor may be manipulated, along with the position of the resulting corner MZMs. Indeed, a first-order topological SC may be promoted to the second order by the application of a magnetic field, with the location of the corner modes determined by the orientation of the field [39,50,63]. It has also been theoretically shown that second-order topological superconductivity can emerge in Josephson junctions, in which case the phase difference between the SCs provides additional means of manipulation [64,65].
We here consider a different scenario in which a variety of MZMs can be generated solely by manipulating the parent p-wave superconducting order in a mixed parity p + id 2D SC. We assume that the p-wave superconducting order, in the absence of any other pairing, hosts a first-order topological state, and exists in either a chiral or helical phase, referring to whether it breaks or preserves time reversal symmetry, respectively. We show that when a d x 2 −y 2 wave superconducting order is induced e.g. via the proximity effect, the system exhibits quite different behavior depending on the parent p-wave phase. In particular, we find that in the helical phase, a localized Majorana mode appears at each of the four corners, as shown in Fig. 1. In the chiral phase, on the other hand, no corner modes appear. Instead, a gap emerges in two out of the four edge modes (Fig. 2). Therefore, the behavior of the Majorana modes can be tuned solely by manipulating the pairing symmetry of the parent topological SC. As we show, the edge geometry can also be relevant in this regard (see Fig. 3), where we display the Majorana states for a circular edge geometry. Finally, we demonstrate that the application of strain may drive a topological phase transition when the parent phase is chiral. As a result, the strain gaps out two out of the four edges, as displayed in Fig. 4.
The rest of the paper is organized as follows. In Section II we introduce the model for the p + id topological SC we consider. Next, we present analytical arguments for the behavior of the resulting edge states in Section III, before moving on to discuss numerical results in Section IV. Finally, we analyze the effect of strain in Section V and present our conclusions together with an outlook in Section VI.

II. MODEL
We employ the standard Bogoliubov-de Gennes formalism to study the system, with the Hamiltonian given as where the corresponding Nambu spinor is , with c k↑ (c † k↓ ) as the annihilation (creation) operator for the quasiparticle with spin up (down) and momentum k . Here, with the blocks describing the normal (nonsuperconducting) state and the pairing, respectively, given by where µ is the chemical potential, m is the quasiparticle mass, the Pauli matrices σ and the unit 2 × 2 matrix, σ 0 , act in the spin space, k F is the Fermi momentum.
Here, ∆ p and ∆ d are the amplitudes of the p-and dwave superconducting orders. For the latter we choose as it features domain walls along the diagonals in the momentum space, located at k x = ±k y , where it changes the sign, and thereby partially gapping out the edge states [71]. We also include a relative phase of π/2 between the two order parameters, implying that the mixed pairing state breaks timereversal symmetry. Notice, however, that the d−wave component preserves the product of the C 4 rotational and time-reversal (T ) symmetries, C 4T = C 4 T implying that the resulting p + id superconducting state may feature the same composite symmetry. This, indeed, occurs when the p−wave component is helical (see below), in which case also a Z 2 topological invariant protects the resulting second-order topological SC [35]. The vector g(k) parametrizes the triplet p−wave superconducting paring, and takes the form g(k) = cos θ k ×ẑ + sin θ (k x + ik y )ẑ, (5) whereẑ is the unit vector pointing in the z direction, assumed to be normal to the 2D plane. The helical phase is found by setting θ = 0, and represents a time-reversal invariant SC, featuring a pair of counterpropagating gapless Majorana edge modes. On the other hand, the chiral phase, which breaks the time-reversal symmetry, is found for θ = π/2. In the following we are interested only in these two special cases, and also refer to the mixed p + id state as either chiral or helical depending upon the phase of the parent p-wave component. Since we only consider the topological regime, we do not include the s−wave pairing in the Hamiltonian in Eq. (2). The bulk spectrum of the Hamiltonian in Eqs.
(2)-(4), is given as with a ± ≡ a ± (k) = k 2 The spectrum for the helical (θ = 0) and the chiral (θ = π/2) phases thus reduces to It is clear that for θ = 0, the bulk band structure always features a gap, as long as ∆ p = 0. The same holds for θ = π 2 , except for the critical value of ∆ d = ∆ p , where the gap closes. We furthermore note that both of the gapped regions ∆ d < ∆ p and ∆ d > ∆ p are topologically nontrivial and feature gapless edge states, as will be evident from the following.

III. ANALYTICAL RESULTS
To understand the effects of the d-wave superconducting order on the edge states present in this system, we introduce an interface which is oriented with an angle α with respect to the horizontal (x−) principal crystalline axis. The Hamiltonian in this rotated coordinate system, whereR is the rotation operator in the Nambu basis defined in Eq. (1), and R = exp[iασ 3 /2] represents the operator of the rotation by the angle α about the z axis in the spin-1/2 representation. The momenta k , and k are related by a rotation about the z−axis with the same angle α, k = R(α)k, with the rotational matrix given by The corresponding edge Hamiltonian is then obtained by projecting a part of the rotated bulk HamiltonianĤ (k ) (Eq. 7) not used to obtain the zero modes onto the subspace spanned by the zero-energy edge modes.
In the following, we assume that both ∆ p and ∆ d are much smaller than the Fermi energy (weak-pairing limit), which we set equal to the chemical potential, µ = 2 k 2 F /2m. In addition, we assume that the wave vector parallel to the interface is much smaller than the Fermi momentum, k || k F . To find the form of the edge Hamiltonian, as previously announced, we first perform the rotation in Eq. (7) and then solve for the zero-energy modes in the absence of the d−wave pairing. To this end, in the rotated Hamiltonian (2), we isolate a part that is proportional to ∆ p and depends on k ⊥ , the wave vector orthogonal to the interface. The total Hamiltonian therefore acquires the form withĤ Here, τ a σ b ≡ τ a ⊗ σ b is a Kronecker product between Pauli matrices in Nambu and spin space, respectively. Since the edge breaks translation invariance, we make the substitution k ⊥ → −i∂ r ⊥ , where r ⊥ is the direction orthogonal to the edge, and solve the resulting differential equation with the boundary conditions ψ(r ⊥ = 0) = ψ(r ⊥ = ∞) = 0. The unperturbed Hamiltonian, given by Eq. (11), admits two degenerate zero-energy eigenstates of the form with a = 1, 2, the spinors and κ = m∆ p /( 2 k F ). The effective Hamiltonian for the edge states is thus found, to first order in the perturbation expansion as the projection which gives where we have kept only the leading term in Eq. (13). The term ∼ k in Eq. (13) vanishes identically after taking the projection (15), and we have neglected the term ∼ k 2 as being less relevant by power counting than the leading one in the low-energy edge Hamiltonian. Note that this term can be, in principle, included in this Hamiltonian following the outlined procedure. The first term in Eq. (16) Clearly, a mass term proportional to cos 2α appears, having domain walls along the two diagonals, located at α = ±π/4. Hence, the edge modes are gapped out, and only the corner modes remain, in agreement with Refs. [35,36]. These corner MZMs are protected by the composite C 4T symmetry and a Z 2 topological invariant [see also the discussion after Eq. (4)].
In the chiral phase, for θ = π/2, we do not have a band crossing which may turn into an anti-crossing and open up a gap at the edge. In this case Eq. (16) takes the form which suggests that the presence of the d wave order parameter only has a trivial effect on the edge states in the regime ∆ d < ∆ p . However, the closing of the bulk gap at ∆ d = ∆ p signals a possible phase transition and therefore something interesting may occur also in this system. It turns out that a selective gapping of the edge states is possible for ∆ d > ∆ p , as can already be seen from Eq. (18). Namely, by taking the values of the angle α = 0 (α = π/2), corresponding to horizontal (vertical) edge, we obtain that the horizontal and vertical edge should be gapped and gapless, respectively. This is exactly what we find by numerical means, as shown in the next section.

IV. NUMERICAL ANALYSIS
In addition to the analytical analysis presented in the previous section, we also numerically study the p + id SC described by Eqs. (2)-(4). The corresponding square lattice Hamiltonian reads with Nambu vector at the lattice site j, ψ j = In the above, t = 2 /2ma 2 , and we use the shorthand notation δn ≡ δ j+n,l . We furthermore remark that the momentum space representation of the above lattice Hamiltonian with periodic boundary conditions is found from Eqs. (2)-(4) by replacing and substituting in the d-wave component the superconducting order parameter. This has the effect of shifting the location at which the bulk gap closes, which is a characteristic feature of the chiral phase, to Furthermore, with µ/t = k 2 F a 2 , it is manifest that the above discretized model becomes equivalent to its continuum counterpart in the limit k F a 1. We now consider the edge states in a square-lattice system. We compute the local density of states at site position j as where E n are the eigenvalues of the Hamiltonian, and v n,j the values of the corresponding eigenvectors at j. The broadening parameter λ is set to 5 × 10 −3 . In the helical phase, as shown in Fig. 1, the edge states quickly vanish with increasing ∆ d , being completely gapped out at ∆ d = 0.2∆ p , and thus leaving only the corner modes.
We turn to the chiral phase, where the gap closing at ∆ d = ∆ c d separates two regions of interest. The region with ∆ d < ∆ c d is topologically equivalent to the case where ∆ d = 0, and we thus expect that the results from Eq. (16) apply here, implying that the d-wave order parameter does not gap out the chiral edge states. This is indeed found to be the case in our numerical analysis, as illustrated in Fig. 2(a), in which the local density of states at ∆ d = t/2 √ 2 = 0.5∆ c d is shown. In contrast, Fig. 2(b), the behavior is different. In that case the horizontal edges are gapped out, but the vertical edge states remain gapless, in agreement with Eq. (18). We furthermore note that a phase shift of π/2 in the relative phase between ∆ d and ∆ p would amount to a π/2 rotation of the result in Fig. 2(b). Therefore, the relative phase between the two pairing order parameters translates into the pattern of the gap at the edge of the system.
We now investigate the edge states in the case of a disk geometry. This is relevant because the behavior of the edge states for any polygonal geometry may be immediately deduced by comparing the corner opening angles with corresponding points on the circle. The modeling of the disk is performed by creating a square grid, and discarding all nodes which fall outside a selected radius, here chosen to be 60 nodes. The results are shown in Fig. 3 for increasing values of ∆ d above the critical value, ∆ c d . Below ∆ c d (not shown), the edge states are uniformly distributed around the entire edge, as expected. Immediately after crossing the critical value, a gap is opened up in the edge states at angles around γ = {π/2, 3π/2}, as shown in Fig. 3(a) for ∆ d = 1.5∆ c d , and γ is the polar angle of the circle. This is consistent with the partial gapping of the edge states observed in Fig. 2(b), which probes the same angles, along with the angles {0, π}, which are gapless. A further increase in ∆ d increases the modulation of the density of states along the edge, and produces additional gapped regions, as can be seen in Fig. 3(b)-(d), which correspond to ∆ d /∆ c d = 3, 4.5, and 6, respectively. Furthermore, the gapped circle sector surrounding γ = {π/2, 3π/2} is seen to narrow as ∆ d becomes larger, but never closes completely, consistent with the domain wall structure of the d x 2 −y 2 -wave pairing.

V. THE EFFECT OF STRAIN
We investigate strain as a potential means to manipulate the edge states. We model its effects by introducing a small strain field to the system, Here, ε xx and ε yy represent axial strain, defined as positive for tensile strain, and ε xy = ε yx represents shear strain. In the presence of such a strain field, the spatial coordinates transform as r i = (δ ij + ε ij ) r j , which implies that, to linear order in the strain tensor, the momentum transforms as By replacing k → k in Eqs.
(2)-(4), then inserting Eq. (26) and retaining only terms up to first order in ε, we find that additional strain-dependent terms are introduced in the Hamiltonian, which may have an effect on the edge states. For an edge with an arbitrary orientation α, the corresponding Hamiltonian, after incorporating the effects of strain, read with Γ p (ε, α) = 1 −ε + 1 2 δε cos 2α + ε xy sin 2α (28) whereε = (ε xx + ε yy )/2 and δε = ε xx − ε yy . We see that strain produces an effective renormalization of the p and d wave order parameters in the edge Hamiltonian, by Γ p and Γ d respectively. Furthermore, for anisotropic strain, Γ d acquires a contribution independent of α. This implies that the corner modes in the helical phase, given by Eq. (17) can be gapped out by strain. This is also consistent with the breaking of the C 4T symmetry by strain, which protects the topological state. However, the gapping only occurs once the mass domain wall at the corners is removed, which requires rather large strains. For instance, with a uniaxial tensile strain of ε xx = ε, the corner modes are gapped out for ε > 50 %. This is much larger than the capacity of any known material, with graphene as the closest contender, reported to sustain strains of up to 25 % [72]. In any case, these levels of strain are certainly far beyond the linear strain regime considered herein, implying the stability of the corner modes against strains in the experimentally realizable range.
In the chiral phase, we consider the effect of strain by numerical means. To this end, we replace the lattice parameter a with directionally dependent equivalents, a x and a y , satisfying a i = (1 + ε ii )a. We ignore shear strain. In the regime ∆ d > ∆ c d , we can conclude from Eq. (27) that ∆ d and ∆ p are renormalized slightly differently by strain, which agrees with the numerical analyses for ∆ d < ∆ c d . Hence, strain may be used to change the ratio between the amplitudes of the d and p wave superconducting order, and thus cause a transition between the two phases exhibiting a different behavior of the edge states, which agrees with our numerical analysis. Setting ∆ d = ∆ c d , which in the unstrained case is gapless, as shown in Fig. 4(a), we perturb the system by applying a uniaxial strain of 10 % along the x and y directions, respectively shown in Fig. 4(b) and (c). In the former case, it can be seen that the system is pushed into the state with uniform gapless edge modes, whereas in the latter case, the selectively gapped state is entered. Therefore, in the chiral phase, strain can be used to tune the form and the localization of the edge states.

VI. CONCLUSIONS AND OUTLOOK
In conclusion, we have studied the Majorana zero modes in a 2D SC with mixed p-and d-wave pairings, considering both the helical and chiral p-wave phases. For the parent helical p-wave order, we found that the effect of the d-wave order parameter is to gap out the edge states, leaving only zero-energy Majorana corner modes. In the case of the parent chiral p-wave pairing, we showed that the edge states can be partially gapped out above a critical value of the d-wave order parameter. Therefore, the parent p-wave phase can control the form of the Majorana modes in the p + id topological SC. We found that the localization of the Majorana modes can also be tuned by the geometry of the edges, as implied by our results in the circular edge geometry. Moreover, the localization of the MZMs in a polygonal geometry may be inferred by identifying the opening angle of a given corner with a corresponding point on the circle. We have also investigated the effect of strain and found that the higher-order topological SC produced in the parent helical phase is robust against strain up to the experimentally reachable values. On the other hand, in the chiral phase, we showed that strain can be used to induce a transition between topologically distinct phases with gapless and partially gapped edge states. Hence, the application of strain, through the strain-induced topological phase transition, may indeed provide a direct means of manipulating the Majorana modes.
The experimental realization of such superconducting system is a challenging issue at present, but we will nevertheless indicate a few promising avenues. First, 2D p + id second-order topological SCs may be hosted in doped second-order topological insulators, as has been recently advocated in Ref. [48]. Another available route is to make use of the proximity effect between a topological p-wave SC and a d-wave SC. Typical examples of the latter are cuprates such as YBa 2 Cu 3 O 7−δ [73], while the former is more difficult to find. Indeed, the leading candidate for topological p-wave superconductivity, Sr 2 RuO 4 , has recently been shown most likely not to be of p−wave nature [74,75]. Intensive research is, however, still ongoing and other potential candidates are the uranium-based heavy-fermion compounds, of which UTe 2 seems to be most promising one [76].
As we are concerned here only with 2D effects, we do not have to rely on bulk superconducting properties. This broadens the choice of materials somewhat. For instance, topological superconductivity in surface states can be achieved via the Fu-Kane construction [77], whereby a conventional s-wave SC is placed in proximity with a topological insulator. Such a system can host edge states upon breaking of time reversal symmetry, e.g. by introducing a magnetic vortex. This behavior has recently been discovered to occur at surfaces of some materials, such as the iron-based SC FeTe 1−x Se x [78][79][80][81][82][83], and is anticipated in PbSaTe 2 [84], which thus places them as promising candidates for our purposes. Hence, our Hamiltonian may effectively describe the Majorana bound states on the surface of such materials, when they are proximitized with a d-wave SC, either as a bilayer, or as a Josephson coupling, as was suggested in Ref. [35]. We also point out that a proximity effect between niobium and graphene, through a single layer of chiral molecules, has been demonstrated to display signatures which may be indicative of chiral p-wave superconductivity [85], and might provide additional means of realizing in the future the model we studied here.
Finally, we remark that a cuprate d-wave SC typically has a significantly larger critical temperature than the pwave SCs discussed above. For instance, YBa 2 Cu 3 O 7−δ may reach a T c ∼ 100 K [86], whereas FeTe 0.55 Se 0.45 has T c ∼ 15 K [83], and T c is doping-dependent in both cases. Hence, the relative size of the two superconducting order parameters, ∆ p and ∆ d , might be tunable by changing the temperature and the doping. In any case, we hope that our findings will motivate further experimental efforts to demonstrate the tunability of the MZMs in mixed-pairing topological SCs by both the parent pwave state and an external nonthermal tuning parameter, such as strain. Our results should also motivate further searches for material platforms where the proposed scenario can be relevant.
In closing, we point out that our mechanism may also be relevant for three-dimensional SCs which can be realized in doped octupolar (third-order) Dirac insulators, hosting corner MZMs [68]. This problem is, however, left for future investigation.