{'fraction_non_alphanumeric': 0.05329757814464091, 'fraction_numerical': 0.014665214112126037, 'mean_word_length': 4.357677342113139, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quasiparticle trapping, Andreev level population dynamics, and charge imbalance in superconducting weak links

We present a comprehensive theoretical framework for the Andreev bound state population dynamics in superconducting weak links. Contrary to previous works, our approach takes into account the generated nonequilibrium distribution of the continuum quasiparticle states in a self-consistent way. As application of our theory, we show that the coupling of the superconducting contact to environmental phase fluctuations induces a charge imbalance of the continuum quasiparticle population. This imbalance is due to the breaking of the left-right symmetry in the rates connecting continuum quasiparticles and the Andreev bound state system, and causes a quasiparticle current on top of the Josephson current in a ring geometry. We evaluate the phase dependence of the quasiparticle current for realistic choices of the model parameters. Our theory also allows one to analyze the quantum coherent evolution of the system from an arbitrary initial state.


I. INTRODUCTION
Quantum coherent superconducting circuits are among the most promising candidates for future large-scale quantum information processing devices, and the last few years have seen an enormous increase in research activity in this direction [1,2]. Their functioning is often limited by the presence of residual nonequilibrium quasiparticles, whose uncontrolled tunneling provide a severe decoherence mechanism [3][4][5][6]. Remarkably, in some cases where the parity of the quantum state matters, the presence of a single extra quasiparticle can determine the macroscopic response of the device [7,8]. The trapping of single quasiparticles in superconducting islands is known as "poisoning" [9]. Although at temperatures well below the superconducting gap ∆, such states have an exponentially small chance to exist in thermal equilibrium, they can have very long lifetimes once they are generated in a nonequilibrium process. Quasiparticle poisoning was also observed in recent experiments [10][11][12] for devices containing a short superconducting weak link with only a few transport channels. Those experiments reported the existence of long-lived nonequilibrium quasiparticles trapped in the Andreev bound states [13] formed near the weak link. Here, we refer to such a superconducting constriction as a "superconducting atomic contact" (SAC) [14].
In this paper, we provide a comprehensive theoretical framework for the understanding of the Andreev bound state population dynamics in a single-channel SAC. In our theory, transitions between different Andreev bound state configurations and their interplay with continuum quasiparticles are fully taken into account. Such transitions can be triggered, for instance, by environmental phase fluctuations or by phonon-induced processes, and also allow for a change of the fermion number parity in the Andreev levels. While the ground state has even parity, there are two spin-degenerate odd-parity Andreev bound state configurations with excitation energy relative to the ground state, where ϕ 0 is the superconducting phase difference across the contact and T the normal-state transmission probability of the contact.
The occupation of such an odd-parity state causes quasiparticle poisoning since the Josephson current is then completely blocked. Due to its long lifetime, the spin degree of freedom corresponding to the two odd-parity states has also been proposed as qubit platform [15,16].
On the other hand, the occupation of odd-parity states severely limits the operation of the "Andreev qubit" [17,18], which is built from the Andreev ground state configuration and the excited even-parity state of energy 2E A , cf. also Ref. [19]; for T → 1 and ϕ 0 ≈ π, these two states are nearly degenerate. Similar superconducting devices are also discussed in the context of Majorana fermion physics [20,21], and questions pertaining to quasiparticle poisoning and the interplay between Andreev (or Majorana) and continuum quasiparticle distribution functions are important in that direction as well. The phase-dependent energy E A in Eq. (1.1) also determines the transition frequencies between different Andreev configurations, which have recently been studied by microwave absorption and supercurrent spectroscopy [11,12], where the odd-parity states can be excited together with a continuum quasiparticle. For the theoretical description of such "artifical atom" spectra in microwave irradiated contacts, see Refs. [19,22,23] and references therein. The above discussion shows that it is an important issue to understand the Andreev bound state population dynamics and its interplay with continuum quasiparticles in SACs. We here study this problem for the simplest single-channel case, where the SAC quasiparticles interact with the phase fluctuations of an electromagnetic environment. To be concrete, we shall consider the plasma mode of another Josephson junction in the ring geometry of Fig. 1, but our subsequent discussion is formulated for a general environment. Previous work has already addressed some aspects of this problem [24], but to the best of our knowledge, the generated nonequilibrium distribution of continuum quasiparticles and its interplay with quasiparticle poisoning has not been discussed in a self-consistent manner up to now.
The structure of the remainder of this paper is as follows. In Sec. II, we introduce the model and provide a second-quantized formulation, where the fermionic (Andreev bound state and continuum) quasiparticles are weakly coupled to the environmental phase fluctuations. In Sec. III, assuming that the electromagnetic environment remains in thermal equilibrium, we provide the master equation description of this system. The density matrix for the quasiparticles can then be factorized into an Andreev part, ρ A (t), and a diagonal density matrix describing the quasiparticles belonging to the continuum spectrum. Importantly, the off-diagonal part of ρ A (t) decouples from the diagonal part, where the latter is determined by Andreev level occupation probabilities. Including quasiparticle relaxation by phonons in a phenomenological way, the resulting steady-state solution can be obtained by self-consistently solving just two coupled nonlinear equations. In Sec. IV, as an application of our approach, we describe a charge imbalance effect caused by an asymmetry in the transition rates between Andreev and continuum quasiparticles. The self-generated nonequilibrium distribution of continuum quasiparticles (we note that no external forces drive the system out of equilibrium in our model) causes a phase-dependent quasiparticle current, and an asymmetric charge profile around the weak link. The predicted charge imbalance effect could be measured by superconductingnormal tunnel junction spectroscopy [25][26][27][28]. Finally, we conclude in Sec. V. Details about our calculations can be found in several Appendices. We often use units with = k B = e = 1.

A. Model
In this work, we consider a single-channel SAC embedded in the asymmetric SQUID geometry sketched in Fig. 1, where the ring contains both the SAC and a conventional Josephson junction. This setup has also been studied in other theoretical works, for instance, see Refs. [23,24]. Denoting the superconducting phase differences across the SAC and the Josephson junction by ϕ and χ, respectively, both phases are linked by where the dimensionless parameter ϕ 0 is related to the magnetic flux (in units of the flux quantum h/2e) threading the ring. Assuming that the Josephson energy, E J , is much bigger than the charging energy, E C = (2e) 2 /2C with the capacitance C, of the Josephson junction, the environmental electromagnetic modes seen by the SAC are well described by an effective LC circuit Hamiltonian, corresponding to an undamped harmonic oscillator. In Sec. III, we shall also include the effects of an additional shunt resistance R, leading to the damping parameter η d = 1/(RC). Note that in the regime E J E C of interest here, fluctuations of χ are small, χ 2 1. Turning to the single-channel SAC, the BCS Hamiltonian is written in terms of a two-component Nambu spinor, ψ(x) = (ψ ↑ (x), ψ * ↓ (x)) T , describing electrons in the left (x < 0) or right (x > 0) superconducting bank, with the contact at x = 0. Using the standard quasiclassical Andreev approximation [2], we introduce slowly varying envelope functions, ψ(x) = e ik F x ψ R (x) + e −ik F x ψ L (x), with Fermi momentum k F . Combining the right-and left-moving envelopes into Ψ(x) = (ψ R , ψ L ) T , where each entry still carries the Nambu spinor structure, the time-dependent wave function satisfies the Bogoliubov-de Gennes (BdG) equation [2], with Fermi velocity v F , the BCS gap ∆, and Pauli matrices σ x,y,z and τ x,y,z in Nambu and right/left-mover space, respectively; the corresponding unit matrices σ 0 and τ 0 are often kept implicit below. As shown in Ref. [18] and references therein, the BdG solutions on both sides of the contact have to be matched at x = 0 by a transfer matrix, which is 4π-periodic in ϕ. For simplicity, the transmission probability, 0 < T ≤ 1, which characterizes the transparency of the constriction in the normal phase, is assumed energy-independent. In what follows, it is convenient to remove the time dependence from Eq. (2.4) by a gauge transformation, with χ(t) in Eq. (2.1). The phase factor in Eq. (2.4) thereby becomes time-independent, with ϕ(t) → ϕ 0 , and H BdG → H BdG + V , where the interaction term is, witḣ χ = ∂ t χ, given by Since χ 2 1, the linearized expression in Eq. (2.6) is sufficient, which now couples the quasiparticle dynamics to the phase χ(t). Using the Josephson plasma frequency, Ω = √ 2E C E J , where we assume Ω < ∆ throughout this paper, the Lagrangian of the coupled system, with Ψ = (Ψ * ) T , is (2.7) Employing the momentum P χ canonically conjugate to the phase χ, the corresponding Hamiltonian is (2.8)

B. Second-quantized formulation
We now switch to a second-quantized language by letting Ψ(x) →Ψ(x), where the electron field operator,Ψ(x), is expanded in terms of the stationary solutions, Ψ ν (x), with energy E ν , of the BdG equation for time-independent matching condition (2.4), i.e., for ϕ(t) = ϕ 0 . The wave functions Ψ ν (x) thus represent the noninteracting SAC eigenstates. Introducing the corresponding quasiparticle creation (annihilation) operators γ † ν (γ ν ), with the standard fermionic anticommutator algebra {γ ν , γ † ν } = δ νν , we arrive atΨ(x) = ν Ψ ν (x)γ ν . The noninteracting SAC Hamiltonian then reads The quantum numbers ν include (i) a pair of Andreev bound states, ν = η = ±, where the energy E η = ηE A , with E A (ϕ 0 ) in Eq. (1.1), is within the BCS gap and Ψ η (x) stays localized near the contact at x = 0, and (ii) delocalized scattering states in the continuum, ν = p = (E, s), where |E| ≥ ∆ and the index s (with s = 1, 2, 3, 4) refers to the four possible types of incoming states (from the left or right side, and of electron-or holelike character). The wave functions Ψ ν (x) are provided in analytical form in Appendix A, see also Ref. [24]. We here employ a semiconductor representation, where quasiparticles are effectively spinless but can have either positive or negative energy. In the ground state of H SAC , all E ν < 0 states are occupied, including the η = − Andreev bound state. Using standard occupation number operators,n ν = γ † ν γ ν , with eigenvalues n ν = 0, 1, the four possibilities for the occupation of the Andreev bound state sector are indexed by (n + , n − ). The ground state, with energy −E A , corresponds to the (0, 1) configuration, which we also denote by the Andreev state |− A . This state carries the equilibrium Josephson supercurrent In the parlance of Refs. [10,24], |− A is an even-parity state, while the odd-parity sector corresponds to the spin-degenerate (0, 0) and (1, 1) states, with excitation energy E A relative to the ground state. The odd-parity states with n + + n − = 0 and 2, resp., are denoted by and imply a vanishing Andreev supercurrent, consistent with the "quasiparticle poisoning" scenario. The lifetime of these states can reach the millisecond regime for high transparency T , and their decay rate exhibits nearly universal scaling as a function of E A /∆ [10,24]. Finally, the (1, 0) even-parity state, denoted as represents an excited "Andreev Cooper pair" localized at the contact, with excitation energy 2E A above the ground state. The |+ A state carries the Josephson current −I A , with opposite sign as compared to |− A , but rather quickly relaxes to the ground state [24]. The second-quantized form of the interacting Hamiltonian (2.8) is thus given by 13) are discussed below and in App. B. For convenience, we now shift P χ → P χ −Ā by means of a unitary transformation, H → U HU † with U = e iĀχ , and represent the (unitarily transformed) phase χ and its momentum P χ by a standard boson operator, b, with commutator After some algebra, we thereby arrive at the Hamiltonian in its final form (up to an irrelevant constant), (2.14) describing fermionic (Andreev level and continuum) quasiparticles coupled to an oscillator mode with the plasma frequency Ω. For E J E C , we are effectively in the weak-coupling regime, λ 1, with the dimensionless coupling strength λ = E C /4Ω. Finally, the Josephson current operator in Eq. (2.14) iŝ where the matrix elements I νν are discussed in App. B and in the next subsection.

