Introduction

The inverse Doppler effect where an observer measures a higher frequency from a source moving away, is such an intriguing phenomenon that it has attracted substantial theoretical and experimental attention since it was first proposed1,2,3,4,5. Our team observed the inverse Doppler effect at optical frequencies for the first time in 20116 using left-hand photonic crystal (LH-PhC)7. However, in that article, the equation

where c is the vacuum velocity of light, f0 is the original frequency of the light source (for instance, CO2 laser in ref. 6) and f1 is the Doppler frequency at the second interface of the PhC, was used to calculate the frequency shift and was based on phenomenological theory. In that research, we considered the contribution of the LH-PhC by treating it as a homogenous medium with a negative np refractive index and used the classic Doppler formula of a mechanical wave8,9,10, which makes no sense with np < 0. Hence, we were inspired to elucidate the physical mechanism underlying the inverse Doppler Effect in PhC and determine the rationale for using equation (1) to obtain the inverse frequency shift6.

We first address the controversy surrounding the explanation of the inverse Doppler effect in magnetic, nonlinear transmission lines2,11. N. Seddon and T. Bearpark attribute it to vg · vp < 0 (vg is the group velocity and vp is the phase velocity) in the second Brillouin zone (BZ), whereas Evan J. Reed considered a time-dependent reflection coefficient in the first BZ, with vg · vp > 012,13. By contrast, Alexander Kozyrev demonstrated that the shift arises from the complex spatial structure of propagating waves, some of which have vg · vp < 0 resulting in the inverse Doppler Effect14. However, for light traveling through a LH-PhC, controversy exists regarding whether vp is opposite to the wave vector k along the propagation direction15,16. Based on observations of Bloch harmonics and negative phase velocities in PhC waveguides resulting from the typical dispersion17, we believe that some spatial harmonics in LH-PhC have vp·k < 0 in terms of the decomposition of Bloch waves, which leads to the inverse Doppler effect in moving PhC.

In this article, we investigate the main spatial harmonics components of the complex field distribution in a LH-PhC and describe out the backward-propagating wave. Subsequently, by analyzing the phase evolution of the entire optical path using the static finite different time domain (FDTD) method, we obtain a Doppler shift that is in accordance with the experimental results. Finally, an improved dynamic FDTD method is proposed to reproduce the continuous motion of the PhC in the experiment6 to the greatest extent and the rationality of np, which was used in Eq. (1)6 is verified.

Results

Analysis of experimental setup and corresponding model

To start our computation, first we analyze the experimental setup used in ref. 6. As it is shown in Fig. 1(a), a transverse magnetic (TM) polarized laser beam at the desired wavelength (λ = 10.6 μm) is divided into two beams by the first Ge splitter. The transmitted beam travels through the LH-PhC which is fixed on a moving platform and the reflected beam is then reflected by a Cu reflector and the second Ge splitters, one of which is fixed on the moving platform. Finally, the two beams interfere on the detecting surface. The LH-PhC under investigation, which was achieved by electron beam Lithography, is a triangular silicon rod array prism with a vertex angle of rhombus geometry of 60°, its effective refractive index (neff) was measured to be −0.50626. The radius of the silicon rods was estimated to be 0.226a (where a is the lattice constant, which was designed to be as 5 μm) from the equi-frequency surface (EFS), which is close to the design requirement6.

Figure 1
figure 1

Schematic representation of the experimental setup used in ref.6.

The blue arrow points in the direction of the moving platform, which holds the PhC and splitter.

Full size image

Phase evolution of light traveling through PhC

To analyze the phase evolution in the LH-PhC as light passes through it, we employed the FDTD method to simulate the electrical field distribution and set the spacing grid to a/40 to ensure accuracy. We recorded the 2-D electrical field data (the size of the recorded data is 2000*7000) of light traveling in the LH-PhC and showed the corresponding spatial field distribution in Fig. 2(a). The entire electrical field distribution of light traveling through the LH-PhC and air is manifested in the inset of Fig. 3(b). In Fig. 2(a), the field distribution was periodic in both the PhC and air along the propagation direction. The periods were designated as λ and λ0 and the horizontal coordinate is the position along the propagation direction. We averaged the recorded data along the ΓK direction (perpendicular to the propagation direction) and obtained a 1*7000 matrix: this intensity-averaged field is shown in Fig. 2(b) as a blue-labeled line with periodic peaks that can be clearly observed.

