1 Introduction

Zitterbewegung (ZB) oscillation represents a high frequency trembling motion of free relativistic Dirac particles described by the Dirac equation, which is caused by the interference between the positive and negative energy states [1], [2]. This phenomenon exhibits an oscillatory motion of a propagating wavepacket transverse to its ballistic trajectory in the absence of external fields. However, a direct experimental observation of ZB effect is difficult owing to its extremely high frequency and small amplitude [2], [3]. Due to the similarity between the relativistic electron and two interacting bands in solid physics, ZB phenomena have so far predicted in a wide variety of quantum and classical physical systems, such as semiconductors [4], [5], superconductors [6], [7], topological insulators [8], Bose–Einstein condensate [9], ultracold atoms [10], plasma [11], as well as optical systems [12], [13], [14], [15], [16], [18], [19]. The experimental observation of quantum phenomena in the electron case is severely hindered by the complex interactions and the limited scale. No such difficulties exist in the optical context, the use of photonic systems to mimic quantum phenomena have attracted continuous and increasing interest from both fundamental and applied viewpoints. More specifically, optical analogs of the relativistic ZB effect have recently been investigated in modulated photonic lattices, which provide a powerful tool for the control of the fundamental aspects of light propagation. In binary waveguide array, the two-component (spinor) wave function dynamics of the Dirac equation can be simulated by paraxial light propagation and the spinor is represented by two interleaved sublattices [12], [13], [14]. In two-dimensional photonic crystals, photonic analogs of the ZB can be observed by measuring the time dependence of the transmission coefficient through photonic crystal slabs [15], [16]. ZB effect is also demonstrated in a honeycomb lattice of coupled rods with slowing varying radii along the direction of propagation, which possesses photonic Weyl points and become Dirac points when projected down to a two-dimensional momentum space [17]. Furthermore, in photonic microcavities, an analogue ZB resulting from the photonic spin–orbit coupling on the real space propagation of polariton wavepackets is also experimentally demonstrated [18]. Another optical simulation of ZB is based on frequency-conversion process of optical pulses in quadratic media, the mismatch of the group velocities of signal and sum-frequency fields leads to an oscillating motion [19].

As mentioned above, optical ZB are most frequently considered in periodic structures with constant transverse modulation by changing the waveguide strength and spacing or the background refractive index. Nevertheless, the periodically driving-induced behaviors of light in space–time periodic optical system can be understood by directly inspecting the band structure in the Floquet picture [20], [21]. Floquet engineering provides additional possibilities for control of beam propagation [22]. Among them, photonic lattice composed of evanescently coupled single mode helical waveguides, fabricated using femtosecond direct laser writing, has received considerable attention for its topological properties and band collapse. The optical analogues of some quantum phenomena originally predicted for electrons have been realized in array of helical waveguides, including Bloch oscillations [23], [24], Zener tunneling [25], conical diffraction and dynamic localization [26]. However, the study of the coexistence and interplay between the different modes remains largely unexplored in the helical waveguide arrays. It is interesting to investigate the interplay between the Floquet mechanism and the inner state interference.

In this work, we show both theoretically and numerically that the optical ZB effect can be realized in Floquet photonic lattices composed of coupled helical waveguides. By employing an alternation of distances between adjacent waveguides, the quasi-energy bands of the superlattice can be effectively controlled, complete collapse and band curvature inversion are obtained. That affects the amplitude, frequency and direction of ZB effect. Even they can be completely arrested or their direction can be reversed. ZB-like motion with the lateral shift is also demonstrated in two-dimensional Floquet photonic lattices composed of coupled vertical and horizontal stacks of Su-Schrieffer-Heeger (SSH) [27] chains. Both completely arrest and inversion of lateral shift are also observed.

2 Optical ZB effect in one-dimensional helical waveguides

An example of scheme of one-dimensional setup is shown in Figure 1(a), which is somewhat similar to SSH lattice. A period array composed of coupled dimers. d₁ is the spacing between sites inside dimers, and d₂ is the spacing between two neighboured dimers. The governing equation, for the complex amplitude of the field ψ(x, y, z) of a light beam propagating in a helical waveguide array with the helix spatial period Z and radius R, can be derived from the Maxwell equations. By using the slowly amplitude approximation, the following two-dimensional Schrödinger-type equation can be arrived:

