Topological edge states of interacting photon pairs emulated in a topolectrical circuit

TL;DR

This paper demonstrates a classical topolectrical circuit that simulates topological states of interacting photon pairs, enabling the study of quantum topological phenomena through accessible measurements.

Contribution

The authors design a topolectrical circuit that emulates a quantum system with topological two-photon states, providing a new platform for exploring quantum topological physics.

Findings

Successfully simulated two-photon topological states in a circuit

Measured bulk and edge state frequencies and topological invariants

Reconstructed two-photon probability distribution from circuit data

Abstract

Topological physics opens up a plethora of exciting phenomena allowing to engineer disorder-robust unidirectional flows of light. Recent advances in topological protection of electromagnetic waves suggest that even richer functionalities can be achieved by realizing topological states of quantum light. This area, however, remains largely uncharted due to the number of experimental challenges. Here, we take an alternative route and design a classical structure based on topolectrical circuits which serves as a simulator of a quantum-optical one-dimensional system featuring the topological state of two photons induced by the effective photon-photon interaction. Employing the correspondence between the eigenstates of the original problem and circuit modes, we use the designed simulator to extract the frequencies of bulk and edge two-photon bound states and evaluate the topological invariant directly from the measurements. Furthermore, we perform a reconstruction of the two-photon probability distribution for the topological state associated with one of the circuit eigenmodes.

Equations

$$(\epsilon -2U){\beta }_{2m,2m}=-2J\left[{\beta }_{2m+1,2m}+{\beta }_{2m,2m-1}\right]+P{\beta }_{2m-1,2m-1},$$

$$(\epsilon -2U){\beta }_{2m+1,2m+1}=-2J\left[{\beta }_{2m+2,2m+1}+{\beta }_{2m+1,2m}\right]+P{\beta }_{2m+2,2m+2},$$

$$\epsilon {\beta }_{m,n}=-J\left[{\beta }_{m+1,n}+{\beta }_{m-1,n}+{\beta }_{m,n+1}+{\beta }_{m,n-1}\right],\text{mspace}8.0mu(m\ne n)$$

$$(\epsilon -2U-2j){\beta }_{2m,2m}=(j+P){\beta }_{2m-1,2m-1}+j{\beta }_{2m+1,2m+1},$$

$$(\epsilon -2U-2j){\beta }_{2m+1,2m+1}=j{\beta }_{2m,2m}+(j+P){\beta }_{2m+2,2m+2}.$$

$$(\epsilon -2U-j){\beta }_{11}=(j+P){\beta }_{22},$$

$$(\epsilon -2U-j){\beta }_{NN}=j{\beta }_{N-1,N-1},$$

$$z_{1,2}=\frac{j+P}{2j^3}\left[2jP+P^2\pm \sqrt{(2jP+P^2{)}^2+4j^4}\right].$$

$${\epsilon }_{1,2}=2U+j-\frac{1}{2j}\left[2jP+P^2\pm \sqrt{(2jP+P^2{)}^2+4j^4}\right].$$

$$P>0\text{mspace}6.0mu\text{or}\text{mspace}6.0muP<-2J^2/U.$$

$$\omega c_{mn}=\sum _{m^{\prime },n^{\prime }}{H}_{mn,m^{\prime }{n}^{\prime }}{c}_{m^{\prime }{n}^{\prime }}-i\gamma c_{mn}-i\kappa E_{mn}.$$

$${\omega }_k{\beta }_{mn}^{(k)}=\sum _{m^{\prime },n^{\prime }}{H}_{mn,m^{\prime }{n}^{\prime }}{\beta }_{m^{\prime }{n}^{\prime }}^{(k)}.$$

$$\sum _{m,n}{\beta }_{mn}^{(k)}{\left[{\beta }_{mn}^{(k^{\prime })}\right]}^{\ast }={\delta }_{k{k}^{\prime }}.$$

$$c_{mn}=\sum _k{\alpha }_k{\beta }_{mn}^{(k)}.$$

$$\omega \sum _k{\alpha }_k{\beta }_{mn}^{(k)}=\sum _{m^{\prime },n^{\prime },k}{H}_{mn,m^{\prime }{n}^{\prime }}{\alpha }_k{\beta }_{m^{\prime }{n}^{\prime }}^{(k)}-i\gamma \sum _k{\alpha }_k{\beta }_{mn}^{(k)}-i\kappa E_{mn}=$$

$$\stackrel{\text{eqref}Eigenmode}{=}\sum _k{\alpha }_k({\omega }_k-i\gamma ){\beta }_{mn}^{(k)}-i\kappa E_{mn}.$$

