Dual Doppler Effect in Wedge-Type Photonic Crystals

The dual Doppler effect, in the simultaneous occurrence of both normal and inverse Doppler effect in one moving two-dimensional wedge-type photonic crystal, is proposed. An improved finite-different time-domain algorithm is used to verify this phenomenon. The spatial Fourier Transformation has been applied to complex electrical field data to reveal the mechanism. The harmonics with negative spatial frequencies show a lagging phase evolution, while those with positive spatial frequencies show a front phase evolution. Different wedge-type photonic crystals are designed to filter out the required harmonics based on the systematic study of spatial Fourier Transformation and wave vector diagram. Our work paves a new way for Doppler cooling of atomic gases, radar deception, invisibility cloaks, microstructure dual frequency interferometer and so on.

The Doppler effect can be observed in every walk of life. Applications of the effect has been widely established, including radar, laser vibrometry, blood flow measurement, and the search for new astronomical objects. It can be described as the effect produced by a moving source of waves in which there is an apparent upward shift in the frequency for observers towards which the source is approaching. However, the inverse Doppler effect performs in an opposite way. In 1968, G. Veselago predicted that when electromagnetic waves propagated through materials with negative permittivity ε and negative permeability μ, the frequency should drop as a source moving towards an observer but increase as it moving away 1 . The inverse Doppler effect was first observed at radio frequency by N. Seddon with an experimental magnetic, nonlinear transmission line in 2003 2 . After that, this anomalous effect was observed in more than one kind of special dispersive media, like left-handed material (LHM) whose effective permittivity ε and permittivity μ are both negative, and photonic crystals (PhC) using shock discontinuities 3 .
LHM possesses the feature of negative refractive index (NRI) 4 . When a light travels through it, the phase velocity v p is antiparallel to the Poynting vector  S . Interestingly, the NRI feature can also be observed in some PhCs whose equi-frequency contour (EFC) is shrinking towards center. Differing from the homogenous LHM, the NRI for the PhC is caused by complex dispersion of the periodic structure. The unusual dispersion leads to many extraordinary phenomena, such as self-collimation 5 , super-prism 6 , zero phase delay 7 , anomalous and negative refraction 8,9 . With one incident wave, the refracted wave in homogenous LHM is unique and acts as a backward wave, of which the energy flows forward while the equiphase front moves backward. Nevertheless, refracted light in the PhC will be no more unique. S. Foteinopoulou had made a detailed illustration how to determine the directions of all refracted lights in the PhC by using wave vector diagram 10 . Propagating waves in the PhC consist of both forward waves and backward waves 11,12 , which has been demonstrated experimentally by analyzing the phase evolutions for the retrieved harmonics of the measured Bloch wave 13 .
The inverse Doppler shift at optical frequency was observed in the PhC in our previous research in 2011 4 . Considering the explanations of N. Seddon's experiment 2,14,15 , we have revealed that the inverse Doppler shift in the PhC is caused by some spatial harmonics in the PhC which have ⋅ < v k 0 g (the backward wave) in terms of the decomposition of Bloch waves, where v g is the group velocity. As to the rest spatial harmonics which have v k 0 g ⋅ > (the forward wave), they cannot couple to the vacuum due to total inner reflection (TIR) 11 . Hence, it is significant to investigate the Doppler effect phenomenon in the PhC if both backward waves and forward waves couple to the vacuum by appropriate design of the wedge type of the photonic crystal.
In this work, by using the improved finite-different time-domain (FDTD) method, we utilize a moving wedge-type PhC to realize a novel phenomenon, dual Doppler effect, in which both normal Doppler shift and inverse Doppler shift occur simultaneously. For the purpose of illustrating the mechanism, we use spatial Fourier Transformation (SFT) on the complex electrical field data of Bloch wave in the propagation direction. Based ScIenTIfIc REPORts | (2018) 8:6527 | DOI:10.1038/s41598-018-24941-8 on our developed algorithm, the phase evolutions of the spatial harmonics could be illustrated in a comprehensible way by the retrieved field and the wave vector diagram. The multiple outgoing waves from the slope of wedge-shape PhC could be determined in use of multiple reciprocal unit cells. It is firstly observed that the dual Doppler shift is governed by the forward wave and backward wave in the PhC, and the shifts can be calculated with the effective indexes of the PhC corresponding to each harmonic.