where ψ(x, y, z) is the envelope of the electric field E(x, y, z) = ψ(x, y, z)exp(ik₀z − iω₀t), k₀ = 2πn₀/λ is the wavenumber, λ is the wavelength, n₀ is the background refractive index, ω₀ = 2πc _v /λ is the frequency with c _v being the light speed in the vacuum, Δn(x, y, z) is the refractive change that define the helical waveguides. To change the z-dependency of the waveguide system into a stationary and straight one, we do a coordinate transform operation through the relations x′ = x + R cos(Ωz), y′ = y + R sin(Ωz) and z′ = z, where Ω = 2π/Z is longitudinal frequency of the helix corresponding to periodicity Z. Now the Eq. (1) is written as [21], [24], [28]:

where ψ′ = ψ′(x′, y′, z′), ▽ = (∂x′, ∂y′); A(z′) = A₀[sin(Ωz′), −cos(Ωz′), 0] is the artificial vector potential induced by the helicity of the waveguides with amplitude A₀ = k₀RΩ. The propagation of light in such weakly coupled waveguides with the nearest-neighbor interaction can be described by the tight-binding model:

Figure 1:

The quasi-energy spectrum of one-dimensional array of helical waveguides. (a) Schematic of one-dimensional array of helical waveguides. Each dimer is distinguished by the index n. (b) Floquet quasi-energy spectrum (red curve) and dispersion relation obtained by the analytic treatment (blue curve) at d₁/d₂ = 0.8, black dotted curves show the corresponding hyperbolic dispersion curves. (c) Effective coupling constant c_eff1 as a function of the value of d₁/d₂. (d) Dispersion relations of the lattice with d₁/d₂ = 0.8 (red curve), 1, 1.66 (blue curve) and 2 (green curve). Other parameters are Z = 0.5?cm, R = 8?μm, n₀ = 1.45, d₂ = 10?μm, λ = 633?nm, c₀ = 60?m⁻¹ and b = 0.03.

To simplify notation, we have omitted the superscript ′. ψ_2n and ψ_2n+1 denote the amplitude functions linked with left and right waveguides inside dimers. r is the displacement vector between the neighboring waveguides, c(d_1,2) is the coupling strength, and assumed to be decay exponentially with the distance c = c₀?exp(−bd) [29], [30], where b characters the exponential decay rate.

It is instructive to analyze the band structure of one-dimensional helical waveguides. Due to the periodicity of the Hamiltonian H(q _x , z) = H(q _x , z + Z) [25] where q _x is transverse Bloch momentum, Floquet theory can be applied to drive a band structure of so-called quasi-energies E. A solution of such z-dependent coupled mode equations is the Floquet states ψ _n (z) = ? _n (z) exp(iEz + iq _x nd) with d = d₁ + d₂ and ? _n (z) = ? _n (z + Z). The quasi-energies can be calculated directly by . A graphic presentation of the Floquet quasi-energy spectrum is shown in Figure 1(b). Two minibands are formed, which intersect at the edges of the Brillouin zone for d₁ = d₂, whereas the two minibands are well separated energetically when c(d₁) ≠ c(d₂). It is shown that, near the Brillouin zone edge, the dispersion relations approximate two opposite hyperbolas, as shown in Figure 1(b), which is similar to that for positive-energy and negative-energy branches of a freely moving relativistic massive particle, therefore beam propagation at incidence angles near the Bragg angle can mimic the temporal dynamic of a one-dimensional massive Dirac electron, that is to say, our model enables an easy visualization in space of the photonic ZB phenomenon.

To verify the above analysis, an analytic treatment of the light dynamics is developed, which can also provide some valuable insight. Discrete diffraction in arrays of helical waveguides occurs like in a straight array with an effective coupling after every period, which can be computed using averaging methods and is given by c_eff1 = c(d₁)J₀(k₀RΩd₁) and c_eff2 = c(d₂)J₀(k₀RΩd₂) for the spacing between adjacent lattice sites d₁ and d₂, respectively [24], where J₀ is the Bessel function of the first kind. Coupled-mode equations then take the following simplified form:

The corresponding band structure reads

