Propagating and Evanescent Bessel Light Beams in \( \mathcal{PT} \)-Symmetric Systems

1. Introduction

Parity-time ($\mathcal{PT}$)-symmetric system as a sort of a non-Hermitian system was first introduced more than two decades ago within the framework of quantum mechanics [1,2] and have since enabled the exploration of numerous remarkable effects in optics, such as unidirectional invisibility [3,4], simultaneous coherent perfect absorption and lasing (CPA laser) [5,6,7,8], asymmetric optical pulling and pushing forces [9], phase transitions between regimes of coherent amplification and absorption [10], enhanced photonic spin Hall effect near exceptional points [11], and so on.

$\mathcal{PT}$-symmetric systems are inherently non-Hermitian, what manifests itself through the properties of either Hamiltonian $\widehat{H}$ or scattering matrix $\widehat{S}$ [7,12,13]. In particular, for nonmagnetic optical media, the $\mathcal{PT}$-symmetry dictates accomplishment of the condition $\epsilon \left(\mathbf{r}\right)={\epsilon }^{*}(-\mathbf{r})$ [3], which is most commonly implemented in structures composed of alternating layers of equal thickness $d_1=d_2=d$ with opposite signs of the imaginary part of the dielectric permittivity as $\mathrm{Im}\left({\epsilon }_1\right)=-\mathrm{Im}\left({\epsilon }_2\right)$. A distinctive feature of such systems is a possibility of spontaneous $\mathcal{PT}$ symmetry breaking at an exceptional point (EP) defined by coalescing eigenvalues and eigenvectors of the scattering matrix or Hamiltonian [13,14,15,16]. The topological nature of these transitions is intimately related to EPs inducing unique phenomena such as eigenmode exchange and phase jumps in reflection coefficients [17]. For propagating waves, the eigenvalues of the scattering matrix S must be unimodular as $|s_1|=|s_2|=1$ in the $\mathcal{PT}$-symmetric regime, whereas the eigenvalues are reversed as $|s_1|=1/|s_2|$ in the $\mathcal{PT}$ symmetry broken regime, the eigenstates being amplified and attenuated [7]. Evanescent waves exhibit an alternative means for insightful analysis of the $\mathcal{PT}$ symmetry through tailoring the transverse wave-vector component as an additional parameter governing the $\mathcal{PT}$-symmetric state of the system [18,19,20]. It must be noted, that no unique definition of the scattering matrix exists when it is introduced independently of system’s Hamiltonian [21]. Since the scattering matrix relates the amplitudes of incoming and outgoing waves, there is a certain freedom in choosing an order of these amplitudes and, accordingly, in defining the scattering matrix itself. In this work, we consider two such definitions, $S_1$ and $S_2$, where the first of them is more convenient for predicting the lasing threshold, while the second one better corresponds to the PT-symmetric Hamiltonian [15,21].

A plane wave has traditionally been a regular object for studying $\mathcal{PT}$ symmetry and has enabled the investigation of numerous physical phenomena [3,7,22], with many more yet to be discovered. However, another tool in the context of the $\mathcal{PT}$ symmetry, a Bessel beam [23,24], has received comparatively little attention, particularly evanescent beams [19,25,26]. Bessel beams are nondiffracting solutions of the wave equation in cylindrical coordinates, formed by a superposition of plane waves whose wave vectors lie on the surface of a cone with half-cone angle $\theta$ [23]. Due to their wavefront, they inherently exhibit rich polarization structure, including $s-$, $p-$, transverse, and longitudinal polarizations [27]. The transverse wave-vector component $k_t=kq$ determines the beam profile as shown in Figure 1. For $q<1$ the beam is propagating; for $q>1$, the wave-vector cone becomes imaginary, and the beam transitions into the evanescent regime, characterized by oscillations in the transverse direction and exponential decay along the propagation axis [19]. Although numerous studies have investigated the properties and applications of paraxial and nonparaxial Bessel beams [27,28,29,30,31], the relatively unexplored combination of evanescent Bessel beams with $\mathcal{PT}$-symmetry as a scientific research is missing.

