Topology and interactions in the photonic Creutz and Creutz-Hubbard ladders

The latest advances in the field of photonics have enabled the simulation of an increasing number of quantum models in photonic systems, turning them into an important tool for realizing exotic quantum phenomena. In this paper we suggest different ways in which these systems can be used to study the interplay between flat band dynamics, topology and interactions in a well-known quasi-1D topological insulator: the Creutz ladder. Firstly, a simple experimental protocol is proposed to observe the Aharonov-Bohm localization in the noninteracting system, and the different experimental setups that might be used for this are reviewed. We then consider the inclusion of a repulsive Hubbard-type interaction term, which can give rise to repulsively bound pairs termed doublons. The dynamics of these quasiparticles are studied for different points of the phase diagram, including a regime in which pairs are localized and particles are free to move. Finally, a scheme for the photonic implementation of a two-particle bosonic Creutz-Hubbard model is presented.


I. INTRODUCTION
Recent years have seen vast advances in the field of quantum simulation, an approach first proposed by Feynman, [1] in which the dynamics of some complicated quantum system is simulated by another, highlycontrollable, quantum system. Analog, purpose-built, non-universal quantum simulators using photonic systems and cold atoms in optical lattices have allowed experimentalists to physically realize a wealth of exotic quantum phases and phenomena that were previously only theoretical proposals. In principle, this technique can provide an enormous speed advantage over the conventional approach of simulating such systems using digital (classical) computers, a feature known as "quantum supremacy". [2] Its analog nature eliminates the need for error correction protocols or high-fidelity quantum gates, making their implementation easier than that of digital quantum computers. [3] It's been estimated that a few tens of controllable qubits would suffice for achieving a significant quantum advantage over a classical computer. [4] A particularly important application of quantum simulation is the investigation of topological insulators; materials which are insulating in the bulk, but possess conductive surface states protected by topological order. [5] The simplest examples of this form of material occur in one-dimensional, [6] or quasi-one-dimensional systems. A notable example of the latter type of system is the Creutz ladder, [7] in which particles (bosonic or fermionic) move in a ladder-like lattice in the presence of a magnetic field, as shown in Figure 1. For some special values of this magnetic field, and in the absence of the runglike links, the energy bands of the Creutz ladder become flat. This means that the group velocity of the particles, v g = ∂ k ω(k), vanishes, and that free propagation of single particles is thus impossible. The mechanism behind this localization is the Aharonov-Bohm (AB) effect, the phase acquired by the wavefunction of a charged particle by moving through a region with a non-zero magnetic vector potential. This effect is also seen in some bipartite lattices such as the T 3 lattice [8] and the rhombus chain, [9] and it has been termed "Aharonov-Bohm caging". This localization effect has been observed experimentally in superconducting wire networks, [10] metal wire lattices, [11] ultracold atom lattices [12] and, most recently, in a photonic system. [13] As the single-particle kinetic energy is zero in a flat band, even extremely weak interactions can have an enormous effect on the system's properties. Introducing a contact interaction of Hubbard-type lifts the degeneracy of the flat bands and breaks AB-caging, significantly complicating the system's properties and producing a rich filling-factor-dependent phase space which can include a quantum pair liquid, a Wigner solid, a condensate, and a supersolid phase. [14,15] For strong interactions, particles can form repulsively-bound pairs called "doublons" which can themselves experience AB-caging at a different value of the magnetic flux. [16] This rich phenomenology makes the Creutz ladder an excellent arena to study the effects and interplay of topology, AB localization and interactions.
Photonic quantum simulators use the properties of quantum or classical states of light to analyze the behavior of different quantum systems. In the last few years, an increasing collection of topological materials have been studied in photonic setups, [17] both using quantum [18][19][20] and classical light. [21][22][23][24] The latter type of simulators are easier to implement but can only simulate first-quantized dynamics, and thus the realization of systems of interacting particles has been traditionally challenging. Nevertheless, interesting advances have been made in the field of nonlinear classical photonics in general, and in strongly correlated systems in particular. Especially relevant for our work are the nonlinear photonic setups which also include topological [25,26] or AB localization effects. [27] Additionally, a quantum topolectrical system with interactions has been recently realized experimentally. [28] For our current purpose, we will focus on photonic waveguide lattices along which classical light propagates in the paraxial regime, where the direction along the waveguide takes the role of the time coordinate in the quantum system. Recently, we have seen the first steps in the use of photonic waveguide lattice simulators for the study of both strongly correlated systems [29][30][31][32][33] and topological insulators [34][35][36]. We will propose a setup to include both interacting and topological effects.
In this paper, we first study the effects of AB interference on one-particle states in the Creutz ladder, and show how they can be observed experimentally in a photonic waveguide lattice. We then proceed to investigate how interactions modify these results, by introducing the minimal model which incorporates interparticle interactions: the two-particle Creutz-Hubbard ladder. Finally, we propose a way to implement the latter model in a photonic waveguide system, focusing on the geometry of the lattice and the necessary use of synthetic dimensions. [37] II. CREUTZ LADDER The Creutz ladder (Figure 1a) consists of two chains of sites connected with horizontal, vertical and diagonal hopping amplitudes. Additionally, a magnetic field is applied perpendicularly to the plane of the ladder. This causes some hopping amplitudes to acquire a so-called Peierls phase and become complex. We choose the gauge in which the horizontal hoppings cause a change of phase of e ±iφ/2 in the wavefunction, and the rest are real. The sign of the phase change is positive if the particle moves to the right (left) on the upper (lower) leg of the ladder, and negative if it moves in the opposite direction. For m = 0 (no "rung" hopping terms), the lattice is bipartite, given that all first neighbours of sites with an odd jcoordinate correspond to an even j, and vice versa.
We will first consider the single-particle case, and go on to consider interactions in Sec. III. Let c ( †) j,α be the destruction (creation) operator for a boson in site |j, α of the ladder, where j = 1, ..., L labels the rungs along the ladder and α = A, B labels the two legs. The Creutz Hamiltonian then takes the form: where the horizontal hopping terms are J α = Je iσαφ/2 , with σ A/B = ±1 and A = B and vice versa. The model is thus specified by two hopping amplitudes (J and m), and the magnetic flux, represented by φ.
FIG. 1. a) Geometry of the Creutz ladder. The coordinates j, α label the sites, and L is the length of the ladder. J and m are the corresponding hopping amplitudes. When a magnetic flux of Φ threads each plaquette, the complex horizontal amplitudes cause a phase change in the wavefunction of ±φ/2 = ±πΦ/Φ0. b) Topological phase diagram for the Creutz ladder. The phase φ is proportional to the magnetic flux through a plaquette and m is the vertical hopping amplitude. The left and right edges of the diagram represent the same systems. Two trivial (ν = 0) and two nontrivial (ν = ±1) topological phases can be distinguished. c) Image of the function (dx(k), dz(k)) [Equation (3)] along the Brillouin zone for some of the points in the phase diagram. Its winding number around the origin corresponds to the topological invariant ν of the phase the point belongs to. For O and T, ν is not well defined, and the image of the function passes through the origin.
The physical magnetic flux that threads each plaquette, Φ, is proportional to the total phase φ acquired by a particle moving along its border in a clockwise sense: φ = 2π Φ Φ0 = qΦ , where Φ 0 is the magnetic flux quantum, and q is the electric charge of a particle. In the m = 0 case, there are no rung-like bonds and plaquette-like paths cannot be defined. However, any 4-step paralelogram-shaped or triangular-shaped path is threaded by the same flux as a plaquette, and can be used instead (compare Figure 1a and Figure 2c).
Since the unit cell of the system is composed of two sites, its band Hamiltonian is a two-dimensional Hermitian matrix. Consequently, in momentum space it can be written in terms of the Pauli matrices as H(k) = d 0 1 + d · σ, where: In the abstract space formed byd = (d z , d x , d 0 ), the d 0 axis corresponds to the possible band degeneracies (i.e. gap closings) in the bands of the Creutz Hamiltonian. This means that all Hamiltonians that satisfy d x (k 0 ) = d z (k 0 ) = 0 for some k 0 in the first Brillouin zone (BZ) are metallic, with their two bands touching at all possible values of k 0 . For the rest of the systems, the winding number of their directed pathd(k) (which is closed because of the periodicity of the BZ) around the d 0 axis determines the topological invariant of the Hamiltonian, of type Z, which serves to distinguish the different topological phases of the system. In Figure 1c, all paths have been projected down to the (d z (k), d x (k)) plane for simplicity, and thus the topological invariant corresponds to its winding number around the origin.
The relevant symmetries of the system, which protect the topological phases, [38,39] are the spatial inversion symmetry I : H(k) → σ x H(−k) (for all values of φ), and the chiral symmetry X = σ y for φ = ±π. Some authors consider the π-flux Creutz ladder to be in the symmetry class BDI, [40] arguing that, in addition to the chiral symmetry, there are particle-hole and time reversal symmetries that are independently conserved [15,[41][42][43]. Other authors classify it as AIII, claiming that the magnetic flux breaks time reversal symmetry. [38,39] For other values of φ, however, only inversion symmetry protects the topological phases, and is not one of the traditional symmetries considered in the classification of topological insulators. [40] The topological systems in which it plays a significant role cannot be sorted into the usual symmetry classes, and do not usually exhibit edge states (something unusual for a topological insulator). [44] The Creutz ladder is an unusual example, as it does present edge states protected by I, making it an excellent tool to probe these exotic symmetry classes. It is also an outlier in the sense that all 1D topological insulators with a Z-type invariant in the traditional classes require chiral symmetry, while it is absent in the Creutz ladder for most values of φ.
We show the topological phase diagram of the system as a function of the phase induced by the magnetic flux, φ, and the vertical hopping amplitude, m, in Figure 1b. [38] Four topological phases can be distinguished. The FIG. 2. a) AB caging. In the system with φ = π, the time evolution of the state |j, α never allows the particle to leave the AB cage centered on that site (shaded in orange). The numbers inside each site and the color code used show the phase acquired by the wavefunction when tunneling to the second neighbors of j, α along the two simplest (i.e. two-step) paths. Each pair of paths interfere destructively, and the particle never leaves the cage. The color code used to represent the phases of the wavefunction is specified by the color wheel on the left. This code is used throughout this work to represent complex quantities. b) Rhombus chain. A magnetic flux represented by the phase φ threads each rhombus. The AB cage of a spinal site is shaded in orange. c) A Creutz ladder with m = 0 can be obtained by identifying the corresponding top and bottom sites of two rhombus chains, represented in red and blue. The purple sites mark the sites where the identification has taken place. The AB cages for this lattice are inherited from the rhombus chains that form it. d) Numerical simulation of the time evolution for state |3, A (see Figure  1a). Half a period is shown. The wavefunction exhibits the oscillatory behaviour predicted in (6). The color code in a) is also used here.
system is topologically trivial for |m| > 2J, with a vanishing winding number, and thus a Zak phase of zero. For |m| < 2J, two different topological phases exist, with winding numbers of ν = ±1, depending on the value of φ. Both phases are associated with a Zak phase of Z = π, given that this quantity is defined modulo 2π (see Appendix A). The boundaries of these regions correspond to metallic systems, in which the two bands of the model touch. The phase diagram is periodic in φ, with period 4π. This means that the left and right sides of Figure 1b represent the same points in phase space. As a simple consistency check, we calculated the Zak phase numerically for different points of the phase space using Equation A3.

