Negative Refraction of Weyl Phonons at Twin Quartz Interfaces

In Nature, α-quartz crystals frequently form contact twins, which are two adjacent crystals with the same chemical structure but different crystallographic orientation, sharing a common lattice plane. As α-quartz crystallizes in a chiral space group, such twinning can occur between enantiomorphs with the same handedness or with opposite handedness. Here, we use first-principles methods to investigate the effect of twinning and chirality on the bulk and surface phonon spectra, as well as on the topological properties of phonons in α-quartz. We demonstrate that, even though the dispersion appears identical for all twins along all high-symmetry lines and at all high-symmetry points in the Brillouin zone, the dispersions can be distinct at generic momenta for some twin structures. Furthermore, when the twinning occurs between different enantiomorphs, the charges of all Weyl nodal points flip, which leads to mirror symmetric isofrequency contours of the surface arcs on certain surfaces. We show that this allows negative refraction to occur at interfaces between certain twins of α-quartz.


I. INTRODUCTION
Negative refraction is a counter-intuitive phenomenon in which incident and refracted waves emerge on the same side of the interface normal 1 , providing potential applications in superlens and sub-wavelength imaging 2 .One prominent strategy to obtain negative refraction uses open isofrequency contours, which has been realised in hyperbolic metamaterials [3][4][5][6] .Recent advances in topological materials offer a new platform to manipulate isofrequency contours arising from topological surface states [7][8][9][10][11][12][13][14][15][16] .For example, the surface arcs of Weyl points can form distinct isofrequency contours for both positive and negative refraction 17,18 , and all-angle reflectionless negative refraction has been observed in Weyl metamaterials 19 .Intuitively, negative refraction can take place at the interface between chiral crystals with left-and righthanded screw symmetries, as the isofrequency contours of the surface Weyl arcs are mirror images of each other for a specific choice of orientation.
One of the most well-known chiral crystal is α-quartz, which exists in nature in two enantiomorphs that belong to a space-group pair and are thus handed: the righthanded screw P 3 1 21 (No. 152) and the left-handed screw P 3 2 21 (No. 154).The crystal structures have opposite chirality, which can be distinguished either by measuring their optical activity (the rotation of the plane of polarisation of plane-polarised light), as first observed in quartz crystals in 1811 by François Arago 20,21 , or by circularly polarised resonant x-ray diffraction [22][23][24] .As α-quartz is an insulator under ambient conditions, any potential electronic topology away from the Fermi level is not easily accessible.To remedy this, we instead consider the topology of the intrinsic lattice vibrations (phonons) in α-quartz, which are not constrained by the Fermi level.The topological properties of phonons have been studied extensively  . In cotrast to metamaterials, phonons are intrinsic quasiparticles in real materials similar to electrons.Additionally, typical phonon frequencies range in 0−50 THz, and negative refraction in the terahertz frequency range has been well studied experimentally 5 .Therefore, it is expected that negative refraction can be straightforwardly measured in topological phonons using existing apparatus.The phonon modes in the two enantiomorphs exhibit chiral behaviours such as opposite pseudo-angular momenta, selective optical transitions and opposite transport direction [50][51][52][53][54][55][56][57] .Furthermore, band crossings between modes in a single enantiomorph of α-quartz form Weyl points because of the lack of inversion symmetry in the space groups P 3 1 21 or P 3 2 21 35 , and it is well known that the Weyl points of two enantiomorphs carry opposite Chern number [58][59][60][61][62][63] .It is therefore expected that the Weyl phonons in α-quartz with left-and right-handed screw carry opposite Chern number.
In nature, quartz naturally forms contact twin structures -two crystals of the same mineral but different orientations, touching along a common plane.Twinned quartz crystals are much more common than untwinned ones on earth 64 , and quartz crystals are therefore generally racemic 65,66 .In twinned crystals, many different structures with various orientation of the unit cells are possible, and have been generally classified by twinning laws 67 .As such, twinned quartz crystals offer a natural and versatile platform to study the interplay between chirality and topology, as the twinning boundaries of αquartz should be easily accessible.It is worth mentioning that the relationship between chirality and topology is a very active area of research 61,[68][69][70][71][72][73][74][75][76][77][78] , as chiral space groups can host a plethora of interesting topological phe-nomena.Hence, a careful comparison of enantiomorphic structures is of significant current interest, with the potential application of negative refraction occurring at the twin boundary in α-quartz.
In this work, we explore the relationship between the twinning type of the chiral crystal structures, their phonon band structure and their associated Weyl points.We show that, for the three most common types of quartz contact twinning, the bulk phonon band structures coincide along all high-symmetry lines in the Brillouin zone.However, depending on the twinning choice, the band structures differ at generic momenta.Furthermore, even when the bulk band structure agrees for enantiomorphic twins, the surface isofrequency contours differ.This allows negative refraction to occur at the twinning interface between two enantiomorphs.

