Electron-boson-interaction induced particle-hole symmetry breaking of conductance into subgap states in superconductors

Particle-hole symmetry (PHS) of conductance into subgap states in superconductors is a fundamental consequence of a noninteracting mean-field theory of superconductivity. The breaking of this PHS has been attributed to a noninteracting mechanism, i.e., quasiparticle poisoning (QP), a process detrimental to the coherence of superconductor-based qubits.Here, we show that the ubiquitous electron-boson interactions in superconductors can also break the PHS of subgap conductances. We study the effect of such couplings on the PHS of subgap conductances in superconductors using both the rate equation and Keldysh formalism, which have different regimes of validity. In both regimes, we found that such couplings give rise to a particle-hole $asymmetry$ in subgap conductances which increases with increasing coupling strength, increasing subgap-state particle-hole content imbalance and decreasing temperature. Our proposed mechanism is general and applies even for experiments where the subgap-conductance PHS breaking cannot be attributed to QP.


I. INTRODUCTION
Subgap states in superconductors are key features of topological superconducting phases [1-13] which offer great promise for quantum information processing [14,15]. Tunneling transport into such Andreev bound states (ABSs) provides the most direct and commonly employed method to detect them [16][17][18][19] (Hereafter ABS refers to any subgap state in superconductors.) Most of our understanding of tunneling into superconductors is based on the celebrated Blonder-Tinkham-Klapwijk (BTK) formalism [20]. One universal consequence of this theory is a precise particle-hole symmetry (PHS) of the conductance into any ABS in a superconductor [18,21,22]. Specifically, this theory predicts that the differential conductance at a positive voltage V inside the superconducting gap precisely match its counterpart value at −V . This symmetry has been shown to be a consequence of the PHS of the mean-field Hamiltonian used in the BTK formalism. However, numerous experiments over two decades [3-8, 11, 23-34] have often observed particle-hole (PH) asymmetric subgap conductances. One way to reconcile this PH asymmetry with the BTK theory is to introduce quasiparticle poisoning induced either by coupling the ABS to a fermionic bath [22,35,36] or through a relaxation process from the ABS to the superconductor's quasiparticle continuum [37].
Quasiparticle poisoning (QP) [38][39][40] refers to a process where an electron tunnels from the bulk of the superconductor to an ABS which changes the occupation (parity) of the ABS. Since the parity is used as the qubit state, QP then introduces bit-flip errors [41][42][43]. Moreover, as QP breaks the PHS of subgap conductances [22,35,36], * setiawan@uchicago.edu one may be tempted to associate the PH asymmetry to short qubit lifetime. We will show that this correlation is not true in general as contrary to commonly held beliefs, the PH asymmetry can also arise without QP.
In this paper, we propose a generic mechanism for PHS breaking of subgap conductances without changing the superconductor's parity state, namely, the coupling between ABSs and bosonic modes. While quantum tunneling in dissipative systems has been widely studied [44,45], previous works consider coupling between bosonic baths and superconductors without ABSs. Motivated by tunneling experiments into ABSs [3-8, 11, 23-34], here we study tunneling transport from a normal lead into an ABS coupled to bosonic modes, e.g., phonons [46,47], plasmons [48], or electromagnetic fields [49], in the superconductor. Our system has a local fermion parity analogous to the spin-boson model [50] with a caveat that our ABSs can participate in transport. Crucially, our study of transport into an ABS coupled to bosonic modes and its relation to PHS breaking of subgap conductances has not been undertaken before. To this end, we present ways to enforce fermion-parity conservation in treating interaction effects on transport into ABSs. We consider two different limits: weak and strong tunneling regimes where the ABS-lead tunnel strength is smaller and larger than the thermal broadening ∼ k B T , respectively. The weak tunneling limit is studied using the rate equation [51,52], which is valid for all values of ABS-boson coupling strength where the tunneling rates are calculated using Fermi's Golden Rule (FGR). In the strong tunneling limit, we study the transport using the Keldysh formalism and treat the ABS-boson coupling within the mean-field approximation. We begin by using FGR to show that while subgap conductances in gapped superconductors (superconductors without baths) preserve PHS even with interactions (including strongly correlated superconductors), the PHS is broken for superconductors with gapless excitations (e.g., phonons, quasiparticles, etc.). The simplest application of FGR [23,53,54] considers the conductance into an ABS at positive [ Fig. 1(a)] and negative subgap energies [ Fig. 1(b)] to arise from the tunneling of electrons and holes, respectively, into the ABS (changing the ABS occupancy n from 0 → 1). The tunneling rates of electrons [R 0→1;e in Fig. 1(a)] and holes [R 0→1;h in Fig. 1(b)] can be calculated from FGR to be proportional to the particle and hole component of the ABS wavefunction, respectively. This suggests that the tunneling conductance into an ABS with different weights of particle and hole component is PH asymmetric. However, this simple argument implicitly assumes the presence of QP [22], which empties out the electron from the ABS after each tunneling event such that its occupancy n returns to n = 0. This implicit assumption can be avoided by taking into account the change in the ABS occupancy n = 0, 1 after each tunneling. (1)]. The first tunneling flips the ABS occupancy n from 0 → 1 and occurs with rates (a) R 0→1;e or (b) R 0→1;h . The second tunneling, which flips n from 1 → 0, occurs with rates (a) R 1→0;e or (b) R 1→0;h . Without bosonic baths, R 0→1;e = R 1→0;h and R 1→0;e = R 0→1;h giving a PH symmetric conductance. However, in the presence of bosonic baths, the second tunneling occurs with a higher rate since it can transfer lead electrons within a larger energy range near the ABS energy (shaded blue region) where the energy difference can be dumped by emitting bosons (green squiggly lines). Therefore, R 0→1;e = R 1→0;h and R 1→0;e = R 0→1;h resulting in a PH asymmetric conductance.
As seen in Fig. 1(a), the electron tunneling flips n either from 0 → 1 (with a rate R 0→1;e ) or vice-versa (with a rate R 1→0;e ). Since each tunneling event flips n →n ≡ 1 − n, a full cycle of transferring a pair of electrons returns the occupancy to the initial |n = 0 occupancy state. The total time for this process that transfers a charge of 2e is τ = (R 0→1;e ) −1 + (R 1→0;e ) −1 leading to a current I = 2e/τ . Combining this result with the analogous argument for negative voltages [ Fig. 1(b)] leads to the expression for the tunneling current (we give a more detailed derivation later): R 0→1;e R 1→0;e R 0→1;e + R 1→0;e for eV |ε A | + k B T , whereε A is the interaction-renormalized ABS energy, k B is the Boltzmann constant, and T is the temperature. The constraints on the voltage V in Eq.
(1) are needed to separate the electron and hole tunneling shown in Fig. 1.
Using FGR, we calculate the electron and hole tunneling rates as R n→n;e ∝ | n|d † A |n | 2 and R n→n;h ∝ | n|d A |n | 2 , whered † A andd A are the electron and hole creation operators in the ABS, respectively. Since R 0→1;e = R 1→0;h and R 1→0;e = R 0→1;h , the current [Eq. (1)] is antisymmetric I(V ) = −I(−V ) and the corresponding subgap conductance shows PHS for a gapped superconductor even with interactions (including strongly correlated superconductors). However, as shown below, this PHS is broken in the presence of bosonic baths.
Here, ε A is the ABS energy,γ (γ † ) is the Bogoliubov annihilation (creation) operator of the ABS, λ is the ABSboson coupling strength,b (b † ) is the boson annihilation (creation) operator, and Ω is the boson frequency. The operatorĉ L,k (ĉ † L,k ) annihilates (creates) the lead electron with momentum k and energy ε L,k . The electron tunneling, represented by the HamiltonianĤ T [37,54], occurs with a strength t and involves the electron operator of the lead [ĉ † L = dkĉ † L,k /(2π)] and ABS (d A = uγ + vγ † [55]) where u ≡ u(x = 0) and v ≡ v(x = 0) are the particle and hole component of the ABS wavefunction at the junction (x = 0). We renormalize the ABS wavefunction such that |u| 2 +|v| 2 = 1. The ABS-boson coupling term can be derived from the microscopic electron-boson interaction by projecting it onto the lowest-energy (ABS) sector (see Sec. I of Ref. [56]). This term can be eliminated using the Lang-Firsov canonical transformationĤ = eŜĤe −Ŝ , whereŜ = λ [56]). The op-eratorŶ is analogous to the operator e −iφ in Ref. [45], through the identification λ(b † −b)/Ω = iφ whereφ is the phase operator of the electromagnetic field used in Ref. [45]. Therefore, our results apply generally to all bosonic modes including electromagnetic fields and plasmons.
The current operator isÎ whereṄ L is the time derivative of the lead electron number. The current is proportional to the tunnel coupling strength Γ ≡ 2πt 2 ν 0 where ν 0 is the density of states at the lead Fermi energy. The ratio Γ/k B T determines two different transport regimes: weak (Γ/k B T < 1) and strong (Γ/k B T > 1) tunneling regimes.

