Asymmetric topological pumping in nonparaxial photonics

Topological photonics was initially inspired by the quantum-optical analogy between the Schrödinger equation for an electron wavefunction and the paraxial equation for a light beam. Here, we reveal an unexpected phenomenon in topological pumping observed in arrays of nonparaxial optical waveguides where the quantum-optical analogy becomes invalid. We predict theoretically and demonstrate experimentally an asymmetric topological pumping when the injected field transfers from one side of the waveguide array to the other side whereas the reverse process is unexpectedly forbidden. Our finding could open an avenue for exploring topological photonics that enables nontrivial topological phenomena and designs in photonics driven by nonparaxiality.

1 Different designs of waveguide structures

Control of coupling coefficient
We should address the waveguide array configuration to mimic a tight-binding Schrodinger equation in optically coupled waveguides in two aspects: the coupling coefficients (κ) between waveguides and the on-site propagation constants (β) of each waveguide. First, compared with the control of the propagation constant (see the next section), the control of the coupling coefficient profiles is straightforward. Since the coupling strength between two guiding modes exponentially decays with the distance between two waveguides, one can change the coupling coefficients by varying the waveguide spacing. Fig. S1 shows the coupling coefficient as a function of the distance, which indicates that the coupling coefficient can be tuned in a wide range. Figure S1: Coupling coefficient as a function of spacing between two waveguides.

Control of propagation constant
Second, to control the propagation constant, we proposed five types of flexible designs based on "H-bar" units by varying height structure (h), periodic structure (p), width structure (w ), up-down shifting structure (s), and rotating structure (a). The existing spoof surface plasmon polaritons (SPP) on the ultrathin corrugated metallic waveguide hold similar dispersion characteristics as SPP in optical frequency, whose dispersion curves are significantly deviated from the light line as shown in Figs. S2(a,b,c,g,h). The dispersion curves are calculated by the finite-difference time-domain CST studio method. Especially in these five types of designs, the dependence of the parameters height (h), width (w ), and period (p) in the "H-bar" structure strongly affects the propagation constant (on-site potential) of the eigenmode, as shown in Figs. S2(d,e,f). Hence, we only consider the varying propagation constant via modulation of the height, width, or period parameters in the topological pumping simulation experiments. Figure S2: Dispersion curves and propagation constants in different designs of waveguide structures. (a,b,c,g,h) Dispersion curves for varying height structure, periodic structure, width structure, up-down shifting structure, and rotating structure, respectively. (d,e,f) propagation constants as functions of height, period, and width, respectively.
While the waveguide spacings are maintained the same to modulate the onsite potentials, we also arranged the structure by varying the width and period of "H-bar" units, as shown in Fig. S3(b) and S3(c), respectively. In the modulation of the period, although the propagation constant departs from the sine function, it only affects the micromotion dynamics of the pumping process. However, it maintains the topological nature of pumping process owing to the topological invariant staying the same.

Precise coupling coefficient in Rice-Mele model
Strictly speaking, the coupling coefficients depend on not only the spacings, but also relate to the geometry of waveguide structures. But in the practical design, the coupling coefficients are usually considered as only depending on the spacing because of the nearly exponential decay with spacing and the slow change with structural parameters. To precisely describe the coupling coefficient, we numerically calculate the coupling coefficients between the first and second waveguides with the same parameters at different eight positions of ϕ=π/8, 2π/8,...7π/8 and π in the Rice-Mele model waveguide array. When injecting microwave from one waveguide, the microwave will oscillate between the two waveguides and the coupling coefficient (κ) can be derived from the length (L) of the oscillations by the relation κ × L = π 2 . By comparing the obtained value with the sine function considering only the spacings, we find that there is small difference between the coupling coefficient used in the manuscript and the updated coupling coefficient, as shown by the star marks on the curve of Fig. S4. Figure S4: Comparison of coupling coefficients with considering spacings between waveguides (red solid line) and with simultaneously considering both spacings and structure parameters (marked as star).

Experimental platform
We carried out experiments utilizing ultrathin bending metallic arrays of coupled "H-shape" waveguides, which support SSPPs at microwave to infrared wavelength. There are three designs for the "H-shaped" waveguide structure to obtain the modulation of the on-site potential. The sample fabricated on a printed circuit board (F4BK) with dielectric constant 2.65, loss tangent 0.001, and thickness 0.2 mm is shown in Fig. S5. We conducted the near-field measurements using a near electric-field platform [see Fig. S5], which is composed of a Vector Network Analyzer (VNA) (E5063A), a monopole antenna as the detector, and a translation stages that can move in the x -and y-directions automatically controlled by a stepper motor. The metallic waveguide is connected to port 1 of the VNA through the SMA to feed the energy, and the end is adhered with an absorber to eliminate reflections. To probe the vertical electric component (E z ) of the SSPPs, we fixed the monopole antenna on top of the metallic waveguide around 1.4 mm and connected it to port 2 of the VNA for recording the field distributions. Figure S5: Experimental platform that consists of stages controller, VNA (E5063A), SMA port, sample and translation stages.