II. METHODOLOGY
Density functional theory (DFT) calculations are carried out using the Vienna ab-initio simulation package (vasp) 79,80 .The generalised gradient approximation (GGA) calculations are performed using the Perdew-Burke-Ernzerhof exchange-correlation functional as revised for solids (PBEsol) 81 , with 4 valence electrons (3s 2 3p 2 ) for silicon atoms and 6 valence electrons (2s 2 2p 4 ) for oxygen atoms.The plane-wave basis set has an upper kinetic energy limit of 800 eV and the k-mesh has a size of 7 × 7 × 7, with the self-consistent field loop stopped when energy differences between steps are below 10 −6 eV.Structural relaxations are carried out until the Hellman-Feynman forces are less than 10 −2 eV/ Å.
Density functional perturbation theory is used in calculating Hessian matrices and phonon frequencies 82,83 , implemented on a 3×3×2 supercell with a 3×3×3 k-mesh.phonopy is used to build the matrix of force constants, diagonalise the dynamical matrix, and obtain the phonon dispersion curves 84,85 .Convergence of the calculations is assessed by varying both the supercell and k-mesh sizes and noting no discrepant results.The calculations include the splitting between transverse and longitudinal optical phonon modes (LO-TO splitting) 86 , but we note that LO-TO splitting plays a minor role on the topological properties of phonons away from the Brillouin centre.WannierTools 87 is used to locate every single band crossing point on a phonon q-mesh of size 51×51×51, to calculate the chiralities of the Weyl nodes, and to compute the phonon surface states via the surface Green's function.

III. RESULTS AND DISCUSSION
A. Crystal structures and lattice dynamics α-quartz is the most stable phase of silica under ambient conditions, and crystallises in the trigonal crystal sys- .455Å, which are in good agreement with previous measured and calculated data 24,[88][89][90] .We first focus on the conventional unit cell choice, where the (x,y) position of all atoms in the unit cell agree, and the difference between the enantiomorphs is solely determined by the relative atomic z coordinates, where the enantiomorphs are related by a mirror symmetry, m z , with respect to the z = 0 plane.This corresponds to Leydolt twinning, as explored below.Using this unit cell, we show the phonon spectra for α-quartz with both space groups in Fig. 1(b), agreeing well with previous calculations and measurements 35,[90][91][92][93] .No imaginary mode is found in the entire Brillouin zone, indicating their dynamic stability.For this choice of twinning, we find that the phonon dispersions of the two chiral structures agrees along all high-symmetry momenta and lines but not at general momenta as shown in Fig. 2(a) and explored below.The phonon dispersions for other twinning choices are shown in the Supporting Information.

B. Weyl points
We focus on the band crossing points formed by phonon branches 18 and 19 in the frequency range 20.5−24.0THz, as they are relatively isolated from the other bands and show the most interesting topological features.On the q z = 0 plane, no band crossing points are formed between these two phonon branches as the two bands are far away from each other.In contrast, the two phonon branches tend to touch on the q z = 0.5 plane (in units of the reciprocal lattice vector 2π/c).
Along the A-L line, the two bands form an avoided crossing (i.e. they are gapped), as shown in Fig. 2(a).The point-group along this line is 2 , so there is no non- trivial unitary point-group symmetry, and hence there is only one irreducible representation (IRREP), GP 1 , to which both bands belong, as shown in Fig. 2(a).On the other hand, the point group of the high-symmetry line Q = H-A is 2, which gives rise to two different IRREPs Q 1 and Q 2 .We find that bands 18 and 19 belong to different IRREPs along Q and as such, these bands form a stable Weyl point.By threefold rotation and time-reversal symmetry, there are thus a total of six Weyl points on the q z = ±0.5 plane, all carrying the same Chern number C within the same space group.However, for different space groups, the charges C of the Weyl points on the q z = ±0.5 plane are opposite, i.e.C = −1 for P 3 1 21 α-quartz and C = +1 for P 3 2 21 α-quartz respectively.
In addition to the Weyl points on the q z = ±0.5 plane, there are also six Weyl points at generic momenta: C = +1 for P The positions of all the Weyl points in the Brillouin zone are demonstrated in Fig. 2(b).The six Weyl points at general q in a single enantiomorph are related to each other by the threefold (non-symmorphic) rotation symmetries and time-reversal, neither of which flip the chirality.

