Topological Josephson Bifurcation Amplifier: Semiclassical theory

We describe the topological Josephson bifurcation amplifier using the resistively and capacitively shunted (RCSJ) model of a Josephson junctionlikharev1986 ; blackburn2016 . In this model (Fig. 1), the Kirchhoff law for the electric current traversing the junction reads

where $\delta$ is the gauge-invariant superconducting phase difference between the two superconducting electrodes, $C$ and $R$ are the effective capacitance of the junction and the characteristic impedance of the microwave source generator, $\varphi_{0}=\hbar/2e$ is the reduced flux quantum, $I_{d}$ is the bias current and $U$ is the grand potential of the junction at temperature $T$. The first, second and third terms in the left hand side of Eq. (1) correspond to the displacement current, the dissipative ohmic current and the dissipationless (Josephson and Majorana) current, respectively. The critical current of the junction is given by $\varphi_{0}^{-1}{\rm max}(\partial U/\partial\delta)$.

In the RCSJ model, $\delta$ is treated as a classical variable, whereas $U(\delta)$ is computed by diagonalizing a quantum mechanical Hamiltonian for each value of $\delta$. Charging energy is not included in the Hamiltonian; instead, the phase fluctuations produced by the charging energy are incorporated via the displacement current in Eq. (1). This semiclassical approach is justified in the "transmon regime" that is relevant to the operation of the JBAvijay . In such regime, the Josephson energy of the junction (approximately given by $(U(\pi)-U(0))/2$) greatly exceeds the charging energy $E_{C}=e^{2}/(2C)$.

In thermodynamic equilibrium, fermion parity constraints are absent for timescales that are long compared to quasiparticle poisoning times. Then, the grand potentialbeenakker2013 takes the form

where $\epsilon_{n}(\delta)$ are the single-particle energies of the system, with $n=\pm 1,\pm 2,...$ and $\epsilon_{n}=-\epsilon_{-n}$ due to particle-hole symmetry (by convention, we take $\epsilon_{n}>0$ for $n>0$). For $T\rightarrow 0$, Eq. (2) reduces to the ground state energy

The eigenvalues $\epsilon_{n}(\delta)$ can be obtained by diagonalizing a tight-binding Hamiltonian describing a single-channel, one-dimensional superconducting/normal/superconducting heterostructure of finite length, where two long superconducting electrodes are separated by a weak link in the normal state. More details of the model are presented in Sec. II.3.

Equation (1) describes the dynamics of a fictitious particle of mass $C$ and position $\varphi_{0}\delta$ moving in a potential $U(\delta)$ under the action of an external force $I_{d}$. We consider a noisy sinusoidal AC-bias current,

where $\omega_{d}$ and $i_{d}$ are, respectively, the drive frequency and amplitude, and $I_{N}(t)$ is a white-noise current produced by thermal fluctuations. We will focus in the case in which $i_{d}$ is small compared to the critical current of the junction. Then, the dynamics of the particle is akin to that of a nonlinear harmonic oscillator centered in one of the minima of $U(\delta)$, which are located at $\delta_{m}=2\pi m$ ($m\in\mathbb{Z}$). Because $U(\delta)=U(\delta+2\pi)$ in the absence of parity constraints, we can concentrate on oscillations near $\delta=0$ without loss of generality.

Following the small oscillation approximation, Eq. (1) can be rewritten as

where the damping rate $\kappa=1/RC$, the plasma frequency $\omega_{p}$ in the harmonic approximation and the leading anharmonicity parameter $\lambda$ play a central role in the functioning of the JBA. They can be extracted from $U$ via

where $U^{(n)}(0)=(\partial^{n}U/\partial\delta^{n})|_{\delta=0}$.

Figure 2 illustrates its solution in the absence of current noise. The various curves in this figure are well-known in the literaturevijay ; here, we highlight the most important points for completeness and later reference. In a frame rotating at the drive frequency, i.e. taking $\delta(t)=\mathrm{Re}[\delta_{z}e^{-i\omega_{d}t}]$, we solve for the steady-state solution of Eq. (5). Ignoring rapidly oscillating terms (rotating-wave approximation), Fig. 2a displays the modulus of $\delta_{z}$ as a function of $\omega_{d}$, for different values of $i_{d}$. For small $i_{d}$, the response of the oscillator is a Lorentzian peaked at $\omega_{d}=\omega_{p}$. As the driving amplitude increases, the response of the oscillator becomes peaked at frequencies lower than $\omega_{p}$ (this is due to the fact that the anharmonicity parameter in Eq. (6) is negative). Beyond a critical value of $i_{d}$ (but still well below the critical current of the junction), there emerges a finite interval of $\omega_{d}$ for which the response of the oscillator is multi-valued, with three values of $|\delta_{z}|$ for each $\omega_{d}$. In such an interval, the intermediate value of $|\delta_{z}|$ is not stable and the oscillator is thus said to be in a bistable regime, with a low- and a large-amplitude oscillation steady state.