$$\left(\omega -{\omega }_k+i\gamma \right){\alpha }_k=-i\kappa \sum _{m,n}{E}_{mn}{\left[{\beta }_{mn}^{(k)}\right]}^{\ast }.$$

$$c_{mn}\stackrel{\text{eqref}Expansion}{=}\sum _k{\alpha }_k{\beta }_{mn}^{(k)}=-i\kappa \sum _{m^{\prime },n^{\prime },k}{E}_{m^{\prime }{n}^{\prime }}\frac{{\left[{\beta }_{m^{\prime }{n}^{\prime }}^{(k)}\right]}^{\ast }{\beta }_{mn}^{(k)}}{\omega -{\omega }_k+i\gamma }.$$

$$E_{mn}=\frac{E_0}{2}\left[{\delta }_{mp}{\delta }_{nq}+{\delta }_{mq}{\delta }_{np}\right],$$

$$c_{mn}=-i\kappa E_0\sum _k\frac{{\beta }_{mn}^{(k)}{\left[{\beta }_{pq}^{(k)}\right]}^{\ast }}{\omega -{\omega }_k+i\gamma }.$$

$$\mathfrak{I}(p,q)=\sum _{m,n}|c_{mn}{|}^2={\kappa }^2|E_0{|}^2\sum _{m,n,k,k^{\prime }}\frac{{\beta }_{mn}^{(k)}{\left[{\beta }_{mn}^{(k^{\prime })}\right]}^{\ast }{\left[{\beta }_{pq}^{(k)}\right]}^{\ast }{\beta }_{pq}^{(k^{\prime })}}{(\omega -{\omega }_k+i\gamma )(\omega -{\omega }_{k^{\prime }}-i\gamma )}$$

$$\stackrel{\text{eqref}Orthogonality}{=}{\kappa }^2|E_0{|}^2\sum _{k,k^{\prime }}\frac{{\delta }_{k{k}^{\prime }}{\left[{\beta }_{pq}^{(k)}\right]}^{\ast }{\beta }_{pq}^{(k^{\prime })}}{(\omega -{\omega }_k+i\gamma )(\omega -{\omega }_k^{\prime }-i\gamma )}={\kappa }^2|E_0{|}^2\sum _k\frac{|{\beta }_{pq}^{(k)}{|}^2}{(\omega -{\omega }_k{)}^2+{\gamma }^2}.$$

$$\mathfrak{I}(p,q)\equiv \sum _{m,n}|c_{mn}{|}^2={\kappa }^2|E_0{|}^2\sum _k\frac{|{\beta }_{pq}^{(k)}{|}^2}{(\omega -{\omega }_k{)}^2+{\gamma }^2},$$

$$\left(\epsilon -2U\right){\beta }_{2m+1,2m+1}=-2J{\beta }_{2m,2m+1}-2J{\beta }_{2m+1,2m+2}+P{\beta }_{2m+2,2m+2},$$

$$\left(\epsilon -2U\right){\beta }_{2m,2m}=-2J{\beta }_{2m-1,2m}-2J{\beta }_{2m,2m+1}+P{\beta }_{2m-1,2m-1},$$

$$\epsilon {\beta }_{mn}=-J{\beta }_{m-1,n}-J{\beta }_{m+1,n}-J{\beta }_{m,n-1}-J{\beta }_{m,n+1}\text{mspace}10.0mu(m\ne n).$$

$$(\epsilon -2U){\beta }_{11}=-2J{\beta }_{12}+P{\beta }_{22},$$

$$\epsilon {\beta }_{1m}=-J\left[{\beta }_{1,m-1}+{\beta }_{1,m+1}+{\beta }_{2,m}\right],$$

$$(\epsilon -2U){\beta }_{NN}=-2J{\beta }_{N,N-1},$$

$${\omega }_0^2=\frac{1}{L{C}_j},\text{mspace}10.0mu\epsilon =\frac{{\omega }_0^2}{{\omega }^2}-4,$$

$$U=\frac{C_P+C_U}{2C_j},\text{mspace}10.0muP=-\frac{C_P}{C_j}.$$

$$\frac{U_{mn}}{L}-{\omega }^2{C}_j\left[\left(U_{mn}-U_{m,n+1}\right)+\left(U_{mn}-U_{m,n-1}\right)+\left(U_{mn}-U_{m+1,n}\right)+\left(U_{mn}-U_{m-1,n}\right)\right]=0$$

$$\epsilon U_{mn}=-U_{m,n+1}-U_{m,n-1}-U_{m+1,n}-U_{m-1,n},$$

$$\frac{U_{1n}}{L}-{\omega }^2{C}_j\left[\left(U_{1n}-U_{1,n-1}\right)+\left(U_{1n}-U_{1,n+1}\right)+\left(U_{1n}-U_{2n}\right)+U_{1n}\right]=0,$$

