A I. Dispersion of two-photon excitations in an infinite array

# Topological edge states of bound photon pairs

## Abstract

We predict the existence of interaction-driven edge states of bound two-photon quasiparticles in a dimer periodic array of nonlinear optical cavities. Energy spectrum of photon pairs is dramatically richer than in the noninteracting case or in a simple lattice, featuring collapse and revival of multiple edge and bulk modes as well as edge states in continuum. Despite the unexpected breakdown of the Zak phase technique and the edge mixing of internal and center-of-mass motion we link the edge state existence to the two-photon quantum walk graph connectivity, thus uncovering the topological nature of the many-body problem in complex lattices.

Nonlinear and many-body phenomena in condensed matter physics and optics are currently in the focus of research interest due to the wide range of opportunities including realization of strongly correlated photon gases, implementation of polariton superfluidity and formation of solitons and vortices [Carusotto and Ciuti, 2013]. One of such striking interaction-induced effects has recently been experimentally observed by Winkler et al. [Winkler et al., 2006]. In the presence of a repulsive interaction two bosons can form a bound pair propagating as a single quasiparticle in an optical lattice Winkler et al. (2006); Valiente and Petrosyan (2008). Such repulsively bound pairs further termed as doublons do not have direct analogues in traditional condensed matter systems being of fundamental interest for the many-body physics and quantum information.

Search for the doublon edge states [Pinto et al., 2009a, b; Longhi and Valle, 2013] has revealed their absence in a simple periodic lattice [Pinto et al., 2009a]. Tamm edge states of doublons may arise similarly to the single-particle case if the lattice has an edge defect [Longhi and Valle, 2013]. However, the existence of topological doublon edge states as well as the calculation of topological invariants for composite particles remain open problems so far.

Inspired by the recent advances in topological photonics [Lu et al., 2014; Khanikaev et al., 2013; Ma et al., 2015; Slobozhanyuk et al., 2015] and quantum optics Greentree et al. (2006); Hartmann et al. (2006), we investigate the edge states of doublons in a dimer lattice of identical optical cavities with a Kerr-type nonlinearity and two alternating tunneling constants, , see Fig. 1. This is a many-body generalization of the Su-Schrieffer-Heeger (SSH) model, considered a simplest example for the topological edge states of photons Schomerus (2013) and plasmons Ling et al. (2015); Cheng et al. (2015). Topological transitions in somewhat similar classical nonlinear systems have been recently predicted in Refs. Hadad et al. (2016); Solnyshkov et al. (2016). Ground state of the many-body SSH model has been studied in Ref. Grusdt et al. (2013). To the best of our knowledge, the two-particle excitations of the interacting SSH model are fully uncharted. Here, we present a rigorous numerical diagonalization for the finite system, accompanied by an exact analytical solution for the bulk two-photon states and by a condition for the edge states, based on the topology of the two-photon quantum walks in the bulk. While the quantum walks have been extensively studied Aharonov et al. (2001); Omar et al. (2006); Lahini et al. (2012); Manouchehri and Wang (2013) and even successfully used to describe the noninteracting topological model Kitagawa et al. (2010), we are not aware of their direct applications to characterize the edge states in interacting systems.

We employ the Bose-Hubbard-type Hamiltonian Essler et al. (2005); Grusdt et al. (2013)

 ^H=ω0∑m^nm+U∑m^nm(^nm−1)−J1∑m(^a†2m−1^a2m+^a†2m^a2m−1)−J2∑m(^a†2m^a2m+1+^a†2m+1^a2m), (1)

where and are the photon creation and annihilation operators for the -th cavity, is the photon number operator, are the tunneling constants, is the interaction strength and is the cavity eigenfrequency (). We assume fixed photon polarization.

