Yu-Shiba-Rusinov qubit

Magnetic impurities in $s$-wave superconductors lead to spin-polarized Yu-Shiba-Rusinov (YSR) in-gap states. Chains of magnetic impurities offer one of the most viable routes for the realization of Majorana bound states which hold a promise for topological quantum computing. However, this ambitious goal looks distant since no quantum coherent degrees of freedom have yet been identified in these systems. To fill this gap we propose an effective two-level system, a YSR qubit, stemming from two nearby impurities. Using a time-dependent wave-function approach, we derive an effective Hamiltonian describing the YSR qubit evolution as a function of distance between the impurity spins, their relative orientations, and their dynamics. We show that the YSR qubit can be controlled and read out using the state-of-the-art experimental techniques for manipulation of the spins. Finally, we address the effect of the spin noises on the coherence properties of the YSR qubit, and show a robust behaviour for a wide range of experimentally relevant parameters. Looking forward, the YSR qubit could facilitate the implementation of a universal set of quantum gates in hybrid systems where they are coupled to topological Majorana qubits.


I. INTRODUCTION
The goal of building a fault-tolerant quantum computer has allowed the understanding of the quantum realm in a plethora of systems to be deepened, as well as leading to an advancement in the development of novel quantum technologies. Trapped ions, semiconductor quantum dots, superconducting circuits, and hybrid semiconductorsuperconductor platforms are some of the examples which have played crucial roles in developing the field of quantum computing [1][2][3][4][5][6][7][8][9][10][11][12][13][14].
While superconducting-circuit-based qubits have been at the forefront of the immense recent progress, proposals that utilize the low-energy bound states in superconductors (SCs), i.e., the Andreev levels, have also been under intense scrutiny for quantum computing [15][16][17][18][19][20][21]. The reasons are twofold: (i) the dimensions of Andreev-statesbased qubits (approximately micrometers) are typically much smaller than the sizes of the conventional superconducting qubits (approximately millimeters), which facilitates the design of quantum registers with higher qubit densities; and (ii) they constitute the building blocks of topological quantum computers based on Majorana zero modes, which have experienced significant theoretical and experimental research efforts [22][23][24][25][26][27].
Magnetic impurities in superconductors lead to localized Yu-Shiba-Rusinov (YSR) in-gap Andreev states [28][29][30][31][32][33][34], with chains and lattices of impurities being viable setups to realize topological superconductors hosting the Majorana modes . The advantage of these implementations is rooted in the ability to pattern superconducting surfaces with magnetic impurities and possibly engineer (topological) quantum processors in a controlled fashion. Moreover, through the use of scanning tunneling microscopy (STM) techniques, they can be interrogated locally, with high spatial resolution. A drawback, however, is that the system parameters are hard to tune, making it difficult to control the topological regime of the system or to manipulate the emerging Majorana modes. Several solutions have been put forward, among which are exploiting the dynamics of the magnetic impurities [59,60], driving the YSR states with microwave fields [61], varying the orientation of the external magnetic fields [62,63], or tuning the Josephson effect through a superconducting tip coupled to the YSR states [64].
The realization of the Majorana-based topological quantum computer in Shiba chains looks distant, as no experimental evidence of quantum degrees of freedom exists as yet in these systems. For this purpose, it would be necessary to experimentally demonstrate that it is possible to coherently manipulate the Majorana qubits before they decohere. In this paper, we show that the minimal system for the demonstration of the quantumness of these systems is a new type of superconducting qubit, the YSR qubit, stemming from two nearby impurities. We demonstrate that the dynamics of the magnetic impurities can be used for controlling the quantum state of the YSR qubit and we uncover the requirements for experimentally observing Rabi oscillations in this system. The precession of the magnetic impurities also leads to a feedback torque acting on the impurities due to the YSR states [65] and we show that this effect can be utilized for the read out of the YSRqubit states. We also address the effect of the spin noises on the coherence properties of the YSR qubit and show robust behavior for a wide range of experimentally relevant parameters. Our proposal is feasible with state-of-the-art experimental techniques, because controlled coupling of YSR states in impurity dimers has already been experimentally demonstrated [66][67][68][69][70] and the manipulation of the impurity spins is possible through STM electron spin resonance (STM-ESR) techniques [71][72][73][74][75]. Finally, we discuss the possibilities of utilizing the YSR qubits in hybrid systems where they are coupled to Majorana qubits.
The paper is organized as follows. In Sec. II, we introduce the model Hamiltonian describing the dynamical spin dimer. Using a time-dependent wave-function approach, in Sec. III we derive the effective YSR-qubit Hamiltonian in the presence of the precessing spins. In Sec. IV, we discuss how to implement coherent Rabi oscillations of the YSR qubit and provide a specific manipulation protocol. Then, in Sec. V, we demonstrate that the spin dynamics can be utilized for the readout of the YSR qubit. In Sec. VI, we introduce a hybrid YSR qubit -Majorana (topological) qubit that can be operated to achieve a universal set of quantum gates. We conclude with a discussion in Sec. VII.

II. MODEL HAMILTONIAN
The time-dependent Bogolioubov de Gennes (BdG) Hamiltonian describing the spin dimer system in Fig. 1 can be written in the Nambu basis FIG. 1. The YSR qubit: two classical spins, target (red) and test (blue), are placed on top of a two-dimensional (2D) s-wave superconductor at a distance R, inducing a double-well potential that accommodates two in-gap YSR states for a given parity. The odd-parity states |0 and |1 define the YSR-qubit states. The asymmetry of the potential stems from the slightly different coupling parameters at the two sites. The hybridization of the two YSR states is quantified by a tunneling Hamiltonian H T . Driving the test spin effectively tunes the potential bias and, when at resonance with the qubit splitting 2 = q , it allows for coherent rotations of the qubit. The target spin is interrogated offresonantly at a frequency 1 = q for quantum non-demolition detection of the YSR-qubit state. Alternatively, the qubit can also be operated with just a single tip.
where H 0 is the superconductor Hamiltonian and V j (t) describes the coupling of electrons to the classical spins S j (t) = S[sin θ j (t) cos φ j (t), sin θ j (t) sin φ j (t), cos θ j (t)] with time-dependent polar θ j (t) and azimuthal φ j (t) angles (j = 1, 2). Here, R j = 0(R) is the position of the spin j = 1 (j = 2), J j are the coupling strengths, σ = (σ x , σ y , σ z ) and τ = (τ x , τ y , τ z ) are the Pauli matrices in the spin and particle-hole spaces, is the superconducting order parameter, and p = p 2 /2m − μ is the kinetic energy of the electrons with effective mass m, momentum p, and chemical potential μ. For simplicity, we neglect the scalar potentials [56] generated by the magnetic impurities, as they do not affect directly the dynamics. The target spin S 1 and test spin S 2 can be addressed and driven individually through STM ESR and they are used for read out and manipulation, respectively.
Before proceeding with a detailed description of the dynamics, let us provide some physical insights into the spin dimer in Fig. 1 based on the recent findings in Ref. [65] concerning the dynamics of a single magnetic impurity in an s-wave SC. We first note that for the static case, the Shiba energy is given by E S = (1 − α 2 )/(1 + α 2 ), with α = πν 0 JS and ν 0 being the density of states at the Fermi level in the normal state. For a spin precessing with frequency (adiabatic limit) at an angle θ around the z axis, the effective Shiba energy has been found to be E S ( ) E S − ( /2) cos θ , i.e., it is shifted by the Berry-phase contribution. Moreover, this dynamical YSR state has been found to act back on the classical spin via a universal torque τ S (t) = −(n S − 1/2)F Sṅ (t), where n(t) = S(t)/S, n S is the YSR state occupation number and F S is the radial Berry curvature. In the absence of spin-orbit interaction, F S = 1/2. This torque modifies the bare resonance frequency 0 of the classical spin as δ / 0 ≈ (1/S)(n S − 1/2)F S , being a direct measurement of the occupation n S . The energy shift (on the Shiba side) and the frequency shift (on the classical spin side) are at the core of our proposal, which is depicted in Fig. 1: the former allows us to control the bias of the double-well potential by driving one of the spins, analogously to tuning the voltage bias in double quantum dots [76], while the latter facilitates extraction of the occupation of the in-gap states, in analogy to quantum non-demolition qubit readouts in cavity quantum electrodynamics setups [77]. The hybridization between the two YSR states will modify the single-impurity findings and in the following we proceed to describe the dynamical YSR dimer system in detail.

