A New Lattice Boltzmann Scheme for Photonic Bandgap and Defect Mode Simulation in One-Dimensional Plasma Photonic Crystals

: A novel lattice Boltzmann method (LBM) with a pseudo-equilibrium potential is proposed for electromagnetic wave propagation in one-dimensional (1D) plasma photonic crystals. The ﬁnal form of the LBM incorporates the dispersive effect of plasma media with a pseudo-equilibrium potential in the equilibrium distribution functions. The consistency between the proposed lattice Boltzmann scheme and Maxwell’s equations was rigorously proven based on the Chapman–Enskog expansion technique. Based on the proposed LBM scheme, we investigated the effects of the thickness and relative dielectric constant of a defect layer on the EM wave propagation and defect modes of 1D plasma photonic crystals. We have illustrated that several defect modes can be tuned to appear within the photonic bandgaps. Both the frequency and number of the defect modes could be tuned by changing the relative dielectric constant and thickness of the defect modes. These strategies would assist in the design of narrowband ﬁlters.


Introduction
Photonic crystals (PhCs) are artificial, periodically modulated structures of dielectric materials that can control the propagation of electromagnetic (EM) waves with a range of forbidden frequencies called photonic bandgaps (PBGs). The first PhC devices were fabricated using an array of alumina lattice elements [1]. Since then, PhCs have been extensively studied for various applications [2][3][4][5]. By introducing a defect mode (DM) inside a PBG, the combination of EM waves and PBGs opens a promising perspective to control spontaneous emission [6]. Spontaneous emission is a fundamental characteristic that limits device performance in various fields, including photonics [7], illumination [8], and imaging [9]. It also hinders the realization of large-scale optical integrated circuits [8][9][10]. However, owing to the fixed physical structure of conventional PhCs, PBGs are not actively tunable, which limits their applicability.
In recent years, plasma media have been widely used as photonic crystals [11][12][13][14] due to the tunability of the permittivity via various parameters, such as the plasma frequency and the plasma density, as well as the collision frequency at both the micrometer scale using the microplasma technique [15][16][17][18] and the millimeter scale using the low-temperature plasma technique [19]. Plasma PhCs provide additional freedom to reconfigure the DM and PBGs through customization of the plasma parameters. Nevertheless, their complex design and manufacturing process limit the potential applications of topological photonics in twodimensional (2D) and three-dimensional (3D) PhCs. Consequently, 1D PhCs are preferred due to their flexibility and ease of manufacturing [20,21]. One-dimensional plasma PhCs are employed to realize surface impedance and the bulk band. Furthermore, 1D plasma PhCs have been employed to manipulate light-matter interactions by tuning the plasma layer thickness in PhCs.
The unique EM properties of 1D plasma media have attracted much attention from researchers, leading to the investigation of the EM responses of plasma media using various numerical simulation methods [22][23][24][25][26]. Several groups have proposed various numerical methods to model the EM waves in plasma, such as the finite-difference time-domain (FDTD) method [27,28], the auxiliary FDTD method [29,30], and the recursive convolution FDTD method [31]. As a general numerical simulation method, the lattice Boltzmann method (LBM) has achieved considerable success in various applications in the last two decades, since it divides all simulators into two decoupling parts: streaming and colliding parts [32]. The former, streaming parts, can easily be implemented with a memory swap function, while the latter, colliding parts, are totally reliant on the data of the local cell without using the data from the adjacent cells; which enables them to exhibit excellent scalability and easy implementation [33], as well as to be widely used in solving a vast range of partial differential equations, including fluid flow [34], wave [35], quantum mechanics [36][37][38], and heat transfer [39][40][41][42][43] problems. Due to the LBM's excellent performance in previous areas, there is considerable interest in investigating the applicability and capability of such a method in solving EM waves in 1D plasma PhCs. Several groups have developed relevant codes to solve the EM propagation problem in nondispersive media based on the LBM [44][45][46]. Chen et al. [25] implemented an LBM scheme for dispersive media, but with pseudo-permittivity and complex force terms. To the best of our knowledge, few LBM schemes have been applied to investigate the EM propagation behaviors in 1D plasma PhCs.
This study introduces a new LBM scheme with a pseudo-equilibrium potential to solve Maxwell's equation in isotropic plasma media. Two validation simulation cases were carried out via comparison with existing analytical and FDTD solutions. Based on the validated method, we studied a 1D PhC to investigate the propagation behaviors of its EM waves and the properties of its PBGs.