Since the Hamiltonian Eq. (1) commutes with the operator , the total number of photons is conserved and the two-photon state can be searched as with . Substituting the wave function into the Schrödinger equation we obtain a linear system of equations to determine the coefficients and the energy . Thus, the one-dimensional two-photon problem is equivalent to a two-dimensional (2D) single-particle problem Longhi and Valle (2013). The corresponding 2D lattice is illustrated in Fig. 2(a) for and in Figs. 2(b,d) for . The and coordinates are just the coordinates of first and second photon. The links represent the tunneling amplitudes ; the energies for the diagonal sites with are equal to (two photons in the -th cavity), the off-diagonal sites have zero energies (space-separated photon pair). Next, we show how the existence of the doublon edge states follows from the connectivity of the two-photon tunneling pathways in this 2D lattice and verify the result by the full diagonalization.

We first reexamine the case , where the doublon edge states are absent Pinto et al. (2009a, b), in the limit of strong interaction . In the zeroth order approximation in all the doublon states have the energy , they are degenerate and the two photons are pinned to the same cavity [dotted diagonal in Fig. 2(a)]. The tunneling couples doublons at different cavities as well as shifts their energies. The resulting effective Hamiltonian for doublons has the general form

 ^Heff=N∑l=1(ε0+δl)^d†l^dl+tN−1∑l=1(^d†l^dl+1+H.c.), (2)

where and . The effective two-photon tunneling amplitude and the two-photon energy blueshift can be expanded in the perturbation series in . We resort to the nearest-neighbor coupling for doublons. In this case the tunneling constant can be calculated by means of the second order perturbation theory,  Pinto et al. (2009a). This corresponds to the two photons tunneling to the adjacent cavity directly one after another. The corresponding tunneling pathways are shown in the two-photon quantum walk graph of Fig. 2(a) by the red arrows. The loops linking the site to itself, i.e. the path in Fig. 2(a), contribute to the energy shifts .

The existence of edge states in the model Eq. (2) depends on the relation between and that, in turn, follows from the topology of the quantum walk graph. Two important identities hold for : and . The first identity stems from the translational symmetry and the fact that the first and last cavities have twice less neighbors. The condition means that the self-induced nonlinear blueshift of the site in the bulk can not exceed the energy shift induced by the left and right neighbors. This condition arises because of the mirror symmetry. Namely, the mirror reflection with respect to the vertical lines 1–3 or 2–4 in Fig. 2(a) maps each two-step path coupling nearest neighbors to the path coupling the site to itself, e.g. to . There exist two symmetry lines, so . In terms of the two-photon quantum walks topology the argument above means that the local vertex connectivity , corresponding to the eigenmode, is equal to . The vertex connectivity is by definition the minimal number of vertices that have to be cut from the graph to make it disconnected West (2001). The relevant two-photon tunneling pathways in Fig. 2(a) are schematically shown in Fig. 2(c) with the points 1–4 forming the minimal set of vertices to be cut. Given the conditions above, the Hamiltonian Eq. (2) describes a simple periodic lattice, where the first and the last sites are detuned by the energy from the middle ones. The detuning is exactly equal to the tunneling amplitude, . Hence, the condition for the edge states in the 1D tight-binding model with the edge defect [Sup. Mat. V] is not satisfied. The edge states do not exist for , although they can be enabled by an arbitrary weak perturbation Pinto et al. (2009b); Longhi and Valle (2013).

The above analysis of the two-photon quantum walks can be generalized to the much less trivial dimer lattice with . We now assume that , the number of cavities is even, , and the lattice ends with a strong tunneling link at both edges. In the perturbation scheme we first take into account the strong links and then the weak links . Due to the dimerization the correct doublon operators in the zeroth order in are the odd and even combinations . The effective doublon Hamiltonian still has the structure Eq. (2) with and being replaced by . The energies and the tunneling amplitude are modified. Since now each doublon mode occupies two cavities instead of one, the consistent derivation of and requires fourth-order perturbation theory [Supplemental Materials, Sec. IV], i.e. including all two-step and four-step pathways in Figs. 2(b,d). We find that contrary to the case, the condition no longer holds in general because the local vertex connectivity of the eigenmode can exceed 4. Indeed, for the even mode as illustrated in Fig. 2(d). The corresponding quantum walk graph [Fig. 2(e)] is inherently irreducible to the graph for in Fig. 2(c). Consequently, more pathways contribute to the blueshift than to the tunneling, and one has [Sup. Mat. IV]. Since the detuning of the bulk sites from the edge exceeds , the edge state does appear. As such, the condition can be used to predict the presence of the edge state.

