Surface potential-adjusted surface states in 3D topological photonic crystals

Surface potential in a topological matter could unprecedentedly localize the waves. However, this surface potential is yet to be exploited in topological photonic systems. Here, we demonstrate that photonic surface states can be induced and controlled by the surface potential in a dielectric double gyroid (DG) photonic crystal. The basis translation in a unit cell enables tuning of the surface potential, which in turn regulates the degree of wave localization. The gradual modulation of DG photonic crystals enables the generation of a pseudomagnetic field. Overall, this study shows the interplay between surface potential and pseudomagnetic field regarding the surface states. The physical consequences outlined herein not only widen the scope of surface states in 3D photonic crystals but also highlight the importance of surface treatments in a photonic system.

The DG photonic crystals are found to exhibit Weyl points due to the geometrical perturbation, which appears as a defect-like shape.To realize a pseudomagnetic field, we constitute a photonic array of the DG unit cells with a spatial gradient of the perturbation, i.e., the degree of the defect varies with the position.We then compute quantized Landau levels and eigenstates along the zeroth Landau level to quantitate the asymmetric localization of photonic waves on the surfaces.This observed asymmetry in wave localization evidences the existence of the surface potential.Then, we tune the translation of the basis of the unit cell to tame the surface termination.This tuning exquisitely controls the surface potential so that the degree of wave localization varies with respect to the tuning.Finally, we implement such results to the evasion behaviors of a photonic wave to observe the interplay between the surface potential and pseudomagnetic field.

Results
The attraction of waves by the surface potential in the Weyl system First, let us describe the pseudomagnetic field's effects using the Weyl Hamiltonian.For this, we consider a system periodic along the x 1 -and x 2 -directions for an orthogonal coordinate system.The system consists of N unit cells along the x 3 -direction (the horizontal rightward direction in Fig. 1a-b), and surface boundaries are parallel to the x 1 -and x 2 -directions.We assume that this system is governed by that describes a Weyl point at k w , where σ 0 is the 2 × 2 identity matrix, and σ i ( i = 1, 2, 3 ) are the Pauli matrices.We decompose the Weyl point's location k w into k w = k w,0 + A n w , as introduced in Refs. 13,37,71,85.k w,0 is a constant while A n w depends on the unit cell's index n .We assume that A n w linearly varies with a specific value p (henceforth, this value is referred to as perturbation strength) which also linearly varies with n: The systems are finite along x 3 -direction (the horizontal direction) and periodic along x 1 -and x 2 -direction.In each panel, the varying p is schematically represented as the size of yellow circles.Both systems are described by Eq. ( 1).(c) Comparisons of their eigenstates' distributions along x 3 -direction, exhibiting the mutual symmetric dispersions.The eigenstates are calculated by Eq. ( 1).(d) Schematic plot of V s , the scalar coefficient of the surface potential.(e) Comparisons of the eigenstates' distributions of (a and b) when V s is applied.Here, the eigenstates are calculated by Eq. ( 4).In (c and e) gray dotted lines are the symmetric curves with respect to the center for comparisons of the red and blue plots.
where p s is a proportional constant (refer to Fig. 1a,b).Then, the pseudomagnetic field can be written as B = ∇ × A n w , which is proportional to p s .Let us suppose two systems with different signs of slopes, i.e., one system described by p s = −p s,0 (Fig. 1a) and the other described by p s = +p s,0 (Fig. 1b) where p s,0 is positive.Their p(n) are schematically illustrated as the size of yellow circles in Fig. 1a,b.The pseudomagnetic fields B for these two are in the same magnitudes and in opposite directions.Their wave localization of the zeroth Landau level at a specific point in the momentum space shows mutually symmetric characteristics for these opposite fields, as shown in Fig. 1c.Now, let us implement a surface potential V as follows: where V s is nonzero only around the boundaries 36,62,[65][66][67] , like Fig. 1d.When we include the surface potential term with V s , further detailed in Methods, Eq. ( 1) is rewritten as: By applying the same V for the two systems in Fig. 1a,b, we can observe mutually asymmetric localization, as shown in Fig. 1e.The two localized states commonly show a slight shift in the left direction, implying that the surface potential attracts them along that direction.
Derivation of a pseudomagnetic field with Eqs.(1) and ( 2) and the relevant data in Fig. 1a,c are not new as they already have been introduced in several studies 13,37,71,85 .The same explanation is applied to the concept of the surface potential 36,62,[65][66][67] .However, to our knowledge, adopting a surface potential into the pseudomagnetic field has not been tried.Thus, the next several sections are discussions about how to interplay them.