$$\epsilon U_{1n}=-U_{1,n-1}-U_{1,n+1}-U_{2,n},$$

$$\begin{array}{rl}& \frac{U_{2m+1,2m+1}}{L}-{\omega }^2{C}_j[\left(U_{2m+1,2m+1}-U_{2m,2m+1}\right)+\left(U_{2m+1,2m+1}-U_{2m+2,2m+1}\right)\\ & +\left(U_{2m+1,2m+1}-U_{2m+1,2m}\right)+\left(U_{2m+1,2m+1}-U_{2m+1,2m+2}\right)]-\\ & -{\omega }^2{C}_p\left(U_{2m+1,2m+1}-U_{2m+2,2m+2}\right)-{\omega }^2{C}_U{U}_{2m+1,2m+1}=0,\end{array}$$

Full text

1. Supplementary Note 1 – Analysis of the strong interaction limit $U\gg J$: Su-Schrieffer-Heeger model for doublons
2. Supplementary Note 2 – Eigenmode tomography
3. Supplementary Note 3 – Mapping of tight-binding model onto the topolectrical circuit: ideal lossless case
4. Supplementary Note 4 – Analysis of topolectrical circuit with realistic losses
5. Supplementary Note 5 – Distribution of elements in the experimental setup
6. Supplementary Note 6 – Mutual inductance
7. Supplementary Note 7 – Circuit impedance and doublon spectroscopy
8. Supplementary Note 8 – Scattering states
9. Supplementary Note 9 – Effects of disorder
10. Supplementary Note 10 – Evaluation of topological invariant from experimental data

Supplementary Note 1 – Analysis of the strong interaction limit $U\gg J$: Su-Schrieffer-Heeger model for doublons

To provide a simple proof of the topological origin of our model, we analyze the strong interaction limit, when $U\gg J$. In such a situation, doublon bands are well-separated from the continuum of two-photon scattering states so that the mixing between doublons and scattering states is negligible. As such, we introduce an effective doublon Hamiltonian which captures the dynamics of doublons excluding other redundant degrees of freedom.

To derive the effective doublon Hamiltonian, we start from the eigenvalue equations provided in the article main text (Methods section):

[TABLE]

In the limiting case $U\gg J$, $\beta_{mm}$ are the dominant coefficients of the doublon wave function with the rest of coefficients $\beta_{m+s,m}$ decaying with the index $s$. Therefore, we neglect all terms proportional to $\beta_{m+2,m}$, $\beta_{m+3,m}$, etc. in Supplementary Equations (1)-(3), treating terms proportional to $\beta_{m+1,m}$ as a perturbation.

Using the approximation $\beta_{m+1,m}\approx-\left(\beta_{mm}+\beta_{m+1,m+1}\right)/(2\,U)$, we derive the following approximate eigenvalue equation for doublon bands:

[TABLE]

The equations for $\beta_{11}$ and $\beta_{NN}$ corresponding to the physical edges of the array are modified if compared to the bulk sites:

[TABLE]

where $j=J^{2}/U$ is the effective doublon hopping rate associated with two consecutive single-particle tunnelings to the neighboring cavity. Here, the array length $N$ is assumed to be odd.

The set of Supplementary Equations (4)-(5) corresponds to the Su-Schrieffer-Heeger (SSH) model [1] with two alternating tunneling amplitudes $j$ and $j+P$, which is known to be the simplest one-dimensional topological model. The only difference from the canonical SSH model is the interaction-induced detuning of the edge sites by $j$ captured by Supplementary Equations (6), (7).

Solving the system of Supplementary Equations (4)-(5) together with boundary conditions Supplementary Equations (6)-(7), we find two states with the localization ratio $z=\beta_{mm}/\beta_{m-2,m-2}$ given by

[TABLE]

Edge-localized states correspond to $|z|<1$. The energies of these states read:

[TABLE]

Supplementary Equation (8) shows that the higher-energy edge state $\varepsilon_{2}$ is localized for any $P\not=0$, while the lower-energy state $\varepsilon_{1}$ disappears when

[TABLE]

In our case with $U/J=7.09$ and $P/J=-4.18$ the latter condition is fulfilled, which means that only the higher-energy state $\varepsilon_{2}$ persists. Note also that the condition Supplementary Equation (10) is equivalent to $|j+P|>j$ which ensures that the site $(N,N)$ is the weak link edge and the respective edge state can be interpreted as the topological state inherent to SSH model.

Supplementary Note 2 – Eigenmode tomography

In an experimental situation, we can examine the excitation of the system applying voltage to various sites. However, we have no direct access to the system eigenmodes and associated probability distributions, which eventually limits our possibilities to observe topological physics in the designed setup.

