Majorana-Weyl cones in ferroelectric superconductors

Topological superconductors are predicted to exhibit outstanding phenomena, including non-abelian anyon excitations, heat-carrying edge states, and topological nodes in the Bogoliubov spectra. Nonetheless, and despite major experimental efforts, we are still lacking unambiguous signatures of such exotic phenomena. In this context, the recent discovery of coexisting superconductivity and ferroelectricity in lightly doped and ultra clean SrTiO$_3$ opens new opportunities. Indeed, a promising route to engineer topological superconductivity is the combination of strong spin-orbit coupling and inversion-symmetry breaking. Here we study a three-dimensional parabolic band minimum with Rashba spin-orbit coupling, whose axis is aligned by the direction of a ferroelectric moment. We show that all of the aforementioned phenomena naturally emerge in this model when a magnetic field is applied. Above a critical Zeeman field, Majorana-Weyl cones emerge regardless of the electronic density. These cones manifest themselves as Majorana arcs states appearing on surfaces and tetragonal domain walls. Rotating the magnetic field with respect to the direction of the ferroelectric moment tilts the Majorana-Weyl cones, eventually driving them into the type-II state with Bogoliubov Fermi surfaces. We then consider the consequences of the orbital magnetic field. First, the single vortex is found to be surrounded by a topological halo, and is characterized by two Majorana zero modes: One localized in the vortex core and the other on the boundary of the topological halo. Based on a semiclassical argument we show that upon increasing the field above a critical value the halos overlap and eventually percolate through the system, causing a bulk topological transition that always precedes the normal state. Finally, we propose concrete experiments to test our predictions.


I. INTRODUCTION
Finding robust experimental realizations of topological superconductivity is an important goal, both for fundamental research of topological matter and for possible applications to quantum technology [1][2][3]. However, materials which naturally host such exotic ground states are scarce. Moreover, measuring non equivocal signatures of topological superconductivity is an outstanding experimental challenge [4][5][6], because such signatures are often obscured by imperfections in the sample or probe. Most candidate materials also realize low-dimensional topological superconducting states. Thus, new candidate bulk superconductors might help overcome such challenges.
The coexistence of superconductivity and ferroelectricity in low-density systems [26][27][28][29][30][31] opens new opportunities in this context. A ferroelectric crystal breaks inversion symmetry spontaneously and therefore can be easily manipulated. Moreover, such systems are often close to their ferroelectric transition, where the dielectric constant is huge [32,33]. As a consequence, the influence of disorder is dramatically suppressed [34,35]. These properties make low-density superconductors close to a ferroelectric quantum critical point prime candidates for engineering unconventional superconducting states.
Motivated by the physics in ferroelectric STO, we revisit the problem of a Rashba spin-orbit coupled superconductor subject to a magnetic field [47,48], where we focus on the case of three spatial dimensions. Rashba spin-orbit coupling originates from the combination of inversion breaking by the ferroelectric moment and atomic spin-orbit coupling [49,50]. Therefore, the axis of the Rashba spin-orbit coupling can vary in space and may also be externally manipulated.
In the absence of superconductivity and magnetic fields, the Fermi surfaces are spin split everywhere in momentum except for two pinching points, which lie along the axis of the polar vector (Fig. 1a). Consequently, in the superconducting state, pair breaking is strongest at the vicinity of these points. When the magnetic field exceeds a critical threshold, the gap closes along this polar axis causing four Majorana-Weyl points to emerge, accompanied by surface Majorana arcs. We then show that the Majorana-Weyl cones can be tilted by tuning the angle between the polar moment and field, such that the superconductor becomes type-II-Weyl with Bogoliubov Fermi surfaces [51,52] above a certain critical angle. We also study the Fermi arcs forming on domain walls between different polarization directions. We find that chiral surface states do not appear for all angles of the magnetic field.
Finally, we turn to the more realistic scenario, where the field is non-homogeneous and penetrates the sample through line vortices. We first study the single vortex problem, where we show that the magnetic field always exceeds the critical threshold close enough to the center, forming a topological halo surrounding the vortex. We show that each such vortex has a single zero mode in its core with a counterpart at the boundary of the halo, yielding corresponding signatures in the tunneling density of states. Then, as the magnetic field is increased towards H c2 the density of vortices increases and the halos begin to overlap, forming larger topological regions. As a consequence we predict that the trivial-superconducting and normal states are always separated by a putative topological phase in any polar superconductor. The topological phase is characterized by percolation of the halos, akin to a transition between integer quantum Hall states.
The rest of this paper is structured as follows. In Section II, we describe the model and show in the meanfield picture, neglecting the orbital effects of the magnetic field, that Majorana-Weyl superconductivity develops when the magnetic field exceeds certain threshold. In Section III, we discuss the Fermi arcs on surfaces and interfaces between the ferroelectric domains. In Section IV, we consider a more realistic model taking into account orbital effects of the magnetic field. We show that in addition to a Majorana string located in the core, an isolated vortex is surrounded by a chiral Majorana mode with the wavefunction peaked at a finite distance from the core. Based on semiclassical considerations, we propose that with the increase of the magnetic field towards H c2 , there is always a percolation-type phase transition to a bulk Majorana-Weyl superconductivity, at which chiral modes going around each vortex overlap. In this section we also calculate contribution from the Majorana modes to the tunneling density of states. Finally, we give our conclusions with emphasis on experimental consequences caused by the physics considered in Section V. Throughout the paper we work in units in which = k B = 1.