Applying surface potential in a finite-sized array
Although there can be several methods to drive a surface potential into a photonic crystal, we herein adjust the surface terminations in a finite-sized array.In a fully periodic crystal, the fields or eigenvectors of a propagating mode are invariant under discrete translation, not depending on a specific position of the basis in a unit cell.On the contrary, a finite-sized array has terminations at the boundaries.A propagating wave profile depends on the terminations according to the variations of the basis's positions in the unit cells 86 , as shown in Fig. 2a.Thus, the effect of surface potential can be quantitated with respect to adjusting the surface terminations.
Here, we remark on the followings: (i) Eqs.(1) to (4) do not consider a detailed geometry in each unit cell.Thus, observation of the relation between the surface termination and surface potential should be phenomenologically carried out using a real array structure.(ii) A photonic band structure for a finite-sized array displays projected bands, called folding of bands.Without the pseudomagnetic field, the zeroth Landau level adheres to the bulk bands so that surface states cannot be obtained.Thus, we should use the pseudomagnetic field, as shown in Fig. 2b.Our study shows the interplay between the surface potential and pseudomagnetic field; the pseudomagnetic field gives the surface localization of a photonic wave, and the surface potential tunes the degree of the localization.
To build the array in Fig. 2b and to induce both the surface potential and pseudomagnetic field, we use a DG photonic crystal that exhibits Weyl points 17 .We adjust the surface terminations by tuning the translation of the basis along the normal direction to the boundaries.

Double gyroid photonic array exhibiting pseudomagnetic field
To drive the surface potential by the surface termination, we use a photonic array that exhibits a pseudomagnetic field by k w = k w,0 + A n w in Eq. ( 4).We consider the DG, as shown in Fig. 3a.The yellow and blue single gyroids (SGs) are given by a set of respectively, where f D 2 and f O are the level-set values that determine the volume fraction of each SG.The yellow SG has selectively the perturbation term pf p (x) to break the inversion symmetry and to generate Weyl points in momentum space.Due to this perturbation, the yellow SG exhibits a defect-like narrow region on the arm passing the a 1 a 2 surface, as shown in Fig. 3a.The higher perturbation strength p induces the deeper defect-like shape; this corresponds to the yellow circle's size in Fig. 1a,b and Fig. 2. (See Methods for the detailed explanations about the DG crystal.)When this DG photonic crystal is periodic along all three lattice vector directions, it can exhibit four Weyl points, as shown in Fig. 3b,c.The Weyl points H 0 and N 0 (marked with blue and red solid points in Fig. 3c, respectively) have positive and negative topological charges, respectively.
Varying p of the DG shifts the positions of the Weyl points, as shown in Fig. 3d.We denote the moved positions of the Weyl points N 0 and H 0 in momentum space by k N w and k H w , respectively.Increasing p makes the Weyl points shift away from Γ-point on the single plane ( (001)-plane).For the narrow range of p around the central value p 0 , the traces of all Weyl points exhibit linear shapes so that we can write their positions as where k H s and k N s are overall directions of the traces (see Fig. 3d).The deviations of shifted Weyl points, δk N and δk H , also show linear relations with p , as shown in Fig. 3e.
To generate the pseudomagnetic field, we design a photonic array made of DGs whose p linearly varies with respect to the position in the array 13,71,72,87,88 .First, we assume a DG-array that consists of DGs with N unit cells  along the a ⊥ -direction between two boundaries, as shown in Fig. 4a.The array is periodic along the a = -and a -directions, but does not preserve the translational symmetry along the a ⊥ -direction.Next, we apply the non- uniform geometry on this array using the perturbation strength p that linearly varies along the a ⊥ : like the inset linear plot in Fig. 4a.Here, x ⊥ is a coordinate along the a ⊥ , d 0 is a distance between the planes at n = 0 , and n = N/2 , and p 0 is a central value of p(x) .Thus, the location of the Weyl point varies with the a ⊥ -directional coordinate.From the information in Fig. 4a,b, the pseudomagnetic field has the form B = ∇ x × A = p s B = a = + B � a � , which are parallel to the array boundaries (see the orange planes in Fig. 4c,d).
As a result, the waves can be localized at the boundaries (see Methods for detailed explanations about the array and pseudomagnetic field calculations).