A. Aharonov-Bohm localization and topological edge states
The Creutz model with φ = ±π, m = 0 is of particular interest because it presents non-dispersive energy bands. Such flat bands have a constant dispersion relation E ± (k) = ±2J, and thus a vanishing group velocity, meaning that particles stay localized in space. This effect is caused by interference arising from Aharonov-Bohm phases, and has been termed AB caging. This effect is most clearly seen in the rhombus chain lattice, shown in Figure 2b). Consider a particle initially localized on one of the central sites with a coordination number of four, which we can term a spinal site. It will propagate to a neighboring spinal site in two successive tunneling processes, and two different routes are available for this; either via the upper side of the chain, or the lower side. Due to the Peierls phases acquired along each path, there is a phase difference of φ between these paths, and consequently if φ = π, the two paths interfere destructively. In that case the occupation of the neighboring spinal sites is frozen at zero, and the particle only undergoes a breathing motion from the initial site to its four neighbours. These five sites constitute the particle's "cage".
Being a quantum interference effect, this phenomenon requires no disorder, and indeed, has been shown to be robust to moderate levels of it. [45] As well as the rhombus chain lattice, [9,16,46] this caging occurs in a variety of other flat-band bipartite lattices, including the T 3 or dice lattice [8] and the Lieb lattice. [47] The most direct way to observe AB caging is to measure the ocupation numbers for all sites in the lattice during the time evolution, and check that a particle initially localized in a single site only visits a small number of them (its cage). This approach is convenient for waveguide lattice setups, given that the initial state requires only pumping light into a single site. Experimental observation of AB caging has indeed been recently observed in this way for a waveguide rhombus chain using classical light. [13] The Creutz ladder differs from the rhombus chain, in that every site has the same coordination number. Their geometries, however, are closely related. A Creutz ladder with no rungs (m = 0) can be obtained by identifying the top and bottom sites of two rhombus chains (see Figure 2c). Accordingly we would expect the resultant Creutz ladder to inherit the AB cages from the constituent rhombus chains. This expectation can be verified explicitly by analyzing the Creutz ladder using a basis of Wannier states |j± with support on the four sites (j, A; j, B; j + 1, A; j + 1, B) that form each one of the plaquettes in the ladder, [15,41,45] defined by the following expression: Their energies in the flat band limit are E ± = ±2J. Any localized state with no support on the ends of the system can be expressed as a superposition of these plaquette states. In particular, Thus, in this limit the time evolution of state |j, α will be periodic, with a frequency equal to the absolute value of the energy of the bands. A similar analysis can be done for the eigenstates of the rhombus chain model [9], although it is less intuitive than the discussion in terms of AB phases we give above. In the Creutz laddder, using the definition of the Wannier basis [Equation (4)] in the time-dependent state yields: which describes an oscillatory behaviour between the original site and its first neighbours. Here, J is the diagonal hopping amplitude and t is time. As expected, this behaviour is reminiscent of the one found in the rhombus chain. [9] The initial site and its first neighbours constitute the AB cage of the particle. Site j, α, in the same rung of the ladder as the initial site, is one of its second neighbours if m = 0 and, as such, is not included in the cage. An exact numerical simulation was carried out for a system with L = 6. The time step used was ∆t = 0.001, where time is measured in units of J −1 (taking = 1). The results are shown in Figure  2d. The obtained half-period was T 1 /2 = 1.572. This result agrees with the analytical calculation in Equation (6), where T 1 = 2π/2J = π. A similar analysis shows that the cage for a particle localized at a single end site (e.g. 1, A, see Figure 1a) is constituted by the four sites closer to that end of the ladder (for 1, A, these are 1, A; 1, B; 2, A; 2, B).
In the case with m = 0, the connectivity of the ladder changes, preventing it from being locally bipartite (the rungless ladder is not globally bipartite with periodic boundary conditions if L is odd) and allowing new paths that do not interfere destructively and release the particle from its cage. This prevents AB caging from taking place.
Another interesting feature of the Creutz model is its edge states, only present in the topological phases and protected by the aforementioned symmetries. A simple yet partially unexplored question so far is the appearance of these edge states for different values of the parameters of the system. The form of the edge states depends drastically on the presence or absence of the vertical hopping  Figure 1b) for a ladder of length 12, obtained using numerical exact diagonalization on the Creutz Hamiltonian. In systems P and Q, AB interference keeps the particle localized in the two end sites. Diagrams similar to Figure 2a are included (left). In system S, localization is exponential. Bar graphs have been included in case S to illustrate the exponential decay of the wavefunctions. In all diagrams, the color of the sites and bars correspond to the phase and modulus of the wavefunction, following the color code specified in Figure 2a. amplitude. For m = 0, AB interference localizes the particle entirely in the first or last sites of the ladder, as was already pointed out by Creutz for the π-flux case. [7] As an example, in the states |L P and |L Q , the paths |1A → |2A and |1B → |2A cancel each other out (see Figures 3a,b). Any other, more complicated, path to site |2A from one of the end sites |1α gets canceled out with the corresponding one from site |1α . All edge states for m = 0 obey the formula: In particular, the edge states of the systems corresponding to the points P and Q in the topological phase diagram (Figure 1b) have the expression: This is confirmed by the numerical results shown in Figure 3a,b.
For all values of φ, the particle is AB-localized in only two sites as long as m = 0, a remarkable result given that the AB caging for single-particle bulk states only occurs for the flat-band systems m = 0, φ = ±π. At this point, it is important to have in mind that the details of the phase differences of the wavefunction between the different sites are gauge dependent, but the physical conclusions obtained (such as localization) are not. Adding rungs to the ladder creates many new paths with amplitudes that do not vanish when added together. Similarly to the case of AB caging, these changes cause the wavefunction to leak out of the end sites, and present an exponential profile instead, as is frequently seen in physical boundary states. A bar graph illustrates the exponential nature of the states in Figure 3c. The edge states represented in Figure 3a-c correspond to the systems labeled P, Q and S in the phase diagram from Figure 1b, and were obtained diagonalizing the Creutz Hamiltonian numerically for a length of L = 12.

