Hidden-symmetry-enforced nexus points of nodal lines in layer-stacked dielectric photonic crystals

It was recently demonstrated that the connectivities of bands emerging from zero frequency in dielectric photonic crystals are distinct from their electronic counterparts with the same space groups. We discover that in an AB-layer-stacked photonic crystal composed of anisotropic dielectrics, the unique photonic band connectivity leads to a new kind of symmetry-enforced triply degenerate points at the nexuses of two nodal rings and a Kramers-like nodal line. The emergence and intersection of the line nodes are guaranteed by a generalized 1/4-period screw rotation symmetry of Maxwell’s equations. The bands with a constant kz and iso-frequency surfaces near a nexus point both disperse as a spin-1 Dirac-like cone, giving rise to exotic transport features of light at the nexus point. We show that spin-1 conical diffraction occurs at the nexus point, which can be used to manipulate the charges of optical vortices. Our work reveals that Maxwell’s equations can have hidden symmetries induced by the fractional periodicity of the material tensor components and hence paves the way to finding novel topological nodal structures unique to photonic systems.


S1. Calculating Band structure using transfer matrix approach
In this section, we derive the analytical expressions of the band structures in the k y = 0 plane for the AB-layer-stacked PhC. For a given angular frequency ω and a wavevector κ = (κ x , 0, k z ) in the k y = 0 plane, there are two plane wave eigensolutions in the homogeneous anisotropic medium described by Eq. (2) in the main text. The two plane wave solutions can be labeled by theirM y parities, and their dispersion relations are given bŷ where k 0 = ω/c, ε 1 ε 3 ≡ ε xx ε zz − ε 2 xz , and we note that ε xz = ±g in layer A and layer B respectively. The corresponding eigenvectors areM M y − even ∶ ϕ even The Bloch eigenfunctions with certainM y parity in the first period of the PhC (x ∈ [−L/2, L/2] ) can be expressed as the superpositions of the plane wave fields in Eqs. (3) and (4): M y − even ∶ ψ even α = c α ϕ even +,α + d α ϕ even −,α .
where α = A, B labels the fields in layer A (x ∈ [−L/2, 0] ) and layer B (x ∈ [0, L/2] ) respectively. And according to the Bloch condition, the field in the m th period is given by ψ α (x + mL) = ψ α (x)e imL kx with k x denoting the x component of the Bloch wavevector inside the 1 st Brillouin zone (BZ).
From the continuity conditions of (E y , H z ) and (E z , H y ) at the intracell interface x = 0, we obtain which shows that wherek ± = κ 2 ± gk z /ε xx . As a result, the transformer matrices forM y -odd and even modes can be written aŝ Solving the two equations, we find that the Bloch wavevector in the x direction takes the simple expression k x = (±κ 1 mod 2π) for odd modes and k x = (±κ 2 mod 2π) for even modes. Therefore, the dispersions ofM y -odd and even bands in the k y = 0 plane readM with m ∈ Z numbering the bands. And the corresponding normalized Bloch states arê Eqs. (11) and (12) exhibit that allM y -odd (even) bands have identical conical dispersions up to a translation along the x-axis in the extended BZ. Therefore, any pair of bands with opposite mirror parities will cross each other along a nodal ling in the k y = 0 plane. We note that the nodal line can be any kinds of conic sections, including ellipse, parabola, and hyperbola, depending on the parameters of the PhC. In particular, the two red nodal rings of our interest in the main text correspond to the intersections of the even band of m even = 0 and the odd bands of m odd = ±1, given by the following implicit equation: And from the Eqs. (11) and (12) , we can also obtain the frequency and wavevectors of the pair of triply degenerate crossing points of the two nodal rings, i.e. the triply degenerate nexus points: ω NP = 2πc L √ εyy−εxx and

S2. Hidden symmetries in the ky = 0 plane
In this section, we discuss the hidden symmetry in the k y = 0 plane induced by the fractional periodicity of the components of constitutive tensors. In this subsystem, the Maxwell's equations in the layer-stacked dielectric PhC can be written as where ε ↔ r takes the form of Eq. (2) in the main text. And hereinafter, we adopt the natural units with ε 0 = µ 0 = 1 for convenience. In general, the Maxwell's equations are invariant under a symmetry transformationÃ, as long as whereĈ can be an arbitrary invertible operator. However, for space group symmetries of the structure,Ĉ is fixed as identity. If the constitutive tensorM is invertible, which is always true for dielectric PhCs, Eq. (17) is equivalent to the invariance of the effective HamiltonianĤ =M −1N , i.e.ÃĤÃ −1 =Ĥ. In the k y = 0 plane, the effective Hamiltonian readsĤ In what follows, we show that the fractional periodicity of the elements of ε ↔ r can give rise to a hidden symmetry of H(k y = 0) beyond space groups.

