Frozen Mode Regime in Finite Periodic Structures
Huanan Li Ilya Vitebskiy Tsampikos Kottos





Department of Physics, Wesleyan University, Middletown CT-06459, USA

Air Force Research Laboratory, Sensors Directorate, Wright-Patterson Air Force Base, OH- 45433, USA

Department of Physics, Wesleyan University, Middletown CT-06459, USA
KBRwyle, Dayton, OH 45431



Periodic structures with Bloch dispersion relation supporting a stationary inflection point (SIP) can display a unique scattering feature, 
the frozen mode regime (FMR). The FMR is much more robust than common cavity resonances; it is much less sensitive to the boundary 
conditions, structural imperfections, and losses. Using perturbation theory, we analyze the FMR in the realistic case of a finite fragment 
of a periodic structure. We show that in close proximity of SIP frequency, the character of the FMR is qualitatively different from the
known case of a semi-infinite structure. 


Frozen Mode Regime in Finite Periodic Structures 
    Tsampikos Kottos
    December 30, 2023
=================================================



Introduction –The ability to engineer composite structures with predefined wave dispersion relation is one of the greatest 
achievements of the last thirty years <cit.>. An outgrowth of this technological development was the realization of 
photonic and phononic band-gap materials which are now used routinely to control the propagation of light and sound by creating 
stop bands and by adjusting and reducing the wave's group velocity. Direct consequences of the latter accomplishment is the 
enhancement of linear (like absorption or amplification) and nonlinear (like Kerr) effects, and the efficient manipulation of 
direction of electromagnetic or acoustic signals.

There are many ways that one can achieve a slow wave propagation, where the group velocity v_g=∂ω/∂ k
|_k=k_0.≈ 0, via dispersion management. Example cases include degenerate or regular band edges 
corresponding to frequency dispersion relations ω(k)-ω_0∝ (k-k_0)^2m with an integer m>1 or m=1 
respectively. Another type of slow wave is associated with stationary inflection points (SIP) singularities ω(k)-ω_0
∝ (k-k_0)^2m+1. The latter results in the formation of the so-called frozen mode regime (FMR) <cit.>, whose 
most prominent feature is the nearly total conversion of an input signal into a slow (frozen) mode and dramatically enhanced amplitude. 
This feature has to be contrasted with the situations encountered in (degenerate) band-edges where the incident wave is typically 
reflected at the interface between free space and a slow-wave medium. 

The FMR is not a resonance phenomenon; it is not particularly sensitive to the size and shape of the 
photonic/phononic structure and it can tolerate all kinds of structural imperfections. On top of it, the FMR 
can withstand much stronger losses than any known cavity resonance. The above features make the FMR very 
attractive for a variety of applications in optics, microwave, RF and acoustics. In a combination with non-reciprocity, the FMR 
can lead to the phenomenon of electromagnetic unidirectionality <cit.>.  In a nonreciprocal photonic 
structure with gain, the FMR can also result in a cavity-less unidirectional lasing <cit.>. 
Due to the underlying mathematical complexity, the FMR has been fully analyzed only in semi-infinite 
periodic multilayered structures <cit.> and multimode waveguide arrays <cit.>. 


Here we are investigating the scattering problem for a finite multi-mode structure, whose periodic counterpart display 
SIPs. Using an abstract transfer matrix formalism 
together with a matrix perturbation approach, we studied the transmission characteristics 
of such set-ups in the  FMR. We derived theoretical expressions for the energy flux carried by the slow propagating 
mode(s) and identify a new scaling behavior with respect to the frequency detuning. Specifically, we find that the energy fluxes, 
associated with the slow propagating mode(s), undergoes a transition at critical sample lengths L_ C∝ |ω-
ω_ SIP|^-1/3 from an S_ p∼ O(1) behavior (characteristic of semi-infinite structures) 
to an S_p∝|ω-ω_ SIP|^-2/3 law. The latter divergence is balanced by a simultaneous development 
of an energy flux carried by pairs of evanescent modes. Our results are confirmed via detailed simulations for set-ups with both 
symmetric and asymmetric spectrum that support SIPs.





Transfer Matrix Formalism near SIPs– We consider (finite) periodic composite structures whose infinite counterpart has a 
dispersion relation ω(k) which supports a SIP at some frequency ω_ SIP. In the absence of any non-reciprocal 
elements, the dispersion relation is reciprocal ω(k)=ω(-k) and therefore an  ω_ SIP is associated with 
two counter-propagating slow modes at ± k_ SIP. In such structures, two sets of three modes are responsible for SIPs at 
± k_ SIP.