Figure 2b shows the dependence of the bistable regime on the drive frequency. Bistability takes place when $\omega_{d}\leq\omega_{c}\equiv\omega_{p}-\sqrt{3}\kappa/2$ and $i_{d}\in(i_{-},i_{+})$, where $i_{-}$ and $i_{+}$ are the lower and upper bifurcation currents, respectively. These currents obey the relations

where $r(\omega_{d})=(\omega_{d}^{2}-\omega_{p}^{2})^{2}-3\omega_{d}^{2}\kappa^{2}$, $\omega_{\delta}=\omega_{p}-\omega_{d}$ and $\omega_{\Sigma}=\omega_{d}+\omega_{p}$. At the bifurcation currents, the response of the oscillator switches between single-valued and two-valued. Thus, for $i_{d}<i_{-}$ ($i_{d}>i_{+}$), the oscillation amplitude is single-valued and low (high). The fact that bifurcation currents depend on $\omega_{p}$ and $\lambda$ will be important below.

Figure 2c presents the dependence of the steady-state solution $|\delta_{z}|$ on the drive amplitude, for a fixed value of the drive frequency ($\omega_{d}<\omega_{c}$, red dashed vertical line in Fig. 2b). In the bistable regime (blue shaded central region), which of the two stable steady-states is reached depends on the initial state of the system and how the drive is turned on. In the case where the drive is turned on sufficiently slowly (compared to the timescale $1/\kappa$) and in the absence of current noise, the steady state follows an hysteresis curve. For an initial state where $|\delta_{z}|\sim 0$, the low-amplitude steady-state is realized for any $i_{d}<i_{+}$, and the system suddenly jumps to the high-amplitude solution at $i_{d}=i_{+}$. Driving the the Josephson junction close to $i_{+}$ thus allows to detect small changes in the upper bifurcation current, independently from the width of the bistable region $|i_{+}-i_{-}|$.

Thus far we have reviewed the phenomena of bifurcation and bistability in the absence of current noise ($I_{N}=0$). In the presence of noise, due to either thermal or quantum fluctuations, the steady-state reached by the system is not deterministic. We define the bifurcation probability $P$ as the probability of the system reaching the high-amplitude state after evolution of the system for a time $T_{0}$ vijay ; bertet2011 . To estimate $P$, we integrate Eq. (5) for a time $T_{0}>1/\kappa$ and $M\gg 1$ realizations of the noise $I_{N}(t)$. We consider uncorrelated white noise by taking $I_{N}(t)$ from a zero-mean Gaussian distribution of variance $\sigma_{N}^{2}=4k_{B}T/R$. We then estimate the bifurcation probability by $P\approx n_{H}/M$, where $n_{H}$ is the number of high amplitude steady-state solutions. This probability $P$ is experimentally measurable.

Figure 3 illustrates the dependence of $P$ on the drive amplitude in the presence of current noise. When $i_{d}<i_{-}$ and $i_{d}>i_{+}$, the oscillator is with certainty in the low and high amplitude states, respectively. Due to the hysteretic nature of the bifurcation dynamics, $P$ is a sharp step function at $i_{d}=i_{+}$ in the absence of noise (blue curve). Current noise increases the width of the step function in the region $i_{d}\lesssim i_{+}$. As we consider a low amplitude initial state and a smooth turn on of the drive, $P$ does not depend on the lower bifurcation current ($i_{-}\approx 2.25$ nA in Fig. 3). Because the bifurcation currents are sensitive to the intrinsic parameters of the junction (such as $\omega_{p}$ and $\lambda$), small variations of the latter lead to significant changes in $P$ when $i_{d}$ lies in the vicinity of $i_{+}$. As we explain below, this high sensitivity may allow for new ways to measure MBS properties in topological Josephson junctions.

Before turning to a more microscopic description of a Josephson junction by calculating the grand potential $U(\delta)$ from a tight-binding model, it is helpful to consider two approximate limiting forms of $U(\delta)$ in ideal tunnel junctions. These limiting cases will be useful to interpret some of the numerical results of the following sections.

In the case of conventional tunnel junction at low temperature, the ground state energy is $U(\delta)=-E_{J}\cos\delta+{\rm const}$, where $E_{J}$ is the Josephson energy associated to the tunneling of Cooper pairs and "const" is a $\delta-$independent number. From Eq. (6), the plasmon frequency and anharmonicity parameter are then $\hbar\omega_{p}=\sqrt{8E_{J}E_{C}}$ and $\lambda=-\omega_{p}^{2}/6$. In this case, $\lambda/\omega_{p}^{2}$ is a constant because $U(\delta)$ has a single harmonic ($\cos\delta$).

