Intertwined Weyl phases emergent from higher-order topology and unconventional Weyl fermions via crystalline symmetry

We discover three-dimensional intertwined Weyl phases, by developing a theory to create topological phases. The theory is based on intertwining existing topological gapped and gapless phases protected by the same crystalline symmetry. The intertwined Weyl phases feature both unconventional Weyl semimetallic (monopole charge>1) and higher-order topological phases, and more importantly, their exotic intertwining. While the two phases are independently stabilized by the same symmetry, their intertwining results in the specific distribution of them in the bulk. The construction mechanism allows us to combine different kinds of unconventional Weyl semimetallic and higher-order topological phases to generate distinct phases. Remarkably, on 2D surfaces, the intertwining causes the Fermi-arc topology to change in a periodic pattern against surface orientation. This feature provides a characteristic and feasible signature to probe the intertwining Weyl phases. Moreover, we provide guidelines for searching candidate materials, and elaborate on emulating the intertwined double-Weyl phase in cold-atom experiments.

In this work, we develop a theory to generate topological phases of matter, by intertwining existing ones protected by a crystalline symmetry. We discover intertwined Weyl phases, which feature the exotic interplay of unconventional Weyl fermions and higherorder topology.
In the bulk, the two phases are independently stabilized by the crystalline symmetry: the Weyl fermions have a monopole charge larger than 1 by the symmetry, and a higher-order topological phase was further superposed by the symmetry. The intertwining results in a specific distribution of the two phases in the bulk, i.e., the higher-order topology exists in the region outside pairs of unconventional Weyl points with opposite charges. Due to their independence, the two phases are separately tunable. The combination of different unconventional Weyl semimetallic and higherorder topological phases results in distinct intertwined Weyl phases. Note that this mechanism is different from the higher-order Weyl semimetal which is determined by the higher-order topology as discussed above. Intertwining topological phases by crystalline symmetry. Red and blue dots represent Weyl points, and green sold lines denote hinge states from higherorder topology. (a) The fourfold rotation symmetry (C4) stabilizes double-Weyl points and higher-order topology at the same time. The interplay between the two phases results in an intertwined double-Weyl phase. (b) The sixfold rotation symmetry (C6) stabilizes triple-Weyl points and higher-order topology at the same time, resulting in an intertwined triple-Weyl phase.
The intertwined Weyl phases are characterized by the the drastic change of Fermi-arc topology on 2D surfaces upon rotating surface termination. This phenomenon comes from the intertwining between the two constituent phases.
Even though higher-order topological phases feature boundary states on 1D hinges, they can intertwine with the unconventional Weyl phase on 2D surfaces of the system: Fermi arcs form due to the underlying unconventional Weyl phases, while the topology of these Fermi arcs is drastically changed by the higher-order topological phase. Thus, we identify a prominent topological feature of the intertwined phase: the change of Fermi-arc topology against surface orientation in a periodic pattern.
The period is determined by rotation symmetry, e.g., a 2n-fold rotation symmetry leads to a period of π/n. Specifically, we discuss intertwined double-Weyl phases, where double-Weyl fermions and higher-order topology are intertwined due to a fourfold rotation symmetry. The topological phase is characterized by the periodic change of Fermiarc topology with period π/2. The periodic behavior can be perfectly explained by an effective boundary Hamiltonian, in accordance with our theory. Finally, we show that intertwined double-Weyl phases are realizable in cold-atom experiments, which can serve as a promising platform to implement our theory in experiment.