which is entirely consistent with Floquet quasi-energy spectrum of Figure 1(b) (the red and blue curves coincide). Let us now analyze the spectrum of the photonic lattice in detail. Equation (5) indicates that the effective couplings can be tuned from positive to negative by varying the waveguide separation, and flat band appears at d_1cr = 1.66d₂, corresponding to zero of c_eff1, as shown in Figure 1(c), the band structure experiences a complete band collapse in the first Brillouin zone (see Figure 1(d)) and would support compact localized modes. A further increase of the distance d₁, results in inversion of the respective dispersion curves. Obviously, the effect of band collapse would be repeat as d₁ increases. Around |q _x | = π/d, the two minibands are separated by a gap β = 2|c_eff1 − c_eff2|, and exhibit two opposite hyperbolas (see Figure 1(b)), where Q = dq _x − π. Let us assume that the array is excited by a broad beam G(x) tilted at the Bragg angle, the modes are excited with a nearly equal amplitude, but the phase difference between adjacent waveguides is equal to π/2 across the array. By setting ψ_2n(z) = (−1) ⁿ φ₁(n, z), ψ_2n−1(z) = −i(−1) ⁿ φ₂(n, z) and introducing the continuous transverse coordinate ξ ↔ n = x/d, the two-component spinor satisfies the following equation:

where σ_x,y are the Pauli matrices. Equation (6) corresponds to the one-dimensional Dirac equation for an electron of mass m₀ provided that the formal change c_eff1 ↔ c _v , , z ↔ t is made, where ? is the reduced Planck constant. Obviously, the Dirac equation in its discretized form can be mapped on an array of helical waveguides with periodically alternating transverse lattice spacing. The Dirac spinor wave function can also be mapped onto the spatial evolution of the field amplitudes φ₁ and φ₂. As a result, ZB effect is observed as a quivering spatial oscillatory motion of the beam center of mass. By solving this equation in the momentum k space, one gets:

where . The expectation value of position for the relativistic particle is defined by , which can be used to determine the spatial motion of the beam center through the relation ⟨n⟩(z) ? 2⟨ξ⟩(z) + 1/2. ⟨ξ⟩(z) is calculated as

where .

The main result of our analysis is summarized by dividing the evolution of the expectation value of the position into the drift and ZB component. The third term denotes the ZB oscillation. The frequency of ZB, which can be estimated as ω = 2|?| = 2|c_eff2 − c_eff1|, grows linearly with increasing a gap β = ω. It can be seen clearly that the motion of the wave packet is accompanied by a transverse drift at a constant velocity v₀ and the direction of drift depends on the polarity of c_eff1, which can be effectively regulated through properly adjusting the lattice distance d₁.

Now, going back to one-dimensional helical waveguides (Figure 1(a)), and examining the mode interference of the system. In order to reduce the influence of the boundary effect, the model consists of more than 400 dimers. The wave transmission for a broad Gaussian beam excitation (corresponding to Bragg angle λ/(2dn₀)) is investigated by directly solving Eq. (3), where w₀ = 10d. Beam dynamics reveals that a clear trembling motion of the beam – corresponding to the optical ZB oscillation – is observed, as shown in Figure 2(a1). One immediately notes the damping of ZB oscillations is also visible, which is ascribable to the angular spectral broadening of the launching broad Gaussian beam. The corresponding oscillation period along the propagation is T = 0.64?m illustrated in Figure 2(b1), which is in good agreement with the theoretical predictions (T = 2π/β). The influence of d₁ on the effective coupling strength c_eff1 is summarized in Figure 1(c). With increasing the distance d₁ (d₁ > d₂), the gap is extended, and the oscillation period increases and the amplitude decreases. In fact, the initial oscillation amplitude can be estimated in the limit at k = 0. Based on Eq. (8), the frequency and amplitude of ZB effect are given by ω_ZB = 2|c_eff1 − c_eff2| and R_ZB = πc_eff1/ω_ZB, respectively. The analytical and numerical results agree with each other very well.

Figure 2:

Optical ZB supported by one-dimensional array of helical waveguides with d₁/d₂ = 0.8, 1.66 and 2, respectively (from top to bottom). Left and right panels show the beam propagation and the evolution of center of energy, respectively. The black solid and red dotted curves correspond to the direct numerical simulated of Eq. (3) and analytical results based on Eq. (8), respectively.