Comparison between simulation and experiment results
By solving the original non-paraxial Helmholtz equation, we give a direct numerical software (CST Studio) simulations of the light propagation in the three different modulated waveguide arrays, without relying on coupled mode theory under a paraxial approximation. The numerical simulations of light propagation injected from left boundary and right boundary waveguides are respectively shown in left-top and right-top panels of Fig. S6 (a)-(c), corresponding to heightvarying, width-varying and period-varying designs. The parameters are chosen as the same as those in the experimental designs. For comparison, we also present the corresponding experiment results in the bottom panels of Fig. S6 (a)-(c). All the numerical simulations agree well with the experimental ones. Comparison with the symmetric topological pumping of light with a much higher frequency (λ = 1.55µm, the corresponding frequency f = 194THz) in Fig. 4, these generic asymmetric topological pumping of microwave with a lower frequency 17GHz indicate that the asymmetric behaviors come from non-paraxial effect of microwave.

Direct drivation of negative NNN coupling coefficient
From the experimental results, it is not difficult to find that most of the electromagnetic field tunnels to the third waveguide at the position of ϕ=π/2 when excited from the first waveguide. But most of the electromagnetic field is still on the tenth waveguide at the position of ϕ=π/2 when excited from the tenth waveguide. We respectively take out the structural parameters of the first, second and third waveguide, as well as the eighth, ninth and tenth waveguide at the position of ϕ=π/2, and build three straight waveguides structure related to these parameters. Based on the simulation results in three coupled waveguides, we can acquire the negative NNN coupling (κ NNN ) induced by nonparaxiality, as shown in Fig. S7. κ ef f =κ 1,3 + κ1,2κ2,3 β1−β2 or κ ef f =κ 8,10 + κ8,9κ9,10 β8−β9 represents the effective coupling coefficient in the three coupled waveguide array, whereκ 1,3 = κ 1,3 +κ NNN ,κ 8,10 = κ 8,10 +κ NNN , κ 1,3 ≃ 0.0019mm −1 (the spacing between the first and third waveguides is 9.67 mm), κ 8,10 ≃ 0.0028mm −1 (the spacing between the eight and tenth waveguides is 8.7 mm), κ 1,2 = κ 2,3 = κ 8,9 = κ 9,10 ≃ 0.021mm −1 , β 1 = β 3 = β 9 ≃ 0.44mm −1 and β 2 = β 8 = β 10 ≃ 0.5mm −1 . The effective coupling length L ef f of the two kinds of three coupled waveguides are extracted as 135mm and 440mm, respectively. From the relation κ ef f × L ef f = π 2 , we can easily calculate κ ef f to be -0.01164m −1 and 0.00357m −1 , respectively. Attention should be payed to the sign of κ ef f , otherwise the opposite results will be obtained, that are, coupling phenomenon is suppressed in the first-second-third waveguide and enhanced in the eighth-ninth-tenth waveguide structures. By using one of the κ ef f ,κ NNN turns out to bẽ κ NNN = −0.00619mm −1 orκ NNN = −0.00658mm −1 , in agreement with the our fitting value.

Deviation of Eq. (3) from Eq. (2) in the main text
The electromagnetic waves are confined on the surface (y direction) of the microwave waveguides, hence we can start from the two-dimensional Helmholtz equation for the (x, z)-plane monochromatic electromagnetic waves, Here, k is the wave vector, and the refractive index n(x, z) is designed according to the Rice-Mele type model. Substituting the in-plane electromagnetic field E(x, z; t) = ψ(x, z)e iωt−ikn0z into the Helmholtz equation, the above equation can be simplified as The above formula can be written as Eq.
(2) in the main text, . Next, we replace the Helmholtz-Hamiltonian operator by the structurally modelled Rice-Mele Hamiltinian H RM in our system. Along the transverse direction, we can expand the electromagnetic field in the basis of Wannier functions w j (x, z) as, where the Wannier functions w j (x, z) are real localized functions centered at the jth waveguides, satisfying the orthogonal relation Substituting Eq. (7) into Eq.
(3), we can obtain Projecting the above equation into Wannier functions w l (x, z) and integrating over x, we can obtain Here, we need to evaluate the magnitudes of w l (x, z) ∂ 2 ∂z 2 w j (x, z) and w l (x, z) ∂ ∂z w j (x, z). For the waveguides designed as Rice-Mele type model, we can approximate Wannier functions as localized Gaussian functions Here, X j (x, z) = jd + νd cos(πj + 2πz/Λ) are the positions of the jth waveguides. ν is the modulation strength of waveguide spacings with an average value d. σ is the width of Gaussian function whose dependence on z can be effectively neglected. For the w l (x, z) ∂ ∂z w j (x, z) term, we have For the w l (x, z) ∂ ∂z w j (x, z) term, we have This term does not depend on the waveguides and contributes to an overall shift which can be neglected. At last, we can obtain a model that describes the coupling between φ j , Here, dxw l (x, z)Hw j (x, z) are the matrix elements of the Helmholtz-Hamiltonian operator by the structurally modelled Rice-Mele Hamiltinian H RM . The above equation can be further rewritten as Then we derive the effective Schrödinger-type form in the main text This equation can be iteratively solved to yield the nonparaxial beam propagation of the electromagnetic field.