Surface potential by surface termination and resulting photonic wave localization
Plugging the design in the previous section to the array of N = 48 primitive cells allows us to generate Landau spectra.The detailed explanations of the Landau levels related to the pseudomagnetic field are in Sect. 1, Supplementary Information.
Here, we extract and compare the zeroth Landau levels by two arrays with p s = −2p s,0 and p s = +2p s,0 , as illustrated in Fig. 5a.Although the two arrays have different internal geometries due to the different p s values, their overall translation status are identical, i.e., they generally use the formulae denoted in Fig. 3a.Thus, we consider that their surface terminations are identical.The photonic band structure in Fig. 5a shows that the resulting zeroth Landau levels accessible from the two arrays are not equal.Furthermore, the eigenstate intensity distributions for the two cases reveal mutual asymmetric (Fig. 5b,c).The distributions in Fig. 5d,e are biased toward the n = 48 from the symmetric curve (the gray dotted lines).From the fact that the results in Fig. 5d,e exhibit the same tendency as with Fig. 1e, we can conclude that there exists a surface potential in these arrays.Now, we replace x as x − ha ⊥ to apply overall translation h of the DG by h|a ⊥ | along the a ⊥ -direction to verify the surface termination's effect, as illustrated in Fig. 6a.Note that the profile of p does not move (see the upper right plot in Fig. 6a).For a specific point on the zeroth Landau levels (marked in Fig. 6b), we calculate the field intensity for several h values, as shown in Fig. 6c.The discussions so far are about that the surface potential can , B H ′ , B N ′′ , and B H ′′ , respectively.These are calculated from the traces around N 0 , H 0 , −N 0 , and −H 0 marked in Fig. 3d, respectively.adjust the degree of localization of surface states generated by the pseudomagnetic field.The results in Fig. 6c indicate that the localization degree can be adjusted by the surface termination.The zeroth Landau level in the photonic band structure also can be moved to another position by changing h .Due to the relatively small h , the overall configuration of the zeroth Landau level and the surrounding bulk bands remains almost intact, as shown in Fig. 6b.

Evasion behavior of photonic wave
We embody the evasion behavior of a photonic wave in a DG structure using oppositely graded p .Adjusting the propagation path of a photonic or phononic wave has been an interest in photonics or wave mechanics.These have been proved using materials with negative-refractive indices 16,[89][90][91] or Dirac/Weyl crystals exhibiting oneway propagations of waves 3,12,13,44,92 .Such proofs are related to invisibility-cloaking.Here, we perform the behavior by pulling a photonic wave that was originally propagating on one boundary toward the opposite boundary.
DG arrays that consist of 8 cells with p s = −2p s,0 , 0 , p s,0 , 3p s,0 , and 5p s,0 are prepared, based on the array in Fig. 4a, ( p s,0 > 0 ).number of cells is counted based on the body-centered cubic primitive cell.)The reason for using a smaller number of cells than in the previous case is to broaden intervals between Landau levels.We classify these arrays into groups with p s = −2p s,0 and the others with p s ≥ 0 .The directions of ∇ x p for these two groups are opposite to each other (see the inset of Fig. 7a,b).The band structures for p s < 0 and p s ≥ 0 (plotted in Fig. 7a,b, respectively) exhibit Landau levels along ŴH ′ -direction.From these plots, we consider ωa/2πc = 0.49 as an optimized frequency for the evasion behavior because this frequency between the bulk bands and all zeroth Landau levels commonly meets this frequency only once.We then constitute a DG array, as shown in Fig. 7c, using the 8-cells arrays depicted in the inset of Fig. 7a,b.Each system consists of 8 and 32 cells along a ⊥ -and The DG array reveals the evasion behavior of photonic waves.At ωa/2πc = 0.49 , a photonic wave is localized in the section of p s = −2p s,0 around the boundary at n = 0 .We then input an incident photonic wave around the localization region, as marked by the star symbols in Fig. 7d-g.If the blocks m = [8, 24] are the same as the other blocks, the localization would continue in the central region.We simulate here the situations that the central region's p s is respectively 0 , p s,0 , 3p s,0 , and 5p s,0 .The results show the gradual attraction of the photonic waves onto the boundary at n = 8 with increasing p s (see Fig. 7d-g).
In addition, we observe the attraction of waves towards the n = 0 boundaries.In the blocks m = [0, 8] and m = [24, 32] , the waves are localized on the boundary.On the contrary, the localized waves in the blocks m = [8, 24] exhibit biased behaviors even when the central blocks' p s is 3p s,0 or 5p s,0 .This coincides to the mutual asymmetric distributions of the waves shown in Fig. 5e, and we see here the existence of a surface potential in the blocks.Then, we can conclude that this arises from surface potential and surface termination.If we adjust the surface termination, we may observe the bias towards the opposite direction.All these show the interplay between surface potential and pseudomagnetic field.