The odd mode still has and no edge states for strong link termination. The contributions from all the pathways crossing the points and cancel each other due to the odd mirror symmetry, so these points are to be excluded from the graph, see Fig. 2(b). More detailed analysis, presented in the Supplementary Materials, confirms that the doublon edge state is formed for the even band but not for the odd one. Its localization length is given by and the energy is close to the lower edge of the even band . The crude argument for the emergence of the edge state only for the even band of doublons is that its blueshift due to the interaction is stronger in general than that for the odd mode. Hence, the energy detuning between the edge and bulk sites is larger and this facilitates the edge state formation. We stress that the two-photon state is formed at the edge with strong tunneling link, where the single-photon edge states are absent.

The analysis above has been performed in the limit of strong interaction and focused on the doublon edge states only. Now we will discuss the whole energy spectrum for an arbitrary interaction strength. In the infinite lattice two types of excitations are possible: (i) pairs of quasi-independent photons that move along the lattice almost without interaction and have the energies , where is the single-photon energy, being the real Bloch wave numbers varying from to and (ii) doublons with complex and of the form , , where is a real number describing the motion of photon pair as a whole, while a complex number describes the relative motion of photons confined to each other. As shown in the Supplementary Materials, the bulk doublon states can be sought in the form of modified Bethe anzatz that yields exact analytical equations for the doublon dispersion [Sup. Mat. I].

New types of two-photon excitations emerge in a finite array, namely, the single-photon edge states when one photon is localized at the edge of the lattice, while the other one moves along the lattice, and also the doublon edge states where two bound photons are localized at the edge. The single-photon edge states exist only at the edge with weak tunneling link and have the energy , where is a real wave number of the delocalized photon. Quite surprisingly, contrary to the case when , the internal and center-of-mass degrees of freedom of doublons are mixed at the edge [Sup. Mat. VI]. We have the dimerization brings massive difficulties into the application of the Bethe ansatz and instead we resorted to the numerical diagonalization of the Hamiltonian in the finite system.

The results for an array with cavities with and a given interaction strength are shown in Fig. 3. Panel (a) presents the bands of quasi-independent photons (green), one-photon edge states (blue), bulk doublons (red). The first Brillouin zone for doublons corresponds to the range , while for quasi-independent photons. Horizontal lines show doublon edge states localized at the edge with strong tunneling link (black) and weak link (purple). The panels (c-f) show the 2D color maps of the corresponding wavefunctions . Doublon localization at the edge with strong link for the uppermost band of even symmetry is in full agreement with our analysis above. We also observe doublon states localized at the opposite edge with the weak link that are facilitated by the interaction with the single-photon edge states.

A quite interesting phenomenon is the decay of both bulk and edge doublon states due to their interaction with the quasi-independent photons and the single-photon edge states, respectively. There exist four bulk doublon bands in total, two of which are unstable against the decay into a pair of quasi-independent photons. We refer to this situation as doublon collapse: no bound solution with complex can be found for a certain range of . For instance, the third from top doublon band in Fig. 3(a) is unstable and exists only for the wave vectors close to . Such effect is illustrated in Fig. 3(b) in a larger scale. Additional analysis of doublon collapse is provided in Supplemental Materials, Sec. II.