C. Surface states
The surface states of the Weyl points projected along the [010] direction are shown in Fig. 6(a), for the same choice of unit cell as in Fig. 1(a).The surface local densities of states (LDOS) is calculated from the imaginary part of the surface Green's function 87 .The projections of the bulk Weyl points are connected via surface arcs.The surface arcs along the high-symmetry lines M-Γ-Ā-L are exactly the same for both enantiomorphs for this choice of twinning, whereas those at generic momenta along the L-Γ are different from each other.For the surface states for other twinning choices, see the Supporting Information.
To have a better view of the distribution of the surface arcs in the reciprocal space, the isofrequency surface arcs at 22.1 THz are plotted in Fig. 6(b).The surface arcs at a fixed frequency, with the choice of unit cells shown in Fig. 1(a) for P 3 1 21 and P 3 2 21 α-quartz, are related by a reflection symmetry along q c .This is analysed further below.

D. Phonon dispersion and twinning
The space groups P 3 1 21 (No.152) and P 3 2 21 (No.154) are chiral and form an enantiomorphic pair 94,95 .When the same compound crystallises in two enantiomorphic structures, the unit cells of these structures are related by an operation of the second kind (changing handedness, e.g.mirror, inversion, rotoinversion or glide symmetries).
We will be interested in the interface between these two enantiomorphic structures.We assume that they touch at a crystal plane, e.g. that they from a contact twin along a contact plane.The study of such twins has a long history in mineraology.One distinguishes between merohedral twins (where the lattices are parallel) and non-merohedral twins (where the lattices are not parallel) 96 .We will only discuss merohedral twinnig here.The three most important merohedral twins of quartz are Dauphiné, Brazil, and Leyoldt (also known as combinedlaw or Liebisch) twins.Of these, Dauphiné twins occur as twinning between two crystals with the same handedness (same enantiomorphs), whereas Brazil and Leydolt twinning occurs between different enantiomorphs.Each of these is characterised by a different twinning operation.Dauphiné twinning occurs between two crystals whose crystallographic axis are related by a C 2z symmetry (twofold rotation around c-axis), Brazil twinning occurs between crystals related by P (inversion through the origin), whereas Leydolt twinning occurs between crystals related by m z = PC 2z (mirroring in the plane normal to the c-axis).In nature, Dauphiné and Brazil twinnings are common, whereas Leydolt twinning is rare 96 .
The effect of these symmetry operations on the bulk phonon dispersion is given by: Dauphiné : ω(q x , q y , q z ) C2z −−→ ω(−q x , −q y , q z ) (1a) Brazil : ω(q x , q y , q z ) P − → ω(−q x , −q y , −q z ) (1b) Leydolt : ω(q x , q y , q z ) mz − − → ω(q x , q y , −q z ).(1c) For the standard unit cell choice in Fig. 1(a), the crystal structures are related by Leydolt twinning which explains the relative positions of the Weyl points between the enantiomorphs.As the Berry curvature behaves as a pseudovector, the sign of the Chern number reverses under mirroring.This will also reverse the propagation direction of the topological surface states.We next analyse how the bulk and surface phonon band structures behave under twinning.

Twinned bulk band structures
By definition, none of these twinning operations are symmetries of the unitary part of the space groups.However, in non-magnetic systems, time-reversal symmetry T additionally enforces ω(q) = ω(−q).From this, we directly conclude that crystals related by Brazil twinning display the same dispersion for all q.
The same is not true for Dauphiné and Leydolt twinning.We therefore expect the band structure for this twinning to disagree at generic points.To understand the structure at high-symmetry points, we start by noting that both Leydolt twinning and Dauphiné twinning lead to the same constraints on the dispersion, as by timereversal symmetry, ω(−q x , −q y , q z ) = ω(q x , q y , −q z ).We therefore consider only Leydolt twinning in what follows.We see immediately that on the q z = 0 and q z = 0.5 planes, all twinning types will display the same dispersion.Thus, the only high-symmetry lines left to discuss are ∆ =A -Γ-A and P = H -K-H. Along both of these lines, 2 110 symmetry (twofold rotation with respect to the [110] axis) guarantees that the dispersion is symmetric around Γ or K respectively, and therefore flipping q z → −q z does not change the dispersion.Thus, the band structures look identical at all high-symmetry points and along all high-symmetry lines for all twins, as confirmed in Fig. S1 in the Supporting Information.This is not true, however, at generic momenta, as shown in Fig. 2(a).We note in passing that the lower-symmetry chiral space-group pairs P 3 1 and P 3 2 do not have 2 110 rotation symmetry, so that in crystals belonging to this chiral space-group pair (twinned by the same operations) we would expect the dispersions to disagree even along the high-symmetry line P (∆ remains symmetric around Γ by time-reversal symmetry).