Pade approximations
In the conventional Schrödinger-type equation, i ∂ ∂z is directly given by H RM . However, in our effective Schrödinger-type formula [Eq. (3) in the main text], the operation i ∂ ∂z is related to itself. Such equation can be solved in a self-consistent way, that is, we iteratively substitute i ∂ ∂z in the denominator by i ∂ ∂z and use the initial condition i ∂ ∂z −1 = 0. For example, the zeroth-order of i ∂ ∂z is given by the first-order of i ∂ ∂z is given by and the second-order of i ∂ ∂z is given by The right-hand side of the above equations can be written as an effective Hamiltonian H n,m ef f = N (n)/D(m) where N (n) and D(m) are polynomials of H RM , which are called as Padé approximants. The first few orders of Padé approximant are given in the Table 1 for details.
The effective Hamiltonians H n,m ef f are still Hermitian. Here, we give a simple mathematical proof. We can expand the Hermitian Rice-Mele Hamiltonian as H RM = j κ j |ψ j ⟩⟨ψ j |, where κ j are real and |ψ j ⟩⟨ψ j | are Hermitian. The effective Hamiltonian as a real function of Rice-Mele Hamiltonian, H ef f = f (H RM ), can be written as H ef f = j f (κ j )|ψ j ⟩⟨ψ j |. Because f (κ j ) are real and |ψ j ⟩⟨ψ j | are Hermitian, the effective Hamiltonian H ef f is also Hermitian.
By numerical examination of the matrix element of the (1,0) Padé approximant, we find that the next-nearest-neighboring (NNN) coupling strengths are always negative. To give a comprehensive understanding, for the (1,0) Padé approximant, we make a Taylor expansion of 1/(1 + H RM /(2kn 0 )) around H RM = 0 up to the first order, and give an approximation of the effective Hamiltonian, The term − 1 2kn0 H 2 RM contributes an effective negative NNN coupling which is approximated as κ N N N ≈ −κ 2 N N /2kn 0 . From Eq. (19), we can immediately know that H ef f is Hermitian, consistent with the general conclusion mentioned above.
(m, n) Figure S8: Planar waveguide array and electromagnetic component of TM0 distributions in the coupled waveguide array. The spatial TM0 distribution of the input and output portions indicate that the left-to-right pumping is symmetric to the right-to-left pumping.

Decreasing the strength of κ NN
We increase the spacing between the waveguides to decrease the nearest-neighbor coupling κ NN . Four cases were compared in the simulation results: (a) g min =1mm and g max =4.2 mm, (b) g min =2mm and g max =5.2mm, (c) g min =3 mm and (d) g max =6.2mm, g min =4 mm and g max =7.2mm. We find that the asymmetry behavior becomes weaker as the waveguide spacing increases. This is because the decrease of nearest-neighbor coupling leads to the dramatical decrease of next nearest-neighbor coupling, consistent with the theoretical prediction (κ N N N ≈ −κ 2 N N /(2kn 0 )). Figure S9: With the decrease of the nearest-neighbor coupling coefficient, the simulation results compare the light beam propagation under the four parameters. (a) g min =1mm and g max =4.2 mm. (b) g min =2mm and g max =5.2mm. (c) g min =3 mm and g max =6.2mm. (d) g min =4 mm and g max =7.2mm.

Different total lengths of waveguide arrays
To exclude the nonadiabatic driving as the reason for the asymmetry, we have performed simulations of three Rice-Mele waveguide arrays with increased total lengths, 600mm, and 800mm, respectively. Since the total length plays the role of the modulation period T in the Rice-Mele model, a larger total length corresponds to a smaller modulation frequency, thus rendering the system more "adiabatic". As shown in Fig. S10, asymmetric pumping can be clearly seen in all of the two simulation results. As a result, we can reasonably rule out the nonadiabatic driving as an underlying factor for the asymmetry. Figure S10: Rice-Mele waveguides arrays with increased total length in the propagation direction (z), namely, (a) 60cm, and (b) 80cm. Note that asymmetric pumping can be clearly seen in all two cases.