For this reason, the main aim of this work is to construct scattering matrices for both propagating and evanescent nonparaxial Bessel beams, thereby providing a quantitative investigation of the $\mathcal{PT}$-symmetry with this class of fields.

Accordingly, in Section 2 we introduce the system under study and the method for constructing the scattering matrices of a Bessel beam; in Section 3 we discuss the obtained results for the two types of scattering matrices as functions of the tunable system parameters; in Section 4 we sum up the obtained results.

2. Methods

Basing on the matrix method proposed in [28], here we will establish the form of the scattering matrix of a Bessel beam interacting with a planar layered system. To this end, we consider the multilayer nonmagnetic structure shown in Figure 1 placed in an air with the dielectric permittivity ${\epsilon }_0=1$. The PT-symmetric structure have to consist of alternating layers of the same thickness $d_G=d_L=d$ with gain and loss described by the complex dielectric permittivity ${\epsilon }_G={\epsilon }^{\prime }-i{\epsilon }^{″}$ and ${\epsilon }_L={\epsilon }^{\prime }+i{\epsilon }^{″}$, respectively. A couple of such slabs forms a unit cell, the number of unit cells N in the system being a variable. Let us consider a monochromatic Bessel beam of the order m propagating along the z-axis incident onto the multilayer system. According to [28], its electric and magnetic field strengths in a medium of permittivity $\epsilon$ read as

$$\begin{array}{c}\mathbf{E}=e^{im\phi +i\beta z}\left(J_m\left(kqr\right)c_2{\mathbf{e}}_z-{\displaystyle \frac{1}{q}}{c}_1({\mathbf{e}}_z\times \mathbf{b})+{\displaystyle \frac{\beta }{kq}}{c}_2\mathbf{b}\right),\\ \mathbf{H}=e^{im\phi +i\beta z}\left(J_m\left(kqr\right)c_1{\mathbf{e}}_z+{\displaystyle \frac{\beta }{kq}}{c}_1\mathbf{b}+{\displaystyle \frac{\epsilon }{q}}{c}_2({\mathbf{e}}_z\times \mathbf{b})\right),\end{array}$$

(1)

where z is the beam propagation direction, $(r,\phi ,z)$ are the cylindrical coordinates, $c_1$ and $c_2$ are complex constants, $J_m\left(x\right)$ is the m-th order Bessel function of the first kind, $k=2\pi /\lambda$ is the wavenumber in vacuum, $\lambda$ is the wavelength, $q=k_t/k=sin\theta$ is the normalized transverse wavenumber, $\theta$ is the half-cone angle of the Bessel beam shown in Figure 1, $\beta =k\sqrt{\epsilon -q^2}$ is the longitudinal wave number, and

$$\mathbf{b}=i{J}_m^{\prime }\left(kqr\right){\mathbf{e}}_r-{\displaystyle \frac{m}{kqr}}{J}_m\left(kqr\right){\mathbf{e}}_{\phi },J_m^{\prime }\left(kqr\right)={\displaystyle \frac{d{J}_m}{d\left(kqr\right)}}.$$

Further we consider the case of a TE-polarized Bessel beam ($c_2=0$). The total field inside the multilayer structure is formed by the superposition of forward- and backward- propagating waves, characterized by complex amplitudes. Thus, the transverse components (with respect to the z axis) of the electric and magnetic fields can then be derived from Equation 1 as

$$\begin{array}{c}{\mathbf{E}}_{\perp }=-{\displaystyle \frac{1}{q}}{e}^{im\phi }\left(c_1^{\prime }{e}^{i\beta z}+c_1^{″}{e}^{-i\beta z}\right)({\mathbf{e}}_z\times \mathbf{b}),\\ {\mathbf{H}}_{\perp }={\displaystyle \frac{\beta }{kq}}{e}^{im\phi }\left(c_1^{\prime }{e}^{i\beta z}-c_1^{″}{e}^{-i\beta z}\right)\mathbf{b},\end{array}$$