III. EFFECTIVE QUBIT HAMILTONIAN
Next, we derive the low-energy Hamiltonian describing the in-gap "molecular" YSR states stemming from the dynamical spin dimer, using a time-dependent wavefunction approach that allows us to identify the effective two-level system defining the YSR qubit. The system dynamics are described by the time-dependent BdG equation T is the BdG wave function. It is instructive to switch to the Fourier space ψ(r, t) = (1/L d ) k e −ik·r ψ(k, t) (where d is the dimension of the system), which in turn allows us to write Assuming that the Shiba energies are close to the Fermi level (deep Shiba limit α 1,2 ≈ 1) and adiabatic dynamics of the classical spins on the scale of T = / , we can follow the approach described in Refs. [41,60] to derive an effective time-dependent 8 × 8 Schrödinger equation is a four-component spinor at position r = 0 (r = R) in the Nambu and spin space, while the diagonal elements H 11 (22) the interaction of the SC with the spins at the positions r = 0(R), and H 12 (t) = H T (t) represents the tunnelling between the YSR states at different impurities: Here, n i = S i (t)/S i andĨ 0,1 (R) are evaluated from the overlap integrals for two impurities separated by a distance R in the superconductor (Appendix A). We point out that the time dependence of the classical spins generates a Berry-phase contribution [the second term in Eq. (3)] that cannot be captured by only making the effective static theory time dependent. As shown later, while this term does not affect the qubit Hamiltonian, it does drastically change the spin expectation values at each impurity; in particular, the contributions perpendicular to the instantaneous classical spin directions (which are responsible to the torques acting on the latter). This is one of the instances when an effective static theory does not suffice to describe the low energy sector dynamics. Let us first consider the dimer in the absence of dynamics. Projecting the above 8 × 8 Hamiltonian blocks, H i and H T , onto the low energy sector results in an effective 4 × 4 Hamiltonian describing the in-gap states [41,60] (for details, see Appendix A). Assuming θ 1 = 0 (i.e., the first spin defines the z axis), the in-gap energy spectrum of the 4 × 4 Hamiltonian becomes where t h = 4α 1 α 2 e −R / √ 2π k F R quantifies the tunneling strength and k F is the Fermi momentum. Note that all energies are expressed in terms of = 1, while all lengths are expressed in terms of the SC coherence length ξ = v F / , with v F being the Fermi velocity. In Fig. 2(a), we show the corresponding energy spectrum as a function of θ. The inset of Fig. 2(a) depicts the relative maximum deviation in the energy difference between the lowest two energy states as a function of δα = α 1 − α 2 for various separation distances R. From these plots, we can infer that: (i) even for moderate values of δα, the dependence of the energies E i on θ is negligible; and (ii) in general, these energies are not equidistant. We can then encode the YSR qubit in the two lowest-energy states defined by the {−E 1 , −E 2 }, which, for a wide range of parameters, are also well separated from the excited pair {E 1 , E 2 }. The qubit Hamiltonian can be written as Here, 2 is the classical spin precession frequency that matches the qubit splitting. The inset shows the energy width δE = q (0) − q (π ) / q (0) + q (π ) as a function of δα = α 1 − α 2 and for several values of R. (b) The many-body energy spectrum E MB for α 1 = 1.15, α 2 = 1.1 and R = 2.9 in the absence of Coulomb interactions. The dotted lines represent the manybody energy spectrum for α 1 = α 2 = 1.1, showing a crossing at θ = π . Here, {|10 , |01 } and {|00 , |11 } label the odd-and evenparity states, respectively, with the YSR qubit being encoded in the former. In both plots, we use k F = 13.55.
It is useful to describe the YSR qubit using manybody states |n 1 n 2 , where n 1,2 = 0, 1 are the occupancy of the single-quasiparticle states. Specifically, the pair of states {|00 , |11 } ({|01 , |10 }) spans the even-(odd-) parity many-body states with energies ±|E 1 + E 2 |/2 (±|E 1 − E 2 |/2). Note that within the BdG description the two parity sectors are decoupled and the YSR qubit defined above acts within the odd-parity states. This choice for the YSR qubit is further justified by its insensitivity to the Coulomb interaction effects that are present for double occupancy (even parity). The many-body energy spectrum depicting the odd-and even-parity states is shown in Fig. 2(b). While the ground state corresponds to the even-parity state |11 for the chosen parameters, the odd-parity sector can be selected by tuning the offset charge with a gate voltage in the case of a finite superconducting island with a sufficiently large charging energy. Alternatively, one can utilize the spin dynamics for the initialization of the system to the odd-parity state.
The many-body picture also allows us to gain further insight into the origin of the qubit states, which is determined by max(δα, t h ). For δα t h , the qubit states stem from the two individual YSR states formed under each of the impurities, while in the opposite regime δα t h , they correspond to the symmetric and antisymmetric superposition of the individual YSR states, being dictated by the tunneling. The first scenario is more advantageous, as the qubit energies become insensitive to θ, as depicted in Fig. 2, rendering it more robust against fluctuations.
Having defined the YSR qubit, we can now reinstate the dynamics of the classical spins, which we exploit for the manipulation and read out of the qubit states. Without loss of generality, in the following we assume that only one spin precesses. Projecting the 4 × 4 time-dependent Hamiltonian onto the YSR-qubit subspace, we obtain the following qubit Hamiltonian (Appendix B): where Equations (6) and (7) establish the imprints of the classicalspin dynamics on the effective YSR-qubit Hamiltonian and represent one of our main findings. Above, we disregard the terms that act as identity in the qubit space. The first two terms in Eq. (7) induce transitions between the qubit states, while the last term allows us to dynamically control the qubit splitting q → q + 2β z (t).
is independent of any of the microscopic parameters.