A. Rate equation
We first study the weak tunneling limit using the rate equation [51,52], which applies for all values of λ. Without the lead coupling, the eigenstates of the ABS-boson system are |n, q with eigenenergies E n,q = nε A + qΩ, where the indices n = 0, 1 and q ∈ Z ≥0 denote the ABS and boson occupation numbers, respectively. The tunneling of electrons and holes from the lead to the ABS introduces transitions between the eigenstates |n, q . If the boson relaxation rate is faster than the tunneling rate Γ/ (typically true in experiments [58]) such that the bosons acquire the equilibrium distribution P b q = e −qΩ/kBT (1 − e −Ω/kBT ), the probability that the system in the state |n, q can be factorized as P n q = P n P b q . In the steady state, P n satisfies the rate equation ( where n|d † A |n and n|d A |n are the bare tunneling matrix elements for electrons and holes, respectively, [56]). Solving Eq. (3) together with the normalization condition P 0 + P 1 = 1, we obtain P 0 and P 1 . Substituting these probabilities into the current I = e n P n R n→n;e − R n→n;h [51], we have We can show that Eq. (5) reduces to Eq. (1) by noting that the hole tunneling is energetically forbidden at large positive voltages (R 0→1;h , R 1→0;h ≈ 0 for eV |ε A | + k B T ) and so is the electron tunneling at large negative voltages [R 0→1;e , R 1→0;e ≈ 0 for eV −(|ε A | + k B T )]. While Eq. (1) implies PHS for subgap conductances of gapped superconductors, the inclusion of a bosonic bath modifies the tunneling rates in Eq. (1) so as to break the conductance PHS. This PHS breaking can be understood more intuitively in the lowtemperature limit as follows. The first tunneling, occurring with rates R 0→1;e [ Fig. 1(a)] or R 0→1;h [ Fig. 1(b)], transfers only lead electrons or holes near the lead Fermi energy and is accompanied by emission of small number of bosons since there are only a few occupied electrons (holes) above (below) the Fermi level. In contrast, the second tunneling, whose rates are R 1→0;e [ Fig. 1(a)] or R 1→0;h [ Fig. 1(b)], has a higher probability of boson emission since it transfers electrons and holes with energies deep inside the lead Fermi energy. This means that R 0→1;e = R 1→0;h and R 1→0;e = R 0→1;h for tunneling into ABSs in superconductors with gapless excitations (e.g., phonons) unlike the gapped superconductor case. Therefore, I(V 0 ) = −I(−V 0 ) [Eq. (1)] and the conductance becomes PH asymmetric, i.e., dI dV V =V0 = dI dV V =−V0 (see Sec. IV A. of Ref.

B. Keldysh
For strong-tunneling limit (Γ > k B T ), we compute the current using the mean-field Keldysh formalism. We begin by rewriting Eq. (2a) in terms of the boson displace- We calculate the mean-field energy ε A + √ 2λ x by selfconsistently solving for x where · · · is the expectation value with respect to the mean-field eigenfunction. To this end, we solve for ∂ Ĥ A /∂ x = 0 and ∂ Ĥ A /∂ p = 0, giving x = − √ 2λ Ω γ †γ and p = 0. The ABS Green's function in the Lehmann represen- x ) where Φ + = (1, 0) T and Φ − = (0, 1) T are the Nambu spinors written in the Nambu basis (γ,γ † ) T . Following Ref. [37], we use the Green's function to evaluate the current as (see Sec. VII of Ref. [56]) where with Γ u = Γ|u| 2 and Γ v = Γ|v| 2 . The mean-field boson displacement x in Eq. (9) is evaluated self-consistently as where (G < αβ ) A = i Ψ † αA Ψ βA is the ABS lesser Green's function (see Sec. VIII of Ref.
[56]) with Ψ A = (γ,γ † ) T and σ z being the z-Pauli matrix in the Nambu space.  shows that the peak conductance increases with increasing ABS-boson coupling strength λ contrary to the rateequation results. However, similar to the rate-equation, the conductance PH asymmetry ζ increases with increasing λ [61]. Unlike the rate equation, the Keldysh approach shows that in the strong-tunneling regime the PHS breaking holds also for high-frequency bosons (see Sec. IX of Ref. [56]), since it arises from nonperturbative effects of tunneling, i.e., the PH asymmetry of the mean-field boson displacement value x .
Our model of tunneling into boson-coupled ABS [Eq.
(2)] can explain the origin of PH asymmetry for subgap conductance observed in a hard superconducting gap [5, 6, 31, 32, 34] which cannot be accounted for by QP. However, similar to QP this model also results in conductance peak areas which are independent of temperatures (see Sec. IV B. of Ref. [56]). In Sec. IV below, we consider another related model, i.e., a boson-assisted tunneling model. This model can not only give rise to PHS breaking of subgap conductances but also account for experimentally observed conductance features which cannot be attributed to QP, e.g., an increase in the conductance peak area with temperature [31].