B. Photonic implementation
A direct implementation of a physical Creutz-like structure will encounter the problem of designing the diagonal crosslinks. Some efforts have been made in this direction. For example, an ultracold atom setup was proposed in which two overlapping zigzag optical lattices could be used to realize the lattice geometry. [48] An alternative way to overcome this challenge, which we will explore in this work, is to use a synthetic dimension. Synthetic dimensions are widely used in photonic materials [49] and ultracold atom systems [50,51] to increase the dimensionality of the lattice. A synthetic dimension employs a non-spatial discrete degree of freedom of the system (often the frequency mode number in photonic arrays and the different atomic hyperfine levels in ultracold atoms) to increase the number of dimensions by one. In this way the physics of a system with d + 1 spatial dimensions can be simulated by a d-dimensional system with one additional synthetic dimension.
The simplest way to implement the Creutz ladder in this manner in a photonic system would be to use a two-level degree of freedom to simulate the vertical dimension of the ladder (i.e. the two different legs) in a one-dimensional chain of resonators. This way, the diagonal links exist in the synthetic space, and their crossing poses no practical problem. One possible realization using a 1D array of parallel waveguides is to make use of two different electromagnetic modes to represent the two legs of the ladder in each waveguide, as shown in Figure 4a). This same scheme, but using atomic orbitals instead of electromagnetic modes, was recently used in an ultracold ytterbium atom system in an optical lattice to achieve one of the first experimental implementations of the Creutz ladder. [52] A system very similar to the Creutz ladder was also previously realized in cold atoms using the same approach. [53] FIG. 4. a) Creutz ladder in one spatial and one synthetic dimension. In a photonic setup two different electromagnetic modes would represent the two legs. In a cold atom system, s and p atomic orbitals would be used instead. This setup was implemented recently using ytterbium atoms in an optical lattice. [52] b) Simulation of the Creutz ladder using a rhombus chain with a large energy offset in the sites with four first neighbours (shaded in green). The particles are restricted to the top and bottom sites, but they can tunnel from one site to another via virtual processes involving the middle sites. This way, an effective lattice with the geometry of the Creutz ladder is obtained. The synthetic magnetic field is inherited from the one in the rhombus chain. [13] Another way to tackle the problem would be to simulate the Creutz ladder in a different system using the approach suggested in Ref. [13], which would make use of the topological similarities between the Creutz ladder and the rhombus chain. In that work, a photonic rhombus chain was constructed using a waveguide lattice. No synthetic dimensions are necessary to build the rhombus lattice, as no tunneling processes cross, making its physical realization considerably more straightforward. If a sufficiently high potential energy is then imposed on the spinal sites, the edge sites will then behave effectively as the two legs of a Creutz ladder (Figure 4b). A very similar photonic rhombus chain was also realized in another recent work. [54] The synthetic magnetic field can be obtained with auxiliary waveguides that retard the wavefunction appropriately, [54] or with Floquet engineering, combining a periodic modulation of the refraction index with an energy gradient. [13] In these systems, it is especially easy to observe the AB caging. In the scalar-paraxial approximation used, the direction along the waveguides behaves as the time dimension of the simulated model. That way, if light is injected into the waveguide j, α, the oscillatory movement predicted by Equation (6) will be observed as spatially periodic (breathing) movements in and out of its first neighbours. It would thus be possible to obtain the Creutz ladder analogue of these results (i.e. Figure 2d) in the existing photonic implementations of rhombus chain, using the mapping between them detailed above.
The highly localized nature of the edge states of the Creutz ladder implies that they too can be directly observed in a similar way. Distinguishing them, however, is inherently harder than simply observing their localization, because of the phase difference needed between different sites even in the simplest edge states [Equations (7,8)].
The discussion above shows that the noninteracting Creutz ladder can be implemented experimentally in photonic systems, and that interesting new phenomena are already accessible in that framework. In the following section, we will discuss the considerably more complicated question of including on-site interactions in the Creutz model, and its experimental feasibility.