To overcome this limitation, we elaborate on the procedure of eigenmode reconstruction (tomography) following the proposals of our recent work [2]. The key idea of this method is to collect the information on system excitation when the voltage is applied to various pairs of sites $(m,n)$ symmetrically with respect to the diagonal. The eigenmode profile is then recovered as a sum of squares of voltages in all nodes of the circuit plotted as a function of the feeding point $(m,n)$. In this section, we provide a summary of the proposed technique highlighting further applications of the developed method and its applicability to a wide class of systems described by tight-binding equations.

We consider a two-dimensional system depicted in Supplementary Figure 1. We describe the behavior of the system under external coherent driving with the following semi-empirical coupled mode equations:

[TABLE]

Here $H_{mn,m^{\prime}n^{\prime}}$ is the Hamiltonian of a closed system which is assumed Hermitian, $\gamma$ quantifies the dissipation rate, and $E_{mn}$ stands for the applied electromotive force (or injected current, depending on the interpretation of $c_{mn}$ coefficients) at frequency $\omega$. In an actual experimental situation, we can apply different drivings $E_{mn}$ at different frequencies and measure the resultant intensity distribution $|c_{mn}|^{2}$. Our goal is to elaborate a protocol to reconstruct the profile of the system eigenmode $\beta_{mn}^{(k)}$ corresponding to the eigenfrequency $\omega_{k}$ which is the solution of the eigenvalue problem

[TABLE]

Since the Hamiltonian of the closed system is Hermitian, the eigenfrequencies $\omega_{k}$ are purely real and the eigenmodes $\beta_{mn}^{(k)}$ are mutually orthogonal, i.e.

[TABLE]

Furthermore, the modes $\beta_{mn}^{(k)}$ form a complete basis and the coefficients $c_{mn}$ can be expanded as

[TABLE]

Putting the expansion Supplementary Equation (14) into Supplementary Equation (11), we obtain:

[TABLE]

Using the orthogonality property Supplementary Equation (13), we recover that

[TABLE]

Physically, Supplementary Equation (17) expresses an intuitive fact that the larger the overlap of the driving profile with the eigenmode, the larger is the contribution of the eigenmode into the driven-dissipative system stationary state. On the other hand, the closer is the driving frequency to the resonance, the larger is the contribution of a given eigenmode to the driven-dissipative system stationary state. The coefficients $c_{mn}$ are expressed as follows:

[TABLE]

At this point we assume that the driving profile (a) includes only one or two points; (b) is symmetric with respect to the diagonal $m=n$. In the other words,

[TABLE]

where $p$ and $q$ provide the coordinates of the excitation point(s). This assumption immediately simplifies the sum. Since the system is symmetric with respect to the diagonal, all non-degenerate eigenmodes are either even or odd; degenerate modes can also be enforced to satisfy this condition by choosing suitable linear combinations. Next, all odd modes have zero overlaps with symmetric pumping profile and they drop out of the sum. The remaining even modes have $\beta_{pq}^{(k)}=\beta_{qp}^{(k)}$. Hence,

[TABLE]

Next, we construct the following quantity:

[TABLE]

Thus, for symmetric pumping into $(p,q)$ and $(q,p)$ with the same amplitudes $E_{0}/2$ and $E_{0}/2$ and equal phases the sum of squares of field amplitudes in all sites of the system reads:

[TABLE]

where the sum extends only over the symmetric eigenmodes. If the driving frequency is close enough to the eigenfrequency $\omega_{k}$, the obtained distribution $\mathfrak{I}(p,q)$ will closely resemble the distribution $|\beta_{pq}^{(k)}|^{2}$, which is actually the eigenmode intensity (or two-photon probability distribution in the original 1D two-particle problem).

The outlined eigenmode reconstruction protocol works especially well once the spectral distance from the given mode to the rest of the modes exceeds the dissipation rate $\gamma$. In such a case, reconstruction of the eigenmode will be very precise as can be seen from the comparison made in the main text. This technique is also useful for the group of modes, e.g. doublon band, well-separated from the rest of the modes.

Note also that the developed eigenmode reconstruction technique is very general and applicable to a wide range of physical systems. In this work, we apply this protocol to LC circuits. Further identification of driven-dissipative model Supplementary Equation (11) with Kirchhoff’s equations for the electric circuit is provided in Supplementary Notes 3-4.

Supplementary Note 3 – Mapping of tight-binding model onto the topolectrical circuit: ideal lossless case

As indicated in the article main text, the two-particle one-dimensional quantum problem with extended Bose-Hubbard Hamiltonian is described by the following system of tight-binding equations:

[TABLE]

with $\beta_{mn}=\beta_{nm}$ and open boundary conditions at the edges and corners:

[TABLE]

where $N$ is the size of the system.

In this Supplementary Note, we analyze the 2D system in Fig. 2a with Kirchhoff’s circuit laws and express the parameters of the Hamiltonian $U$ and $P$ in terms of the circuit parameters. In our analysis we assume that the voltages in the circuit nodes $U_{mn}\propto e^{-i\omega t}$, and hence $\ddot{U}_{mn}=-\omega^{2}\,U_{mn}$. We also assume $J=1$ and denote

[TABLE]

To establish one-to-one correspondence of the initial tight-binding problem with LC circuit, we analyze several representative situations.

1. Site $(m,n)$ with $m\not=n$, not at the edge of the system. From the first Kirchhoff’s law we get

[TABLE]

or with designations Supplementary Equations (29)-(30)

[TABLE]

which is consistent with Supplementary Equation (25).

2. Site $(1,n)$ with $n\not=1$, at the edge of the system.

[TABLE]

where the last term in the square bracket is associated with extra capacitance $C_{j}$ connected to the ground. Supplementary Equation (33) can be rearranged to yield

[TABLE]

which reproduces open boundary condition for the tight-binding model Supplementary Equation (27). The equations for the sites $(m,1)$, $(m,N)$ and $(N,n)$ (not at the corner) where $N$ is the system size are completely analogous.

3. Site $(2m+1,2m+1)$ at the diagonal of the system, not the corner one.

[TABLE]

which can be rewritten as

[TABLE]

and coincides with Supplementary Equation (23) provided the pattern of voltages is symmetric, i.e. $U_{2m,2m+1}=U_{2m+1,2m}$ and $U_{2m+2,2m+1}=U_{2m+1,2m+2}$.

4. Site $(2m,2m)$ at the diagonal of the system

[TABLE]

which can be recast in the form

[TABLE]

in agreement with Supplementary Equation (24).

5. Site $(1,1)$ at the corner of the system. For this corner site, we include two capacitances $C_{j}$ to the ground to compensate the absence of the two neighbors.

[TABLE]

which is equivalent to

[TABLE]

exactly reproducing an open boundary condition Supplementary Equation (26) for the corner of the system.

6. Site $(N,N)$ at the corner of the system ($N$ is assumed to be odd). In contrast to $(1,1)$ site, besides two capacitances $C_{j}$ connected to the ground, we also include an extra capacitance $C_{P}$ connected to the ground:

[TABLE]

which can be rearranged to yield

[TABLE]

in agreement with open boundary condition for the corner of the system, Supplementary Equation (28).

Hence, the proposed experimental setup allows us to emulate two-body physics in Bose-Hubbard model, but with several constraints on parameters that enter the Bose-Hubbard Hamiltonian: (i) $P$ is always negative, cf. Supplementary Equation (30); (ii) $U>|P|/2$, i.e. we cannot realize the regime of too weak on-site interactions; (iii) to emulate boson states, the pattern of voltages in a circuit should be symmetric with reasonable accuracy.

Supplementary Note 4 – Analysis of topolectrical circuit with realistic losses

As we have demonstrated, the lossless 2D topolectrical circuit corresponds precisely to the Bose-Hubbard model. However, realistic circuits necessarily have losses, and in this section, we examine excitation of the system taking the effect of loss into account. We assume that the system is excited with a current source and current is pumped into one or several lattice sites. We aim to compare the governing equations with the simple driven-dissipative model outlined in Supplementary Note 2, which is described by the equation

[TABLE]

where $\omega$ is a driving frequency, $c_{mn}$ describe the stationary state of the system, $\gamma$ is the dissipation rate and $\varkappa$ is the coupling coefficient.

In circuit analysis, we assume $e^{-i\omega t}$ time dependence of voltages, which yields the following impedances:

[TABLE]

For sufficiently small losses, the ratios of impedances read

[TABLE]

[TABLE]

Similarly to Supplementary Note 3, we analyze several characteristic situations:

1. Site $(m,n)$ with $m\not=n$, not at the edge of the system. First Kirchhoff’s law now yields:

[TABLE]

which can be rearranged as

[TABLE]

where $\gamma$ is a frequency-dependent damping

[TABLE]

Note also the sign in front of $\gamma$ in Supplementary Equation (51): it is different from the sign in Supplementary Equation (43). This happens due to the definition of “energy” variable Supplementary Equation (29) which ensures that imaginary parts of $\omega$ and $\varepsilon$ have opposite signs:

[TABLE]

Since in the dissipative case eigenmode frequency has negative imaginary part, “energy variable” $\varepsilon$ should have positive imaginary part.

