Satellite Dirac cones in phononic crystals

Domain: physics
Published: 2025-12-18
Editorial outcome: minor_revision

Abstract
Satellite Dirac cones in graphene systems are challenging to observe due to limited momentum-space resolution, hindering further exploration of their bulk-edge correspondence and associated transport properties. Here, we tackle this challenge by simplifying the bilayer configuration into a monolayer one with strong long-range couplings—an arrangement unattainable in natural materials but achievable in artificial structures for classical waves. Using phononic crystals for acoustic waves, we experimentally observe the satellite Dirac cones in the presence of third-nearest neighbor couplings and identify the acoustic edge states determined by their distinct non-trivial winding numbers. The satellite valleys emerge by introducing extra on-site potential, leading to multiple valley interface states. We further identify a frequency window enabling anomalous forward transport of the valley interface states in a standard four-channel splitter. Note that the introduction of long-range couplings is a systematic approach to engineering topological properties. Our findings establish a promising framework for multiple Dirac cones and valleys, opening new avenues for wave manipulation and signal processing in topological materials across both classical and quantum systems.

Introduction
Following the rise of graphene, layer-stacked graphene came up to the stage as a platform for exploring a plethora of intriguing physics phenomena1–4, such as unconventional quantum Hall effects5,6 and domain-wall valley transports7–9. These phenomena are closely tied to the quadratic Dirac cone. It is demonstrated that, in the presence of trigonal warping effect, the quadratic Dirac cone is unstable and evolves into a central linear Dirac cone accompanied by three additional linear Dirac cones surrounding it, collectively referred to as satellite Dirac cones1,6,10–13. Despite being proposed nearly two decades ago, satellite Dirac cones are experimentally challenging to identify, and their direct spectral observation remains elusive14.

Meanwhile, the Dirac cones inherently carry topological characteristics, as exemplified by the presence of edge states in graphene via the principle of bulk-edge correspondence15,16. Additionally, introducing gaps in the Dirac cones lead to topological insulating phases, featuring by topologically protected edge or interfacial states17. Dirac cones and the associated topological transport are not confined to electronic systems and can be generalized to classical wave systems18–20, where the flexibility in structural design allows the relaxation of locality constraints in interactions. For example, roton-like dispersions can be driven by dominant long-range couplings in metamaterials21–25, which enable the emergence of larger topological invariants26–30. From this perspective, an intriguing question arises: could symmetry-permitted longer-ranged couplings generate additional Dirac cones in the Brillouin zone (BZ), forming satellite Dirac cones? If so, what would be the associated bulk-edge correspondence, and what novel phenomena might emerge if these satellite Dirac cones are gapped into satellite valleys?

In this work, we achieve the realization of satellite Dirac cones by introducing strong third-nearest-neighbor (TNN) couplings. These couplings enable the symmetry-protected Dirac cones in a monolayer honeycomb structure to generate satellite Dirac cones, accompanied by a band reshaping. The resulting satellite Dirac cones carry Berry phases in an opposite sign to that of the central cone. Consequently, satellite Dirac cones induce a doubling of edge states compared to the original one. Our model is implemented in phononic crystals (PCs), which have long served as versatile platforms for exploring topological phenomena19,20. PCs offer flexibility for performing broadband measurements and facilitate counterintuitive couplings, including the strong TNN couplings required in this study. Moreover, by introducing staggered on-site potentials, the Dirac cones are gapped, forming satellite valleys that support multiple valley interfacial modes. The satellite Dirac cones and their associated transport phenomena are validated through both simulations and experiments. These findings highlight the potential of longer-range couplings in band engineering, enabling the emergence of topologically nontrivial phases.

Results
Lattice model with satellite Dirac cones
We incorporate TNN hoppings into a monolayer graphene lattice to achieve satellite Dirac cones, which significantly simplifies the scenario in Bernal-stacked bilayer graphene (see Supplementary Note 1). As illustrated in Fig. 1a, our model yields the following Bloch Hamiltonian in the (A, B) sublattice basis:

H(k) = [0 t1 γ1(k) + t3 γ3(k); t1 γ1*(k) + t3 γ3*(k) 0], (1)

where t1 and t3 represent the nearest-neighbor and TNN hopping strengths, respectively, both of which are real-valued. The explicit forms of γ1(k) and γ3(k) are given as γ1 = 1 + eik·a2 + eik·(a2 − a1) and γ3 = eik·a1 + eik·(a1 − a2) + e−ik·a2. The energy spectra are given by E± = ±|t1 γ1(k) + t3 γ3(k)|, and the band-touching points are determined by solving E± = 0. The existence of these solutions is guaranteed by the combined symmetries of time-reversal and spatial inversion. In real graphene, the TNN hopping strength (t3) is approximately 40 times weaker than t1, and thus t3 is often neglected when analyzing electronic properties of graphene. This approximation gives rise to the well-known Dirac cones at the corners of the BZ.

To explore the impact of stronger beyond-neighbor couplings (t3 > t1), we show the evolution of the Dirac points as t3 increases, as depicted in Fig. 1b. Notably, this evolution induced by TNN hoppings resembles a Lifshitz transition in non-interacting systems31, and, in twist-bilayer graphene near magic angle, similar Dirac cone evolution and Berry curvature distributions can be observed32. Starting from t3 = 0, where a regular Dirac cone is located at K, increasing t3 gradually reduces the gap at M until it closes at t3 = t1/3. This gap closure corresponds to the formation of a semi-Dirac cone characterized by an anisotropic Fermi velocity. As t3 continues to increase, the semi-Dirac cone splits into two Dirac cones along K-K0. Subsequently, when t3 = t1/2, these Dirac cones merge with the central cone at K, resulting in a quadratic node, which accounts for the low-energy band structure in Bernal-stacked bilayer graphene without trigonal warping. Further increasing t3, the satellite Dirac cones reappear from the quadratic band crossing along the ΓK directions, resembling the trigonal warping effect. A summary of these Dirac cones with their low-energy Hamiltonians is provided in Supplementary Note 2 with Table S1. Hereafter, we concentrate on the regime of stronger t3, using t3 = 0.36t1 as a representative case.