C. Current
We first note that due to the spatial homogeneity of the extended quasiparticle states (plane waves) away from the contact, the matrix elements A pp and W pp , and hence also I pp , between continuum states can be finite only when their energies match, E = E , i.e., phase fluctuations do not induce intraband transitions. As we show in App. B, one finds that W pp = 0 even for E = E , implying that I pp = 0. Therefore, delocalized continuum states can contribute to the supercurrentÎ S only through transitions mixing them with Andreev levels.
The Josephson current operator then contains a part I A , coming from the Andreev sector only, and a partÎ cA , describing the mixing of continuum and Andreev states, I S =Î A +Î cA . We find for the pure Andreev current [17,18] where γ = (γ + , γ − ) T combines the two Andreev level fermion operators, and the Pauli matrices η x,y,z act in the corresponding space. Note that the Andreev current operator (2.16) is written in energy representation, where the Hamiltonian projected to the Andreev sector is diagonal, For non-ideal transparency of the contact, T < 1,Î A does not commute with H A -Andreev level eigenstates are superpositions of current eigenstates implying that strong fluctuations of the supercurrent are generated for ϕ 0 ≈ π [29]. Similarly, the supercurrent contribution caused by the mixing of continuum and Andreev states iŝ where the matrix elements I η,p are specified in App. B.
Finally, the total current flowing through the contact also contains a conventional dissipative quasiparticle contribution due to continuum states, I qp , on top of the supercurrent contribution Î S . We provide the standard scattering theory expression for I qp in App. B.
To study the physics described by the interacting Hamiltonian, H = H 0 + V , with the noninteracting piece H 0 = H SAC +Ωb † b and the interaction contribution V = λ(b + b † )Î S , we now turn to a master equation approach. In this work, we assume that the plasma mode remains in thermal equilibrium with a heat bath of temperature T env at all times, and we thus neglect feedback effects on the phase dynamics.