IV. MODEL II. BOSON-ASSISTED TUNNELING INTO ABS
In this section, we consider boson-assisted tunneling into an ABS via virtual hopping of electrons or holes from the lead into higher-lying states in superconductors which are boson-coupled to the ABS. The higher-lying states can be either higher-energy ABSs or states from the continuum above the gap. By integrating out the higherlying states, we derive the effective low-energy Hamiltonian for the boson-assisted tunneling into the ABS as (see Sec. X of Ref. [56]) Note the extra (b+b † ) term in the above tunneling Hamiltonian as compared to Eq. (2c) in Sec. III. Figure 4 shows the current and conductance calculated using the rate equation within the boson-assisted tunneling model for different temperatures T and ABS-boson coupling strengths λ. Contrary to the tunneling model in Sec. III where the current at V = ±∞ is independent of temperature (see Sec. IV B. of Ref.
[56]), for the boson-assisted tunneling model the current magnitude at V = ±∞ increases with increasing temperature. This is because the boson-assisted tunneling rate (see Sec. X of Ref. [56]) is proportional to q|(b+b † ) 2 |q , which increases with increasing temperature. Crucially, we find that the current I(V = ±∞) or equivalently the peak area of the conductance versus voltage curve has a faster-than-linear increase with temperature [inset of Fig. 4(a)], providing excellent agreement with experimental results [31]. Since QP preserves the conductance peak area under different temperatures and necessarily induces "soft-gap" conductance features, our proposed boson-assisted tunneling process is thus more likely to be responsible for the PH-asymmetric subgap conductances inside a hard superconducting gap observed in Ref. [31]. Contrary to the model in Sec. III, for the bosonassisted tunneling model, the current calculated at large positive and negative voltages need not be perfectly antisymmetric, i.e., increases with increasing ABS-boson coupling strength λ [see inset of Fig. 4(b)]. This current PH asymmetry (or equivalently the asymmetry between the conductance peak area for positive and negative voltages) as well as the dependence of the conductance peak area with temperature can serve as signatures for the boson-assisted tunneling process. Similar to the model in Sec. III, the conductance PH asymmetry ζ calculated using the boson-assisted tunneling model also decreases with in-