2. Site $(1,n)$ with $n\not=1$, at the edge of the system. In a similar way we derive an equation

[TABLE]

which yields

[TABLE]

3. Site $(2m+1,2m+1)$ at the diagonal of the system, not the corner one.

[TABLE]

This yields an equation

[TABLE]

where for the diagonal sites

[TABLE]

Thus, the diagonal elements have extra loss in tunneling $P^{\prime}$ and on-site losses $\gamma^{\prime}$ stemming from the insertion of additional elements. This problem, however, can be circumvented by requiring that

[TABLE]

in which case all elements (in the bulk and at the diagonal) are characterized by the same magnitude of loss, though depending on frequency. However, even if Supplementary Equation (59) is not fulfilled, doublon bands will remain almost unaffected, since the maximal voltages are expected at the diagonal sites only.

4. Site $(2m,2m)$ at the diagonal of the system.

[TABLE]

which yields

[TABLE]

5. Site $(1,1)$ at the corner of the system.

[TABLE]

which yields

[TABLE]

6. Site $(N,N)$ at the corner of the system [$N$ is assumed to be odd].

[TABLE]

which can be transformed as

[TABLE]

To summarize, the designed LC circuit is described by the same driven-dissipative system Supplementary Equations (43) apart from the following differences:

• •

dissipation coefficient $\gamma$ depends on the frequency of driving $\omega$;

• •

the magnitude of dissipation for diagonal and off-diagonal sites is different unless an additional constraint Supplementary Equation (59) is fulfilled;

• •

instead of frequency in Supplementary Equation (43), we deal with the auxiliary “energy variable” $\varepsilon$ which is inversely proportional to frequency and therefore has a positive imaginary part in the dissipative case.

• •

instead of external field $E_{mn}$ we have a combination $Z_{C}I_{mn}$ which is frequency-dependent.

In all other aspects, the standard driven-dissipative model captures all essential features of the proposed topolectrical circuit and hence the tomography technique developed in Supplementary Note 2  ‣ Supplementary Note 2 – Eigenmode tomography can be applied to the designed LC circuit.

Supplementary Note 5 – Distribution of elements in the experimental setup

In contrast to the idealized model of electric circuit discussed above, the actual values of inductances and capacitances of setup elements slightly fluctuate from one element to another. Hence, the experimental circuit possesses an inherent disorder in the tunneling constants $J$ and $P$ (ascribed to capacitors $C_{J}$ and $C_{P}$) as well as in the interaction strength $U$ and on-site resonant frequencies $\omega_{0}$, depending on capacitors $C_{U}$ and inductors $L$, respectively. The elements used in our experimental sample are Murata GRM32RR71H105KA01 for $C_{J}$, Murata GRM31CR61A476ME15L for $C_{P}$, Murata GRM31CR71C106KAC7 for $C_{U}$, and Bourns RLB1314-220KL for $L$. In the process of fabrication, we created a detailed map of elements that allows us to determine the precise value of the given circuit bond. The distributions of the element values are shown in Supplementary Figure 3. As seen from the histograms, actual mean values in the fabricated circuit are $L=22.77\mspace{2.0mu}{\rm\mu H}$, $C_{J}=1.0024\mspace{2.0mu}{\rm\mu F}$, $C_{U}=10.031\mspace{2.0mu}{\rm\mu F}$, and $C_{P}=4.1863\mspace{2.0mu}{\rm\mu F}$ which are slightly different from the ones provided in specifications. Typical fluctuations of parameters are of the order of 1-2%. We use these measured mean values for numerical simulations described in the Article main text.

Supplementary Note 6 – Mutual inductance

Some elements in the experimental setup can have extra parasitic couplings besides the couplings introduced intentionally. This effect is especially pronounced for inductive coils, which have relatively large size compared to the inter-coil distance, as can be seen from Fig. 1 of the article main text. We estimate now the effects of such inductive coupling on the performance of the studied setup.

Incorporating mutual inductance between the coils (grounding elements) into the Kirchhoff’s rule for the site $(2m+1,2m+1)$, which we consider as an example, we get:

[TABLE]

where $U_{\rm ind}$ denotes the voltage induced in the coil connecting the site $(2m+1,2m+1)$ to the ground by all surrounding coils. A diagonal capacitor $C_{P}$ is present between the sites $(2m+1,2m+1)$ and $(2m+2,2m+2)$. Taking into account the interaction of the given coil with its nearest neighbors and also with the coils having both coordinates $m$, $n$ shifted by $1$ (diagonal neighbors), we can express this induced voltage as

[TABLE]