A. Hidden symmetries in theMy-odd subspace
In theM y -odd subspace, since the electric field is polarized in the y direction, theM y -odd band structure is entirely determined by ε yy . Especially, for the AB-layered PhC in Fig. 1 of the main text, ε yy is a global constant, and hence the band structure of the odd modes on the k y = 0 plane is directly obtained by folding the light cone in a homogeneous medium with a constant permittivity. Consequently, allM y -odd bands are twofold degenerate along Γ − Z except for the one connected to ω = |k| = 0 point. In fact, we can relax the condition of a constant ε yy so that the period of ε yy is a fraction 1/N of the primitive period L of the whole PhC, namely ε yy (x + L/N ) = ε yy (x), the width of the genuine BZ of theM y -odd subspace (marked as BZ(odd)) should be 2N π/L which is N times as large as the primitive BZ of the whole system (marked as BZ(whole)). Therefore, the finial band structure in BZ(whole) is obtained by translating the bands in BZ(odd) periodically with spacing ∆x = 2π/L. As a result, twofold degenerate Kramers-like NLs can appear along Γ − Z as long as N ≥ 3. If we also require the PhC respects the space group R 2 ⋊ Rod(22), N should be an even number and thus the minimal value of N is 4. We introduce the 1/4-period translation operator in theM y -odd subspace: whereT x (L/4) denotes the original translation operator, andP ± = 1 2 (Î ±M y ) is the projection operator ontoM yeven(+)/odd(-) subspace. So the fractional periodicity of ε yy is equivalent toT odd x (L/4)Ĥ(k y = 0)T odd x (L/4) −1 = H(k y = 0). At the same time, theM y -odd subsapce is also invariant under a generalized 1/4-period twofold screw operation about the x axis:S namely, the combination of a 1/4-period translation and a twofold rotation about the x axis acting on the odd subspace. We will prove in the next section that the generalized 1/4-period twofold screw symmetry, together with the time reversal symmetry, guarantees the emergence of the Kramers-like NLs along Γ − Z.

B. Generalized fractional translation symmetry in theMy-even subspace
For theM y -even subspace, it is convenient to deal with the sub-Hamiltonian of Eq. (18) acting on ψ even = Though the period of the sub-Hamiltonian is the same as the full system, we will show that it can be converted into a new form with fractional period via a similarity transformation. First, we change the eigen basis from ψ even = and the sub-Hamiltonian is transformed accordingly: It shows that −i∂ x always appears together with an effective gauge potential A x = k z εxz εxx inĤ ′ even . Therefore, we can εxx(ξ) dξ] to remove the effective gauge potential in the Hamiltonian: We note that the gauge transformation is k z dependent. Under the coordinate z representation, it appears as a translation operator with x-dependent translation along the z axis: After the combined similarity transformation (ÛĜ),Ĥ ′′ even is only explicitly dependent on ε xx and ε 1 ε 3 = ε xx ε zz −(ε xz ) 2 as shown in Eq. (24). For the AB-layer-stacked PhC shown in Fig. 1 of the main text, both ε xx and ε 1 ε 3 are constant. Therefore, theM y -even band structure are also obtained by folding the light cone of a homogeneous system, and all bands are twofold generate along Γ − Z except for the lowest one attached at ω = |k| = 0.
Similar to theM y -odd case, if we relax the condition to that the common period of ε xx (x) and ε 1 (x)ε 3 (x) is a fraction 1/N of the full period L of the PhC with N ≥ 3, twofold degeneracies of even bands can still appear along the Γ − Z line. And without changing the space group of the system, the minimal value of N is 4 and the corresponding components of the permittivity should satisfy In the following, we will focus on this special case with minimized constraints on the PhC that are compatible with the space group R 2 ⋊ Rod(22) and support the Kramers-like NLs along Γ − Z.
The L/4 periodicity ofĤ ′′ even can be regarded as a generalized fractional translation symmetry operating on thê M y -even subspace for the original Hamiltonian (18),T even whereĜ is reformulated asĜ Moreover, sinceĤ ′′ even also respects the twofold rotation symmetryĈ 2x about the x axis, we can further introduce the generalized 1/4-period twofold screw operator: and the generalized 1/4-period twofold screw operator in the k y = 0 plane: whereŨ It can be directly checked that the effective Hamiltonian given in Eq. (18) is invariant underT x (L/4) andS L/4 , providing that the permittivity of the PhC satisfies Eq. (26), Consequently, we have demonstrated that the fractional periodicity of the elements of permittivity tensor engenders the hidden symmetries of the Maxwell's equations. And we have the relation between the generalized 1/4-period screw rotation and the generalized fractional translation: and the similar result is established for the generalized screw rotation: with s = ±1, ±i denoting theS L/4 branch index of the Bloch state. We note that theS L/4 branch index is only well defined for the states on the k x -axis. Nevertheless, sinceS 2 L/4 =T x (L/2), we havẽ Therefore, theS 2 L/4 -parity (equal to the square of the branch index s 2 = ±1) is well defined for the whole band on the k y = 0 plane, and is determined by the branch index of the states Ψ (s) (k x , 0, 0) on that band. For convenience, we will assign a "pseudo branch index" for every Bloch state on the k y = 0 plane as Ψ (s) (k x , 0, k z ), while we remind that only its square s 2 denoting theS 2 L/4 -parity of the state is meaningful in general (the sign of the pseudo branch index for a state with k z ≠ 0 is indeterminate, as s 2 = (−s) 2 ), but the sign of s makes practical sense for the states on the k x -axis.