II. MAJORANA-WEYL SUPERCONDUCTIVITY IN THE PRESENCE OF A ZEEMAN FIELD
We now describe the microscopic model. We start with the coupling between the optical phonon displacement P and the conduction electrons [53][54][55][56][57] where ψ k is an annihilation operator for the electron with momentum k andλ is a coupling. This term has its microscopic origin as a consequence of combined effect of spin-orbit coupling and interorbital hybridization allowed by inversion breaking [49].
In the ferroelectric phase, the displacement field P develops a non-zero expectation value. We emphasize however, that this expectation value does not imply the presence of long-ranged electric fields, which are always screened by the itinerant electrons beyond the Thomas-Fermi length scale. The term "ferroelectricity" is commonly used in the literature to describe the polar state even when it is metallic. However, in this case the "ferroelectric" phase transition refers to a structural transition, where inversion symmetry is broken. As a result, the coupling Eq. (1) leads to the celebrated Rashba spin orbit wheren is a unit vector parallel to the ferroelectric order parameter P and λ = | P |λ.
Additionally, we consider a Zeeman coupling to an external magnetic field B, and neglect its orbital effects pro tem [58]. Without loss of generality we align the zaxis with the local polarization P (hencen =ẑ), and obtain the dispersion Hamiltonian where σ = (σ x , σ y , σ z ) is a vector of Pauli matrices in spin space and we have assumed the dispersion k = k 2 2m − µ is spherically symmetric. In the following, we work in units in which gµ B /2 = 1.
We next add an attractive interaction between electrons, which causes a Cooper instability at low temperature. For simplicity we restrict ourselves to s-wave superconductivity [59], which is also reported in the experiments on paraelectric STO [35].
Finally, writing the Hamiltonian in BdG form we ob-tainĤ is the Nambu spinor, ∆ = iσ y ∆ in the s-wave BCS channel, and we choose a gauge in which ∆ is real. The BdG Hamiltonian above enjoys a particlehole symmetry, implemented by P = τ x C, where τ j , j = x, y, z are Pauli matrices in the particle-hole space and C is the complex conjugation operator. Namely, the Hamiltonian obeys PH BdG (k)P † = −H BdG (−k). Additionally, when B P the Hamiltonian has rotational symmetry about the axis parallel to the polarization, where The Fermi surface of the free Rashba gas (in the ferroelectric phase). Blue and green arrows denote spin texture of the outer and inner sheet of the FS, respectively. (b) In the presence of non-zero magnetic field parallel to the polarization, the Fermi surface develops a gap at the points of overlap of the two sheets of the Fermi surface of the Rashba gas. Blue and green arrows denote the direction of the spins corresponding to the momenta on the FS lying on kz-axis. (c) Depairing effect of the Zeeman field is maximal along kz, and sufficiently strong magnetic field destroys superconductivity locally in momentum space, making the spectrum gapless at the specific momenta along kz -the Weyl points denoted by the red and blue spheres, which colors signify positive and negative chiralities, respectively. (d) Bogoliubov Fermi surfaces appear when the magnetic field is perpendicular to the polarization direction.
the rotation includes both spatial and spin rotation. In the presence of higher order terms due to the lattice, the continuous rotational symmetry is reduced to discrete four-fold rotations about the polarization axis. The energy dispersion is determined from the solutions of a quartic equation [see Eq. (A1)], which for a magnetic field parallel to the polarization yields where k ⊥ = (k x , k y ) denotes the projection of the momentum onto the xy-plane. The 3D Fermi surfaces of the free Rashba gas described by Eq. (2) have the shape obtained by rotating two displaced circles around the axis connecting their crossing points (see Fig. 1a). Consequently, the crossings form pinching points along the k z axis, where two Fermi sheets with opposite helicities touch. Upon turning on a magnetic field in the z-direction, the two sheets separate, and the spins at these points becomes co-linear with the field direction. Thus, the depairing effect of the magnetic field in the superconducting phase is expected to be strongest at these pinching points.
Indeed, a sufficiently strong magnetic field closes the gap at the pinching points on the k z axis. From Eq. (4), we see that the gap closes for B 2 > ∆ 2 at momenta p = (0, 0, p z ), where This equation is satisfied at four points with j = 1, . . . , 4 labeled in descending order along the k z -axis (see Fig. 1c). The closing of the gap at these momenta can be viewed as a topological phase transition in the two-dimensional Hamiltonian H BdG (p x , p y , p z ), where p z is a tuning parameter. Indeed, for p 2 < p z < p 1 and p 4 < p z < p 3 the two dimensional Bloch bands have non-zero Chern numbers ±1 (of equal sign), signaling that the Weyl nodes are monopoles of Berry charge. It is worth noting that in the low density limit there are only two Weyl nodes p 1 and p 4 , in accord with the finding of previous literature [18][19][20][21].
Rotation of B with respect to the ferroelectric moment P profoundly changes the quasiparticle spectrum. Due to the rotational symmetry, the dispersion is symmetric for both k → −k and E → −E separately, when B P . However, in the presence of a perpendicular component, the spectrum is invariant only under the combined action of these two operations. This means that when the angle is large enough, the Weyl cones over tilt and become type II [60,61], which is accompanied by the development of the Fermi surface of zero-energy Bogoliubov quasiparticles [51,52] (see Fig. 1d). This mechanism is analogous to the one described in Ref. [62] for the surface of a 3D topological insulator and 2DEG Rashba spin-orbit gases with the proximity induced superconductivity and applied inplane magnetic field. For more details see Appendix A.
To make these observations more concrete, we derive the low-energy effective Hamiltonian in the vicinity of the Weyl nodes by projecting to the low-energy subspace. This yields the 2 × 2 Hamiltonian where ), and all other components of the matrix A are equal to zero. The chiralities of the Weyl nodes are determined by and are controlled by B z , which is the projection of the magnetic field, B, on the polarization vector P . The σ 0 -term in Eq. (7), which is proportional to the components of B that are perpendicular to P , is responsible for tilting the Weyl cones when the magnetic field and polarization are not collinear. This can be seen by noting the energy spectrum of the Hamiltonian Eq. (7) As mentioned above, the system can even be driven into a type-II phase, where the cones tilt is so strong they dip below the Fermi energy and form Bogoliubov Fermi surfaces [51]. The condition for Bogoliubov Fermi surfaces to develop is the existence of non-zero k such that Using the expressions in Eq. (8), we find that this criterion is satisfied when B 2 ⊥ > ∆ 2 . Close to the cone, the Bogoliubov Fermi surface sheet defined by (k) = 0 from Eq. (10) is a cone with the opening angle in k x k y -plane φ = π − 2 arcsin ∆ B ⊥ . However, inspecting the full Hamiltonian Eq. (3) (see Appendix A), we find that, in fact, the Bogoliubov Fermi surfaces form the shape of two bananas touching at the Weyl points (see Fig. 1d).
Before proceeding to the physical consequences of the Weyl nodes, we comment that in our model they appear exactly at zero energy. This is however, not fixed by symmetry, but an artifact of the gap function we chose, which is purely the A 1g representation (s-wave). The inversion symmetry breaking renders this representation indistinguishable from A 2u (p z , which is triplet). Therefore the gap is in general a mixture of the two, which is characterized by nodes shifted from zero energy, where the sign of the shift for each node depends on the sign of the momentum along z. Such a shift will inflate the nodes leading to small Bogoliubov Fermi surfaces (see Appendix B).
We finally note that the angle between P and B can be spatially manipulated, for example across a domain wall separating different ferroelectric domains. This opens a path to control the Weyl nodes, as we discuss in the following section.

