Tunable Photonic Band Gaps in Two-Dimensional Bravais–Moiré Photonic Crystal Composed of High-Tc Superconductors

: In this study, we perform a theoretical study of light propagation properties in two-dimensional square photonic crystals (PCs) following Bravais–Moiré (BM) patterns composed of copper oxide high-temperature superconductors (HTSCs). The BM PCs are made of cylindrical cores formed from the combination of two square Bravais lattices. The Moiré pattern forms due to a commensurable rotation of one of these lattices with respect to the other. The dielectric function of the superconducting material is modeled by the two-ﬂuid Gorter–Casimir theory. We report on the corresponding gap, the mapping as a function of the radius of dielectric cores, as well as the dispersion relations of TM modes for BM PCs and for the waveguide system built of defect lines within such a crystal. The BM PCs were composed of copper oxide HTSCs, which exhibit large tunability in terms of temperature.


Introduction
Photonic crystals (PCs) are structures whose refractive index varies periodically in space. By creating a contrast in the refractive indices of the materials that compose the structure, it is possible to obtain photonic band gaps (PBGs) for electromagnetic radiation. This implies the absence of transmission for signals with frequencies within a given interval of values. Such a property is highly useful in applications involving optical devices [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. For these applications, it is essential to tune the properties of PBGs and to obtain the widest possible band gaps. Typically, the contrast between refractive indices is achieved by changing the filling fraction (the ratio between the radius of the cylinders and the lattice constant) and the type of lattice. Once the geometry is determined, PBGs can only be adjusted by changing the constituent materials. In this sense, using tunable materials in the design and fabrication of PCs allows one to adjust the optical properties with high flexibility. Tunable materials include metals and semiconductors, which are highly dispersive and whose band structure depends strongly on plasma frequencies [21,22]. However, these types of materials have some drawbacks. For example, metals exhibit inherent losses, which are caused by the extinction coefficient. To overcome these problems, a variety of unconventional materials have been used, such as superconductors [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. In superconductors, optical properties are strongly influenced by temperature and magnetic fields. Compared to PCs with metallic components, superconductors have the advantage of low loss and adjustable permittivity. In the superconducting state, the electromagnetic wave can only propagate within the London penetration depth of the materials. In the normal state, the London penetration depth becomes infinite, and the electromagnetic wave can propagate in the material without limitations [7].

Description of the System
In this work, we focus on the analysis of guided modes in BM networks. We consider the special case presented in ref. [26], where a commensurate rotation of φ = 53.13 • with r = 3, s = 1 (BM-R3S1, according to the notation presented in [26]) is considered. Figure 1a shows the unit cell (UC) corresponding to this structure. In this case, the UC is composed of 40 dielectric cores that have a circular cross section (blue circles in the figure), of which 20 correspond to the square Bravais lattice without rotation (with radius r nr and dielectric permittivity ε nr ), and 20 dielectric cores with circular cross section (green circles in the figure), which correspond to the rotated Bravais lattice (all with radius r r and dielectric permittivity ε r ). This rotated-at a commensurate angle-square lattice was superimposed over the first one to form the BM pattern. In the construction, a is the lattice constant and ε b is the dielectric permittivity of the background. This construction is further detailed in refs. [26,33,34].
To design CROWs, the structure is optimized in such a way that the PBGs are as wide as possible to ensure that the guided modes appearing in the PBG region, after introducing defects, are located in the central part of the PBG and that they do not interact with modes in the continuum [40]. For this analysis, a scan is performed over the filling fraction, and the band structure is calculated in the UC of Figure 1a. The parameters for the UC that optimize the PBG are used to construct the CROW. The CROW considered in this work is shown in Figure 1b. It is constructed by removing some cylinders from the UC. Particularly, in Figure 1d, the cylinders that are removed are shown in red. For the analysis of the dispersion relations in the CROW, it must be considered that the waveguide is not a perfectly periodic structure in both directions, as can be seen in the UC marked in yellow in Figure 1b. However, it is still periodic in the waveguide direction. For the calculation, an infinite translation of a part of the PC containing the defect was assumed. Thus, the UC is larger and contains several periods (see the yellow box in Figure 1b). Since the signal propagates in the waveguide direction, radiation confinement occurs in the transverse direction. Hence, the wave vector is separated into two components: the first one is perpendicular to the direction of wave propagation and is connected to the dispersion relation of the waveguide, and the second one, which is in the waveguide direction, is known as the propagation constant and is designated as β [40]. Figure 1c shows the Brillouin zone of the waveguide. The unit cell utilized for the analysis is represented by the yellow rectangle with orange edges. The dashed rectangle indicates the direction along which the wave propagates inside the guide. β is the propagation constant and k x is the x-component of the wavevector. (c) Shows the Brillouin zone and the k-path in waveguide case. (d) (in red) the cores removed from the unit cell in order to build the waveguide in (b). In our model, the cylindrical active cores are assumed to be made of a high critical temperature Bi 1.85 Pb 0.35 Sr 2 Ca 2 Cu 3.1 O y superconductor compound, and the background is air.

Theoretical Framework
Maxwell equations govern the propagation of electromagnetic waves in the PC. Through their solutions, we shall be able to investigate the photonic properties of the system of interest. For a 2D medium, electric and magnetic fields can be written as linear combinations of transverse electric (TE) and transverse magnetic (TM) components. Here, we restrict ourselves to dealing with TM modes. Accordingly, the sourceless Maxwell-Helmholtz equation for the magnetic field reads: where H z ( r) is the field's z-component at the position r. Time-harmonic electric and magnetic fields come from the solution of Equation (1), according to the following: PC periodicity is used in order to obtain the corresponding dispersion relations. It comes from the application of the Bloch theorem at the borders of the unit cell. When the electric (magnetic) field propagates between the PC points separated by a lattice vector, R, the only effect on the field is the change in phase given by the following: Here, R is a photonic lattice vector and k is the signal wavevector. In addition, the dielectric function of the superconductor material (SC) was modeled by the two-fluid model to describe the electrodynamics of the superconductor cores at nonzero temperatures [2,10]: ε ∞ is the high-frequency permittivity of HTSC and γ is a damping term. ω sp = c/(λ L (T) √ ε ∞ ) and ω np = (n n e 2 /mε 0 ε ∞ ) 1/2 are the plasma frequencies of superconducting and normal conducting electrons, respectively. Here, λ L (T) = (m/µ 0 n s e 2 ) 1/2 is the London penetration depth, and the dependency of it on the temperature from the Gorter-Casimir model can be written for HTSC as per the following [7]: where c, n n , n s , m , ε 0 , µ 0 , e, and T c are the speed of light in vacum, concentration of normal electrons, concentration of superconductor electrons , the electron mass, permittivity,permeability of vacuum, charge of the electron, and the critical temperature of the superconductor, respectively.

Simulation Settings
A solution of the wave Equation (1) was accomplished via the finite element method (FEM). The core of FEM consists of mess generation: the partition of geometry into small units of a given shape. Once constructed, dependent variables are approximated by means of known -shape-functions. By inserting this approximation in Equation (1), it is possible to generate a set of algebraic equations, which are suitably solved. Then, the solution of the differential problem is assembled. To carry out our study, we used FEM as implemented in the COMSOL-MULTIPHYSICS licensed software package 5.6 [41]. The description of this tool, including the foundations of FEM, the building of meshes, the discretization of differential equations, optimization, and the convergence criteria can be found in refs. [42,43].

Results and Discussion
In the design of the 2D PC considered in this work, the cylindrical cores (both rotated and non-rotated) were made of the SC compound Bi 1.85 Pb 0.35 Sr 2 Ca 2 Cu 3.1 O y . This compound has a ε ∞ = 12, λ L (0) = 22.452 µm and Tc = 107 K [7]. Figure 2 shows the calculated PBG mapping for the TM modes in the BM PC, which is schematically depicted in Figure 1. This mapping appears plotted in Figure 2a-for the fixed value of temperature, T = 15 K-as a function of the ratio between the core radius and lattice constant, r/a, which considers that both rotated and unrotated elements have the same radius, r r = r nr = r. Figure 2b contains the mapping as a function of T for the value r/a = 0.0334 = r m , which corresponds to the widest PBG obtained from Figure 2a. In this variation, the widest PBG (∆ω = 0.37) appears just for T/T c = 1. The dispersion relation that corresponds to r m and T = 15 K appears in Figure 2c, with the vertical axis in units of ωa/2πc. The forbidden regions, represented by the blue fringes in Figure 2a,b, and the lowest gray stripe in Figure 2c correspond to the interval of frequencies below the T-dependent cut-off value, ω c a/2πc = 0.67, at which the effective dielectric function of the structure vanishes. Notice that, as expected, ω c becomes zero when T = T c , but when T is fixed at 15 K, it has a slight increasing rate with r/a. This indicates the effect of geometry on the setting of the overall dielectric response in the unit cell. The second gray stripe highlighted in the dispersion relation is, precisely, the widest PBG, resulting from setting r m at 15 K, (∆ω = 0.25), which is located between ωa/2πc = 2.01 and ωa/2πc = 2.26. It is worth noting that, in this and the following figures showing PBG mappings, the lowest colored stripes just indicate the frequency interval below the cut-off value, for which no propagation of any electromagnetic signal in the system is allowed.
According to Figure 2, due to the chosen design, our working frequency range lies between 0 and 2 THz, which falls within the far-infrared region of the spectrum. As pointed out in refs. [5,7], if the collision frequency γ = 1/τ for unpaired electrons is much larger than the working frequency, the contribution from the third term in Equation (5) can be neglected, as has already been conducted in previous works [7]. Based on the results presented in ref. [44], the collision time for Bi-based materials can be estimated to be around τ = 1.38 × 10 −14 s. This leads to the condition ω 72.5 THz, which is fulfilled without problems within the declared working range. Without the influence of such a term in ε sc (ω), the absorption losses from the normal state subsystem are not relevant in our calculation.
To work within the mentioned frequency range, the lattice constant of the structure will then be chosen to operate in this region of the spectrum; thus, we set a = 0.45 mm. In order to achieve a greater contrast in the refractive indices of the materials in the unit cell, the cylindrical cores are considered in an air background, where ε b = 1.0. Furthermore, copper oxide HTSCs exhibit strong two-dimensional anisotropy. In this case, to evaluate the TM modes in the two-dimensional photonic crystals considered here, the electric field is considered parallel to the c-axis of the HTSC ( E c) [7].
In Figure 3a, the results for the gap mapping of TM modes in the considered BM PC, the fixing of unrotated cylinders at r nr /a = 0.0334, and the varying in the size of the unrotated cores are shown. The remaining parameters are kept with the same values employed to produce Figure 2. One now observes the opening of a main gap for lower frequencies, which widens toward the smallest value of r r /a, thus reaching the maximum ∆ω = 0.39 when r r /a = 0.01, as depicted in Figure 3c. It is also noticed that the forbidden fringe below the cut-off frequency evolves differently, with a slightly smaller value for the maximum gap configuration when compared with Figure 2c.
Another result related with PBG properties appears from the analysis of gap mapping presented in Figure 4. This time, the calculation has a fixed value of rotated SC cylinders, r r /a = 0.0334, while the radius of unrotated ones, r nr , changes. As before, the remaining parameters are the same used for Figure 2. The main difference revealed with this particular setup is the appearance of a double PBG for the lowest values of r nr /a. In the case of maximum amplitudes (r nr /a = 0.01), both gaps have almost the same widths, ∆ω 1 = 0.26 and ∆ω 1 = 0.24 (in units of ωa/2πc), and are separated by a narrow interval of ωa/2πc = 0.03.
To summarize the effects of geometry manipulation on the presence and features of PBGs in BM PCs made of high-T c Bi 1.85 Pb 0.35 Sr 2 Ca 2 Cu 3.1 O y SC, Table 1 contains the compilation of the above-commented results. In accordance, the particular setup shown in the middle row (corresponding to the situation depicted in Figure 3) results in the most suitable setup for PBG enhancement.  The gap mapping as a function of temperature, for the geometry at which the widest gap in (a) appears: r r /a = 0.0334 and r nr /a = 0.01. Again, when t = T c (i.e., when the maximum phonic gap is obtained ∆ω = 0.33), then T/T c = 1 is obtained. (c) Dispersion relation for the same TM modes, corresponding to the same filling fraction with r r /a = 0.0334 and r nr /a = 0.01, for T = 15 K. There is a double gap with almost the same amplitudes: ∆ω 1 = 0.26, located between ωa/2πc = 1.83 and ωa/2πc = 2.09 and ∆ω 2 = 0.24, located between ωa/2πc = 1.56 and ωa/2πc = 1.80, as well as a cutoff frequency of ω c a/2πc = 0.55. The bottom stripes correspond to the cut-off frequency. The flatness of dispersion relations for some of these coupled resonator modes indicates at a significant reduction in the group velocity, v g , of the associated light signals. For that reason, we performed a detailed analysis of this quantity for the different guided modes revealed in order to identify the temperature-dependent slow light properties in the system. In this sense, Figures 6 and 7 contain the analysis of coupled resonator waveguide modes, which vary as a result of the increase in the temperature of the structure from T = 15 K to the T c value. These four figures, respectively, correspond to the guide mode cases depicted in Figure 5d,e (with separate analysis of the modes within the lower ( Figure 8) and upper (Figure 9) PBGs) and Figure 5f. In general, it is possible to observe a trend showing a lowering in the frequency position of the guided modes with an increase in T. In addition, some of those modes showing greater dispersion at low values of the temperature become flatter when it approaches the critical value. However, the influence of the geometrical setup on their group velocities should be more carefully analyzed. In fact, the situations in which there is a single main PBG with either one or two localized modes-that is, the ones observed in Figure 5a,c-display rather symmetric variations of |v g |/c with respect to the center value of the interval for propagation constant values. Such cases are depicted in Figures 6b and 7b. Their respective normalized maximum values (NMVs) for such a quantity are presented for the different temperatures considered in Tables 2 and 3. From them, it is possible to observe values of the NMVs of group velocity in the order of 10 −3 , which decrease as long as T augments toward T c ; in particular, the mode labeled as A 1 B 1 in Figure 5f, is where this stands out. This mode is one of those that reduces its dispersion and becomes flatter with augmenting T. It bears the smallest calculated value for the NMV of group velocity: |v g |/c = 0.0013 (almost three orders of magnitude of reduction from c). As such, it is possible to speak here of a real slow light phenomena.    Table 2. Maximum normalized group velocities corresponding to the guided modes analyzed in Figure 6.

T(K)
Mode Figure 7. The same as in Figure 6 but for the guided modes corresponding to the configuration leading to Figure 5c,f. Table 3. Maximum normalized group velocities corresponding to the guided modes analyzed in Figure 7.
However, the picture showing the guided modes appearing under the geometrical setup corresponding to Figure 5b,e is significantly different. They are localized within the two close PBGs of almost the same widths that appear in the configuration with a radius of rotated SC cylinders r r /a = 0.0334 and with radius of the rotated ones r nr /a = 0.01. As already stated, the modes confined within the lower (A 0 B 0 and A 1 B 1 ) and upper (A 2 B 2 and A 3 B 3 ) PBGs are analyzed separately in Figures 8 and 9. Here, the NMVs of |v g |/c appear, and they are reported in Tables 4 and 5, respectively.  By first analyzing the guided modes inside the lower PBG, it is possible to observe their tendency to become largely dispersive for the larger values of the propagation constant, and for them to become very close in frequency. This indicates a possible coupling between them, as noticed in the case of T = 75 K. True, such a dispersive tendency reduces with the increment of T toward the superconductor transition value. This phenomenon is reflected in the group velocity by suppressing any kind of symmetric variation with β, as occurs in the above-discussed cases. In fact, the position of the curve maximum shifts to higher values in the horizontal axis of Figure 8b. From Table 4, one may observe that, once again, increasing temperature leads to a reduction in the NMVs of the group velocity for both modes, whereby the upper one achieved the smaller value. Nonetheless, the achieved values remain in the order of 10 −2 , and these results are not better than those discussed in the two previous cases. Table 4. Maximum normalized group velocities corresponding to the guided modes analyzed in Figure 8.

T(K)
Mode With regard to the second (upper) pair of waveguide modes shown in Figure 5b,e, the plots in Figure 9a indicate the already commented shift to smaller frequencies. These are induced by the increment of temperature, as well as due to the trend of the lower one to become dispersionless with T. In contrast, the higher of these two modes keeps a rather strong dispersion. Interestingly, this particular mode, at the time of evolving with the temperature, is located farther away from the other. As a consequence, it was seen that the associated curve of |v g |/c tends to recover the symmetric character with respect to the center of the β region and the group velocity augments, as observed from Table 5 (mode labeled as A 3 B 3 ). On the other hand, the lower of the two waveguide modes under discussion shows a non-symmetric variation of group velocity, which shifts toward the intermediate and higher values of the propagation constant, at which point the greater dispersion occurs. As we have already said, the two modes become closer in that region, and their coupling seems to affect the group velocity of the signal in such a way that |v g |/c curves turn to asymmetry within the region, with a reduction in the NMVs of signal velocity. The only exception for these is what occurred in the picture at T = 45 K. In this case, the A 3 B 3 mode strongly disperses, approaching A 2 B 2 mode at higher values of β. According to Table 5, the MNVs of group velocity remain higher than the values presented for the other two geometrical setups, although the A 2 B 2 mode evolves within small multiples of 10 −2 . Thus, one may also talk about the situation of slow light for it within the range of temperature involved. Table 5. Maximum normalized group velocities corresponding to the guided modes analyzed in Figure 9. Previous works on Cu-based HTSC 2D PCs, such as the one reported in refs. [7,45], where simple square lattices of SC cylinders are considered, our system presents better properties for adjusting PBG properties, as can be seen in Figures 2-4. This is because, by having multiple cores in the UC, the active component of the SC cores can be changed by the filling fraction with a greater flexibility than in the case of the typical square lattice. In fact, our results show that the band structure is strongly affected when considering geometrical changes, particularly in the radius of the non-rotated cores compared to the rotated ones. Compared with more recent works, such as the one reported in [10], we found that-in addition to the flexibility for adjusting the guided modes as a function of temperature-the guided modes in our system exhibit a low group velocity, which is not observed in the referenced work. As such, it is worth highlighting, at this point, that such low values of normalized group velocities for waveguide signals are obtained with the use of a HTSC as constituent of BM PCs. These are, in some cases, almost an order of magnitude smaller than those calculated for the BM PCs made of dielectric cores [26]. It is also possible to compare with the works that consider CROWs in other designs of lattices with multiple atoms, such as the one reported in [23]. There, the authors consider CROWs in Archimedean-like lattices with dielectric cores, and with report group velocities that are smaller than those reported in usual networks, some even as low as 0.030 times the speed of light in a vacuum. In our work with SC BM PCs, we found velocities as small as 0.0013 times the speed of light in a vacuum, which is three orders of magnitude smaller than the speed of light. This represents a true slow-light behavior, with the additional possibility of adjusting such velocity with temperature, thus indicating the advantage of SC-based photonic structures for slow-light applications.

Conclusions
In this work, we have studied the properties of light propagation in a two-dimensional photonic crystal whose unit cell was constructed according a Bravais-Moiré pattern and was composed of high-T c superconducting cylinders. This study also includes the proposal for a design of a line-defect waveguide.In addition, the coupled resonator modes associated to it were also investigated. It was shown that the different geometrical setups arising from the fixation, or the variation of the core cylinder radii, have an impact on the photonic gap mapping of the crystal and also on the appearance and amount of waveguide modes inside such gaps, without a coupling with the continuum.
In addition, the analysis of the group velocity for the guided modes, when considering the variation in the temperature up to the superconductor transition value, allows one to identify particular configurations for which the values of the normalized maximum value of this property reduces to almost three orders of magnitude compared to the vacuum speed of the electromagnetic waves with those frequencies. This fact indicates positive prospects for slow-light applications. Funding: This research was funded by the Mexican CONACYT (grant number A1-S-8218 (MEMR)). CAD is grateful to the Colombian Agencies: CODI-Universidad de Antioquia (Estrategia de Sostenibilidad de la Universidad de Antioquia) and projects "Propiedades magneto-ópticas y óptica no lineal en superredes de Grafeno", "Estudio de propiedades ópticas en sistemas semiconductores de dimensiones nanoscópicas", "Propiedades de transporte, espintrónicas y térmicas en el sistema molecular ZincPorfirina", and "Complejos excitónicos y propiedades de transporte en sistemas nanométricos de semiconductores con simetría axial", as well as to the Facultad de Ciencias Exactas y Naturales-Universidad de Antioquia (CAD exclusive dedication project 2022-2023).

Data Availability Statement:
Data is partially available through direct contact with the corresponding authors.

Conflicts of Interest:
The authors declare no conflict of interest.