The wave propagation can be analyzed using the transfer matrix approach. The transfer matrix T(z, z_0;ω) connects the wave amplitudes (in mode space– see Ref. <cit.> for a coupled mode theory 
implementation) Ψ of a monochromatic wave at two different positions z and z_0 through the relation Ψ(z)= T(z, z_0;ω)Ψ(z_0). For the specific case of periodic structures, the transfer 
matrix of a unit cell 𝒯(ν)≡𝒯(1, 0;ω_0+ν) dictates the transport. 
Here we assumed that the length of the unit cell is L_ uc=1, ω_0=ω_ SIP and ν is the frequency 
detuning. 

We consider a minimal model for which the unit transfer matrix 𝒯(ν) is 6×6 and it is analytic around the 
SIP. Since a symmetric 
spectrum develops two SIPs at ν=0, 𝒯(0) can be represented by its Jordan normal form as

    𝒯(0)=g_0(0)([ J^-   0;   0 J^+ ])g_0^-1(0);
    J^±≡ e^± k_0([ 1 1;   1 1;     1 ]),

where g_0(0)=[j_0^-,j_1^-,j_2^-,j_0^+,j_1^+,j_2^+] is an invertible 6× 6 matrix with 
columns given by the Jordan basis vectors and ± k_0=± k_ SIP. When ν 0 (but still ν→0), 
𝒯(ν) reduces to its normal form <cit.>

    g_0(ν)^-1𝒯(ν)g_0(ν)=    ([ T^-(ν)      0;      0 T^+(ν) ]),

where T^±(ν)=J^±+T_1^±ν+⋯≡ e^± k_0(I_3+Z^±
(ν)) and the matrix g_0(ν) depends analytically on ν in the vicinity of ν=0. 

Next we focus on the eigenvalue problem associated with the individual blocks of the normal form Eq. (<ref>). 
Let us consider, for example, the block matrix T^-(ν) or its equivalent 
problem associated with the matrix Z^-(ν). Simple-minded normal perturbation theory is not useful in cases like 
ours when the leading term of the operation expansion is nilpotent i.e.  (Z^-(0))^3=0. Indeed in such cases
the standard Taylor series assumed for the eigenvalues is not the appropriate expansion; rather one has to develop the eigenvalue
perturbation expansion using a Puiseux series <cit.>. Nevertheless, a singular perturbation theory provides a recipe 
to “re-construct" the appropriate operator expansion after identifying the correct leading order term <cit.>. Using this approach 
we have found that Z^-(ν)=Z_0^-(ν̃)+Z_1^-ν̃+⋯ where  Z_0^-(ν̃)≡([ 0 1 0; 0 0 1; 0 0 0 ])+ν̃([ 0 0 0; 0 0 0; 1 0 0 ]) and ν̃≡[Z^-(ν)]_31=-3!/ω^”'(-k_0)ν+𝒪(ν^2). 

The diagonalization of Z_0^-(ν̃), gives  

    G_0^-1(ν̃)Z_0^-(ν̃)G_0(ν̃)=  ν̃^1/3Λ_0; Λ_0=diag(c_0,c_1,c_2),

where c_n≡ e^2π/3n and the similarity transformation matrix G_0(ν̃) is a 
Vandermonde matrix <cit.> of order 3.
Further, the diagonalization process for Z^-(ν̃)
(or equivalently T^-(ν)) can be continued order-by-order, leading to the following compact form

    e^-SG_0(ν̃)^-1T^-(ν)G_0(ν̃)e^S
    
    =     e^- k_0(I_3+ν̃^1/3Λ_0+ν̃^2/3Λ_1+⋯),

where the matrix S≡ S(ν̃^1/3)=ν̃^1/3S_1+ν̃^2/3S_2+⋯
and Λ_1,⋯ are diagonal matrices.  A similar treatment applies for the eigenvalue problem associated with T^+(ν).