Methods
The LBM has achieved considerable success in various applications as a general numerical simulation method [34][35][36][37][38][39][40][41]. In property scheme designing, the evolution of an interested physical field can be recovered from the particle distribution of the LBM [33][34][35]. In this study, we propose an LBM scheme with a new form term design, and prove that the new proposed scheme is mathematically consistent with Maxwell's equations in plasma media using the Chapman-Enskog technique.
Maxwell's equations in dielectric materials at a particular location, x, and time, t, are given by [25]: where B and H are the magnetic induction and the magnetic field intensity, respectively; and D and E are the electric displacement and the electric field intensity, respectively, where D(ω) = ε(ω)E(ω), with ε(ω) denoting the dielectric constant of a linear, isotropic, and dispersive medium. We also have B = µH, where µ. denotes the permittivity. Using the relative constants µ r and ε r , the relations µ = µ r µ 0 and ε = ε r ε 0 exist. Since, in the time domain, D(ω) = ε(ω)E(ω) can be rewritten as D( (1) and (2)) in one-dimensional dispersive media could be rewritten in the following form: To solve Equations (3) and (4) and to recover the dispersive phenomena, we propose using the well-known LBM-BGK model, which is expressed as per [31]: where e i and f i (x, t) are the lattice vector and the particle distribution function in the i-th direction, respectively. When i = 0, e 0 is the zero vector [31]. ω = 1 τ , where τ and ∆t are the relaxation time and the lattice time step with b = 2 for one-dimensional media, respectively. Additionally, we further propose a new group of definitions for the macroscopic electromagnetic quantities E y (x, t) and H z (x, t), in addition to the equilibrium distribution f eq i (x, t): where ψ(x, t) is the new proposed pseudo-equilibrium potential, which is used to represent the dispersive effects; A 0 , A, and B are constants. In the following section, we will prove that the LBM model of Equation (5) can recover the Maxwell's equations of Equations (3) and (4) using our proposed forms of Equations (6) and (8), and that the proposed pseudo-equilibrium potential ψ(x, t) could be used to recover the dispersive phenomena. Two extra conservation constraints are introduced to specify the distribution weights: Based on Equations (6) and (8) and Equations (9) and (10), we find: To calculate the pseudo-equilibrium potential in plasma media, we take the Drude model in the time domain: where ϑ c is the damping frequency, ω p is the plasma frequency, and U(t) is the unit step function. We use the Chapman-Enskog expansion to construct the distribution weights and the pseudo-equilibrium potential, through which the evolution features of the EM waves controlled by Equations (3) and (4) can be recovered.
The particle distribution, f i (x, t), can be expanded up to the second order with respect to the expansion parameter, θ: where f 1 i (x, t) and f 2 i (x, t) are the formal expansion functions representing the distribution function at the first and second orders, respectively. Combining Equations (9), (10) and (13): Considering that e 0 is a zero vector, and that θ is an infinitesimal value, Equations (14) and (15) imply: ∑ i e i f k By using the Taylor expansion and combining it with Equation (5), we obtain: The operators ∂ t and ∂ xα can then be expanded by using the Chapman-Enskog expansion: and we can rewrite the pseudo-equilibrium potential in Equations (6) and (8): where ψ(x, t) is assumed to be: Grouping terms with the same order of θ leads to: and With Equations (6) and (8), and summing Equation (24) over i, we obtain: Multiplying both sides of Equation (25) by e iα : Summing both sides of Equation (24) over i gives: Multiplying Equation (24) by e iα : Through comparison with the Maxwell's equations of Equations (3) and (4), we note that the following constraints exist: Combining Equations (25)-(31) gives: From Equations (32) and (33), it can be verified that, with the new pseudo-equilibrium potential, the Maxwell's equations can be recast as ∆t, and that θ can go to zero. Considering Equation (21), the pseudo-equilibrium potential in Equation (7) could be written as: Compared with the FDTD method, the LBM scheme proposed in this study does not require the parameter values from adjacent cells when we need to calculate the equilibrium distribution and pseudo-equilibrium potential, so it can be easily parallelized, which is also one of the most significant advantages of the LBM widely used in other applications [32][33][34][35]. Since we focussed on a 1D EM wave simulation problem, the proposed method cannot be directly used to solve 2D or 3D problems. However, we could have derived the pseudoequilibrium potential for 2D and 3D dispersion materials following the Chapman-Enskog expansion procedure, Equations (13)- (34). Another aspect of this study is that one must Photonics 2022, 9, 464 6 of 15 convert between the LBM and physical spaces. Using the speed of light as a critical unit conversion parameter, Table 1 summarizes the conversion rules between the LBM quantities and their corresponding physical values. Table 1. Conversion between the LBM and physical spaces.