We investigate the topological property and bulk-edge correspondence of the satellite Dirac cones. The topological invariants are defined by the winding number W(kx) = (2πi)−1 ∫ dky ∂ky [log(t1 γ1 + t3 γ3)], as shown in Fig. 1e. Beyond W=1 regions, higher W=2 regions are clearly observed, which originates from a twice change in phase induced by additional satellite Dirac cones. The winding number can also be intuitively viewed by counting the winding number of pseudospin vector as ky traverses the BZ. Under a proper gauge choice, one can define a pseudospin vector n^(k) = (⟨σy⟩, −⟨σx⟩) and track its evolution, where ⟨σs⟩ represents the expectation value of the Pauli matrices for the eigenstate |u+(k)⟩. For example, as shown in Fig. 1f, at kx = 0.25, the pseudospin vector completes two full twists before returning to its starting point, indicating W=2. In contrast, at kx = 0.36, the pseudospin vector exhibits a single twist, corresponding to W=1. For kx with trivial W values, the pseudospin vector shows no winding at all. Consequently, a doubled number of edge states occur at kx = 0.25 in Fig. 1g, compared to a single pair at kx = 0.36. We further visualize the distributions of these edge states in Fig. 1h, where edge states with larger winding numbers exhibit shorter decay lengths. By extending the coupling range, additional Dirac points and higher winding numbers become attainable, as summarized in Table S2 of Supplementary Note 8, highlighting that long-range couplings provide a simple and systematic approach for engineering topological band degeneracies and manipulating topological invariants33,34.

Acoustic Realizations
To realize satellite Dirac points experimentally, we design an acoustic analog of graphene with enhanced TNN couplings. Figure 2a presents a photograph of our 3D-printed sample, where cylindrical cavities are connected by pipe to form a honeycomb structure with 14×11 unit cells in the x-y plane, with a lattice constant of a = 25 mm. Figure 2b shows the schematic of the sample design for realizing stronger acoustic TNN couplings (see Supplementary Note 3 with Fig. S3 for geometric details). All cavities in the structure are identical in height and diameter, ensuring that the unit cell belongs to the point group D3, which guarantees the band crossings at BZ corners of PC. In the experiment, the system is driven near the frequency of the first longitudinal cavity mode, which features a wave node at the midplane of the cavity (lower-right inset in Fig. 2a). This mode strengthens the coupling for tubes connected at the cavity ends, while reducing the coupling for tubes connected near the midplane. This design facilitates the required TNN couplings for achieving the acoustic satellite Dirac points.

The measured bulk dispersion is depicted in Fig. 2c, which matches well with the simulated and analytical results represented by orange solid and white dashed lines. These results confirm the appearance of extra Dirac cone dispersions along the Γ-K line, arising from the interplay between strong TNN couplings and trigonal crystalline symmetry. Figure 2d, e present the measured acoustic dispersions along K-Γ1 and K-Γ2, further confirming the presence of satellite structure in band crossings around frequency fD = 4.94 kHz. The measured frequency contours near K and K0 are present in Supplementary Fig. S4. Note that the absence of genuine chiral symmetry in the real system leads to slightly difference for the dispersions above and below fD. Nevertheless, our PC captures the key features and profiles of the satellite Dirac cones with high fidelity. One can further improve the resolution of measured dispersions by using samples with more unit cells.

We further experimentally validate the bulk-edge correspondence in the gapless PC with satellite Dirac cones. Excited by a point-source placed in an edge cavity of the PC, the scanned acoustic pressure field is displayed in Fig. 3a. The measured field decays more rapidly toward the bulk than along the edge, indicating the existence of edge states. The response spectra of the sample in edge and bulk area are given in Fig. S5. Figure 3b shows the projected dispersion as a function of kx, where the colormaps represent experimental results and the dots correspond to simulated dispersions for a strip with bearded edges. Around the Dirac frequency fD, the bright colormaps align with the simulated edge states along kx, covering the wavevectors between satellite Dirac point projections, as illustrated in Fig. 3c. Specifically, at kx = 0.25, four edge modes are present, with the simulated acoustic pressure fields shown in Fig. 3d (only the upper boundary is displayed for clarity). Similarly, at kx = 0.36, two edge modes are observed with the corresponding fields shown in Fig. 3e. Simulated results for acoustic edge states with zigzag boundaries are provided in Supplementary Fig. S6. These simulated and experimental results show good agreement, providing support for the satellite Dirac points in our acoustic system.

In experiments, the response signal near the sample boundary generally results from a collective contribution of all edge modes. This complicates a clear experimental identification of the kx-dependent bulk-edge correspondence solely based on edge responses, as shown in Fig. 3b. To address this subtlety, we measure the variation of pseudospin (⟨σx⟩, ⟨σy⟩) across a ky path at three specific kx, marked by colored solid lines in Fig. 3c, which provides an intuitive representation of the winding numbers that corresponds to the number of edge states, as discussed in the lattice model. Figure 3f–h show the gapped bulk dispersions along ky at three representative kx values, which can alternatively be interpreted as the dispersions of reduced 1D subsystems parameterized by different kx. For each subsystem, the spinor used to construct ⟨σi⟩ is approximated by Fourier transforming the response of a single band from the two sublattices separately (see Methods and Supplementary Note 4 with Fig. S7). The corresponding pseudospin trajectories for the first band over a ky-cycle are depicted in Fig. 3i–k. In Fig. 3i, the pseudospin winds around the origin once, corresponding to the presence of the two edge states at kx = 0.36. At kx = 0.25, where four edge states are present, a double winding of the pseudospin spiral is clearly observed in Fig. 3j. Finally, at kx = 0, the pseudospin trajectory oscillates back and forth without forming a winding, as illustrated in Fig. 3k, which accounts for the absence of edge state responses. To show the extensibility of our approach, we further experimentally identified a winding number of 3 in a similar honeycomb PC with strong eighth-nearest-neighbor couplings (see Supplementary Note 8), further confirming that longer-range couplings enable larger winding numbers.