In the case of a topological tunnel junction, one instead finds (disregarding parity constraints) $U(\delta)\simeq-E_{M}|\cos(\delta/2)|-E_{J}\cos\delta+{\rm const}$, where $E_{M}$ is the Majorana energy associated to the tunneling of single electrons. In particular, if the junction transparency is very low, $E_{J}\ll E_{M}$ and $U(\delta)\simeq-E_{M}|\cos(\delta/2)|$. In this tunnel regime where Cooper pair tunneling across the junction is negligible compared to single-particle tunneling, we have $\omega_{p}\sim\sqrt{2E_{M}E_{C}}/\hbar$ and $\lambda\sim-\omega_{p}^{2}/24$. Once again, $\lambda/\omega_{p}^{2}$ is a negative constant, but of smaller amplitude.

In general, when the $\delta-$dependence of $U(\delta)$ has more than one harmonic (e.g. in topological junctions with comparable $E_{J}$ and $E_{M}$, or even in trivial junctions that are neither in the tunneling nor in the perfectly transparent regimegolubov2004 ), $\lambda/\omega_{p}^{2}$ will not be a constant.

Finally, for later reference, we note that in the literature of superconducting qubitsblais2020 , anharmonicity in the transmon regime is often defined as $\lambda^{\prime}=(E_{2}-E_{1})-(E_{1}-E_{0})$, where $E_{n}$ are the energy levels for the quantum anharmonic oscillator evaluated to first order in the strength of the quartic potential. The relation between $\lambda^{\prime}$ and our anharmonicity parameter is $\lambda^{\prime}=6\lambda E_{C}/\omega_{p}^{2}$. Thus, $\lambda^{\prime}=-E_{C}$ in a conventional tunnel junction.

In order to account for topological effects, we calculate the grand potential $U(\delta)$ defined in Eq. (2) by diagonalizing a tight-binding model of a 1D nanowire junction (see Fig. 4 for a sketch of the model). As this model has been extensively reviewed in the literature, see e.g. Ref. [cayao2018, ], we limit ourselves to providing a brief description.

We consider spinful fermions with annihilation (creation) operator $c_{j,\sigma}^{({\dagger})}$ ($\sigma\in\left\{\uparrow,\downarrow\right\}$). Introducing the spinor $\psi^{\dagger}_{j}=(c_{j,\uparrow}^{\dagger},c_{j,\downarrow}^{\dagger},-c_{j,\downarrow},c_{j,\uparrow})$, the Hamiltonian reads

where $h_{j}$ and $u_{j}$ are (respectively) the onsite and hopping matrices, with site-dependent parameters in order to create a junction. We define these matrices in terms of two sets of Pauli matrices, $\sigma_{\alpha}$ and $\tau_{\alpha}$ ($\alpha\in\left\{x,y,z\right\})$, acting respectively on the spin and the particle-hole sector of the spinor. Electrons can hop between lattice sites with a spin-independent nearest-neighbor tunneling amplitude $t$ and, due to spin-orbit coupling, a spin-flip tunneling amplitude $\alpha$. Then, the hopping matrix is

where $t_{j}=t$ away from the junction interface (see Fig. 4 and below). The onsite matrix reads

where the local mean-field superconducting gap $\Delta_{j}$ is $\Delta$ in the left superconducting lead, 0 in the weak link and $\Delta e^{i\delta}$ in the right superconducting lead. The Zeeman energy $B$ associated to the applied magnetic field is assumed to be spatially homogeneous. Likewise, the chemical potential $\mu$ is assumed to be uniform and identical in the two superconducting electrodes, though we allow for a mismatch in the chemical potentials of the electrodes and the weak link. The value $\mu=0$ corresponds to the chemical potential being in the middle of the Zeeman gap of the normal state electronic bands.

Using parameters relevant to InAs nanowires aguado2017 , we take $t=25$ meV and $\alpha=t/10\sim 2.5$ meV (corresponding to a lattice constant $a=10$ nm and an electronic effective mass $m^{*}=0.016m_{e}$, with $m_{e}$ the bare electron mass). We take throughout the proximity-induced $s$-wave superconducting gap to be $\Delta=0.9$ meV. We also neglect the dependence of $\Delta$ on $B$ and on $\mu$. In reality, this dependence is smooth across a topological phase transition, and thus clearly distinguishable from the sharp features we will be discussing below. It is well-known that variations of $B$ and $\mu$ can drive the superconducting electrodes across a topological phase transitionreviews2 . The superconducting electrodes are in the topological phase (with a MBS in each extremity of the superconductors) for $B>B_{c}\equiv\sqrt{\Delta^{2}+\mu^{2}}$, and in the trivial phase otherwise.

To model the effect of a potential barrier at the superconductor/normal interface, we use a phenomenological transparency parameter $\gamma\in[0,1]$. Both the hopping amplitude $t$ and the spin-orbit coupling parameter $\alpha$ are scaled by $\gamma$ at the interface. Throughout the text, we consider the cases $\gamma=0.6$ and $\gamma=1$ to study respectively the effect of a large potential barrier at the junction interface and the case of a clean junction with no potential barrier. We denote the former as the "tunnel" regime and the latter as the "transparent" regime. However, we note that the junction transparency also depends on the normal state Fermi velocity $v_{F}$ BTK ; heck2018 and is hence a function of $B$, $\alpha$ and most importantly $\mu$. For example, even when $\gamma=1$, the junction’s characteristics (e.g. $\omega_{p}$) can resemble those of a true tunnel junction if $\mu$ is sufficiently negative.

