Parametric amplification of topological interface states in synthetic Andreev bands

A driven-dissipative nonlinear photonic system (e.g. exciton-polaritons) can operate in a gapped superfluid regime. We theoretically demonstrate that the reflection of a linear wave on this superfluid is an analogue of the Andreev reflection of an electron on a superconductor. A normal region surrounded by two superfluids is found to host Andreev-like bound states. These bound states form topological synthetic bands versus the phase difference between the two superfluids. Changing the width of the normal region allows to invert the band topology and to create"interface"states. Instead of demonstrating a linear crossing, synthetic bands are attracted by the non-linear non-Hermitian coupling of bosonic systems which gives rise to a self-amplified strongly occupied topological state.

Topological physics relies on the specific structure of the eigenstates of Hamiltonians in a parameter space. It became one of the most active fields of research of the last decades. Topological invariants were successively used to characterise superfluid excitations (quantum vortices) [1,2], solitons in polyacethylene [3], Landau levels [4], electronic Bloch bands in solids [5], and more generally energy bands in periodic media [6]. A related concept is the bulk-boundary correspondence [7,8], which associates the change of the band topology with the gap closing and therefore with the existence of an interface state between media of different topologies. Such interface states can demonstrate unique properties, such as the chiral (one-way) transport in topological insulators [9].
The recent proposals and implementations of topological media supporting protected edge states in nonlinear systems opened new possibilities [10][11][12][13][14][15][16][17][18]. The most well-known examples are the so-called topological lasers, where lasing occurs at the topological edge states of a photonic lattice [19][20][21][22]. Another recent orientation is synthetic topological matter, where the parameters of the Hamiltonian are not imposed by the system itself (such as the wave vector in Bloch bands), but rather externally, by experimental conditions. This approach allows to explore regimes which cannot be accessed otherwise, such as the 4D quantum Hall effect [23]. One possibility to implement synthetic topological matter relies on using Josephson junctions between conventional superconductors [24]. The dependence of Andreev bound states versus the phase difference between the superconductors which form the junction allows to define synthetic bands with the synthetic dimensionality controlled by the number of terminals. These bands were found to be topological, showing Weyl singularities, and significant efforts are made to implement them [25,26]. Andreev reflection occurs at the interface of a supercon-ducting medium [27], where an incoming electronic wave undergoes an anomalous reflection, where its wavevector is reversed and is interpreted as a hole excitation associated with opposite charge and mass. The analogy between Andreev reflection and the reflection of a wave over a Bose-Einstein Condensate has been studied theoretically [28]. In photonic systems, one can implement a superconductor analog using resonant driving, typical for cavity exciton-polaritons [29,30]. In such a case, a gap opens in the spectrum of elementary excitations of the pumped modes which corresponds to the formation of a "gapped superfluid". It is, moreover, possible to modulate the pump intensity in space to realize two or more gapped regions with controllable phase, separated by normal (non-superfluid) regions [31][32][33].
In this work, we study the analog of the Andreev reflection of an incident wave on a gapped driven-dissipative polariton superfluid. We show the existence of Andreevlike bound states between the two driven superfluids. The topology of the synthetic bands defined versus the phase difference between the two superfluids is found to be non-trivial. The interface states between the regions of inverted topology are found to be self-amplified. They are macroscopically-occupied topologically protected Andreev-like states.
The model. We consider a resonantly-pumped 1D nonlinear optical media described by the driven-dissipative Gross-Pitaevskii equation for scalar particles, which is similar to the Lugiato-Lefever equation [34]: where P = e −iωpt is the quasi-resonant pumping term at normal incidence, α > 0 is the repulsive interaction term, and γ the decay rate. The pumped mode dispersion is parabolic (mass m), typical for photons or excitonpolaritons in a microcavity. All polarisation effects are arXiv:2012.05709v1 [cond-mat.mes-hall] 10 Dec 2020 neglected. The spatially-homogeneous solution ψ s e −iωpt shows a bistable behaviour [35] with a large occupation of the pumped mode above the bistable threshold. To simplify the analytics, we consider the limit γ → 0 (see [36]). The wavefunction describing weak excitations of the macro-occupied mode reads: Here, one plane-wave excitation is coupled by the nonlinear term to its complex conjugate. Inserting this wavefunction into Eq. (1) allows one to derive the Bogoliubovde Gennes equations: where L = ( k − E p + 2αn). The coefficients u, v are the Bogoliubov coefficients, k = 2 k 2 /2m and ψ s = √ ne iφ . The energy E = ω, measured with respect to the laser detuning E p = ω p , is found by cancelling the determinant and reads: A bogolon is formed by the superposition of two plane waves of amplitude u and v * at energies E, −E. Considering E positive, the ratio |u|/|v| is fixed by Eq. (3) and is larger than 1. The specific choice of a normalisation condition |u| 2 − |v| 2 = 1 allows to define a bogolon as a particle of positive total energy (|u| 2 − |v| 2 )E, but containing amplitudes at both E and −E. When αn > E p , the spectrum of allowed energies (containing both positive and negative parts) shows a gap of magnitude 2∆ centered on the pump energy: Next, we consider two semi-infinite regions ( Fig. 1(a)). The left part is not pumped and characterized by a parabolic dispersion. The right part is resonantly pumped and described by the 1D GP equation. We consider an incident wave coming from the normal part to the pumped area at an energy within the gap. We have therefore to consider evanescent bogolon modes, instead of the usual propagative ones. The corresponding wavefunction reads: ψ(x, t) = e −iωpt ψ s + ue −κx e −iωt + v * e −κx e +iω * t .  Importantly, at a given positive E, there exist two different evanescent waves with two different inverse decay lengths and two different eigenvectors: The "minus" state has a longer decay length and |u − | > |v − |. Its dominant component has an energy E > 0 and the state is normalized as |u − | 2 − |v − | 2 = 1. It is quite similar to propagative bogolons, continuing their dispersion within the gap, as shown on Fig. 1(c). We call this type of state, where the positive energy component is larger, a "particle". The "plus" state has a shorter decay length and |u + | < |v + |. Its dominant component has an energy −E < 0 as |u + | 2 − |v + | 2 = −1. This evanescent solution has no propagative counterpart [36]. The κ + branch shown in Fig. 1(c) is disconnected from the propagative states dispersion. This type of solution never described before to our knowledge can be assimilated to a "hole" state. It is associated with a local decrease of the particle density with respect to the homogeneous superfluid. It plays a crucial role in the Andreev-like reflection, as described below. As shown in Fig. 1(b), we consider a plane wave of energy 0 < E < ∆ incident on the pumped region. It excites the two above-mentioned types of evanescent bogolons, which provoke reflection at the energies E and −E, respectively. The wavefunction in the normal (left) region reads: with k 1,2 = 2m (E p ± E)/ , valid for ∆ < E p (see [36] for other cases with evanescent states in the normal region). In the superfluid, the wavefunction combines two evanescent waves: To compute the reflection coefficients, are the normal and Andreev reflection coefficients for an incident "particle" (dominant positive energy component). The group velocities allow to conserve the current. Similarly,