V. CONCLUSIONS
Contrary to widely held belief, we show that the PHS breaking of subgap conductances in superconductors can arise without QP. Specifically, the coupling of ABSs to a bosonic bath (or multimode bosonic baths [62]) can break the PHS of subgap conductances without changing the superconductor's parity state. Therefore, contrary to QP, our mechanism is not detrimental to the coherence of superconductor-based qubits. (Topological qubits are exponentially protected from the bosonic bath dephasing due to the spatial separation of Majoranas [63].) We find that the conductance PH asymmetry increases with increasing ABS's PH content imbalance, increasing ABS-boson coupling strength and decreasing temperature. Our theory is general as it applies to all ABSs, e.g., quasi-Majorana states [64, 65], Yu-Shiba-Rusinov states [66-68], Caroli-de Gennes-Matricon states [69], etc., which couple to bosonic modes such as phonons, plasmons, electromagnetic fields, etc., in superconductors. Contrary to QP, our mechanism applies even for ABSs observed inside a hard superconducting gap [5, 6, 31, 34] and can give rise to an increase in the conductance peak area with temperature as observed in experiments [31].
Our PHS breaking mechanism results from boson emissions or absorptions accompanying the electron/hole tunneling. Since these bosons such as phonons are ubiquitous in superconductors, we expect electronphonon interactions (EPIs) to significantly affect transport in superconductors, particularly the semiconductorsuperconductor heterostructures used to realize topological superconductors [3-5, 7, 8, 11, 28-30]. In fact, measurements of transport in semiconductors have observed features [70-74] associated with EPI that are theoretically understood [75][76][77]. We estimate that for a typical topological superconductor which uses either an InAs or InSb semiconductor with a length of ∼ 1 µm (having a phonon frequency Ω ∼ v s π/ = 7.2 µeV where v s ≈ 3.5 × 10 5 cm/s [78-80] is the sound velocity), EPI can give rise to a conductance PH asymmetry in the tunneling limit for ABSs with energies ε A Ω/2 = 3.6 µeV. Therefore, contrary to QP, EPI does not affect the zerobias Majorana conductance.
Compared to diagrammatic techniques, FGR is a more controlled approach in treating the effect of interactions on transport in superconductors (even for strongly correlated superconductors) for the strict tunneling limit. This is because interaction diagrams can generate an imaginary self-energy [37], resulting in a conductance PH asymmetry similar to QP [22]. Therefore, it is crucial to enforce a fermion parity conservation in the diagrammatic treatment of ABS-boson couplings like our meanfield treatment of interactions in the Keldysh formulation. Our work thus motivates the formulation of the nonequilibrium Green's function beyond the mean field approximation that conserves fermion parity. We note that our mechanism is quite distinct from the subgapconductance PHS breaking due to the bias-voltage dependence of the tunnel barrier [81]. While this mechanism can be treated within the Keldysh approach by moving the interaction term from the ABS to the barrier, it vanishes in the tunneling limit where FGR applies.   [61] Since we ignore the Fock term in the mean-field approximation, the conductance calculated from the Keldysh approach has no boson sidebands.
[62] While here we focus exclusively on single-mode bosonic baths, our proposed mechanism is expected to hold also for multimode bosonic baths since the result for multimode bosonic baths is qualitatively similar to averaging multiple results for different single-bosonic modes. This averaging is justified as bosonic modes do not interact with each other and are therefore independent.  Supplemental Material for "Electron-boson-interaction induced particle-hole symmetry breaking of conductance into subgap states in superconductors"

I. DERIVATION OF THE ABS-BOSON COUPLING FROM THE MICROSCOPIC ELECTRON-BOSON INTERACTION
In this section, we derive the ABS-boson coupling term in Eq. (2a) of the main text from the microscopic electronboson interaction. We begin by writing a generic Hamiltonian for a superconductor with an ABS aŝ where h l,m describes the dynamics of the electrons in the superconductor with an ABS, ∆ l,m is the superconducting pairing potential, d † l,m (d l,m ) is the electron creation (annihilation) operator of the superconductor and the indices l, m represent both the orbital and spin degrees of freedom. We can diagonalize the above Hamiltonian using the Bogoliubov transformationγ where the lowest energy level corresponds to the ABS energy, i.e., ε 1 = ε A . The Hamiltonian of the electron-boson coupling is given bŷ where g lm is the electron-boson coupling strength, andb (b † ) is the boson annihilation (creation) operator. Substituting into Eq. (S-4), we havê where we have definedλ Projecting the above Hamiltonian into the lowest energy sector α = β = 1 which corresponds to the ABS energy sector, we havê where we have definedγ ≡γ 1 as the Bogoliubov operator for the ABS, λ ≡ 2λ (c) 11 as the ABS-boson coupling strength, and χ ≡ lm dx dx g lm (x, x ). Note that in evaluating Eq. (S-8), we have used the anticommutation relation {γ, γ † } = 1, {γ, γ} = 0, and {γ † , γ † } = 0. We can eliminate the term χ(b † +b) in Eq. (S-8) by introducing the shift b →b − χ/Ω andb † →b † − χ/Ω which gives the ABS Hamiltonian aŝ Introducing the shift ε A → ε A + 2λχ/Ω and shifting the overall energy by χ 2 Ω , i.e.,Ĥ A →Ĥ A − χ 2 Ω , we have the Hamiltonian for the boson-coupled ABS as in Eq. (2a) of the main text:

II. LANG-FIRSOV TRANSFORMATION
In this section, we follow Ref.
[53] to derive the matrix elements for the tunneling of electrons (d † A ) and holes (d A ) and the boson absorption or emission matrix elements Y qq in Eq. (4) of the main text. We begin by writing the Hamiltonian of an ABS coupled to a one-dimensional normal lead and bosonic modes, e.g., phonons, plasmons, etc., as the sum of the Hamiltonian of a boson-coupled ABS, lead and tunnel coupling,Ĥ =Ĥ A +Ĥ L +Ĥ T [Eq.
(2) of the main text], whereĤ Here, ε A is the ABS energy,γ (γ † ) is the Bogoliubov annihilation (creation) operator of the ABS, λ is the ABS-boson coupling strength,b (b † ) is the boson annihilation (creation) operator and Ω is the boson frequency. The operator c L,k (ĉ † L,k ) annihilates (creates) the lead electron with momentum k and energy ε L,k . The tunneling Hamiltonian H T [37,54] represents the electron tunneling between the normal lead and ABS, where the electron annihilation operator of the lead and ABS at the junction given byĉ L = dkĉ L,k /(2π) andd A , respectively. The operatord A is obtained by projecting the operatord 1 (x = 0) [Eq. (S-5a)] to the ABS energy sector (α = 1), where we havê d A = uγ + vγ † . For notational simplicity, here we define γ ≡ γ 1 , u ≡ u 11 (x = 0) and v ≡ v 11 (x = 0) where u and v are the particle and hole components of the ABS wave function at the junction (x = 0). In this paper, we renormalize the ABS wave function such that |u| 2 + |v| 2 = 1. Note that since we consider only the subgap state and ignore the above-gap states, the relationd A = uγ + vγ † is only approximate which makesd A nonfermionic. The operatord A becomes fermionic if all the states in the superconductor including the above-gap states are taken into account [see Eq. (S-5a)]. Our conclusion on the PHS breaking of the subgap conductance due to the ABS-boson coupling does not rely on the fermionic properties ofd A .
To eliminate the ABS-boson coupling, we can transform the Hamiltonian [Eq. (S-11)] using a canonical transfor-mationĤ we can write the transformed annihilation and creation operators for the Bogoliubov quasiparticles, electrons and bosonic modes asγ =γŶ , (S-15a) . Under this transformation, the number operator remains the same, i.e.,γ †γ = γ †γŶ †Ŷ =γ †γ and the Hamiltonians [Eq. (S-11a) and Eq. (S-11c)] transform aŝ We can evaluate the matrix elements for the electron and hole tunneling which change the ABS occupancy number n from 0 → 1 and the boson occupancy from q → q as respectively. Using the Baker-Campbell-Haussdorf formula, we havê where |Y qq | 2 is symmetric under the interchange q ↔ q . Note that in going to the second line of Eq. (S-19), we have used the following relations:

III. RATE EQUATION AND TUNNELING RATES
The stationary-state rate equation satisfied by the probability P n q of an ABS-boson system being in the state |n, q , i.e., having an ABS occupation number n and boson occupation number q, is given by [51, 52] 0 = ∂P n q ∂t = q Pn q Rn →n;e q →q + Rn →n;h q →q − P n q q R n→n;e q→q + R n→n;h q→q + P n q+1 η q+1;− + P n q−1 η q−1;+ − P n q (η q;+ + η q;− ) . (S-21) The second line in Eq. (S-21) represents the probability flux due to hopping of an electron (e) or hole (h) from the lead to the ABS which changes the ABS occupation number fromn ≡ 1 − n to n and the boson occupancy from q to q and vice versa. The quantity P n q denotes the probability that the system is in the state |n, q and R n→n q→q denotes the transition rate from the state |n, q to the state |n, q . The third line of Eq. (S-21) represents the boson relaxation where the boson emission and absorption probabilities are η q;+ = A(q + 1) and η q;− = Bq, respectively, with A = Be −Ω/kBT . These probability rates are consistent with the fluctuation-dissipation theorem. If the boson relaxation rate is faster than the tunneling rate Γ/ such that the bosons acquire the equilibrium distribution P b q = e −qΩ/kBT (1 − e −Ω/kBT ), the probability P n q can be factorized as P n q = P n P b q . Summing Eq. (S-21) over q for these factorized probabilities gives 0 = Pn(Rn →n;e + Rn →n;h ) − P n R n→n;e + R n→n;h + P n q P b q+1 η q+1;− + P b q−1 η q−1;+ − P b q (η q;+ + η q;− ) , For the tunneling Hamiltonian in Eq. (S-16b), the rates of the electron and hole tunneling processes can be calculated from Fermi's Golden Rule to be where n|d † A |n and n|d A |n are the bare tunneling matrix elements for electrons and holes, respectively, Y qq = q |e −λ(b † −b)/Ω |q is the boson emission or absorption matrix element, and f (E) = [1 + exp(E/k B T )] −1 is the lead Fermi function.