III. CREUTZ-HUBBARD LADDER
As we noted earlier, interactions between particles play an interesting role in both topological and flat band systems. This motivates the inclusion of an interaction term in the Creutz Hamiltonian. In the literature, nearest neighbor [41,45,55] or on-site interactions, both attractive [56,57] and repulsive [14,15] have been considered, We will consider an on-site repulsive (U > 0) Hubbard interaction term in the Creutz Hamiltonian: Given that c ( †) jα are bosonic operators, there is no limit to the number of particles that a single site can hold. Nevertheless, for reasons that will become apparent soon, we will focus on the case where only two particles populate the ladder, the minimal case for interaction effects to occur. We are able to do this because the Creutz-Hubbard Hamiltonian preserves the number of particles. Additionally, its correlated nature prevents the use of a single-particle approach. We will use a bosonic Fock space basis {|1 i,α 1 j,β } i,α,j,β , defined for the two-particle subspace as: The nature of the photonic implementation we will consider in this work forces us to also consider a first quantization Hilbert basis of the two-particle space {|i, α; j, β } i,α,j,β . The main difference between both bases is that the states in the former are intrinsically symmetric under the exchange of both particles (and so |1 i,α 1 j,β is the same state as |1 j,β 1 i,α ), while in the latter the states |i, α; j, β and |j, β; i, α are distinct. The relationship between both bases is: Using each of these basis, any state |ψ can be expressed as: where we have defined an ordering for the sites i, α: 1A < 1B < · · · < LA < LB. The bosonic nature of the particles imposes the constraint λ i,α;j,β = λ j,β;i,α , and the components in each basis are related by η i,α;j,β = [ We will expand on the need for the two different bases in section III B.

A. Effective Hamiltonian, doublon caging and edge states
For sufficiently large values of U , the repulsive interaction can actually bind both particles together. This can be understood in energetics terms: if U J, m, the states of the system with a doubly-occupied site have a much higher energy than the others. As there is no dissipation in our model, conservation of energy means that a doubly occupied state cannot evolve into a pair of separated particles under time evolution. Both particles thus move together around the lattice as a pair, forming a doublon. Thus, a doublon is a quasiparticle composed of two particles which are bound as a result of a repulsive Hubbard interaction, in which both particles occupy the same site. A model where dissipation is taken into account would see doublons break apart after a certain lifetime, due to different mechanisms as doublon-doublon scattering or phonon scattering. [58] The fact that the doublon bands (formed by all the doubly occupied states of the Fock basis) are decoupled from the rest of the spectrum (the singly occupied states) can be exploited to execute a Schrieffer-Wolff transformation and obtain an effective Hamiltonian for the doublon subspace. [59][60][61] The resulting Hamiltonian, for a finite system, takes the form: where d jα ) 2 are the doublon creation and annihilation operators, and ∆ = [(2m 2 + 8J 2 )/U ]1 is an energy offset of the whole Hamiltonian with respect to the original scale of energy. The first sum is a renormalized Creutz model, with −J → J 2 /U, −J α → J 2 α /U, −m → m 2 /U . The second term is the Hubbard interaction term, which remains unchanged. Additionally, a chemical potential −µ = −2J 2 /U appears in the end sites, because of their different coordination number than the rest of the lattice. This same phenomenon was demonstrated in other 1D and 2D effective models for doublons. [61,62] The rungless (i.e. m = 0) π-flux Creutz-Hubbard ladder with repulsive interactions has been studied in the literature for different values of its filling factor. These studies have shown that, for a dense enough system, a Tomonaga-Luttinger liquid (TLL) of doublons is spontaneously formed, while some fraction of the particles remain unpaired and localized in the background. [14,15] The dynamics in this situation have some similarities and differences with the sparse system with one or two particles that we focus on. The most important common ground is that single particle localization for φ = π is due to AB caging for any filling factor. On the other hand, the interacting nature of the Hamiltonian makes the spectrum of the system fairly sensitive to the filling factor. Indeed, all states with double occupancy remain in the high energy region of the spectrum for sparse systems, but as the number of particles increases they eventually take a part in the ground state of the TLL. This illustrates the complexity and richness that often characterize doublon dynamics in Hubbard models.
Given that the doublon will have an electric charge 2q (with q the charge of one particle), we can expect to find analogues of the behaviour found in the single particle dynamics for half the magnetic flux than in the latter case. In particular, AB caging for doublons is expected to be found for φ = ±π/2. Because of the periodicity of the phase diagram in φ, caging should also be found for φ = ±3π/2. For these values of the flux, the independent particles are not localized, and this paints an interesting picture: an ensemble of localized doublons and free-moving particles. This is the opposite situation to the rungless π-flux regime, in which particles were caged and pairs could propagate. The behaviour for m = 0 is analogous to the single-particle case: including a vertical hopping amplitude breaks AB caging for all values of φ.
In the following, we will focus on the case with m = 0. Doublon edge states, which have been studied for different 1D and 2D systems [61][62][63][64], cannot be treated in the same way as their single particle counterparts in 1D and quasi-1D systems, as the usual bulk-boundary correspondence is not valid and the value of the Zak phase is not correlated to the existence of topological doublon edge states. [63] Consequently a different approach must be employed. In the m = 0 case, the simplicity of the ABlocalized single-particle edge states in Equations (7,8) allows us to probe for analogous doublon edge states. The relative phase between the sites must now be twice that of the single-particle states, to match the AB phase acquired in each hopping by the doublon of charge 2q: Exact numerical time evolution of these states confirm their stationary nature and identifies them as eigenstates of the Hamiltonian (9) (see Figure 5f). Photonic doublon edge states can be used to create maximally entangled multiphoton (so-called NOON ) states, which find applications in high-precision quantum metrology. [65,66] In our system, the state |2L φ + |2R φ constitutes a NOON state. A driving potential could also be used to implement a doublon transfer protocol between both ends of the ladder, which could find applications in the field of quantum information.
[61] Figure 5 shows several numerical simulations of the rungless Creutz-Hubbard ladder for different values of its parameters. To get an intuitive feeling for the dynamics, color was used to encode two important features of the wavefunctions (Figure 5a): the expectation value of the occupation number for each site, n i,α , and a measure for the doublonness of the wavefunction, D i,α . This measure is the probability that, if the site i, α is measured and found not to be empty, two particles (and not only one) are found in it. It can be calculated as a conditional probability. The numerator represents the probability of finding a doublon in site i, α if a measure is performed, and the denominator is the probability of finding that the site is not empty: The initial state for Figures 5b-d is |2 3 J. This state is stationary. g) Each particle stays localized at one end. This can also happen for any value of φ, and the state is also stationary. [see Equation (7)], for 5f it is 2L φ=π/2 [Equation (13)] and for 5g, As shown in the figure, the characteristic times for doublon dynamics are around one order of magnitude greater than the ones for single particle dynamics.