In the numerical simulations, unless mentioned otherwise, the weak link is two sites long ("short junction" regime), whereas each superconducting electrode contains $200$ sites. The main qualitative points are relatively insensitive to the strength of spin-orbit coupling, though larger values of $\alpha$ make it easier to observe the relevant features. Likewise, longer junction lengths do not change the main features discussed below, though they host additional properties that complicate the JBA-based detection of the topological phase transition; we elaborate on this issue below.

Figure 5 shows the dependence of the Josephson plasmon frequency $\omega_{p}$ and the anharmonicity parameter $\lambda$ on the magnetic field, for small oscillations of the superconducting phase difference in the vicinity of $\delta=0$, in the many-body ground state. We show the results both in the "tunneling" regime ($\gamma=0.6$, Fig. 5a-d) and in the "transparent" regime ($\gamma=1$, Fig. 5e-h). In all cases, the most salient feature is the pronounced peak in $\partial\omega_{p}/\partial B$ and $\partial\lambda/\partial B$ at the topological phase transition. As we explain next, these peaks can be attributed to the emergence of MBS.

In the tunneling regimetunnel1 , Fig. 5a evidences an abrupt growth of the plasmon frequency at the onset of the topological phase ($\omega_{p}$ grows by more than $0.1$ GHz as $B/B_{c}$ changes by less than $10\%$). This behavior can be understood by analyzing $\epsilon_{n}(\delta)$ as a function of $\delta$. In our model for a topologically trivial tunnel junction, all single-particle "bands" are rather flat in $\delta$ (they would be completely independent of $\delta$ if $\gamma=0$). Accordingly, $U(\delta)$ depends weakly on $\delta$ and hence $\omega_{p}$ is relatively small. In the topological phase, the single-particle bands remain flat, except for those that result from the hybridization of MBS. The latter disperse parametrically more strongly with $\delta$, because their "bandwidth" scales with the square root of the transparency, rather than with the transparency itself (as is the case for the non-Majorana bands). Hence, this explains why $\omega_{p}$ increases strongly at the onset of the topological phase, thus producing a strong peak in $\partial\omega_{p}/\partial B$ (Fig. 5b).

The aforementioned argument, derived in the context of a single-channel nanowire model, applies also to a less idealized tunnel junction that may host multiple subbands in the normal state. The simplest way to see this is by adding a term of the form $-E_{J0}\cos\delta$ to the ground state energy, where $E_{J0}$ models the contribution from the multiple subbands to the conventional Josephson energy of the junction. As a result of the additional term in the energy, $\omega_{p}$ is shifted upwards in both the trivial and the topological regimes, but the ramp up of $\omega_{p}$ and the peak of $\partial\omega_{p}/\partial B$ at the onset of the topological phase remain unchangedEJ0 .

A related behavior has been identified for the DC critical current in earlier works, with a marked increase at the onset of the topological phase sanjose2014 ; cayao2017 ; huang2017 . At first glance, a similarity in the behavior of $\omega_{p}$ and the critical current is not surprising, the latter being proportional to $\omega_{p}^{2}$ in conventional Josephson junctions. Yet, for the unconventional junctions we are interested in, this relation of proportionality does not apply. Instead, all one can assert is that $\omega_{p}$ depends on the characteristics of the single-particle bands near $\delta=0$, while the critical current depends on the slope of the single-particle bands near large ($\simeq\pi$) values of $\delta$. Moreover, from an experimental point of view, the measurement of $\omega_{p}$ is qualitatively different from the measurement of the critical current; it can therefore offer an alternative way to probe the interpretation of a recent experiment reporting an enhanced critical current at the onset of the topological phasetiira2017 .

Let us now discuss the anharmonicity parameter. Still in the tunneling regime, Fig. 5c reflects a peculiar evolution of $\lambda/\omega_{p}^{2}$ across the topological transition, which can be understood analytically. Deep in the trivial regime, we find $\lambda/\omega_{p}^{2}=-1/6$, as expected for a conventional Josephson junction with $U(\delta)=-E_{J}\cos\delta$. Deep in the topological regime, $\lambda/\omega_{p}^{2}$ tends asymptotically to $-1/24$, as expected for $U(\delta)=-E_{M}|\cos(\delta/2)|$ (cf. Sec. II.2). In other words, in a tunnel junction, the evolution of $\lambda/\omega_{p}^{2}$ across the topological phase transition reveals the emergence of MBS.