Intertwining topological phases by crystalline symmetry
We begin by showing how to intertwine topological phases by crystalline symmetry. The essential physics can be captured by the following simple but generic Hamiltonian in 3D momentum space, where k ± = k x ± ik y , σ ± = (σ 1 ± iσ 2 )/2, and σ i and τ i (i = 1, 2, 3) are Pauli matrices for (pseudo)spin and orbital degrees of freedom, respectively. H Weyl , responsible for generating unconventional Weyl points, is superposed with mΓ h , a term for introducing higherorder topology. m is a real model parameter, and Γ h a 4×4 matrix. After being projected to the boundary, mΓ h acts as a mass term that is constrained by the symmetry, as we will discuss in Eqs. (3) and (4). Each of the two parts is invariant under the 2n-fold rotation symmetry C 2n about the z axis, with 2n ≤ 6 by lattice restriction (n ∈ {2, 3}). As shown in Figure 1, the intertwining can be induced by a crystalline symmetry that superposes the two phases at the same time: 1. The bulk Weyl points with higher monopole charge are stabilized by rotation symmetry. Unconventional Weyl points with charge ±n (n > 1) can be generated by breaking time-reversal and/or inversion symmetry in H Weyl . After the substitution of (k x , k y ) = (k cos θ, k sin θ) [see Figure 2(a)], H Weyl can be rewritten as H Weyl (θ) = k n e +inθ τ 3 σ − + k n e −inθ τ 3 σ + + k z τ 3 σ 3 . As can be seen, these unconventional Weyl points are protected by rotation symmetry [14] C 2n H Weyl (θ)Ĉ −1 2n = H Weyl (R 2n θ).
2. The higher-order topology is further enforced by rotation symmetry. The topological protection by symmetry can be explicitly demonstrated at the boundary. After projecting the Hamiltonian to the boundary subspace, we can get an effective Hamiltonian on the surface in the form of (See Supplementary Note 2 for derivation details.) where a i is a real coefficient, which may depend on k z , θ is the surface orientation, and k the momentum parallel to the surface [see Figure 2(a)].
The higher-order term m(θ)γ h , originated from mΓ h , is still rotation-symmetric,ĉ 2n m(θ)γ hĉ −1 2n = m(R 2n θ)γ h , with the projected rotation operator c 2n = σ 3 . Note that m(θ), which stems from the constant m term in Eq. (1), becomes θ-dependent after the projection by the orientation-dependent boundary wavefunctions. γ h anticommutes with σ 3 , e.g., σ 1 or σ 2 , so that it acts as a mass term in Eq. (3). Thus, we obtain This relation means m(θ) must have a zero value between θ and θ + π/n, which leads to a gapless point in the surface spectrum. The 1D Hamiltonian (3) for a fixed k z is characterized by a Z topological invariant. Note that the same mechanism leads to higher-order topological insulators protected by crystalline symmetry, where gapless points correspond to corner or hinge modes [43,44]. Here, as a manifestation of higher-order topology, the gapless points result in 1D hinge states, as shown by the green solid lines in Figure 1.
The two topological phases in superposition no longer act individually, but intertwine with each other, as we discuss next. Thus, we refer to the resultant phase as "intertwined Weyl phase", which is different from each individual phase.
We emphasize that in the intertwined Weyl phase, the unconventional Weyl points and the higher-order topology are independently stabilized by the symmetry. This is different from the higher-order Weyl semimetals, where the higher-order Weyl points depend on the higherorder topological phases and are typically conventional ones of topological charge-1 [51][52][53][54][55].
Due to the independent nature of the two phases, they can be tuned separately to generate various intertwined Weyl phases.
Specifically, different unconventional Weyl fermion phases can be incorporated to build with higherorder topological phase under the same or different crystalline symmetries, resulting in distinct kinds of intertwined Weyl phases, as we will discuss later.
Intertwining and the change of Fermi-arc topology in a periodic pattern In the bulk, the intertwining results in a specific distribution of the two phases. This is because the higher-order topology cannot survive in regions with non-trivial Chern number. To have vanishing Chern number, it is required that the unconventional Weyl points with opposite charges appear in pairs. In this way, there can exist regions with trivial Chern numbers to host higher-order topology. As shown in Figure 2 by green arrows, the higher-order topological (HOT) phase exists outside pairs of unconventional Weyl points with opposite charges.
Inside pairs of unconventional Weyl points, e.g., the lower and upper two sets of Weyl points in Figure 2 there is no m(θ), the surface Hamiltonian (3) describes Fermi arcs like those in unconventional Weyl semimetals. On 2D surfaces, the existence of these Fermi arcs is isotropic, i.e., they do not depend on surface orientation θ. However, in HOT regions, the higher-order topological term m(θ) exists. The existence of Fermi arcs in this region depends on m(θ) and becomes anisotropic. We find that due to the intertwining of the two parts of Fermi arcs, the Fermi-arc topology on 2D surfaces is drastically changed. We use m(θ) in Eq. (4) to explain the phenomenon in detail below.
First, the gapless point enforced by the higher-order topological phase in Eq. (4) crucially affects the Fermiarc topology. We use the spectrum of the surface Hamiltonian (3), E surface = ± ( n i=0 a i k i ) 2 + m 2 (θ), for illustration. If there is no gapless point [m(θ) = 0] in the HOT phase region for a specific angle θ, no Fermi arc can go through this HOT phase region because E surface is gapped, e.g., as shown by Figures 2 (b) and (d).
In contrast, if there is a gapless point [m(θ) = 0] in HOT region, Fermi arcs can go through this HOT region because E surface is gapless for some k , e.g., as shown by Figure 2 (c). In this regard, the gapless point [m(θ) = 0] determines whether the Fermi arcs can go through the higher-order phase region at an angle θ. Thus, strikingly, the Fermi-arc topology changes with θ.
Second, m(θ) in Eq. (4) is antiperiodic. Remarkably, it becomes periodic in the surface spectrum E surface after taking the square, since the square of m(θ) obeys m 2 (θ) = m 2 (θ + π/n), which has a period of π/n. The periodic function m 2 (θ) indicates that the topology of Fermi arcs not only changes with the surface orientation θ, but also in a periodic pattern. We emphasize that the change of Fermi-arc topology comes from the intertwining of unconventional Weyl fermions and higher-order topology. Thus, this phenomenon is absent in any individual phase alone. This is distinct from Ref. [59], where the edge states in different samples of different geometry are determined by higher-order topology. Thus, it constitutes the characteristic feature of intertwined Weyl phases.
After figuring out how m(θ) influences the behaviour of boundary states on the 2D surfaces, it is also necessary to discuss its impact on the boundary states on the 1D hinges. The hinge states, seen in the HOT phases, are determined by the condition m(θ) = 0. For m(θ) = 0, no localized zero modes can be observed. This is in accordance with the behaviour of the higher-order topological phase. Thus, the hinge states of higher-order topology can still be observed in 1D, as shown by the green lines in Figure 1.
A similar story applies to the intertwined triple-Weyl phase which is enforced by sixfold rotation symmetry C 6 , i.e., n = 3. On 2D surfaces, each triple-Weyl point emits three Fermi arcs on the surface, and the topology of Fermi arcs changes periodically with a period of π/3, as shown by Figures 2(e, f, g) in a full period [0, π/3]. On 1D hinges, the boundary states occur at, e.g., θ = nπ/6 with (n = 0, 2, 4, 6, 8, 10). They are in accordance with the boundary states of higher-order topological phase with C 6 symmetry, as shown in Figure 1(b).
In Figure 3 (a), the Chern number in k x k y -plane against k z is plotted. Here, without loss of generality, we assume k w2 < k w1 , for the chosen parameters given in the caption of Figure 3. We can see that the Chern number changes by 2 when passing a double-Weyl point. This is because the Chern number on the surface that encloses the Weyl point equals to the monopole charge, as it does for conventional Weyl phases. This explains the non-trivial Chern number between the upper pair k z ∈ (k w2 , k w1 ) and lower pair k z ∈ (−k w1 , −k w2 ) of Weyl points. However, the distribution of the Chern number in the intertwined phase is different from conventional ones. In order to host higher-order topological phases, it is required to have a region where the Chern number vanishes. Here, the region is between the lower and upper pairs of Weyl points, i.e., k z ∈ (−k w2 , k w2 ), as shown by the HOT region (green arrow) in Figure 3. In this region, the existence of boundary states depends on the orientation θ, as shown by Figure 3 (c). In contrast, in the region with non-zero Chern number, i.e., inside the two sets (upper and lower) of Weyl points, the Fermi arcs exist regardless of the angle θ.