Discussion
We have demonstrated the asymmetric localization of photonic waves via an interplay between surface potential and pseudomagnetic field using a DG photonic crystal.The pseudomagnetic field has induced the surface states of the photonic waves and the surface potential adjusted the degree of localization.The pseudomagnetic field was formed by the graded location of Weyl points, and the surface potential was applied using surface termination by tuning unit cells' basis translation.We have observed the shift of the surface states as a result of the surface potential.In fact, the interplay by the surface potential and pseudomagnetic field generates natural and predictable results so that this might be underestimated.However, the various methods of applying a surface potential (even though they have been utilized in metals, semimetals, or acoustics) can be interplayed with a pseudomagnetic field in photonics in the near future.Thus, we expect that our current study will be a starting point of these directions.
In the case of graphene 66 , boron nitride 67 , and semimetals 36,65,68 , surface potentials can be applied using materials' surface/edge passivation.Although this study used surface termination by crystal's basis translation, we believe that there could be several types of photonic passivation, for example, doping with a thin material such as a dielectric sheet, imposing perfectly magnetic conductor boundary condition, cutting position of the crystal, or attaching a band-gap material.Then, this study and the follow-up studies will open the various possibility of using topological photonic crystals.Meanwhile, there is no example of the detailed analysis of the surface states along the zeroth Landau level by the Hall effect and pseudomagnetic field in three-dimensional Weyl photonic crystals.Therefore, this study will fill this gap thereby this study will positively affect the other studies on the three-dimensional quantum (spin) Hall effect [93][94][95][96][97][98] .

Eigenstates localization by Weyl Hamiltonian
We apply Eq. ( 1) or (4) to the finite system illustrated in Fig. 1a-b to observe the wave localization like Fig. 1c or e.Instead of deriving the Landau level analytically using ladder operators 33,99 , we use a traveling wave solution ψ(x, t) = ue i(k•x−ωt) to consider the finite array.The traveling wave solution can be rewritten as the product of a state u n that depends on only n and an exponent e i(k 1 x 1 +k 2 x 2 −ωt) that depends on only other variables: ψ(x, t) = u n e i(k 1 x 1 +k 2 x 2 −ωt) .The surface localization is obtained by substituting this into Eq. (1)or (4).To consider the surface potential in Eq. (3), we used the plot in Fig. 8 as V s .Detailed derivations, explanations, and additional results related to Fig. 1 are given in Sect.2, Supplementary Information.

Preparations of DG photonic crystal
In the following, we give detailed explanations of the DG photonic crystal and array.In this study, we consider the DG, reported in ref. 17 .The body-centered cubic primitive cell of this structure is defined by the lattice vectors a i = a/2[1, 1, 1] − a x i , where [ x 1 x 2 x 3 ] = I .The mathematical formulae of f SG,Y (x) , f p (x) , and f SG,B (x) denoted i n F i g .3 a a r e f SG,Y (x) = sinX 1 cosX 2 + sinX 2 cosX 3 + sinX 3 cosX 1 , f p (x) = sin(X 1 + X 2 ) , a n d f SG,B (x) = sin X 1 cos X 2 + sin X 2 cos X 3 + sin X 3 cos X 1 , respectively.Here, the local coordinates are given by X = [X 1 , X 2 , X 3 ] = (2π/a)(x − a s s) and X = X 1 , X 2 , X 3 = (2π/a)x , where a is a lattice constant, s = 0.0578 is the shift coefficient that describes the translation of the yellow SG, and a s = a 1 + a 2 + 2a 3 is the translation

Simulation details
All photonic structure simulations were performed using COMSOL Multiphysics®.To input periodicity, the Floquet periodic boundary condition was imposed on all periodic boundaries.All band structures were obtained using the 'Eigenfrequency' solver.In each band structure, we plot only nine bands above and below the zeroth Landau level so that only 19 bands are displayed in a band structure.