where $I_{m,n}^{(g)}$ is a current through the inductance which connects the site $(m,n)$ to the ground, $Z_{M}$ and $Z_{M}^{\prime}$ are the impedances corresponding to the mutual inductance between the nearest neighbors and diagonal neighbors, respectively.

In the limit of strong interaction $U\gg 1$, which is the case in the considered model, all voltages are mostly concentrated at the diagonal sites of the circuit in the frequency range of interest. Then, the above expression takes the form

[TABLE]

where $I_{2m,2m}$ and $I_{2m+2,2m+2}$ denote external currents applied to the corresponding sites. Then, Supplementary Equation (66) reads

[TABLE]

If only one site of the circuit is driven, which is the case in our measurements, then the second term at the right-hand side vanishes, and we finally obtain the equation

[TABLE]

which has the same form as the equation describing the circuit without inductive couplings up to renormalization of the effective model parameters caused by the extra factor $(1-Z_{M}^{\prime}/Z_{P})$.

Supplementary Note 7 – Circuit impedance and doublon spectroscopy

To demonstrate doublon states experimentally, we need a technique to distinguish doublon modes from the rest of the modes supported by the two-dimensional sample. In this regard, it is useful to consider the quantity [2]

[TABLE]

where $\varphi_{ij}^{mn}$ is the potential at the site $(i,j)$ of the circuit when the external driving voltage at frequency $f$ is applied to the site $(m,m)$. Since doublon modes are characterized by the voltage maxima at the diagonal sites, this quantity exhibits resonant peaks at frequencies of doublon modes. If the external voltage is applied to the site $(N,N)$, where $N=15$ in our case, a single peak corresponding to the doublon edge state is observed. On the other hand, if any other diagonal site is driven, then two peaks centered at frequencies of bulk doublon bands emerge, Supplementary Figure 4a.

We now demonstrate that the frequency dependence of the quantity $S_{m}$ is consistent with the results of the circuit impedance spectroscopy. Indeed, one can introduce Green’s matrix of the circuit $\hat{G}$ [3], which is by definition related to the admittance matrix $\hat{Y}$ introduced in the main text as

[TABLE]

Then, potentials at the nodes of the externally driven circuit are related to the driving current $I_{mn}$ as

[TABLE]

and therefore $\varphi_{nn}^{mm}=G_{nn,mm}\,I_{\rm ext}$ where $I_{\rm ext}$ is a value of the driving current. Hence, Supplementary Equation (71) takes the form

[TABLE]

At the same time, the characteristic impedance $R_{m}$ between the given node $(m,m)$ of the circuit and the ground is simply given by the diagonal element of the Green’s matrix:

[TABLE]

As seen from Supplementary Figure 4b, the characteristic impedance of the circuit $R_{m}$ demonstrates a very similar structure of resonant peaks compared to the quantity $S_{m}$ for various positions of the driven site $(m,m)$, which highlights the dominant role of diagonal entries of the Green’s matrix in our system.

The equivalent scheme of the experimental setup with an external source applied to the diagonal site has the form shown in Supplementary Figure 4c. It represents a series connection of the total circuit impedance with the real part $R_{c}$ and the voltage source equivalent impedance $R_{s}=50\mspace{3.0mu}{\rm\Omega}$. In the experiment, the voltage drop between the given site of the circuit and the ground $U_{c}$ is measured as a function of the driving frequency $f$. The external voltage $U_{\rm ext}$ is fixed and set to the value 1 V. Then, the circuit impedance $R_{c}$ can be obtained as

[TABLE]

while the external current flowing into the circuit is given by the relation $I_{\rm ext}(f)=U_{\rm ext}/(R_{c}(f)+R_{s})$. All calculations in the article main text are carried on for the constant external current $I_{\rm ext}$ flowing into the system, whereas the experimental data is obtained for the fixed external voltage $U_{\rm ext}$. To make a fair comparison of these results, we multiply the experimental values of $S_{m}(f)$ by the factor $I^{2}(f_{0})/I^{2}(f)$ with $f_{0}=6$ kHz defining a calibrating value of the external current. Experimental results on doublon spectroscopy recalculated in this way are presented in Fig. 2c of the article main text.

Supplementary Note 8 – Scattering states

Besides two bulk doublon bands and doublon edge state, there exists a vast set of the system eigenmodes termed as scattering states. These modes shown by the shaded region in the dispersion (Fig. 2a of the article main text) correspond to such a state of two photons when they are typically located at distinct resonators and possess energy equal to the sum of single-photon energies. Supplementary Figure 5 shows the probability distributions $|\beta_{mn}|^{2}$ for such states obtained by the diagonalization of the tight-binding Hamiltonian in the absence of dissipation. Here, panels (a-d) correspond to the eigenmodes of an ideal system without any disorder. Introducing the disorder in coupling constants $J$ and $P$ with the uniform distribution within the range $\pm 10$%, we observe only slight changes in the eigenmode intensity patterns accompanied by weak energy shifts, Supplementary Figure 5(e-h). However, when the strength of disorder is increased further up to 30%, quite strong distortions of eigenmode profiles are observed and the energy shifts become comparable with the spectral distance between the spectrally close modes, Supplementary Figure 5(i-l).