Results
Verification of the dual Doppler effect. The structure of the PhC under investigation is shown in Fig. 1, silicon rods with the radius of 0.226a (a is the lattice period) and ε of 11.96 form themselves as a wedge type 4,11 . This PhC, with the EFC showed in Fig. 1(c) portraying a shrinking contour at the normalized frequency of 0.4717(2π/a), has been verified experimentally to be a negative refraction material with the NRI measured as −0.5 4 . In that experiment, light travels along ΓM direction and the only light exits from the slope of which the vertex angle is 60°. Here, by changing the light propagation direction and the vertical angle, refracted light exiting from the slope can be divided into two beams. Figure 1 shows two cases, as for case 1 in Fig. 1(a), light traveling along ΓM direction (y direction) with the vertex angle of 30°, the silicon rods array with the period of a 3 forms the slope. While for case 2 in Fig. 1(b), light traveling along MK direction (x direction) with the vertex angle of 30°, the slope is formed by silicon rods array with the period of a.
Refracted light beams in these two cases appear both positive and negative refractions. Since the negative refraction beam leads to inverse Doppler effect 4 , we speculate that with both positive and negative refractions happen simultaneously, a dual Doppler phenomenon will occur. We use the dynamic FDTD method mentioned above to verify that. Figure 2 shows the corresponding Sketch map, where the rods in green and purple indicate the PhC before and after moving, the double-sided arrow indicates the moving direction, the point P on the second interface of prism represents the position where the light exits and moves to Q with the motion of PhC. For the sake of saving running time and assuring the accuracy 11 , the moving speeds of the PhCs are set as 0.005c, 0.008c and 0.02c, where c is the speed of light in vacuum. With the PhCs moving at a speed of 0.008c, the corresponding field distributions are detected and exhibited in the insets of Fig. 3.  Three sampling points S 0 , S 1 and S 2 , labeled with red dots in Fig. 3(a,b), are set to record electrical field of light source, negative refracted light and positive refracted light correspondingly. The corresponding time domain frequencies ω s , ω ni , and ω pi (i = 1, 2, corresponding to case 1 and case 2) are then obtained by the fast Fourier Transformation (FFT) method. The spectrum shown in Fig. 3 manifests that the frequencies of positive refracted light and negative reflected light shift in opposite way, with their frequencies locating in both sides of the source frequency. As is shown in Fig. 3(c), for all the velocities, frequencies of negative refracted lights are larger than those of light source, while frequencies of positive refracted lights are smaller than the source frequency. Without loss of generality, the incident light in our model is chosen as CO 2 laser, which possess a wavelength of 10.6 μm. The detecting source frequency ω 0 is 2.827022 × 10 13 Hz, which has a relative error of as low as −0.04% to the theoretical frequency ω s of 2.82823 × 10 13 Hz (calculated by c/λ). This confirms the accuracy of the dynamic FDTD algorithm.