A. Master equation
Within the master equation framework [2], the Liouville-von Neumann equation for the density matrix of the complete system, ρ tot , is expanded to second order in the small interaction parameter λ 1. Writing time-dependent operators in the interaction picture as (3.1) Our assumption of thermal equilibrium for the plasma mode implies a factorized form of the density matrix, where ρ osc ∼ e −(Ω/Tenv)b † b is a thermal density matrix for the plasma mode and ρ(t) describes the time evolution of fermionic quasiparticles. Taking the trace over the oscillator degree of freedom, Eq. (3.1) yields where we have employed the Markov approximation, valid at long times t and not too low temperatures [2,30]. The boson correlator in Eq. (3.3) reads  5) and the environmental spectral density, (3.6) We use Eq. (3.6) below also for ω < 0, and directly include the Ohmic damping parameter, η d , to capture the effects of a shunt resistance, see Sec. II. For η d → 0, the spectral density has the limit J(ω) = λ 2 δ(|ω| − Ω)sgn(ω). For finite η d , Eq. (3.6) exhibits sharp peaks for |ω| = Ω.
The equation of motion (3.3) is still quite cumbersome, and we shall here proceed by making two approximations. First, we neglect entanglement between the Andreev and continuum quasiparticles, which means that the reduced density matrix factorizes into an Andreev part and a continuum part, This approximation is justified in the weak-coupling regime λ 1, since higher-order terms in λ are needed to coherently couple Andreev and continuum states [31]. The factorized density matrix (3.7) is expected to be highly accurate away from the zero-temperature limit, since the thermal energy uncertainty causes a blurring of continuum quasiparticle wavepackets that rapidly destroys entanglement between Andreev and continuum states. Second, we also assume that the density matrix ρ c (t) describing continuum quasiparticles remains diagonal during the time evolution. This approximation is justified by noting that there are no direct matrix elements in H connecting different continuum states, and implies that ρ c (t) is fully determined by specifying the time-dependent occupation probabilities n p (t) of continuum states, where |1 p = γ † p |0 p corresponds to a filled single-particle state p = (E, s). Note that ρ c (t) in Eq. (3.8) is always normalized, Tr c [ρ c (t)] = 1. On the other hand, the density matrix ρ A (t) describing the Andreev sector, with normalization condition Tr A [ρ A (t)] = 1, may have offdiagonal entries reflecting quantum coherence.
Tracing over the Andreev part in Eq. (3.3) then yields an equation of motion for the continuum state occupation numbers n p (t). Similarly, tracing instead over the continuum states, one obtains an equation for the time evolution of the reduced Andreev density matrix ρ A (t). In these equations, the transition rates between different levels follow as Fermi golden rule expressions, 9) with the Bose function n B (ω) in Eq. (3.5) and the spectral density J(ω) in Eq. (3.6). By using Eq. (2.16), we observe that the direct rates connecting different Andreev states are given by Notice that these rates vanish for perfect transparency, T → 1. Recalling now that I pp = 0 for arbitrary T , we see that transition rates between continuum states are always absent, Γ pp = 0. Finally, the supercurrent matrix elements between Andreev and continuum states, I ηp , see Eq. (2.15) and App. B, determine the corresponding transition rates, Γ η,p , for exciting an Andreev quasiparticle into the continuum, plus the reverse process with rate Γ p,η . Such transitions must involve the absorption or emission of an environmental photon. Since |E| ≥ ∆ and the spectral density is sharply peaked around the Josephson plasma frequency Ω, those rates are sizeable only when Ω > ∆ − E A [24].
Performing now the trace over the Andreev sector in Eq. (3.3), we find (3.11) The time-dependent continuum state distribution function, {n p (t)}, thereby couples to the Andreev level occupation probabilities, Tracing instead over the continuum states in Eq. (3.3), we find ({A, B} denotes the anticommutator) This equation has been obtained under the assumption that the coupling to the plasma mode provides the only relaxation mechanism, but in Eq. (3.23) below, we will also include the effect of other mechanisms (e.g., phonons) through a phenomenological relaxation time τ qp . Notice that the terms ∼ γ † η ρ A (t)γ η and ∼ γ η ρ A (t)γ † η in Eq. (3.13) describe "parity jumps" where the fermion number parity of Andreev quasiparticles can change.
Since there are four Andreev configurations (n + , n − ), the Andreev density matrix is a 4 × 4 matrix. We here represent ρ A (t) in the basis spanned by the Andreev ground state |− A , corresponding to the (0, 1) configuration, the spin-degenerate odd-parity states |0 A and |2 A in Eq. (2.10), and the excited even-parity state |+ A in Eq. (2.11). The diagonal elements of ρ A (t) yield the respective occupation probabilities, P 0 (t) = A 0|ρ A (t)|0 A , and likewise for P η=± (t) and P 2 (t). Thereby the normalization condition for ρ A (t) gives P 0 (t) + P 2 (t) + η P η (t) = 1, (3.14) and the n η=± (t) in Eq. (3.11) are expressed as n η (t) = P η (t) + P 2 (t). (3.15) We now observe that the off-diagonal components of ρ A (t) decouple from the equations for the diagonal part in Eq. (3.13); we briefly discuss the dynamics in the off-diagonal sector in Sec. III B. The diagonal part of Eq. (3.13) determines the dynamics of the Andreev state occupation probabilities, where we finḋ Together with Eq. (3.11), we thereby arrive at a set of coupled nonlinear equations determining the timedependent continuum distribution function, {n p (t)}, and the Andreev level probabilities, P ±,0,2 (t). Importantly, despite of the approximations involved in their derivation, these coupled equations automatically satisfy the normalization condition (3.14). The resulting Andreev bound state population dynamics is schematically illustrated in Fig. 2. The rates Γ η,−η [Eq. (3.10)] connect the even-parity Andreev states |η = ± A , without involving continuum quasiparticles. However, processes that populate or depopulate odd-parity Andreev states sensitively depend on the continuum distribution function {n p (t)}.