where $c_1^{\prime }$ and $c_1^{″}$ are the amplitudes of the forward and backward propagating Bessel beams, respectively. The magnitudes of transverse electric ${\mathbf{E}}_{\perp }=E_{\perp }{\mathbf{e}}_z\times \mathbf{b}$ and magnetic ${\mathbf{H}}_{\perp }=H_{\perp }\mathbf{b}$ fields can be expressed in terms of amplitudes $c_1^{\prime }$ and $c_1^{″}$ as

$$\left(\begin{array}{c}{E}_{\perp }\left(z\right)\\ H_{\perp }\left(z\right)\end{array}\right)=P\left(z\right)\left(\begin{array}{c}{c}_1^{\prime }\\ c_1^{″}\end{array}\right),P\left(z\right)=e^{im\varphi }\left(\begin{array}{cc}-{\displaystyle \frac{1}{q}}{e}^{i\beta z}& -{\displaystyle \frac{1}{q}}{e}^{-i\beta z}\\ {\displaystyle \frac{\beta }{kq}}{e}^{i\beta z}& -{\displaystyle \frac{\beta }{kq}}{e}^{-i\beta z}\end{array}\right).$$

(2)

Following its definition, the scattering matrix relates the amplitudes of outgoing waves to those of incoming ones. Therefore, our objective is to relate the tangential fields at the starting point $z=0$ with those at the final point $z=2Nd$ for any number of unit cells N and thickness of the single cell $2d$ using the matrix representation given in Equations 2. The transverse field components at two boundaries of the multilayer are related by the transfer matrix ${\mathsf{\Theta }}_0^z$ as follows [28]

$$\left(\begin{array}{c}{E}_{\perp }\left(z\right)\\ H_{\perp }\left(z\right)\end{array}\right)={\mathsf{\Theta }}_0^z\left(\begin{array}{c}{E}_{\perp }\left(0\right)\\ H_{\perp }\left(0\right)\end{array}\right),{\mathsf{\Theta }}_0^z=P\left(z\right)P^{-1}\left(0\right)=\left(\begin{array}{cc}cos\beta z& -i{\displaystyle \frac{k}{\beta }}sin\beta z\\ -i{\displaystyle \frac{\beta }{k}}sin\beta z& cos\beta z\end{array}\right).$$

In the considered structure, the Bessel beam passes through the gain and loss layers within a single unit cell. Consequently, the transfer matrix of the unit cell is a product of the transfer matrices of each layer of the same thickness:

$${\left[{\mathsf{\Theta }}_0^d\right]}_L{\left[{\mathsf{\Theta }}_0^d\right]}_G=\left(\begin{array}{cc}cos{\beta }_{L}d& -i{\displaystyle \frac{k}{{\beta }_L}}sin{\beta }_{L}d\\ -i{\displaystyle \frac{{\beta }_L}{k}}sin{\beta }_{L}d& cos{\beta }_{L}d\end{array}\right)\left(\begin{array}{cc}cos{\beta }_{G}d& -i{\displaystyle \frac{k}{{\beta }_G}}sin{\beta }_{G}d\\ -i{\displaystyle \frac{{\beta }_G}{k}}sin{\beta }_{G}d& cos{\beta }_{G}d\end{array}\right),$$

(3)

where the subscripts $i=G,L$ denote the gain or loss media and the longitudinal wave number is given by ${\beta }_i=k\sqrt{{\epsilon }_i-q^2}$. The total transfer matrix of N cells is obtained by raising the matrix described by Equation 3 to the N-th power:

$$\left(\begin{array}{c}{E}_{\perp }\left(2Nd\right)\\ H_{\perp }\left(2Nd\right)\end{array}\right)={\mathsf{\Theta }}_{tot}\left(\begin{array}{c}{E}_{\perp }\left(0\right)\\ H_{\perp }\left(0\right)\end{array}\right)={\left({\left[{\mathsf{\Theta }}_0^d\right]}_L{\left[{\mathsf{\Theta }}_0^d\right]}_G\right)}^N\left(\begin{array}{c}{E}_{\perp }\left(0\right)\\ H_{\perp }\left(0\right)\end{array}\right).$$

