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.

Though topological gapless phases characterized by nontrivial band degeneracies, and topological gapped phases featuring anomalous boundary states, can be stabilized by the same crystalline symmetry, so far, their interplay remains to be explored.
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 higher-order 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 HOT 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 HOT phases results in distinct intertwined Weyl phases. Note that this mechanism is different from the higher-order Weyl semimetal which is determined by the higherorder topology as discussed above.
The intertwined Weyl phases are characterized by 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 HOT 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 HOT phase. Thus, we identify a prominent topological feature of the intertwined phase: the change of Fermiarc 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 Fermi-arc 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 coldatom 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 higher-order 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 Fig. 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 θÞ (Fig. 2a), 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 Here,Ĉ 2n ¼ τ 3 σ 3 , and R 2n θ = θ + π/n acts in momentum space.
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 (Fig. 2a). 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ĉ 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 mðθÞ ¼ ÀmðR 2n θÞ ¼ Àmðθ þ π=nÞ: 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 HOT 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 Fig. 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 HOT 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 HOT 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 higherorder topology. As shown in Fig. 2 by green arrows, the 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 Fig. 2b, 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 Red and blue dots represent Weyl points, and green sold lines denote hinge states from higher-order topology. a The fourfold rotation symmetry (C 4 ) stabilizes double-Weyl points and higherorder topology at the same time. The interplay between the two phases results in an intertwined double-Weyl phase. b The sixfold rotation symmetry (C 6 ) stabilizes triple-Weyl points and higher-order topology at the same time, resulting in an intertwined triple-Weyl phase. orientation θ. However, in HOT regions, the HOT 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 HOT phase in Eq. (4) crucially affects the Fermi-arc topology. We use the spectrum of the surface Hamiltonian (3) , 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 Fig. 2b 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 Fig. 2c. 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 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 behavior 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 behavior of the HOT phase. Thus, the hinge states of higher-order topology can still be observed in 1D, as shown by the green lines in Fig. 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 Fig. 2 The intertwining results in the change of Fermi-arc topology in a periodic pattern. a The relation between different coordinate systems. θ indicates the orientation of the surface (black plane), on which the local density of states is computed numerically. b-d Fermi arcs of intertwined double-Weyl phase in a full period of θ ∈ [0, π/2]. Each projected double-Weyl point emits two Fermi arcs on the surface. The topology of Fermi arcs changes when rotating from the (100) surface (θ = 0) to the (110) surface (θ = π/4), and it returns to the initial topology after rotating further to the (010) surface (θ = π/2), completing one period. e-g Fermi arcs of intertwined triple-Weyl phase in a full period of θ ∈ [0, π/3]. Each projected triple-Weyl point emits three Fermi arcs on the surface. The topology of Fermi arcs changes when rotating from the (1000) surface (θ = 0) to the ð1010Þ surface (θ = π/6), and it returns to the initial topology after rotating further to the ð0010Þ surface (θ = π/3), completing one period.
surface, and the topology of Fermi arcs changes periodically with a period of π/3, as shown by Fig. 2e-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 the HOT phase with C 6 symmetry, as shown in Fig. 1b.
In Fig. 3a, 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 Fig. 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 HOT 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 2 Àk w2 ; k w2 ð Þ , as shown by the HOT region (green arrow) in Fig. 3. In this region, the existence of boundary states depends on the orientation θ, as shown by Fig. 3c. 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θ ; ffiffi ffi 2 p þ 1; 0; 0Þ approximately in the region of k z ∈ (−k w2 , k w2 ), and ψ 2 / ð0; 0; e À2iθ ; ffiffi ffi 2 p þ 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= ffiffi ffi 2 p ϵσ 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 2 ðÀk w2 ; k w2 Þ between the middle two Weyl points, where the Fig. 3 The location of Weyl points and the orientation-dependent mass term. a Left: The four double-Weyl points on the k z axis. Red and blue dots denote Weyl points with monopole charge +2 and −2, respectively. Right: Chern number on k x k y -plane against k z , which changes by ± 2 when a double-Weyl point is passed. It is nontrivial within upper and lower two pairs of double-Weyl points. b The higher-order term mðθÞ 2 ¼ m 2 cos 2 ð2θÞ is a periodic function with period π/2. c Angular dependence of Fermi arcs in the range of [0, π/4], which is between Fig. 2b  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 Fig. 1.
On 2D surfaces, clearly, the periodic change of Fermi-arc topology, shown in Fig. 2b-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. Figure 2b-d shows 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 Fig. 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 kdependence 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 2 Z. Commonly, a set of four crossings like in Fig. 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 Fig. 1a 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 Fig. 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 0 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 nearest-neighbor (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 0 2 ðrÞ leads to the NN coupling ðcos k y À cos k x Þσ 1 . After making operator transformations, all the intra-layer terms host opposite signs for different layers. Furthermore, the HOT term τ 1 σ 1 is naturally introduced by the inter-layer hopping. The details of the realization proposal are shown in 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 HOT phase, and exhibits its distinct characteristic topological features. The intertwining results in the drastic change of Fermi-arc 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][29][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. Fig. 4 Illustration of experimental setup. The atoms are confined in a bilayer optical lattice. Couplings between opposite pseudospins are processed via laser fields of modes M 1 (r), M 2 (r), and M 0 2 ðrÞ. The on-site energy offset is prepared as Γ λ (ϕ) which not only depends on the layer index but also is manually controlled by the parameter ϕ. In intertwined Weyl phases, the cross section A s depends on the surface orientation θ. As shown in Fig. 3c, 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.