The lattice spacing affects not only amplitude and frequency of the optical ZB, but it can also arrest this effect and even invert their direction of drift. When the intracell spacing d₁ is small, oscillations occur in the region n > 0, the lateral shift of the center of beam energy decreases with increase of d₁. At d₁ = 1.66d₂, corresponding to the special case c_eff1 = 0, two flat band pinned at E = ±c_eff2 are presented. This means that excitation in such a lattice will propagate without diffraction, which results in complete localization of the wavepacket in the system. This effect is known as dynamic localization. The ZB oscillations are undamped and they are localized, as shown in Figure 2(a2) and (b2). Further increasing the distance d₁, corresponding to negative c_eff1, subsequent reappearance of ZB effect at d₁ > d _1cr is accompanied by the reversal of the direction of oscillations, where ZB effect appears in the region n < 0 (see Figure 2(a3) and (b3)).

We further investigate the robustness of the optical ZB in arrays of helical waveguides, by introducing the random noise in the coupling strength c(d_1,2). A coefficient (1 + δr _n ) is multiplied on coupling strength, where δ is the perturbation strength, r _n is a random number between −1 and 1. Figure 3 shows that, the beams exhibit the similar ZB-like lateral shifts with δ up to 10 % (namely, δ = 0.1). Consequently, they show a great robustness against perturbation of the coupling strength.

Figure 3:

Robustness of optical ZB with d₁/d₂ = 0.8 under the random perturbation. (a) Light evolution at δ = 0.1. (b) Center of the energy with different perturbation.

Here, we focus on the impact of lattice spacing on the ZB effect in arrays of helical waveguides, it is shown that longitudinal rotation of waveguides (helix radius and period) leads to notable variations of Floquet quasi-energy spectrum, that will certainly drastically affect the ZB oscillations. Figure 4(a) shows the effective coupling constants as a function of the helix spatial period. The ZB oscillations can occur in the regions n > 0 and n < 0 at Z = 0.40?cm and 0.23?cm, respectively. ZB effect with undamped amplitude is also observed at Z = 0.31?cm.

Figure 4:

Optical ZB for different helix spatial periods. (a) Effective coupling constants c_eff1 (solid curve) and c_eff2 (dotted curve) as a function of the helix spatial period Z with d₁/d₂ = 0.4. (b) Lateral shifts of the center of energy at Z = 0.40?cm (red curve), 0.31?cm (black curve) and 0.23?cm (blue curve). Other parameters are same as that in Figure 1.

3 Optical ZB effect in two-dimensional helical waveguides

The ZB-like motion is further investigated in two-dimensional geometry shown in Figure 5(a), which is based on two-dimensional SSH model composed of coupled vertical and horizontal stacks of SSH chains [31]. (m, n) are the integers that enumerate the lattice sites. The lattice consists of four sites per unit cell with dimerized nearest neighbor couplings. Light propagation can be described by the following coupled-mode equation for the modal field amplitudes ψ_2n,2m,

Figure 5:

Dispersion relation of two-dimensional helical waveguide array. (a) Schematic of two-dimensional helical waveguide array. (b) Corresponding first Brillouin zone of the square lattice in (a), (c) dispersion relation of two-dimensional model with d₁/d₂ = 0.8 along high symmetry points. Other parameters are same as Figure 1.

After setting ψ_2n,2m = (−1)^n+mφ₁(n, m), ψ_2n−1,2m = −i(−1)^n+mφ₂(n, m), ψ_2n,2m−1 = −i(−1)^n+mφ₃(n, m) and ψ_{2n−1,2m−1} = (−1)^n+mφ₄(n, m), and introducing the two continuous transverse coordinates η ↔ n, ζ ↔ m, the four-component bispinor satisfies the following two-dimensional Dirac-type equation:

where, α, ρ and τ are 4 × 4 matrices with all elements equal to zero, except α_1,2 = α_2,1 = α_3,4 = α_4,3 = 1, ρ_1,3 = ρ_2,4 = ρ_3,1 = ρ_4,2 = 1 and τ_1,2 = τ_1,3 = −τ_2,1 = τ_2,4 = −τ_3,1 = τ_3,4 = −τ_4,2 = −τ_4,3 = 1. Correspondingly, two-dimensional optical ZB effect can also be realized in arrays of helical waveguides.