S3. Kramers degeneracies induced by the generalized 1/4-period twofold screw symmetry
Here, we express the time reversal operator in the coordinate-momentum mixed representation (x, k y , k z ): Here, (Ũ †Ũ )-pseudo-antiunitarity means that where ψ, ϕ represent two arbitrary states.
Θ L/4 operating on a Bloch state Ψ (s) (k x , 0, k z ) on k y = 0 plane yields a new Bloch stateΨ(−k x , 0, k z ) =Θ L/4 Ψ (s) (k x , 0, k z ) of the same frequency at (−k x , 0, k z ). In addition, sincẽ the new Bloch stateΨ(−k x , 0, k z ) =Ψ (s * ) (−k x , 0, k z ) has pseudo branch index s * and has the sameS 2 L/4 -parity (s * ) 2 = s 2 as Ψ (s) (k x , 0, k z ). In the mean time, In particular, on the Γ − Z line (k where the second equality of the first line is due to the pseudo-antiunitarity ofΘ L/4 . SinceŨ †Ũ is positive definite, On the other hand, it is well known that the ordinary twofold screw symmetryŜ 2x together with the time reversal symmetry T protects all Bloch states on the k x = π/L plane are doubly degenerate. And the two degenerate states Ψ(π/L, k y , k z ) andΨ(π/L, k y , k z ) are correlated by the combined antiunitary operatorΘ L/2 = TŜ 2x : In particular, at X point (k = (π/L, 0, 0)), if Ψ (s) (π/L, 0, 0) has branch index s, we havẽ where the second equality of the first line is due to the fact thatS L/4 andΘ L/2 are commutable. This result shows thatΨ(π/L, 0, 0) =Ψ (−is * ) (π/L, 0, 0) has branch index −is * . In other words, the pair of Bloch bands (along the k x -axis) degenerate at X must have branch indices of either s = 1, s = −i or s = −1, s = i. S4. Determining the branch indices of the two bands connected to |k| = ω = 0 In this section, we figure out the branch indices of the two bands attached at |k| = ω = 0. The expansion of the Bloch states Ψ(k x , 0, 0) at k x = ω = 0 reads with u 0 = (e 0 , h 0 ) ⊺ . Substitution of the expansion into the Maxwell's equation (16) yieldŝ where d 0 x , b 0 x denote the x components of the 0-order D and B fields respectively. When k z = 0,Û is reduced to the identity matrix, thusŨ (k Therefore, both the two bands stemming from |k| = ω = 0 have the same branch index s = −1.

S5. Asymptotic dispersion of bands at infinity
For a general permitivitiy given by Eq. (1) of the main text, the wave equations ofM y -odd and even modes on the k y = 0 plane are, respectively, M y -odd: M y -even: where {⋅, ⋅} denotes the anticommutator. By introducing a new functionH εxx(ξ) dξ)H y , Eq.(53) can be transformed into a standard Sturm-Liouville equation: Therefore, the eigen-equations forM y even and odd bands can be uniformly expressed as follows: with Bloch boundary condition ψ (0) = e ikxL ψ (L), d dx ψ (0) = e ikxL d dx ψ (L), whereK is a positive semidefinite operator arising from the fact that the permittivity tensor of dielectrics is positive definite, and V (x) = V (x + L) is a positive piecewise smooth function (so it has infimum inf (V ) = V min ). Theorem 1. All bands ω n (k z ) of Eq.(55) tend to a unique asymptotic linear dispersion, as k z → ∞, and the asymptotic slope is determined by the infimum of V (x): Proof. 1) Lower bound: left product of ψ n (x) to Eq.(55) yields ω 2 n /c 2 = ⟨ψ n |K|ψ n ⟩ + k 2 z ⟨ψ n |V |ψ n ⟩ ≥ k 2 z V min , sinceK is positive semidefinite. Hence, we have ωn kz ≥ c √ V min .
On the other hand, ∀ϵ > 0, we can find n-D function space F n satisfying the boundary conditions, such that ⟨u|V |u⟩ As k z → ∞, both the lower and upper bounds of ω n /k z tend to inf ϵ>0 √ V min + ϵ = √ V min , so we obtain the limit lim kz→∞ ω n /k z = c √ V min . And due to L'Hôpital's rule, we also have lim kz→∞ dω/dk z = c √ V min .