Denomination LBM Context Physical Context
Space step c is the speed of light, DL is the domain length, L is the cell number, and θ V = 1 kV/m.

Results and Discussion
To validate the accuracy of the proposed method, numerical and analytical solutions for three typical cases were considered. The results show that the proposed LBM method was able to correctly recover electromagnetic propagation behaviors in nondispersive and plasma media.

Electromagnetic Waves in Nondispersive Media
For the first validation case, we investigated an electromagnetic Gaussian pulse propagated through a 1D nondispersive medium with 10 4 cells and a dielectric interface at the center. The length of the entire simulation domain is 2.3 mm. The left half of the simulation domain is a vacuum with a relative dielectric constant of ε 0 r = 1.0, while the right half is occupied with a nondispersive dielectric medium with a relative dielectric constant of ε r = 2.0. The electromagnetic Gaussian pulse could be described as [26,46]: where H 0 M and E 0 M are the amplitudes of the magnetic and electric fields, respectively. The constant, α = 400 (in units of lattice cells), fixes the pulse width, and x cen is in the quarter.
The theoretical amplitudes of the reflected, E M,R , and transmitted, E M,T , pulses can be calculated based on the incident pulse, E 0 M , and the relative dielectric constants, ε 0 r and ε r , as [7,47]: from which we calculate that the theoretical ratio of

PBGs in 1D Plasma PhCs Based on LBM-SEF
We further investigated the PBGs of 1D plasma PhCs to demonstrate the adaptability of the proposed LBM method. The dielectric material, D, and the plasma material, P, were placed in alternating layers to form 1D plasma PhCs, as shown in Figure 2  We investigated an incident Gaussian pulse propagating in the abovementioned 1D plasma PhCs, and the amplitude of the pulse is described as Equation (35) with = 10 (in units of lattice cells). The incident wave propagated from the right air medium, traversed the 1D plasma PhCs, and returned to the air medium. The transmitted and reflected electric fields of the incident Gaussian pulse are distributed as shown in Figure 3.

PBGs in 1D Plasma PhCs Based on LBM-SEF
We further investigated the PBGs of 1D plasma PhCs to demonstrate the adaptability of the proposed LBM method. The dielectric material, D, and the plasma material, P, were placed in alternating layers to form 1D plasma PhCs, as shown in Figure 2