B. Off-diagonal part of Andreev density matrix
Within our approach, off-diagonal components of the Andreev density matrix ρ A (t) decouple from the diagonal ones and obey their own set of dynamical equations. Although later on, we restrict ourselves to diagonal initial density matrices, where off-diagonal components do not appear at all, future experiments may test the corresponding quantum-coherent correlations in the Andreev sector. Anticipating such experiments, which would probe an NMR-like quantum dynamics within the Andreev sector after careful choice of the initial conditions, we provide the relevant equations now. Taking into account Hermiticity of the Andreev density matrix ρ A (t), we find from Eq. (3.13) that ρ 0,2 (t) ≡ A 0|ρ A (t)|2 A = ρ 0,2 (0), while the dynamics of the even-parity matrix elements is determined by The dynamics of matrix elements connecting states with different parity follows from In the remainder of the present work, however, we assume that the initial state (at t = 0) is diagonal. In that case, the decoupled off-diagonal density matrix elements remain zero during the entire time evolution.

C. Steady-state distribution of quasiparticles
We now proceed under the assumption that the initial Andreev density matrix, ρ A (0), is diagonal in the basis {|± A , |0 A , |2 A }. In the long-time limit, the system will then reach a time-independent steady-state distribution, which is fully characterized by the probabilities P ±,0,2 together with the continuum quasiparticle distribution function {n p }. In order to determine these quantities, we first observe that P 2 = P 0 due to the spin degeneracy of the two odd-parity states. Moreover, by using the normalization condition in Eq. (3.14), P 0 can be expressed in terms of P ± alone, For the Andreev level occupations, we thus find 22) and the steady-state version of Eq. (3.11) yields (3.23) where we added a phenomenological relaxation term for continuum quasiparticles describing, for instance, the effect of phonons [24]. According to the estimates detailed in Ref. [24] for SACs made of aluminum, we expect τ qp ∆ ≈ 10 4 in the phonon-dominated regime (given by E A < ∆ − Ω, see below). For simplicity, we here assume an energy-independent relaxation time, τ qp , and a Fermi distribution function for the noninteracting continuum quasiparticles, where the temperature T qp may differ from the temperature T env governing environmental phase fluctuations. We mention in passing that the theory in Ref. [24] corresponds to the fast equilibration case with Γ p,η τ qp 1. Taking into account Eq. (3.21), the rate equation (3.16) then yields the steady-state relation (3.25) and Eq. (3.17) is automatically fulfilled.
It is now a simple matter to solve Eq. (3.23) for the continuum quasiparticle distribution function, (3.26) which is thereby expressed by the P ± -dependent effective ratesΓ To obtain the Andreev level probabilities P ± , we then insert Eq. (3.26) back into Eq. (3.25). After some algebra, we arrive at two coupled nonlinear equations, with the auxiliary functions 29) where p = (E, s) and ν 0 = L/(π v F ) is the normal density of states at the Fermi level. The nonlinear system in Eq. (3.28) can be solved by numerical iteration, where a relative accuracy of 10 −6 was ensured by using a Newton-Raphson algorithm. This is necessary because the continuum quasiparticle distribution {n p }, which follows by virtue of Eq. (3.26) from the self-consistent solution for P + and P − , strongly responds even to tiny changes in the P ± . Below, it will be useful to consider the rate Γ in for transitions from the even-parity to the odd-parity sector (i.e., P 0 increases), as well as the escape rate, Γ out , out of the odd-parity state (i.e., P 0 decreases). Assuming an equilibrium quasiparticle distribution function {n p }, those rates were considered in Ref. [24]. Here, by employing the self-consistent continuum quasiparticle distribution function, both rates can be read off from Eq. (3.17), (3.30) As observable of primary interest, we will discuss the quasiparticle current I qp , which follows with our selfconsistent solution for {n p } by using standard scattering theory expressions. We summarize these for the convenience of the reader in App. B.