B. Bidimensional model and photonic implementation
The classical light propagating in a photonic lattice, in the paraxial approximation, follows a wave equation that can be mapped to the Schrödinger equation, and thus the amplitude of the electric field can be regarded as analogous to a wavefunction in (first quantization) quantum mechanics. This is the analogy on which quantum simulation in photonic lattices is based. To implement the first quantization Hilbert space corresponding to the two-particle 1D system, a 2D system is needed. This is a well-known method of treating the degrees of freedom in low-dimensional systems, and has been employed for both conventional [29,33,67] and topological materials. [63,68].
Using the Fock basis, the Schrödinger equation for the two-particle photonic Creutz-Hubbard system, H CH |ψ = E |ψ , turns into a system of linear equations that can be interpreted as describing a single bosonic particle hopping around a quasi-2D lattice. This system of equations is calculated in Appendix B. The resulting lattice can be obtained as the Cartesian product of two Creutz ladders. The sites in this new lattice are labeled with the four indices i, α, j, β, where i, α would be the position of the first particle in the original model, and j, β that of the second one (see Figure 6a). In this way, the diagonal sites i, α, i, α correspond to states with a doubly occupied site in the original Creutz lattice. If the particle never leaves the diagonal in the 2D lattice, that implies that the corresponding 1D doublon would not decay into two separate particles. The interaction energy U has to be reinterpreted as an on-site potential in the diagonal sites of the 2D model. For the photonic lattice to correctly emulate a bosonic system, all initial states considered must be symmetric under exchange of the two particles (this is equivalent to the aforementioned constraint λ i,α;j,β = λ j,β;i,α ). The exchange symmetry of the equations of motion, inherited from the simulated bosonic system, ensures a symmetric state will always remain symmetric.
Although the resulting quasi-2D lattice might be thought of as four-dimensional, the small size of two of its dimensions can be exploited to rearrange the lattice in a 3D geometry with some nonlocal hopping amplitudes, as shown in Figure 6b). To do this, it is useful to relabel the states with a new index ζ = 1, 2, 3, 4, corresponding to the values (α, β) = (A, A), (A, B), (B, A), (B, B). This small perspective shift gives the key to a possible scheme for an experimental implementation of the model. In order to construct the equivalent of a 3D lattice evolving in time in a photonic system, at least one dimension must FIG. 7. Quasi-2D plots for the simulations in Figure 5. a) Different time evolutions for the initial state |23,A , depending on the parameters of the system (see Figure 5b-d). b) Single particle edge state and caged particle (see Figure 5e). c) Doublon edge state (Figure 5f). d) One particle at each end (Figure 5g). Note that colour corresponds to the first quantization components (see Equation 11), and hence the apparent preference for the diagonal in the doublon collapse plot, which actually corresponds to a more even distribution in the Fock basis (Figure 5d.). be synthetic. In the Creutz-Hubbard model, the natural choice is the ζ dimension described above. This is especially convenient, because synthetic gauge fields and nonlocal hopping amplitudes are easier to implement in synthetic space than real space. [49,69] Three nonlocal hoppings would be necessary per unit cell (two for the rungless ladder), represented with dotted lines in Figure 6b. It is worth noting that, if the synthetic dimension is periodic, no nonlocal hoppings are needed.
A specific setup for the photonic Creutz-Hubbard ladder could be similar to the one proposed in a recent work. [37] In it, oscillating columns of waveguides act as sites in a 2D lattice, and the synthetic dimension is represented by the modes of oscillation of the column. Four modes in a square lattice would suffice to implement our system. In Figure 7, the quasi-2D representations of the simulations in Figure 5 are shown. This way of presenting them carries more information than its counterpart used in Figure 5, but it is also much less intuitive. The lattice sites are arranged as indicated by Figure 6a. In an experiment, the initial state |2 iα would be easy to prepare, allowing for the observation of the different phenomena in Figure 7a.