PBGs in 1D Plasma PhCs Based on LBM-SEF
We further investigated the PBGs of 1D plasma PhCs to demonstrate the adaptabil of the proposed LBM method. The dielectric material, D, and the plasma material, P, w placed in alternating layers to form 1D plasma PhCs, as shown in Figure 2  We investigated an incident Gaussian pulse propagating in the abovementioned plasma PhCs, and the amplitude of the pulse is described as Equation (35) with = (in units of lattice cells). The incident wave propagated from the right air medium, tra ersed the 1D plasma PhCs, and returned to the air medium. The transmitted and reflect electric fields of the incident Gaussian pulse are distributed as shown in Figure 3. We investigated an incident Gaussian pulse propagating in the abovementioned 1D plasma PhCs, and the amplitude of the pulse is described as Equation (35) with α = 10 (in units of lattice cells). The incident wave propagated from the right air medium, traversed the 1D plasma PhCs, and returned to the air medium. The transmitted and reflected electric fields of the incident Gaussian pulse are distributed as shown in Figure 3.  The simulation results for the 1D plasma PhCs are presented in Figure 3. Figure 3a shows that, at time t = 0 fs and t = 15.3 fs, the Gaussian pulse is in the air zone. Additional after passing through four dielectric layers and four plasma layers from the left, a signi cant proportion of the incoming EM waves is reflected back into the air domain, while t remaining proportion penetrates the plasma. At t = 30.6 and 45.9 fs, Figure 3c,d shows th most of the EM waves propagate through the leftmost plasma layer and still move towa the right. Additionally, there are still tens of small-amplitude pulses propagating in t 1D plasma PhCs. A comparison of the results obtained by our proposed method and t FDTD method is shown as Figure 3a-d, which demonstrates that the electric files obtain by our LBM model agree well with those calculated using the FDTD method. Figure 4 shows the transmission curves of 1D plasma PhCs. There are clearly tw PBGs in the frequency domain region of 2.0-10.0 THz. In this case, the calculated tran mission coefficients, obtained both with the present model and the FDTD method, a plotted as functions of frequency, Figure 4. The two curves are again in good agreeme This is because the PBG frequency domain lies above /2 , and the relative permittiv of plasma in the PBGs region is near that of a vacuum, according to the Drude model. T theoretical most significant bandgap of 1D plasma PhCs is around × = 4.6 THz, where d is the layer thickness of the 1D plasma PhCs, c is light speed, and k is positive integer number. Figure 4 shows that the first gap at 3.8-4.8 THz and the seco gap between 8.5 THz and 9.5 THz are in very good agreement with the theoretical mod The simulation results for the 1D plasma PhCs are presented in Figure 3. Figure 3a,b shows that, at time t = 0 fs and t = 15.3 fs, the Gaussian pulse is in the air zone. Additionally, after passing through four dielectric layers and four plasma layers from the left, a significant proportion of the incoming EM waves is reflected back into the air domain, while the remaining proportion penetrates the plasma. At t = 30.6 and 45.9 fs, Figure 3c,d shows that most of the EM waves propagate through the leftmost plasma layer and still move toward the right. Additionally, there are still tens of small-amplitude pulses propagating in the 1D plasma PhCs. A comparison of the results obtained by our proposed method and the FDTD method is shown as Figure 3a-d, which demonstrates that the electric files obtained by our LBM model agree well with those calculated using the FDTD method. Figure 4 shows the transmission curves of 1D plasma PhCs. There are clearly two PBGs in the frequency domain region of 2.0-10.0 THz. In this case, the calculated transmission coefficients, obtained both with the present model and the FDTD method, are plotted as functions of frequency, Figure 4. The two curves are again in good agreement. This is because the PBG frequency domain lies above ω p /2π, and the relative permittivity of plasma in the PBGs region is near that of a vacuum, according to the Drude model. The theoretical most significant bandgap of 1D plasma PhCs is around f ≈ k × c 5d = k × 4.6 THz, where d is the layer thickness of the 1D plasma PhCs, c is light speed, and k is a positive integer number. Figure 4 shows that the first gap at 3.8-4.8 THz and the second gap between 8.5 THz and 9.5 THz are in very good agreement with the theoretical model.  To calculate the convergence order of our proposed method, Richardson's method [48] is further used to estimate the exact solution of the EM energy density E as: where Er( ) is the estimation error with an order of n+1, δ = , and L and N are the computational domain length and total cell number, respectively. Figure 5 shows the numerical error as a function of the LBM grid cell number, ranging from 10 3 to 10 5 cells, which implies that relative error decreases with an increasing LBM grid cell number, as δ . . The decrease in error verifies that the present scheme has a second-order convergence.

