Relativistic Zitterbewegung in non-Hermitian photonic waveguide systems

Zitterbewegung (ZB) is a phenomenon in relativistic quantum systems where the electron wave packet exhibits a trembling or oscillating behavior during its motion, caused by its interaction or coupling with the negative energy state. To directly observe ZB in electronic systems is difficult, due to the challenges associated with the atomic scale wavelength of the electron. Photonic systems offer an alternative paradigm. We exploit the concept of pseudo parity-time (pseudo PT ) symmetry to study ZB in non-Hermitian quantum systems implemented as an experimentally feasible optical waveguide array. In particular, the non-Hermitian Hamiltonian is realized through evanescent coupling among the waveguides to form a one-dimensional lattice with periodic modulations in gain and loss along the guiding direction. As the modulation frequency is changed, we obtain a number of phenomena including periodically suppressed ZB trembling, spatial energy localization, and Hermitian-like ZB oscillations. We calculate phase diagrams indicating the emergence of different types of dynamical behaviors of the relativistic non-Hermitian quantum system in an experimentally justified parameter space. We provide numerical results and a physical analysis to explain the distinct dynamical behaviors revealed by the phase diagrams. Our findings provide a deeper understanding of both the relativistic ZB phenomenon and non-Hermitian pseudo-PT systems, with potential applications in controlling/harnessing light propagation in waveguide-based optical systems.