Supplementary Note 9 – Effects of disorder

To examine the impact of disorder and losses on the results of tomography procedure (Supplementary Note 2), we simulate eigenmode tomography for the designed circuit taking into account Ohmic losses as well as fluctuations in the values of all lumped elements $L$, $C_{J}$, $C_{U}$, and $C_{P}$. For the sake of simplicity, we consider uniformly distributed fluctuations, and Ohmic losses are assumed to be constant and equal to their maximal possible values within the considered spectral range according to elements specifications.

The evolution of low-energy doublon mode, doublon edge state, and high-energy doublon mode with the increase of disorder is shown in panels (a-c), (d-f), and (g-i) of Supplementary Figure 6, respectively. It is seen that enhanced localization of doublon modes can be observed even at 10% disorder, which finally results in the formation of localized states in the bulk of the circuit at 30% disorder. At the same time, the edge state easily survives 10% disorder, since it is spectrally well-separated from the bulk doublon bands having a relatively small spatial overlap with them. However, strong 30% disorder can mix it with bulk doublon states, as seen in Supplementary Figure 6f. It should be stressed that the considered levels of disorder in the values of circuit elements exceed those expected for the experimental circuit and the corresponding results are calculated for the illustrative purpose only.

Supplementary Note 10 – Evaluation of topological invariant from experimental data

The definition of topological invariant for interacting many-body systems is currently an open problem which is being actively investigated. In our specific case, however, the clue to topological characterization is provided by the fact that the effective doublon Hamiltonian corresponds to that of the Su-Schrieffer-Heeger model once strong interaction regime $U\gg J$ is realized. This is the case for our experimental sample with $U/J=7.09$ and $P/J=-4.18$.

The effective Hamiltonian for doublons written in the basis of isolated resonators’ eigenstates takes the form (see Supplementary Note 1 for details):

[TABLE]

where $j_{1}=J^{2}/U$ and $j_{2}=j_{1}+P$ are effective tunneling amplitudes for doublons. Hence, for the given Bloch eigenmode the ratio of two components of the doublon wave function $\left|\psi\right>=(\psi_{\rm A},\psi_{\rm B})^{T}$ is given by:

[TABLE]

We notice that the effective doublon Hamiltonian Supplementary Equation (77) is represented in chiral basis and therefore the winding number is determined by plotting its off-diagonal block $q(k)=j_{1}+j_{2}\,e^{-ik}$ on the complex plane [4]. Furthermore, according to Supplementary Equation (78) $q(k)=\varepsilon(k)\,\psi_{\rm A}/\psi_{\rm B}$, where doublon energy $\varepsilon(k)$ is purely real. Hence, the winding number can be found by plotting the ratio $\psi_{\rm A}/\psi_{\rm B}$ for a particular Bloch eigenmode.

In an experimental situation, we measure the distribution of voltages at the diagonal of the sample: $U_{\rm A}(n)$ and $U_{\rm B}(n)$. To get a result, corresponding to the given value of wave number $k$, we perform a discrete Fourier transform of those voltages extracting

[TABLE]

and then plotting the ratio $U_{\rm A}(k)/U_{\rm B}(k)$ on the complex plane.

Calculating the winding number from experimental data, we assume that the pattern of voltages excited in a topolectrical circuit by the source inserted in the middle of the diagonal resembles the pattern of voltages expected for the eigenmode, which is justified in the case of doublon bands well-separated from the scattering continuum.

Supplementary References

• [1]

Su, W. P., Schrieffer, J. R. & Heeger, A. J.

Solitons in Polyacetylene.

Phys. Rev. Lett. 42, 1698–1701

(1979).

• [2]

Gorlach, M. A. et al.

Simulation of two-boson bound states using arrays of driven-dissipative coupled linear optical resonators.

Phys. Rev. A 98, 063625

(2018).

• [3]

Chung, F. & Yau, S.-T.

Discrete Green’s Functions.

Journal of Combinatorial Theory, Series A 91, 191–214

(2000).

• [4]

Ryu, S., Schnyder, A. P., Furusaki, A. & Ludwig, A. W. W.

Topological insulators and superconductors: tenfold way and dimensional hierarchy.

New J. Phys. 12, 065010

(2010).