Effects of the Defect Layer Thickness on DMs and PBGs in 1D Plasma PhCs
To demonstrate the adaptability of the proposed method in the investigation of defect modes, we used the proposed lattice Boltzmann scheme to analyze the EM wave propagation behaviors in 1D plasma PhCs with one defect layer, and this presented an interesting phenomenon displayed by the plasma medium in response to incoming EM waves. By tuning the defect layer thickness, we can change the EM wave propagation behavior and the defect modes of 1D plasma PhCs.
The schematic diagram in Figure 6 shows the proposed topological PhCs. One-dimensional PhCs comprised alternating layers of plasma, P, and dielectric, D, with one To calculate the convergence order of our proposed method, Richardson's method [48] is further used to estimate the exact solution of the EM energy density E as: where Er δx n+1 is the estimation error with an order of n + 1, δx = L N , and L and N are the computational domain length and total cell number, respectively. Figure 5 shows the numerical error as a function of the LBM grid cell number, ranging from 10 3 to 10 5 cells, which implies that relative error decreases with an increasing LBM grid cell number, as δx 1.92 . The decrease in error verifies that the present scheme has a secondorder convergence. To calculate the convergence order of our proposed method, Richardson's m [48] is further used to estimate the exact solution of the EM energy density E as: where Er( ) is the estimation error with an order of n+1, δ = , and L and the computational domain length and total cell number, respectively. Figure 5 sho numerical error as a function of the LBM grid cell number, ranging from 10 3 to 10 which implies that relative error decreases with an increasing LBM grid cell num δ . . The decrease in error verifies that the present scheme has a second-order c gence.

Effects of the Defect Layer Thickness on DMs and PBGs in 1D Plasma PhCs
To demonstrate the adaptability of the proposed method in the investigation fect modes, we used the proposed lattice Boltzmann scheme to analyze the EM wave agation behaviors in 1D plasma PhCs with one defect layer, and this presented an esting phenomenon displayed by the plasma medium in response to incoming EM w By tuning the defect layer thickness, we can change the EM wave propagation be and the defect modes of 1D plasma PhCs.
The schematic diagram in Figure 6 shows the proposed topological PhCs. O mensional PhCs comprised alternating layers of plasma, P, and dielectric, D, wi

Effects of the Defect Layer Thickness on DMs and PBGs in 1D Plasma PhCs
To demonstrate the adaptability of the proposed method in the investigation of defect modes, we used the proposed lattice Boltzmann scheme to analyze the EM wave propagation behaviors in 1D plasma PhCs with one defect layer, and this presented an interesting phenomenon displayed by the plasma medium in response to incoming EM waves. By tuning the defect layer thickness, we can change the EM wave propagation behavior and the defect modes of 1D plasma PhCs.
The schematic diagram in Figure 6 shows the proposed topological PhCs. Onedimensional PhCs comprised alternating layers of plasma, P, and dielectric, D, with one embedded defect layer of thickness, DL = 2d (Figure 7) and DL = 4d (Figure 8), in the region bounded by x = 0.28L + 8d and x = 0.28L + 10d, as well as by x = 0.28L + 6d and x = 0.28L + 10d (d = 1.13 µm in this case). We focussed on the effect of the DL on the DM and PBGs in 1D plasma PhCs, where the permeability of the defect layer is set to be the same as that of the dielectric layer D.
Photonics 2022, 9, x FOR PEER REVIEW 10 of 15 region bounded by = 0.28 + 8 and = 0.28 + 10 , as well as by = 0.28 + 6 and = 0.28 + 10 ( = 1.13 μm in this case). We focussed on the effect of the DL on the DM and PBGs in 1D plasma PhCs, where the permeability of the defect layer is set to be the same as that of the dielectric layer D.   region bounded by = 0.28 + 8 and = 0.28 + 10 , as well as by = 0.28 + 6 and = 0.28 + 10 ( = 1.13 μm in this case). We focussed on the effect of the DL on the DM and PBGs in 1D plasma PhCs, where the permeability of the defect layer is set to be the same as that of the dielectric layer D.   shows that the rightward-propagating EM waves are slightly slower for a thicker defect layer, and new reflected signals are observed around the left side of the 1D plasma PhCs as the DL increases from 2d to 4d. Figure 9 shows the DMs of 1D PhCs with different DLs. As the DL increases, new DMs emerge from the bands below and above the bandgaps, which are consistent with the previous research conclusions [49,50]. This figure indicates that the number of DMs can be controlled by varying the DL. Thus, a range of defect thicknesses is suitable for designing multichannel filters.  Figure 8b,c this shows that the rightward-propagating EM waves are slightly slower for a thicker defect layer, and new reflected signals are observed around the left side of the 1D plasma PhCs as the DL increases from 2d to 4d. Figure 9 shows the DMs of 1D PhCs with different DLs. As the DL increases, new DMs emerge from the bands below and above the bandgaps, which are consistent with the previous research conclusions [49,50]. This figure indicates that the number of DMs can be controlled by varying the DL. Thus, a range of defect thicknesses is suitable for designing multichannel filters.  Figures 7b,c and 8b,c this shows that the rightward-propagating EM waves are slightly slower for a thicker defect layer, and new reflected signals are observed around the left side of the 1D plasma PhCs as the DL increases from 2d to 4d. Figure 9 shows the DMs of 1D PhCs with different DLs. As the DL increases, new DMs emerge from the bands below and above the bandgaps, which are consistent with the previous research conclusions [49,50]. This figure indicates that the number of DMs can be controlled by varying the DL. Thus, a range of defect thicknesses is suitable for designing multichannel filters.