Remark:
The theorem can also be proved using min-max theorem of eigenvalues [2]. Indeed, some procedures in the proof of min-max theorem have also been used in the above proof. In addition, Ref. [3] offers a different proof.
As a consequence, the asymptotic group velocities of even bands and order bands are, respectively, This asymptotic group velocities are numerically verified, as shown in Fig. 2d of the main text. Therefore, along Γ − Z direction, all bands tend to a linear dispersion as k z → ∞. All even and all odd bands have identical asymptotic group velocities respectively. For almost any dielectric photonic crystal respecting the symmetries, there exits a pair of triple points as the nexus of the nodal lines (intersection either of 1st even and 2nd,3rd odd bands or of 1st odd and 2nd,3rd even bands). The only exception (zero measure) occurs when the asymptotic group velocities of even and odd bands are identical, namely ε max xx = ε max yy .

S6. Robustness of nexus points against the variation of materials
In the main text, we have expounded that the special photonic band connectivity induced by the hidden symmetry guarantees that the triply degenerate NPs can almost always emerge on the 4 lowest bands along Γ − Z. In Fig. S1, we use the biaxial dielectrics to construct the AB-layer-stacked PhCs and fixed the values of ε xx , ε zz , and ε xz for the PhCs in all the panels, while we change the value of ε yy in each panel to show how the nodal structure changes with the parameter of the PhC and to demonstrate the robustness of NPs. We note that there are infinite NLs in the band structures, while only those connecting with the lowest triple NPs are plotted. In Fig. S1a, we let ε yy < ε xx . As such, the asymptotic group velocity ofM y -odd bands (light blue) is larger than that ofM y -even bands (light red) in the Γ − Z direction. So the lowest triply degenerate NPs (orange dots) are formed by the 1 st odd and 2 nd , 3 rd even bands at the intersections of two nodal rings (red and green) and the lowest Kramers-like NL (blue).
If we increase the value of ε yy while keeping it less than ε xx , the two nodal rings will grow bigger and will eventually intersect with their transnational counterparts at the fourfold degenerate NPs (purple dots) on the two boundaries of the BZ as shown in Fig. S1b. As a result, the nodal rings in the extended BZ connect together and form a nodal chain in the k y = 0 plane.
As ε yy increases, the eccentricity of the two nodal rings grows accordingly, and the two triple NPs move outward along the z axis from the origin. When ε yy reaches the critical value ε yy = ε xx , the asymptotic group velocities of even and odd bands are accidentally identical in the z direction. Then, the two NPs move to the infinity and vanish, and the two nodal rings are reduced to two straight lines parallel to Γ − Z, as shown in Fig. S1c. As we have emphasized in the main text, the two triple NPs in the lowest 4 bands can only disappear in this accidental condition of ε yy = ε xx . However, if we arbitrarily select the parameters of the PhCs respecting the symmetries, the probability of encountering these exceptional cases is zero, since they are restricted to a subset of measure zero for all possible parameters.
Once ε yy is larger than ε xx , Fig. S1d demonstrates that the two triple NPs will reappear immediately. However, in this case, the two NPs are on the band crossings of 1 st even and 2 nd ,3 rd odd bands, as the asymptotic group velocity of even bands surpasses that of odd bands. At the same time, Fig. S1d shows that the two previous ring shaped nodal lines convert to two hyperbolas. And the nodal hyperbolas intersect with their translational counterparts infinite times on the Γ − Z line as well as on the two BZ boundaries, forming fourfold degenerate NPs (magenta dots and purple dots).
In Fig. S2, we plot the band structure near a typical fourfold degenerate NP on the BZ boundary (X − M ). As shown in Fig. S2a,b, there are 4 nodal rings intersect on the NPs. In the section of k x = π/L, there are two doubly degenerate bands intersecting at the NP as displayed in Fig. S2c, resulting from the combined symmetry of time reversal and twofold screw rotationΘ L/2 = TŜ 2x . Furthermore, Fig. S2d,e show that both the band structure in the k z = k NP z section and the iso-frequency surfaces at ω NP disperse as 2D double Dirac cones [4]. Figure S2. Band structures near a fourfold degenerate nexus point in a AB-layer-stacked PhC with parameters ε1 = ε2 = 1, ε3 = 12 and θ = π/12. a. Band structure near the nexus of 4 nodal lines on the ky = 0 plane. b-d, Zoomed in band structure near the NP in the sections of ky = 0, kx = π/L, and kz = k NP z , respectively. e. Iso-frequency surfaces near the NP at the frequency of the NP.