The above approach allows us to evaluate perturbatively the eigenvalues θ_n^∓(ν) and the eigenvectors 
f_n^∓(ν) of the unit transfer matrix 𝒯(ν). We get

    𝒯(ν)f_n^∓(ν)=    θ_n^∓(ν)f_n^∓(ν), n=0,1,2
    
    θ_n^∓(ν)≈    e^(∓ k_0+λ_n^∓); λ_n^∓(ν)≡α_0^∓c_nν^1/3
    
    f_n^∓(ν)≈   [1-σ_2^∓λ_n^∓+η^∓(λ_n^∓)^2]j_0^∓
    
    +    [λ_n^∓-σ_1^∓(λ_n^∓)^2]j_1^∓-
    (λ_n^∓)^2j_2^∓

where α_0^∓=(3!/ω^”'(∓ k_0))^1/3,
η^∓=γ_3^∓-γ_1^∓+1/2((σ_1^∓)^2+σ_1^∓σ_2^∓-(σ_2^∓)^2),
 σ_l^∓=1/3[T_1^∓]_l+1,l/[T_1^∓]_31,
γ_l^∓=1/3[T_1^∓]_l,l/[T_1^∓]_31 and 
j_n^∓ is the Jordan basis of T(0). 

We assume that ν→0^+ and that the incident wave is entering the finite structure 
from the left interface at z=0; for an example see the dispersion relation in Fig. <ref>b. We are now ready to decompose 
any wave inside the structure to the forward (backward) propagating f_0^- (f_0^+) and evanescent f_1^-,f_2^+ (f_1^+,f_2^-) 
modes and thus evaluate the associated conversion coefficients. We shall also analyze the energy flux carried from these modes and 
determine their scaling with respect to detuning ν. 


Conversion Coefficients – We consider that the finite structure consists of N periods of the unit cell. In contrast to the semi-infinite 
case <cit.>, finite structures involve two interfaces at z=0 and z=N and therefore both forward and backward modes can 
participate in the scattering process. When ν→0^+, the eigenmodes Eq. (<ref>) associated with different blocks 
in Eq. (<ref>) become degenerate within the specific
block. This observation forces us to construct a new “well-behaved" basis ℬ_fb={ f_0^-, f̃_1^-,
 f̃_2^-, f_0^+, f̃_1^+, f̃_2^+}, where the new basis vectors f̃_1^∓=
f_1^∓-f_0^∓/α_0^∓ν^1/3(c_1-1) and f̃_2^∓=-1/3(α_0^∓)^2ν^2/3(c_2f_2^∓+c_1f_1^∓+f_0^∓) together with
f_0^±, are independent in the limit of ν→0^+. 

Next we introduce semi-infinite leads and coupled them to the left and right of the structure. We shall assume that the leads do not develop 
any spectral singularity around ω_0. We then request continuity of  Ψ(z) at the interfaces at z=0,N together with the scattering condition 
that the incident wave enters the structure from the left i.e. that 
the coefficients of the backward modes on the right leads are zero. Finally the identification of the appropriate (non-degenerate in 
the ν→0 limit) basis guarantee that the scattering problem has unique solution and that the expansion coefficients {φ_1^+,φ_2^+,φ_3^-,φ_1^-,φ_2^-,φ_3^+} of Ψ(z=0^+) in the basis ℬ_fb exists for any incident wave. 

Obviously the specific values of the expansion coefficients depend on the particular form of the incident wave. Nevertheless 
some features are independent of the incident waveform; we find, for example, that φ_l^±(ν)=
φ_l^±(0)+𝒪(ν^1/3) while the envelopes scale as |φ_3^±(0)|∼𝒪(N^-1) and |φ_j^±(0)|∼𝒪(N^0), j=1,2., in the large N limit. Correspondingly, in terms of the eigenmodes of T(ν), the 
expansion of Ψ(z=0^+) is given as

    Ψ(z=0^+)=∑_σ=±[-φ_3^-σ/3(λ_0^-σ)^2+-φ_2^σ/(c_1-1)λ_0^-σ+φ_1^σ]f_0^-σ
       +[-c_1φ_3^-σ/3(λ_0^-σ)^2+φ_2^σ/(c_1-1)λ_0^-σ]f_1^-σ+[-c_2φ_3^σ/3(λ_0^σ)^2]f_2^σ,

where σ=+/- correspond to forward/backward modes. Substitution of the scaling expressions for the expansion 
coefficients φ_l^± together with λ_0^-σ (see Eq. (<ref>)) in Eq. (<ref>), 
allow us to estimate the scaling of the conversion coefficients. Specifically we find that each of the square bracket terms
in Eq. (<ref>) scale as 

    [⋯] β_2/(Nν^2/3)+ β_1/ν^1/3;
 
