Spatiotemporally controllable honeycomb superlattice plasma photonic crystals in dielectric barrier discharge

We present the experimental realization of tunable honeycomb superlattice plasma photonic crystals (PPCs) in dielectric barrier discharge by utilizing mesh-liquid electrodes. Fast reconfiguration among the simple honeycomb lattice, honeycomb superlattice, and honeycomb-snowflake superlattice is achieved. A dynamic control on the sizes of center scattering elements in the honeycomb superlattice has been realized. A phenomenological activator-inhibitor reaction diffusion model is established to demonstrate the formation and reconstruction of the honeycomb superlattice. The simulations reproduce well the experimental observations. The photonic band diagrams of different honeycomb PPCs are studied by using the finite element method. The addition of large center elements in honeycomb superlattice yields remarkable omnidirectional band gaps that are about 2.5 times larger than in the simple honeycomb lattice. We propose an effective scheme to fabricate spatiotemporally controllable honeycomb lattices that enable great improvement in band gap size and dynamic control of microwave radiations for wide applications.


Introduction
Photonic crystals (PCs) are metamaterials that offer an unprecedented opportunity for light manipulation and applications in optical sensing and communication [1][2][3][4]. Since the unique features of PCs rely heavily on the band structures, the exploration of PCs with enhanced band gaps is of great interest for both fundamental and applied research. The honeycomb lattice is a promising candidate for the fabrication of a PC lattice, which possesses large band gaps that extend over a wide range of filling fractions [5,6]. In particular, when another rod with different diameter or geometry is introduced into the center of each honeycomb cell, a honeycomb superlattice can be created, giving rise to even larger band gaps [7]. The sizes and shapes of these center elements have a significant effect on the widths and positions of the band gaps. Crystal symmetry reduction by adding disparate scattering elements in a superlattice lifts degeneracy in the photonic band diagrams and opens new ways for engineering photonic gaps [7,8]. Moreover, honeycomb superlattice brings about new properties for PCs. For instance, valley PCs with a honeycomb superlattice structure support robust and valley-dependent topological states [9,10]. They are friendly to micro-nano fabrication and beneficial for designing compact on-chip optical devices [11]. So far, extensive efforts have been made in honeycomb superlattice PCs [9][10][11]; nevertheless, it remains a challenge to dynamically tune the sizes and shapes of the scattering elements. The possibility of tuning the geometric configurations of scattering elements would invest PCs with considerable flexibility, and could lead to new opportunities for photonic devices.
Recently, there has appeared great interest in plasma photonic crystals (PPCs), which are periodic lattices alternating plasmas and dielectric materials [12,13]. The integration of plasmas leads to advantageous properties for PPCs, such as dynamical tunability [14], negative refraction [15], and rapid reconfiguration [16]. Applications of PPCs in plasma stealth aircraft, microwave filter, multichannel communication, and highly sensitive sensors have been proposed [17][18][19][20][21][22][23]. Dielectric barrier discharge (DBD) is one of the most significant systems to realize PPCs. A variety of plasma lattices such as square, hexagonal, and annular lattices have been obtained in DBD by using very simple device [16,[24][25][26]. The high nonlinearity of the system facilitates the rich diversity of the lattices [27]. Besides, rapid reconfiguration between different PPCs is allowed in DBD by readily modulating the lattice constant, electron density, or structural configurations [16,[24][25][26][27]. Only a few seconds are required for the reconfiguration by virtue of the short ionization and recombination time scales in plasmas. To date, great progress has been made in fabricating and measuring PPCs, which find wide applications in the manipulation of microwave or terahertz waves. Nevertheless, most of previous PPCs are exclusively confined to the primitive lattices comprised of identical elements with the same sizes and shapes. It is highly required to achieve tunable superlattice PPCs with improved band gaps. Particularly, it is interesting to realize an independent control on the sizes or shapes of different scattering elements in plasma superlattices, which may offer new properties and capabilities.
In this paper, we present the realization of tunable honeycomb superlattice PPCs in DBD. Rapid reconfiguration among simple honeycomb lattice, honeycomb superlattice, and honeycomb-snowflake superlattice is realized. Significantly larger band gaps have been obtained in the honeycomb superlattice with the addition of a large center element in each unit cell. Active control on both the sizes and shapes of the center elements in honeycomb superlattice has been realized, allowing for real-time control on the widths and positions of the band gaps. A two-component phenomenological reaction diffusion model is established to demonstrate the formation mechanism of the honeycomb superlattice. The numerical simulations agree well with the experimental observations.

Experimental device
The schematic of our DBD system with mesh-liquid electrodes is illustrated in figure 1. A 10 cm × 10 cm metal mesh comprised of hexagonal cells is covered by a quartz glass layer with thickness of 1.5 mm. The side length of each hexagonal unit cell L = 4 mm. The lattice constant defined as the distance between the centers of two neighboring hexagonal cells a = 6.35 mm. The mesh electrode is connected to an alternating-current (AC) power source with driving frequency f = 50 kHz. It introduces a two-dimensional Laplacian electric field to provide a constrained symmetry and lattice constant of the plasma lattices, which is key for yielding controllable and stable PPCs. If the structure of the mesh electrode has changed, the symmetry of plasma structures will be varied correspondingly [24]. The left liquid electrode is fabricated by filling water into a cylindrical vessel with diameter D = 7.5 cm. The water also serves as a coolant and transparent medium for the observation and measurement of the plasma filaments. The water electrode is sealed with 2.0 mm thick quartz glass plates. A metallic ring connected to the ground is placed into the vessel. As the discharge boundary, a square glass frame with thickness of 1.5 mm and inner side-length of 40 mm is clamped between two parallel glass layers. The whole cell is placed in a large vacuum chamber that can be filled with different gas mixtures and pressures. An intensified charge-coupled device (ICCD: Andor DH334T) is used to make temporal resolved measurements for different plasma lattices.
Diagnostics of microwave transmission characteristics of various PPCs are carried out. The microwave is excited at a vector network analyzer (Rohde & Schwarz ZNB40, 100 kHz-40 GHz) and launched from a pyramidal horn antenna (26.5-40 GHz, XB-GH28-20K). It propagates through an aperture with an opening of 4 cm × 8 mm adjacent to the PPC. A pyramidal antenna receiver (26.5-40 GHz, XB-GH28-20K) is placed on the opposite side. The microwave with the transverse magnetic (TM) mode is utilized. The distance between the transmitting and receiving antennas is about 30 cm. The band characteristics can be evaluated from measurements of the transmittance spectra S 21 . More details about the experimental setup can be found in appendix A. Figure 2 shows reconfiguration between the honeycomb lattice, honeycomb superlattice, and honeycomb snowflake superlattice with changes of the discharge parameters. When the applied voltage U = 4.3 kV, the filaments at the six vertices of the hexagonal cells are ignited and arrange into a well-ordered honeycomb structure (figure 2(a)). The increase in applied voltage U = 4.7 kV leads to appearance of the honeycomb superlattice with addition of large plasma elements in the center of each cell ( figure 2(b)). It is a periodic array of scattering elements with different sizes but similar shapes. By changing the air/argon mixture ratio, a novel honeycomb-snowflake superlattice forms at U = 3.1 kV as shown in figure 2(c). Snowflake-like plasma elements appear in the center of each cell and exhibit the same directions. Hence, the honeycomb-snowflake superlattice consists of two types of plasma elements with both different sizes and shapes. With  reconfiguration from honeycomb lattice to honeycomb-snowflake superlattice, the crystal symmetry is reduced. Symmetry breaking enables great improvement in the band gap size as will be demonstrated in figure 12 below. The regulation efficiency between different honeycomb-like lattices is high in our system. These structures are very stable and can be changed from one to another by readily adjusting the applied voltage and gas compositions. One can also control the lattices by regulating the gas pressure, the frequency and waveform of the applied voltage, the distance of gas gap, and the geometry of mesh electrode additionally. The reconfiguration always takes place quickly with a response time of only several seconds. It is reversible with high repeatability. The high stability and quick response to discharge parameters are crucial for fabricating fast-modulated and highly-integrated optical devices.

Experimental observations
By using fast camera diagnostics, we explore the time-resolved measurements of diverse honeycomb lattices. It is shown that these well-defined plasma lattices are actually time averaged structures that possess intriguing spatio-temporal dynamics. For the simple honeycomb lattice as shown in figure 3, only a single current pulse occurs in each half cycle of the applied voltage. This implies that all filaments are ignited simultaneously in a time window ∆t 1 = 150 ns. They exhibit good temporal periodicity, which repeat the same discharge behavior in each half-period with the jittering of current pulses below 30 ns. For honeycomb superlattice, as shown in figure 4, two distinct current pulses have been observed during each half period of  the applied voltage. The simple honeycomb sublattice (H) emerges during the first pulse in a time window about ∆t 1 = 480 ns, while a triangular sublattice (T) composed of large center elements appears during the second pulse with ∆t 2 = 400 ns. Hence, the honeycomb superlattice results from the spatial-temporal integration of the honeycomb substructure and the triangular substructure with the order of H-T-H-T in each period of the applied voltage. Figure 5 illustrates the spatiotemporal dynamics of the honeycomb-snowflake superlattice, which also has two distinct current pulses in each half-cycle of the applied voltage. Similar to the honeycomb superlattice, honeycomb sublattice emerges during the first pulse ∆t 1 , and triangular sublattice composed of large center elements appears during the second pulse ∆t 2 . These large elements are bigger but without irregular shapes, in comparison with their counterparts shown in figure 4(c). Interestingly, inspection of this triangular sublattice reveals that it actually consists of two different discharge stages (insets in figure 5(c)). The triangular sublattice with circle elements emerges in the former stage ∆t 3 , while a network structure (N) emerges in the latter stage ∆t 4 . The network structure is very weak, and even no current growth can be observed in the current waveform. A higher momentary AC voltage is required to reach its breakdown threshold. Consequently, the honeycomb-snowflake superlattice is a superposition of three sublattices: honeycomb sublattice (H), triangular sublattice (T), and network sublattice (N), respectively, with the emerging sequence H-T-N-H-T-N in each period of the applied voltage. All these honeycomb-like structures exhibit good temporal periodicity. The harmonic response with the applied voltage allows one to temporally control the emergence of plasma lattices by modulating the AC voltage frequency.  In comparison with the simple honeycomb lattice and honeycomb-snowflake superlattice, the honeycomb superlattice is a better candidate for the fabrication of a PC structure. Not only can larger band gaps be realized in the honeycomb superlattice [7], it is also easily manufactured from the fabrication standpoint. In view of these advantages, the studies of the honeycomb superlattice are concentrated on in this work. A dynamic control on the sizes of large center elements has been realized as demonstrated in figure 6. When the working gas is 100% ambient air, the radius of large center elements R increases from 0.81 mm to 0.96 mm with a reduction of the applied voltage and gas pressure (figures 6(a)-(c)). When argon has been introduced into the working gas (figure 6(d)), the radius can be further increased to 1.20 mm. By contrast, minor changes take place for the radii of small edge elements, allowing one to control the sizes of large center elements almost independently.
The honeycomb superlattice is robust and can be produced in a wide range of discharge parameters. As shown in figure 7, it is obtained under argon concentration in a range of 0%−60% and gas pressure 76-342 Torr. In contrast, the honeycomb-snowflake superlattice emerges at a high argon concentration. Figure 7(b) shows a detailed study on the radius R of the large center elements in the honeycomb superlattice. R can be adjusted in a range of 0.72-1.45 mm, which scales with the argon concentration, but decreases with increasing of the gas pressure and the voltage. To the best of our knowledge, we present the first realization of controllable sizes of scattering elements in a plasma honeycomb superlattice, while holding the symmetry and lattice constant of the structures being invariant. Since the photonic bands are directly related to the sizes of scattering elements, the tunable sizes will lead to changes in band gap structures to manipulate the light flow of different frequencies. Figure 8 provides a qualitative explanation on the formation mechanism of the honeycomb superlattice. The combined effects of space charges, surface charges and spatially periodic Laplacian fields produced by the mesh electrode is responsible for the generation of a honeycomb superlattice. At t = t 1 , the Laplacian electric fields E app provided by the applied voltage are the largest at six vertices of each hexagonal cell, giving preference to discharges at these locations ( figure 8(b)). Correspondingly, the simple honeycomb sublattice ignites during ∆t 1 . This process leaves surface charges on the dielectric layers, or 'footprints,' in accordance with the array of discharge filaments. These deposited charges will build up a surface charge electric field E s opposite to the current applied voltage, and thus extinguish the discharge eventually. Meanwhile, surface charges also inhibit the discharges in the vicinity of their positions due to spreading along the dielectric layers [28]. With a further increase in momentary AC voltage at t = t 2 , the electric field in the middle of each hexagonal cell, which is far away from the inhibition regions of surface charges, reaches the breakdown threshold and induces large plasma elements in these regions during ∆t 2 . Consequently, after the extinction of these two discharges at t = t 3 , the surface charges are located at the six vertices and the center of each hexagonal cell. As the polarity of the applied voltage is reversed in next half-period, the maximum electric fields occur first again at the six vertices of hexagonal cells, and thus a similar discharge sequence will take place. The plasma lattice tends to spatially coincide with that of the previous half-period because of the 'memory effect' of surface discharges [29].

Discussion on experimental results
The distributions of E s along the surface of dielectric layer covering the cathode at t = t 2 are presented in figure 8(c), corresponding to the four honeycomb superlattices with increased sizes of the center elements as shown in figures 6(a)-(d). At this instant, the simple honeycomb sublattice that occurs during the first current pulse has been extinguished. The total quantity of the transferred electric charges Q can be estimated by integrating the first current pulse Q =´Idt. In principle, it can be assumed that the total amount of surface charges is roughly equivalent to that of the transferred electric charges [30]. This ansatz enables us to use Q to estimate the effects of surface charges qualitatively. The average quantitiesQ of the transferred electric charges for one single pulse areQ 1 = 23.5 nC,Q 2 = 22.2 nC,Q 3 = 17.8 nC, andQ 4 = 15.1 nC, respectively.Q i decreases, leading to the reduction of inhibition regions for surface charges. To be specific, the inhibition effects of surface charges in the center of each honeycomb cell are weakened with the reduction ofQ i as shown in figure 8(c). As a result, the sizes of center elements occurring beyond the inhibition regions are increased.
The formation of tunable sizes of plasma center elements in the honeycomb superlattice can also be explained by the avalanche mechanism. It is known that the visible avalanche radius of the discharge filament r v has the form [31] r where r A is the avalanche radius of the filament, z is the propagation distance in the discharge gap, α denotes the Townsend ionization coefficient, T e is the electron temperature and E is the electric field. It suggests that r v grows proportionally with z and α 1/2 . Since the gap distance z remains invariant in our experiment, α is the main factor that influences r v , which is expressed as Obviously, α increases with the reduction of the gas pressure p, leading to larger sizes of visible avalanches correspondingly. This is consistent with the experimental observations shown in figures 6(a)-(c). On the other hand, the parameter B is 365 V·cm −1 ·Torr −1 in air and 180 V·cm −1 ·Torr −1 in argon. Therefore, the addition of argon leads to a lower value of B and thus to an increased α. Accordantly, the radius of visible avalanche r v is increased with incorporation of argon gas, which agrees well with the experimental result in figure 6(d). Here, we would like to point out that the radii of small edge elements r in the honeycomb superlattice are also increased. Its growing rate is comparable with that of large center elements. Owing to the fact that r is rather small and that the changes of r have little impact on the band structures (appendices C and D), the radius growing of small elements is not considered in this work.

Reaction-diffusion model
In a sense, the generation of honeycomb superlattices in DBDs is a typical self-organization phenomenon in non-equilibrium systems. The reaction-diffusion model has proven to be a crucial way to characterize the phenomenon of self-organization in nonlinear systems, such as biological, chemical, physical systems, etc [32]. Recently, it is shown that the spatial-temporal plasma lattices in gas discharges can also be depicted by a phenomenological reaction-diffusion model [33,34]. The formation of a honeycomb superlattice results from the interplay of space charges, surface charges and Laplacian fields produced by the mesh electrode as described above. Surface charges can be considered to be the inhibitor that suppresses the discharge. The space charges serve as the activator, which activates not only the inhibitors but also of itself through charge multiplication in avalanche processes. The effects of the mesh electrode can be described by a spatially variable parameter. Based on these facts, we explore the modified Gierer-Meinhardt reaction-diffusion model [35] with a spatially variable parameter. In a dimensionless form, the model equations are as follows: where u, υ represent the concentrations of the activator and inhibitor, and D u , D υ denote their corresponding diffusion coefficients. From standard linear theory, it is known that D u must be smaller than D υ for the well-known diffusion-driven instability [36]. ζ and ψ are positive parameters that refer to the transfer rates of the activator and inhibitor, respectively. µ denotes the reaction rate. A spatially variable parameter w is involved in ψ and D u to denote the effects of the spatially periodic potential. They can be expressed as follows: where A is the modulation amplitude and w has the form The spatial periodic distribution described by w is consistent with the electric potential distribution provided by the hexagonal mesh electrode in experiments. Here, k F is the wave number of the hexagonal cells on the mesh electrode with the number of cycles n = 4 and the system size L s = 2. Equation (3) is solved by an explicit Euler method on a square grid containing 100 × 100 points with the time step ∆t = 0.01 and the grid spacing ∆x = ∆y = 0.02. The other parameters are invariable with µ = 10, ζ = 10, ψ 0 = 20, and D υ = 0.27. Figure 9 shows the honeycomb superlattices with increased sizes of the center elements obtained in our model. They are achieved by increasing the modulation amplitude A of the inhibitor. This scenario can be understood as a consequence of the increased Townsend ionization coefficient α with the decrease in gas pressure p or with the addition of argon gas. The neutralization velocity for the surface charges increases with α, leading to a corresponding increase in A. The formation of the honeycomb superlattice with tunable sizes of the center elements can be well described by the Gierer-Meinhardt reaction-diffusion model with spatially periodic modulations. A similar honeycomb superlattice was observed in a Faraday system with the introduction of a third forcing frequency [37]. A rhombus superlattice alternating of large and small spots was also obtained by using the Gierer-Meinhardt reaction-diffusion model with spatially varying parameters [36]. It is believed that the spatiotemporally heterogeneous parameter introduced by the latticed electrode is the key factor in producing such complex superlattices and enhancing the robustness of superlattice formation. The simulation well reproduces the experimental results.

Photonic band diagrams of honeycomb lattices
The dispersion relations of different honeycomb lattices are calculated by utilizing COMSOL software with the Lorentz-Drude model. Photonic band diagrams under TM mode are considered. The propagation of electromagnetic (EM) waves can be described by the Maxwell equations, with the form of Helmholtz equation [38], where the wave vector k 0 can be expressed as The permeability µ r (r) = 1, the dielectric conductivity σ = 0, and the relative dielectric constant ε r (r) can be expressed as ε g in neutral gas (ε g = 1) and ε p in plasma regions, written as [39] where the plasma frequency ω pe = (e 2 n e /ε 0 m) 1/2 , n e is the electron density. Depending on the discharge parameters, n e is generally of the order of 10 13 cm −3 -10 15 cm −3 , which has a significant effect on the band structures. ν m is the electron elastic collision frequency. The boundary of parallelogram primitive cell is chosen as Floquet periodic boundary condition ( figure 10(a)). The normalized eigenfrequencies (ωɑ/2π) for   figure 10(b)). According to Bloch's theorem, each eigenvector E(r) has the form as: where u k (r) is a periodic function of the lattice, and R l is the lattice vector. Figure 11 illustrates the photonic band diagrams of the simple honeycomb lattice and the honeycomb superlattice. One can see that two omnidirectional band gaps (OBGs) and one unidirectional band gap (UBG) in the K-M direction can be obtained for the simple honeycomb lattice ( figure 11(a)). They locate at the positions of 31.9-36.9 GHz (1 st OBG), 54.7-55.1 GHz (2 nd OBG), and 63.8-67.2 GHz, respectively. For the honeycomb superlattice ( figure 11(b)), an OBG forms in the range of 56.3-61.3 GHz, while two UBGs are created in M-Γ, and K-M directions in the ranges of 38.8-43.1 GHz, and 43.0-47.9 GHz, respectively. It is evident that the band diagrams have changed significantly with the reconstruction from the simple honeycomb lattice to the honeycomb superlattice. Figure 11(c) shows an experimental verification of the microwave transmission characteristics for these two honeycomb lattices. The transmission spectra S 21 in the range of 34.2-36.2 GHz is detected, corresponding to the 1 st OBG. One can see that a band gap with the central frequency f 0 = 34.7 GHz can be observed for the simple honeycomb lattice, while no band gap forms for the honeycomb superlattice in this range. The experimental results are consistent with the theoretical calculations. As indicated by the orange arrows in figures 11(a) and (b), there is an OBG for the simple honeycomb lattice, whereas no band gap occurs for the honeycomb superlattice at this position. Owing to the limitation of our experimental apparatus that the stop frequency is 40 GHz, the diagnosis of the higher band gaps for honeycomb superlattices will be performed in the future.
The active control on the sizes of plasma elements in PPCs offers a new route to modulate the photonic band structures. Figure 12 shows the gap map for the honeycomb lattice as the radius of the large element R increases. Here note that, in order to clarify the effects of large elements exclusively, we set the electron density and the radius r to be average values which are invariable in the simulation. It is seen that the addition of large center elements produces a great band gap that is about 2.5 times larger than the simple honeycomb lattice case (R = 0). More importantly, an increase in R leads to great changes of the positions and widths of band gaps. The width of the first OBG gradually reduces with increasing R, and this OBG disappears at R = 0.33 mm. The width of the second OBG decreases to zero at R = 0.5 mm and then increases substantially to ∆f = 12.9 GHz at R = 1.5 mm. These facts suggest that large band gaps can be obtained only in a certain range of filling fractions, which is consistent with the results presented in [7]. Additional discussions on the effects of the radius of small elements r and the electron density n e upon the band structures can be found in appendices D and E, respectively.

Conclusion
We present experimentally and numerically the realization of controllable honeycomb superlattices in DBD with uniquely designed mesh-liquid electrodes. Rapid reconfiguration among simple honeycomb lattice, honeycomb superlattice, and honeycomb-snowflake superlattice is realized. An active control on both the sizes and shapes of the center elements in honeycomb superlattices is achieved. The addition of large center elements in honeycomb superlattices yields a remarkable band gap, which is about 2.5 times larger than in the simple honeycomb lattice case. With changes of the geometric configurations of different honeycomb lattices, the positions and widths of band gaps can be modulated dynamically. Moreover, the spatiotemporally resolved measurements of different honeycomb lattices are carried out by using fast camera diagnostics. It reveals that the simple honeycomb lattice, honeycomb superlattice, and honeycomb snowflake superlattice result from an integration of one, two, and three distinct sublattices, respectively, which can be controlled spatiotemporally. A qualitative explanation for the formation of honeycomb superlattices is provided in terms of the avalanche mechanism, surface charge field, and Laplacian potential of the mesh electrode. The Gierer-Meinhardt reaction-diffusion model with spatially variable modulations is explored to reveal the formation mechanism of the honeycomb superlattice with increased sizes of center elements. The numerical simulations are well consistent with the experimental observations. Our results may provide inspiration for designing new types of controllable metamaterials and find promising applications in the manipulation of microwaves and terahertz waves.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).  Figure 1A shows images of the real experimental device taken from the front and side, respectively. The discharge cell is placed in a large vacuum chamber that can be pumped and filled with a gas mixture. The left water electrode is made by filling tap water into a cylindrical container with diameter of 7.5 cm. The water also serves as a coolant and transparent medium for the observation and measurement of discharge filaments. The right is a mesh electrode that is fabricated by covering a quartz glass layer on a metal mesh composed of hexagonal units. The metal mesh is connected to a sinusoidal AC power supply. It introduces a two-dimensional Laplacian electric field to provide a constrained symmetry and lattice constant of the plasma lattices, which is key for yielding controllable and stable plasma lattices. If the structure of the mesh electrode has changed, the symmetry of plasma structures will be varied correspondingly. A square glass spacer is placed between two parallel electrodes, serving as a lateral boundary. When the applied voltage is increased to the discharge threshold, the filaments will ignite between the gas gap and construct into different structures as the discharge parameters change. Two pyramidal antennas are located on two opposite sides of the vacuum chamber for emitting and receiving the microwaves.  Figure 2 shows the images of plasma lattices taken by a regular digital camera (Canon EOS 6D). This is the same as the appearance seen by our naked eyes. In these pictures, the color of plasma units is lilac because the working gas is 100% air or a mixture of 80% argon and 20% air. The emission spectra of nitrogen molecules and Ar I range from 330-410 nm and 690-820 nm, respectively, resulting in the lilac color of the discharges. Figure 4 shows the time-resolved honeycomb superlattice obtained by using an intensified charge-coupled device camera (ICCD: Andor DH334T). The raw images are actually in black and white mode. Figure 2A shows the raw images of honeycomb sublattice and triangular sublattice, corresponding to the first and second current pulse in the current waveform, respectively. To demonstrate these structures more clearly and furthermore to distinguish the different discharge stages of two sublattices, the images are dyed in different colors, in accordance with the stripe color marking the current pluses. As shown in figure 4(a), the first current pulse is marked by a red stripe. Since the honeycomb sublattice forms during this current pulse, we dye the plasma units of honeycomb sublattice red correspondingly. Similarly, the second current pulse is indicated by a green strip. The triangular sublattice appears during this current peak, and we thus dye the plasma units of triangular sublattice in green color. In this manner, the discharge behavior of different sublattices can be clearly illustrated.