Figure 2
figure 2

Component analysis of the field when light traveling though the LH-PhC.

(a) Field distributions along the propagation direction in air and in the PhC. (b) The Original averaged field is labeled with a blue line and the red line represents the spatial waveform of the electrical field component E1 extracted from the original field. (c). The blue line is a close look into the blue line in the dash rectangle of (b) and the red line indicates the extracted spatial waveform of the field component E2.

Full size image
Figure 3
figure 3

Phase changes for the main components along the propagation direction.

(a) The spatial Fourier spectrum of the original field data reveals several spatial frequencies: F0, F2 and F3. The insets present the electrical field distributions of the incident wave, E1 and E2, at adjacent times, corresponding to F0, F1 and F2, the black lines shos field distribution at one moment, the blue lines shows the field distribution for the previous time and the red lines shows the field distribution for the later time. (b,c) Spatial field distributions recorded along the light path as light passes through the LH-PhC and RH-PhC, The red arrows indicate the interfaces between the PhC and air. In the LH-PhC, the phase evolutions on either side of the interfaces are opposite. In the RH-PhC, the phase evolutions on either side are the same and the insets show light propagating through the LH- PhC and RH- PhC. The black dotted lines denote the interface between the PhC and air.

Full size image

The electrical field in a PhC can be decomposed into a series of plane waves with different amplitudes and frequencies because the field obeys the Bloch theorem. The Fast Fourier Transform (FFT) method is highly advantageous for analyzing the complex frequencies in a 2-D PhC16,17,18,19. Here, a 1-D FFT was applied to the averaged electrical field data and the spectra are depicted in Fig. 3(a). We denote F0, F1 and F2 for the three peaks and list their corresponding spatial frequencies in the form of wave vectors in Fig. 3(a). As expected, one of the spatial frequencies, 0.472(2π/a), is exactly the same as the actual spatial frequency of a CO2 laser in air as calculated from F0 = a/λ0 (λ0 = 10.6 μm). Other two peaks labeled with F1 and F2 are the spectra in the PhC. These two main components have large amplitudes, while there are also few tiny noises with small amplitude, which are arising from the numerical FFT. Differing with that in the homogenous left hand materials (LHMs), the beam propagation in the PhC takes the form of a Bloch wave, which can be decomposed into a series of complex plane wave; however, it is difficult to physically assign global phase-front velocities and single wave vector to the global Bloch wave. To overcome this issue, B. Lombardet proposed an appropriate way, which is to use the component of which the wave vector dominates the Fourier decomposition of Bloch wave to analyze the Bloch wave19. In Fig. 3(a), we can see that F1 and F2 possess most of the components’ energy and thus, we focus on these two components in the following analysis.

For further study, we extracted the two main electrical field components E1 and E2 from the field in the PhC. More specifically, we filtered F1 and F2 from the spectrum and then applied an inverse fast Fourier Transform (iFFT). The corresponding reconstructed fields were shown in Fig. 2(b,c) with red lines. E2 is a close look into the dashed rectangle of E1 and E1 is in a manner similar to an envelope of E2 by modulating its amplitude.