Figure 1 shows the amplitudes of the Bessel beams on the left and right sides of the multilayer structure, propagating along the z axis $(A$ and $D)$ and in the opposite direction $(B$ and $C)$. Expressing the field components in terms of the introduced amplitudes ${(D,C)}^T$ at $z=2Nd$ and ${(A,B)}^T$ at $z=0$ using Equation 2, we arrive at

$$\left(\begin{array}{c}{E}_{\perp }\left(2Nd\right)\\ H_{\perp }\left(2Nd\right)\end{array}\right)=P_0\left(2Nd\right)\left(\begin{array}{c}D\\ C\end{array}\right),\left(\begin{array}{c}{E}_{\perp }\left(0\right)\\ H_{\perp }\left(0\right)\end{array}\right)=P_0\left(0\right)\left(\begin{array}{c}A\\ B\end{array}\right),$$

where the subscript 0 in the matrices $P_0\left(2Nd\right)$ and $P_0\left(0\right)$ denotes the ambient medium and superscript T stands for the transpose. As a result, the amplitudes on the right side of the multilayer system are related to those on the left side as

$$\left(\begin{array}{c}D\\ C\end{array}\right)=M\left(\begin{array}{c}A\\ B\end{array}\right),M=P_0^{-1}\left(2Nd\right){\mathsf{\Theta }}_{tot}{P}_0\left(0\right).$$

(4)

The scattering matrix relates the outgoing amplitudes to the incoming ones as [15]

$$\left(\begin{array}{c}B\\ D\end{array}\right)=S_1\left(\begin{array}{c}A\\ C\end{array}\right),S_1=\left(\begin{array}{cc}{r}_L& t\\ t& r_R\end{array}\right),$$

(5)

where $r_L=r_{LR}$ and $t=t_{LR}$ ($r_R=r_{RL}$ and $t=t_{RL}$) are the complex field reflection and transmission coefficients for the wave incident from the left (right) to the right (left). Using Equations 4, we determine the amplitudes B and D by means of A and C and get the following expression for the scattering matrix $S_1$:

$$S_1={\displaystyle \frac{1}{M_{22}}}\left(\begin{array}{cc}-M_{21}& M_{11}{M}_{22}-M_{12}{M}_{21}\\ 1& M_{12}\end{array}\right).$$

Condition $detM=1$ [15] allows us to rewrite the matrix $S_1$ in a simpler form

$$S_1={\displaystyle \frac{1}{M_{22}}}\left(\begin{array}{cc}-M_{21}& 1\\ 1& M_{12}\end{array}\right).$$

Scattering matrix $S_2$ relates the switched outcomming amplitudes according to

$$\left(\begin{array}{c}D\\ B\end{array}\right)=S_2\left(\begin{array}{c}A\\ C\end{array}\right),S_2=\left(\begin{array}{cc}t& r_R\\ r_L& t\end{array}\right)$$

(6)

and can be expressed in terms of the matrix $S_1$ as

$$S_2=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)S_1={\displaystyle \frac{1}{M_{22}}}\left(\begin{array}{cc}1& M_{12}\\ -M_{21}& 1\end{array}\right).$$

3. Results

First of all, we are interested in studying eigenvalues of the scattering matrix for both propagating and evanescent Bessel beams to determine the state of the system in the context of available $\mathcal{PT}$ symmetry. The real part of the dielectric permittivity of layers in the unit cell is taken as ${\epsilon }^{\prime }=2$ and the radiation wavelength is set to $\lambda =1000$ nm. The dependence on the following key parameters affecting $\mathcal{PT}$ symmetry is studied: (i) the imaginary part of the dielectric permittivity ${\epsilon }^{″}$ and (ii) the value of the normalized transverse wavevector component q that determines the wave propagation regime (propagating Bessel beam for $q<1$ and evanescent beam for $q>1$). The parameter space of these two quantities is used to determine the regions of $\mathcal{PT}$-symmetric and $\mathcal{PT}$ symmetry broken states and their evolution when the number of unit cells N changes. It should also be noted that the scattering matrices given by Equations 5 and 6 are independent of the integer azimuthal number m of the Bessel beam. In particular, the m-dependent factors of $P_0^{-1}\left(2Nd\right)$ and $P_0\left(0\right)$ matrices cancel out in Equation 4, leaving the $\mathcal{PT}$-state of the system unaffected.