III. FERMI ARCS ON SURFACES AND DOMAIN-WALLS
In this section we discuss the Majorana Fermi arcs, which appear on surfaces and domain walls. We first review the well known case of an interface between a single domain and vacuum. We then turn to the case of internal tetragonal domain walls.

A. Majorana arcs on the surface of a single domain
We first show that Majorana arc states appear on the boundary between a single domain and the vacuum. Assuming that the ferroelectric moment P is tilted with an angle θ to the interface, we pick a coordinate system such that the yz-plane is in the plane of the interface, the z-axis aligns with the projection of P onto the interface, and the x-axis directs into the domain. The Hamiltonian for the domain is given by Eq.
ing open boundary conditions: The Bogoliubov quasiparticles operators are defined as Thus, the reality conditionγ We now look for the solution of Eq. (11a) in the form Ψ k || (x) = Ψ 0,k || e −αx , where Re(α) > 0 implies decaying solutions as x → ∞. Plugging this into Eq. (11a), we obtain the characteristic equation for α [ee Appendix C, Eq. (C1)], the solutions of which for the parallel momentum k || , denoted α k || , obey α k || = α * −k || in accord with the reality condition Eq. (13). Analysis shows (see Appendix C) that for |B y | < ∆ there are four roots with positive real part. In this case, a general decaying solution for Eq. (11a) is a linear combination of four solutions: . Plugging this into the boundary condition Eq. (11b) and requiring vanishing of the determinant of the resulting set of the linear equations with respect to coefficients C i , one obtains k || for which a non-trivial solution, corresponding to the Majorana-Fermi arc, exists. For |B y | > ∆, when the Weyl cones overtilt in x-direction, we do not find Fermi arcs on the x = 0 surface.
It is easy to find analytical solution for B y = 0. In this case, it is expected that Majorana-Fermi arcs are formed at k y = 0. Indeed, in this case Eq. (C1) for α splits into two simpler ones where η = ±1, and we find v kz,s = ηu kz,s , as required by particle-hole symmetry. Thus, the problem separates into two sectors corresponding to η = ±1. The number of roots in the right half-plane depends on the sign of the quantity Defining a new (primed) coordinate system, rotated such that its k z -axis aligns with the ferroelectric moment, k z = k z cos θ, we see that Ξ < 0 is just the condition for the momenta k z to lie between the two Weyl nodes p 1 and p 2 or p 3 and p 4 . Precisely, for Ξ < 0 (Ξ > 0), there are three (two) roots in the right half-plane for η = 1, and one (two) roots for η = −1. The Dirichlet boundary condition Eq. (11a) and the normalization of the wave-function define three conditions to be satisfied. Thus, for η = 1 for momenta on k z -axis lying between the projections of two nearby Weyl nodes, a non-trivial solution corresponding to Majorana-Fermi arc exists, see the dashed lines in Fig. 3a. In Appendix C, we show that for non-zero B x and |B y | < ∆, the Majorana-Fermi arcs remain to be straight lines connecting the projections of the Weyl nodes.

B. Majorana arcs on domain walls
In the previous subsection we showed that Majorana zero modes (MZMs) connecting into Fermi arcs appear on the boundary with vacuum. We now turn to discuss another situation relevant to experiments in STO: Domain walls between different tetragonal domains. To understand the nature of such domain walls, we recall that low-temperature STO has spontaneously broken its cubic symmetry into tetragonal structure. In this phase each oxygen octahedra rotates about one of the three cubic axis, clockwise or anticlockwise, alternating from unit cell to unit cell [35,63], which is known as antiferrodistortive (AFD) order. The rotation axis fixes the polarization direction when tuning into the ferroelectric phase.
The AFD phase is notoriously known to breakout in domains [66,67], which appear in two types, one endows the system with the reflection symmetry about the wall and the other endows the system with the reflection symmetry about the wall combined with a glide [68]. The AFD order parameters in neighbouring domains constitute ± π 2 angle with each other. In turn, the ferroelectric polarization in the neighbouring domains will also differ by direction with a relative angle of π 3 or − 2π 3 , see Fig. 2. We fix the polarization vector in the first domain to be A 4 ( Fig. 2). When the polarization vector in the second domain is B 2 or B 4 , the Weyl nodes coincide when projected onto the momentum plane parallel to the wall (see Fig. 3a). In contrast, if the polarization vector in the second domain is B 1 or B 3 , the projections of the Weyl nodes from the two domains are at different points ( Figs. 3b and 3c). Below we present a qualitative description of the resulting Fermi arcs for these scenarios.
FIG. 2. Interface between the two AFD domains, D1 and D2, with ferroelectric orders. AFD1,2 -direction of the AFD distortion in the left (red) and right (blue)) domains, respectively; Ai, Bi -possible directions of the ferroelectric moments in the left (red) and right (blue)) domains, respectively.
In both scenarios, the effective low-energy Hamiltonian is given by [69,70] where 1,2 (k || ) are the low-energy chiral modes of each of the domains, D 1 and D 2 , and the off-diagonal matrix component a(k || ) are the couplings. The eigenvalues of Eq. (16) are given by x 2 − ( 1 + 2 )x + 1 2 − |a| 2 = 0 and therefore, the Fermi arc states obey the equation (i) The scenario in which the projections of the Weyl points onto the interface of both domains coincide-This happens when the polarization vector in D 1 is A 4 , and the polarization vector in D 2 is B 2 or B 4 . We assume the magnetic field B lies in xz-plane for simplicity. We then identify two cases: case I-The chiralities of the Weyl nodes with coinciding projections are the same. In this case 2 = − 1 , and the condition becomes − 2 1 = |a| 2 , which can be satisfied only when |a| 2 = 0 for k || at which 1 = 0. However, there is no symmetry that fixes a(k) = 0 for k on that line. Therefore, the arcs are gapped out in the general case. In Appendix C, we discuss such unprotected zero energy solutions. case II-The chiralities of the Weyl nodes with coinciding projections are opposite. Here 1 = 2 . Consequently, the arcs are robust and found on the lines for which 1 (k || ) = ± a(k || ) (see Fig. 3a).
An important consequence of the scenario of coinciding Weyl points when projected to the domain wall, is that a rotation of the magnetic field B about the y-axis allows to continuously tune between case I and case II. Then we expect arc states to disappear and reappear as a function of angle.
(ii) The scenario where projections of the Weyl nodes do not coincide-This happens when the polarization vector in D 1 is A 4 , and the polarization vector in D 2 is B 1 or B 3 . For ∆ 2 < B 2 < ∆ 2 + µ 2 the Majorana Weyl arcs will "repel" and "attract" each other as shematically illustrated in Fig. 3b. For B 2 > ∆ 2 + µ 2 , a more signifi-cant reconstruction of the Majorana Fermi arcs happen. For the point close to the crossing point, we can write which defines a hyperbola in the vicinity of the crossing point, now connecting the projections of the Weyl nodes of the same chirality (see Fig. 3c).