The stability regions for the bulk and edge states can be further traced by their energy dependence on the interaction strength, shown in Fig. 4(a). Color of the circles and squares depicting doublon edge states in Fig. 4(a) encodes their degree of localization defined as . The number of stable bulk doublon bands varies from 2 to 3 depending on the value of , see Fig. 4(b). The uppermost band with and the corresponding edge state, localized at the strong link edge, [squares in Fig. 4(a)], always remain stable. For weak nonlinearity this edge state exists in the continuum of quasi-independent photons [Fig. 4(c)]. We have numerically verified that it retains exponential localization in the continuum [Sup. Mat. III] and that it exists for an arbitrary ratio of the tunneling constants . The second from top doublon band, with does not have edge states. The states, localized at the weak link edge [circles in Fig. 4(a,b)], can collapse or revive as function of due to the interaction with one-photon edge states (blue bands).

To summarize, our condition for the edge states of bound photon pairs in terms of the local vertex connectivity of the quantum walks graph has been fully confirmed by the numerical calculation. The predicted edge state of topmost even photon pair band is localized at the edge with greater tunneling constant, where the single-photon edge states are absent. This state exists for all interaction strengths and for arbitrary unequal tunneling constants. Its stability against the interaction with the continuum of quasi-independent photons might be a sign of topological protection Zhen et al. (2014). It remains to be understood if any meaningful topological invariant aside from the local vertex connectivity can be assigned to the doublon bands. We are not aware of approaches that can handle the composite quasiparticle decay, inherently present in the problem. For instance, the Zak phase, typically used as an evidence of the topological character of the noninteracting Su-Schrieffer-Heeger model Shen (2012), is not informative here. Our calculations [Sup. Mat. VII] show that the Zak phase is equal to , if , and it is equal to , when . This result is the same for all the four doublon bands provided that they are stable and hence does not explain why only the state from the topmost band can be localized at the edge with greater tunneling constant. The breakdown of Bethe anzatz at the edge might further hinder the application of traditional bulk-boundary correspondence. On the contrary, the vertex connectivity of the quantum walks can be analyzed for various composite particles (e.g. photon triplets Pinto et al. (2009a)) in complex 1D lattices.

Thus, the simple appearance of the two-particle Su-Schrieffer-Heeger model is deceptive. It uncovers a wide spectrum of fundamental phenomena such as interaction-induced edge localization and decay of the bound bulk and edge quasiparticles into the weakly interacting ones. Our results may be useful for a whole range of quantum systems described by Bose-Hubbard model thus paving a way to nonlinear topological physics. The simplest demonstration might be provided by waveguide lattices where the quantum walks of noninteracting Peruzzo et al. (2010); Schreiber et al. (2012); Solntsev et al. (2014); Mittal et al. (2016) and interacting photons  Corrielli et al. (2013) can be emulated classically.

The authors acknowledge valuable discussions with D. Yudin, G. Zhilin, M. Hafezi, I.S. Sinev, A.K. Samusev, A.V. Poshakinskiy, M.M. Glazov, A.A. Sukhorukov and Yu.S. Kivshar. This work was supported by the “Dynasty” foundation, investigation of bulk doublon dispersion was supported by the Russian Foundation for Basic Research (Grant No. 15-32-20866), calculation of Zak phase was supported by the Russian Science Foundation (Grant No. 16-19-10538). ANP was supported by the Russian President Grant No. MK-8500.2016.2.

Supplementary Materials

## Appendix A I. Dispersion of two-photon excitations in an infinite array

Since the system Hamiltonian [Eq. (1) in the main text] commutes with the operator , the total number of photons is conserved and the two-photon wave function can be represented as

 |ψ⟩=∑n√2βnn|2n⟩+∑m≠nβmn|1m1n⟩, (S1)

where , and denotes the stationary state when one photon is located in -th cavity, and the other one is located in the -th cavity. Without a loss of generality we assume that .

Substituting this wave function and the system Hamiltonian into the eigenvalue equation

 ^H|ψ⟩=(ε+2ω0)|ψ⟩ (S2)