D. Perfect transparency
As application of our theory, we will discuss a charge imbalance effect in Sec. IV. This discussion is simplified when considering a SAC with perfect transparency, T = 1. We specify the explicit form of the Andreev and continuum state wave functions for T = 1 in App. C. The Andreev bound state energies, ηE A with η = ±, then follow from E A (ϕ 0 ) = ∆| cos(ϕ 0 /2)|, see Eq. (1.1), and for ϕ 0 → π, the Andreev levels tend to zero energy. Moreover, Eq. (3.10) shows that transition rates between different Andreev states vanish for T = 1, i.e., Γ η,−η = 0.
We show in App. C that for given energy E with |E| ≥ ∆, there are two decoupled types of scattering states Ψ p=(E,s) , namely s = {1, 4} and s = {2, 3}. Those channels correspond to a net charge transfer across the weak link in opposite directions. Indeed, charge flows from the left to the right side for s = {1, 4}, but from the right to the left when s = {2, 3}, as is directly seen from the definition of the scattering states, see Eqs. (A4) and (A6). This also implies that the supercurrent matrix elements between Andreev and continuum states, I p,η , are nonzero only when η = −sgn(π − ϕ 0 ) for s = {1, 4}, and η = +sgn(π −ϕ 0 ) for s = {2, 3}. In what follows, we take the phase difference across the contact as 0 ≤ ϕ 0 ≤ π.
For n R (E) = n L (E), continuum quasiparticles are driven out of equilibrium. For given steady-state Andreev occupation probabilities P ± , the distribution functions in Eq. (3.32) follow from Eqs. (3.26) and (3.27), taken with the above T = 1 rates.