Figure 5 also displays the results for several values of $\mu/\Delta$. Even though we consider $\mu$ to be homogeneous throughout the system (i.e., $\mu=\mu^{\prime}$ in Fig. 1b), we have verified that the results do not change significantly when we vary $\mu$ in the superconducting electrodes while keeping $\mu^{\prime}$ pinned to zero in the weak link. We note that $\omega_{p}$ increases with $\mu/\Delta$ in Fig. 5. This effect can be understood by observing that, for fixed $\gamma$, the junction transparency is a monotonically decreasing function of the ratio of the barrier amplitude and the normal-state Fermi velocity $v_{F}$ BTK . As $\mu/t\ll 1$ for all the cases considered (Fermi level near the bottom of the band), increasing $\mu$ means a higher $v_{F}$ and hence an increase of the junction transparency. Although Fig. 5 shows only positive values of $\mu$, we have checked that, for $\mu<0$, the junction transparency effectively decreases. For $\mu/\Delta\lesssim-4$, the behavior expected for a tunnel junction is found irrespective of the value of $\gamma$ as the Fermi level in the normal state then lies below the bottom of the band (insulating regime).

Another effect of increasing $\mu/\Delta$ is the appearance of oscillations in $\omega_{p}$ and $\lambda$ as a function of $B$; their amplitude grows with $B$. These oscillations are the microwave counterpart of the ones predicted in the DC critical current by Cayao et al.cayao2017 The oscillations take place only in the topological phase and arise from the hybridization between MBS localized at opposite extremities of each superconducting electrode, as well as from the hybridization between the MBS localized at the two opposite extremities of the entire system. Accordingly, the oscillations are washed out when the length of the superconducting electrodes becomes long enough (we find no evidence of them when doubling the electrodes size to $400$ sites).

Let us now discuss the transparent junction regime $\gamma=1$ (see Fig. 5e-h). In this regime, the behavior of $\omega_{p}$ at the topological phase transition is also reminiscent of that predicted in earlier works for the DC critical current.sanjose2013 ; cayao2017 In those works, a kink of the critical current at the topological phase transition was identified and attributed to a band inversion taking place at the topological phase transition. A similar mechanism is at play in our case, as we explain next.

The underlying explanation begins by recognizing that a long superconducting electrode in the trivial phase hosts two different and quasi-independent $p$-wave superconducting gaps in the energy spectrumsanjose2013 ; cayao2017 ; murthy2020 : one at the inner Fermi points (the "$\Delta_{-}$ gap") and the other at the outer Fermi points (the "$\Delta_{+}$ gap"). When $B\ll\alpha$, $\Delta_{+}\simeq\Delta$. As $B$ increases, $\Delta_{+}$ decreases gradually and gets significantly suppressed when $B\gg\alpha$. The $\Delta_{-}$ gap is far more sensitive to $B$, vanishing when $B=B_{c}$. This gap-closing point marks the topological phase transition. Near the transition, we have $\Delta_{-}\simeq|B-B_{c}|$.

The next part of the explanation is to recall that in a short, transparent and topologically trivial junction, there is a single Andreev bound (or quasi-bound) state associated to each of the two $p$-wave gaps of the electrodestunnel . In this regime, the rest of the single-particle states (the so-called scattering states) are largely dispersionless in $\delta$, and therefore the main contribution to $\omega_{p}$ and $\lambda$ originates from the two Andreev states bound by $\Delta_{+}$ and $\Delta_{-}$. Although there is no simple analytical expression for the Andreev bound state dispersions in the presence of a generic magnetic fieldheck2018 , qualitative insight can be gained by neglecting the coupling between the two $p$-wave gaps, adopting the Andreev approximation, and assuming perfect transparency. Then, we posit $\epsilon_{\pm}(\delta)\simeq\Delta_{\pm}\sqrt{1-\sin^{2}(\delta/2)}$ for the two ABS in the trivial phase, and we get

Across a topological phase transition from the trivial to the topological phase, $\Delta_{-}$ crosses zero and inverts its sign, thereby changing its character from a superconducting gap to a magnetic gap. Accordingly, in the topological phase, $\Delta_{-}$ no longer binds an Andreev state. Therefore, the contribution from $\Delta_{-}$ to $\omega_{p}$ and $\lambda$ in Eq. (III.1) is turned off for $B\geq B_{c}$. This results in a kink for $\omega_{p}$ and $\lambda$ as a function of $B$ at $B=B_{c}$ (and corresponding discontinuities in $\partial\omega_{p}/\partial B$ and $\partial\lambda/\partial B$). Of course, such reasoning of "suddenly turning off" the contribution from $\Delta_{-}$ at $B=B_{c}$ makes sense only at zero temperature; finite-temperature effects will be discussed in Sec. III.2.

Some of the features in Fig. 5 match qualitatively with the behavior expected from Eq. (III.1). First, $\omega_{p}$ is decreasing as we approach the topological phase transition. Second, a kink is visible at the transition, especially for larger values of $\mu/\Delta$, and $\omega_{p}$ continues to decrease as we get deeper in the topological phase (because $\Delta_{+}$ keeps decreasing as $B$ is made larger). On the other hand, the results from Fig. 5g do not match qualitatively with Eq. (III.1). For one thing, $\lambda/\omega_{p}^{2}$ is not a constant. Also, a clear kink in $\omega_{p}$ at the topological phase transition is present only if $\mu/\Delta$ is sufficiently large. There are various possible reasons for these discrepancies.