In Figure 2, we show the behavior of the absolute difference in the moduli of $S_1$-matrix eigenvalues $|s_{1,2}|$ while varying the aforementioned parameters q and ${\epsilon }^{″}$ for the fixed number of unit cells $N=1$. In the density plots, the colored regions exhibit non-zero values of the absolute difference and define the $\mathcal{PT}$ symmetry broken regions. The blue lines are the EP lines separating the $\mathcal{PT}$-symmetric and $\mathcal{PT}$ symmetry broken phases. The black regions attached to the EP lines are $\mathcal{PT}$-symmetric with equal moduli of the eigenvalues $|s_1|=|s_2|$. For the layer thicknesses $d=75$ nm in Figure 2a and $d=125$ nm in Figure 2b, the EP lines between the colored and black regions are absent at $q<1$. This indicates, on the one hand, that the system is in the $\mathcal{PT}$ symmetry broken state for any ${\epsilon }^{″}$, and, on the other hand, that the layer thickness d has no clear effect on the emergence of $\mathcal{PT}$-symmetric states in the region of propagating beams. The increase of the layer thickness results in broadening the region of the great values of the difference $\left|\right|s_1|-|s_2\left|\right|$ in the range from ${\epsilon }^{″}\approx 2.5$ to ${\epsilon }^{″}\approx 3.7$, thus, narrowing the black-colored domain. This means that the black-colored region is not characterized by equal absolute eigenvalues, but rather demonstrates the gradual decrease of $\left|\right|s_1|-|s_2\left|\right|$. For $q>1$, there exist $\mathcal{PT}$-symmetric regions in which the eigenvalue moduli are equal as $|s_1|=|s_2|$, but not equal to the unity. Increasing the layer thickness d leads to reduction of the $\mathcal{PT}$ symmetry broken region and, consequently, allows one to avoid a precise tuning of parameters ${\epsilon }^{″}$ and $q>1$ to observe $\mathcal{PT}$-symmetric states for evanescent Bessel beams experimentally. It should also be noted that the white region in both plots in the range from $q=1$ to $q\approx 1.17$, associated with the excitation of an eigenmode upon transition to the evanescent regime, corresponds to the growth of the only eigenvalue.

The change in the absolute difference $\left|\right|s_1|-|s_2\left|\right|$ when varying the number of unit cells N at fixed $d=125$ nm is shown in Figure 3. The plot for the number of cells $N=1$ depicted in panel (a) has been already discussed above as Figure 2b showing no EPs in the case of propagating Bessel beams. Upon the transition to a couple of unit cells $N=2$, we notice the emergence of $\mathcal{PT}$-symmetric and $\mathcal{PT}$ symmetry broken phases sharply detached one from another for both propagating and evanescent waves. We can conclude that only starting from $N=2$ it is possible to identify a $\mathcal{PT}$-symmetric phase for propagating beams in the case of the scattering matrix $S_1$. In general, $\mathcal{PT}$ symmetry broken ($\mathcal{PT}$-symmetric) phase is observed for the great non-Hermiticity parameters ${\epsilon }^{″}$ and propagating (evanescent) Bessel beams independently of N. The number of $\mathcal{PT}$ symmetric regions M changes at every even number of unit cells $N_e=2,4,6,8$ for propagating and evanescent Bessel beams according to the relationship $M=N_e/2$. For the odd number of unit cells and evanescent Bessel beams, it always exist a region of $\mathcal{PT}$-symmetric phase at small ${\epsilon }^{″}$ and q near 1.

So far, we have considered the eigenvalues of the scattering matrix $S_1$. However, based on the fact that there is a certain freedom in defining the scattering matrix, we also investigate the eigenvalues of the matrix $S_2$, defined by Equation 6.