S7. Derivation of k ⋅ p Hamiltonian for layer-stacked photonic crystals
Here we introduce a general k ⋅ p framework for a generic layer-stacked photonic crystal that is periodic in x direction and homogeneous in the other two dimensions r ⊥ = (y, z). In comparison with using wave equation (a 2-order PDE), it is more convenient to perform the derivation directly from the Maxwell's equations (1-order PDEs): Expending of the Bloch state Ψ n,k (r) = e ik⋅r u n,k (x) = e ik⊥⋅r⊥Ψ n,k (x) using the states at k 0 (usually a degenerate point), we obtain Substitution of the expansion into Eq. (58) yields Since the Brillouin zone of a layer-stacked PhC is infinite in the transverse plane, the orthogonality of the Bloch states should be written as where k and r are transverse wavevector and position vector respectively, δ(k ′ ⊥ − k ⊥ ) represents Dirac delta function, while the other two delta functions are Kronecker delta. Performing inner product between Ψ m ′ ,k0 and Eq. (60), we obtain the k ⋅ p eigen-equation of the 1 st -order which can be rewritten as with the element of the 1 st -order k ⋅ p Hamiltonian

A. k ⋅ p Hamiltonian near the nexus points
At the triply degenerate NPs with k NP± = (0, 0, ± 2π L √ εxx εyy−εxx ) and ω NP = 2cπ L √ εyy−εxx for the AB-layer-stacked PhC in Fig. 1 of the main text, the three degenerate eigenstates can be obtained from Eqs. (13) and (14): Indeed, the equality of the charge of the generated optical vortex and the difference of the spin quantum number of finial and initial states reflects the conservation of the generalized total angular momentum during the diffraction process. Since the effective 2D spin-1 Hamiltonian (75) has anisotropic Fermi velocity in the xy plane, the total angular momentum is not a conserved quantity of the Hamiltonian. Nevertheless, we can rewrite the Hamiltonian aŝ H(δk xy ) = v xŜ z δk x +ṽ yŜ y δk y = e a iŜa δk i , (a = z, y, i = x, y) with regarding the anisotropic Fermi velocity tensor as the tetard (e a i ) = v ↔ f = diag(v x ,ṽ y ) of an anisotropic space. The corresponding metric tensor of the space is given by (g ij ) = (δ ab e a i e b j ) = diag(v 2 x ,ṽ 2 y ). Then we can define a generalized orbital angular momentum operator for the anisotropic space: where ϵ zij denotes the antisymmetric symbol, (r i ) = (x, y). In terms of the coordinate transformation δke iϕ = v x δk x + iṽ y δk y wherek,ϕ can be viewed as the generalized polar coordinates in the momentum space, we have Therefore, in the momentum representation, the generalized orbital angular momentum operator can be alternatively expressed asL Thus the eigenstates of the generalized orbital angular momentum operator are exactly given by ψ l = exp [ilϕ(δk xy )] with l denoting the orbital angular quantum number: In addition, the generalized total angular momentum operator is given bỹ where the negative sign in font ofŜ x is utilized for the sake of the consistent chirality of the spin frame (Ŝ z ,Ŝ y , −Ŝ x ) and of the coordinate frame (x, y, z). Then it can be directly verified that the generalized total angular momentum commutes with the Hamiltonian: [J,Ĥ(δk xy )] = 0.
Therefore, the generalized total angular momentum quantum number j = l−s is conserved during the wave propagation along the z direction. For an incident state with l i = 0 orbital angular momentum and s i spin momentum, the output state satisfies j = l f − s f = 0 − s i , thus the final orbital quantum number l f = s f − s i is determined by the difference of the finial and initial spin quantum numbers.