First, the analytical expressions for $\epsilon_{\pm}(\delta)$ are approximately valid at $B=0$, but may fail qualitatively at larger $B$. Even at $B=0$, $\gamma=1$ does not guarantee a perfect transparency of the junction; the actual transparency depends on the value of $\mu$ as well. If the transparency is not unity, one expects various competing harmonics in $\epsilon_{\pm}(\delta)$ near $\delta=0$ (rather than just $\cos(\delta/2)$), which in turn leads to $\lambda/\omega_{p}^{2}$ not being constant. Second, the Andreev approximation upon which Eq. (III.1) is based is not satisfied in the low electron density regime required to have a topological Josephson junction; this is evidenced by the single-particle energy spectra we have calculated (not shown), wherein the ABS in the trivial phase have significant gaps at $\delta=\pi$. Third, near $\delta=0$, the Andreev state bound by $\Delta_{+}$ is buried in the continuum of scattering statessanjose2013 ; murthy2020 ; this is due to the fact that $\Delta_{+}\gg\Delta_{-}$ in the vicinity of the phase transition. In other words, the ABS bound by $\Delta_{+}$ undergoes multiple avoided crossings with scattering states near $\delta=0$. Consequently, it is not accurate to neglect the $\delta$-dependence of those scattering states, nor is it to ignore their contribution to $\omega_{p}$ and $\lambda$.

Next, we compare our results to relevant works in the literature. Our finding that $\partial\omega_{p}/\partial B$ and $\partial\lambda/\partial B$ display prominent peaks at $B=B_{c}$ is seemingly more optimistic than that of a recent numerical study by Keselman et al.keselman2019 , who conclude that there are "no strong signatures of the topological transition itself on the simulated frequency spectra of the junction." This discrepancy in viewpoint might be explained by the different parameter regimes considered. In particular, the exact numerical treatment of charging energy in Ref. [keselman2019, ] comes at the expense of studying smaller systems ($\sim 40$ sites), where finite-size effects smoothen the phase transition into a crossover. Below, we show that finite-temperature effects can likewise broaden and eventually erase the sharp signatures in $\partial\omega_{p}/\partial B$. In another recent workavila2020 , Avila et al. have, to the contrary, reported an abrupt dip in the anharmonicity parameter $\lambda^{\prime}$ at $B=B_{c}$, stating it to be a "precise smoking gun" for the topological phase transition. We recall (see discussion at the end of Sec. II.2) that $\lambda^{\prime}=6E_{C}\lambda/\omega_{p}^{2}$.

The theoretical treatments of Avila et al. and Keselman et al. differ from ours in that they treat the charging energy quantum mechanically, while we do so semiclassically. Yet, both approaches ought to yield similar results in the transmon regime. In Fig. 5g, we do find that $\lambda/\omega_{p}^{2}$ can have a dip at $B=B_{c}$ in the transparent regime, though it is nowhere as pronounced as in Fig. 9b of Ref. [avila2020, ]. Moreover, the dip is not a generic feature when we vary parameters such as $\gamma$ and $\mu$ (e.g., it is absent in Fig. 5c). What appears more generic in our case is the peak in $\partial\lambda/\partial B$ at $B=B_{c}$, whose origin we have explained above.

On a related note, while Avila et al. investigate the behavior of the nonlinear Josephson inductance in the topological and trivial phases, they do not point out specific signatures at the phase transition (at least Fig. 8 of Ref. [avila2020, ] does not show enough values of $B$ to ascertain any signatures). Yet, the inverse of this nonlinear inductance, close to $\delta=0$, is proportional to $\omega_{p}^{2}$. It would be interesting to find out whether the theory of Ref. [avila2020, ] agrees with our prediction for $\partial\omega_{p}/\partial B$ near $B_{c}$.

Thus far, we have concentrated on short junctions either in the tunneling or in the transparent regimes. In junctions of intermediate transparency, $\omega_{p}$ does not display a kink or a sharp upturn at the topological phase transition; the peak of $\partial\omega_{p}/\partial B$ at $B=B_{c}$ can then be less pronounced. Even in those cases, we still observe a clear feature in $\partial\lambda/\partial B$ at $B=B_{c}$.