Introduction
There has been a great deal of recent interest in investigating the role of parity-time reversal (PT ) symmetry in wave propagation at different scales with the discoveries of phenomena such as non-reciprocal beam propagation [1][2][3][4][5], and uni-directionally transparent invisibility [6][7][8]. The phenomenon of topologically protected states [9,10] was originally discovered in condensed matter physics associated with electronic transport, but recently it has been demonstrated in optics [11][12][13][14][15][16][17], due naturally to the correspondence between matter and optical waves. For example, topological photonics/acoustics were demonstrated by exploiting the analogy between electronic and synthetic photonic crystals, where tunable and topologically protected excitations were observed [18,19]. There were also efforts in exploring new wave features in artificial photonic crystals with/without the PT -symmetry [20][21][22][23][24][25][26][27]. All these have led to the emergence of a forefront area of research in optics: light propagation in non-Hermitian PT -symmetric media with balanced loss and gain profiles. The new field offers the possibility to engineer light propagation, potentially revolutionizing optics with unconventional applications.
In conventional quantum mechanics, the observable operators are required to be Hermitian to ensure real eigenvalues. This requirement is the result of one of the fundamental postulates in quantum mechanics: the physically observable or measurable quantities are the eigenvalues of the corresponding operators. However, the seminal works of Bender et al. [28][29][30] demonstrated that non-Hermitian Hamiltonians are also capable of generating a purely real eigenvalue spectrum and therefore are physically meaningful, if the underlying system possesses a PT symmetry. This opened a new research field called non-Hermitian PT symmetric quantum mechanics. In fact, in physical systems, non-Hermitian Hamiltonians exist in contexts such as electronic transport (in open Hamiltonian systems) and gain/loss materials in optics [31]. Especially, in optics, materials with a complex index of refraction can effectively be a non-Hermitian PT symmetric system [32]. However, in such a case, the PT symmetry constrictions require that the real and the imaginary parts of the refractive index be an even and odd function in space, respectively [33][34][35], which may be challenging to be fabricated for experimental study. An alternative configuration that does not require even/odd spatial functions is the pseudo-PT symmetric system with periodic modulations [36]. While such a system does not conserve the energy at any given instant of time because of gain/loss, the energy does not diverge within any practically long time.
An outstanding problem is how relativistic quantum effects manifest themselves in non-Hermitian photonic systems with full or pseudo PT symmetry. This is motivated by the tremendous development of 2D Dirac materials in the past decade such as graphene [37][38][39][40][41][42][43], topological insulators [9,10], molybdenum disulfide (MoS 2 ) [44,45], HITP [Ni 3 (HITP) 2 ] [46], and topological Dirac semimetals [47,48], where the underlying quantum physics is governed by the relativistic Dirac equation. It is thus of interest to investigate experimentally more realizable optical systems hosting relativistic excitations and to study the manifestations of the fundamental phenomena [49] such as the Klein tunneling, Zitterbewegung (ZB), and pseudo spin to exploit their topological origin. In this regard, including gain/loss in synthetic optical systems can generate non-Hermitian relativistic PT -symmetric excitations through engineering the gain and loss in a balanced manner. For example, arranging PT -symmetric couplers periodically provides a unique platform to realize the analogy of the non-Hermitian relativistic quantum systems in optics [50][51][52][53].
In this paper, we focus on a fundamental phenomenon in relativistic quantum mechanics -ZB oscillations. While the optical analog of the relativistic ZB effect in Hermitian quantum systems was previously observed [54], we investigate the manifestations of ZB in non-Hermitian photonic systems. In particular, we consider a binary waveguide array with periodically modulated imaginary refraction index in the guiding direction. This kind of time-periodic gain and loss photonic system is widely used in many theoretical and experimental works [36,[55][56][57][58][59]. The optical waveguide array system, due to its controllable degrees of freedom, has been a paradigm to study a host of fundamental physical phenomena [60] such as the neutrino oscillations [61], Bloch oscillation [62], Zener tunneling [63,64], and Klein tunneling [65]. We configure the system to have a pseudo-PT symmetry, so that it exhibits quasi-stationary light propagation with slowly varying time-averaged total intensity [36]. By calculating a detailed phase diagram in an experimentally meaningful parameter space, we uncover a number of phenomena in the system. In particular, we find that, in certain parameter regime, the wave amplitude tends to diverge. A remarkable phenomenon is that there are parameter regimes in which the dependence of ZB oscillations on the modulation frequency is non-monotonic, where the system exhibits a striking, periodically suppressed ZB effect with revival or intermittent ZB oscillations for low frequencies and a Hermitian-like ZB effect in the high frequency regime. In the intermediate frequency regime, a surprising spatial energy localization behavior emerges. By solving the corresponding time-dependent Dirac equation, we obtain analytic results that provide explanations for the numerically observed, ZB manifested phenomena. The findings have implications. For example, the wave divergence phenomenon may be exploited for applications in optical amplifier and lasing. Intermittent BZ oscillations can potentially lead to a new mechanism to manipulate/control light propagation.
Before describing our results in detail, we remark on the two unique aspects of our work.
First, what results can be considered as constituting a highly nontrivial or unexpected general fundamental insight into the ZB phenomenon for researchers working in this field? Photonic crystals allow researchers to build up classical simulators of quantum systems. For example, through engineering the gain/loss media, non-Hermitian, PTsymmetric and non-relativistic quantum physics has been analogously realized using two sets of evanescently coupled waveguides [32] through non-reciprocal beam dynamics. Most previous studies on the optical analogies of non-Hermitian PT-symmetric quantum systems are nonrelativistic. Quite recently, optical analogies of the Hermitian relativistic Dirac equation have been articulated, such as synthetic photonic graphene and photonic topological insulators. An outstanding issue concerns the PT-symmetry and the related physical effects in synthetic relativistic quantum systems. Existing theoretical proposals [50,52] rely on an exact design of the gain-loss profile or a sophisticated strainbased control scheme for the given gain-loss configuration, leading to stationary, non-Hermitian PT-symmetric systems. Distinct from the existing works, our work reports a dynamical scheme to optically realize time-dependent, pseudo-PT symmetric relativistic Dirac equation, and we find the striking phenomenon of periodically suppressed and revival ZB oscillations. This "intermittent" ZB phenomenon is in sharp contrast to conventional ZB oscillations with a constant amplitude. Another aspect of our work lies in observing a resonance boundary between the pseudo-PT and PT-breaking regions. This kind of resonance boundary has been reported in previous works, between the conventional PT-symmetric region and PT-broken region [55][56][57][58][59]. Here, our results show that the boundary between the pseudo-PT region and PT-breaking region follows the similar behavior.
Second, are there potentially feasible experimental schemes? The answer is affirmative. Experimentally, conventional non-Hermitian PT-symmetry systems require that the real part of the effective potential to be an exactly even while the imaginary part be an exactly odd function in space. Our system design relaxes these requirements and significantly expands the configuration range of the non-Hermitian systems through introducing spatial modulation. The width of ZB trembling is on the order of micrometer and the length of the waveguides is millimeters. Concretely, our system can be realized in a waveguide array configuration similar to that in Ref. [54], where the gain effect can be realized, for example, through doping of dye molecules such as Rhodamine B (emission at 627 nm, excitation at 554 nm -green laser). Loss can be introduced using metals evaporated into the waveguides during the fabrication process.
In Sec. 2, we provide a detailed description of our photonic waveguide array system with pseudo PT symmetry, and derive the underlying Dirac equation. In Sec. 3, we present our main results: phase diagrams, amplitude divergence, regular and revival ZB oscillations, and spatial energy localization. In Sec. 4, we summarize the main results and offer a discussion. Certain details of the analytical derivation based on the relativistic quantum Dirac equation are given in Appendix.

Optical waveguide system with pseudo PT symmetry and description in terms of the Dirac equation
We consider a two-dimensional binary photonic superlattice consisting of two kinds of interleaved waveguides, A and B, as shown in Fig. 1. The distance between two adjacent waveguides is a and the coupling strength between them in the effective tightbinding Hamiltonian is κ. Due to the similarity between the optical and electronic waves, an optical waveguide system can simulate the time evolution of the electronic wavefunction subject to an equivalent potential. If the real parts of the refractive index of the nearest neighboring waveguides have the same mismatch, the waveguide superlattice can simulate the behaviors of the finite-mass Dirac equation. In particular, the dispersion relation of the lattice is a hyperbola about the edge of the first Brillouin zone [54,65], and the A and B sublattices can be regarded as corresponding to the two components of the spinor wavefunction underlying the Dirac equation. When the incident angle of a wave packet is approximately the first Bragg angle, for an initial group velocity close to that of the edge of the first Brillouin zone, the optical system obeys the relativistic Dirac equation. Experimentally, it is difficult to realize a waveguide system with time-dependent refractive index. However, the z direction represents effectively the time dimension. The two-dimensional binary photonic superlattice system is then completely equivalent to a time dependent, one-dimensional relativistic Dirac system, making possible experimental studies of various relativistic quantum effects. We employ the tight-binding approximation to analyze the waveguide superlattice.
The amplitude of the optical field is described by where a n is the amplitude of the optical field in the nth waveguide, κ is the effective coupling constant, and σ is the propagation mismatch. We assume that the real part of the mismatch σ r between the nearest neighboring waveguides is constant. The space dependence of σ(z) is embedded in its imaginary part: σ i (z) = rσ r sin (ωz), where σ r and σ i denote its real and imaginary parts, respectively, with r being the ratio between them that characterizes the strength of gain/loss. The quantity ω is the frequency of spatial gain/loss modulation in the z-direction. Equivalently, ω can be regarded as the frequency of the imaginary part of the time-varying potential for the corresponding quantum system that is fundamentally non-Hermitian due to the imaginary potential. We note that, for a Hermitian system, i.e., a superlattice system with a purely real refractive index, previous works [54,65] showed the emergence of two minibands with the dispersion relation given by ω ± (q) = ± σ 2 + 4κ 2 cos 2 (qa). For the wavevector along the x-direction, q is close to π/(2a) and the quantities ω ± (q) exhibit a hyperbolalike behavior with a gap of 2σ, which is similar to the electron and positron dispersion curves from the finite-mass Dirac equation. For a general non-Hermitian system, it is straightforward to show that the dispersion relation has the same form as that for a Hermitian system (see Appendix), but there are two differences: (a) the quantity σ becomes complex and (b) the system is time dependent due to the external driving.
In the limit of small imaginary part, the complex frequency in the dispersion relation possesses the following real part ± (σ 2 r − σ 2 i (z)) + 4κ 2 cos 2 (qa), guaranteeing that our non-Hermitian system can still retain the equivalence to the relativistic quantum system for a massive particle at the boundaries of the Brillouin zone. This condition constricts the incident angle of the wave packet to be about θ B ≈ λ/(4n s a), where λ is the wavelength and n s is the refractive index of the substrate. We emphasize that this analogy holds only when the wave packet is located near the Brillouin zone boundary. Analogy in other regions would be affected by the dispersive effect of the waveguides. In fact, many experiments have shown that the analogy can hold up to hundred millimeters [65][66][67].
To demonstrate the equivalence of the waveguide equation to the Dirac equation, we make the following substitutions [54,65,67]: From the periodicity and symmetry of the system, we see that a unit cell from sublattice A combining with one from sublattice B is effectively a unit cell that forms a Dirac spinor. It is necessary to use the dimensionless continuous transverse coordinate ξ ↔ n = x/(2a) to form a spatial derivative operator. With these considerations, we obtain the following equivalent one-dimensional Dirac equation where are the Pauli matrices. These considerations hold when σ is complex. Comparing with the standard Dirac equation, we note the following equivalence: i.e., the coupling constant between the waveguides corresponds to the speed of light while the mismatch between the waveguides is equivalent to the mass of the underlying particle. Indeed, as we can see from the dispersion relation of the superlattice, the band gap for the Brillouin zone is also proportional to σ, corresponding to the dispersion relation of a massive relativistic quantum particle.
In our simulations, we assume that the initial optical field has the form E( The superlattice is made of 200 waveguides separated from each other with the distance a = 16µm, and the effective coupling rate is κ = 0.14mm −1 . The Gaussian beam has the wavelength of λ = 633nm with a spot size about 105µm. The modulation frequency is scaled by the factor ω 0 = 1/a and the lengths are normalized by a.

Emergence of pseudo-PT symmetry
To ascertain the emergence of the pseudo-PT symmetry, we probe into the parameter space of the photonic superlattice system systematically by calculating the phase diagrams. The results are shown in Fig. 2, where the same initial wave packet with a group velocity determined by the relativistic Dirac point is used for different parameter combinations. The evolution is simulated for a relatively long time to ensure that the system has settled into a steady state. Figures 2(a) and 2(b) show the phase diagrams in the r −ω and σ r −ω parameter planes, respectively. The color is coded in terms of the maximum value of the intensity of the optical field for the corresponding parameters. Considering the fact that, in the pure gain regime, the intensity grows exponentially, we use a logarithmic function to rescale the intensity for better visualization. We set a cut-off intensity level beyond which the optical field is deemed to have diverged, corresponding to the PT breaking case. The cut-off criterion we used is E tot /E 0 tot ∼ 10 9 , where E tot is the maximum of the total intensity and E 0 tot is the total intensity when there is no modulation. Even for this relatively high cut-off intensity, there is little change in the shape of the boundary. In most parts of the boundary of the PT breaking phase and the pseudo-PT phase, the change in the color from red to blue are quite rapid. In Figs. 2(a) and 2(b), the divergent regions are coded as dark red. For a small modulation frequency, a number of waveguides absorb energy within the simulation time so that their intensity can exceed the cut-off value. For such waveguides, a pseudo-PT behavior cannot be numerically detected. Because of this, there is a persistently red region in the phase diagrams when the modulation frequency is close to zero. From Fig. 2(a), we see that a pseudo-PT behavior exists either in the parameter region with small gain/loss strength or in the region with a high modulation frequency. The reason is that, when the gain/loss of the waveguides is small, energy absorption is insignificant so that the waveguides will have sufficient time to transfer the energy within the superlattice and dissipate the absorbed energy into the adjacent waveguides. When the modulation frequency is high, the period becomes small so that the waveguides are able to dissipate the absorbed energy through itself when the imaginary part of the refractive index changes its sign during light propagation. In general, the stronger the gain/loss strength, the higher the modulation frequency is needed to balance the gain and loss, providing an explanation for the shape of the boundary of the two regions in the high modulation frequency regime in the phase diagrams. Figure 2(b) exhibits a similar behavior, i.e., the pseudo-PT behavior exists in the larger real refractive index (corresponding to weaker gain and loss effects) and higher modulation frequency regions. A general feature is that higher modulation frequency is favorable for the emergence of the pseudo-PT symmetry. A phenomenon present in both Figs. 2(a,b) is that the boundary of the pseudo-PT symmetric and the non pseudo-PT regions exhibits oscillatory variations as the modulation frequency is changed. For ω/ω 0 ∼ 4 − 6, the system exhibits a divergent behavior even for the small gain/loss strength. To understand the origin of the boundary variations, we examine the details of the time evolutions of the system. In Fig. 2(a), we mark the positions of the valleys in the variation with white lines, which are ω/ω 0 = 0.6, 0.9, 1.5 and 4.2, leading to the ratio of the period, i.e., 1/ω, to be about 7 : 5 : 3 : 1, indicating that the system tends to diverge when the modulation period is odd times of the rightmost line. When the ratio is even, the system tends to stay in the pseudo-PT phase. Here we want to emphasis that this kind of resonance boundary has been extensively studied in the conventional PT-symmetric system [55][56][57][58][59]. However, when the ratio becomes large, it is hard to observe the oscillatory behavior, due mainly to the finite resolution of the simulation. Another reason is that the system tends to shift into the high energy regime when the modulation frequency is low, so the oscillations are buried in the regime that does not correspond to PT breaking. The boundary in Fig. 2(b) exhibits a similar variational behavior. Particularly, as σ r becomes larger, the imaginary part becomes relatively less significant. As a result, the band gap and thus the ZB trembling frequency is given by 2σ r (to be described in Sec. 3.2). When the modulation frequency is equal to the ZB trembling frequency, a resonant effect appears, forming a pseudo-PT phase boundary that starts from ω/ω 0 = 2 in Fig. 2(b). The specific ratio relation suggests the occurrence of resonances in the system that can enhance or suppress the gain and loss of the superlattice. A more detailed understanding can be obtained through the dynamical oscillations in the system, i.e., the ZB effect.

Relativistic Zitterbewegung in photonic superlattice
ZB is a purely relativistic quantum effect resulting from the interference between the positive and negative energy states of the Dirac fermion. To experimentally observe ZB in electronic systems is challenging due to the small wavelength and the extremely high oscillating frequency. Photonic superlattice systems with the underlying equation having the same mathematical form as the Dirac equation provide an alternative paradigm for detecting and characterizing the ZB phenomenon [54]. Our goal is to investigate whether ZB can emerge in pseudo-PT symmetric photonic systems. Figure 4 shows a number of representative time series of the beam center of mass of the wave packet for different parameters. There are apparent oscillations in the time series. Since oscillations are a generic feature of the wave equations, the issue is whether these oscillations are true manifestations of ZB. We address this issue by analyzing the equivalent Dirac equation to determine if the oscillations have a relativistic quantum origin. The basic idea is to write the Dirac equation in the momentum space to obtain the time evolution of the wavefunction [68]. In order to arrive at an analytical form, we make the assumption that the total Hamiltonian at any two given instants of times are commutative so that a sophisticated treatment of the time ordering operator is not necessary. Note that, while for time independent quantum systems time ordering is not necessary, for a time dependent Hamiltonian this may not be the case. Nonetheless, we expect that, in the parameter regime where the non-Hermitian effect is weak, i.e., small σ i , our analytic results would agree with the direct simulations results.
The main results of our analysis can be summarized by dividing the time evolution of the expectation value of the position operator into three components: the drift ξ d (t), the ZB trembling component ξ ZB (t), and the purely imaginary component ξ Im (t). We have (detailed derivation can be found in Appendix) where is the normalization constant associated with the Euler formula for the Pauli matrices, and is the angular spectrum of the Gaussian beam envelop. The modulus of the wavefunction is We note that, for a relativistic quantum system described by the Dirac equation, the time evolution of the expectation value of the position operator contains the drift and the ZB trembling components only. For a stationary Hermitian system, we have A = A * ∝ t, the drift component can be simplified as which describes the motion of the wave packet at constant velocity. The ZB trembling component ξ ZB (t) can be simplified as an oscillatory term described by sin(A) cos(A) ∝ sin(2A), and the purely imaginary component ξ Im (t) is simply zero. The analytical results can be used to interpret the numerically observed dynamical behaviors of the system associated with the boundary variations in the phase diagrams in Fig. 2. Specifically, Fig. 3 shows the analytical results for r = 0.5 and ω/ω 0 = 4.2, 1.5 (the blue solid and red dashed curves, respectively), which correspond to the two rightmost white lines marked in Fig. 2(a). Figure 3(a) shows the time evolution of the normalized modulus of the two sublattices A and B. As expected, the modulus of the red dashed curve is about three times slower than that of the blue solid curve, which is the ratio between the modulation frequencies. However, the time evolution of the ZB trembling term ξ ZB (t) has the same frequency for the two sublattices, as shown in Fig. 3(b), indicating that ZB trembling has little dependence on the modulation frequency. A comparison of the time series in Figs. 3(a,b) indicates that the frequency of ZB oscillations is close to that of the modulus oscillation for ω/ω 0 = 4.2.
Note that, for r = 0 the band gap is 2σ r = 4.2. These results imply that the variations in the boundary between the pseudo-PT symmetric and the pseudo-PT symmetry breaking regions are due to the resonant interaction between the modulation to the waveguides and ZB trembling. The energy absorption process is enhanced when the ratio between the modulation period and the ZB trembling period is an odd integer. For the high frequency region where the modulation period is smaller than the ZB trembling period, the system exhibits a relatively simple behavior, where the boundary can be described by a linear relation between r and ω/ω 0 . Figure 4 shows the simulation (blue) and analytic (red) results for the ZB effect, where the behavior of the real part of the quantity (ξ d (t) + ξ ZB (t))/|ψ(t)| 2 is illustrated for six different parameter combinations indicated by the red stars in Fig. 2(a). When the imaginary part of the refractive index is relatively small, i.e., when the system is only weakly non-Hermitian, the analytical and simulation results agree well with each other, as shown in Figs. 4(a-c) for r = 0.2. Note that the scale of the x axis in Fig. 4(a) is about three times of that in Fig. 4(b-c), and the periods of the oscillations for these cases are the   same, indicating that the oscillation periods have little dependence on the modulation frequency. Our analytical results provide a base to attribute the oscillations as resulting from the ZB trembling term ξ ZB (t), demonstrating their relativistic quantum origin. We note that, for a Hermitian system, the dynamical evolution of the wave packet can be represented as the superposition of a constant drift behavior and simple sinusoidal oscillations, but the dynamics of ZB oscillations in the relativistic non-Hermitian system are richer. Specifically, from Fig. 4(a), we see that, for low modulation frequencies, ZB oscillations can be enhanced or suppressed in different time intervals. However, for high modulation frequencies [e.g., Fig. 4(c)], ZB oscillations are conventional in the sense that there are no apparent enhancement or reduction effects. An intermediate situation arises between these two cases, e.g., for ω/ω 0 = 1.0. These results indicate that the modulation frequency of the refractive index can affect the amplitude of ZB oscillations.
The role of the modulation frequency in the amplitude evolution of ZB oscillations can be assessed in a more detailed manner through a close examination of the oscillation term in ξ ZB (t), i.e., sin A cos * A, for different values of the modulation frequency, as shown in Figs. 3(b-c), where the oscillating patterns have the same modulation frequency, implying that the frequency of ZB oscillations is entirely determined by its sine and cosine components. Further, the ZB oscillations for the two different values of the modulation frequency are in pace with each other in time, indicating that the modulation frequency has little effect on the frequency of ZB oscillations. This is expected because the quantity A(t) is an integration of σ(t) so that the amplitude of σ(t) is a key factor. Nonetheless, the modulation frequency does affect the oscillating intensity, as shown in Fig. 3(a). Thus, the envelope behavior of ZB oscillations in Figs. 4(a-c) is determined by the modulation of the periodic refractive index: for high and low modulation frequencies, ZB oscillations are suppressed and enhanced, respectively. A conclusion is that, for high modulation frequency and low modulation intensity, the effect of refractive-index modulation is weak, as exemplified in Fig. 4(c). For r = 0.5, the analytically predicted amplitude of the oscillation behavior [Figs. 4(df)] deviates from the simulation results (except for the initial phase of the dynamical evolution). In spite of the disagreement, the predicted phase behavior of the oscillations agrees with the simulation results. The failure of the analytic theory to predict correctly the amplitude behavior of ZB oscillations stems from the hypothesis used in our analysis: the Hamiltonians at different times are commutative. This hypothesis is violated when the non-Hermitian effect becomes pronounced, i.e., as the imaginary component of the refractive index and the time dependent modulation are relatively large. Indeed, as we increase the value of r from 0.2 to 0.5, the agreement becomes increasingly worse, especially for large time. However, for the parameter region on the right side of the phase diagrams in Fig. 2 where the imaginary part of the refractive index is small, a reasonably good agreement between the analytical and simulation results is obtained, due to the relatively weak non-Hermitian effect. For the high modulation frequency region (the rightmost part of Fig. 2), ZB oscillations are similar to those of the conventional case where the modulation in the refractive index is absent. Physically, in the high modulation frequency region, the waveguides do not have time to absorb and dissipate energy, generating a mean-field like effect.
In the intermediate frequency regime, the analytical and simulation results do not agree with each other. Numerically, we find the phenomenon of spatial energy localization, as shown in Fig. 5. In particular, in the parameter space near the boundary of the phase diagram, e.g., for r = 0.5 and ω/ω 0 = 4.0, spatially the wave energy is localized within a small region about the center of the waveguide array, as shown in Fig. 5(a). This phenomenon of energy localization is sensitive to the value of ω. Figure 5(b) shows a contrary example for ω/ω 0 = 3.5, where the wave packet spreads Figure 6. Phase diagrams for different values of σ r /κ. (a-e) Phase diagrams for σ r /κ = 1.1, 1.5, 2.1, 2.5, 3.1, respectively, where ω 1 is the base frequency, i.e., the largest resonant frequency. The white dashed lines on (c) correspond to the cases of ω 1 /ω = 1, 3, 5, 7, 9, respectively. There is a general relation: ω base ∝ σ r /κ. after a short time. We find that the wave spread becomes faster as the value of ω is reduced. Figures 5(c,d) show the total wave intensity corresponding to the cases in Figs. 5(a,b), respectively. We observe an apparent envelope modulation behavior in the case of energy localization. The strong reduction in the spread of the wave packet is similar to that reported in Ref. [69], where a localization behavior in quantum diffraction was found to result from the linearization of the quasienergy spectrum close to the PTsymmetry breaking boundary. Generally, this is a dynamical localization phenomenon in one-dimensional driving systems [70][71][72][73]. The width of the localized wave packet depends on the width of the initial wave packet.
Phase diagram in 1/ω axis. As can be seen from the phase diagrams in Fig. 2, the boundary between the pseudo-PT breaking and the pseudo-PT phases exhibits an oscillation pattern. A detailed investigation indicates that the inverse of the peak frequencies are odd multiples of the inverse of the rightmost peak frequency. To better understand this behavior, especially in the limit ω → 0, we present a series of phase diagrams in the (1/ω, r) plane for different values of σ r /κ, as shown in Figs. 6(a-e), where the parameter values in Fig. 6(c) are the same as those in Fig. 2(a). For convenience, in Fig. 6(c) the 1/ω axis is rescaled as 1/(ω 1 /ω 0 ), where ω 1 /ω 0 = 4.2. Other panels are scaled using the base frequency of Fig. 6(c). In Fig. 6(c), the white dashed lines are for ω/ω 1 = 1, 3, 5, 7, 9, which is indicative that, in the vicinity of these lines, the system is in the pseudo-PT breaking phase. As noted, we have ω 1 /ω 0 ∼ 2σ r /κ, so as σ r /κ is varied, we expect to see a similar resonance phenomena but with a different base frequency, as shown in Figs. 6 (a,b,d,e). Physically, the quantity σ r /κ measures the band gap of the massive relativistic dispersion relation, and it is known for the Hermitian case that the ZB trembling frequency is proportional to the band gap. As a result, the modulation frequency, which is at resonance with the ZB trembling frequency, is also modified. Similar results were obtained recently [55][56][57][58][59].

Conclusion
Novel electronic materials obeying relativistic quantum mechanics are a current focus of condensed matter physics and materials science, which provide us with a platform to uncover, understand, and exploit unusual physical phenomena. It is of great interest to use optical systems to gain insights into the fundamentals of relativistic quantum solid state devices, and vice versa, i.e., to exploit relativistic quantum electronic behaviors to generate revolutionary ideas/methodologies in manipulating light propagation. The main reason for such an interest is that, electronic materials face many experimental and technological challenges due to the extremely small wavelength issue. The conventional method is to engineer, e.g., through strain, doping and voltage control, some specific materials to realize the desired configuration, constraining researchers to some limited kinds of lattice structures and interactions. Compounded on the small wavelength issue are various kinds of disorders and impurities. Photonic crystals and waveguides can be used to mimic the electronic properties, e.g., the dispersion relations, of solid state materials, with the advantage that the optical wavelength is orders of magnitude larger than the electronic wavelength. Substantially larger optical devices can then be fabricated to investigate a host of electronic behaviors in relativistic quantum solid devices. These are synthetic photonic materials and devices.
This paper is a case study of exploiting the equivalence between coupled wave equations and the Dirac equation to model/simulate relativistic quantum effects in a non-Hermitian setting. Standard quantum mechanics requires observable operators to be Hermitian to ensure real eigenvalues. However, non-Hermitian Hamiltonians are also capable of generating a purely real eigenvalue spectrum insofar as the system has a PT symmetry. In optics, non-Hermitian Hamiltonian systems can be realized using materials with a complex index of refraction, rendering experimentally feasible such systems [2,4,23,24]. We studied a class of non-Hermitian waveguide systems with periodic modulated imaginary refraction index in order to uncover the optical counterpart of the relativistic ZB effect, which exhibit PT symmetry breaking at any instant of time but the symmetry is preserved on a larger time scale (herewith the term pseudo PT symmetry). We generated phase diagrams of the pseudo-PT behavior, which provide a general picture for controlling and harnessing light propagation in the waveguide system. The phenomenon of oscillatory boundary variations in the phase diagrams is explained. In the low modulation frequency region, we observed periodically enhanced and suppressed ZB oscillations, as predicted by an analysis of the equivalent Dirac equation. The phenomenon that there are time periods where the ZB trembling is absent, interspersed by time intervals where the relativistic quantum oscillations emerge and disappear, a kind of "revival" phenomenon, can potentially lead to a new mechanism to manipulate/control light propagation. In the high frequency regime, a conventional ZB trembling behavior arises due to a mean-field effect. In the intermediate region, the equivalent description based on the Dirac equation no longer holds. However, numerically we uncovered a spatial energy localization phenomenon. Derivation of the main analytical result Eq. (6). We start from the Dirac equation in the momentum space: The wavefunction at any given time can be written as where T is the time ordering operator. Since our system is a time dependent non-Hermitian system, i.e., [H(t 1 ), H(t 2 )] = 0, the expansion of the exponential factor under the time ordering operator can be quite sophisticated. In order to gain analytical insights, we assume [H(t 1 ), H(t 2 )] = 0. This does not mean that the non-Hermitian property of our system is lost. In fact, non-Hermitian characteristics still exist in E ± (t).
We first discuss the validity of the approximations. The effective Hamiltonian of the system is H(k, t) = κkα 1 + σ(t)α 3 , where time modulation occurs in the imaginary part of σ: σ(t) = σ r + iσ i sin(ωt), while the real part σ r is time independent. In the limit σ i = 0, system becomes time independent, rendering irrelevant the time ordering operator, as lim σ i →0 [H(t 1 ), H(t 2 )] ∝ κkσ i → 0. The time dependent weight on the α 1 term comes from the integral t 0 H(t ′ )dt ′ , which is responsible for the dynamical phase in the limit.
For a small value of σ i , the time ordering operator can be written as where, for simplicity, we have ignored . Take the second order term ] as an example. The difference induced by the time ordering operator is on the order of [H(t 1 ), H(t 2 )] ∝ κkσ i , while the main part of the second order term is dominated by κ 2 k 2 + σ 2 r . Since we focus on the excitations about the k ∼ 0 points, i.e., the Brillouin zone boundary in the original waveguide system, the ratio of the two terms, κkσ i /(κ 2 k 2 + σ 2 r ), is negligibly small, providing a justification that the effect of the time ordering operator can be neglected in the regime σ i /σ r → 0. For higher order terms, the effect of the time ordering operator will be even more negligibly small.
Taken together, for σ r = 0, our results are exact for σ i = 0, as the system reduces to a time independent system. For σ i << σ r (e.g., σ i /σ r = 0.2), the effects of the time ordering operator can be neglected. From the comparison with the numerical results, we can see the approximation works well when σ i /σ r = 0.2. This agreement demonstrates that the oscillations we observed are indeed Zitterbewegung, which possesses a modulated amplitude resulting from the non-Hermitian modulation of the system. When σ i /σ r increases to 0.5, our analytical results cannot quantitatively reproduce the numerical results. However, the positions which are of tremendous oscillations are well predicted.
Since σ(t) = σ r + iσ i sin(ωt), we have Next we consider the initial condition. At the input plane z = 0, the initial wave packet can be chosen to have the shape of a slowly varying Gaussian form. The initial amplitude of ψ 1,2 (ξ, 0) is then proportional to G(2na) and G((2n − 1)a) ≈ G(2na), so the two components of the Dirac spinor have the same initial condition, i.e., ψ k (0) = G(k) [1,1] T . The Pauli matrices operate on the initial Dirac spinor, and we get the time evolution of the two components of the Dirac spinor as The expectation value of the position operator can be calculated through A lengthy algebra leads to the main analytic result: the formulas in Eq. (6) as well as the time evolution of the modulus of the wavefunction. The convergence of the integrals in Eq. (6) is guaranteed by the exponential decay of the Gaussian wave spectrum G(k). In particular, for small t values, | cos(A)| 2 is about unity. For large values of t, we have |ψ(t)| 2 ≈ | cos(A)| 2 + | sin(A)| 2 1. Thus, |ψ(t)| 2 can never approach 0, so Eq. (6) converges.
The real part of the dispersion relation in the limit of small imaginary part is then given by ± σ 2 r − σ 2 i (z) + 4κ 2 ] cos 2 (qa). (.7)