IV. CONCLUSIONS
In this work, we have studied numerically several exotic phenomena that arise in the photonic Creutz and Creutz-Hubbard ladders. The former has already been implemented experimentally, and we propose a 3D waveguide lattice setup to implement the latter, possible thanks to the new possibilities that the use of synthetic dimensions can offer. It is our understanding that the current state of the art in the rapidly evolving field of photonics allows this proposal to be implemented in a real system. As far as we are aware, this would mark the first instance of a quasi-1D topological insulator with interactions in a waveguide-lattice-type quantum simulator, representing an increase in complexity compared to the strongly correlated systems that have been implemented so far in waveguide lattices (which were both 1D and topologically trivial). [29][30][31][32][33] In the single-particle model, we have focused on the caging of particles and the topologically protected edge states, which are AB-localized in the absence of rung-like hoppings, even when no localized bulk energy eigenstates exist. We then considered the two-particle model, the minimal model with interactions, and studied the doublon collapse, the different caging regimes (in particular, the doublon caging for m = 0, φ = π/2) and the doublon edge states. These results show that the phase diagram of these systems have many interesting phenomena to observe away from the usually-studied φ = π limit. In particular, analogous caging regimes will be present for multiplons with different numbers of particles at different values of the magnetic flux, giving rise to complex physical situations that require further study.
Once realized, these photonic systems could be the first to observe many of the interesting features of the Creutz and Creutz-Hubbard models, marking a milestone in the research of low-dimensional topological systems and flat band models. Moreover, in the interacting case, it would mark a step forward in the understanding of the complicated interplay that occurs between interactions and topology, the details of which are still not completely understood. The Creutz and Creutz-Hubbard models constitute an ideal playground in which flat band dynamics, particle interactions and nontrivial topology can be studied in depth, separately or at the same time, allowing for a complete study of the rich interplay between them.