Longer junctions show the same phenomenology as the short junctions when it comes to the behavior of $\omega_{p}$ and $\lambda$ at the topological phase transition (we have simulated weak links with up to $50$ sites). However, one additional aspect of the longer junctions is that they may host accidental (imposed neither by electronic topology nor by symmetry), parameter-sensitive, zero-energy touchings at $\delta=0$ in the topologically trivial phasevuik2019 . Because of particle-hole symmetry inherent to the single-particle spectrum, the positive-energy and negative-energy Andreev states that cross at zero energy have opposite curvatures and anharmonicities. As a result, the curvature of the occupied (negative-energy) single-particle states changes abruptly as a function of $B$ at the accidental zero-energy touchings. A corresponding discontinuity takes place in $\omega_{p}$ and $\lambda$ (see Fig. 6 for an example). The magnitude of the discontinuities is diminished when the $\delta$-dependence of the ABS at the zero-energy touching is smaller; hence, the effect is generally less pronounced in tunnel junctions. We note that no such discontinuities occur at the true topological phase transition, where positive- and negative-energy Andreev bound states coalesce at zero energy (as opposed to crossing it) as a function of $B$. Hence, measuring $\omega_{p}$ and $\lambda$ as a function of $B$ may allow to distinguish an accidental gap closing in the trivial phase (with a discontinuity in $\omega_{p}$) from a gap closing associated to a topological phase transition (with no discontinuity in $\omega_{p}$ but a peak in $\partial\omega_{p}/\partial B$).

Unfortunately, the aforementioned distinction is blurred at finite temperatures, where the discontinuities are transformed into rapid but continuous variations. Hence, at nonzero temperature, the behavior of $\omega_{p}$ and $\lambda$ at accidental zero-energy touchings in the trivial phase can mimic the signatures of the topological phase transition. In this case, one way to disentangle the true MBS signatures is to analyze the dependence of the peaks in $\mu$; experiments of this type are feasiblefrolov2017 . The MBS-related peak at the phase transition will be displaced in $B$ as $B_{c}=\sqrt{\mu^{2}+\Delta^{2}}$, while the peaks of trivial origin will not generically show such dependence, and can even disappear completely as $\mu$ is varied (see Fig. 6).

To conclude this subsection, we remark that our main results do not rely on the fact that a magnetic field has been used to induce a topological phase transition. Similar conclusions are reached when $B$ is held fixed and the topological phase transition is driven by variations in $\mu$. In that case, we find that $\partial\omega_{p}/\partial\mu$ and $\partial\lambda/\partial\mu$ show strong peaks at the topological phase transition, for the same reasons as the ones alluded to above.

In the preceding subsection, we have investigated $\omega_{p}$ and $\lambda$ at zero temperature, and have identified signatures of the topological phase transition therein. In this subsection, we study the effect of finite temperature on those features. We concentrate on temperatures that are low compared to the proximity-induced $s$-wave superconducting gap $\Delta$. We find that the peaks of $\partial\omega_{p}/\partial B$ and $\partial\lambda/\partial B$ at the topological phase transition are thermally broadened, though they remain noticeable at experimentally attainable temperatures. We also find a nontrivial temperature dependence of $\omega_{p}$ in the topological phase of tunnel junctions, which can be ascribed to the MBS localized on the opposite extremities of the weak link.

Our results are summarized in Fig. 7. The main aspects of this figure can be understood from the finite temperature expressions (derived from Eqs. (2) and (6)) for the plasmon frequency and the anharmonicity parameter near $\delta=0$,

where $\epsilon_{n}(0)$ is the energy of the $n$-th single particle state at $\delta=0$ (recall that $\epsilon_{n}(\delta)\geq 0$ for $n>0$), and $\epsilon_{n}^{(k)}(0)$ denotes the $k$-th derivative of $\epsilon_{n}$ with respect to $\delta$ (evaluated at $\delta=0$).

It follows from Eq. (III.2) that, whenever the energy of the lowest-lying $n>0$ state at $\delta=0$ largely exceeds $k_{B}T$, $\omega_{p}$ and $\lambda$ will have a negligible temperature-dependence. In other words, the condition for a significant temperature-dependence of $\omega_{p}$ and $\lambda$ is that $k_{B}T\gtrsim\epsilon_{1}(0)$ (with $\epsilon_{1}(0)$ the smallest single-particle eigenenergy) Another useful piece of information in Eq. (III.2) is that, for a fixed $n$, the contribution of the given single-particle level to $\omega_{p}$ is suppressed when $k_{B}T\gg\epsilon_{n}(0)$.

Let us now analyze in detail the content of Fig. 7. We divide the discussion in three parts, according to three regimes: (i) the trivial phase, (ii) the onset of the topological phase, and (iii) deep in the topological phase.

In the trivial phase of short junctions, there are no subgap states in our model at $\delta=0$ (regardless of the transparency) vuik2019 . Therefore, $\epsilon_{1}(0)=\Delta_{-}$ and the temperature-dependence of $\omega_{p}$ and $\lambda$ is strong only if $k_{B}T\gtrsim\Delta_{-}$. At low temperatures, such condition is satisfied only close to the topological phase transition (because $\Delta_{-}\to 0$ when $B\to B_{c}$). This explains why, in Figs. 7a, c, e and g, the curves for different temperatures are superimposed deep in the trivial phase, but begin to diverge closer to the topological phase transition. The fact that $k_{B}T\gg\Delta_{-}$ is easily attained at the topological phase transition also explains why the zero-temperature kinks in $\omega_{p}$ and $\lambda$, which take place at the topological phase transition in highly transparent junctions (Fig. 5e-h), are washed out by thermal effects. Still, remnants of the kinks survive at low ($k_{B}T\ll\Delta$) but finite temperatures, in the form of broadened peaks for $\partial\omega_{p}/\partial B$ and $\partial\lambda/\partial B$. These peaks are evident in Figs. 7b, d, f and h.