IV. WEYL-SUPERCONDUCTIVITY IN THE PRESENCE OF VORTICES
Up to this point we have only considered the Zeeman coupling to the magnetic field. We now turn to consider the consequence of the orbital coupling. In a type-II superconductor, the field can induce vortices when it exceeds the value H c1 . We distinguish two limits of interest. In the small magnetic field limit H c1 < B H c2 the distance between vortices is much greater than the coherence length and each vortex can be treated independently. In the opposite limit, B H c2 the vortices become densely packed, overlap and significantly reduce the global average value of the order parameter.
In what follows, we focus on these two limits. We start with the single vortex problem. Using the results of Section II, we show that individual vortices in ferroelectric superconductors can contain non-trivial Majorana bound states, even when the bulk superconducting state is trivial. Then in the next step, based on semiclassical considerations (namely, assuming strong localization of the Majorana states on a scale much smaller than the coherence length), we find that there is always a critical magnetic field B * < H c2 , marking a percolation transition to a putative topological state with Majorana-Weyl nodes in the bulk.

A. The single vortex problem -Non-trivial bound states
In the solution of the Ginzburg-Landau equations for a single vortex, the superconducting order parameter ∆(r) and magnetic field B(r) both depend on the radial distance from the vortex core. Starting from the core and moving outwards, the order parameter is initially zero, and adjusts back to its bulk value at a distance of the order of the coherence length ξ. The magnetic field, on the other hand, is maximal at the core and gradually decays to zero at a distance given by the penetration depth λ L (we assume that λ L ξ). The dependence of these two fields is schematically plotted in Fig. 4.
In light of the discussion in Section II, this implies that somewhere between the vortex core and r → ∞ there is a "halo" radius r h , where the critical threshold for creating Majorana-Weyl nodes B(r h ) = ∆(r h ) is satisfied (see Fig. 4). Majorana-arc states then appear on a cylinder of radius r h and at the core of the vortex. Clearly, such states can only be observed if their localization length l M is significantly smaller than r h .
To obtain these states we solve the BdG equation explicitly (see Appendix D). We consider two models. First we consider a toy model, which we solve analytically. In this model B is taken to be constant and we mimic the spatial dependence of the gap near the vortex core by breaking it into two steps (see dashed lines in Fig. 4). Namely, the core region is defined to be in the region r < r 1 , where the gap is zero. The second region is the topological "halo" defined in the region r 1 < r < r 2 (where r 2 is the halo radius r h in this model). In this region the gap takes a non-zero value ∆ 1 , which is smaller than the field, such that the topological criterion ∆ 1 < B is satisfied and there are Weyl nodes. The third region is r > r 2 , where we assume ∆(r > r 2 ) = ∆ 2 such that ∆ 2 > B and therefore the superconducting state is trivial and fully gapped.
The explicit solution shows there are two exponentially localized Majorana bound states, which are slightly split in energy due to the finite spatial separation between the boundaries at r 1 and r 2 . The key result we obtain from the toy model is an estimate of the localization length of these states where ξ 0 = v F π∆2 is the Pippard coherence length estimated at the momentum k z located between the Weyl nodes.
As can be seen, the length scale Eq. (19) appears in units of ξ 0 and is proportional to the parameter λ/v F . Close to H c2 , the halo size becomes of the order of the correlation length. Therefore the Majorana arc states on the edge of the Halo can be resolved from the core Majorana states in the limit λ v F . Recent theoretical results estimate the electron coupling to the transverse optical phonon mode in STO [57]. Using the average displacement in the ferroelectric phase [39,45], this coupling constant can be estimated to be λ = 254 meV· • A. Using this value of λ, we find the concentration at which v F becomes greater than λ (which happens when the Fermi surface crosses the Dirac point) is n ≈ 2.3·10 19 cm −3 . For higher densities, the ratio λ/v F diminishes. For reference, this parameter diminishes to λ/v F = 1/5 at n ≈ 1.3 · 10 21 cm −3 . It is worth noting however, that other estimates of λ are smaller [54]. To confirm the results of the toy model we also solve the BdG problem numerically using a more realistic profile of the gap, ∆(r) = ∆ 0 tanh(r/ξ), where ∆ 0 = exp(iφ)|∆ 0 | and |∆ 0 | > B. We solve the BdG problem inside the interior of a cylinder of radius R. As before, the topological criterion ∆(r) < B is only satisfied within a finite halo radius r h surrounding the core. The resulting amplitude of one of the two BdG wave functions with nearly zero energy is shown in Fig. 5. We observe two peaks, corresponding to location of the core and the critical radius r h .
An interesting aspect of the halo is that it realizes a local pseudo magnetic field [71]. The continuous variation of |B(r) − ∆(r)|, which controls the distance between the Weyl nodes, therefore acts as a pseudo gauge field in the z direction A z (r). The resulting pseudo magnetic field looks like a vortex circulating the core of the halo. An important physical consequence of this field is the emergence of a whole spectrum of Landau levels, which in this case are labeled by angular momentum. These states are plotted in Fig. 10 in Appendix D. For more details regarding the analytic and numeric solutions of the BdG problem we refer the reader to Appendix D.