Effective boundary theory
To understand the periodic behavior of Fermi-arc topology, a boundary theory applicable to any surface orientation is constructive. We can achieve this goal by firstly deriving two boundary states for each θ in the absence of the higher-order term mτ 1 σ 1 in Eq. (6). The spinor part of two boundary states takes the form of ψ 1 ∝ (e −2iθ , √ 2 + 1, 0, 0) approximately in the region of k z ∈ (−k w2 , k w2 ), and ψ 2 ∝ (0, 0, e −2iθ , √ 2 + 1) in the region of k z ∈ (−k w1 , k w1 ). Note that the contribution from spatial part of the boundary states does not affect the main results, and is neglected for simplicity. By projecting the whole Hamiltonian into the subspace spanned by the boundary states, we can obtain the effective boundary Hamiltonian as (See Supplementary Note 1 and 2 for derivation details.) up to a constant term −1/ √ 2 σ 0 , and k 2 c = 2t cos k z − M . Clearly, the boundary Hamiltonian is in the form of Eq. (3). The effective Hamiltonian is valid in the region where ψ 1 and ψ 2 overlap, i.e., k z ∈ (−k w2 , k w2 ) between the middle two Weyl points, where the Chern number is trivial. On a 1D hinge, the gapless points at m cos(2θ) = 0 result in hinge Fermi arcs, as shown by the green solid lines in Figure 1.
On 2D surfaces, clearly, the periodic change of Fermiarc topology, shown in Figures 2 (b)-(d), is caused by the higher-order term m cos(2θ). The period is π/2, as determined by m 2 cos 2 (2θ) in the spectrum, in accordance with the general theory of Eq. (5). Within a single period, the higher-order topology enforces the appearance of gapless point [Eq. (4)], which is located at m cos(2θ) = 0. The gapless point drastically changes the topology of Fermi arcs, because it determines whether the arcs can go through the region of k z ∈ (−k w2 , k w2 ) between the middle two Weyl points or not. Figures 2  (b)-(d) show three representative orientations in a full period of θ ∈ [0, π/2]. The Fermi arcs cannot go through the region between the two middle Weyl points at θ = 0. After rotating to θ = π/4 at the gapless point, the Fermi arcs are allowed to go through. Finally, at θ = π/2, the Fermi arcs return to their initial topology, completing one period.