we obtain the system of linear equations with respect to the expansion coefficients :

 (ε−2U)β2m,2m=−2J1β2m−1,2m−2J2β2m,2m+1, (S3) (ε−2U)β2m+1,2m+1=−2J1β2m+1,2m+2−2J2β2m,2m+1, εβ2m,2m+2n=−J1(β2m−1,2m+2n+β2m,2m+2n−1)−J2(β2m+1,2m+2n+β2m,2m+2n+1), εβ2m,2m+2n−1=−J1(β2m−1,2m+2n−1+β2m,2m+2n)−J2(β2m+1,2m+2n−1+β2m,2m+2n−2), εβ2m+1,2m+2n=−J1(β2m+2,2m+2n+β2m+1,2m+2n−1)−J2(β2m,2m+2n+β2m+1,2m+2n+1), εβ2m−1,2m−1+2n=−J1(β2m,2m−1+2n+β2m−1,2m+2n)−J2(β2m−2,2m−1+2n+β2m−1,2m−2+2n).

In Eqs. (S3) one has . The system of equations can be interpreted as a two-dimensional problem for a single particle. This equivalent two-dimensional problem is illustrated in Fig. S1.

The expansion coefficients can be searched in the form of the standard Bethe ansatz Essler et al. (2005):

 βmn=Cj(m,n)eik1m+ik2n, (S4)

where for the even-even, even-odd, odd-even and odd-odd pairs of indices , respectively. We assume that in Eq. (S4). In the case when , the relation is used.

Analysis of Eqs. (S3) with the ansatz Eq. (S4) yields the law of dispersion of two-photon excitations

 ε4−4ε2[J21+J22+2J1J2cos(k1+k2)cos(k1−k2)]+16J21J22sin2(k1+k2)sin2(k1−k2)=0. (S5)

Importantly, Eq. (S5) describes all considered types of two-photon excitations: bound pairs and quasi-independent pairs, bulk pairs and edge pairs. Equation (S5) can be rearranged in the equivalent form

 ε=±√J21+J22+2J1J2cos2k1±√J21+J22+2J1J2cos2k2. (S6)

In the case of real and this equation describes the states of quasi-independent photons. The energy of such state is represented as a sum of single-photon energies. Both Bloch wave numbers and are varying in the range from to .

In addition to the quasi-independent photon states, the photon-photon interactions determined by the term of the Bose-Hubbard model give rise to the other type of two-photon states, namely, bound photon pairs (doublons) [Winkler et al., 2006]. Such excitations are characterized by complex and of the form , , where is a real number describing the motion of the photon pair as a whole, while a complex number describes the relative motion of photons. The latter number is assumed to have a positive imaginary part, , that captures the effect of photons binding. We note that the first Brillouin zone for doublon is one-dimensional and spans from to .

To obtain the ansatz in the form applicable to doublons, we rewrite Eq. (S5) as follows:

 2J21J22sin2k⋅ξ2+ε2J1J2cosk⋅ξ−18[ε4−4ε2(J21+J22)+16J21J22sin2k]=0, (S7)

where . Importantly, for the fixed values of doublon energy and Bloch wave number Eq. (S7) is a quadratic equation with respect to . This means, that in the general case there exist two values of with the positive imaginary part corresponding to the given Bloch wave number and energy ; these two roots are denoted further as and . Therefore, to describe the dispersion of doublons, we use a modified Bethe ansatz:

 βmn=eik/2(m+n)[Cj(m,n)eiϰ/2(n−m)+¯Cj(m,n)ei¯ϰ/2(n−m)] (S8)

with for even-even, even-odd, odd-even and odd-odd pairs of indices and . This ansatz yields the following system of linear equations:

 εC1+(J1e−ik2+J2eik2)C2+(J1e−ik1+J2eik1)C3=0, (S9) (J1eik2+J2e−ik2)C1+εC2+(J1e−ik1+J2eik1)C4=0, (S10) (J1eik1+J2e−ik1)C1+εC3+(J1e−ik2+J2eik2)C4=0, (S11) (J1eik1+J2e−ik1)C2+(J1eik2+J2e−ik2)C3+εC4=0, (S12) ε¯C1+(J1e−i¯k2+J2ei¯k2)¯C2+(J1e−i¯k1+J2ei¯k1)¯C3=0, (S13) (J1ei¯k2+J2e−i¯k2)¯C1+ε¯C2+(J1e−i¯k1+J2ei¯k1)¯C4=0, (S14) (J1ei¯k1+J2e−i¯k1)¯C1+ε¯C3+(J1e−i¯k2+J2ei¯k2)¯C4=0, (S15) (J1ei¯k1+J2e−i¯k1)¯C2+(J1ei¯k2+J2e−i¯k2)¯C3+ε¯C4=0, (S16) Missing or unrecognized delimiter for \left (S17) Missing or unrecognized delimiter for \left (S18)