Appendix C. The changing law of r and R in the honeycomb superlattice
Both the sizes of small elements and large elements in the honeycomb superlattice change as the discharge parameters are varied. Figure 3A shows a detailed study on the changes of r and R with an increase of gas pressure under different gas compositions. It can see that r and R can be adjusted in the range of 0.33-0.90 mm and 0.72-1.45 mm, respectively. To further demonstrate the relative changes of two radii, we define the ratio k r = r/r min and k R = R/R min, where r min and R min denote the minimum radii of small and large elements. It can be seen that k r is in the range of 1-2.67, while k R is 1-2.01. To be specific, the growing rates of r and R are comparable. Nevertheless, r is rather small, the changes of r can hardly be recognized, which have little impact on the band structures (appendix D).  Figure 4A shows the photonic band diagrams of the two honeycomb superlattices presented in figures 6(b) and (d), respectively. Here, we keep the radii of large center elements R fixed, while changing the radii of small elements. One is set as the average value as used in the manuscript, and the other is set as the real value. To clarify the effects of the radii of small elements, we set average electron densityn e = 2.3 × 10 14 cm −3 in the simulation. As shown in figures 4A(b) and (c), for R = 0.88 mm, the position and width of OBG varies little when r is changed from 0.47 mm (the real value) to 0.50 mm (the average value). Similarly, for R = 1.2 mm as shown in figures 4A(e) and (f), the position and width of OBG also changes slightly when r is increased from 0.50 mm (the average value) to 0.59 mm (the real value). Consequently, the small changes of r have little influence on the band structures. Remarkably, compared to the effect of r, the radius R has a greater impact on the band structures. As shown in figures 4A(c) and (e), significant changes of band structures can be observed when r = 0.5 mm is fixed, while the radii of the large elements R have changed from 0.88 mm to 1.2 mm. Based on above results, the growing of the radius of small elements is not considered in our manuscript. For simplicity, an average radiusr = 0.5 mm is adopted to exclusively clarify the effects induced by the addition of large elements in figure 12.