To verify the propagation characteristics, the proposed structure was modeled in a 3D grid by employing commercially available software42 based on the finite-difference time-domain (FDTD) method. Perfectly matched layer (PML) boundary conditions were used to terminate the computational domain and a TM polarized (electric field is in the y- direction) Gaussian pulse with a pulse length equal to ct/a = 30 in normalized time units was launched. The source was placed at a distance of 0.2a in front of the structure. Fig. 3(a–c) show the electric field intensities, calculated along the propagation (z-) direction, at the center point of the xy- plane at each time frame for normalized central frequencies a/λ = 0.450, 0.468 and 0.496, respectively. Figures 3(a–c) suggest that the PhC exhibits distinct localization regions for each frequency. More precisely, the incoming wave is slowed down until it reaches its localization region (or turning point) and is then trapped for a finite time interval before it starts propagating backwards. If the trapping time is defined as the time duration of the electric field intensity decaying from its maximum value until its 1/e 19, then the normalized trapping times for Figs 3(a–c) would be obtained as ct/a = 18.9, 20.3 and 24.1; respectively. Further inferences can be made from the time evolution of the pulses, if one pays attention to the change of the magnitude of the field intensity. At the localization region the intensity increases, since the local energy density is expected to be inversely proportional to the group velocity20, 21, 43. Furthermore, to obtain a broadband analysis of the field enhancement, the steady-state electric field intensity was numerically calculated at the same region (along the z- propagation direction, at the center of the xy- plane), and is shown in Fig. 3(d). If the turning point is defined as the spatial position where the intensity becomes maximum, then a nearly linear increase of the position of the turning point is observed, which originates from both facts that the upper cut-off region is nearly linearly dependent on the layer-to-layer distance (see Fig. 1(b)) and that the chirp function is varying linearly. One interesting observation from Fig. 3(d) is that the intensity enhancement is highly sensitive to the operation frequency. In particular, for the normalized operation frequency range a/λ = {0.42, 0.50} the maximum normalized intensities occurring at the localization regions deviate between 8 and 132. We attribute this deviation to the oscillations occurring between the turning point and the entrance of the structure. For instance, Fig. 3(a–c) suggest that the intensity enhancement is nearly equal for these frequencies, which is expected since the nearly linear cut-off variation will induce similar slowdowns for different frequencies, even though their steady-state intensity enhancements are much diversified, as indicated with white arrows in Fig. 3(d). The main difference between these two cases is that the steady state analysis is the response of a continuous source excitement, which comprises multiple enhancement factors such as slowing down and multiple interferences of the reflected waves from the turning point and at the entrance of the structure; whereas the time evolution analysis is the response of a source with a finite temporal pulse width, which includes only the slowing mechanism in terms of the intensity enhancement. We therefore believe that the oscillatory behavior of the intensity enhancement stems from the peculiar relation between the optical path and the local wave vector distribution. In this regard, for a specific frequency, the local wave vector distribution leads to such a local wavelength distribution that it will encounter constructive interferences within the optical path (which is also unique to the operational frequency). In other words, the slowing down of the wave will be accompanied by a structural resonance, which in result give rise to an additional enhancement factor. Such frequencies correspond to the narrowband peaks in Fig. 3(d) and one can infer from this figure that the spectral distance between these peaks decreases as the operational frequency increases. This is also reasonable, since increasing the frequency will cause the wave to localize and reflect further away from the entrance of the structure and to enter the structure with a wave vector away from the cut-off region, compared to lower operation frequencies. This will have the consequence that the wavelength distribution will initiate with lower spatial wavelengths as the wave enters the structure21, in which case it will be enough to increase the operational frequency less compared to lower frequencies, to match up the next constructively interfered frequency.Figure 3Instantaneous and steady-state field analysis. Time evolution of a Gaussian pulse propagating inside the 3D PhC (along the z- propagation direction, at the center of the xy- plane), for three frequencies: (a) a/λ = 0.450, (b) a/λ = 0.468 and (c) a/λ = 0.496. (d) Broadband steady-state field intensity distributions at the same propagation line. Steady-state field distributions obtained at the xz- plane of the center of the y- direction for the normalized operational frequencies equal to (e) a/λ = 0.450, (f) a/λ = 0.468 and (g) a/λ = 0.496. The magnitude of the field intensities are normalized by dividing the calculated intensities to the intensity produced by the source at the specific frequencies. The distances in the z- propagation direction are given with respect to the position of the source.