the reflection coefficients
for an incident "hole" having a dominant negative energy component. The continuity of the wavefunctions and of their derivatives at the interface gives an analytical expression for these reflection coefficients (see [36]) which are plotted in Fig. 1(c) for parameters m = 5 × 10 −5 m 0 , E p = 0.5 meV, and αn ≈ 1.11E p , characteristic for GaAs exciton-polaritons. This yields a gap ∆ ≈ 0.5E p . The Andreev-like reflection is comparable in amplitude to the normal reflection. Ultimately, the phenomenon occurring here is very similar to the Andreev reflection, but can also be interpreted as a non-linear frequency conversion. An incoming wave at the frequency ω p ± δω is partially reflected both at ω p ± δω and at ω p ∓ δω. In the case of in-gap energies, the reflection (normal and Andreev together) is total, since |r N p,h | 2 + |r Ap,h | 2 = 1. Optical phase conjugation, discovered and studied in the 70's in nonlinear optics [37] also shows a strong analogy with Andreev reflection [38,39].
Andreev bound states analog. The next step is to consider a Superfluid -Normal -Superfluid (SNS) junction ( Fig. 2(a)). The only difference between the two superfluids is their phase, φ L = 0 and φ R ≡ φ, respectively. The width of the normal region is a. This structure exhibits trapped states (similar to quantum well eigenstates), determined by the density profile even without the Andreev reflection. The latter provides a correction and generates a second energy component (bogolon image) for each of these states.
The energy components of the Andreev bound states are computed using the scattering matrix formalism. The scattering matrices for the reflection phenomena at each interface read: The eigenvectors of S T determine the wavefunction via Eqs. (S12), (10) and allow one to determine if a state is either particle-like or hole-like. Depending on parameters, the eigenenergies can be either real (stationary Andreev-like bound states) or imaginary (self-amplified bound states). Figures 2(b) and (c) show two examples of hole-like bound states. In Fig. 2(b), φ = 0 and both energy components have the same parity (s-like state). In Fig. 2(c), φ = π/2. The main component keeps its parity, whereas the bogolon image becomes p-like, because of the phase shift. We compute the probability current J ± in the normal region [36], similar to the Josephson current in superconducting junctions: The total probability current of one component at the energy +E is fully compensated by the current at the energy −E. Topological synthetic bands. Both positive and negative energy components of a bound state form synthetic energy bands with respect to φ. For a given energy band, the wavefunction is a superposition of two counterpropagating plane waves, as defined by Eq. (S12). The amplitudes of these two plane waves define a pseudospinor X + ∼ (A, B) T for positive energies and X − ∼ (C, D) T for negative energies. The evolution of these pseudospinors along the band allows to compute the Zak phase [41,42]: Figure 3(a-c) shows the synthetic bands and their Zak phases for three thicknesses a. In Fig. 3(a,c) all energies are real. The bands shown in blue correspond to a particle-like state with a dominant positive energy part (solid lines). Its negative energy counterpart (bogolon image) is smaller in amplitude (dashed lines) and shows a non-zero Zak phase. Indeed, at φ = 0, C = ±D. The current is zero, and the associated pseudospin lies in the plane of a Bloch sphere representation. At φ ≈ π/4, D = 0. The current is maximal, and the pseudospin points towards the pole. Finally, at φ = π/2, C = ∓D: the symmetry of the state has changed. Between φ = −π/2 and φ = π/2, the pseudospin covers a full great circle of the Bloch sphere (constrained by the mirror symmetry of the problem, see [36]), and the accumulated Zak phase is π. On the other hand, the pseudospin of the majority component slightly moves towards the pole at φ ≈ π/4 to go back to its original position at φ = π/2. In general, the band associated with the dominant energy component of the bogolon (the original trapped state) shows a null Zak phase, whereas the minority component (the bogolon image) is topologically non-trivial. Fig. 3(a,c) show the band topology inversion between two values of a which implies the gap closing for a critical thickness a 0 . This topological band crossing normally gives rise to Dirac, or Weyl points, depending on the system's dimensionality. It is also at the heart of topologically protected edge or interface states. In our bosonic system of interacting particles, crossing bands interact through a non-linear non-Hermitian coupling [Eq. (3)]. Fig. 3(b) shows the states at the critical thickness a 0 , where the gaps are closing. Instead of simply crossing, the bands merge: their real parts become flat, while the imaginary parts are opposite. The positive imaginary part means that the topological state corresponding to the crossing is amplified and becomes strongly populated. This amplification occurs when the two bogolon components are resonant with the linear eigenstates of the potential trap. Fig. 3 (d) and (e) show for φ = π/2 and φ = 0, respectively, the mode energy versus a. For φ = 0, amplification occurs when the majority and minority component have the same parity [(n, n)], when the potential well states n = 1, 2, . . . are resonant with the laser. For φ = π/2, the amplification occurs when the bogolon image of a state of given parity becomes resonant with an original trapped state of different parity.
Numerical simulations. To confirm the analytical theory, we perform numerical simulations, solving Eq. (1) over time with a weak probe exciting the bogolon states, and a finite lifetime τ = /2γ = 30 ps. The width a varies versus y as: a(y) = a 0 + y/l, where l = 64 µm. The calculated spectra are presented in Fig. 4 for a < a 0 for φ = 0 (a) and φ = π/2 (b). The two energy components of each of the Andreev states are marked in the Figure (e.g. 1 is the "original" trapped state and 1 is its bogolon image). The change of the symmetry of the bogolon images 1 and 2 is clearly visible. The gap closes at a = a 0 , where each original state is resonant with the bogolon image of the other: E 2 = E 1 and E 1 = E 2 , and their symmetries also coincide for φ = π/2. As a result, the Andreev states are amplified and dominate the spectrum (Fig. 4(d)). Figure 4(c) shows the spatial distribution of emission, with the amplified Andreev state visible at y = 0. We stress that the normal part of the junction is a non-interacting medium, otherwise solitons are formed there [31][32][33], and the approximations used to make our calculations are not valid anymore.
To summarize, we predict an analogue of the Andreev reflection in a photonic driven-dissipative gapped superfluid. Such systems can be used to form bosonic SNS junctions hosting Andreev-like bound states. These bound states form topologically nontrivial synthetic bands as a function of the phase difference between the pumping lasers. By changing the width of the normal region, one can invert their topology, and the associated topologically protected interface states are found to be self-amplified, giving rise to strongly emitting topologically protected photonic modes.
In the main text, the spectrum of elementary excitations of bogolons is proven to present a gap. This gap directly comes from the dispersion relation: (S1) At first glance, in the propagative case, this provides two solutions for the norm of the wavevector |k| (that is, four solutions for k): However, the solution k − is actually imaginary, since in the propagative case E > ∆, which yields (αn) 2 + E 2 > 2αn − E p . Thus, there are only two solutions with the same norm, but opposite propagation direction: This is obviously different from the case of evanescent bogolons considered in the main text, where the two inverse decay lengths are different.
For comparison, the energy E = ω of a bogolon in the evanescent case can be expressed with respect to the inverse decay length κ: This expression should be compared with Eq. (S1).
Finally, we note that the descriptions of the normal region based on the Schrödinger equation and on the Bogoliubov-de Gennes equations (BDG) with exactly zero interactions are equivalent. While it may seem that the BDG equations have two solutions with opposite energies, zero interactions mean that the BdG matrix is already diagonal and the negative energy enters the wavefunction only with the minus sign. Thus, these two components simply correspond to two waves propagating in opposite directions with the same energy. This energy is positive when measured from the bottom of the band.