IV. DETAILS ON MODEL I. TUNNELING INTO BOSON-COUPLED ABS
A. Proof for the particle-hole asymmetry of boson-coupled-ABS conductance In this section, we will prove that, unless |u| = |v|, the current into a boson-coupled ABS [Eq. (5) of the main text] is in general not PH antisymmetric, i.e., I(V 0 ) = I(−V 0 ) resulting in a PH asymmetric conductance, i.e., We will prove below that the function F (x) in the denominator of Eq. (S-25) is an increasing function of x and hence the denominator in Eq. (S-25) is asymmetric with respect to the interchange V ↔ −V unless |u| = |v|. By rewriting W (x) in Eq. (S-26a) as In the following, we will prove that F (x) is a monotonic function of x. We first begin by noting that Q(ω)−Q(−ω) ≤ 0 for ω ≥ 0. The proof is as follows where in the third line we interchange q with q for the second sum and use |Y qq | 2 = |Y q q | 2 . For ω ≥ 0, the delta function forces q ≥ q implying that (P To prove that the boson-assisted tunneling model (model II) can also break the PHS of subgap conductances, we simply replace |Y qq | by |X qq − λY qq /Ω| in the above derivation, where X qq ≡ e − λ 2 2Ω 2 q |e λ Ωb † (b † +b)e − λ Ωb |q [Eq. (S-82b)]. Even though the conductance is not PH symmetric, under a simultaneous interchange of V ↔ −V and |u| ↔ |v|, the current is antisymmetric (I → −I) resulting in a symmetric conductance. This means that the conductance for |v| 2 > |u| 2 can be obtained from the conductance for |v| 2 < |u| 2 (shown in Figs. 2 and 3 of the main text) by interchanging both |u| ↔ |v| and V ↔ −V simultaneously. As a result, the higher and lower peaks switch sides when |v| ↔ |u|, which changes the sign of the PH asymmetry ζ.
The PH asymmetry of the conductance can be understood more intuitively in the limit of large positive and negative voltages |eV | |ε A | + k B T . In the large-positive-voltage regime (eV |ε A | + k B T ), hole tunneling processes are energetically forbidden [R 0→1;h , R 1→0;h ≈ 0 since W (ε A,+ ), W (−ε A,− ) ≈ 0]. On the other hand, in the large-negativevoltage regime where eV −(|ε A |+k B T ), electron tunneling processes are not energetically allowed [R 0→1;e , R 1→0;e ≈ 0 since W (ε A,− ), W (−ε A,+ ) ≈ 0]. In this limit, the current [Eq. (S-25)] thus reduces to Eq. (1) of the main text: implying that the current at large positive and negative voltages are due to sequential tunnelings of electrons and holes, respectively (see Fig. 1 of the main text). For both the boson-coupled ABS model and boson-assisted tunneling into ABS model, the current is in general not PH antisymmetric, i.e., I(−V 0 ) = I(V 0 ) or the conductance is PH asymmetric ( dI dV V =V0 = dI dV V =−V0 ) because of the rate asymmetry between the first and second tunneling processes of electrons and holes (i.e., R 0→1;e = R 1→0;h and R 1→0;e = R 0→1;h ). This rate asymmetry arises because the second tunneling process which happens at energy deep inside the Fermi level is energetically allowed to emit more bosons hence occurs with a larger rate than the first tunneling process. Without the ABS-boson coupling (λ = 0),

V. CURRENT CALCULATED FROM THE RATE EQUATION AND KELDYSH APPROACH
In this section, we show the current calculated from the rate equation (Fig. S1) and mean-field Keldysh approach (Fig. S2) corresponding to the conductance shown in Figs. 2 and 3 of the main text, respectively. As shown in Figs. S1(a) and S2(a), the current decreases with increasing ABS's PH content imbalance ||u| 2 − |v| 2 | where I = 0 when ||u| 2 −|v| 2 | = 1. This is due to the fact that the terms R 0→1;e R 1→0;e and R 0→1;h R 1→0;h in the current expression [Eq. (5) of the main text] are ∝ |uv| 2 = [1 − (|u| 2 − |v| 2 ) 2 ]/4.  Figs. S3(g,h)]. Figure S3(b) shows that the magnitude of the conductance PH asymmetry ζ has a nonmonotonic dependence on the ABS-boson coupling strength λ. In the limit ε A − λ 2 /Ω k B T (where the two conductance peaks are well separated), the conductance PH asymmetry ζ increases with increasing λ; this corresponds to the results shown in Fig. 2(b) of the main text. As λ keeps increasing, the two conductance peaks approach each other and in the regime where ε A − λ 2 /Ω < k B T , the two peaks start to overlap with each other and ζ decreases with increasing λ. Note that in the regime where ε A − λ 2 /Ω > 0, the higher peak is at positive voltage for the case where |u| 2 > |v| 2 while for the case where |v| 2 > |u| 2 , the higher peak is at negative voltage. When ε A − λ 2 /Ω = 0, the two conductance peaks merge at the zero voltage which gives a zero conductance PH asymmetry (ζ = 0). Increasing λ beyond this point splits the peaks but with the low and high peaks now switching sides which in turn changes the sign of ζ. As λ increases further, the two peaks move away from each other and the PH asymmetry ζ increases in magnitude; beyond a certain value of λ, ζ becomes weakly dependent on λ as shown in the inset of Fig. S3(b). Note that for large enough λ, the position of the conductance peaks are no longer PH symmetric [see green curve in Fig. S3(b)]. shows that the subgap conductances exhibit PH asymmetry only forε A /Ω 0.5 whereε A ≡ |εA − λ 2 /Ω| + kBT /2. The parameters used for all panels are: |u| 2 /|v| 2 = 1/9. Figure S3(d) shows that the conductance PH asymmetry ζ decreases with increasing temperature T . This is due to the fact that temperature broadens the conductance peaks. The dependence of the ABS conductance on the boson frequency Ω is shown in Fig. S3(f). The PH asymmetry ζ has a nonmonotonic behavior with the boson frequency Ω where it first increases with increasing Ω and then after reaching its maximum, it decreases with increasing Ω. The initial increase of ζ with increasing Ω can be attributed to the fact that the two conductance peaks move away from each other as Ω increases (ε A = ε A − λ 2 /Ω increases with increasing Ω). The decrease of ζ for large Ω is due to the fact that the effective ABS-boson coupling strength λ/Ω decreases with increasing Ω. Figure S3(h) shows the dependence of the conductance on the ABS energy. As shown in the inset of panel (h), the peak conductance only exhibits the PH asymmetry for Ω 2|ε A | + k B T . This can be understood from the fact that the second tunneling process [whose rate is R 1→0;e in Fig. 1(a) or R 1→0;h in Fig. 1(b) of the main text] can transfer lead electrons or holes with an energy difference up to ∼ 2|ε A | + k B T from the subgap state, where this energy difference is transferred in form of boson energy Ω. Even though in this paper, we focus only on the regime k B T λ where the boson sidebands vanish due to the thermal broadening [51], the dependence of the ABS conductance peak on the above parameters also hold true in the case where there are boson sidebands. Moreover, the PH asymmetry of the boson sidebands also have similar dependences on the above parameters as that of the ABS conductance peak.

VII. DERIVATION OF THE CURRENT IN THE KELDYSH FORMALISM
In this section, we derive the current [Eq. (8) of the main text] following Refs. [37,82,83]. We begin by writing the Hamiltonian asĤ In Eq. (S-35a),Ĥ A is the mean-field Hamiltonian of the ABS-boson system obtained by replacingx in Eq. to the HamiltonianĤ in Eq. (S-34), whereN L =ĉ † Lĉ L and N S =d † Ad A are the lead and substrate electron number, respectively, withd A = uγ + vγ † . With this transformation, the single-particle energies in the lead and substrate are measured from the chemical potential of the lead (µ L ) and substrate (µ S ), respectively, where the transformed Hamiltonian isĤ with the tunneling Hamiltonian transformed aŝ The current operator is given bŷ By taking the expectation value of the current operator, we have where σ z is the z-Pauli Matrix in the Nambu basis. In Eq. (S-40), we have introduced the hopping matrix -41) and the lesser Green's function in the Nambu space [(G < αβ ) ij = i Ψ † βj Ψ αi ] with i, j = L, A denoting the quantities for the lead and ABS, respectively, where Ψ L = (ĉ L ,ĉ † L ) T and Ψ A = (γ,γ † ) T . We can Fourier-expand the current and Green's functions in terms of the frequency ω 0 = eV / , where we have Let us denote G mn (ω) ≡ G(ω + mω 0 , ω + nω 0 ) for which we have G mn (ω) = G m−n,0 (ω + nω 0 ). The dc current which is the zeroth order (I 0 ) in the Fourier expansion of the current [Eq. (S-42a)] is given by where the superscripts ee, eh, he, and hh denote the matrix elements in the Nambu space. Using the Langreth rule [85] where g L = diag g ee L , g hh L , we have G <,ee AL,10 = t G r,ee A,11 u * + G r,eh A,11 v * e −iω0τ g <,ee L,00 + G <,ee A,11 u * + G <,eh A,11 v * e −iω0τ g a,ee L,00 , (S-45a) G <,he AL,10 = t G r,he A,11 u * + G r,hh A,11 v * e −iω0τ g <,ee L,00 + G <,he A,11 u * + G <,hh A,11 v * e −iω0τ g a,ee L,00 , (S-45c) G <,ee LA,01 = t g <,ee L,00 e iω0τ uG a,ee A,11 + vG a,he A,11 + g r,ee L,00 e iω0τ uG <,ee A,11 + vG <,he G <,eh LA,01 = t g <,ee L,00 e iω0τ uG a,eh A,11 + vG a,hh A,11 + g r,ee L,00 e iω0τ uG <,eh A,11 + vG <,hh A,11 . (S-45h) Substituting Eq. (S-45) into Eq. (S-43) and using g ee L,00 = g hh L,00 = g L , we have Furthermore, by using the relation G < − G > = G a − G r , we obtain The current can be written more compactly as with Γ u = Γ|u| 2 and Γ v = Γ|v| 2 . Note that the term x is evaluated self-consistently using or Eq. (10) of the main text: In this section, we evaluate the expressions for the ABS lesser and greater Green's functions G <,> A (ω) which are used to calculate the current [Eq.
To evaluate G <,> A (ω), we begin by writing the ABS Green's function in the Lehmann representation as where Φ + = (1, 0) T and Φ − = (0, 1) T are the positive-and negative-energy eigenfunction of the ABS written in the Nambu basis (γ,γ † ) T . The ABS's self energy due to the lead coupling is Σ r A (ω) =ť † diag(g r L (ω − ), g r L (ω + ))ť wherě t = t u v −v * −u * is the hopping matrix, and g r L (ω − ) = g r L (ω + ) = −iπν 0 is the lead retarded Green's function with ω ± = ω ± eV . Similar relations apply for Σ a,<,> A (ω). The ABS lesser Green's function is [85] where G r,a A = g r,a A (1 − g r,a A Σ r,a A ) −1 , g < j (ω) = f (ω)(g a j − g r j ) and g > j (ω) = −(1 − f (ω))(g a j − g r j ) with j = L, A. The explicit expressions of the matrix elements of G < A (ω) = can be evaluated as where ω ± = ω ± eV and with Γ u = Γ|u| 2 and Γ v = Γ|v| 2 . The expressions for the matrix elements of the ABS greater Green's function G > A (ω) can be obtained from Eq. (S-57) by using the substitutions:  Fig. 3 of the main text, we have shown that the PHS breaking holds for the case of low-frequency bosons, here we will show that it also holds for the case of high-frequency bosons, i.e., Ω > 2ε A + k B T . Unlike the perturbative calculation in the rate equation, the PHS breaking calculated from the Keldysh approach arises due to non-perturbative effects of tunneling, i.e., the PH asymmetry of the mean-field boson displacement value x . In the non-perturbative regime, electrons can tunnel from the lead into virtual states in the superconductor by emitting or absorbing bosons with high frequencies where energy violation is allowed for sufficiently large tunnel coupling (Γ Ω), resulting in PHS breaking of subgap conductances. This energy violation is allowed as long as the energy violation in the first tunneling process is negated by the second tunneling process which the conserves the total energy of a full cycle of transferring a pair of electrons in the two-step tunneling process. Figure S4(b) shows that the magnitude of the conductance PH asymmetry ζ has a nonmonotonic dependence on the ABS-boson coupling strength λ. The conductance PH asymmetry ζ first increases with increasing λ where the two peaks approach each other until they reach a certain minimum distance. Note that for this range of λ, the higher peak is at positive voltage for the case where |u| 2 > |v| 2 while for the case where |v| 2 > |u| 2 , the higher peak is at negative voltage. After the PH asymmetry reaches a maximum, it decreases to zero and stays there for a range of λ where the peaks remain more or less at the same place. As λ increases and reaches a certain value, the high and low peaks switch positions, i.e., from negative to positive voltage and vice versa. As λ keeps increasing, the two peaks move away from each other and the magnitude of the PH asymmetry increases. For large enough λ, the positions of the conductance peaks are no longer PH symmetric [see purple curve in Fig. S4(b)]. Note that the results for large λ may not be reliable as our mean-field treatment of interactions may break down in this regime. Figure S4(d) shows that the conductance PH asymmetry ζ decreases with increasing temperature T due the temperature broadening of the conductance peaks. The dependence of the ABS conductance on the boson frequency Ω is shown in Fig. S4(f). The PH asymmetry ζ has a nonmonotonic behavior with the boson frequency Ω where its magnitude first decreases to zero with increasing Ω. This corresponds to the two conductance peaks moving towards each other as Ω increases which is due to the decrease in the effective ABS-boson coupling strength λ/Ω. After the PH asymmetry reaches zero, it switches sign which corresponds to the high and low peaks switching sides. As Ω increases, the two peaks move towards each other and the PH asymmetry increases to a certain maximum value. Having reached its maximum, the PH asymmetry ζ decreases with increasing Ω which corresponds to the decrease in the effective ABS-boson coupling strength λ/Ω. Figure S4(h) shows the dependence of the conductance on the lead tunnel coupling Γ. As shown in the inset of panel (h), the PH asymmetry ζ has a non-monotonic dependence on the lead-tunnel coupling Γ where it first decreases as Γ increases. After the PH asymmetry reaches zero, it changes sign and increases in magnitude to a certain maximum value as Γ increases. Having reached its maximum, the PH asymmetry then decreases as Γ increases. Note that unlike the rate equation, our mean-field Keldysh approach shows that the conductance in the tunneling limit (Γ/Ω 1) still exhibits PH asymmetry even for high-frequency bosons. Since the treatment of interactions within the rate equation is exact in the tunneling limit, our Keldysh results obtained using the mean-field treatment of interactions may not be correct in this tunneling limit. This is because the meanfield approximation breaks down in this limit due to the singularity in the tunneling density of states. For the case where the tunnel coupling is not too small, the mean-field approximation is valid and we can see from Fig. S4(d) that unlike the rate equation, that subgap conductance calculated from the Keldysh approach can still be PH asymmetric for high-frequency boson case. Finally we note that since we ignore the Fock term in the mean-field approximation, the conductance calculated from the Keldysh approach has no boson sidebands. This tunneling Hamiltonian can be obtained by first projecting the microscopic Hamiltonian [Eq. (S-6)] onto the lowest and second-lowest energy sector α, β = 1, 2 and followed by integrating out the second-lowest Bogoliubov operator γ 2 from the total Hamiltonian of the system. Projecting the ABS and tunneling Hamiltonian onto the lowest and second-lowest energy state giveŝ αβ γ α γ β + H.c. (b † +b) + Ω(b † +b), (S-60a) For simplicity, we will choose parameters such that the tunneling term into the lowest Bogoliubov (ABS) operator (γ 1 and γ † 1 ) vanishes where we only have the tunneling into the second lowest Bogoliubov operator (γ 2 and γ † 2 ), i.e., H T =tĉ † L ũγ 2 +ṽγ † 2 + H.c. (S-61) Note that in the above, we have definedtũ ≡ t 2 u 22 (0) + t 1 u 12 (0),tṽ ≡ t 2 v 22 (0) + t 1 v 12 (0), and we have also chosen parameters such that the tunneling term into the ABS vanishes, i.e., t 2 v 21 (0)+t 1 v 11 (0) = 0 and t 2 u 21 (0)+t 1 u 11 (0) = 0. We note that we choose these parameters only for simplicity and our results on PHS breaking hold in general even without this simplification.