Twinned surface states
The relationship between the surface states for different crystal twins is more subtle, as this depends on both the surface symmetries and on the crystal termination.We focus on the [010] surface with termination for Leydolt twinned crystals in Fig. 6.We find, as shown in Fig. 6(a), that the surface band structures agree along the high-symmetry lines, but differs at generic points.This may be a result of remnant time-reversal symmetry, or of the antiunitary symmetry 2 100 = 2 100 T , which leaves q y invariant.We show the surface band structures for the other twinning operations in Fig. S2 in the Supporting Information.More generally, it is not easy to predict the surface symmetries, as these will generally also be influenced by surface termination and reconstruction.
The surface arcs shown in Fig. 6(b) however, arise as a consequence of the bulk Weyl points, and should therefore display some topological protection.We expect these surface states to curve oppositely in the different enantiomorphs, as they are related by a handedness-inverting symmetry, flipping the chirality.This should be the case for both Leydolt and Brazil twinning, as inversion symmetry is generically broken at the surface.Such an inversion of direction in the surface arcs can give rise to negative refraction at the surface of twins.

E. Negative refraction
Benefiting from the isofrequency contours of P 3 1 21 and P 3 2 21 α-quartz, negative refraction can take place at the interface between these enantiomorphs.We focus on the [010] surface for Leydolt-twinned quartz [i.e.crystal structures shown in Fig. 1(a)].
We plot the schematics of negative refraction of surface Weyl arcs along the [010] direction that takes place at the interface between the two enantiomorphs in Fig. 4. For surface phonons at 22.1 THz in P 3 2 21 α-quartz with wave vector q 1 and incidence angle θ 1 , the tangential wave vector is conserved 5 , i.e. q 1 sinθ 1 = q 2 sinθ 2 (where q 2 and θ 2 are the wave vector and refraction angle in P 3 1 21 α-quartz), exhibiting positive refraction for the wave vector.However, the Poynting vector S, which is normal to the isofrequency surface arcs and directed along the energy flow, can exhibit negative refraction, as demonstrated by the Poynting vectors S 1 and S 2 with incidence and refraction angles of ϕ 1 and ϕ 2 respectively.As a result of negative refraction, the Poynting vectors S 1 and S 2 of the surface phonons are on the same side of the interface normal.Most interestingly, the negative refraction is tunable by varying the surface phonon frequency (as shown in Fig. S3 the Supporting Information).Such negative refraction can transform the linear interface between P 3 2 21 and P 3 1 21 α-quartz into a lens capable of focusing/defocusing phonon waves depending on the frequency-tunable incidence and refraction angles.
In terms of feasibility for synthesising such interfaces, we note that the [010] surface for twinned quartz can be easily found in nature, as quartz twin crystals are much more abundant than the untwinned ones 64 .In virtually every natural quartz, the two morphologically distinct natural crystals are internally twinned, i.e., the enantiomorphs can be found within the same crystal with an interface 65,66 .Moreover, the interfaces are experimentally realisable, as it has been reported that an atomically sharp internal interface between two enantiomorphs can be synthesised 97 .
In terms of measuring negative refraction, previous experiments have used an illumination frequency of about 27 THz by fabricating a gold antenna on one side of the interface to measure the propagation wave 5 , which is similar to the frequency in our work.The refractive behaviour can be measured by a tunable quantum cascade laser in a scattering-type scanning near-field optical microscope 98,99 .
In terms of potential applications, negative refraction in the terahertz region holds significant implications for thermal emission by controlling the flow of thermal energy.This discovery can also facilitate the development of thermal imaging techniques for medical imaging, aerospace and manufacturing.

IV. CONCLUSION
We explore the relationship between chiral crystal structures, twinning and topological charges of Weyl points.We find that the, depending on the choice of twinning operation, the bulk band structure for different twins can agree or disagree at generic points.We find further that the Weyl points in opposite chiral structures carry opposite Chern number, which can lead to negative refraction at the twinning surface between different enantiomorphs.

FIG. 1 .
FIG. 1.(a) Crystal structures of α-quartz with space group P 3121 (No. 152) and P 3221 (No. 154).(b) Phonon dispersion of α-quartz.The two structures are related by a mirror operation, leading to Leydolt twinning as explored below.

FIG. 3 .
FIG. 3. (a) Topological surface states along the [010] direction of P 3121 (No. 152) and P 3221 (No. 154) α-quartz with unit cell related by mz symmetry (Leydolt twinning), and (b) their isofrequency surface arcs at 22.1 THz.Note that the bulk Weyl points are at 22.5 THz so there is no projection of the bulk Weyl points in (b).

FIG. 4 .
FIG. 4. Schematics of negative refraction of surface Weyl arcs along the [010] direction at 22.1 THz, which takes place at the interface between P 3221 (No. 154) and P 3121 (No. 152) α-quartz, related by a mirroring operation.