IV. YSR-QUBIT MANIPULATION
The YSR qubit can be manipulated by utilizing the second term in Eq. (6). The pulse sequence for introducing Rabi oscillations is shown schematically in Fig. 3(a). The logical states of the qubit are defined in a parallel classicalspin alignment and the resonant oscillations between the states of the qubit are induced in the antiparallel configuration. Before describing the details of the sequence, let us underline the physical reasons for this choice. The β i (t) terms in Eq. (7) are much weaker for deviations δθ around θ = 0 [β x,z ∝ (δθ ) 2 and β y ∝ δθ] than when the same deviations occur in the proximity of θ = π (β x ∝ δθ and β y,z ∝ constant), which makes them rather inefficient in the parallel configuration. In the idle phase, on the other hand, this is beneficial, since the qubit will be more robust against random fluctuations in the angles θ and φ (discussed below). Nevertheless, the qubit can also be operated fully in the antiparallel geometry, at the expense of shorter coherence times.
Let |0 and |1 be the eigenstates of the static qubit Hamiltonian at θ = 0 and let us assume that the qubit is initialized in state |0 at time t = 0. Then, at time t, the qubit state becomes |ψ(t) = U q (t, 0)|0 , where the evolution operator is U q (t, 0) = Te −(i/ ) t 0 dt H q (t ) , with T being the time-ordering operator. In step 1 of the protocol in Fig. 3(a), the right classical spin is rotated from 040347-4 the parallel to the antiparallel configuration via a pulse θ(t) = π tanh(2π t/T a ), where T a is the pulse length, and the evolution operator is U q,1 ≡ U q (T a , 0) [78]. The amplitude of the Rabi oscillations is largest if the qubit remains in state |0 during this pulse. Thus, ideal results are obtained if 0 − π transition is adiabatic, i.e., T a T q with T q = / q , but almost ideal Rabi oscillations can also be achieved for fast pulses (Appendix E). In the second part of the sequence, the classical spin is driven into circular precession around the z axis, so that the precession frequency which first stabilizes the precession of the spin to a cone angle θ = π − δθ in a time T b (step 2 ), then causes a precession of the spin for a duration T l − 2T b (step 3 ), and finally restores the classical spin back to θ = π in time T b (step 4 ). Assuming that T b T q implies that the evolution induced byφ(t) during the ramping periods is practically frozen and we can write U q,2 ≈ U −1 q,2g U q,2r U q,2g , where U q,2g and U q,2r correspond to the evolution from π to π − δθ atφ = 0 and the evolution induced by the circular precession at fixed δθ (θ = 0) during the time T l − 2T b , respectively. Finally, the classical spin is rotated back to the parallel configuration using U q,3 ≡ U −1 q,1 , shown by step 5 in Fig. 3(a).
The amplitude and the period of the Rabi oscillations can be determined by calculating how the probability for the qubit to be in state |1 after a pulse, , depends on the precession time T l . We implement numerically the evolution operator pertaining to U q (t) and in Fig. 3 showing the Rabi oscillations of the qubit for the parameters δα = 0.05, δθ = 0.1, R = 2.9. An increase in the precession angle δθ increases the Rabi oscillation frequency as R ∝ t h δθ but in turn reduces its amplitude, as depicted in Figs. 3(c) and 3(d). The latter is a consequence of the transformation U q,2g , which generates a finite weight c 1 ∝ (t h /δα)δθ on the state |1 for t h δα. Therefore, for a given R , the requirement for P 1m ≡ max[P 1 (T l )] ≈ 1 is R q , which, coincidentally, is similar to the adiabaticity condition in the first part of the protocol. Moreover, in the limit T b T q , the transformation U q,2g is purely geometrical (Appendix E) and thus independent of the details of the pulse that tilts the classical spin away from the z axis by an angle δθ. As stressed above, the manipulation can be fully performed in the antiparallel configuration, in which case U q (T l , 0) ≡ U q,2 . Both the parallel and the antiparallel configurations have been observed experimentally, their realization depending on the specific implementation and the distance between the impurities [68,79].
To give some estimates for the time scale of the Rabi oscillations, let us assume that δθ = 0.1, R = 2.9 and δα = 0.05. These rather conservative parameter values result in the Rabi oscillation period T R ≈ 8.5 ns, which is comparable to the Rabi times observed in implementations of the Andreev qubits [17]. For the YSR qubit to be useful, the Rabi time should be much shorter than the time scales over which it loses its coherence, namely the relaxation (T 1 ) and pure dephasing (T φ ) times, which, to the best of our knowledge, are largely unknown for the YSR states. Nevertheless, we can readily identify several possible sources of decoherence: (i) quasiparticle poisoning [80][81][82], (ii) thermal fluctuations in the magnetic moments (magnons) that define the YSR states, and (iii) phonon or photon coupling to the Shiba electrons [83]. Decoherence induced by nonequilibrium quasiparticle poisoning is highly specific to the system and thus it is difficult 040347-5 to provide precise scalings and estimates. Recent studies, both experimental and theoretical, show that the relaxation times pertaining to this mechanism can range from milliseconds to even seconds [80][81][82]. The general consensus is that their effect can be minimized by improving the samples and that it can be accounted for by a phenomenological linewidth of the isolated YSR states, which in principle can be extracted from STM-ESR measurements in the limit of weak tunnel coupling [84,85]. The last two mechanisms, on the other hand, have not been discussed in the literature for the YSR molecule. In the following, we give a short account of the magnon-induced decoherence, while the details of the phonon (and photon) mechanism is described in Appendices F and G. The Hamiltonian describing the coupling of the qubit to the magnetization fluctuations of spins k = 1, 2 reads as follows: where the tensor χ k ≡ [χ μν k ] quantifies the coupling of the two orthogonal fluctuations μ = 1, 2 [δn k (t) ⊥ n k ] of each classical spin k to the qubit Pauli matrices ν = x, y, z. The elements of the tensor χ μν k can be found by projecting ∂ n k H BdG onto the qubit basis (Appendix F). Within the Bloch-Redfield framework [85], we find the following expressions for the dephasing and relaxation times, respectively: 1 where χ μσ is the noise spectrum pertaining to the fluctuations δn k,μ (t). Above, the σ = +(−) terms represent the emission (absorption) rates, which are related by the detailed balance condition at equilibrium. The spectrum of the fluctuations δn k (t) is determined by the specific form of the classical spin free energy F S (n) and in order to give estimates for the above decoherence times, we consider the following form (assuming the free energies of the two spins to be identical): where κ measures the crystal anisotropy (intrinsic or induced by the surface), B is the externally applied magnetic field, and γ is the gyromagnetic ratio. Considering κ > 0, this free energy per spin is consistent with the perpendicular to the surface configurations observed in experiments. At finite temperatures B → B + δB(t), with δB(t) being the stochastic contribution the Fourier components of which satisfy the fluctuation-dissipation relations δB μ (ω)δB ν (ω ) = [87], in which α g and γ are the Gilbert damping and gyromagnetic coefficient, respectively. By utilizing the Landau-Liftshitz-Gilbert (LLG) equation that describes the dynamics of classical magnets in the presence of stochastic magnetic fields δB(t), we can evaluate the correlators S k μν (ω) in terms of δB μ (ω)δB ν (ω ) (for more details, see Appendix F). For simplicity, we focus only on the static (idle) parallel-and antiparallel-spin configurations, assuming a spin S = 5/2 at each site [88,89]. In both cases, we find that the pure dephasing rate is zero and, furthermore, for θ = 0, the relaxation rate is also zero, quantitatively justifying our choice for the qubit basis in the idle phase. However, at θ = π , the longitudinal relaxation rate is nonzero and assuming that α g = 0.001, κ = 0.1 meV, α 1 = 1.15, α 2 = 1.1, R = 2.9, and temperature T 0 = 100 mK < q , we obtain T 1,m ≈ 3.5 μs [90]. Comparing that to the Rabi oscillation period, we estimate that the YSR qubit can undergo a large number of Rabi oscillations before it decoheres due to magnons.
We find that both the phonon and photon couplings vanish in the antiparallel configuration (where the YSRQ is operated), in stark contrast to the magnons, which have their maximal effect. That is because neither phonons nor photons can induce spin flips, which are required for quasiparticle tunneling between the two YSR states in this configuration. In the parallel arrangement, the phononinduced relaxation is instead maximal and we evaluate it to be T 1,ph ≈ 5.8 μs. This is a slightly longer time than the coherence time induced by the noise in the magnetic moments. However, all these sources of decoherence seem to be of similar magnitude and they are also similar to the coherence time observed in the Andreev qubits [17].