Hence, we can assign two different refractive indexes for these two components when they propagate in the LH-PhC because of the unusual dispersion character of the LH-PhC and the refractive index can be calculated as |neff| = λ0/λi = Fi/F0 (i = 1, 2), with two results of 0.5016 for E1 and 1.9406 for E2. Notably, the value 0.5016 is close to the absolute value of the measured value of the negative refraction index, 0.5062. If F1 corresponds to negative refraction, the sign for neff should be negative, which implied the component E1 propagates backward. E2 should be propagating forward. To confirm our speculation, we analyzed the phase evolution Δφ20 of the electrical field by recording the electrical field at adjacent time slots with an interval ΔT of 4.067fs, which is smaller than the light period, T = 35.3fs, to ensure the validity of the inferred propagation direction. The field distributions of the incident wave, components E1 and E2 are illustrated in the insets of Fig. 3(a), beside three peaks. Here, black lines represent the field distribution at one moment, blue lines represent the field distribution at the previous time and the red lines represent the field distribution at the next time. As expected, for the incident wave and E2, the red line always leads the black line, indicating a phase lead, Δφ > 0. By contrast, for E1, Δφ is negative, suggesting a phase lag21. Considering the phase velocity vp = Δφ/(k · ΔT), for the incident wave and E2, vp > 0, which means that the wave moves forward, neff = 1.9406 for E2. However, as for E1, vp < 0, then vp · k < 0, i.e., the phase velocity and the wave vector are antiparallel. Therefore, the phase of E1 is delayed and neff is −0.5016 instead of 0.5016. E1 and E2 correspond to plane waves propagating with opposite phase velocities but they are coupled because of the periodicity and travel forward at the same group velocity.

As for E1, the corresponding wave vector is 0.237(2π/a), which is in the first BZ, because |neff| = 0.5016 is smaller than that of air (nair = 1), E1 can exit the interface of the PhC. Meanwhile, each point on the equi-phase surface of E1 undergoes different phase reduction when it exits obliquely from the interface of the LH-PhC to the normal medium. Thereby, under the phase conservation law, these points undergo different optical path and form a new equi-phase surface which is on the same side of normal as the incident beam, yielding the negative refraction phenomenon. However, for component E2, the corresponding wave vector is 0.915(2π/a), which is out of the first BZ, neff = 1.9406 is larger than nair and because the critical angle (31°) is smaller than the incident angle (60°), total reflection occurs, thus it cannot couple to free-space, as illustrated by the black arrow in the inset of Fig. 3(b).

Doppler shift obtained using the static FDTD method

The spatial distributions of the field inside and outside of the PhC were independently discussed above. Now let us consider the entire process of the transmission of light through a stationary PhC after the field stabilizes. Since E1 is the only component that can couple into free-space from the LH-PhC in an incident angle of 60° at the second interface of the LH-PhC, we extracted E1 from the recorded data of the PhC’s area at one time and treated it as the only contributing field in the PhC because vp < 0. Then, we expanded the spatial field distributions in front of, inside and behind the PhC in the same optical path. Using the same process, the spatial field distributions at the previous and later times were all analyzed, as depicted in Fig. 3(b). The phase velocity directions are opposite on either sides of the interface between the air and PhC, which echoes the analysis above. For comparison, we designed an RH-PhC with a neff of 1.12 at the incident wavelength and a = 5 μm. At an incident wavelength of 31.25 μm, the corresponding propagation schedule inserted in Fig. 3(c) exhibited a normal refraction behavior. By performing the same analysis, the spatial field distribution was obtained, as shown in Fig. 3(c). Under these conditions, the phase change remained positive in both the air and the PhC, because the red line always appears in front of the blue line. Therefore, the phase evolution of E1 was indeed abnormal because of its inverse phase velocity, which caused the phase on the detecting surface to vary in an unusual way.

In the experiment6, in which the platform was moving, the optical path changed continuously in both branches of the optical setup, which caused a phase variation on the detecting surface. Those two branches can be simplified as a signal beam with light traveling throughout the LH-PhC and a reference beam with light reflected by the moving splitter. As displayed in Fig. 4(a), at time t1, the phase of the light source is φl1 and the phases of the signal and reference beams on the detecting surface are φs1 and φr1, respectively. After the platform moves distance Δx at velocity v, at time t2, the phases became φl2, φs2and φr2, respectively. During the interval Δt = t2 − t1, the phase increments for each item are φs = φs2 − φs1 and φr = φr2 − φr1. From this perspective, the frequency shift can be expressed as Δf = |φs − φr|/(2π · Δt).