Here, and . Note that the determinant of Eqs. (S9)-(S12) as well as Eqs. (S13)-(S16) is zero due to the same dispersion equation Eq. (S7). Thus, the calculation of the doublon dispersion can be accomplished as follows:

• Express and as functions of energy and wave number from Eq. (S7).

• Express , and via , and from Eqs. (S9)-(S12).

• Express , and via , and from Eqs. (S13)-(S16).

• Obtain a system of homogeneous linear equations with the unknown variables and from Eqs. (S17), (S18). The coefficients in the obtained system depend on the doublon energy and the wave number .

• Doublon dispersion equation is deduced setting to zero the determinant of the derived system of equations.

## Appendix B II. Doublon collapse

We analyze the dispersion equation Eq. (S7) in the special case denoting by . It is straightforward to show that

 Missing or unrecognized delimiter for \right (S19)

Photons will be bound together if Im , i.e. . In the opposite scenario , is purely real, photons are no longer confined to each other and doublon thus decays into the pair of quasi-independent photons. Further we refer to this phenomenon as doublon collapse. It can be shown that the condition for doublon collapse reads

 2|J2−J1|<|ε0|<2|J2+J1|. (S20)

On the other hand, the energy intervals , and for correspond to the energy bands of quasi-independent photons. Thus, we conclude that if the doublon energy falls into the upper or lower energy band of quasi-independent photons, the doublon becomes unstable and collapses. Besides that, our calculations show that the doublon is always stable for the values of close to . The phenomenon of doublon collapse is illustrated in more detail in Fig. S2 where two panels correspond to different strengths of interaction.

## Appendix C III. Two-photon states in a semi-infinite array

New types of two-photon excitations emerge, if the array has an edge. In this section, we analyze various types of edge states in a semi-infinite geometry, when the structure starts from the cavity with the number . In such situation system of equations Eq. (S3) should be supplemented by the boundary conditions of the form:

 (ε−2U)β11=−2J1β12, (S21) εβ1,2n=−J1(β1,2n−1+β2,2n)−J2β1,2n+1, (S22) εβ1,2n+1=−J1(β2,2n+1+β1,2n+2)−J2β1,2n, (S23)

where .

First we consider a situation when one photon is localized at the edge, while the other one moves freely along the lattice. We assume ansatz Eq. (S4) for coefficients. Then Eqs. (S3) are compatible with Eqs. (S21)–(S23) only under the condition . Simple analysis reveals that edge states exist if . The energy of the edge state and its localization parameter read

 ε=±√J21+J22+2J1J2cos2k2, (S24) e2ik1=−J1/J2. (S25)

In fact, the edge state described by Eqs. (S24)-(S25) is a direct analog of that existing in a one-dimensional single-particle Su-Schrieffer-Heeger model.

Another type of edge modes is represented by the edge states of bound photon pairs. However, their analytical description is extremely cumbersome (if possible at all, see Sec. VI), and therefore we resorted to the full numerical diagonalization of the finite system Hamiltonian. In particular, we demonstrate that doublon edge states can arise in continuum. The results of numerical diagonalization of the system Hamiltonian are presented in Fig. S3.