B. Local tunneling density of states in the vicinity of a single vortex
Using our results from the previous subsection, we now compute the local tunneling density of states in the vicinity of a vortex. The resolution of a typical scanning tunneling microscope is much smaller than the size of the vortex, and therefore it may be capable of distinguishing the core and edge states described above. The local √ 2mµ (such that pz = 0) as a function of energy and distance from the vortex core; b) νp z (E, r) (normalized by νp z ,0, the peak value of νp z (E, r) in the plotted region) for pz = 1.05 √ 2mµ as a function of energy and distance from the vortex core. The parameters used for the simulations are: m = 1, ∆0 = 2, Bz = 1.71, λ = 0.1, µ = 10, ξ = 100, R = 700, and the temperature T = 5 · 10 −5 ∆0. density of states, which is often proportional to the differential conductance [72,73], is given by ν(E, r) = dp z ν pz (E, r) Here Ψ σ,i (r) is the radial part of the σ-component of the Nambu wavefunction corresponding to the i-th eigenmode of energy E i , ν pz (E, r) stands for the contribution to the local density of states from eigenmodes corresponding to a particular p z , and we substituted deltafunctions with the negative derivatives of the Fermi-Dirac distribution at low temperature. In Fig. 6a, we plot ν pz (E, r) for p z = √ 2mµ, as a function of energy and distance from the vortex core r. One can clearly see peaks at zero bias for r = 0 and r ≈ r h . For p z 's away from √ 2mµ (but for which MZMs still exist), the distance between the peaks of the MZMs deacreases, while the localization length of MZMs increases. This results in broadening of the peaks in r-direction and further separation in Edirection; see Fig. 6b [74] in which ν pz (E, r) is plotted for p z = 1.05 √ 2mµ. Consequently, the full density of states ν(E, r) (and the differential conductance) will have smeared zero-bias peaks.
C. The many-vortex problem -Percolation of the topological phase The picture presented above, where each isolated vortex is surrounded by a topological halo, suggests the possibility of a percolation transition, where the halos overlap and the topological phase percolates through the system. At large field, B ∼ H c2 , the ground state of the system is expected to be a vortex lattice. Therefore, let us assume the magnetic field is large enough such that the lattice state is formed, yet the halos are still separated and each vortex is encircled by chiral Majorana zero modes. Upon increasing the field even further, the halos grow, and eventually touch, creating a connected sea of the topological phase. Below we develop a crude estimate for this percolation threshold B * and find that it is always smaller than H c2 . Full microscopic calculations are required to verify this scenario more rigorously, especially considering that the proof for the existence of Majorana zero modes presented in this manuscript is strictly valid only at low magnetic fields.
We use Abrikosov's theory [75], applicable for magnetic fields close to the upper critical field H c2 . The harmonic approximation solution for the first Ginzburg-Landau equation can be written in the form where ∆ 0 is the gap function at zero magnetic field and where ξ(T ) is the coherence length at temperature T , x n is a position of nth vortex core on x-axis, 2π/q is the periodicity in the y-direction, and D n are dimensionless coefficients. Substituting Eq. (21) into the second Ginzburg-Landau equation, one finds an expression for the magnetic field where κ = λ L /ξ is the Ginzburg-Landau parameter. We note that Abrikosov's theory is valid for small f 2 , therefore in what follows we will consider B 0 0.9H c2 .
To find the next order correction to f in the small parameter 1 − B 0 /H c2 one requires that [75] where O stands for the averaging O over one unit cell of the vortex lattice. Then, using a parameter which characterizes a lattice structure (for the square lattice β = 1.18, and for the triangular one β = 1.16), the non-zero solution for f 2 is The spatial profile f (r) of the order parameter is defined by the coefficients D n in Eq. (21). For simplicity, in the following we consider the case of the square lattice, for which D n are constants denoted by D, and x n = nqξ 2 .
Percolation of the topological phase will occur when the magnetic field at the half distance between the neighboring vortices' cores, d/2, reaches the critical value for the topological phase transition (see Fig. 7), B(d/2) = ∆(d/2), i.e., where f 0 is given by Eq. (22) in which all D n set to one. An important remark here is that we use the topological criterion B > ∆ derived in Section II for a uniform case and which remains valid for the vortex problem in a small magnetic field when the vector potential terms in BdG Hamiltonian can be neglected. In principle, in high magnetic fields, the contribution from these terms might modify the topological criterion. We leave the investigation of this question for future work. Combining this equation with Eq. (26), we find the value of B * needed to be applied to reach the percolation point where The field B * at which the topological phase percolates is therefore controlled by two phenomenological parameters. The first is K = ρ s /2∆ 2 0 N (0)ξ 2 0 , where ρ s is the superfluid stiffness and N (0) is the density of states of the underlying metal. Assuming a full volume fraction, a parabolic band dispersion ρ s = 2 n/4m [76] and ∆ 0 ≈ 1.76T c [77], we have K ≈ 0.1(µ/T c ) 2 /(k F ξ 0 ) 2 , which can be estimated directly from experiment (at n = 10 18 cm −3 µ = 2 meV, T c = 200 mK [78]) to (b) At the critical value B * , previously bounded topologically non-trivial "puddles" touch with the neighbouring "puddles" at one point.
(c) For large values of the magnetic field, topologically non-trivial "puddles" around each vortex overlap, creating a topologically non-trivial "sea".
(d) In the more realistic disordered network of vortices, the percolation will have a disordered character as well. Red lines are Majorana "halos" surrounding the topological phase.
FIG. 7. Schematic illustration of the percolation of the topological phase. The black line structure is a contour plot for ∆((r)), orange color denotes topologically trivial regions, and blue color -topologically non-trivial regions. be between 1 and 5 depending on the value of ξ 0 between 100 and 50 nm, respectively. It is interesting to compare this result with the prediction of BCS theory K = 0.5 [79]. The second parameter controlling B * is the ratio ∆ 0 /H c2 . Comparing with the experimental data of Ref. [80] we find that this parameter can be on the order of (and even larger than) 1. We plot the resulting schematic [81] phase diagrams in the space of magnetic field B and temperature T for different values of ∆ 0 /H c2 and K in Fig. 8, where we naively extended the use of Eq. (28) beyond the region of validity of the Ginzburg-Landau theory. The value of δ in Eq. (28) does not depend much on κ for κ > 3, and we fix it to equal to 10. As can be seen, for all values of K there is a topological phase separating between the trivial superconducting and normal state. This result is much more generic than our particular model. We predict that any noncentrosymmetric superconductor where inversion is broken by a vector [53] will develop such topological halos above H c1 . Consequently, all such superconductors may undergo a percolation transition to a bulk topological phase before giving way to the normal state.