Figure 4
figure 4

Analysis of the inverse Doppler effect based on the static FDTD method.

(a) Simplified graph of the two interference branches when the platform moves. (b) The spatial field distribution along the signal beam path at different times. The terminal point denoted by the red circle is the detection point. (c) Comparison between Δf, Δfe and Δfc. The insets present the frequency spectra of the experimentally measured beat waves.

Full size image

Because the software we used (FDTD solutions) weakly treats moving objects, we have to divide the continuous movement into a series of discrete motions at different times. We assumed that v = 0.0244 mm/s, as was one of the speed used in the experiment. As the platform moved away from the detector step by step over Δx (Δx = 1 μm), the initial phases φli for each time ti were adjusted to their appropriate values in the software using the equation (λ/c)/(i · Δx/v) = 2π/φli. Therefore, by employing the analysis mentioned above, the spatial electrical field distributions of the signal beam for three different moments along optical path are shown in Fig. 4(b), where the red circle indicates the detection point. The phase differences Δφs between each adjacent time can be obtained by subtracting the former value of the detection point from the latter. For reference beam, the change in light path was Δl = v · Δt · sin(ϕ/2 + θ1 + θ2)/sin(ϕ/2), which can be inferred from the geometrical relationship shown in Fig. 3(a) by using the sine theorem. Using the same subtraction operation, the phase difference Δφr = 2π · Δl/λ could also be obtained. Consequently, the frequency difference Δf of the two arms was calculated to be 1.949 Hz. Other frequency difference data corresponding to the velocities used in the experiment, namely, 0.0122 mm/s, 0.0244 mm/s and 0.0732 mm/s, were also be obtained, as shown in Fig. 4(c) and the measured frequency shifts Δfe and theoretical frequency shifts Δfc calculated from the expression are also presented.

The calculated Δf data were unambiguously consistent with the measured Δfe data and the theoretical Δfc data, with a certain error between Δf and Δfc that does not affect the judgment of the inverse Doppler effect. It illustrates that the above analysis of the physical mechanism is correct. Part of the discrepancy can be clarified by the spectra of measured beat waves, which is depicted in the insets of Fig. 4(c). Because they are influenced by the measurement noise, the spectra by no means resemble idealized Dirac functions but are expanded in a small range.

Doppler shift obtained by the dynamic FDTD method

The above analysis was based on the variation of phase change obtained by slicing the continuous motion into a suite of static movements at different time and the relationship between each movement only serves to add an increased source phase source to the former one. This analysis cannot sufficiently reproduce such a continuous process. In fact, unlike a mechanical wave, the propagation of electromagnetic waves is influenced by the electromagnetic fields at previous times when relative movement occurs because they obey Maxwell’s equations. Specifically, the existing electrical field and magnetic field data from the former time will be added to the new field data at the later time to create a new field distribution and another field distribution will be produced for the next time using Maxwell’s equations. To this end, we developed an improved dynamic FDTD method. However, this method is limited in handling low-speed motion because the movement at each step time of the FDTD is very small. To ensure appropriate resolution, the grid should be set to a nanometer scale, which would require thousands of hours of simulation time with a general-purpose computer. In this study, after evaluating the time consumption, resolution and analysis requirements, we set the grid size of the FDTD to 0.05 μm and simulated the inverse Doppler Effect with velocities of 0.01 c, 0.02 c and 0.05 c, in that order. In the simulation, a specified point along the outgoing beam was treated as the detector, at which the electrical field data were recorded at each time. The electrical field distribution for a velocity of 0.05c is depicted in Fig. 5(a), the red dot denotes the position of the detector, the red arrow points in the direction of emitted light and the blue arrow points in the direction in which the PhC moves. Taking the FFT of the record field data from the detector, the detected frequency fd was found to be 2.854 × 1013 Hz. Both the record field and spectrum are shown as insets in Fig. 5(a). Because this frequency shift could not obtain experimentally since such high speeds were not attainable in the experiments, for comparison, a calculated frequency fc of 2.855 × 1013 Hz was obtained by using the expression in6, where f0 is the original frequency of the CO2 laser source, 2.827 × 1013 Hz and the effective refractive index np is −0.5062, which was the same as that in the experiment6. Using the same approach, the detected frequencies and calculated frequencies for the remaining two speeds were also obtained, as shown in Fig. 5(b).