At the onset of the topological phase ($B\simeq B_{c}$), four MBS emerge. Two of them, localized on the outer extremities of the superconducting electrodes, have near-zero energies that are essentially independent of $\delta$ in the vicinity of $\delta=0$; these dispersionless states do not contribute to $\omega_{p}$ and $\lambda$. For this reason, we will hereafter omit the two "outer" MBS from the discussion. In contrast, the two "inner" MBS localized on the opposite extremities of the weak link have a non-negligible bandwidth in $\delta$, given by the hybridization energy $E_{M}$; these dispersive states can contribute appreciably to $\omega_{p}$ and $\lambda$. In a junction of high transparency, $\epsilon_{1}(0)=\Delta_{-}\ll E_{M}\simeq\Delta_{+}$. Hence, in such junction, the leading temperature-dependence of $\omega_{p}$ and $\lambda$ at $k_{B}T\gtrsim\Delta_{-}$ originates mainly from the non-Majorana states. In contrast, in a junction of low transparency, $\epsilon_{1}(0)=E_{M}<\Delta_{-}\ll\Delta_{+}$ and the leading temperature-dependence at $k_{B}T\gtrsim E_{M}$ originates from the inner MBS.

As the system is driven deeper into the topological phase, Fig. 7 reflects two possible scenarios. In the scenario of Fig. 7e-h, the junction is highly transparent and there are no dispersive subgap states at $\delta=0$. In this case, $\epsilon_{1}(0)={\rm min}(\Delta_{-},\Delta_{+})$, where we once again ignore the outer MBS. If $\Delta_{-}<\Delta_{+}$ (which is common not far from the topological phase transition), $\epsilon_{1}(0)=\Delta_{-}$ grows linearly with the magnetic field. Therefore, the temperature-dependence of $\omega_{p}$ and $\lambda$ weakens gradually as $B$ grows, and becomes once again negligible when ${\rm min}(\Delta_{-},\Delta_{+})\gg k_{B}T$. This is why, in Fig. 7e and g, the curves for different temperatures tend to converge deep in the topological phase.

Contrastingly, in the scenario of Fig. 7a-d, the junction has low transparency. In this case, there are dispersive subgap MBS of energy $E_{M}$ at $\delta=0$, so that $\epsilon_{1}(0)=E_{M}$. Moreover, the low junction transparency implies $E_{M}\ll\Delta_{-}$. Then, as temperature is raised well above $E_{M}$, the contribution of MBS to $\omega_{p}$ is suppressed via the $\tanh(E_{M}/2k_{B}T)\ll 1$ factor in Eq. (III.2). At the same time, the energies of the higher-energy (non-Majorana) scattering states are relatively dispersionless in $\delta$ and therefore make a relatively modest contribution to $\omega_{p}$, much like in the trivial phase. As a result, $\omega_{p}$ is strongly reduced in the topological phase when $k_{B}T\gg E_{M}$. This explains the traits of Fig. 7a in the topological phase.

One intriguing outcome from the foregoing discussion is that it appears to be possible to thermally switch on and off the contribution from MBS to $\omega_{p}$ in short tunnel junctions, by going from $k_{B}T\ll E_{M}$ to $k_{B}T\gg E_{M}$. This statement applies to longer junctions of low transparency as well. Although such junctions host low-energy Andreev-bound states in the trivial phase, their $\delta$-dependence in the tunneling regime is relatively small; therefore, changing their thermal occupancy does not make a substantial difference to $\omega_{p}$. We have verified this point for junctions with 20 and 50 sites in the weak link, and $|\mu|\ll\Delta$.

Thus far, we have been concerned with the prediction and understanding of peaks in the derivatives of $\omega_{p}$ and $\lambda$ across a topological phase transition. The connection between $\omega_{p}$ and $\lambda$ with the operation of a JBA is made by Eq. (7). From this equation, it is clear that the bifurcation currents of the TJBA will display signatures of a topological phase transition.

The results are summarized in Fig. 8, which shows the lower and upper bifurcation currents as a function of the magnetic field, for short junctions of low (Fig. 8a-d) and high (Fig. 8e-h) transparency. The drive frequency $\omega_{d}$ is chosen to be below $\omega_{p}$, as otherwise there is no bifurcation (cf. Sec. II).

The bifurcation currents are much smaller in the tunneling regime than in the transparent regime, because the critical current of the junction is much lower in the former. At the topological phase transition, the behaviour of bifurcation currents as function of magnetic field parallels that of $\omega_{p}$, and thus their derivatives display peaks at $B=B_{c}$. These peaks are broadened at finite temperature.