Spectral response of Josephson junctions with low-energy quasiparticles

Anna Keselman, 2 Chaitanya Murthy, Bernard van Heck, 4 and Bela Bauer Station Q, Microsoft Corporation, Santa Barbara, California 93106 USA Kavli Institute for Theoretical Physics, University of California, Santa Barbara California 93106 USA Department of Physics, University of California, Santa Barbara, CA 93106 Microsoft Quantum Lab Delft, Delft University of Technology, 2600 GA Delft, The Netherlands (Dated: May 10, 2019)


Introduction
Due to the macroscopic coherence of the superconducting state and their non-linearity as circuit elements, Josephson junctions are a workhorse of quantum state engineering [1]. They are the fundamental building block of superconducting qubits [2,3] like the transmon [4] and the fluxonium [5]. In these quantum engineering applications, the superconducting circuit embedding the Josephson junction [6] is operated at frequencies ω ∼ 5-10 GHz, which are much smaller than the superconducting gap of the electrodes, ∆/h 50 GHz for aluminum. As a consequence, quasiparticles are not involved in the coherent dynamics of the circuit, although their presence influences the relaxation and dephasing of superconducting qubits [7][8][9][10][11][12][13][14][15][16][17].
New frontiers in superconducting devices force us to reconsider the role of quasiparticles in Josephson junction dynamics. Most notably, the presence of Majorana zero modes (MZMs)-topologically protected zero-energy quasiparticles that emerge at the ends of topological superconducting wires [18,19]-can drastically affect the behavior of a superconducting circuit. MZMs are able to non-locally encode qubits [18] and, via non-Abelian braiding [20], allow fault-tolerant processing of quantum information. Therefore, they form a potential platform for topological quantum computation [21] which is actively being pursued [22]. A junction between two topological superconductors exhibits a 4πperiodic Josephson effect [18,23,24], a hallmark feature of topological superconductivity which has been experimentally sought [25,26]. Several practical schemes for topological quantum computation rely on the coupling of MZMs across a Josephson junction and on the use of microwave circuits for control and readout of topological qubits [27][28][29][30]. In conjunction with the growing interest in MZMs, different research groups have developed and studied superconducting devices with semiconductor-based Josephson junctions either in nanowires [31,32] or 2DEGs [33], as well as in graphene-based heterostructures [34,35]. These junctions, characterized by few conducting channels and potentially high transparency, can also be used for the development of qubits based on conventional Andreev bound states (ABSs) [36][37][38][39][40]. The presence of low-energy ABSs in nanowire-based junctions has been directly measured via microwave spectroscopy [41][42][43].
Understanding present and future experimental developments in this direction calls for adequate theoretical approaches that can fully incorporate the role of quasiparticles in the circuit dynamics. In most theoretical descriptions of Josephson-junction dynamics, quasiparticles are included as a fermionic bath [44] (see Ref. [45] for a recent exception). In the case of Andreev qubits, detailed models which are amenable to an analytical approach are available in simple limits, such as that of a short Josephson junction with a single conducting channel [37,[46][47][48]. On the other hand, most of the theory literature treating the presence of MZMs in superconducting circuits [49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66] relies on simple toy-models with phenomenological terms representing Majorana couplings, bypassing a microscopic description of the topological phase.
In this paper, we carefully examine this problem using accurate numerical simulations of a microscopic model for a one-dimensional topological superconductor. While we confirm the applicability of simplified models in certain limits, we find that in other, experimentally relevant limits they are insufficient to describe the system's behavior. We focus on the topological phase transition and on the case of additional subgap Andreev bound states in the junction, which we expect to generically appear in wires with large spin-orbit coupling and external magnetic field. We find that the phase transition does not lead to strong signatures in the response of the capacitively shunted junction, while additional subgap Andreev states exhibit complex interplay with the plasma modes and significantly alter the response.
Our numerical simulations are based on matrix product states (MPS) [67][68][69]. Specifically, we use the density matrix renormalization group (DMRG) [70][71][72] and time-evolving block decimation (TEBD) [73][74][75][76] to compute the time-dependent charge correlation function of a nanowire Josephson junction shunted by a large capacitor (i.e. a transmon circuit) across the topological phase transition. In the frequency domain, the correlation function determines the observed spectra in a typical circuit QED (cQED) experiment, making our method suitable for direct comparison with experimental measurements. This approach allows us to determine the expected frequency spectra even close to the critical point-a regime which cannot be captured by existing toy models-and to easily include additional Andreev bound states.
In order to interpret the results of the MPS simulations, and extending previous studies [48], we also develop a simple perturbative approach which is valid in the harmonic limit, i.e. when the charging energy is small compared to the Josephson energy. The method allows one to derive an effective Hamiltonian for the capacitively shunted junction, starting from an arbitrary quadratic Hamiltonian describing the quasiparticles. The effective Hamiltonian takes the form of a generalized Jaynes-Cummings model describing the interaction between Josephson plasma modes and quasiparticle excitations. This model describes the energy spectra obtained from the MPS simulations deep in the harmonic limit quite well, but cannot reproduce non-perturbative effects that arise away from this limit (and are fully captured by the MPS simulations), such as the charge dispersion of energy levels and certain couplings between plasma modes and fermionic modes.
The paper is structured as follows. In Sec. 2 we present the setup and the general model used to describe a nanowire-based Josephson junction shunted by a capacitor, incorporating the fermionic degrees of freedom. In addition, we discuss the experimentally relevant probes and the parameter regimes we will be addressing in this study. As a simple application of the general model, in Sec. 3 we discuss and review a minimal model with a single low-energy fermionic mode on each side of the junction, which captures the essential ingredients of a nanowire in the topological phase. In Sec. 4 we discuss how MPS-based techniques can be used to calculate the experimentally relevant quantities that probe the response of the system. We then introduce the microscopic model for the nanowire that we use in our numerical study, and present the results. In Sec. 5 we consider the limit of small charging energy, and derive an effective theory based on a perturbative expansion which successfully captures the coupling between the plasma mode and the fermionic quasiparticles in this limit. We then discuss the effect of non-perturbative corrections. Figure 1: Schematic setup: a semiconducting nanowire proximity coupled to a superconductor, forming a Josephson junction shunted by a capacitor and controlled by the gate voltage V g . The right half of the wire is connected to a superconducting ground. The dashed line represents the partitioning of the system used in our theoretical model, see main text for details.

Setup and model
The setup we consider is schematically depicted in Fig. 1. A semiconducting nanowire is proximity-coupled to a grounded superconductor on its right half, and to a floating superconducting island on its left half. A short segment in the middle of the wire, which is not in direct contact to any superconductor, forms a junction between the two superconductors. The conductance of the junction can be tuned by a gate underneath this middle region. The voltage on this gate is denoted by V b and determines the strength of the Josephson coupling between the floating superconducting island and the grounded superconductor. The island is shunted by a large capacitance C to the ground. The charge induced on the island is controlled by a gate with voltage V g .
The Hamiltonian describing the system is given by The first term above is the electrostatic energy of the island, with charging energy E c = e 2 /2C and dimensionless gate charge N g = C V g /e, where e is the electron charge. The electrostatic energy is determined by the total number of electrons on the island (counted from the neutrality point), which is the sum of the number of paired electrons in the superconductor, N , and of the number of electrons in the left segment of the semiconducting wire, N f,L (to be better specified below). The second term in Eq. (1) describes the dynamics of the fermionic degrees of freedom in the semiconducting wire. The dynamics are prescribed by a Bogoliubov-de Gennes (BdG) Hamiltonian, H BdG , which includes the coupling between the semiconductor wire and the superconductors as well as the coupling between the two wire segments across the junction. For concreteness, we consider a lattice description of the system, such that H BdG is written in the Nambu basis where c † i,σ (c i,σ ) is the creation (annihilation) operator of an electron on site i with spin σ.
To simplify the treatment of the charging energy, we consider a sharp boundary between the left part of the wire, which is coupled to the floating superconducting island, and the right part of the wire, which is coupled to the grounded superconductor. We choose to place this sharp boundary at the left end of the junction, as indicated by the dashed black line in Fig. 1. Although this choice can have a quantitative effect on the spectrum of the full Hamiltonian (with other parameters held fixed), we expect the qualitative physics to be insensitive to a specific (but generic) choice for the position of the boundary.
The sites i in the lattice description of the wire are divided into two sets, I L and I R , depending on whether they belong to the left or to the right part of the wire. The number of electrons in the left part of the wire is defined as The induced s-wave superconductivity is included in H BdG via pairing terms of the form ∆e iφ c i,↑ c i,↓ +h.c. (if i belongs to the part of the wire coupled to the floating island) or ∆c i,↑ c i,↓ +h.c. (if i belongs to the part of the wire coupled to the grounded superconductor), where ∆ is the induced superconducting gap. The pairing vanishes in the junction region. In what follows we will consider junctions of finite extent as well as junctions consisting of a single weak link. The operator e iφ (e −iφ ) adds (removes) a Cooper pair to (from) the left superconductor, and is canonically conjugate to the charge operator N , i.e.
The pairing terms in H BdG thus commute with the total charge of the floating island, N + N f,L , which enters in the charging energy in Eq. (1). On the other hand, hopping terms in H BdG which connect I L and I R do not commute with the charging energy term. More general forms of the induced pairing, e.g. a spatial variation of the pairing term strength or different pairing symmetries, could easily be included. If the fermionic quasiparticles are gapped and one is interested in the behavior of the system at frequencies ω far below the excitation energy for quasiparticles, i.e. ω ∆, then one may replace the Hamiltonian H BdG with its phase-dependent ground state energy, where n (φ) are the positive energy eigenvalues of H BdG 1 . For a weakly transparent junction, such as the tunnel oxide junctions used in Al-based superconducting devices, the energy E GS (φ) is well-approximated by the form −E J cos φ. In this case, one recovers the canonical superconducting qubit Hamiltonian, H = E c (N − N g ) 2 − E J cos(φ). In the "transmon" limit E J E c , the low-energy excitations of the junction are quantized charge oscillations-due to Cooper-pair tunneling across the junction-with a characteristic plasma frequency, and (crucially for qubit applications) a slightly anharmonic spectrum. If, on the other hand, H BdG has low-energy quasiparticle excitations with energies n ω p , this description is no longer valid. Any correct description must include both the bosonic and the fermionic low-energy degrees of freedom present in the Hamiltonian of Eq. (1). Low-energy fermionic excitations naturally appear if the system is driven into a topological superconducting phase, where zero-energy Majorana modes emerge at the spatial boundaries between trivial and topological regions in the system. Moreover, as the system undergoes the topological phase transition between the conventional and the topological superconducting phase, the gap in the bulk of the system closes, giving rise to a continuum of states at energies below the plasma frequency. In general, one may expect sub-gap quasiparticles to arise quite generically in nanowire-based Josephson junctions as well as in junctions based on other systems engineered to support topological superconductivity, such as proximitized surfaces of topological insulators [77]. In these systems, many competing effects are involved, such as large magnetic fields, spin-orbit coupling, and interfaces between materials with very different properties, which can lead to complicated junction spectra even when the system is not tuned to the topological phase.
In a realistic model of a nanowire, the number of fermionic degrees of freedom appearing in H BdG (φ) may be too large to treat exactly. However, we expect that the effect of highenergy quasiparticles at energies n ω p can be captured by their contribution to the phase-dependent ground state energy. Assuming the latter can be approximated by a cosine dispersion, we rewrite the Hamiltonian in Eq. (1) as where now the fermionic degrees of freedom c correspond to the low-energy degrees of freedom, and E 0 J accounts for the high-energy degrees of freedom.

Experimental probes and parameters
The dynamics of the junction can be experimentally probed in a circuit quantum electrodynamics (cQED) setup [6], in which the junction is coupled to a microwave resonator cavity. The system is driven by an AC field, which corresponds to a time-dependent voltage V g in Fig. 1, V g (t) = V g + δV g (t). This time-dependent voltage couples to the total charge operator on the left island, N tot = N + N f,L , leading to a small time-dependent contribution to the Hamiltonian, δH(t) = E c N tot (C/e) δV g (t). We will characterize the response of the system to this perturbation by the spectral function of the charge operator, Here, |α denotes an eigenstate of the system with energy E α , with α = 0 the ground state. The spectral function exhibits a peak at each transition energy of the system, with the intensity of the peak related to the matrix element of the total charge operator between the initial and final states of the transition. We now discuss interesting parameter regimes for typical cQED circuits which we consider in our simulations. Transmon qubits typically operate in the regime E J /E c ≈ 20 or higher in order to suppress charge noise. E c typically varies in the 200-500 MHz range, yielding plasma frequencies in the 5-10 GHz range. Gate-controlled nanowire junctions allow for a tunable E J and thus allow a device to be operated not only in the transmon regime but also in a regime with lower E J /E c . The lower the ratio E J /E c , the larger the charge dispersion, i.e. the dependence of the energy levels on the dimensionless charge N g . An intermediate ratio E J /E c ≈ 5 is interesting since, as shown in Refs. [60,61], the presence of Majorana zero modes coupled across the junction is associated with distinct features in the charge dispersion.

Gauge transformation
To study the model (1), it is useful to perform a unitary gauge transformation, H → U HU † , with [78] U = e iφ N f,L /2 .
Under this gauge transformation, the number operator N = −2i∂ φ transforms as N → N − N f,L . This simplifies the first term in Eq. (1), which after the transformation is given by E c (N − N g ) 2 . In addition, the fermion operators that belong to the left island transform as c i,σ → e −iφ/2 c i,σ , while the fermions on the right island are unaffected by the transformation. Hence, after the transformation, the operators c i,σ with i ∈ I L are charge-neutral, while the operator N counts the total charge of the superconducting island and of the left part of the wire. The effect of the transformation on the terms in H BdG (φ) is the following. The pairing terms on the left part of the wire, which initially take the form ∆e iφ c i,↑ c i,↓ , lose their phase dependence and are given by ∆c i,↑ c i,↓ . Meanwhile, terms describing hopping between the left and the right parts of the wire acquire phase dependence: c † i,σ c j,σ with i ∈ I L and j ∈ I R becomes e iφ/2 c † i,σ c j,σ . We must also discuss the effect of the gauge transformation on the wave functions. A complete basis for the Hilbert space of the model (1) is given by |n, n L , n R , where n L(R) denotes the vector of occupation numbers for the fermionic degrees of freedom on the left (right) part of the wire, and n ∈ Z is the number of Cooper pairs (counted from charge-neutrality) in the left superconductor, i.e. N |n, n L , n R = 2n |n, n L , n R . A generic many-body state is thus given (in the original gauge) as Although we will eventually use this "number basis" in the numerics, it is convenient to temporarily work in the "phase basis", formed by the states |φ, n L , n R = (2π) −1 n e iφn |n, n L , n R .
The wave function coefficients in the phase basis are Ψ(φ, n L , n R ) = n e iφn Ψ(n, n L , n R ).
They are 2π-periodic: Ψ(φ + 2π, n L , n R ) = Ψ(φ, n L , n R ). The action of the gauge operator U on a phase basis state is simply where n L is the number of electron in the left island. Thus, the gauge transformation maps the state |Ψ to a new state, |Ψ = U |Ψ , with wavefunction coefficients Because of the extra phase factor, the boundary conditions forΨ are periodic or antiperiodic depending on the parity of the number of fermions in the left island [79]: With this new boundary condition, the spectrum of the operator N = −2i∂ φ changes from 2Z (in the original gauge) to Z (in the new gauge). This is in agreement with the fact that, as mentioned above, N must now account for the total charge of the left island and not only for the part due to the paired electrons. However, this accounting is only consistent if the parity of N , which we refer to as the bosonic parity and denote by P b = e iπN , is the same as the parity of the number of (now charge-neutral) fermions on the left island, Figure 2: Minimal model describing a nanowire in the topological phase, with Majorana zero modes, γ i=1,..,4 , at the ends of the topological superconducting regions.
That is, in the new gauge every physical state must obey the following constraint: Note that, since P b = e iπN is the operator that translates φ by 2π, the constraint (16) is simply a rewriting of (15) in a basis-independent form. Finally, we comment on the effect of the gauge transformation on the calculation of the spectral function (8). Since, after the gauge transformation, N is the total charge on the left island, the correlation function that appears in the integral in Eq. (8) is now simply N (t)N (0) .

Minimal model for a nanowire in the topological superconducting phase
We now turn to a minimal realization of Eq. (7), namely the case in which there is exactly one low-energy fermionic mode on each side of the junction. The setup is chosen to capture the essential ingredients of a nanowire in the topological phase, hosting one Majorana mode at each end of each topological segment. We denote the Majorana modes at the ends of the left (right) island by γ 1,2 (γ 3,4 ) as depicted in Fig. 2. Assuming the absence of additional low-energy fermionic quasiparticles, the effective Hamiltonian, written in the gauge where the boundary conditions of Eq. (15) hold, is given by: with The term H δ couples the Majorana modes within each island. Such a coupling can arise due to the finite length of each island, and it is chosen to be identical on both islands for simplicity. In the topological phase, δ vanishes exponentially with the length of each island. The phase-dependent coupling E M cos(φ/2) originates from single-electron tunneling across the junction. This model was introduced and studied by Ginossar and Grosfeld [60].
In the present work we discuss it in the context of the more general model (1), and address additional points that were not discussed in detail in Ref. [60]. To solve the model (17), we first deal with H J and H δ separately, and then consider the effect of H M , which couples the bosonic and fermionic excitations. The eigenfunctions of H J are known exactly in terms of Mathieu functions (approximate eigenfunctions are also immediate to find numerically). Anticipating the role of the boundary conditions (15), we will find eigenfunctions in the space of 4π-periodic functions. The 2π-periodicity of H J implies that [H J , P b ] = 0, and hence we can choose eigenstates to have well-defined bosonic parity. We denote the m-th energy eigenstate with even/odd bosonic parity by |m, ± b , and the corresponding energy by E m,± . The wave functions ψ m,± (φ) = φ|m, ± b are given by where me ν (z, q) is the Mathieu function with characteristic exponent ν, and where The wave functions satisfy the boundary condition ψ m,± (φ+2π) = ±ψ m,± (φ). In the limit E J E c , they reduce to normalized plane waves, ψ m,± (φ) ∼ e ik m,± φ / √ 4π, with k m,+ ∈ Z and 2k m,− ∈ Z. In the opposite limit E J E c , the wavefunctions are localized near the minima of the potential −E J cos φ, i.e. near φ = 0 and 2π.
To describe the fermionic sector, we observe that two Majorana modes together form a fermionic mode, and define its occupation as n ij = (1 + iγ i γ j )/2. A basis for the fermionic Hilbert space can thus be obtained by arranging the Majorana modes into pairs and specifying the occupations of the pairs. We focus on the sector with even total fermion parity, and take the pairs of Majorana modes to be those in the same superconducting region, such that the basis states also have well-defined fermion parity in each island. Then, the two possible states are |n 12 = n 34 = 0 and |n 12 = n 34 = 1 .
We now define a basis for the full Hilbert space describing both the bosonic and the fermionic sectors, using states which obey the constraint in Eq. (16): These states have equal boson parity and fermion parity on each island; for the remainder of this section, we will refer to this as island parity. The Hamiltonian H J + H δ is diagonal in this basis; the state |m, p has energy E m,p − 2pδ, where p = ±1 is the island parity. The four lowest-energy states in the regime with δ ω p (m = 0, 1 and p = ±) are shown schematically in Fig. 3(a), illustrating the parity constraint.
The coupling term H M is entirely off-diagonal in this basis, since it couples states of opposite island parity. Its nonzero matrix elements are In the limit E J E c , an asymptotic analysis of the integral in (23) shows that the dominant matrix elements are the diagonal ones, η m,m = 1−O(z), where z = 2E c /E J 1. The matrix elements in which m and m differ by an even integer 2 are subdominant, η m− ,m+ = O(z ). The matrix elements in which m and m differ by an odd integer, η m,m+2 +1 , are exponentially small in the ratio E J /E c . They vanish identically when N g is an integer and are largest when N g is half-integer.
The low-energy spectrum as a function of N g is shown in Fig. 3(b) in the two limits of large δ ω p (left panel) and δ = 0 (right panel). The colored dashed lines correspond to the spectrum for E M = 0, while the solid black lines correspond to the spectrum in the presence of a finite E M . The dispersion of energy with N g can be understood in (c) Excitation spectrum as a function of δ/ω p for N g = 0, and the spectral function at δ = 5ω p (blue) and δ = 0 (red). The spectral function is convolved with a Gaussian with σ/ω p = 10 −2 .
terms of instanton tunneling processes between the two minima of the potential at φ = 0 and 2π (quantum phase slips), and has a magnitude proportional to e − √ 8E J /Ec . The dispersion of levels with opposite island parity is shifted by one unit along N g , leading to level crossings at half-integer values of N g at δ = 0, see right panel of Fig. 3b. H M introduces a coupling proportional to E M η m,m (N g ) between states of oppposite parity, leading to avoided crossings at the degeneracy points.
We now study the behavior of the excitation spectrum for a choice of parameters that mimics driving the system from a trivial and fully gapped superconducting phase into a topological phase with well-separated Majorana zero modes. Fixing N g = 0 and E J /E c = 5, we vary the coupling between the Majorana modes on the same island from a large value δ ω p to a small value δ ω p , and compare the excitation spectrum with E M = 0 to the one with a small but finite E M ω p . The resulting spectra are shown in Fig. 3(c).
For large δ ω p , the ground state is |0, + and all the states with odd island parity, |m, − , are at high energy. The spectrum of the junction resembles that of a trivial Josephson junction, with level spacing set by ω p , up to anharmonic corrections. When δ ≈ ω p /4, there are degeneracies in the spectrum between states |m, − and |m + 1, + . At finite E M the coupling between these states is proportional to η m+1,m , and hence vanishes for N g = 0. For δ ω p , the two states in each doublet, |m, ± , get closer in energy. In the limit δ = 0 the energy splitting between these states is set by the larger of the two energy scales E M and the splitting due to charge dispersion, |E m,+ − E m,− |. Note that for the doublets with odd m, one has E m,− < E m,+ at N g = 0, and thus the state with odd island parity crosses the one with even parity in energy as δ is decreased. In Fig. 3(c) this can be seen as the crossing of the lines corresponding to |1, + (red) and |1, − (green). At finite E M , the coupling between these states gives rise to an avoided crossing between them, with size of order η m,m E M .
Using the excitation spectrum, we calculate the spectral function (8) for δ ω p and δ = 0 (see right panel of Fig. 3(c)). When δ ω p , the ground state is simply |0, + . Since the operator N can only couple states with the same parity, only the spectral lines within the even parity sector are observed in the spectral response. Once different parity states couple due to finite E M , additional spectral lines can be observed. In particular, for δ = 0, both spectral lines in the m = 1 manifold can be observed. Deep in the harmonic limit, when |E m,+ − E m,− | → 0, the eigenstates are superpositions of states with well-defined fermionic parity of the junction (i.e. occupation of n 23 ), and due to the total fermionic parity constraint also of n 14 . As the charge operator N acts locally at the junction, it cannot change the occupation of n 14 , and only a single spectral line is visible again. For a further discussion of this limit, see Sec. 5.2.

Numerical simulations
Building on the picture developed in the previous sections, we now turn to a numerical study, using DMRG and TEBD, of the spectral function defined in Eq. (8). We will perform the study on a microscopic model that allows us to treat the behavior both in the topological phase and in the vicinity of the topological phase transition, as well as to examine the effect of additional Andreev subgap states in the junction both in the topological and non-topological regime.
Numerical simulations presented in this work were performed using the ITensor library [80].

Lattice model and numerical techniques
In our numerical simulations we use a lattice description of the system. To this end, we introduce the tight-binding parameters corresponding to the hopping amplitude, t = 1/(2ma 2 ), and the spin-orbit coupling, v = α/a, where a is the lattice constant (see Appendix A for the explicit tight-binding Hamiltonian). Unless otherwise specificed, all energies hereafter are measured in units of ∆. We describe the system in the gauge introduced in Sec. 2.2, where the Hilbert space consists of a lattice of fermionic degrees of freedom and a bosonic degree of freedom describing the total charge on the floating island. We represent this bosonic mode in the charge basis and as a single additional site. The occupation of this charge mode depends on the ratio E J /E c , but the number fluctuations grow slowly, δN 2 ∼ (E J /E c ) 1/2 for E J > E c [4], so in practice we find that truncating to 10−20 charge states yields sufficiently small error. The charging energy term acts only on the bosonic charge state, which is coupled to the fermionic modes via the hopping terms between the left and the right parts of the wire, i.e. via e iφ/2 c † i,σ c j,σ with i ∈ I L , j ∈ I R (recall that in the charge basis the operator e iφ/2 simply changes the charge by one). Hence, if the charge site is placed at the boundary between the left and the right part of the wire, the Hamiltonian is local, simplifying the numerical study. The constraint introduced by the boundary conditions (15) is handled by defining a conserved Z 2 quantum number P b P f,L = 1.
To obtain the spectral function S N (ω) as defined in Eq. (8) we first obtain the ground state |0 of the system using DMRG. We then perform time evolution of the state |0 ≡ N |0 in order to obtain the real-time correlation function 0|N (t)N (0)|0 = e −iE 0 t 0 |e iHt |0 . For the time evolution, we use TEBD with a 4th order Suzuki-Trotter decomposition with a time step of dt = 0.02, truncating the MPS wavefunction to bond dimension of M max = 50 and truncation error tr = 10 −7 . We note that similar numerical techniques could be used to characterize the correlation function for some other initial state, such as a low-lying excited state or even a mixed state. This may be of interest to describe experimental situations where finite temperature or non-equilibrium effects lead to a finite occupation of excited states. The computational cost of this time evolution is heavily dominated by the terms of the Hamiltonian that involve the bosonic degree of freedom, which has a large on-site Hilbert space dimension and thus limits the bond dimension we can treat. To obtain the spectral function in real frequencies, we typically compute the real-time correlation function up to times t f 400 and use linear prediction methods [83][84][85] to extrapolate the real-time correlation function to times ∼ 2t f . We then apply a Gaussian In the trivial phase, for B/∆ B c /∆ = 1, a single spectral line is seen, corresponding to the plasma frequency. As the system enters the topological phase, a second spectral line appears, and an avoided crossing between the two can be observed (see main text for further discussion). The red dashed lines correspond, for B < B c , to the spectral lines expected for a conventional Josephson junction, and for B > B c , to those expected for a junction deep in the topological phase. The latter are obtained using the minimal model given by the Hamiltonian (17) (see main text for more details).
windowing function and finally perform the Fourier transformation. This allows us to reach a frequency resolution of order 10 −2 ω p in the parameter regime we consider.
We consider two different models for the junction, corresponding to the limits of a short and a finite-length junction, that will be described in detail below.

Short junction (weak-link model)
We start from the short junction limit, valid for junctions much shorter than the superconducting coherence length, in which the junction can be modeled as a single weak link between the islands. All hopping terms on this link are reduced by a factor κ < 1 compared to their values in the bulk. The transmission of the junction, and hence E M , is determined by κ; for small κ, E M is proportional to κ. This setup is shown schematically in Fig. 4(a). The exact tight-binding Hamiltonian is given in Appendix A.
In Fig. 4(b), we plot the spectral function obtained for this model, as a function of the Zeeman energy B, for N g = 0 and E 0 J /E c = 5. In the trivial phase, for B B c , a single spectral line is present (in the frequency range shown) at a frequency ω ≈ ω p , as expected. As the Zeeman energy is increased and the wire is driven into the topological phase, a second spectral line appears due to finite E M as discussed in Sec. 3. The avoided crossing between the two spectral lines in Fig. 4(b), which can be observed close to the topological phase transition at B c /∆ = 1, is reminiscent of the avoided crossing between the states |1, + and |1, − in Fig. 3(b).
To validate our numerical results, we overlay on the spectral function the spectral lines expected deep in the trivial and the topological phase, obtained using the minimal model discussed in Sec. 3. These are plotted as red dashed lines in Fig. 4(b). For the trivial phase (B < B c ), the position of the spectral line is calculated from H J of Eq. (18a) with Josephson energy set to E 0 J . For the topological phase (B > B c ), we use the the Hamiltonian (17) with δ = 0 and a value of E M determined numerically from the 4π-periodic component of the ground state energy deep in the topological phase (more specifically, at B/∆ = 2), as explained in detail in Appendix B. As can be seen, our numerical results indeed agree very well with these values in both limits.

Topological phase transition
The finite-size gap at the transition is determined by the spin-orbit coupling strength as ∼ α∆/(BL) (see Ref. [86] for details). For the parameters used to obtain Fig. 4, we have ∼ ω p . To probe the response of the system closer to the continuum limit we consider a smaller spin-orbit coupling, such that many states cross the plasma frequency close to the topological phase transition, but still large enough to have a sizable gap in the topological region. However, we still do not observe any significant features in the spectral response associated with the gap closing, as can be seen in Fig. 5, where we also plot the energies of the two-quasiparticle excitations on top of the spectral function. We attribute this to the fact that the critical states have weak dependence on the phase difference at the junction, and thus small matrix elements of the charge operator N which determines the spectral response via Eq. (8). An intuitive picture for this observation is that the critical states

Finite-length junction
We now consider a finite-length junction, which is modeled as a finite segment of length W of normal (non-superconducting) wire. The hopping and the spin-orbit coupling between the left (right) superconducting region of the wire and the junction is reduced by a factor κ L(R) compared to their values in the bulk. The exact tight-binding description of the model is given in Appendix A and is shown schematically in Fig. 6(a).
In Fig. 6(b) we plot the spectral function obtained for this model. We find that in this case the structure of the spectral function is more complicated with additional spectral lines appearing both in the trivial and in the topological phase. To understand this spectral function, we first obtain the spectrum expected on the basis of the minimal model of Sec. 3 deep in the trivial and topological phase, as we did for the short junction case (see previous sub-section for more details). The red dashed lines plotted on top of the spectral function for B < B c (B > B c ) correspond to this spectrum. In addition, we plot as white solid lines the energies of two-quasiparticle excitations, obtained from the tight-binding BdG Hamiltonian at a phase difference of φ = 0 between the superconducting islands, as a function of B. These are the fermionic excitations (in the even parity sector) originating from the second term in the Hamiltonian (1), neglecting the dynamics of the field φ and the finite charging energy. It can be clearly seen that the lines in the spectral response originate from avoided crossings between these two types of excitations. In Sec. 5 below, we will present a perturbative approach that will allow us to understand this coupling and to calculate the avoided crossings in the harmonic limit, i.e. for E J /E c → ∞.

Effective theory for the coupling between bosonic and fermionic modes
We now examine the observed avoided crossing betweens the bosonic plasma mode and fermionic subgap states in more detail. To this end, in Sec. 5.1 we develop a perturbative theory in the harmonic limit. In Sec. 5.2 we discuss non-perturbative couplings that cannot be captured by the perturbative expansion, and show their effect on the spectrum of the problem using a simple model.

Harmonic expansion
We start from the Hamiltonian (7), where we assumed that the effect of the high-energy quasiparticles is captured by their contribution to the phase dependent ground state energy that can be approximated by a cosine dispersion −E 0 J cos(φ). In the limit E 0 J E c , phase fluctuations are small and centered around φ = 0. Therefore, we can ignore the fact that φ is a compact variable, and thus also the periodic or anti-periodic boundary conditions (15). In this approximation, the N g -dependence of the Hamiltonian can be "gauged away" by generalizing the gauge transformation (9) to U = e iφ(N f,L −Ng)/2 and the dependence of the eigenvalues on N g is lost.
Expanding the cosine dispersion to lowest order in φ, we denote by H 0,b the resulting harmonic oscillator Hamiltonian acting on the bosonic degree of freedom, Here ω 0 p = 8E 0 J E c , and a † and a are harmonic oscillator raising and lowering operators, respectively. We will denote the eigenstates of H 0,b as |m , where m is the occupation level of the harmonic oscillator.
We now expand H BdG (φ) around φ = 0, We denote the zeroth order term in this expansion by H 0,f and diagonalize it to obtain Here, the operators Γ † i (Γ i ) create (annihilate) quasiparticle excitations with non-negative energies i . We define a corresponding Nambu spinor Γ † = (Γ † i , Γ i ) and denote the singleparticle unitary that diagonalizes H BdG (0) by V , i.e. c = V Γ. We will denote the  Figure 7: (a) Excitation spectrum and (b) spectral function obtained using the harmonic expansion for a finite-length junction with the same parameters as in Fig. 6. The black dashed lines in (a) correspond to the excitation energies of the unperturbed Hamiltonian, while the blue solid lines are obtained including the perturbation to lowest order in φ. The truncated basis used here consists of two lowest energy harmonic oscillator states and two-quasiparticle fermionic states constructed from four lowest energy single-particle states. For presentation purposes we plot the spectral function in (b) convolved with a Gaussian with σ = 5 × 10 −3 .
eigenstates of H 0,f as | ν , where ν = (ν i ), with ν i ∈ {0, 1}, is the vector of occupations of the quasiparticle levels. In particular, the ground state will be denoted as | 0 . The excited states are then given explicitly as | ν = i (Γ † i ) ν i | 0 . We take the sum H 0 = H 0,b + H 0,f as the unperturbed Hamiltonian for the problem. Note that since H 0,b(f ) acts solely on the bosonic (fermionic) degrees of freedom, the eigenstates of H 0 are simply the tensor products |m, ν ≡ |m ⊗ | ν . Their unperturbed energies are given by The second term in the expansion (26) acts as a perturbation to H 0 , which we denote as δH. This term introduces a coupling between the bosonic and fermionic degrees of freedom. Noting that the phase φ has the representation φ = √ z(a + a † ), where z = (2E c /E 0 J ) 1/2 is a small parameter, and using the basis | ν for the fermionic states, δH can be written as The matrix elements of H BdG (0) are related to the dispersion of the quasiparticle levels with the phase across the junction, and hence to the current carried by these levels.
The perturbation δH introduces a coupling between the states |m, ν and |m , ν when m − m = ±1, with a matrix element proportional to ν|Γ † [V † H BdG (0)V ]Γ| ν . The problem can now be easily tackled numerically, solving the Hamiltonian H 0 + δH with a truncated basis consisting of low-energy many-body states. The spectral function S N (ω) of Eq. (8) can also be computed numerically using that, in the harmonic limit, N = i(a † − a)/ √ z. The spectrum and the spectral function calculated for the finite-length junction model (see Sec. 4.3) using this approach are plotted in Fig. 7, for the same model parameters as in Fig. 6. Comparing Figs. 6 and 7, we find that many of the qualitative features appearing in the spectral function are reproduced within the harmonic expansion. In particular, the characteristic behavior of the spectral function near avoided level crossings can be easily understood. For concreteness, consider a crossing between the unperturbed levels |1, 0 and |0, ν as the magnetic field B is tuned through some value B 0 . Near the crossing, we may approximate the unperturbed energies as E 1, 0 (B) ≈ ω 0 p and E 0, ν (B) ≈ ω 0 p + 2λ(B − B 0 ), where 2λ ≡ E 0, ν (B 0 ). Including the perturbation δH and solving the resulting two-level problem, we obtain hybridized eigenstates of the form |± = a ± (B)|1, 0 + b ± (B)|0, ν with energies The brightness of the spectral line corresponding to the transition |0 → |± , where |0 denotes the ground state, is proportional to | ±|N |0 | 2 . This matrix element is just |a ± (B)| 2 /z: The intensity of the two spectral lines is thus identical at the degeneracy point, B = B 0 , and becomes more and more asymmetrical away from the degeneracy point. At the same time, since Fig. 6 was obtained for E 0 J /E c = 5, i.e. not very deep in the harmonic limit, some features are not captured correctly. First, in Fig. 7 the plasma frequency is higher than Fig. 6. This discrepancy, which also shifts the exact positions of some of the avoided crossings, can be explained by the anharmonic correction to the plasma frequency that is brought by the fourth-order expansion of −E 0 J cos φ, δω 0 p ≈ −E c [4], and which is automatically included in the MPS simulations of Sec. 4. Second, we note the level crossing at B/∆ ≈ 2.7 and ω/∆ ≈ 0.6 in Fig. 7, which is avoided in Fig. 6 (near B/∆ ≈ 2.3 and ω/∆ ≈ 0.5). The fact that this crossing is not avoided within the harmonic expansion can be understood as follows. Since φ appears in H BdG (φ) only in the phases of the terms that hop fermions between the left and right parts of the wire, the perturbative couplings in δH only involve fermion quasiparticle levels with support at the junction and are thus local. On the other hand, we find that the two states involved in the crossing differ in the occupation of fermion modes far away from the junction and thus cannot be coupled in the harmonic expansion. In Sec. 5.2, we will discuss how the avoided crossing can be understood in terms of non-perturbative effects. To show that both discrepancies are due to the relatively low value of E 0 J /E c , in Appendix C we compare the spectral function obtained using the exact MPS simulation and the one obtained using the harmonic expansion for a higher value E 0 J /E c = 20, and show that there is an excellent agreement between the two.
Finally, we note that in principle it should be possible to systematically improve the perturbative expansion in z by including the higher-order terms which were so far neglected in Eq. (26) as well as in the expansion of the cosine potential. In general, the n-th order term introduces a coupling between the states |m, ν and |m , ν with |m − m| ≤ n. In addition, it can cause a shift in the eigenvalues, even away from level crossings. However, in order to be consistent, the expansion would have to take into account the contribution of the low-energy fermionic degrees of freedom to the plasma frequency as well as the coupling between the low-and high-energy degrees of freedom which were separated in Eq. (7); thus, it appears to be a challenging approach.

Non-perturbative couplings
We now address the non-perturbative couplings that are not captured within the harmonic expansion. As mentioned in passing in Sec. 3, non-perturbative corrections arise due to quantum phase slips between the minima of the potential energy at φ = 0 and φ = 2π. In the case of the standard qubit Hamiltonian Figure 8: Spectrum of the Hamiltonian (32) projected onto the low-energy subspace, for E J /E c = 5, E M /E c = 1. We take (B) =ω + λ(B/B 0 − 1) with λ/E c = 3 and the coupling ζ = 0.1. Black dashed lines correspond to the spectrum ofH, while the blue solid lines correspond to the spectrum ofH +H ζ . In (a) we artificially set the magnitude of the charge dispersion δ 1 to zero, thus neglecting all non-perturbative corrections, while in (b) the spectrum with finite δ 1 is plotted. Once non-perturbative corrections are taken into account, hybridization between |e + and |e − (due to finite charge dispersion) results in the avoided crossing between the hybridized states |ẽ ± and |f ± . that, in the limit E J E c , quantum phase slips give rise to the charge dispersion of the energy levels with magnitude ∝ e − √ 8E J /Ec [4]. This effect is, for instance, visible in the energy spectra of Fig. 3(b). These corrections are non-perturbative, as manifested by their singular dependence ∝ e −c/z on the expansion parameter z of the previous Section.
In this Section, we examine in detail an example for how these non-perturbative corrections can prominently affect the microwave response in the topological phase. In particular, we show that non-perturbative effects can resolve level crossings that are not resolved perturbatively, leading to the appearance of avoided crossings in the spectral response. The basic physics of this phenomenon is as follows: consider a level crossing between two given states |e + and |f + , as in Fig. 8(a), which are not coupled by perturbative terms local at the junction. If non-perturbative effects cause |e + to hybridize with a third state |e − that does couple perturbatively to |f + , this will lead to a non-perturbative avoided crossing, as seen in Fig. 8(b).
As a concrete example, we extend the minimal model presented in Sec. 3 to include an additional fermionic mode localized (for simplicity) on the right island. We assume that this mode couples to the Majorana mode on the left island via single-particle tunneling across the junction. The Hamiltonian is then given by where and we have implicitly set N g = 0. Here, c † (c) is the creation (annihilation) operator for the extra fermionic mode, is its energy, and ζ is its coupling strength to the Majorana mode on the left island. We assume that the energy is comparable to the plasma frequency ω p , and that it can be tuned by an external parameter (such as the Zeeman energy B, as in the previous section). Writing the additional fermion operator in terms of Majorana operators, c = (γ 5 + iγ 6 )/2, we rewrite H ζ as Analogously to the basis defined in Eq. (21) that was used to study the four-Majorana model in Sec. 3, we can define a basis for the full Hilbert space of this model that satisfies the parity constraint (16) and diagonalizes H J + H : |m; n 12 , n 34 , n c = |m, (−1) n 12 b ⊗ |n 12 , n 34 , n c .
Here, |m, ± b are the bosonic states defined in Sec. 3 and correspond to the m-th energy eigenstate of H J with even/odd bosonic parity and energy E m,± ; n 12 and n 34 are the occupation numbers encoded in the MZMs of the left and right island, n ij = (1 + iγ i γ j )/2; and n c is the occupation number of the extra fermionic mode.
In the harmonic limit E J /E c → ∞, δ m ≡ E m,+ − E m,− → 0 and thus, for each m, the states with different n 12 but equal n c become degenerate eigenstates of H J + H . Furthermore, in this limit, H M has non-vanishing matrix elements only between states with the same m. Restricting our discussion to low energies (namely, to states with energies up to order ω p ), in the harmonic limit the eigenstates of H J + H + H M are given by: At a large but finite E J /E c , the energy levels of H J acquire a finite charge dispersion, In the harmonic limitω → ω p . We see that |g ± and |f ± are still eigenstates ofH, but H e mixes |e ± ; its eigenstates, to lowest order in q ≡ δ 1 /(2η 1,1 E M ), are given by At this point, we can understand why, as pointed out in Sec. 3, the two spectral lines in the topological phase are only visible away from the harmonic limit. To this end, note that the charge operator N couples |g + to |e + but not to |e − ; on the other hand, it couples |g + to both |ẽ + and |ẽ − .
Finally, we consider a finite coupling ζ. Projecting H ζ onto the basis (36), we obtaiñ where, similarly to Eq. (23), χ m,m is the overlap integral with ψ m,± (φ) ≡ φ|m, ± b . Note thatH ζ only couples |e + to |f − and |e − to |f + . Thus, in the harmonic limit, or more specifically for vanishing δ 1 , there is no avoided crossing between |e + and |f + , as can be seen in Fig. 8(a). Away from the harmonic limit, the magnitude of the avoided crossing between |ẽ + and |f + is This avoided crossing is clearly visible in Fig. 8(b). It is manifestly non-perturbative, and vanishes in the harmonic limit.

Conclusions and Outlook
We have carefully studied the response of a capacitively shunted topological Josephson junction, using a combination of accurate numerical techniques and theoretical approaches that allow us to incorporate a microscopic description of the topological phase. Our results indicate that detecting the signatures of the topological phase in the transmon-like setup of Fig. 1, as proposed for instance in Ref. [60] and also as described in Sec. 3, can be complicated by two factors. First, we do not find strong signatures of the topological transition itself on the simulated frequency spectra of the junction. Second, the presence of non-topological Andreev bound states can hinder the signatures of coupled Majorana zero-modes at the Josephson junction in the measurable low-frequency spectrum. Our simulations show that this second problem becomes particularly relevant in parameter regimes that tend to have a higher density of sub-gap states that couple to the phase degree of freedom. For a quantitative understanding of experiments performed on nanowire Josephson junctions, including the electrostatic environment and the presence of multiple sub-bands can be very important. While the numerical calculations presented here are based on an effective one-band model of a Majorana nanowire, the methods themselves can in principle be applied to more realistic models. This is particularly true for the perturbative approach described in Sec. 5, which can take as input low-energy BdG spectra obtained from largescale microscopic simulations of realistic devices. For future research, it would also be valuable to find a way to systematically include the contribution of quantum phase slips without resorting to an intensive MPS calculation, and to consider the case in which the nanowire Josephson junction hosts a quantum dot with finite charging energy.  Figure 9: Energies of the ground state (E 0 ; blue dashed line), and the first excited state (E 1 ; green dashed line), obtained using the BdG spectrum, calculated for the weak-link tight binding model (see Appendix A). Tight binding parameters used are the same as in Fig. 4, with a Zeeman energy B = 2B c = 2. The approximate crossings in the spectrum at φ = π, 3π allow us to define a 4π-periodic ground state energy E 4π , plotted as the solid red line.
For the finite-length junction model studied in 4.3, the Hamiltonian for the junction is In order to extract an effective E M from the tight binding model of the junction, we diagonalize the BdG Hamiltonian and calculate the energies of the ground state and first excited state as functions of the phase φ defined on a 4π-periodic domain (see Fig. 9). In a Josephson junction geometry and for an infinite wire, the model of Eq. (24) would give an Andreev level crossing of at φ = π + 2πn (where n ∈ Z), leading to the 4π-periodicity of the ground state. For a finite system, the crossings will be avoided, with the splitting determined by the coupling between the Majorana modes at the junction and those at the far ends of the wire. Nevertheless, If the system is deep in the topological phase (in practice we consider a Zeeman energy of magnitude B = 2B c ), the avoided crossings at φ = π, 3π will be small enough to allow us to define a 4π-periodic ground state energy, E 4π (φ). An effective E M is then defined as . C Numerical results in the harmonic limit In Fig. 10 we plot the spectral function obtained from the MPS simulations in the harmonic limit, E 0 J /E c = 20, and the spectral function obtained using the harmonic expansion for the same parameters. To highlight the agreement between the two for the strength of the couplings between the plasma mode and the fermionic quasiparticle excitations, we take into account the anharmonic corrections that shift the energy of the first excited state with respect to the plasma frequency in the harmonic limit ω p = √ 8E J E c , by δω p ≈ −E c . This is achieved by including the next order term in the expansion of the cosine potential, namely −E J φ 4 /24 = −E J z 2 a + a † /24, as part of the bosonic Hamiltonian H 0,b , given in Eq. (25).