Andreev reflection
In this section, we present explicit expressions for the reflection coefficients for normal and Andreev reflection. We start by commenting the limit of vanishing lifetime used for the analytical calculations. A finite γ reduces both reflection coefficients. Its effect is stronger for the case of large penetration length 1/κ ± , discussed in the main text. This occurs especially for E → ∆. Therefore, the analytical results that we obtain should not be applied for energies close to the edge of the gap and for particularly narrow gaps, ∆ → 0.
As said in the main text, both the wavefunctions and their derivatives have to be continuous on the interface. This imposes: Matching the wavefunctions and their derivatives at the interface gives an analytical expression for the reflection coefficients: The precedent expressions can be reformulated to make the phase appear explicitly by considering the group velocities v 1,2 = k 1,2 /m and the relation be-tween the Bogoliubov coefficients u + = −v − e 2iφ and v + = u − e −2iφ : (S7) This form makes appear explicitly the role played by φ, and more specifically the change of sign of the Andreev reflection coefficient when φ = π/2. The reflection coefficients for holes are given as r N h (E) = r Np (−E) and The dependence on the different energies at stake is here hidden in the complexity of the formulas. However, for E E p , an approximate expression can be given for both coefficients: With these expressions, we clearly notice that the three crucial energies to consider are E p , αn−E p and 3α n −E p . It allows to find the asymptotic values of the reflection coefficients for ∆ → 0. The maximal value for Andreev reflection coefficients, |r A | 2 = 2/3 = 2|r N | 2 , is achieved for αn → E p (∆ → 0). However, we note that these values are beyond the domain of the validity of the theory, since for E = ∆ = 0 any finite decay γ plays a nonnegligible role. For realistic ∆, |r A | 2 < |r N | 2 .
In the main text, we also consider the case of a SNS junction with a normal region of width a. In the derivation of the energy components of a bound state, based on the scattering matrices of the interfaces formed by the reflection coefficients discussed above, we also need the expression of the scattering matrix describing the propagation of the wave in the normal region S N . Regardless of the direction of propagation, this matrix reads: This matrix, together with the ones describing the reflection processes on both interfaces, allow one to find a condition for the existence of a bound state which takes the form of an cancellation of a determinant (see main text). This condition is equivalent to: (S10) where k ± = k 1 ± k 2 , which is more compact and makes appear the reflection coefficients explicitly.