Guidelines for material search
Two key ingredients of the discussed intertwined Weyl phases are fourfold (sixfold) rotation symmetry and the unconventional Weyl points on the rotation axis. Both requirements can be fulfilled in tetragonal, cubic, or hexagonal space groups, where the desired symmetries emerge. There is a correspondence between the absolute value of the chirality of a crossing and the rotation eigenvalues of the involved bands [62], such that a crossing between certain rotation eigenvalues leads to the required unconventional Weyl points. The total phase of eigenvalues accumulated by such crossings for any pair of bands is restricted by the periodicity of Brillouin zone. To be more specific, to obtain crossings of the required charge of ±2 (±3), the phase of the eigenvalues of fourfold (sixfold) rotation must change by π. By looking at Figure 1 it is evident that intertwined Weyl points require an even number of such changes of the eigenvalue phase. Thus, the necessary crossings can only occur if the total accumulated phase, which is consistent with the periodicity of the Brillouin zone, is equal to 0 mod 2π, like it is fulfilled in Eq. (6). This excludes for nonsymmorphic systems certain filling factors, where the k-dependence of rotation eigenvalues would, for example, require an odd number of crossings. An appropriate filling b for a 2n-fold rotation with a fractional translation of m 2n has to fulfill bm 2n ∈ Z. Commonly, a set of four crossings like in Figure 1 would be accidental. Nevertheless, enforced crossings can be found within fourfold degenerate points, e.g., in space groups 106 and 133 on the path M-A. These comprise unconventional Weyl points related by mirror symmetry [63]. Hereby, mirror pairs of enforced unconventional Weyl points coincide. If a weak time-reversal and mirror symmetry breaking is introduced to such a system, one may obtain Weyl points in the configuration discussed in Figure 1(a) with no other enforced bands at the Fermi energy. A comprehensive material search is left for future works.