IV. CHARGE IMBALANCE EFFECT
We now discuss a charge imbalance effect which is predicted to be observable in high-transparency SACs. We shall discuss this effect for a perfectly transmitting SAC, T = 1, and by assuming ϕ 0 ∈ [0, π]; for ϕ 0 ∈ (π, 2π), the sign of the induced quasiparticle current discussed below is reversed. Noting that our theory allows for arbitrary 0 < T ≤ 1, we find very similar results also for (not too small) T < 1 and observables taken as function of E A /∆. We then put T = 1 from now on.
Using Eqs. (3.31) and (B4), the quasiparticle current flowing through the SAC is given by with the energy-resolved dimensionless quasiparticle current (|E| ≥ ∆), and the self-consistent distribution functions n L,R (E) in Eq. (3.32). Evidently, if a charge imbalance is present, n L (E) = n R (E), one generally expects a finite quasiparticle current I qp from Eq. (4.1). We also define the total accumulated quasiparticle charge, (4.3) Since the density of states ν 0 ∝ L/ξ 0 , where L is the channel length and ξ 0 = v F /∆ the BCS coherence length, Q qp tends to vanish for a very short channel, L/ξ 0 → 0, while the induced quasiparticle current remains finite in that limit.
Let us first address the steady-state Andreev populations, P ±,0,2 , where P + and P − follow from the selfconsistent solution of Eq. (3.28), and the occupation probability of the degenerate odd-parity state, P 0 = P 2 , is then given by Eq. (3.21). Representative results for P ±,0 vs E A /∆ are shown for experimentally relevant parameters in Fig. 3. Since E A = E A (ϕ 0 ), see Eq. (1.1), Fig. 3 essentially shows the phase dependence of the Andreev state probabilities for ϕ 0 ∈ [0, π]. The charge imbalance turns out to be absent in the strong relaxation regime τ qp ∆ < 1 (see also below), where our the- ory reduces to the approach of Ref. [24] and thus selfconsistency plays no role. We therefore focus on the weak relaxation regime τ qp ∆ 1 in this section. The main panel in Fig. 3 is for T env = T qp , while the inset studies a case where T env substantially exceeds T qp . From Fig. 3, we can distinguish two qualitatively different regimes, E A > ∆ − Ω (with P 0 → 0) and E A < ∆ − Ω (with P 0 = 0), respectively. For E A > ∆ − Ω, the system remains quite close to the ground state, |− A , since environmental photons can rapidly excite quasiparticles from an odd-parity state into the continuum. On the other hand, for E A < ∆ − Ω, the frequency Ω is too low to achieve such a transition.
The corresponding rates Γ in and Γ out , see Eq. (3.30), for populating and depopulating the odd parity states, resp., are shown in Fig. 4, again as a function of E A /∆. We observe now more clearly that E A > ∆ − Ω and E A < ∆−Ω correspond to qualitatively different physical regimes. For E A ≈ ∆ − Ω, the rates increase over several orders in magnitude with very small ϕ 0 variation, and one enters a regime where the odd-parity state quickly decays. This regime, E A > ∆ − Ω, has been termed "fast relaxation regime" in Refs. [10,24].
Next, in Fig. 5 we show the induced quasiparticle current I qp , see Eq. (4.1), vs E A for Ω = 0.5∆. This quantity clearly demonstrates that there is a significant charge imbalance effect throughout the regime E A > ∆ − Ω, but not for E A < ∆−Ω. The induced current gets reduced as the quasiparticle relaxation rate 1/τ qp increases, and is only significant for τ qp ∆ 1, which is the typical regime for SACs made of aluminum [10]. Further insight on the generated charge imbalance is obtained by analyzing the distribution functions n R,L (E) for right-and left-moving quasiparticles, see Eq. (3.32). As illustrated by the insets in Fig. 5, the generated imbalance is maximal for E A → ∆ (i.e., for ϕ 0 → 0, where the supercurrent Î S can be vanishingly small), and becomes smaller as E A de- creases. For the present case with T qp = T env , the smaller n(E) curve (i.e., the n R component for E > 0, and 1−n L for E < 0) coincides with the Fermi distribution at the corresponding temperature, indicated by dotted curves in the insets. Noting that the Josephson current for a fully transparent SAC is of order Î S = Î A ≈ e∆/ , the induced quasiparticle current is a few percent of this value for the parameters in Fig. 5. In the ring geometry of Fig. 1, the quasiparticle current I qp flows in opposite direction to the Josephson current Î S , which can be rationalized as follows. The rate from |+ A to the leftmoving s = 2 continuum states with E > 0 -carrying negative current -is much larger than the one from |− A to the (E > 0, s = 1) states carrying positive current, due to the much shorter distance in energy. This difference is able to outweigh the fact that P + < P − favors the same sign of I qp and Î S . The parameters considered up to now were inspired by those realized in available experimental reports [10,24]. However, as we show next, it is also interesting to consider an alternative scenario, where the temperature of the environmental modes is so high to put them into a classical regime, T env Ω. Experimentally, such a situation can be realized by replacing the electromagnetic environment by an external microwave radiation source of frequency Ω. We here consider the case Ω = 0.2∆, with T env = 0.5∆ = 2.5Ω. For the quasiparticle temperature T qp , we take T qp = 0.01∆, significantly smaller than T env . The Andreev state populations for this case were shown in the inset of Fig. 3, and the corresponding Γ in/out rates in the inset of Fig. 4  regimes (in the parlance of Ref. [10]) can be identified for E A > ∆ − Ω and E A < ∆ − Ω, respectively. However, in this case, the generated quasiparticle populations differ more strongly from the Fermi distributions (see insets of Fig. 6), and a significant quasiparticle current is induced throughout the whole E A range. This is illustrated in the main panel of Fig. 6. It is also interesting to notice that the induced quasiparticle current exhibits a sign change for E A Ω.