Zak phase
To compute the Zak phase, we use: As in the main text, the wavefunction in the normal (central) region is written as: and the associated vectors for each energy component ±E are written as: where the coefficients are the ones defined in (S12) but normalized to one (|A | 2 +|B | 2 = 1 and |C | 2 +|D | 2 = 1). The coefficients themselves are computed numerically. Such vectors can indeed be plotted on the Bloch sphere.
The vector corresponding to the dominant energy always gives Φ Zak = 0 while the other one gives Φ Zak = π.
We note that both spinors are constrained to the great circle of the Bloch sphere by the symmetry of the probability density distribution. Indeed, since the problem is completely symmetric with respect to x = 0, the probability density has to exhibit mirror symmetry with respect to this point. For this, the relative phase between A and B and also between C and D has to be either 0 or π, which means that the pseudospin can only make a circle through the constant longitude plane (azimuthal angles 0 • and 180 • ). If one allows an arbitrary phase between these coefficients, the probability density is shifted and becomes asymmetric. Indeed, (S14) If the phase difference between A and B is 0 or π, we can assume that both are real. In this case, the probability density simply writes which is symmetric. We can now introduce the phase difference between the two coefficients explicitly: which gives an asymmetric probability density |Ae ik1x +Be iϕ e −ik1x | 2 = 1+2AB cos 2 (k 1 x−ϕ/2). (S17) We conclude that the azimuthal angle on the Bloch sphere ϕ has to be zero, and that the pseudospin has to follow the great circle. The probability current J ± of each energy component ±E can be computed numerically from the expression: These currents can be plotted with respect to the phase difference (see Fig. S1). Furthermore, by plotting on the same graph ±J 0 sin(2φ), where J 0 (which is phase constant) is the maximum absolute value of J ± , one can notice that there is an excellent match between J ± and ±J 0 sin(2φ): This trend can be traced back to the definition of the current on the right interface. The probability current of each energy component measures the exchange between particles at +E and particles at −E. Thus, considering the positive energy component for instance, it follows: Regarding the expressions of these coefficients given in Eq. (S7), one can deduce that: Finally, we retrieve the expression: where J 0 has a small phase dependence, contrary to J 0 . This phase dependence comes from the dependence of the energy of the components of a bound state on the phase difference (via the norm of the reflection coefficients).