In order to provide more insight into the properties of different two-photon states, we plot the dispersion of various types of two-photon excitations for a fixed tunneling constants ratio and a fixed interaction strength (Fig. S4a, cf. with Fig. 4 in the main text). Such values of parameters correspond to the partial overlap of quasi-independent photons energy band with the energy band of single-photon edge states. Therefore, some types of doublon edge states present for larger tunneling constants ratio, e.g. for , disappear. This tendency is further illustrated by Fig. S4b. Figure S4b shows also that only two from four doublon energy bands are strongly affected by the nonlinearity, exhibiting almost linear growth with the interaction strength . The reason for such behavior is that photons are predominantly localized in the same cavity for two upper doublon bands, while for two lower doublon bands photons are mainly localized in the neighboring cavities.

## Appendix D IV. Effective Hamiltonian for two upper doublon bands

As it is pointed out in the previous section, energies of the two from four doublon bands depend almost linearly on . Thus, in the limit of strong interaction these two bands are well separated from the remaining states (see Fig. S4b) and the two photons are well confined to each other. Therefore, it is possible to develop an effective one-dimensional model describing the properties of the upper doublon bands.

Examining the calculated probability distributions for the two upper doublon bands, we notice that the coefficients with are always negligible. Therefore, for this special case the system of equations Eqs. (S3) together with the boundary conditions Eqs. (S21)-(S23) can be truncated. Excluding from these truncated equations , and leaving only the terms , we obtain:

 εα1=(2U+j1+τ)α1+(j1+2τ)α2+τα3, (S26) εα2=(2U+j1+j2+5τ)α2+(j1+2τ)α1+(j2+4τ)α3+τα4, εα3=(2U+j1+j2+6τ)α3+(j2+4τ)α2+(j1+4τ)α4+τα1+τα5, εα4=(2U+j1+j2+6τ)α4+(j1+4τ)α3+(j2+4τ)α5+τα2+τα6, …

where and . Equations (S26) describe a one-dimensional lattice with the next-nearest neighbor hopping.

In order to further simplify the analysis we consider the limit of strong interaction and strongly different tunneling constants . In such case, in the zero-order approximation the array consists of isolated dimers and its eigenstates are

 u(±)=(1/√2±1/√2). (S27)

We rewrite Eqs. (S26) in terms of the vectors and project these equations onto the eigenvectors . As a result, we derive the following equations for the projections :

 εψ1 =(ε(±)0+δ±/2)ψ1+t±ψ2, (S28) εψm =(ε(±)0+δ±)ψm+t±(ψm−1+ψm+1), (S29)

where , and . The “” sign choice corresponds to the antisymmetric eigenstates and describes the second from top doublon energy band. On the contrary, the “” sign choice corresponds to the symmetric eigenstates and describes the upper doublon energy band. For both signs, the equations Eq. (S28) and Eq. (S29) formally correspond to the one-dimensional array with the detuned resonator at the edge. The magnitude of detuning is equal to . Thus, the well-known condition for edge states emergence reads (see Sec. V). From the derived results it is straightforward to show that and . Using the correspondence with the one-dimensional lattice, it is also easy to estimate localization length of doublon edge state associated with the upper band:

 l=1ln(δ+2|t+|)≈2|t+|δ+−2|t+|≈(UJ1)2. (S30)

These conclusions are in perfect agreement with the general topological arguments provided in the article main text.

## Appendix E V. Edge states in a simple lattice with the detuned edge

In the previous section we have reduced the problem for both symmetric and antisymmetric doublon modes to the following one:

 εψ1=Δψ1+tψ2, (S31) Missing dimension or its units for \mspace (S32)

Equations (S31),(S32) correspond formally to the case of one-dimensional semi-infinite lattice with the edge cavity detuning equal to . We search edge states in the system assuming with having positive imaginary part. It is straightforward to show that . Since due to nonzero imaginary part of , the condition for edge state existence reads . Localization length of edge state can be calculated as

 l=1Imk=1ln|t/Δ|. (S33)

## Appendix F VI. Bethe ansatz breakdown

In this section we discuss a general ansatz for that allows one to investigate analytically all types of two-photon excitations supported by a finite cavity array. Due to the property , it is sufficient to consider only the coefficients with .