Effects of the Relative Dielectric Constant on the DMs and PBGs in 1D Plasma PhCs
Using the same parameters that were used in Figure 6 (except for the relative dielectric constant of the defect layer), the influence of the relative dielectric constant on the electric field and transmission curves is given in Figures 10 and 11, respectively.

Effects of the Relative Dielectric Constant on the DMs and PBGs in 1D Plasma PhCs
Using the same parameters that were used in Figure 6 (except for the relative dielectric constant of the defect layer), the influence of the relative dielectric constant on the electric field and transmission curves is given in Figures 10 and 11, respectively.  Simulation results for different permeabilities are presented in Figure 10a-d. Figure  10a shows that, at t = 0, the initial Gaussian-modulated sinusoidal pulse is in the air zone.

Effects of the Relative Dielectric Constant on the DMs and PBGs in 1D Plasma PhCs
Using the same parameters that were used in Figure 6 (except for the relative tric constant of the defect layer), the influence of the relative dielectric constant electric field and transmission curves is given in Figures 10 and 11, respectively.  Simulation results for different permeabilities are presented in Figure 10a-d. 10a shows that, at t = 0, the initial Gaussian-modulated sinusoidal pulse is in the air Simulation results for different permeabilities are presented in Figure 10a-d. Figure 10a shows that, at t = 0, the initial Gaussian-modulated sinusoidal pulse is in the air zone. It then reaches the right side of the defect layer at approximately t = 15.3 fs, as shown in Figure 10b. Compared with the situation in Figure 7b, the EM waves propagate significantly slower here, owing to the higher permeability of the defect layer in Figure 10b. Figure 10c shows that at t = 30.6 fs most of the EM waves have passed through the 1D plasma PhCs, while there remains a significant reflected signal on the left side of the defect layer compared with Figure 7c. Figure 11 shows the transmission curves for different permeabilities of the defect layer with ε DL = 4.5, 5.0, and 6.0. This shows that, as ε DL varies, the DMs of 1D plasma PhCs are modified accordingly. The defect frequency decreases until it merges with the low edge of the bandgap, and a different DM is generated from the upper edge of the bandgap. The DMs in the third bandgap (8.5-9.5 THz) are more dependent on the permeability of the defect materials than those in the second gap . In addition, a range of defective materials can be used to design multichannel filters.

Conclusions
In summary, a novel LBM method with a pseudo-equilibrium potential was introduced in this study to investigate' the propagating behaviors of EM waves in 1D plasma PhCs. The Drude model was used to describe the dispersive effects of plasma media, and the Chapman-Enskog expansion method was used to prove that the proposed method is mathematically consistent with Maxwell's equations of EM waves in plasma media. The accuracy of the proposed method was demonstrated by comparing it to the analytical and FDTD methods. Two typical cases were considered in order to analyze the characteristics of EM waves that propagate through 1D plasma PhCs involving a defect layer of varying thickness and relative dielectric constants. The results showed that both the number and frequency of DMs are controlled by the thickness and relative dielectric constant of the defect layer. These features not only meet the requirements for designing tunable narrowband filters, but also demonstrate the suitability of the proposed LBM-SEF method for 1D plasma PhCs.

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