Cold-atom experimental realization
Owing to technical advances, cold atoms have been widely applied in quantum simulations of topological matter [64][65][66][67], and now are also readily available for realizing our intertwined double-Weyl semimetal described by Eq. (6). Here we present the realization proposal using fermionic atoms. We choose two hyperfine states as the pseudo-spins ↑↓ for the σ degrees of freedom. For our model Hamiltonian, which requires two extra degrees of freedom, we consider the atomic gases loaded in a 2D bilayer optical lattice, and thus, the layer index λ is used to represent the τ subspace. The setup is illustrated in Figure 4. The energy offset of the hyperfine states is prepared as Γ λ (φ) = M 0 + (−1) λ − 2t cos φ that is manually controllable by the cyclical parameter φ. It not only generates the on-site energy of the model Hamiltonian, but also introduces the additional parameterized space represented by φ. In this way, the band physics of the Weyl Hamiltonian can be captured in the 3D parameterized space (k x , k y , φ) [68]. Here, φ corresponds to k z in Eq. (6).
In order to engineer the intra-layer spin-flipped hopping, we use laser fields of three modes to couple the pseudo-spins. The spatial modulations of the field modes are prepared as M 1 (r) = iM 1 sin(k L x) sin(k L y), M 2 (r) = M 2 cos(k L x) cos(3k L y) cos(k L z), and M 2 (r) = −M 2 cos(3k L x) cos(k L y) cos(k L z), where k L = π/d and d denotes the lattice constant. Due to the odd parity of M 1 (r) in the xy-plane [69], the on-site and nearestneighbor (NN) couplings vanish, while the next-NN coupling dominates, resulting in sin k x sin k y σ 1 . Due to the crystalline symmetry, the combination of M 2 (r) and M 2 (r) leads to the NN coupling (cos k y − cos k x )σ 1 . After making operator transformations, all the intralayer terms host opposite signs for different layers. Furthermore, the higher-order topological term τ 1 σ 1 is naturally introduced by the inter-layer hopping. The details of the realization proposal are shown in the Supplementary Note 3.

DISCUSSION
Based on the theory to create topological phases by intertwining existing ones, we have discovered the intertwined Weyl phase. The intertwined Weyl phase is different from the individual unconventional Weyl semimetallic phase or higher-order topological phase, and exhibits its distinct characteristic topological features.
The intertwining results in the drastic change of Fermiarc topology in a periodic pattern, which constitutes the characteristic feature of the intertwined phase. We have proposed a feasible cold-atom experiment to verify our theory and to realize the intertwined Weyl phases.
Our theory could serve as a guiding principle to generate topological phases based on existing ones. A direct application would be to investigate the intertwining between topological gapless phases with emergent particles other than Weyl fermions and topological crystalline phases, that are protected by the same symmetry [28-30, 70, 71].
Finally, we note that Fermi-surface topology is crucial for electronic properties of material. The change of Fermi-arc topology in the intertwined Weyl phases offers a direction of tuning Fermi-surface topology. It would inspire further research on electrical, magnetic, thermodynamic, and transport properties, that are determined by Fermi-surface topology. For instance, a potential application is to investigate the quantum oscillations of Fermi-arc surface states [72]. These oscillations are periodic against the inverse of magnetic field 1/B. Their frequency F is determined by the Onsager relation F = Φ 0 /(2π 2 )A s . Here Φ 0 = h/2e is the magnetic flux quantum, and A s is the Fermi surface cross section. In intertwined Weyl phases, the cross section A s depends on the surface orientation θ. As shown in Figure 3 (c), by changing θ, one may interpolate between small Fermi surfaces (localized around k w1 and k w2 ) and extended Fermi surfaces (connecting k w1 (2) to −k w1(2) ), which must yield a substantial change in the observed quantum oscillation spectra. Thus, the intertwined Weyl phases are expected to have a significant change of quantum oscillations upon rotating surface termination.

Analytical derivation and symmetry analysis
The derivation of the analytical results makes use of the k · p methods and symmetry analysis. The numerical calculation is based on tight-binding model and Green's function method. Details of these derivations are given in the Supplementary Notes 1, 2, and 4. This includes the calculation of Fermi arcs in intertwined double-Weyl semimetals, and the analysis of symmetry for intertwined double-Weyl and triple-Weyl semimetals. The Supplementary Note 3 also contains a detailed discussion of the cold-atom realization of intertwined double-Weyl semimetals.

DATA AVAILABILITY
The numerical data of the plots within this paper are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The numerical codes that support the findings of this paper are available from the corresponding author upon