Evanescent states in the normal region
In the main text, the case E p > ∆ is considered because it leads to propagative states in the normal region for both positive and negative values of E, which is the most interesting case. However, the case with E p < ∆ is possible as well. Then, for incident particles with energies E > E p , the reflected part at the energy symmetric with respect to the pump detuning is evanescent. We have solved the reflection problem in this case and obtained a non-zero amplitude of the reflected evanescent wave. Our calculations show that SNS junctions with this type of states can exist as well. They present the same global behaviour as for the propagative states (see Fig. S2(a,b) for the two configurations with a different phase), with the wavefunction in the normal region being a linear combination of hyperbolic functions. However, the bands they form no longer cross. Indeed, the crossing of the bands in the main text occurred when an original state of the quantum well had the same energy as the bogolon image of another state. This is not possible when the original states are propagative and the images are evanescent, since they are always at the opposite sides of zero. Thus, the topologically protected self-amplified interface states discussed in the main text cannot be observed for these bands.

Non-topological configuration of Andreev bound states
To confirm that the band crossing and amplification observed in the main text are indeed due to the topology of the bands, we consider an alternative configuration where the bands do not exhibit any inversion of the topology with the variation of the parameter, and thus they anticross instead of crossing each other. We introduce an additional potential barrier at x = 0, splitting the trap into two parts with two distinct trapped states (left-and right-localized) having the same symme-try (s for the lowest state). We also introduce a potential step, responsible for the detuning of these two states. Instead of varying the width of the trap a, we vary the height of the step as a function of the second coordinate y. The total potential therefore reads where U 0 = 2 meV is the barrier height, σ = 0.6 µm is its width, U 1 = 0.2 meV is the characteristic step height, and l = 50 µm is the characteristic variation length of the step height.
The results of numerical simulations of this configuration are shown in Fig. S3. Panel (a) presents an example of the calculated spectrum of the Andreev bound states for a particular value of U 1 y/l = 0.15 meV (the total detuning between the left and right states is 0.3 meV). The original (s-type) states are clearly visible, as well as their bogolon images exhibiting p-symmetry due to the laser phase φ = π/2. Each of the four visible states belongs to a band (as a function of the synthetic variable φ). The Zak phases of the bands are calculated as in the main text. They are shown in Fig. S3(b), together with the energies of the band extrema at φ = π/2 plotted as a function of the step height. The two lowest bands, formed from the original s-symmetric states, have a zero Zak phase. Their symmetry is the same. Thus, when the step height changes sign and the detuning inversion leads to the state inversion (the lowest state changes localization from left to right), the topology of the system does not change. There are no topological reasons for the crossing of the bands, and indeed, it does not occur: their anticrossing is controlled by the tunneling t across the barrier in the center (controlled by its height U 0 ). The same concerns the two upper bands, sharing the same topology (Zak phase π, different from the two lowest bands). Finally, panel (c) confirms that no amplification due to a band crossing occurs in this case (since the crossing is actually avoided): as the step height changes with y, the detuning of the states with respect to the laser changes, and we observe the transfer of maximal intensity from left to right, but no signs of mascroscopically populated oscillating modes are visible. This confirms that the band crossing and the resulting mode amplification discussed in the main text are indeed of a topological origin.