V. CONCLUSIONS AND DISCUSSION
We studied Majorana-Weyl superconductivity emerging in systems with intertwined superconducting and ferroelectric orders due to the application of a magnetic field. First, we considered the effect of a uniform Zeeman field. We confirmed that above the Clogston-Chandrasekhar threshold gµ B B > 2∆, Weyl cones appear in the Bogoliubov quasiparticle spectrum along the axis of the polarization moment, regardless of the charge density. We also showed that rotating the magnetic field with respect to the polarization tilts the Weyl cones and eventually causes Bogoliubov Fermi surfaces shaped as bananas to appear.
However, the magnetic field is not expected to be uniform in the superconducting state. Instead it threads through the sample in the form of vortices. Due to the vanishing of the gap at the core of each vortex, the critical threshold gµ B B > 2∆ is always fulfilled in some area surrounding it, which we dub the "halo". Such halos are characterized by Majorana strings at their core and chiral Majorana arc states going around them. When the magnetic field is increased towards H c2 the vortices become denser, the halos merge and the system undergoes a percolation type phase transition to a bulk Majorana-Weyl superconductivity. This transition always precedes H c2 .
Our predictions have a number of sharp experimental consequences. The first is the emergence of topological halos surrounding vortices at small magnetic fields above H c1 . These can be observed in the local tunneling density of states using an STM. However, we expect a clear separation of scales between the size of the halo and the arc state's localization length, only close to H c2 . This is because the magnetic field at the center of an isolated vortex is of order H c1 , which is much smaller than the critical threshold. Therefore, the halo radius is very small when the magnetic field is far from H c2 . In addition to the zero modes, the nodes also modify the tunneling density of states away from zero energy. Namely, due to the bulk nodes there will be a quadratic dependence on bias. The arc and nodal states can also be observed in the heat conductivity. For example, we anticipate that close to H c2 , in the topological phase, the system will become heat conducting albeit still superconducting. Finally, when tilting the magnetic field to be perpendicular to the polarization direction we expect Bogoliubov Fermi surfaces to emerge. Close to H c2 these surfaces will contribute a T -linear term to the specific heat and a constant tunneling density of states. Finally, it is also possible that the existence of Majorana zero modes surrounding vortices will contribute a constant term to the specific heat close to H c2 , which will manifest itself as a Schottky anomaly at low temperatures. The size of the anomaly should diminish by a factor of 1/2 when crossing to the topological phase.
Full quantum calculations are needed to verify the proposed scenario of the percolation of the topological phase.
In the presence of slowly varying disorder, naively, we may anticipate a scenario similar to the transitions between integer quantum Hall states [82], where the topological nature of the phases is manifested as long as a network of edge-states percolates through the bulk. Furthermore, it is also interesting to consider the transition between the topological state considered here and the FFLO state, which is also a relevant ground state when the magnetic field is perpendicualr to the polarization [12][13][14]. To that end, one needs to solve selfconsistently for the lowest energy ground state. We postpone the study of such questions to future work. The energy dispersion of Eq. (3) is determined from the equation For B x = B y = 0, its solutions are easily found and given in Eq. (4). Here we analyze its zero-energy solution for the case when the magnetic field is not parallel to the polarization. We first show that the conditions for the gap closure is essentially the same as for the case of perpendicular magnetic field. For a generic quartic equation a product of its roots Considering Eq. (A1) for k ⊥ = 0, e = 0 is satisfied at B 2 = ∆ 2 + 2 k , signifying that there is a zero root. In addition, this root is double, as it can be readily seen form Eq. (A1): the free term is zero, and the linear term is zero at k ⊥ as well. Thus, the gap closes at k ⊥ = 0 for B 2 > ∆ 2 at k z determined from the equation B 2 = ∆ 2 + 2 k . Other non-degenerate zero solutions might be determined from equation e = 0, where e is the free term in Eq. (A1). Without loss of generality, choosing the direction of the magnetic field such that B y = 0, we recast this equation in a form which determines the dependence of k y on k ⊥ and k z for the momenta satisfying the condition e = 0. It is indeed a solution if |k y | ≤ k ⊥ , which may happen only if B 2 x > ∆ 2 . Also, note that these roots are non-degenerate, and there are roots of different sign amongst those four corresponding to the solution of Eq. (4). It is easy to see that the solution with k ⊥ = 0, ∆ 2 + 2 k −B 2 = 0 (i.e. away from the Weyl nodes) is impossible. Thus, for B 2 x > ∆ 2 , we infer that the momenta at which E = 0 form closed surface(s) defining a 3D Bogoliubov Fermi surface. Numerical investigation shows that these surfaces connect the Weyl nodes at p 1,2 and p 3,4 , respectively (see Fig. 1d). In particular, for B 2 x = ∆ 2 , B z = 0, zero solution can exist only for k x = 0, and from Eq. (A3) we obtain which defines two intersecting circles We emphasize that this result is obtained under the assumption that the superconducting order parameter remains s-wave order. Here we illustrate that the presence of a triplet component (k z dependent) leads to the inflation of the Weyl nodes into Bogoliubov Fermi surfaces. We consider Eq. (3) with ∆ = iσ y (∆ + ∆ 1 k z σ z ) treating ∆ and ∆ 1 as parameters. In Fig. 9, we plot zero-energy surfaces in k-space in the parallel (a) and not parallel (and overtilted) (b) to P magnetic field for the case of B > ∆ ∆ 1 . We note that overtilting of the magnetic field produces large Bogoliubov Fermi surfaces.
Appendix C: Additions to the "Fermi arcs on surfaces and Domain-walls" section of the main text The characteristic equation for α in the ansatz solution )λ cos θk y − iB z λ(sin θk z − i cos θα) − λ 2 sin 2 θk 2 y (B x + B y ) 2 = 0.
One may view this equation as an equation with real coefficients with respect to iα. Thus, its roots are sym- metric with respect to the imaginary axis, i.e., if α is a root, then −α * is also. Therefore, Eq. (C1) can have four roots with positive real part.
In the main text, we showed the existence of the Majorna-Fermi arcs for the case of B || = 0. Here, we present a solution for an arbitrary B. To find a locus of Majorana zero modes in the k y k z -plane by substituting the general solution of Eq. (11a) into boundary conditions Eq. (11b) without any assumption is quite difficult. Instead, we check if k y = 0, k z ∈ (p z2 , p z1 ) (p z4 , p z3 ) is the locus of zero-energy solutions.
For k y = 0 at artibitrary B, as for the case of B y = 0, the characteristic equation for α, Eq. (C1), splits into two simpler equations where η = ±1, and we find Again, the problem separates into two sectors corresponding to η = ±1, and the further analysis proceeds in analogy to the presented one in the main text.
In the main text, based on the low-energy theory, we pointed out that non-protected Fermi arcs still may exist in the case II of scenario (i), where the Weyl nodes of the same chiralities in two domains project onto same points on the interface. Here we show that such solution exists in our continuous model.
We choose the coordinate system as described in the main text for the case of a boundary between the single domain and vacuum with x−axis pointing into the first domain, D 1 . The boundary problem to be solved is where H 1(2) and Ψ 1(2),k || (x) are the Hamiltonian and the wavefunction for the first (second) domain.
Again, we are looking the solutions for k y = 0 in the form Ψ 1(2),k || (x) = Ψ 01(2),k || e −α 1 (2) x , where α 1 (2) are determined from Eq. (C2). We consider λ > 0 in the first domain, and the flip of the Weyl node chiralities in the second domain correspond to the flip of the sign of λ, i.e., λ < 0 in D 2 . The decaying solutions in D 1(2) imply α 1(2) > (<)0. For k z ∈ (p z2 , p z1 ) (p z4 , p z3 ), there are three α 1 with Re(α 1 ) > 0 and three α 2 with Re(α 2 ) < 0 in case of η = 1. For η = −1, there are one α 1 with Re(α 1 ) > 0 and one α 2 with Re(α 2 ) < 0. For k z ∈ (p z2 , p z1 ) (p z4 , p z3 ) there are two α 1 (2) with Re(α 1(2) ) > (< 0). Again, boundary condition Eq. (C4c) imply that we can stitch solutions in D 1 and D 2 corresponding to the same η = ±1 only, and that the problem separates into two sectors η = ±1. This results in that that the boundary conditions effectively give us four constraints, and together with the normalization condition there are five constraints. For k z ∈ (p z2 , p z1 ) (p z4 , p z3 ), the general solutions for Eqs. (C4a) and (C4b) are linear combinations of three functions, which gives us six unknown coefficients to be found. This is one more then the number of constraints we have, which implies that we get a family of solutions parametrized by one parameter, which might be thought of as an angle in two-dimensional vector space. Thus, this set of solutions can be thought as a linear combination of two orthogonal solutions. To obtain MZM in the presence of vortices in the low field regime, we solve the BdG problem in the vicinity of a single vortex. The corresponding BdG Hamiltonian in cylindrical coordinates assumes the form where the z-component of the momentum remains a good quantum number and Here pz = p 2 z 2m − µ and we have neglected the coupling to the vector potential [83]. The phase of the order parameter winds by 2π around the vortex origin, ∆(r) = ∆ 0 (r)e iφ . The cylindrical form of Eq. (D1) suggests to look for the energy eigenstates in the form (D3) Following Ref. [8], in searching for the Majorana modes, we focus on the l = 0 channel, which is also justified by our numerical calculations. Substitution of Eq. (D3) into Eq. (D1) leads to a system of ordinary differential equations (ODE) with real coefficients, and thus the functions Ψ iEl (r), i = 1..4 in Eq. (D3) are real. For l = 0, the particle-hole symmetry implies σ x ⊗ σ 0 Ψ pzE0 (r) * = ηΨ −pz−E0 (r), where η is a phasefactor. Given that the Hamiltonian Eq. (D1) is even in p z , we obtain σ x ⊗ σ 0 Ψ pzE0 (r) * = ηΨ pz−E0 (r). Combining this with the statement about the reality of Ψ i (r), we conclude η = ±1 and Ψ 3(4)E0 = ηΨ 1(2)−E0 . Although we anticipate the splitting in energy due to overlapping of Majorana states at r = 0 and r = r h , we start with seeking the zero-energy solution. In the following, we drop the subindices for the putative E = 0, l = 0 state and use Ψ i ≡ Ψ i00 , and thus we have Ψ 3(4) = ηΨ 1 (2) . Then the zero-energy eigenstate equation for BDG Eq. (D1) reduces to the system of two ODEs where Ψ(r) = (Ψ 1 (r), Ψ 2 (r)) T .
In what follows, we first present an analytic analysis of Eq. (D4) for a simplified piece-wise constant model. Then in the next step we present numerical analysis for more realistic profiles of the gap and magnetic field, which continuously vary in space.
The gap structure in the simplified model constitutes of three regions We also assume the magnetic field is uniform (justified by the type II condition λ L ξ). We then focus on the limit ∆ 1 < B z , in which case the intermediate region is "topological". In the region 0 < r < r 1 , where ∆(r) ≡ 0, we look for the solution in the form [8] where J n (z) are Bessel functions. Substituting this into Eq. (D4), we find a characteristic equation for α which has four solutions: ±α 1 and ±α 2 . Thus, the general solution in this region is For the regions r 1 < r < r 2 and r > r 2 , where ∆ 0 (r) = 0, we look for the solution in the form [8] and get a set of algebraic equations for the coefficients (a n , b n ). For n = 0, we obtain where ∆ stands for either ∆ 1 or ∆ 2 depending on the region under consideration. This gives the following equation forq = −iq The rootsq i of this equation satisfy the condition For the region r > r 2 , the decaying solutions correspond to suchq that Re(q) > 0. For 4 i=1q i > 0, there are two such roots for either η; for 4 i=1q i < 0, there are three such roots for η = −1, and one such root for η = 1.
Two boundaries (at r 1 and r 2 ) with smooth continuity conditions for a two-component vector and one normalization condition bring nine conditions, in total. Now we count the number of yet unknown coefficients in the constructed solution to be obtained from these conditions focusing on the case 4 i=1q i > 0 for the region r > r 2 for all p z (which correspond to ∆ 2 > B z ), where the middle region (the halo) r 1 < r < r 2 is in the topological phase, while the outer and inner regions are trivial. In the region 0 ≤ r < r 1 , there are two coefficients; in the region r 1 ≤ r < r 2 , there are four coefficients; and in the region r > r 2 , there are two coefficients. This brings in total eight coefficients, which is not enough to satisfy nine conditions. Thus, there is no zero-energy solutions. In fact, this is anticipated and corresponds to the overlapping of two Majorana states at r 1 and r 2 . Specifically, removing the "domain" wall at r 2 (or moving it to infinity), at the boundary r = r 1 we have to satisfy only five conditions at r = r 1 . In this case, for r > r 1 we have to single out only decaying at infinity solutions, which gives three coefficients in this region. And it total we have five coefficients to satisfy five conditions. Analogously, moving r 2 to infinity, and requiring that physical solutions decay far away from r 2 in the topological phase, we look for suchq in Eq. (D11) that Re(q) < 0. There is one such root for η = −1 sector, and three such roots for η = 1. Then, in η = 1 sector we again have equal number of constraints and coefficients. As we move r 2 from large distance closer to r 1 , the overlapping of the two Majorana modes leads to splitting in energy of the states constructed out of the linear combinations of these Majoranas. Alternatively, we can think in the following way. For the case B z > ∆ 2 , we have one Majorana zero mode localized around the vortex core, which is essentially the case considered in Ref. [8]. But as we decrease B z to values just below ∆ 2 , we lose the zero-energy solution. The only way it can happen is via pairing the Majorana zero mode at the vortex core with another one at r = r 2 . We also note that while here we considered a crudely discretized model, the argument presented extends to the arbitrary fine discretization. Indeed, the introduction of a new segment brings in four new boundary conditions and, at the same time, four new constants to be found, thus leaving the balance between the number of conditions and the number of coefficients untouched.
However, because the halo has finite size it is essential to estimate the splitting of the zero modes due to their overlap. To this end, at zero temperature, we assume the separation between the two boundaries r 2 − r 1 ∼ ξ 0 is on the order of the coherence length ξ 0 = v F π∆2 [84]. This length should be compared with the localization length of the zero modes, l M , which can be estimated from the low-energy effective Hamiltonian Eq. (7) (for B = Bẑ) H ef f = vk y σ x + vk y σ y + E g σ z . (D12) where E g is the gap at a given k z away from the Weyl point. We then find that l M ∼ v Eg . Comparing with Eq. (7), we find that v = λ∆ B , E g = pz pkz 2mB yielding We then evaluate l M for k z located at the middle point between the two Weyl points under assumption µ √ B 2 − ∆ 2 . Focusing on the region r 1 < r < r 2 , we find . Thus, the ratio of the length scales is controlled by the small parameter λ/v F and is therefore expected to be very small except for very close to the nodes or close to the transition point.
To confirm our analytical considerations, we perform numerical calculations. For simplicity, we consider a cylinder of a radius R with a single vortex located at the axis of the cylinder and impose zero boundary conditions at r = R. For r < R, we assume a radial dependence of the order parameter given by the function ∆(r) = ∆ 0 tanh(r/ξ), which reflects a typical behavior in a vortex core center. Also, we assume the magnetic field is uniform, and smaller than the bulk threshold B z < ∆ 0 .
We represent the radial part of the spinor Ψ pzEl (r) in the Bessel-Fourier series form where µ l i is the set of roots of the equation J l (µ l i ) = 0, which guarantees that the boundary conditions are satisfied. Substituting this representation into Eq. (D1) and projecting onto J ν (µ i r R ), ν = l − 1, l, l + 1, we obtain an infinite system of algebraic equations, which is solved approximately by truncation. In the calculations used for producing plots in this article, we cut the system of algebraic equations at size 600 × 600.
We plot eigenenergies corresponding to the wavefunctions Ψ El (r, θ) in Fig. 10. Under PHS, l → −l and E → −E. Thus, in fact, for l = 0, there are two near-zero energy solutions (for the parameters considered, these energies are on the order of 10 −4 · ∆ 0 ) that are indistinguishable in the plot and correspond to the states that are linear combinations of Majorana zero modes.