The satellite-valley interface states
We further explore the acoustic valley interface states related to the satellite Dirac cones. A typical sample with a zigzag-type interface is shown in Fig. 4a. The staggered sublattice potential is introduced by varying the cavity heights, hB − hA = ±δ with δ = 0.05h0 for the two distinct phases I and II. As such, the satellite Dirac cones split into satellite valleys located around K(K0) with opposite Berry curvatures, resulting in multi-valley insulating phases, as shown in Fig. 4b. Consequently, as shown in Fig. 4c, the valley interface channel holds four guiding modes, each exhibiting two valley polarizations propagating in the same direction (see Supplementary Note 9 for the lattice model results). Distinct from approaches rely on high-order modes35, these interface states arise from satellite valleys induced by long-range couplings, enabling a wider frequency window. The robustness of these valley guiding modes against sharp path bends is demonstrated in Supplementary Note 9.

As is well known, the preservation of valley pseudospin enables valley-locked one-way topological transport, which plays a key role in governing wave behavior in a prototypical four-channel beam splitter. These splitters serve as fundamental building block for functional topological valley-based networks. Note that in the multi-valley cases, the Berry curvature becomes more diffuse, facilitating intervalley transitions. We exploit these transitions to enable transmission through channels otherwise forbidden by valley locking. As illustrated in Fig. 4d, we construct such a beam splitter by combining two distinct phases (see Supplementary Note 5 with Figs. S8 and S9 for the configuration details of the splitter and eigenfields in the channels), with the input channel as shown in Fig. 4a. We then simulate the acoustic interface state carrying valley index q at the input channel and examine how it bifurcates at the intersection. Remarkably, within a frequency window of approximately 4.73–4.78 kHz, the acoustic interface state can transport forward into channel 4, as demonstrated by both simulations and experiments shown in Fig. 4d, e. Note that such forward transport is previously heavily suppressed due to the mismatched valley polarization, while here the satellite valleys enable this anomalous forward transport. Outside this frequency window, the system exhibits normal wave-splitting behavior (see Supplementary Note 6 for the channel transmissions at different frequencies). The measured Fourier spectra at 4.75 kHz for the channels are presented in Fig. 4f–i. For a q-valley polarized input acoustic wave, the output waves from channels 2 and 3 retain the valley polarization as expected (Fig. 4g, h). However, the acoustic wave in channel 4 is converted to be q0-valley polarized (Fig. 4i), breaking the valley conservation. This is probably caused by intervalley scattering, which occurs due to reduced momentum-space separation between the valleys36. (see Supplementary Note 7 with Fig. S13 for simulated results of the splitter with different interface configurations). This phenomenon may add fresh ingredients for valleytronics and valley-based transportation.

Methods
Experiments
The PC samples were fabricated by 3D printing, using photosensitive resin. The wall thickness of all cavities and connecting tubes was set at 2 mm. Due to the significant impedance mismatch, the walls act as rigid boundaries for sound waves. Small holes (3 mm in diameter) were drilled on top of each cavity to facilitate the excitation and probing of sound signals; these holes were sealed with stoppers when not in use.
In experiments, a network analyzer (B&K Type 3560B) was employed to generate and analyze the signals in the frequency domain, with a frequency step of 4 Hz. Sound generation was achieved using loudspeakers (about 2 mm in diameter) inserted into designed cavities, while sound detection was carried out using two acoustic probes. One probe (B&K Type 4187 1/4-inch microphone) was fixed near the sound source to collect the reference signals, while the other (B&K Type 4182) was manually swept in the sampling area of desired cavities to measure the field distributions (Fig. 4e). Spatial Fourier transformations of the acoustic fields were performed to extract the corresponding dispersions for bulk states (Fig. 2c–e), edge states (Fig. 3b), and interface states (Fig. 4c and Fig. 4f–i).

The experimental method for extracting pseudospin windings in Fig. 3i–k is detailed as follows. A point source was placed in a bulk cavity of the sample shown in Fig. 2a to generate sound waves (see Supplementary Fig. S5 for the measured response spectrum). The acoustic fields in sublattice cavities A and B were separately recorded and subjected to Fourier transformations. From the transformed data, the one-dimensional bulk band versus ky was extracted at kx = kxi as shown in Fig. 3f–h. The response intensity of a single band at wavevector kyi and frequency fi was approximated by a peak, which was subsequently Fourier transformed to yield extended waves in real space y. Using the transformed data, the field amplitudes at the same y from the two sublattices were adopted to construct the spinors, thereby obtaining the pseudospin vector at kyi. This procedure is described in Supplementary Note 4 with Fig. S7. This process was repeated for each kyi and yielded the pseudospin windings shown in Fig. 3i–k.