Figure 4a shows the domains of $\mathcal{PT}$-symmetric and $\mathcal{PT}$ symmetry broken states of the scattering matrix $S_2$ for the layer thickness $d=75$ nm and number of unit cells $N=1$. It is obvious that, in contrast to the eigenvalues of the matrix $S_1$ depicted in Figure 2a at the same system parameters, there exists well-defined $\mathcal{PT}$-symmetric phase with the unimodular eigenvalues $|s_1|=|s_2|=1$ for the propagating Bessel beams, while the region of evanescent beams, on the contrary, corresponds to the $\mathcal{PT}$ symmetry broken phase. As it can be noticed in the case of the $S_1$ matrix, the layer thickness variation only slightly shifts the boundaries of the parameters at which one regime gives the way to another. That is why we will look closer at the response caused by the increasing number of unit cells. First noticeable changes occur at $N=3$. As it is shown in Figure 4b, for $q<1$, an additional line of eigenvalues’ degeneracy points (the green line number 3) appears in the parameter space of ${\epsilon }^{″}$ and q, where no $\mathcal{PT}$ symmetry phase transition takes place. That is, the line 3 is composed of diabolic points. At the same time, the evanescent beams remain completely in the $\mathcal{PT}$ symmetry broken phase. With the further increase in the number of unit cells, the new degeneracy lines (green lines) emerge, and the previous lines shift. It is important to note that this rule applies only to the eigenvalue degeneracy lines, but not to the blue EP lines that do not change their position in the considered parameter space. This tendency is shown for $N=5$ in Figure 4c, where the eigenvalues’ degeneracy line number 3 present in Figure 4b extends into the evanescent region along with the appearance of another line 4 in the propagating region. Despite the appearance of a black region at $q>1$, which may look like $\mathcal{PT}$-symmetric one, the evanescent region still remains in the $\mathcal{PT}$ symmetry broken phase. This assertion is confirmed by the absence of blue lines of the $\mathcal{PT}$ symmetry phase transition. The absolute difference in the eigenvalue moduli tends to zero in this region, but it is not equal to zero. Additionally, it should be noted that the matrix $S_2$ does not exhibit the behavior inherent to the matrix $S_1$ related to the qualitative distinctions between even and odd number of cells.

To clarify the eigenvalue behavior in black regions of Figure 4c, we examine the eigenvalue moduli of the scattering matrix $S_2$ for the same layer thickness $d=75$ nm and number of unit cells $N=5$, but different non-Hermiticity parameters ${\epsilon }^{″}$ as shown in Figure 5.

Figure 5a depicts the absolute eigenvalues for ${\epsilon }^{″}=2$ and propagating Bessel beams ($q<1$). In this case, the system stays in the $\mathcal{PT}$ symmetry broken phase up to the EP at $q\approx 0.780$, beyond which it undergoes a phase transition. In Figure 4c, the EP is recognized by the abrupt color change when transition through the blue EP line occurs. The point $q=0.365$ is a point of eigenvalue degeneracy, where eigenvalues exchange their behaviors. However, this point does not alter the phase, since eigenvectors are not degenerate (diabolic point). A similar situation is observed for evanescent beams ($q>1$), where, as mentioned earlier, the eigenvalues in the $\mathcal{PT}$ symmetry broken phase are unequal. One more diabolic point that does not vary the phase can be noticed at $q\approx 1.215$. With further increase of the normalized transverse wave-vector component q, the eigenvalue moduli tend toward each other but cannot become equal at finite parameters q. Thus, for evanescent Bessel beams, the system keeps the $\mathcal{PT}$ symmetry broken phase despite the presence of the eigenvalue degeneracy point. In Figure 5b, we investigate the behavior of the eigenvalue moduli $|s_{1,2}|$ for the greater value of the parameter ${\epsilon }^{″}=7.5$. In the region of propagating Bessel beams, the point $q\approx 0.261$ is an EP, ensuring the phase transition from the $\mathcal{PT}$-broken to symmetric state. The eigenvalue moduli for evanescent beams, as in the case of the lower non-Hermiticity parameter in Figure 5a, asymptotically tend to merge, indicating a $\mathcal{PT}$ symmetry broken state of the system, which cannot be clearly grasped from Figure 4c.