where β_1,2 are some constants independent of N and ν.
Equation (<ref>) signifies a scaling transition from  1/ν^2/3 (small samples) to 1/ν^1/3 (large samples)
at some critical sample length 

    L_ C= L_ uc N_ C∝ L_ uc/ν^1/3.

While the latter scaling law for the conversion coefficients is already known from the case of semi-infinite structures, the 
former one is completely new and a trademark of the finite length nature of the scattering setting.  

Modal energy flux –  We now turn our focus on the consequences of the scaling (<ref>) in the modal energy 
flux. First we recall that near a SIP the Bloch dispersion relation takes the form ω -ω _0∝( k-k_0) ^3. The group velocity of the slow propagating mode(s) is

    v_g=∂ω/∂ k∝( k-k_0) ^2∝( ω -ω _0) ^2/3

while the associated energy flux contribution S_p is

    S_p=W_pv_g∝ W_p ν ^2/3

where W_p is the energy density of the slow propagating mode. An  estimation for
the scaling of W_p is provided from the behavior of the conversion coefficients associated with f_0^±, see 
Eqs. (<ref>,<ref>), i.e. W_p∝|β_2/Nν^2/3+β_1/ν^1/3|^2. In other words, W_p 
undergoes a transition from an 1/ν^4/3 (for N<N_ C) to an 1/ν^2/3 (for N>N_ C) scaling 
with respect to the detuning ν. 

In the latter limit of “semi-infinite” samples the sole contribution to the energy flux comes from the slow mode and thus 
S=S_p=W_pv_g∝ 1, as expected also from previous studies <cit.> (see also Appendix). 
In contrast, in finite scattering set-ups, the contribution S_p from the slow propagating mode(s) to the 
total energy flux S is  

    S_p=W_pv_g∝ W_pν^2/3∝ν ^-2/3

where we have used Eq. (<ref>) together with the scaling behavior of W_p for short samples.

The anomalous scaling Eq. (<ref>) of the modal energy flux of the propagating modes near the SIP 
can be balanced only by the same type (but different in sign) of divergence of modal energy flux S_ ev carried 
by the two pairs of forward and backward evanescent modes. This is necessary in order to get a total energy flux S∼ O(1) and it is a new feature associated with the fact that the scattering set-up is finite.
In the remaining of this paper we will be checking these theoretical predictions using 
some simple numerical examples.


Tight-binding model – We first consider a tight-binding (TB) model supporting a symmetric dispersion relation
with two symmetric SIPs, see Fig. <ref>a,b. This system can be realized as a quasi-one-dimensional 
array of coupled resonator
<cit.>. The system consists of M=3 chains of coupled resonators where the resonators 
of each chain have equal nearest-neighbor coupling (set to be 1 as coupling unit). The vertical inter-chain 
coupling between the nearest chains is κ_1. In addition, the resonators at the first two chains have an on-site 
potential contrast κ_0 (with respect to the resonators of the third chain) and they are also coupled via an inter-chain 
diagonal coupling κ_2. In this TB model a monochromatic electromagnetic wave is described by  

    ω E_l^(1)=     E_l-1^(1)+E_l+1^(1)+κ_1E_l^(2)+κ_2E_l+1^(2)
    +κ_0E_l^(1)
    ω E_l^(2)=    
    E_l-1^(2)+E_l+1^(2)+κ_1(E_l^(1)+E_l^(3))+
    κ_2E_l-1^(1)+κ_0E_l^(2)
    ω E_l^(3)=     E_l-1^(3)+E_l+1^(3)+κ_1E_l^(2),

where E_l^(m) is the field amplitude at the site l of the chain m. Substituting E_l^(m)=
A^(m)e^ kl in Eq. (<ref>), we get 

    ω u_A=Du_A;  
    D≡([   ϵ(k)   v(k)      0; v^*(k)   ϵ(k)    κ_1;      0    κ_1 2cos k ])
 
where ϵ(k)=2cos k+κ_0,   v(k)=κ_1+κ_2e^ k and u_A=(A^(1), A^(2), 
A^(3))^T. Then the dispersion relation ω(k) is obtained by setting (D-ω I_3)=0. 
Generally there are three bands for this model and we mainly focus on the band supporting SIPs characterized by ω'
(± k_0)=ω”(± k_0)=0 and ω”'(± k_0)≠0. An example is given in
Fig. <ref>b, where for the parameter values κ_0=κ_1=κ_2=5 and SIPs at ± k_0= ±π/2 and ω_0=-5.  