Figure 1 .
Figure1.Effects of surface potential on the wave localization by pseudomagnetic fields.(a, b) Schematic illustrations of two systems with opposite p s .The systems are finite along x 3 -direction (the horizontal direction) and periodic along x 1 -and x 2 -direction.In each panel, the varying p is schematically represented as the size of yellow circles.Both systems are described by Eq. (1).(c) Comparisons of their eigenstates' distributions along x 3 -direction, exhibiting the mutual symmetric dispersions.The eigenstates are calculated by Eq. (1).(d) Schematic plot of V s , the scalar coefficient of the surface potential.(e) Comparisons of the eigenstates' distributions of (a and b) when V s is applied.Here, the eigenstates are calculated by Eq. (4).In (c and e) gray dotted lines are the symmetric curves with respect to the center for comparisons of the red and blue plots.

Figure 2 .
Figure 2. Schematics on the realization of surface potential.(a) Inducing different surface potentials by different surface terminations.The surface terminations are adjusted by the basis's positions in a unit cell.(b) Adding pseudomagnetic field on the systems in (a).Generating surface localization of a photonic wave is performed by the pseudomagnetic field, and the degree of localization depends on the surface potential.

Figure 3 .
Figure3.Design of 3D photonic crystal for pseudomagnetic field.(a) A DG photonic crystal.The SGs are defined by x that satisfies the inequalities denoted in the figure.The yellow SG's inequality has an additional term pf p (x) to impose a defect-like shape on this structure (see the bottom-right arm). (b) Simulated photonic band structure exhibiting four Weyl points, the band degeneracies.(c) Weyl points marked in the threedimensional first Brillouin zone.The four Weyl points in (b) are marked on the green plane.(d) Simulation result of Weyl points' movements on the plane with the perturbation strength p .Weyl points in (c) are also marked here with red and blue solid points and denoted by N 0 and H 0 , respectively.(e) Deviations of the Weyl points' movements with p from the positions of N 0 and H 0 , shown in the right and left enlargements in (d), respectively.

Figure 4 .
Figure 4. Design of a photonic array and resulting pseudomagnetic fields.(a) Schematic geometry which consists of N cells with spatially linearly varying p along a ⊥ -direction.This gradient is imposed only on the yellow structure.(b) Bulk and surface Brillouin zones.(c, d) Direction vectors of the pseudomagnetic field around each Weyl point by the above design, denoted by B N ′, B H ′ , B N ′′ , and B H ′′ , respectively.These are calculated from the traces around N 0 , H 0 , −N 0 , and −H 0 marked in Fig.3d, respectively.

eFigure 5 .
Figure 5. Photonic wave localization and evidence of surface potential.(a) Comparison of Landau levels by the opposite perturbation fields.The left and right cases correspond to p s = −2p s,0 and p s = 2p s,0 , respectively. (b, c) Normalized eigenstates respectively along the blue-and red-colored Landau levels shown in (a) to see surface states localized on the boundaries.(d, e) Comparisons of the wave localization along the vertical lines in (b andc).Data with the same symbols are overlapped in the same plot.Gray dotted lines are the symmetric curves with respect to n = 24 for comparisons of the red and blue plots.All the plots are simulation results.

Figure 6 .
Figure 6.Adjustment of surface potential by surface termination.(a), Schematic illustration of three surface terminations.Observation of the surface potentials is carried out by adjusting the basis's positions in each unit cell.The upper left three structures are enlargements of the bottom right arrays.The lower left three figures are respectively the schematics of the upper left structures, and they clearly show the translation of the structures along a ⊥ -direction.The right upper plot indicates that the perturbation p versus the position n is the same for all three cases.(b, c) Simulation results: zeroth Landau levels (b) and eigenstate intensities (c) at the point marked in (b) for three different surface terminations tuned by translation h along the a ⊥ -direction when N = 48 and p s = 2p s,0 are used.

Figure 7 .
Figure 7. Evasion behavior of photonic waves using heterogeneous blocks.(a, b) Zeroth Landau levels by the 8-block system with several perturbation fields, where p s,0 = 7.0711 × 10 −3 a −1 .The perturbation fields used in (a, b) are opposite.(c) Heterogeneous DG system.The PEC is applied to the pink and green boundaries.In pink and green colored bounded blocks, the perturbation fields' directions are opposite as marked.The system is periodic only along a -direction.(d-g) Evasion behavior of photonic waves with several p s > 0 values of the central blocks with fixing p s < 0 values of the blocks around both ends.Incident points are marked as the star symbols.

Figure 9 .
Figure 9. DG primitive cells array that was input to get Fig. 7a,b in the main text.The array consists of 8 primitive cells with spatially linearly decreasing or increasing p along a ⊥ -direction.