V. YSR-QUBIT READOUT
The ability to measure the outcome of a computation efficiently and rapidly is a prerequisite for a practical qubit. Furthermore, it allows us to initialize the qubit state at the beginning of the computation. Here, we show that the qubit state can be measured using STM-ESR techniques via the torques induced by the YSR states on the classical spins. In the following, we focus on the case when the measurement is performed in the parallel spin configuration and the left spin (target) is interrogated off-resonantly with the qubit splitting, as shown in Fig. 1. The former condition is considered in order to minimize the decoherence effects, while the latter allows us to physically separate the manipulation and detection.
The dynamics of the left spin, S 1 , are governed by the LLG equatioṅ 040347-6 where τ = −J 1 S 1 × σ δ(r) is the torque pertaining to the electrons in the SC that act on spin S 1 , including the YSR-qubit contribution, while B(t) is the time-dependent external magnetic field utilized to drive the precession. We employ a Green's function approach that describes the 8 × 8 Hamiltonian [60] to evaluate the total torque τ σ for the two YSR-qubit states σ = 0, 1 (for more details, see Appendix C). Considering S 1 to precess with frequency 1 in the adiabatic limit 1 , we can write τ σ ≈ τ σ s + τ σ d , where the first term (τ σ s ) originates from the misalignment of the two classical spins and describes the contribution of the in-gap states to the Ruderman-Kittel-Kasuya-Yoshida (RKKY) interaction, while the latter (τ σ d ∝ 1 ) has been unraveled recently in Ref. [65] and has been found to have a geometrical (Berry-phase) origin. In Fig. 4(a), we show the magnitudes of the total torque τ σ , as well as the two individual contributions τ σ s and τ σ d , as a function of δα for each of the two qubit states. We see that the torques are determined by the static term τ σ s in the limit 1 |δα| t 2 h , while in the opposite regime, 1 |δα| t 2 h , the dynamical contribution τ σ d dominates and, moreover, it reaches a universal value associated with an isolated impurity [65]. We need to mention that throughout the section we neglect the effect of the bulk states on both the static and dynamical torques. In Ref. [91], it has been shown that the (static) bulk contribution, which represents the conventional RKKY interaction, becomes negligible compared to that of the YSR in-gap states for separations R ≥ 1. Furthermore, a full nonequilibrium calculation for a single impurity has shown that the YSR states dominate the dynamical torque in the deep Shiba adiabatic regime and, moreover, that a finite YSR linewidth imprints onto the magnetic impurity linewidth [65], which could be utilized to measure the coherence times of the YSR qubit.
In the limit of small cone-angle precession (θ ∼ 0), we can linearize the LLG equation and extract the renormalized resonance frequency of spin S 1 for each qubit state in terms of the torques as (Appendix D) where , and 0 is the bare resonance frequency. The difference δ = r,0 − r,1 discriminates between the two qubit states in the STM-ESR measurements and represents one of our main findings. In Fig. 4(b), we plot δ as a function of δα for various distances R. For separations R such that t 2 h 1 |δα|, the difference δ saturates to a constant value, which we find to be δ ≈ 8S 0 /(16S 2 − 1) [65]. That is because the two impurities become practically decoupled, resulting in τ σ s → 0, and only the dynamical torque from the isolated impurity contributes to the signal. We see again here that the optimal regime for operating the YSR qubit is when tunneling between the two isolated YSR states is smaller than their energy difference, in which case δ is almost invariable for wide range of system parameters.
To enrich the understanding of the above results, we present a heuristic derivation of the YSR-qubit torques from basic energy considerations. In the readout regimė θ = 0 and 1 =φ = q , so that the effective qubit splitting is eff q = q + 2β z and we can neglect the β x,y terms in Eq. (7). Then, the magnitude of the torque acting on the spin S 1 by the YSR qubit in state σ = 0, 1 can be expressed as We see that for t 2 h 1 |δα|, the dynamical torque dominates, reaching a universal value τ d ≈ ( 1 /4) sin θ, consistent with the findings in Fig. 4(a). On the other hand, for t 2 h 1 |δα|, the torque is controlled by the static contribution τ s and reaches the asymptotic value τ s ≈ −(t h /4α 1 α 2 ) sin(k F R + π/4) sin (θ/2), which depends strongly on the separation R between the impurities. This is again consistent with the findings in Fig. 4(a). Interestingly, while for t 2 h 1 |δα| the torque τ d behaves similarly around θ = π , for α 1 = α 2 ≡ α, we obtain τ s ≈ −(t h /4α 2 ) sin(k F R + π/4), which is a consequence of the two qubit levels crossing each other: even though the spins 040347-7 are antiparallel, a torque is exerted between the two, which is to be contrasted with the RKKY interactions mediated by the bulk [91]. This can also be interpreted as a fractional spin Josephson effect that is protected by the presence of the inversion symmetry.
To give estimates for the possible frequency shifts δ , let us consider the following experimentally pertinent values: α 1 = 1.15, α 2 = 1.1, R = 2.9, = 1 meV ≈ 241 GHz, and S = 5/2. Interrogating the classical spin with frequencies 0 ∼ 25 GHz results in δ ≈ 4.9 GHz, which is well within the state-of-the-art experimental resolution [74]. Note that with the above parameters, the qubit splitting q ≈ 9.5 GHz and thus the target spin precesses off-resonantly, which is essential for noninvasive readout of the qubit.