The scattering sample is attached to the left and to the right with semi-infinite leads, which are composed of three decoupled 
chains with constant nearest-neighbor coupling κ_L in each chain. Thus the leads support a traveling wave whenever
its frequency is within the band ω(k_L)=2κ_Lcos k_L where -π≤ k_L< π. The field amplitude in each 
lead-chain can be written as a sum of two counter-propagating waves, i.e., E_l^(m)=a^(m)e^|k_L|l+
b^(m)e^-|k_L|l. In the simulations, we assume κ_L=4 so that b^(m) 
represents the amplitude of incoming waves since v_g≡∂ω/∂ k_L|_-|k_L|.>0. 





Finally the energy flux through a section l in the scattering domain can be defined using the continuity equation, 

    d/dt[∑_mE_l^(m)*E_l^(m)]=F_l-1→ l-F_l→ l+1

where F_l-1→ l≡2ℐm[∑_mE_l-1^(m)E_l^(m)*+
E_l-1^(1)κ_2E_l^(2)*] denotes the flux flowing from section l-1 to l. At the
same time the field amplitudes can be parametrized as E_l-1^(m)=a_l-1^(m)+b_l-1^(m) and E_l^(m)
=a_l-1^(m)e^ k+b_l-1^(m)e^- k, where ω=2cos k and k is in general a complex number. 
The self-consistency requirements impose E_l^(m)≡ a_l-1^(m)e^ k+b_l-1^(m)e^- k=
a_l^(m)+b_l^(m), which together with Eq. (<ref>) allows us to calculate the unit transfer matrix 𝒯(ν) such that Ψ_l=𝒯(ν) Ψ_l-1 where  
Ψ_l≡(a_l^(1),b_l^(1),a_l^(2),b_l^(2),a_l^(3),b_l^(3))^T.  

We are now ready to analyze numerically the scaling of the modal energy fluxes of the TB model Eq. (<ref>). First 
we have verified that Eq. (<ref>) is valid for 𝒯(0) using the aforementioned parameters. In 
Fig. <ref>c we report our numerical findings for the modal energy flux associated with the slow propagatings S_p 
and evanescent S_ ev modes for three different system sizes N. We find that while for ν→ 0 these 
quantities scale according to the new scaling law Eq. (<ref>), the modal fluxes saturate to a constant value at
different ν_ C∝ 1/N^3 in accordance to our theoretical prediction, see Eq. (<ref>). In Fig. <ref>d  
we report the data for one of the N values referring to a linear-linear plot. We find that  
S_ ev^∓, associated with each of the two pairs of evanescent modes (corresponding to the T^∓ blocks
in Eq. (<ref>)), balances the divergent of the S_p contribution so that the total flux S∼ O(1).

Non-reciprocal layered structures – It is straightforward to reproduce Eqs. (<ref>,<ref>,
<ref>) for the case of finite set-ups with spectral non-reciprocity i.e ω(k)≠ω(-k). Here, instead, we confirm 
numerically the validity of these equations for the example case of a multilayer periodic magnetic photonic crystal (PC) with 
proper spatial arrangement, see Fig. <ref>a <cit.>. When analyzing the modal energy 
flux, we find that the forward slow propagating mode carries an energy flux which scales according to Eq. (<ref>) 
while the pair of the associated evanescent modes balanced this divergence in a similar manner, see Fig. <ref>c,d. The 
energy flux of the remaining (fast) backward propagating mode does not show any divergence and has minimal contribution 
to the total energy flux, see Fig. <ref>d.



Conclusions - We developed a theory of FMR for realistic finite structures. We find that the character of the FMR undergoes a 
transition which is reflected in a dramatic change of the scaling behaviour (with respect to detuning ν) of the modal energy flux 
of the slow propagating modes at critical lengths L_C∝ 1/ν^1/3. As opposed to the semi-infinite case, below this length
-scale, the energy flux is carried even by (pairs of) evanescent modes. Our results might have important applications to non-reciprocal 
transport.

Acknowledgments - (I.V.) acknowledges support from AFOSR via LRIR14RY14COR.




1
D13P. A. Deymier, Acoustic Metamaterials and Phononic Crystals, (Springer, 2013)

JJWM08J. D. Joannopoulos, S. G. Johnson, J. N., Winn, & R. D. Meade, Photonic Crystals: Molding the Flow of Light, 2nd ed. (Princeton University Press, Princeton, NJ, 2008).