Analysis of the dual Doppler effect.
To analyze the spatial harmonics in the PhC and identify the positive spatial frequencies and negative spatial frequencies, the SFT is applied for the complex electrical field along the propagation direction in the PhC 16 . As for case 1, the spatial spectrum for complex field data along ΓM direction, as shown in Fig. 4(a), has a series of peaks which have a certain periodic property. This can be illustrated from the perspective of Bloch theorem. Waves traveling in the PhC consist of a series of plane waves of which the wave vectors are periodic in respect to the reciprocal lattice vector. For conciseness, only several reciprocal unit cells are shown in Fig. 4(b), where the blue circle and green circles indicate the EFCs of air and PhC in the frequency of 0.4717(2π/a) correspondingly. With the phase match condition 17 , countless plane waves along ΓM direction are excited by the normal incident wave vector (bold blue arrow). Under the guidance of the principle that the Poynting vector should point to the shrinking direction and should be away from the source 10 , the directions of Poynting vectors S m and the wave vectors K m (m is the order, m  ∈ ) of the excited plane waves can be determined, as shown in Fig. 4(b), labeled with black and pink arrows respectively. Wave vectors K m can be written as . This periodicity is obvious in the spectra of Fig. 4(a). Furthermore, the phase evolution for each spatial harmonics is obtained by using FFT and iFFT methods 7 , as is shown in the inset of Fig. 4(a). It shows the phase of harmonics with m ≤ 0 is lagging, which corresponds to the backward waves. While for harmonics with m > 0, the phase is front, which corresponds to forward waves. These are consistent with the analysis in the wave vector diagram in Fig. 4(b).
For each spatial harmonic, the corresponding effective refractive index of the PhC can be assigned as n f f / m s = , where f s = 0.4717/a is the spatial frequency of light source. Considering the TIR condition, with regard to the reflection happened on the slope of the PhC, only the corresponding spatial harmonics with a corresponding |n| smaller than 1/sin(θ 0 ) can couple into the vacuum. Among all these harmonics, only f 0 (n = −0.498) and f 1 (n = 1.946) are satisfactory, with the emergent angles being 0.08π and 0.43π in terms of Snell's law. Actually, the outgoing directions of reflected lights can also be predicted precisely from the wave vector diagram. As is shown in Fig. 4(b), the dotted lines represent the normals of the interface, also known as the conservation lines. The conservation line of K 0 intersects the EFC of air at two points A and B, while another line for K 1 intersects the EFC of air at points C and D. In consideration of the causality, the correct directions of the outgoing waves can be determined as OA and OC, as labeled by red arrows, which agrees with the propagation in Fig. 1(a). As for other plane waves, the corresponding conservation lines have no intersection with EFC of air, which means they cannot couple to the vacuum. Likewise, the only refractive light for our previous work 3 is also clarified in Fig. 4(b). As the conservation line (represented by the solid black line) intersects the EFC of air at points E and F. As expected, only OE is the right emergent direction.
In case 2, the Bloch waves traveling along MK direction possess a period of G x = 4π/a 18 . Wave vector diagram in Fig. 5(b) shows three excited plane waves in one reciprocal unit (as shown in the red circle). Hence, the corresponding wave vectors can be expressed as is the wave vector in first reciprocal unit, and i = 1, 2, 3 represent those three excited waves. Likewise, the spatial frequencies in the spectra can be . However, the components of f m,2 vanish in the spectra. A possible explanation can be given from the perspective of the definition of group velocity, ω = ∇ v g k . That the equi-frequency contour is sharp when the f m,2 components are excited, which makes the derivative for wave vector of the sharp point discontinuous. Thus, it leads a failure to define a certain direction of group velocity. Nonetheless, peaks in Fig. 5(a) can be divided into three groups from the perspective of order number m, with a spatial interval of G /2 x π. Phase evolutions for these harmonics depicted in the inset of Fig. 5(a) show that the harmonics are forward for m > 0 while backward for m ≤ 0. Similarly, the directions of outgoing light beams can also be easily determined from the wave vector diagram. As is shown in Fig. 5(b), two conservation lines intersect the EFC of air at four points. To preserve causality, only OH and OG point to the reasonable direction. The exact angles can be obtained from the spectra with the method mentioned above, as θ 1 of 14.18° for negative refraction and θ 2 of 54.48° for positive refraction.
Waves traveling in the PhC consist of both forward waves and backward waves in these two cases. Taking one of them to analyze the Doppler effect, as showed in Fig. 2, the PhC prism moves with a velocity of v while light traveling through it. Therefore, the Doppler frequency shift ω′ at the second interface can be calculated simply as 4 where n 1 and n 2 are effective refractive indexes of the PhC corresponding to those two spatial harmonics which show negative refraction and positive refraction. The calculated dual Doppler frequencies (ω cn , ω cp ) and the detected values (ω n , ω p ) for those three velocities are showed in Table 1. The Doppler frequencies deduced from equation (2) and equation (3) are verified precisely by the dynamic FDTD method with a relative error of less than 0.2% between detecting values and theoretical values. With ω ω > n s and ω ω < p s , the dual Doppler effect that both inverse Doppler effect and normal effect appeared simultaneously are proved successfully.

Discussion
In conclusion, we have found a novel phenomenon, the dual Doppler effect, with both normal and inverse Doppler effect occurring simultaneously in one moving two-dimensional wedge type PhC. It was obtained by a dynamic two-dimensional FDTD method which can deal with the situation of electromagnetic wave propagation in moving objects. The theoretical values showed a high consistence with the detecting values along with a relative error less than 0.2%. In addition, we have shown that the outgoing direction of the complex harmonic coupling to the vacuum can be predicted precisely by using multiple reciprocal unit cells. The novel phenomenon could benefit any technique currently exploiting the conventional Doppler effect, such as the Doppler cooling of atomic gases, radar deception, invisibility cloaks, microstructure dual frequency interferometer 19,20 and so on. The technique is also applicable to systems with tunable sources of terahertz radiation and other frequency regions that are difficult to access.

Methods
Spatial Fourier Transform. The obtained optical field in the structure is a complex field A·exp(iϕ). A complex Spatial Fourier Transformation (SFT) is used to retrieve both the wavevector (k) and amplitude (A) of the field propagating. The SFT of the complex optical field also reveals the sign of k, i.e. both positive and negative wavevectors. As shown in the following equation, the convolution operation shifts the spatial frequencies to the right positions along the axis in the spectrum.

Dynamic FDTD algorithm.
To investigate the performance of the Doppler effect, a dynamic algorithm is developed based on the traditional FDTD method. In the general FDTD algorithm, the position of object is represented by assigning special ε and μ to the corresponding grids. It cannot be changed during the whole iterative operation, making it useless in dealing with the case of wave propagation in a moving object. The dynamic FDTD method can solve this problem. As shown in Fig. 6, the improved FDTD algorithm can be simply described as: the parameters of ε and μ are reassigned to the adjacent space grids before each iterative operation, so as to reappear the movement of the object. For example, the object moves in a speed of v, the time step is Δt and the grid size is Δx, then the parameters of ε and μ will be reassign to the grids vΔt/Δx away from the present girds. Before each reassigned operation, the electromagnetic field data have been saved, and after the reassignment, the saved field data is used as the initial field data for the coming iterative operation.