To get a general idea of constructing ansatz for , we consider first a situation with , previously studied in the literature in Refs. [Pinto et al., 2009a, b; Longhi and Valle, 2013]. We map the two-photon problem onto the corresponding two-dimensional single-particle problem. Two-photon state is represented as a superposition of the wave and various scattered waves which arise from the scattering at the diagonal and at the edge of the two-dimensional sample. Such “diagonal” and “edge” scattering is described by the operators (Fig. S5).

 d(k1,k2)=(k2,k1), (S34) e(k1,k2)=(−k1,k2). (S35)

We note further that the elements and describing two types of scattering generate a group of the eighth order, i.e. ansatz for will contain eight terms

 βnm=C(1)j(m,n)eik1n+ik2m+C(2)j(m,n)e−ik1n+ik2m+C(3)j(m,n)eik1n−ik2m+C(4)j(m,n)e−ik1n−ik2m++C(5)j(m,n)eik2n+ik1m+C(6)j(m,n)e−ik2n+ik1m+C(7)j(m,n)eik2n−ik1m+C(8)j(m,n)e−ik2n−ik1m, (S36)

It is this ansatz that was employed in the analysis of Tamm-Hubbard states in an array with equal tunneling constants  [Longhi and Valle, 2013].

However, the situation becomes much more complicated if . Besides the elements and describing scattering at the diagonal and the edges, one also needs to take into account the element that acts on wave vector as follows:

 x(k1,k2)=(¯k1,¯k2), (S37) Missing or unrecognized delimiter for \left (S38) Missing or unrecognized delimiter for \left (S39)

Indeed, this scattering is responsible for doublon formation. It turns out that the group generated by the elements , and is infinite in the general case. Thus, it is problematic to find the ansatz analogous to Eq. S36 and valid for the description of all two-photon states in a finite array. For that reason, we opted to study doublon edge states using full numerical diagonalization of the finite system Hamiltonian.

## Appendix G VII. Calculation of the Zak phase for doublons

In this section we calculate the Zak phase for doublons propagating in an infinite lattice. To this end, we extract the periodic part from the entire doublon wave function calculated in Sec. I. The choice of the two-photon states comprising the periodic part is illustrated in Fig. S1, so reads:

 |uk⟩={√2β00|20⟩+2∞∑n=1β−2n,2n|1−2n12n⟩+2∞∑n=1β2−2n,2n|12−2n12n⟩}+{2∞∑n=1β2−2n,2n−1|12−2n12n−1⟩+∞∑n=1β−2n,2n−1|1−2n12n−1⟩+∞∑n=1β2−2n,2n+1|12−2n12n+1⟩}+{2∞∑n=1β1−2n,2n|11−2n12n⟩+∞∑n=1β1−2n,2n−2|11−2n12n−2⟩+∞∑n=1β3−2n,2n|13−2n12n⟩}+{√2β11|21⟩+2∞∑n=1β1−2n,2n+1|11−2n12n+1⟩+2∞∑n=1β1−2n,2n−1|11−2n12n−1⟩}. (S40)

Note that the chosen unit cell possesses inversion symmetry. Therefore, we expect the Zak phase to be quantized [Zak, 1989]. Using the modified Bethe ansatz Eq. (S8) for the coefficients , we obtain:

 |uk⟩=4∑j=1∣∣uj⟩, (S41) ∣∣uj⟩=Cj∣∣vj⟩+¯Cj∣∣¯vj⟩, (S42) |v1⟩=√2|20⟩+2∞∑n=1e2iϰn|1−2n12n⟩+2eik∞∑n=1eiϰ(2n−1)|12−2n12n⟩, (S43) |v2⟩=2eik/2∞∑n=1eiϰ(2n−3/2)|12−2n12n−1⟩+ e−ik/2∞∑n=1eiϰ(2n−1/2)|1−2n12n−1⟩+e3ik/2∞∑n=1eiϰ(2n−1/2)|12−2n12n+1⟩, (S44)
 |v3⟩=2eik/2∞∑n=1eiϰ(2n−1/2)|11−2n12n⟩+ e−ik/