FV06A. Figotin and I. Vitebskiy, Waves Random Complex Media 16, 293 (2006).

FV11A. Figotin and I. Vitebskiy. Slow Wave Phenomena in Photonic Crystals (Review article), 
Laser & Photonic Reviews 5, 201 (2011).

FV13A. Figotin and I. Vitebskiy, Electromagnetic Unidirectionality in Magnetic Photonic Crystals; 
Book chapter in: Magnetophotonics: From Theory to Applications. Springer Series in Materials Science, 2013.

RKVK14H. Ramezani, S. Kalish, I. Vitebskiy, and T. Kottos, Phys. Rev. Lett. 112, 043904 (2014). 

FV01A. Figotin, I. Vitebskiy, Phys. Rev. E 63, 066609 (2001).

FV06a A. Figotin, I. Vitebskiy, J. Magn. Mater 300, 117 (2006)

FV06bA. Figotin, I. Vitebskiy, Phys. Rev. E |bf 74, 066613 (2006)

GSSB12 N. Gutman, C. Martijn de Sterke, A. A. Sukhorukov, L. C. Botten, Phys. Rev. A 85, 033804 (2012).

GDSSS12 N. Gutman, W. Hugo Dupree, Y. Sun, A. S. Sukhorukov, C. M. de Sterke, Opt. Express 20, 3519 (2012).

SLCPK08A. Sukhorukov, A. Lavrinenko, D. Chigrin, D. Pelinovsky, Y. Kivshar, J. Opt. Soc. Am. B 25, C65 (2008).

CMTW-P Huang, J. Opt. Soc. Am. A 11, 963 (1994).

K95T. Kato, Perturbation theory of  linear operators, Springer (1995).

MPT V. N. Bogaevski and A. Povzner, Algebraic Methods in Nonlinear Perturbation Theory (Springer-Verlag, Berlin, 1991). 

LT85P. Lancaster, M. Tismenetsky, The theory of matrices, Academic Press (2985).

KMMVK16U. Kuhl, F. Mortessagne, E. Makri, I. Vitebskiy, T. Kottos, submitted (2016).

SFKMS17T. Stegmann, J. A. Franco-Villafane, U. Kuhl, F. Mortessagne, T. H. Seligman, Phys. Rev. B 95, 
035413 (2017).

BKMM13M. Bellec, U. Kuhl, G. Montambaux, F. Mortessagne, Phys. Rev. B 88, 115437 (2013).









Supplemental Materials

















§ SEMI-INFINITE STRUCTURES


We can also consider semi-infinite scattering set-ups for the case of coupled waveguide arrays, see Fig. 1a,b. 
In this case the semi-infinite scattering domain is attached
to one sem-infinite lead (say to the left). Again we assume that the lead does not develop any singularity in its 
spectrum around ω=ω_0. We consider that the incident wave is sent from the lead towards the 
scattering domain, thus exciting the forward mode(s). In this case the forward mode(s) consist of one forward 
slow propagating mode f_0^-(ν) and two forward evanescent modes f_1^-(ν) 
and f_2^+(ν). 

It is important to notice that in the limit ν→0^+, the modes f_0^-(ν) and 
f_1^-(ν) are degenerate because f_0^-(ν)-f_1^-(ν)→0,
see Eq. (5). Following the same strategy as in the case of finite structures, we construct a “well behaved”
basis for the forward modes as ℬ_f={ f_0^-(ν), f̃_1^-(ν), f_2^+(ν)} with f̃_1^-(ν)≡f_1^-(ν)-f_0^-(ν)/α_0^-ν^1/3(c_1-1), which are linearly independent as ν→0^+. 

Next we decompose the propagating waves in this basis. At the interface z=0, the wave state Ψ is assumed to 
be continuous. By matching the boundary condition, the expansion coefficients {φ_1, φ_2, φ_3} of Ψ(z=0^+) in the basis ℬ_f are obtained. Although these coefficients 
depend on the specific form of the incident waves, they, nevertheless, are characterized by some general features i.e. φ_l(ν)=φ_l(0)+𝒪(ν^1/3).
As results, in terms of forward Bloch modes, we obtain the expansion of Ψ(z=0^+) as 

    Ψ(z=0^+)=    [φ_1-φ_2/α_0^-(c_1-1)ν^1/3]f_0^-
    
    +    [φ_2/α_0^-(c_1-1)ν^1/3]f_1^-+φ_3f_2^+.