Appearance of the phase transition between $\mathcal{PT}$-symmetric and $\mathcal{PT}$ symmetry broken phases can be better understood from Figure 6 showing the logarithms of the transmittance $lnT=ln{\left|t\right|}^2$ in Figure 6a, a’, reflectance to the right $ln{R}_R=ln|r_R{|}^2$ in Figure 6b, b’, and reflectance to the left $ln{R}_L=ln|r_L{|}^2$ in Figure 6c, c’ for the layer thickness $d=75$ nm.

Despite the fact that the transmission t and reflection $r_R$, $r_L$ coefficients define both $S_1$ and $S_2$ matrices, they are able to predict the $\mathcal{PT}$-symmetric phase transitions following exclusively from the definition of matrix $S_2$. The dark lines in Figure 6 represent the zero value of the corresponding quantity. The only dark lines at $q=1$ in the plots of $lnT$ [panels (a) and (a’)] are related to the transition from the propagating to evanescent regime, which is also distinguishable in Figure 4a for the $S_2$ matrix. Figure 6a, b, c show the logarithms of the transmittance and reflectances for the number of cells $N=1$. The only bright line in all three plots for the transmittance and reflectances [panels (a), (b) and (c)] defines the excitation of an eigenmode at $q>1$. The horizontal lines ${\epsilon }^{″}=2$ and ${\epsilon }^{″}=7.5$ [panels (b), (c), (b’), and (c’)] correspond to Figure 5a and b, respectively. Vertical blue lines in Figure 5 mark the values of q at which the reflections are minimal. Thus, the dark lines 1 and 2 in Figure 6b and c are the EP lines in the parameter space of ${\epsilon }^{″}$ and q. They can be also noticed in Figure 4a. In the case of the number of cells $N=5$ [panels (a’), (b’) and (c’) in Figure 6], the dark lines numbered 1 and 2 are again the EP lines describing the $\mathcal{PT}$ symmetry phase transition. However, the lines 3 and 4 emerging simultaneously in two plots are the lines of eigenvalues’ degeneracy that are not accompanied by the phase transition. Thus, the vanishing of the single reflectance, $R_L=0$ or $R_R=0$, is the criterion of the $\mathcal{PT}$ symmetry phase transition for the $S_2$ matrix, whereas when $R_L=R_R=0$, the eigenvalues of the scattering matrix degenerate, but the eigenvectors do not, that is, an EP is not formed.

4. Conclusion

We have derived expressions for two types of nonparaxial Bessel beam’s scattering matrices to investigate $\mathcal{PT}$ symmetry phase transitions in multilayer structures composed of alternating gain and loss layers of equal thickness. To comprehensively explore these systems, we have varied key parameters including the imaginary part of layers’ permittivity ${\epsilon }^{″}$, layer thickness d, number of unit cells N, and normalized transverse wavenumber $q=sin\theta$, which is determined by the half-cone angle of the propagating Bessel beam $\theta$ (in the case of evanescent beams, the wavenumber $q>1$).

We have demonstrated that varying the layer thickness does not qualitatively alter the $\mathcal{PT}$ symmetry of the system and can eliminate the need for fine-tuning parameters to observe $\mathcal{PT}$-symmetric states. In contrast, additional layers in the system dramatically change its behavior both for propagating and evanescent Bessel beams. For the scattering matrix $S_1$, this leads to the emergence of new regions with opposite $\mathcal{PT}$-symmetric phases in the parameter space of the non-Hermiticity ${\epsilon }^{″}$ and wavenumber q, the distinct behaviors being observed for even and odd numbers of layers. For the matrix $S_2$, however, increasing the number of layers results in the appearance of additional eigenvalues’ degeneracy (diabolic) lines with weak distinction between even and odd numbers of unit cells. We have conducted a comparative analysis of the transmittance T and reflectances $R_L$ and $R_R$ and associated reflectances’ minima with the EP and eigenvalues’ degeneracy lines of the scattering matrix $S_2$. This identifies a $\mathcal{PT}$ symmetry phase transition when a minimum in the single reflection coefficient is observed and a diabolic line when both left and right reflectances vanish. The results of this work may pave a way for future $\mathcal{PT}$ symmetry research exploiting propagating and evanescent optical Bessel beams.