Figure 5
figure 5

Analysis of the inverse Doppler effect based on the dynamic FDTD method.

(a) Simulation of light transmitted through the moving PhC at a velocity of 0.05 c. The red dot denotes the position of the detector, the red arrow points in the direction of the emitted light and the blue arrow points in the direction of which the PhC moves. The inset shows the recorded field obtained at the detection point and its Fourier spectrum. (b) Comparison between the detected frequencies fd and the calculated frequencies fc at three speeds: 0.01 c, 0.02 c and 0.05 c.

Full size image

The detected values of fd obtained using the improved dynamic FDTD method at three velocities namely, 0.01 c, 0.02 c and 0.05 c were very close to the calculated values of fc, with a relative error that did not exceeded 5%, which is less than that of the static FDTD method. This finding supports our assumption that the electromagnetic field of the former time affected the later field and influenced the frequency shift when there was relative movement, possibly contributing to the discrepancy mentioned in the discussion of the static FDTD method.

Discussion

We analyzed the Bloch wave in the PhC by treating it as a summation of a series of spatial harmonics and extracted the main spatial harmonic components. The comparison between the phase changes of the electrical fields in the PhC and those in air identified an unusual component in the PhC, with vp·k < 0, producing a negative refraction phenomenon. Subsequently, in terms of the phase evolution, by using the electrical field data obtained from the static FDTD method, we calculated the frequency difference between the signal and the reference beams at detection point. The result demonstrates that the experimental data were identical to the calculated values obtained with the static FDTD method, which supposes that the component E1 that yields to negative refraction leads to the inverse Doppler effect. To address the insufficiency of the static FDTD method for handling moving objects and to further verify our analysis, we improved on the traditional FDTD method and simulated the inverse Doppler effect with the same PhC parameters used in the experiment. Furthermore, the simulated result generated by the dynamic FDTD method was consistent with the calculated value by using the expression in ref. 6, which strongly supports the rationality of considering the contribution of the LH-PhC by simply using a negative np. Hence, the phenomenological result of our experiment can be explained by a clearly physical mechanism. Additionally, our study provides a new way to analyze the behavior of EM waves travelling in PhC, particularly, in moving objects.

Methods

EFS

The EFS varies with the parameters of the PhC, such as the radius of the rod (hole), refractive index and shape. From the EFS, we can derive the expression np = k/(2π · ω). Thus, the relationship between the radius of the rod and the refractive index is known and we can use it to infer the radius of the rod from the measured refractive index.

FFT and filter

We extracted the main components of the complex spatial field distribution in the PhC. We then applied the FFT to the original data and identified the desired peak in the spectrum. Specifically, we retained the amplitude of the required frequency and set other amplitudes to zero and then used the iFFT to reconstruct the required component without losing the phase information.

Dynamic FDTD method

In contrast to the static FDTD method in which objects remain stationary, for moving objects, the dynamic FDTD method considers the fact that the existing electrical fields and magnetic fields at earlier times should be added to the new field data at later time to create a new field distribution; this new field can then be used to produce another field distribution at the next time according to Maxwell’s equations. In the dynamic FDTD method, at each sliced time, the position of PhC was changed by varying its parameters (namely, the dielectric constant and magnetic permeability) in each grid. The electrical field and magnetic field were stored as the initial field for the next calculation time.

Additional Information

How to cite this article: Jiang, Q. et al. Mechanism Analysis of the Inverse Doppler Effect in Two-Dimensional Photonic Crystal based on Phase Evolution. Sci. Rep. 6, 24790; doi: 10.1038/srep24790 (2016).