Simulations
All simulations of PCs were implemented with COMSOL Multiphysics, a finite element solver package. The speed of sound in air was set to 344 m/s at room temperature and the air density was 1.29 kg/m3. The solid walls of the sample were treated as rigid boundaries. To obtain the bulk dispersions (orange line in Fig. 2c), periodic boundary conditions were applied to the unit cell in both directions. For the projected bands corresponding to edge states (dotted lines in Fig. 3b), a strip-like supercell consisting of 11 unit cells along the y direction was used, with solid boundary conditions on the edges and periodic boundary conditions along the x direction. The calculation of valley interface states in Fig. 4c was carried out in a similar way, with states localized at the external edges filtered out for clarity. To simulate the frequency-domain field distributions shown in Fig. 4d, two identical point sources were placed in two interface cavities at x = −10. Sound damping in simulation was incorporated by assigning an imaginary part of the sound velocity with a value of 0.5i m/s. The height of the point sources was set to 0.9(hA − δ) for exciting the acoustic field.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The work entitled “Satellite Dirac cones in phononic crystals” reports the design and realization of a two-dimensional phononic crystal hosting multiple spinless Dirac points. The Dirac points are induced by a long-range hopping that preserves the chiral symmetry. These Dirac points, together with the induced zero-energy modes, were observed through direct dispersion measurements. The authors also explored the gapped phase where multiple valley interface states were discovered. This paper is easy to read, with sufficient supporting numerical results and experimental details. In general, I believe this is a nice work that would contribute to the areas of phononic crystal and topological acoustics.\nHowever, I am concerned about the novelty and general interest of this work. Firstly, systems containing multiple Dirac points are not new, even in experiments (e.g., Ref. [14] in the present paper and Science 361, 782 (2018)). Also, the observed phenomena are standard consequences of Dirac points and can be similarly found in the literature. For example, mid-gap edge states with high winding numbers have been observed in photonics (Adv. Photon. 7, 046005 (2025)). Interfaces supporting multiple valley interface states have also been realized (e.g., Laser Photon. Rev. 18, 2300686 (2024)). Therefore, I feel that simply a new design for these known phenomena and a demonstration of the splitting function cannot justify publication in Nature Communications.\nA technical issue: the origin of the valley interface states is not thoroughly discussed. Can they be characterized by the valley Chern number? Or can they be obtained from the low-energy theory? I notice that the projections of the valleys are quite close. Will this cause any potential trouble or effect?", "reviewer_2": "This work studies Satellite Dirac cones by simplifying the bilayer configuration into a monolayer one, which is realized by introducing third-nearest neighbor couplings. By using phononic crystals for acoustic waves, they experimentally observe the satellite Dirac cones, and identify related properties like the acoustic edge states determined by distinct non-trivial winding numbers. This is an interesting work, trying to avoid the pitfall of complex calculation in bilayer graphene. I will recommend this paper to be published after solving following concerns.\nIn Fig.2, are the acoustic dispersions gapless? please give the partial enlarged graph around K point and explain the reason.\nIn Fig.3d, the acoustic eigenfields are different with a same k_x. Why?\nWhy the authors employ the Bernal-stacked bilayer graphene but not the turbostratic stacked (AA) bilayer graphene? Is the developed method also applicable to the latter?\nHow to quantitatively demonstrate the equivalence between BBG and the TNN model? From Fig. S1, they are not perfectly equivalent. Will this difference lead to distinct properties in Dirac cones? Furthermore, can this simplified TNN model replace Bilayer models in other physical phenomenon? Or they are only equivalent in Dirac cones prediction?\nIn the supplementary information Note 6, please give the specific value of the group velocity V_g and other physical quantities, which can make the readers to reproduce relevant data. Also, the author should also check the main text and give the specific value of related physical quantities.\nHow to accurately define this word “neighbor” ? In one periodic cell, there are six carbon nodes (if numbered 1-6 in anticlockwise direction). For node 1, it is related to node 2 and node 6 by conventional interaction and related to node 4 by the third neighbor effect. Then what about node 3 and node 5? Are they negligible?\nSome additional minor questions:\nIn Fig.1d, the author obtains the projected band spectrum through the path “q-K-q’-q-K’-q”, but I do not find any illustration about the q and q’ point. please explain it or mark them in the figure.\nIn Fig.4d and Fig.4e, I would recommend the author use a same color bar to facilitate the comparison between simulation and experimental results.\nIn recent years, the magic-twisted bilayer graphene is a research hotspot. Can the researchers do some related discussions? Note that twisting may generate flat band, which will cause some exotic phenomenon. Is this TNN model applicable to flat band or superconductivity study?"}, "author_response": "RESPONSE TO REVIEWER #1:\nThe work entitled “Satellite Dirac cones in phononic crystals” reports the design and realization\nof a two-dimensional phononic crystal hosting multiple spinless Dirac points. The Dirac points\nare induced by a long-range hopping that preserves the chiral symmetry. These Dirac points,\ntogether with the induced zero-energy modes, were observed through direct dispersion\nmeasurements. The authors also explored the gapped phase where multiple valley interface\nstates were discovered. This paper is easy to read, with sufficient supporting numerical results\nand experimental details. In general, I believe this is a nice work that would contribute to the\nareas of phononic crystal and topological acoustics.\nReply: We thank Reviewer #1 for the constructive summary and the positive assessment of our\nwork as “a nice work that would contribute to the areas of phononic crystal and topological\nacoustics”, supported by “sufficient supporting numerical results and experimental details”. We\ngreatly appreciate the recognition of both theoretical and experimental contributions. Below,\nwe provide a point-to-point response to the concerns and questions.\nHowever, I am concerned about the novelty and general interest of this work. Firstly, systems\ncontaining multiple Dirac points are not new, even in experiments (e.g., Ref. [14] in the present\npaper and Science 361, 782 (2018)). Also, the observed phenomena are standard consequences\nof Dirac points and can be similarly found in the literature. For example, mid-gap edge states\nwith high winding numbers have been observed in photonics (Adv. Photon. 7, 046005 (2025)).\nInterfaces supporting multiple valley interface states have also been realized (e.g., Laser Photon.\nRev. 18, 2300686 (2024)). Therefore, I feel that simply a new design for these known\nphenomena and a demonstration of the splitting function cannot justify publication in Nature\nCommunications.\nReply: We thank Reviewer #1 for the opportunity to clarify the novelty and general interest of\nour work. Our contribution goes beyond a new design, we introduce a novel systematic\nmechanism to generate multiple topological degeneracies and large winding numbers, which is\nfundamentally distinct from, and offers significant advantages over, the references raised by\nReviewer #1. Our core novelty lies in the introduction of strong long-range couplings to\ngenerate topological phases. The long-range couplings can induce roton-like effect that reverts\nband dispersions [cf. e.g., Nat. Rev. Phys. 7, 299 (2025)]; however, existing works focused on\ngapped insulator phases, our work is the first to apply this concept to semimetal cases,\nunlocking a simple, general, and highly tunable method to engineer topological properties.\nTo demonstrate the advantage of our approach, the main text reports 3rd NN coupled\nacoustic graphene with eight Dirac points in one BZ and an experimentally measured winding\nnumber of 2. As summarized in Table R1, extending the coupling range systematically\nincreases both the number of Dirac points and the associated winding numbers in a graphene\nlattice. Specifically, introducing 8th NN couplings is predicted to produce up to 32 Dirac points\nwith a winding number of 3, and 11th NN couplings up to 50 Dirac points with a winding\nnumber of 4. To be more concrete, we performed experiments to demonstrate our approach.\nWe constructed an acoustic graphene lattice with 8th NN coupling. The measured dispersions\nare shown in Fig. R1, together with high winding number of 3 experimentally identified. The\nsimulated eigenfield distributions of the sixfold edge states—equally distributed at the upper\nand lower edges—are further provided in Fig. R2.\nIn Refs. Nat. Commun. 15, 3133 (2024) and Science 361, 782 (2018), multiple Dirac\npoints arise from the trigonal warping effect in a high-purity bilayer graphene or from Dirac\nquasiperiodic order, both posing substantial materials/structural challenges that hinder direct\nobservation of associated topological edge phenomena. In contrast, our long-range-coupled\nacoustic platform is simple, symmetry-preserving, and tunable. We measure the satellite Dirac\ncones symmetrically orbiting a central one and experimentally identified the high winding\nnumbers. With longer-range couplings, we systematically scale up the number of Dirac points\n(as summarized in Table R1).\nIn Ref. Adv. Photon. 7, 046005 (2025), high winding numbers and the related zero modes\nare constructed by displacing Dirac points via 𝐶 -symmetry breaking, which greatly limits its\n3\nability to achieve higher winding numbers. Besides, the existence of the zero modes strictly\nrequires specific twig-type boundaries. In contrast, our method preserves 𝐶 -symmetry and\n3\ndoes not rely on special boundaries, enabling higher winding numbers and more mid-gap edge\nstates in the presence of longer-range couplings.\nIn Ref. Laser Photon. Rev. 18, 2300686 (2024), multiple valley interface states are\nobtained by increasing the number of scatterers in each photonic unit cell, which inevitably\ninvolves high-order modes and narrows the bandwidth. In contrast, our interface states originate\nfrom satellite valleys generated by long-range couplings, avoiding higher-order modes and\noffering a wider frequency window. Furthermore, extending the coupling range can readily\nfurther increase the number of interface channels while maintaining nearly the same bandwidth\n[as exemplified in Fig. R3(c)].\nIn summary, we establish a broadly applicable and highly tunable framework, which\nrequires only long-range couplings, to realize and scale multiple Dirac points and large winding\nnumbers in a symmetry-preserving and easy-to-implement platform. We hope the added data\nand the above clarifications address the reviewer’s concerns about novelty and general interest.\nThe above discussion has been added to the revised manuscript, and the previously missing\nreferences have been added accordingly. A new section has been added to the Supplementary\nMaterials.\nCoupling range Maximal number of Dirac points Largest winding number\n3rd NN 8 2\n4th NN 14 2\n7th NN 14 3\n8th NN 32 3\n9th NN 26 3\n11th NN 50 4\nTable R1. Maximum number of Dirac points and largest winding number achievable in a\nchiral-symmetric graphene lattice with 𝑛th nearest-neighbor (NN) couplings.\nFig. R1. Dirac cones and winding numbers for an acoustic lattice with 8th NN couplings. (a)\nSchematic of the acoustic unit cell. (b) Photograph of the sample. (c) Measured bulk spectrum\nnear the Dirac point frequency. The hexagon outlines the first BZ. (d) and (e) Measured\n(colormaps) and simulated (dotted lines) dispersion curves along the solid horizontal lines in\n(c). (f)-(j) Measured pseudospin trajectories (color points) along the five dashed lines in (c).\nThe gray solid curves are obtained from simulation.\nFig. R2. Eigenfields of acoustic edge states in the lattice with 8th NN couplings. (a) Projected\ndispersions for the acoustic lattice with periodic boundaries along the 𝑥 direction and 28 unit-\ncells along the 𝑦 direction. The green and blue dashed lines mark 𝑘 = 0 and 𝑘 = 0.13 ,\n𝑥 𝑥\nrespectively. (b) Bulk winding number for different 𝑘 .\n𝑥\n(c) and (d) Eigenfrequency spectra at\n𝑘 = 0 and 𝑘 = 0.13.\n𝑥 𝑥\n(e) Corresponding eigenfield distributions for the edge states.\nFig. R3. (a) Schematic illustrating the origin of valley interface states in the main text. (b) Berry\ncurvatures and valley interface states from a lattice model with 3rd NN couplings. (c) Same to\n(b) for 8th NN couplings. The parameters are 𝑡 = 𝑡 = 5𝑡 and 𝑚 = 𝑡 , where 𝑡 represent the\n3 8 1 1 𝑖\ncoupling strength of the 𝑖th NN couplings and 𝑚 denotes onsite potential.\nA technical issue: the origin of the valley interface states is not thoroughly discussed. Can they\nbe characterized by the valley Chern number? Or can they be obtained from the low-energy\ntheory? I notice that the projections of the valleys are quite close. Will this cause any potential\ntrouble or effect?\nReply: We thank Reviewer #1 for the questions and apologize for not explaining the origin of\nthe valley interface states more clearly in the original manuscript. Our system follows the\nestablished valley-contrasting mechanism, even when the projections of the valleys become\ncloser. Nontrivial valley phases are characterized by valley Chern numbers, which, via bulk-\nedge correspondence within a low-energy effective theory, predict the valley interface states.\nA valley-specific low-energy Hamiltonian can be constructed at each valley to compute the\nvalley Chern numbers (approaching half-integer quantization as the gap closes) and to infer the\ncorresponding interface modes. In Fig. R4, we present numerical valley Chern numbers for a\nlattice model with 3rd NN couplings and on-site potential bias 𝑚.\nIn practice, closer valley projections are accompanied by relatively more delocalized Berry\ncurvature distributions in the BZ, which manifests in the valley interface states as their\ndispersion not well traverse the band gap, and meanwhile facilitates intervalley transitions of\nthe valley interface states. In our work, we exploit these intervalley transitions to enable\ntransmission through channels that are otherwise forbidden by valley locking. As shown in Fig.\nR5, for a small bias 𝛿 = 0.02 ℎ , each band of valley interface state traverses well the gap [Fig.\nconventional four-channel splitting consistent with valley locking [Fig. R5(c)]. For a larger bias\n𝛿 = 0.05 ℎ , the Berry curvature becomes more diffuse [Fig. R5(e)], and forward transmission,\n0\nwhich is previously suppressed by valley locking, appears at 4.81 kHz [Fig. R5(f)].\nThe above discussion has been added to the revised manuscript and the Supplementary\nMaterials.\nFig. R4. Valley Chern numbers for a model with 3rd NN couplings and 𝑡 = 3𝑡 . (a) Schematic\n3 1\nof integration regions (shading boxes) used to compute the valley Chern numbers at the 𝑞′ and\n1\n𝐾 valleys. (b) and (c) Valley Chern number versus on-site potential bias 𝑚 at 𝑞′ and 𝐾 ,\n1\nrespectively.\nFig. R5. (a) Projected band structure of an acoustic lattice in strip geometry with a bearded\ninterface and on-site bias 𝛿 = 0.02 ℎ . The orange circles denote simulated results, agreeing\n0\nwell with the dispersions from lattice model results (green lines) with the fitting parameters\n𝑡 = 132 Hz , 𝑡 = 50 Hz and mass term 𝑚 = 25 Hz . (b) Berry curvature distribution\n3 1\ncalculated from lattice model. (c) Acoustic pressure field distribution in a four-channel beam\nsplitter with frequency domain excitation at 4.81 kHz . The red star denotes the excitation\nsource. (d) Similar as (a) for 𝛿 = 0.05 ℎ and 𝑚 = 100 Hz. (e) and (f) Corresponding Berry\n0\ncurvature and field distribution for 𝛿 = 0.05 ℎ .\nThis work studies Satellite Dirac cones by simplifying the bilayer configuration into a\nmonolayer one, which is realized by introducing third-nearest neighbor couplings. By using\nphononic crystals for acoustic waves, they experimentally observe the satellite Dirac cones, and\nidentify related properties like the acoustic edge states determined by distinct non-trivial\nwinding numbers. This is an interesting work, trying to avoid the pitfall of complex calculation\nin bilayer graphene. I will recommend this paper to be published after solving following\nconcerns.\nReply: We thank Reviewer #2 for the positive assessment and recommendation. We appreciate\nthe reviewer’s recognition of our approach, that is, simplifying the bilayer configuration into\nan effective monolayer via longer-range couplings. We agree that this strategy offers a\ntransparent alternative to the computational complexity of bilayer graphene models. Below, we\naddress each concern point-by-point to further strengthen the manuscript, and the manuscript\nhas been accordingly revised.\nIn Fig.2, are the acoustic dispersions gapless? please give the partial enlarged graph around K\npoint and explain the reason.\nReply: We thank Reviewer #2 for this question. In response, we have included a partially\nenlarged view around the 𝐾 point (see Fig. R6). Experimentally, the acoustic branches are\nexpected to be gapless (protected by crystalline symmetry of 𝐷 ); however, due to the finite\n3\nsize of the phononic crystal sample, achieving a perfectly gapless acoustic bulk dispersion is\nchallenging. A larger sample would provide more spatial sampling points and thus higher\nFourier transformation resolution, which would more clearly reveal the crossing at 𝐾.\nWe have added the above discussion to the revised manuscript.\nFig. R6. Acoustic dispersion [Fig. 2(c) of the main text] (a), and an enlarged view (b) of the\nregion enclosed by the gray rectangle.\nIn Fig.3d, the acoustic eigenfields are different with a same k_x. Why?\nReply: We thank Reviewer #2 for this question. At 𝑘 = 0.25, a bulk winding number of 2\n𝑥\nimplies four edge eigenfields, two localized at upper boundary and two at the lower boundary.\nBecause these eigenfields are mutually orthogonal, their acoustic field distribution differ, as\nexpected. In Fig. 3(d) of the main text, we present the two simulated eigenfields near the upper\nboundary at 𝑘 = 0.25.\n𝑥\nWhy the authors employ the Bernal-stacked bilayer graphene but not the turbostratic stacked\n(AA) bilayer graphene? Is the developed method also applicable to the latter?\nReply: We thank Reviewer #2 for these insightful questions. The Bernal-stacked bilayer\ngraphene is a system with multiple Dirac points in low energy [e.g., Nat. Commun. 15, 3133\n(2024)], and therefore is employed to compare with our model that has multiple Dirac points\ninduced by long-range couplings. For the AA-stacked bilayer graphene, the effective\nHamiltonian reads,\nℎ = 𝜎 ⊗(𝑝 𝜎 +𝑝 𝜎 )+𝑡𝜎 ⊗𝜎 , (R1)\n𝐴𝐴 0 𝑥 𝑥 𝑦 𝑦 𝑥 0\nwhere 𝜎 are Pauli matrices and 𝑡 denote interlayer AA couplings. This Hamiltonian is known\nto host a nodal ring at zero energy, in contrast to multiple point-degeneracies in Bernal-stacked\nbilayer graphene or in our study. Because the four bands generated by Eq. R1 can not be\npartitioned into a separated pair of bands near zero energy, the effective two-band reduction\nprocedure is in principle invalid. (By contrast, in Bernal-stacked bilayer graphene the two inner\nbands are spectrally separable from the others near zero energy.) If one nevertheless enforces a\ntwo-band reduction, one obtains (𝑡−𝑝)2[(𝑝 −𝑖𝑝 )𝜎 +(𝑝 +𝑖𝑝 )𝜎 ], which indeed yields\n𝑥 𝑦 𝑥 𝑥 𝑦 𝑦\na nodal ring but also produces a Dirac cone. The latter should not appear in the low-energy\nsolutions of Eq. R1. Overall, since the AA-stacked bilayer graphene departs from the case of\nmultiple Dirac cones, we exclude the AA stacking configuration and compare only with the\nBernal-stacked bilayer graphene.\nThe above discussions have been added to Supplementary Note 1.\nHow to quantitatively demonstrate the equivalence between BBG and the TNN model? From\nFig. S1, they are not perfectly equivalent. Will this difference lead to distinct properties in Dirac\ncones? Furthermore, can this simplified TNN model replace Bilayer models in other physical\nphenomenon? Or they are only equivalent in Dirac cones prediction?\nReply: We thank Reviewer #2 for these insightful questions. At low energy, the effective\nHamiltonian of the TNN model near valley 𝐾 is given by\n2\n0 𝑣 𝜋+𝑣 (𝜋†)\nℎ = [ 1 2 ], (R2)\n𝑇𝑁𝑁\n𝑣 𝜋†+𝑣 𝜋2 0\n1 2\nwhere 𝜋 = 𝑝 +𝑖𝑝 , 𝑣 = √3(2𝑡 −𝑡 )⁄2 , and 𝑣 = (𝑡 +4𝑡 )⁄8 with 𝑡 and 𝑡 being the\n𝑥 𝑦 1 3 1 2 1 3 1 3\nnearest- and 3rd-neighbor-nearest couplings, respectively. For nonzero 𝑡 , the linear terms 𝑣 𝜋\n3 1\nand 𝑣 𝜋† vanish at 𝑡 = 𝑡 ⁄2 , yielding a purely quadratic Dirac cone. Away from this\n1 3 1\ncondition, the spectrum exhibits four linear Dirac cones with a central one orbiting by three\nsatellite cones. Noted that Eq. R2 is identical to the well-known effective Hamiltonian of BBG\nwith trigonal warping [e.g., Rep. Prog. Phys. 76, 056503 (2013)], implying that the low-energy\ndispersions and Berry-curvature distributions of the TNN model and BBG coincide, as\nconfirmed numerically in Fig. R7. We apologize for any confusion caused by Fig. S1. In Fig.\nS1(c), the spectrum was calculated from the TNN Hamiltonian with 𝑡 = 𝑡 ⁄2, while in Fig.\n3 1\nS1(d), it was calculated from the same Hamiltonian with 𝑡 = 𝑡 ⁄4. As sharing the common\n3 1\neffective Hamiltonian in Eq. R2, the simplified single-layer TNN model can replace the bilayer\nmodels at low energy. Consequently, similar phenomena can be realized in both systems such\nas engineering complex Landau levels via tunable Dirac material [Nat. Commun. 15, 3133\n(2024)]. Moreover, since the TNN model is of a simple single-layer structure, it offers greater\nflexibility for exploring effects associated with larger topological numbers.\nWe have added the above discussions in Note 1 of the Supplementary Materials.\nFig. R7. Band structures and Berry curvature distributions (color hotspots) near 𝐾 point for the\nTNN model and Bernal-stacked bilayer graphene. In the TNN model, 𝑡 and 𝑡 denote the\n1 3\nnearest and the 3rd-nearest-neighbor couplings. In bilayer graphene, 𝑡 and 𝑡 denote the\n⊥ 𝑤\nvertical and skew interlayer couplings, respectively.\nIn the supplementary information Note 6, please give the specific value of the group velocity\nV_g and other physical quantities, which can make the readers to reproduce relevant data. Also,\nthe author should also check the main text and give the specific value of related physical\nquantities.\nReply: We apologize for the ambiguity in Supplementary Note 6 and thank Reviewer #2 for\npointing this out. Following Ref. Phys. Rev. Appl. 21.034026 (2024), the transmission ratio at\nport 𝑗 is defined as the ratio between the output and the input acoustic energy:\n∫ 𝐼 𝑑𝑣/∫ 𝐼 𝑑𝑣 , where 𝐼 = |𝑝|2/𝜌 𝑣 . The term 𝜌 𝑣 is a constant existing in both\n𝑉𝑗 𝑗 𝑉in in air g air g\ndenominator and numerator and canceled out, leaving the integrand only the square of acoustic\npressure field |𝑝|2 . Figure R8 shows the geometry used to calculate the transmissions in\nSupplementary Note 6, the integration area 𝑉 and 𝑉 are manually selected (blue regions).\nin 𝑗\nTwo identical point sources are placed at the input channel with a sweeping step of 5 Hz. Global\ndamping is included via a complex sound speed (344+0.3 𝑖) m/s. To suppress reflections,\ncladding layers with stepped damping 344+𝑙 0.8 𝑖 (𝑙 = 1,2,…11) are attached to each port.\nWe have revised Supplementary Note 6 to incorporate these clarifications.\nFig. R8. Geometry used for transmission calculation. Blue regions indicate the integration areas\n𝑉. Red stars mark the point sources. Gray arrows indicate the increasing direction of stepping\n𝑗\ndamping.\nHow to accurately define this word “neighbor” ? In one periodic cell, there are six carbon nodes\n(if numbered 1-6 in anticlockwise direction). For node 1, it is related to node 2 and node 6 by\nconventional interaction and related to node 4 by the third neighbor effect. Then what about\nnode 3 and node 5? Are they negligible?\nReply: We thank Reviewer #2 for these important questions. The word “neighbor” in this work\ndefined by direct geometric distance between atomic (cavity) nodes, as schematically illustrated\nin Fig. R9. For an arbitrary carbon node, there are three nearest neighbors (NN), six second-\nnearest neighbors (2NN), three third-nearest neighbors (3NN), six fourth-nearest neighbors\n(4NN), etc. In our phononic crystals, couplings arise only through physical connection tubes\nbetween cavities. Since cavities at nodes 3 and 5 are not connected to the cavity at node 1 by\ntubes in the unit cell, the couplings between them are negligible and can be safely omitted.\nFig. R9. Schematic of the 𝑛th-nearest neighbor (nNN) up to 8NN for then node labeled by 1.\nSome additional minor questions:\nIn Fig.1d, the author obtains the projected band spectrum through the path “q-K-q’-q-K’-q”,\nbut I do not find any illustration about the q and q’ point. please explain it or mark them in the\nfigure.\nReply: We thank Reviewer #2 for this careful remark. We apologize for omitting the bar\nnotation for the projections of Dirac points. Consistent with the conventions in Figs. 3 and 4 of\nthe main text, we denote the projections of 𝐾 and 𝐾′ as by 𝐾̅ and 𝐾̅′, respectively. Likewise,\nthe projections of 𝑞 and 𝑞′ are labeled 𝑞̅ and 𝑞̅′, respectively. We have updated Fig. 1(d) and\nassociated Fig. 1(e) to include these labels.\nIn Fig.4d and Fig.4e, I would recommend the author use a same color bar to facilitate the\ncomparison between simulation and experimental results.\nReply: We thank Reviewer #2 for this helpful suggestion. We have revised Fig. 4d and Fig. 4e\nin the main text so that they share the same color bar to facilitate direct comparison between\nsimulation and experiment.\nIn recent years, the magic-twisted bilayer graphene is a research hotspot. Can the researchers\ndo some related discussions? Note that twisting may generate flat band, which will cause some\nexotic phenomenon. Is this TNN model applicable to flat band or superconductivity study?\nReply: We thank Reviewer #2 for these insightful questions. The low-energy bands of twisted\nbilayer graphene near the magic angle can be significantly compressed, forming nearly flat\nband dispersions. In Ref. Nature 556, 80 (2018), the evolution of Dirac cones in twisted bilayer\ngraphene versus the twisting angle close to a magic angle is calculated with the results shown\nin the left panel of Fig. R10, where the color shows the hotspots of Berry curvature for each\nband. In our TNN model, by tuning the 3rd-nearest neighbor coupling, a similar evolution of\nDirac cones and Berry curvature distributions are identified, as shown in the right panel of Fig.\nR10. This suggests that our TNN model can effectively capture the low-energy bands profiles\nand topological characteristics of twisted bilayer graphene near the magic angle.\nFurthermore, as highlighted in our work and Ref. Nat. Rev. Phys. 7, 299 (2025), the long-\nrange couplings can induce roton-like effect, reflected in the reversion of the band slopes.\nThis implies that the incorporation of longer-ranged couplings could revert the band slopes multiple\ntimes and potentially squeeze the bands into a narrow energy range, as demonstrated in Fig.\nR11. Such behavior may suggest a unique and fruitful platform for exploring flat bands.\nHowever, we note that our TNN model is formulated within a single-particle framework and,\nas such, does not capture many-body correlation effects necessary to phenomena such as\nsuperconductivity. Incorporating such effects would require more extensions.\nWe have added the above discussions in the revised manuscript.\nFig. R10. Left panel [adapted from Nature 556, 80 (2018)]: Evolution of the low-energy band\nstructures and the Berry curvature in twisted bilayer graphene. Right panel: Evolution of the\nBerry curvature as a function of 3rd-neighbor coupling 𝑡 in the TNN graphene model.\n9th NN couplings, and (c) simultaneous 3rd NN and 9th NN couplings."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the authors for their comprehensive explanation regarding the novelty of this work and the additional numerical and experimental results. I am now convinced that this paper can meet the criteria of Nature Communications. However, before making a recommendation, I would like to invite the authors to further address the technical issue raised in my original report.\nIn the revised version, the authors show that the interface states are characterized by the valley Chern number and follow a bulk-boundary correspondence similar to the standard valley interface states. Then I think it is necessary to check that if these interface states can go through sharp bends (like 60-degree or 120-degree path bends) or not, as this is the most important property of valley transport in phononic/photonic valley systems.", "reviewer_2": "The authors have properly replied to my comments. Can be published."}, "author_response": "RESPONSE TO REVIEWER #1:\nI thank the authors for their comprehensive explanation regarding the novelty of this work and\nthe additional numerical and experimental results. I am now convinced that this paper can meet\nthe criteria of Nature Communications. However, before making a recommendation, I would\nlike to invite the authors to further address the technical issue raised in my original report.\nReply: We thank Reviewer #1 for the positive comments on our manuscript and are pleased\nthat our work is now considered to meet the criteria for publication in Nature Communications.\nIn the revised version, the authors show that the interface states are characterized by the valley\nChern number and follow a bulk-boundary correspondence similar to the standard valley\ninterface states. Then I think it is necessary to check that if these interface states can go through\nsharp bends (like 60-degree or 120-degree path bends) or not, as this is the most important\nproperty of valley transport in phononic/photonic valley systems.\nReply: We sincerely thank Reviewer #1 for their valuable suggestion. We have conducted\nadditional simulations to examine the robustness of the interface states when traversing sharp\nbends, such as 60-degree and 120-degree path bends. The results, shown in Fig. R1, confirm\nthat the valley interface states can robustly propagate through these bends without\nbackscattering, thereby demonstrating the key property of valley transport.\nThe above discussion and Fig. R1 have been added to revised manuscript and\nSupplementary Information.\nFig. R1. Simulation for the valley interface states propagating through (a) 60-degree and (b)\n120-degree path bends. The stars indicate the locations of the point sources.\nRESPONSE TO REVIEWER #2:\nThe authors have properly replied to my comments. Can be published.\nReply: We sincerely thank Reviewer #2 for the recommendation to publish our work."}]

Decision letter