This result is similar to the one obtained in layered structures. The only differences is the additional forward evanescent mode originating
from the T^+ block in Eq. (2). This mode has relatively negligible contribution to the field when ν→0^+.

Following the same argumentation as in layered structures we can estimate, for the semi-infinite systems, the scaling behaviour of the modal
energy flux of the slow propagating mode. Specifically, using Eq. (<ref>) we can evaluate the wave Ψ( z) 
transmitted to lossless semi-infinite periodic structure. The latter is composed of propagating Ψ_p( z) and evanescent 
Ψ _ev( z) contributions. An estimation of the scaling of these components with the detuning ν is given from the 
conversion coefficients associated with f_0^- (propagating amplitude) and f_1^- (diverging evanescent amplitude), see Eq. (<ref>).
We get

    |Ψ _pr| ^2∝|Ψ
    _ev| ^2∝ν^-2/3.

At the interface z=0 of the semi-infinite structure, the two diverging Bloch components interfere destructively, thereby satisfying the continuous
boundary conditions at z=0. As the distance z from the interface increases, the evanescent contribution Ψ _ev( z)
vanishes, while the diverging propagating component Ψ _p(z) provides the sole contribution to the transmitted field Ψ( z) and determines its saturation value. In short, in the vicinity of a SIP, the wave transmitted to the semi-infinite periodic
structure has diverging saturation amplitude (<ref>) and vanishing group velocity. The latter can be easily calculated
 from the dispersion relation which in the vicinity of a SIP is

    ω (k)-ω _SIP∝ (k-k_SIP)^3.

resulting in a vanishing group velocity which scales as

    v_g∝( ω -ω _SIP) ^2/3

Using Eqs. (<ref>,<ref>) we get that the energy flux is

    S∝ v_g|Ψ _pr| ^2

provided by the frozen mode remains finite even at ω =ω _SIP. Under certain conditions, the energy flux of the transmitted frozen mode can
be close to that of the incident wave, implying effective coupling with the incident wave.





§ SIMULATION DETAILS OF MULTI-LAYERED STRUCTURE


The basic unit of the PC shown in Fig. 2a, contains three components consisting of two nonmagnetic misaligned anisotropic layers (red and blue) 
with a magnetic layer (grey) in-between. 
The magnetic layer – in the presence of static magnetic field or spontaneous magnetization – guaranties the violation of time-reversal 
symmetry while the anisotropic layers break the mirror reflection symmetry; a necessary condition for achieving spectral non-reciprocity 
i.e ω(k)≠ω(-k). We can find appropriate parameters for which the spectrum supports one SIP (see Fig. 2b). 

The numerical analysis for the multilayered structure was performed using the transfer matrix approach (for a detail 
presentation see Ref. [6]). The permeability tensor for the gyrotropic magnetic layer has the form

    μ̂_f = ( [ μ_xx  i β    0; -i β μ_xx    0;    0    0    1 ])
 
while the permittivity tensor is

    ϵ̂_f = ( [  ϵ_f  i α    0; -i α  ϵ_f    0;    0    0    1 ])
 
where α is the gyrotropic parameter responsible for Faraday rotation. The width of this layer is L_f.

Similarly the permeability tensor for the birefrigent layers is a 3×3 identity matrix 
μ̂=  1 while the  permittivity tensor takes the form

    ϵ̂_1,2 = ( [ ϵ_A + δcos (2 ϕ_1,2)       δsin (2 ϕ_1,2)                    0;       δsin (2 ϕ_1,2) ϵ_A - δcos (2 ϕ_1,2)                    0;                    0                    0                    1 ]),

where ϵ_A is the permittivity of the medium. Moreover δ is the magnitude of 
the in-plane anisotropy and ϕ_1,2 is the angular orientation of the principal axes in the xy-plane.  The width of these layers is L_1,2=L.  

In our simulations, we have placed the multilayered structure in air and we have used the following parameters for ϵ_f≈ 3.9, α≈ 0.9, μ_xx≈0.814,
β≈ 0.015, ϵ_A≈ 7.3, δ=0.54, L_f=1/2 L and L=1. Finally the we have used ϕ_1 = π/4 
and ϕ_2=0.  The SIP (see Fig. 1c) was found at ω_0=2.5c/L_ uc where L_ uc=2.5.