Substituting the periodic solution ψ_n,m(z) = ? _n (z) exp[iEz + i(q _x n + q _y m)d] in the form of Floquet–Bloch-like wave in Eq. (3), where q _x and q _y are the wave numbers along n and m, the dispersion relation of the two-dimensional helical waveguide array along high symmetric points (Figure 5(b)) is shown in Figure 5(c). At the high symmetry point M, there are three eigenvalues E₁ = −E₃ = 2|c_eff2 − c_eff1|, and E₂ = 0. Similar to the case of one-dimensional geometry, the waveguide array is also excited by a broad circular Gaussian beam, with w₀ = 6d. Figure 6 shows the lateral ZB-like motion of the Gaussian beam. Although there contain two oscillation frequencies generated by interferences between the upper and zero modes (E₁ − E₂), and between the upper and low modes (E₁ − E₃). The former plays the leading role in determined the ZB-like motion. The oscillation period is determined by the T = 2π/E₃. Figure 6(a) show the lateral shifts of the Gaussian beams with d₁/d₂ = 1.5, 1.66 and 2.2. At c_eff1 = 0 (d₁/d₂ = 1.66), corresponding to dynamic localization, the beams oscillation without lateral shift shown in Figure 6(c). When d₁/d₂ = 1.5 (for c_eff1 > 0), the wave packet expands in both transverse dimensions, optical ZB-like motion occurs with the lateral shift along the bisectrix of the (n, m) plane in the first quadrant(n, m) > 0, as shown in Figure 6(b) The inversion of the direction of this optical ZB-like motion is also confirmed with the lateral shifts in the third quadrant (n, m) < 0 when d₁/d₂ = 2.2 (for c_eff1 < 0), as shown in Figure 6(d). The lateral shifts of the two beams are opposite shown in Figure 6(a). The corresponding oscillation periods are T = 0.15?m, 0.13?m and 0.10?m, respectively, which are agreement with the theoretical predicated.

Figure 6:

Optical ZB effect in two-dimensional helical waveguides. (a) Lateral shifts of the center of the Gaussian beam in two-dimensional geometries at d₁/d₂ = 1.5 (red line), 1.66 (black line) and 2.2 (blue line). (b)–(d) Field modulus distributions at z = 0.92?m for d₁/d₂ = 1.5, 1.66 and 2.2, respectively.

4 Conclusions

In conclusion, propagation of light waves in one-dimensional and two-dimensional Floquet photonic system made of SSH-like helical arrays of evanescently coupled optical waveguides is theoretically and numerically investigated and shown to provide a classic wave optics analog of the ZB effect, where a spatial oscillatory motion of an optical beam is observed. In one-dimensional case, Floquet quasi-energy spectrum consists of two minibands. Near the edges of the first Brillouin zone, the dispersion curve looks like that for positive-energy and negative-energy branches of a freely moving relativistic massive particle and light transport simulates the temporal dynamics (ZB effect) of the relativistic Dirac equation. The lattice spacing determines the effective coupling strength, and further drastically not only affects the oscillation amplitude and frequency, but also invert their direction of drift when the effective coupling strength is tuned from positive to negative. The ZB-like motion is also demonstrated in two-dimensional geometry. At the high symmetric point M, there are three eigenvalues, the interference between the upper and zero modes plays the leading role in determined the ZB-like motion. Their oscillation and direction strongly depend on waveguide spacing (coupling strength). Complete arrest and inversion of the drift direction of ZB effect are reported. Robustness of the optical ZB phenomenon to the perturbations of coupling strengths is also demonstrated. Here, we focus on the impact of lattice spacing on the ZB effect in arrays of helical waveguides, it is shown that longitudinal rotation of waveguides (helix radius and period) leads to notable variations of Floquet quasi-energy spectrum, that can certainly drastically affect the ZB oscillations. Floquet engineering provides additional possibilities for control of ZB effect. Our results will not only enrich the ZB research in Floquet waveguide arrays but also provide a new opportunity to manipulate light propagation. Furthermore, due to its easy accessibility and high controllability, we hope the Floquet photonic lattices composed of coupled helical waveguides enables an easy visualization of other optical analogues of relativistic quantum phenomena (such as Klein tunneling).