VI. YSR-MAJORANA HYBRID QUBIT
The YSR qubits described here could be utilized for quantum information tasks on their own; for example, by creating a network of weakly interacting spin dimers on top of superconductors that can be addressed individually. More importantly, they could be integrated with Majorana zero modes hosted at the ends of spin chains in superconductors and exploited for performing universal quantum computation. Indeed, the braiding statistics of the Majorana zero modes alone are not sufficient for implementing the universal set of topological gate operations necessary for quantum computation and additional nontopological gates are needed to achieve universality. A viable way to implement the missing π/8 phase gate is to control the couplings of the Majorana zero modes [92], e.g., by varying the magnetic fluxes in the transmon geometries [23,93,94], and extremely robust geometric [95] and distillation [96] protocols can be utilized if sufficiently accurate control of the couplings is possible. However, it is not easy to realize suitable pulses to control the couplings of the Majorana modes in Shiba chains and to our knowledge no proposals currently exist for robust protocols to implement the π/8 gate in these systems. An alternative proposal to implement the π/8 gate is to integrate the Majorana zero modes with quantum-dot-based qubits [24,[97][98][99] and here we show that this idea can be transferred to the context of Shiba chains by utilizing the YSR qubit.
In Fig. 5, we sketch our proposal for generating a universal set of gates in a hybrid YSR-qubit (YSRQ) -Majorana-qubit (MQ) system. A T-junction formed by three Shiba chains generated by magnetic impurities (or adatoms) placed on top of an s-wave superconductor interacts via tunnelings t l,r with a spin dimer that encodes a YSRQ described in the previous sections. The left and right Shiba chains are topological, hosting Majorana zero modes γ † i = γ i (i = 1, . . . , 4), while the lower one is in the trivial regime. Two Majorana zero modes on the left (right) chain define a fermionic state, which can be occupied and fusing the Majorana zero modes. Driving the dimer locally rotates the YSRQ coherently (see Sec. IV) and by combining this rotation with the SWAP gates facilitated by t l,r (described in detail in Ref. [99]), the missing π/8 phase gate can be implemented. The SWAP gate can also be used for the measurement of the MQ via the readout of the YSRQ. or empty. Thus, the states of the MQ can be encoded either as {|0 l 0 r , |1 l 1 r } (in the case of even parity) or as {|0 l 1 r , |1 l 0 r } (odd parity).
The hybrid qubit can be operated as suggested in Ref. [99]. The Clifford gates can be implemented in a topologically protected fashion by fusing Majorana zero modes and by braiding them in a T-junction shown in Fig. 5 [22]. Additionally, the nontopological π/8 gate can be performed by first swapping the MQ state to YSRQ, then performing the π/8 gate on the YSRQ, and finally swapping the YSRQ state back onto the MQ. The SWAP gate can be constructed by utilizing the general ideas presented in Ref. [99]. The tunable couplings t l,r lead to an interaction Hamiltonian between the qubits where M j (j = x, y, z) are the Pauli matrices acting on the MQ and J ij (t) are time-dependent coupling strengths originating from t l , t r = 0. In the case of q = 0 (a degenerate YSRQ qubit), a specific sequence of operations for implementing the SWAP gate has been provided in Ref. [99]. In our case, q = 0 can be achieved easily in the antiparallel configuration by applying a magnetic field B z along the 040347-8 z axis. Indeed, for θ = π , the tunneling vanishes and we obtain that q = 0 at gμ B B z = 2δα, with g and μ B being the g factor and the Bohr magneton, respectively. Moreover, since we are assuming the deep Shiba limit (δα 1), the condition gμ B B z 1 is satisfied, which means that the bulk SC remains unaffected. Alternatively, it might be possible to utilize the dynamics of the classical spin in the implementation of the SWAP gate. The SWAP gate can also be used for the initialization and measurement of the MQ via the readout of the YSRQ.

VII. DISCUSSION AND CONCLUSIONS
In this work, we introduce and study a novel type of quantum bit, the YSR qubit, that is encoded in the energy states of a spin dimer coupled to an s-wave superconductor. We demonstrate theoretically that both the coherent manipulation and the readout of the YSR qubit can be efficiently implemented by harnessing the dynamics of the spins that engender it. Furthermore, we scrutinize the effect of the classical spin fluctuations on the coherence of the YSR qubit and show robust behavior compared to the manipulation times. Given the ability to manipulate magnetic adatoms on superconducting substrates with a high degree of control, the YSR qubit could be utilized together with Majorana topological qubits to facilitate the performance of universal quantum computation. We propose one such hybrid implementation that is based on topological Shiba chains and a spin dimer hosting the YSR qubit.
There are several avenues for future studies. An immediate objective would be to generalize the time-dependent formalism described here to account for the spin-orbit effects originating from both the substrate [69] and the anisotropy of the exchange coupling between the classical spins and the superconducting electrons [68]. The spin-orbit coupling effects can provide the microscopic mechanism for the easy-axis anisotropy and therefore they are potentially useful for the operation of the YSR qubit but, additionally, they can stabilize the ferromagnetic order in adatom chains and facilitate the realization of the Majorana modes. Thus, these effects are crucial for the operation of the hybrid qubit.
Another important direction is to establish the Hamiltonian of the hybrid qubit from the microscopics in order to study the possible quantum gates and to optimally engineer the adatom deposition. Additionally, coherent transfer of the information between the two types of qubits might be beneficial for entangling MQs that are separated by a large distance. The YSRQs could be entangled, for example, by utilizing cavity quantum electrodynamics in a similar fashion to the way it has been employed in various other solid-state qubits [77]. While the coupling of the YSR qubit to the magnetic field of a microwave cavity should be weak (<kHz), the stronger electric field component could instead couple to the qubit by affecting the tunneling t h or via the spin-orbit coupling. Alternatively, these interactions could be ignited indirectly, via the coupling of the quantum fluctuations of the classical spins to the microwave photons [51,100], which could also be used to drive them into precession. Moreover, such setups would naturally allow for the YSRQ to interact with other types of qubits and enhance their functionality.
Further down the road, it would also be interesting to extend the dynamical framework describing the spin dimers to Shiba chains and 2D Shiba islands that can host Majorana end modes and chiral Majorana edges [58,101], respectively. We believe that by triggering the magnetic dynamics, it should be possible to both manipulate and detect the Majorana edge modes, the latter by leaving their fingerprints on the STM-ESR signals. Moreover, such approaches should be advantageous, as they would allow us to interrogate Shiba systems with wellestablished methods from spintronics [102]. As a catalyst for another direction of future work, we speculate that the in-gap Shiba states, either in dimers or chains, could mediate out-of-equilibrium spin interactions in the presence of external magnetic drives that have no counterparts in the static situations. Given the nonperturbative nature of the coupling between the electrons in the superconductor and the spins, which converts into torques as in Eq. (15), that would require the self-consistent solution of the combined dynamics of the two systems. This might result in novel spin configurations that are stabilized dynamically and, on the electronic side, induce new types of (possibly dissipative) phases [103].
In conclusion, the YSR qubit proposed in this work operates well within the current experimental capabilities and we expect it to open up new possibilities for future studies on superconducting systems patterned with spins. We hope that our findings help to create a road map toward a functional MQ in these systems.

ACKNOWLEDGMENTS
We would like to thank Thore Posske for the interesting and fruitful discussions. The work is supported by the Foundation for Polish Science through the International Research Agendas (IRA) program cofinanced by the European Union within Smart Growth Programme.

APPENDIX A: DERIVATION OF THE LOW-ENERGY HAMILTONIAN
In this appendix, we show the derivation of the Hamiltonian describing the YSR qubit given by Eq. (6). The total SC Hamiltonian written in the Nambu basis (r) = 040347-9 where, as described in the main text [Eq. (1)], the BdG Hamiltonian is  ψ(r, t). By using the Fourier decomposition ψ(r, ψ(k, t), we can recast the BdG equation in the following form: which, in the static limit, pertains to the substitution i∂ t → E and coincides with the equation for the spectrum presented in Ref. [41]. In the frequency domain, and retaining only the leading-order terms in ω, we obtain where In the above expressions, we retain only the leading-order corrections in ω, assuming that the time dynamics of the classical spins as well as those of the emerging Shiba energies are such that ω (adiabatic regime). For a 2D superconductor, the integrals I 0,1 (r) can be written as I 0 (r) = (2ν 0 / )Ĩ 0 (r) and I 1 (r) = 2ν 0Ĩ 1 (r) [41], wherẽ ν 0 is the density of states, k F is the Fermi momentum, v F is the Fermi velocity, and K 0 is a modified Bessel function of the second kind. At r = 0, this givesĨ 0 = π/2 andĨ 1 = 0, while for k F r 1 (a limit utilized throughout our work), the asymptotic expressions are as follows [60]: Next, we switch back to the time domain and we obtain t), we can manipulate this expression further by writing the combined evolution as where and n i ≡ n i (t) = S i (t)/S i . As mentioned in the main text, the dynamics of the spins induce an extra term in the local Hamiltonian matrix element of Berry-phase origin. Without this term, the transverse spin expectation values at the positions of the impurities would have the wrong sign.
To help distinguish the low-and high-energy sectors, which in turn allows us to eliminate perturbatively the terms that couple them, it is instructive first to perform a unitary transformation U 0 = exp (iπτ y /4) that converts τ x ↔ τ z , followed by a (time-dependent) U i (t) = τ 0 ⊗ U i (t) that acts on site i = 0(R) and diagonalizes the terms ∝ n i · σ : These rotations affect the terms in Eq. (A11) and they become The low (high) 4 × 4 energy sector is spanned by the σ z τ z = −1(1) and the corresponding energies of the isolated Shiba states are

040347-10
Consequently, we can then project the remaining terms, i.e., the tunneling and the velocity contributions ∝φ i ,θ i , onto the low-energy sector to obtain an effective timedependent Hamiltonian.
To simplify the discussion, from here onward we assume the left spin (1) is static and aligned along the z direction, or θ 1 , φ 1 = 0, and that θ 2 ≡ θ and φ 2 ≡ φ. Furthermore, we also consider that all lengths are expressed in terms of ξ SC = v F / and set = 1. Then, the projected 4 × 4 lowenergy Hamiltonian can be written as where P l is the corresponding projector, while represents the instantaneous projected Hamiltonian, with A l,φ (t) and A l,θ (t) being the projected gauge field terms associated with theφ andθ contributions in Eq. (A14).

APPENDIX B: EFFECTIVE QUBIT HAMILTONIAN
We can further diagonalize the instantaneous 4 × 4 Hamiltonian H l,0 (t) in order to identify the effective 2 × 2 YSRqubit Hamiltonian presented in the main text. This is achieved by another time-dependent unitary transformation, with s 1 = sign[sin(k F R + π/4)] and s 2 = sign[cos(k F R + π/4)]. Its effect on H l (t) can be formally written as with s = φ, θ . Note that whileH l,0 (t) is now diagonal, with energies ±E 1,2 , where E 1,2 = (B ± C)/(2α 1 α 2 ), the gauge field terms can induce transitions between its eigenstates. More importantly, these terms have two separate contributions: one from the initial gauge fields, originating from the first unitary transformations U 1,2 , and one from the diagonalization of the 4 × 4 effective (instantaneous) time-dependent Hamiltonian H l,0 (t). Both are required to correctly capture the low-energy-sector dynamics and 040347-11 starting from the static effective theory by making the parameters θ and φ time dependent would lead to erroneous results.
In order to establish a qubit that is well separated from the excited states, we assume that both α 1 , α 2 = 1, as well as t h , δα 2 min |1 − α 1,2 |, with δα = |α 1 − α 2 |. Then, we can further projectH l (t) to the two lowest-energy states, resulting in the qubit Hamiltonian presented in the main text (up to terms that act as identity in this subspace): where q = C/α 1 α 2 is the qubit splitting energy.

APPENDIX C: DETAILS ON THE READOUT VIA TORQUES
Here, we provide details on the calculation of the torque τ = −J 1 Sn 1 × σ (0) acting on the precessing spin S 1 by the SC electrons and its effects on the STM-ESR signal. We first note that at the operator level, the torque can be written asτ = −n 1 ×ĥ, where we introduce the magnetic field operatorĥ(t) = ∂ n 1 H tot (t). Then, for a given manybody state | (t) that acts in the occupation number basis, we have For a static BdG Hamiltonian, we can write H tot = being the BdG eigenvalues and γ i (γ † i ) being the Bogoliubov annihilation (creation) operator found from diagonalization. Then, in such a case we obtain the average field: where f i = γ † i γ i is the occupation of state i. This encodes both the well-known RKKY interaction mediated by the bulk states as well as the (static) YSR contribution [91]. The dynamics can induce transitions between different instantaneous energy levels and, in general, a full diagonal form for the BdG Hamiltonian might not be found.
However, in our perturbative scheme in the dynamics, when q = 2β z , we can neglect the transitions caused by β x and β y . Moreover, since 1 , the bulk states are also unaffected. Then, the many-body Hamiltonian is still diagonal and the magnetic field reads where the first and second terms determine the bulk contribution (all levels i empty or f i = 0) and YSR in-gap states contributions, respectively. Importantly, eff i are the full single-particle energies that include the shifts induced by the dynamics (which, in more formal language, corresponds to Berry-phase effects [65]). The YSR states that define the qubit states correspond in the many-body picture to the configurations f 1(2) = 0(1) and f 1(2) = 1(0). Thus, the field for each qubit state is where σ = 0, 1. Consequently, one can find the corresponding torques from τ σ = −n 1 × h σ , τ s = −n 1 × ∂ n 1 q and τ d = −2n 1 × ∂ n 1 β z , as presented in the main text (note that these are more general, as they assume arbitrary changes in the angles θ and φ).
For the numerical evaluation of the torques, we employ a Green's function approach that describes the dimer when the target spin precesses circularly. In this case, an exact solution can be found, assuming φ(t) = t and θ = const, where is the precession frequency. Indeed, the dynamical problem in Eq. (A2) can be made static by rotating it with the time-dependent unitary transformation U(t) = e −i( /2)σ z t . In this frame, the stationary Schrödinger equation from Eq. (A2) can be written as where the second term acts as a fictitious magnetic field on the superconductor. Following Ref. [60], the wave function at any point r can be written as ψ(r) = r j ∈0,R G 0 (r − r j , E)V j ψ(r j ), where 040347-12 with (C7) We then obtain the following set of eigenvalue equations: and the in-gap spectrum can be found numerically for arbitrary frequencies /2 < 1 from the 8 × 8 determinant Then, the associated torques (that include all orders in /2) can be evaluated as in the previous subsection. The plots depicted in Fig. 4 in the main text are obtained assuming the deep Shiba limit (α 1,2 ∼ 1) and /2 1 (adiabatic driving), which is the relevant regime in this work. Nevertheless, this approach can be readily employed to study the effects of the dynamics beyond the adiabatic realm.

APPENDIX D: LINEARIZATION OF LLG EQUATION AND RESONANCE-FREQUENCY RENORMALIZATION
The LLG equation describing the dynamics of the classical spin S 1 in the presence of the torque τ σ torque pertaining to the YSR qubit in state σ = 0, 1 can be written asṠ is the external magnetic field, being the sum of a constant term along z, which defines the bare resonance frequency 0 = γ B 0 , and a weak inplane rf component. Specifically, we consider B ⊥ (t) = B ⊥ [cos( 1 t), sin( 1 t), 0] and B ⊥ B 0 . In the stationary limit, the impurity spin can be written as S 1 (t) = S z z + δS(t), with δS(t) = S ⊥ [(cos( 1 t + φ) , sin ( 1 t + φ), 0], where S ⊥ = S sin θ ≈ Sθ and S z = S cos θ ≈ S and φ quantifies the lagging of the spin with respect to the driving field. In this limit, we can also expand the torque τ σ in terms of the small parameter θ , which in turn gives [65] where From the above equations, we can readily evaluate both the amplitude S ⊥ and the phase lag φ, respectively: The resonance frequency r,σ of the precessing spin is shifted depending on the qubit state σ = 0, 1 as Note that each type of torque will also contain a constant contribution, independent of the qubit state, that originates from the (occupied) bulk states. Hence, we can write τ σ s → τ bs + τ σ s and τ σ d → τ bd + τ σ d , where the index b labels the bulk contribution. Nevertheless, as shown in Ref. [91], the static bulk contribution is negligible for R ≥ 1, whilem as argued in Ref. [65], the dynamical contribution of the bulk states is negligible in the adiabatic regime. We can then extract the resonance-frequency difference as which reflects only the in-gap state effects. The Rabi frequency, and hence the time period of the Rabi oscillation, is determined by β x . Notably, the β i terms are much weaker for deviation δθ around θ = 0 as compared to such deviation around θ = π , making it inefficient for manipulation in the parallel configuration. In Fig. 6(a The probability amplitude of the qubit initialized in state |0 to remain in that state after time T a , 0|ψ(T a ) , in terms of T q = / q for δα = 0.05 or t h /δα = 0.34 (in the inset, δα = 0.01 or t h /δα = 1.75). The blue curve is the result obtained from the full numerical implementation of the evolution, while the red-dashed line is 0|ψ(T a ) evaluated in the geometrical limit using Eq. (E2). As T a is increased, the full curve starts to deviate from the geometric limit that is valid for T a T q and already almost reaches unity for T a ∼ T q , consistent with the adiabatic result. All plots are obtained using k F = 13.55.

Numerical approach for qubit-state evolution
The time-evolution operator corresponding to the qubit Hamiltonian can be written as U(t, t 0 ) = Te where T represents the time-ordering operator. We implement the evolution of the qubit state by performing time slicing with a small increment δt, so that the evolution operator during one slice can be expanded as U(t, t − δt) = 1 − (i/ )H q (t)δt. Then, starting from the initial state |0 , the qubit state at time t can then be written as which we evaluate numerically for |ψ(t) by evolving the state under the sequence of pulses described in the main text.

Analytical approach for qubit-state evolution in the geometric regime
The qubit-state evolution subjected to the pulse θ(t) = π tanh(2π t/T a ) can be studied analytically in two extreme limits: the adiabatic (T a T q ) and geometric (T a T q ) limit, respectively. In the adiabatic limit, the qubit evolution is trivial, as it remains in state |0 during the pulse. In the geometric limit, the energy splitting q becomes unimportant (and thus can be neglected) and the qubit evolution is solely determined by β y y ∝θ . Then, the evolution operator under an arbitrary rotation of the qubit from θ = 0 to a final θ 0 reads and thus the first pulse in Fig. 3(a) corresponds to θ 0 = π in the geometric limits. Figure 6(b) shows the probability amplitude of the qubit state to be in |0 as a function of time T a (scaled with T q ), c 0 (T a ), starting from |0 . The blue solid line represents c 0 evaluated numerically, while the red dashed line corresponds to the geometric limit evaluated as c 0 = cos A(π ). The deviation of the geometric amplitude from the adiabatic result increases with decreasing δα, this being due to the effect of tunneling, which makes it easier for the qubit to explore the Bloch sphere. For δα > q , the deviation is negligible and c 0 (T a ) depends weakly on the pulse length T a . Thus, this situation is preferable for the qubit manipulation.

Geometric effects around θ = π
The Rabi oscillation amplitude is also reduced because the geometric pulse θ(t) = (π − δθ) tanh(2π t/T b ) leads to a probability amplitude c 1 to excite the qubit from state to state |1 . Assuming that t h < δα and δθ 1, we obtain This linear increase of c 1 with δθ shows good agreement with the full numerical results, as shown in Fig. 7 for various values of δα. In the limit t h δα, c 1 → 0 and the amplitude of the Rabi oscillations approaches unity, as argued in the main text.

APPENDIX F: DECOHERENCE OF YSR QUBIT
Below, we give the detailed analysis of the decoherence in the YSR qubit induced due to the fluctuations in the magnetic moments and phonon coupling to the Shiba electrons.

Magnon-induced decoherence
Here, we provide details on the decoherence of the YSR qubit by the stochastic fluctuations in the magnetic moment orientations. The fluctuations of the spin k = 1, 2 can be accounted for by performing the substitution n k → n k + δn k (t), where δn k (t) = e k,1 δn k,1 + e k,2 δn k,2 ⊥ n k describe the induced fluctuations of the magnetic moment perpendicular to the deterministic orientations n k .
Here, e k,i and δn k,i , with i = 1, 2, label the orthogonal fluctuation directions and the corresponding magnitudes, respectively, with e k,1 = z × n k /|z × n k | and e k,2 = n k × (z × n k )/|n k × (z × n k )|.
Then, the coupling between the qubit and the classical spins changes accordingly, V → V + δV(t), with whereĥ k is the magnetic field operator acting on the electrons. Projecting the above Hamiltonian onto the qubit subspace leads to the following extra contribution: where χ μν k represent the components of the tensor coupling between the fluctuations and the qubit, which can be extracted from above: with the trace being taken over the qubit states. In the following, we assume that the external driving is absent and only focus on the static coherence properties. Then, we can evaluate explicitly the matrix elements of the field: σ |ĥ k |σ = −(−1) σ 1 2 ∂ n k q = (−1) σ +k 1 2 ∂ θ q e k,2 , (F4) σ |ĥ k |σ = (−1) σ q σ |∂ n kσ = −(−1) σ +k q σ |∂ θσ e k,2 + 1 sin θ σ |∂ φσ e k,1 , where the factor (−1) k reflects the fact that the magnetic fields are opposite for given relative angles. Note that for δα = 0, the diagonal terms vanish at both θ = 0 and θ = π . From the above expressions, we can write the total magnetic field operator acting in the qubit subspace aŝ h q k = (−1) k 1 2 ∂ θ q z e k,2 − q β y x e k,2 + β x sin θ y e k,1 , where β x = β x /φ and β y = β y /θ in the qubit Hamiltonian in Eq. (B6). From here, the matrix χ μν k can be readily identified. Let us evaluate the above field for the two cases of interest θ = 0 (parallel) and θ = π (antiparallel) configurations, respectively. In the former case,ĥ q k ≡ 0, meaning 040347-15 that no dephasing or relaxation occurs because of the coupling to the magnetic fluctuations, while in the latter casê h q k = (−1) k+1 g c x e k,2 + y e k,1 , where g c = t h sin(k F R + π/4) 4α 1 α 2 (F7) is the effective coupling strength of the qubit to the fluctuations, the magnitude of which is dictated by the tunneling t h . For this specific orientation, Eq. (F2) becomes δH q (t) = g c k (−1) k+1 δn k,2 (t) x + δn k,1 (t) y . (F8) We are now in a position to calculate the decoherence rates engendered by this coupling. We first introduce the noise power spectrum pertaining to the fluctuations δn k (t) in the Fourier space: where the averages are taken over the thermal equilibrium and we assume that the fluctuations of the two spins are not correlated. Within the Bloch-Redfield framework [85], the dephasing and the longitudinal relaxation rates read, respectively, φ,m = k=1,2 (|χ 1z k | 2 + |χ 2z k | 2 )S k 11 (0), 1,m = μ,ν k=1,2;σ =± χ μσ k χ νσ k S k μν (σ q ), where χ μ± k = χ μx k ± iχ μy k . The pure dephasing rate φ vanishes at both θ = 0 and π and 1,m = 0 at θ = 0. The relaxation rate at θ = π is 1,m = 1/T 1,m = 2g 2 c k=1,2 while the dephasing time satisfies T 2,m = 2T 1,m . In order to provide estimates, we need to describe the noise spectrum of the magnetic fluctuations. To do that, we start by employing the stochastic LLG equation describing the magnets in the presence of magnetic noises (here, we disregard the effect of the qubit on the dynamics, as it would only manifest in higher orders in the coupling): n k = −γ n k × B k,eff + δB k (t) + α g n k ×ṅ k . (F13) Here, B k,eff = −γ −1 δF S (n k )/δn k is the effective magnetic field acting on the impurity, with F S (n k ) being the k th classical spin free energy, and δB k (t) is the stochastic magnetic field the Fourier components δB k,μ (ω) of which, with μ = e k,1 , e k,2 , satisfy the fluctuation-dissipation relation [87]: In order to describe the experimental observations [88,90], we assume z to be an easy axis (the spin orients perpendicular to the surface), so that the free energy can be written as in the main text: where B = B z z is external magnetic field along z and κ is the strength of the anisotropy, which is assumed to be identical for the two spins. Consequently, the effective magnetic field that determines the dynamics can then be written as B k,eff = B k,eff z, with the magnitude B k,eff = B z + (κ/γ )n k,z . For the antiparallel alignment and considering the deterministic direction of the spins to be along the z axis, inserting the effective field in the LLG equation, we can extract the noise spectrum, for which we find where k = 1, 2. To provide some estimates, we assume that α 1 = 1.15, α 2 = 1.1, R = 2.9, and α g = 0.001. Considering the magnetization anisotropy energy κ = 0.1 meV and 0 = 25 GHz, we find T 2,m ≈ 7 μs, allowing around 800 Rabi oscillations to be experimentally observable before the qubit is hampered by the decoherence stemming from magnetic fluctuations. In the presence of a small applied magnetic field-say, B z = 0.2 T, which corresponds to one. This is in stark contrast to the relaxation induced by the impurities fluctuations, 1,m ∝ sin 2 (θ/2), which vanish in the parallel configuration. Consequently, the two mechanisms do not compete with each other in the two qubit-operation configurations, allowing us to extract their effects separately. The reason for such behavior is that phonons cannot cause spin-flip transitions during the tunneling processes. In the antiparallel configuration, the tunneling of the YSR quasiparticles involves spin flips and thus results in zero coupling. The pure dephasing rate φ,ph = 0, as a consequence of the phonon power spectrum J ph (ω) ∝ ω 2 in the current 2D setup [85]. Then, similarly to magnons, the phonon-induced dephasing entirely originates from longitudinal relaxation, or T 2,ph = 2T 1,ph .
In order to provide estimates, let us focus on a 2D Pb SC slab. We assume that x = 10 nm, i.e., much smaller than the coherence length ξ , c s ≈ 10 4 m/s, and g ph ≈ 1.2 × 10 −8 μeV cm 3 . This leads to 1,ph (θ ) = 1.72 × 10 5 cos 2 (θ/2) s −1 , reaching its maximum at θ = 0. The corresponding phonon-induced relaxation time in the parallel configuration is then T 1,ph ≈ 5.8 μs, which is comparable in magnitude to that stemming from impurity fluctuations.