V. CONCLUDING REMARKS
In this work, we have formulated and applied a theoretical framework for the Andreev bound state population dynamics in single-channel superconducting weak links. Taking into account phase fluctuations by an electromagnetic environment, we have developed a master equation approach for the quasiparticle dynamics, capturing the rich interplay between Andreev states and continuum states. In particular, the role of odd-parity Andreev states and the need for a self-consistent treatment of the generated nonequilibrium continuum quasiparticle distribution has been emphasized. As an application of our theory, we have demonstrated that the asymmetry in Andreev-continuum quasiparticle transition rates causes an intriguing charge transfer across the weak link, reflected in a quasiparticle current. Using established experimental techniques, this charge imbalance should be measurable in a SAC. Our theory could be also applied for the study of the quantum coherent dynamics of this system, including the effect of parity mixing processes. This is of relevance for the various proposals of using Andreev levels as qubits [15,16,18,19]. As another ex-tension of our formalism, it would be very interesting to study the Andreev-and Majorana bound state dynamics in topological superconductor weak links, or to study the interaction-induced effects (see also Ref. [33]) on Andreev bound state dynamics when the constriction contains a quantum dot with sizeable charging energy, or couples to local phonon modes.
Note: During the preparation of this manuscript, we became aware of related work [34], where the charge imbalance effect described here has also been pointed out. Where there is overlap, our results match theirs. However, in contrast to our work, Ref. [34] considers a monochromatic external microwave driving and the continuum quasiparticle distribution was not treated in a self-consistent manner.