Order-by-disorder from bond-dependent exchange and intensity signature of nodal quasiparticles in a honeycomb cobaltate

Recent theoretical proposals have argued that cobaltates with edge-sharing octahedral coordination can have significant bond-dependent exchange couplings thus offering a platform in 3d ions for such physics beyond the much-explored realisations in 4d and 5d materials. Here we present high-resolution inelastic neutron scattering data within the magnetically ordered phase of the stacked honeycomb magnet CoTiO3 revealing the presence of a finite energy gap and demonstrate that this implies the presence of bond-dependent anisotropic couplings. We also show through an extensive theoretical analysis that the gap further implies the existence of a quantum order-by-disorder mechanism that, in this material, crucially involves virtual crystal field fluctuations. Our data also provide an experimental observation of a universal winding of the scattering intensity in angular scans around linear band-touching points for both magnons and dispersive spin-orbit excitons, which is directly related to the non-trivial topology of the quasiparticle wavefunction in momentum space near nodal points.


INTRODUCTION
Spin-orbit coupling is at the origin of many remarkable properties of condensed matter uncovered in recent years [1][2][3][4][5]. It is central to the appearance of nontrivial topological invariants in electronic band structures and underlies the existence of bond-dependent exchange couplings that have been shown to bring about exotic features in many quantum magnets [6][7][8]. In the latter case much of the effort in materials discovery has focussed on heavy 5d and 4d ions in which the spin-orbit coupling is one of the dominant energy scales. Notable are the honeycomb iridates A 2 IrO 3 (A=Na,Li) and related materials, and α-RuCl 3 , which displayed a range of many novel exotic magnetic properties including spin-momentum locking [9], incommensurate orders with counter-rotating spin spirals [6], broad scattering continua in the spectrum of spin excitations [10] or unconventional field-dependent thermal Hall effect [11]. The origin of these exotic forms of behaviour is the presence of significant anisotropic, bond-dependent exchange, which in extreme cases has been predicted to stabilize quantum spin liquids, such as the celebrated Kitaev honeycomb model with Ising ex-changes along orthogonal directions for the three bonds that meet at each site [12]. The path to the discovery of the unusual magnetic properties of those materials has been a fruitful one starting with theoretical proposals that bond-dependent exchange couplings can arise in certain iridates and ruthenates with edge-sharing octahedra [13,14]. The octahedra supply a crystal field environment that leads to an effective low-energy spin one-half degree of freedom for the magnetic ions and the edge-sharing provides the local exchange pathway that, in conjunction with the spin-orbit coupling, produces anisotropic bond-dependent exchange. There is now evidence for significant such exchanges in honeycomb iridates and ruthenates [6][7][8].
More recent theoretical work has argued that significant bond-dependent exchange in the form of Kitaev and related couplings may also arise between Co 2+ ions in edge-sharing octahedral coordination [15][16][17] thus extending the original proposals into a surprising new setting. To investigate such effects we report here inelastic neutron scattering (INS) measurements of the spin dynamics in the stacked honeycomb magnet CoTiO 3 . Our data show propagating spin wave excitations with a clear low energy spectral gap, which was inferred but could not be resolved by previous studies [18]. We show that the spin wave spectrum is not merely compatible with the presence of bond-dependent exchange, but that such couplings must be present in the low energy pseudo-spin one-half theory in order to explain the origin of the gap. Moreover, we show that the gap opening must occur via a quantum order-by-disorder mechanism [19][20][21][22][23][24][25] as a consequence of unusually strong constraints on the possible mechanisms that can open the spectral gap. In view of the low-lying crystal field excitations in this material compared to the exchange coupling, we provide compelling evidence that virtual crystal field excitations are the driving mechanism for order-by-disorder [26,27] assisted by spin-orbital exchange and supply a calculation of the spin wave spectrum including this effect that captures the principal features of the data.
CoTiO 3 is part of a growing list of materials [18,28,29] explored as candidates displaying Dirac magnons. Earlier studies established the presence of Dirac nodal lines [18], which make this material ideal for the exploration of a recently predicted [30] fingerprint of a topologically non-trivial magnon band structure, namely a universal azimuthal modulation in the dynamical structure factor around linear band touching points, not probed exper-arXiv:2007.04199v2 [cond-mat.str-el] 26 Jun 2021 (11) (01) imentally before and which originates from the special topological features in the wavefunction of nodal quasiparticles. We indeed observe clear evidence for the predicted intensity winding around the nodal points, thus providing a direct measurement of the non-trivial topology of the Dirac magnon wavefunctions and establishing that there are meaningful features in the momentumand-energy dependent dynamical structure factor beyond simply revealing the quasiparticle dispersion relations. Furthermore, we observe analogous features in the dispersive spin-orbital excitations at higher energy, highlighting the universal properties of Dirac bosonic quasiparticles. Finally, we investigate the effect of the bonddependent exchange on the Dirac nodal lines arguing that they are robust to gap opening and likely appear as 'double helices' winding around each zone corner. We show that the same type of bond-dependent anisotropic exchange that opens up the spectral gap provides a natural explanation for a 'double-peak' structure in energy scans near the nodal points.

RESULTS
Magnon dispersions -The magnon dispersions along high-symmetry directions in the honeycomb plane obtained using inelastic neutron scattering (INS) measurements on single crystals of CoTiO 3 (for details see Supplementary Note 6A) are summarized in Fig. 1a). Wavevectors are indexed in reciprocal lattice units of the hexagonal structural unit cell. Near the (1,1,3/2) magnetic Bragg peak the lowest mode has a near-linear in-plane dispersion. As the honeycomb layers are ferromagnetically ordered with moments confined to the crystallographic ab plane, the linear dispersion indicates predominant easy-plane-type exchange couplings for inplane neighbors. Fig. 1c) observes a finite dispersion at low energies in the direction normal to the layers, in-dicating finite inter-layer couplings, and a small but finite spectral gap ∆ = 1.0(1) meV, clearly resolved above the magnetic Bragg peak. Ref. [18] proposed that a finite spin gap would be needed to account for the observed non-linear magnetization curve in small in-plane fields [31], but it was not possible to directly resolve the gap excitation in the earlier lower-resolution INS data [18]. Apart from the finite gap, the main features of the magnon spectrum can be accounted for by a minimal exchange Hamiltonian H XXZ for the stacked honeycomb geometry in CoTiO 3 , allowing for each bond a different exchange coupling between the moment components along the c-axis, and between the components in the ab plane. For a single ferromagnetic honeycomb layer, two magnon bands (acoustic/optic) would be expected with linear crossings at the corners (K-points) of the hexagonal Brillouin zone. For finite interlayer couplings that stabilize antiferromagnetic stacking of layers, the number of bands doubles and inter-layer resolved lower bands are expected with almost degenerate higher bands, as observed in Figs. 1a,c). H XXZ has a gapless (Goldstone) mode corresponding to moments rotating freely in the ab plane, so to capture the observed gap we assume that the physical mechanism responsible for gap generation only modifies the dispersion relations ω(k) of H XXZ by adding a gap in quadrature, i.e. experimental dispersion points are compared withω(k) = ω 2 (k) + ∆ 2 . We call this parameterization the XXZ∆ model to emphasize that the gap ∆ is not intrinsic, but is an additional, empirical fitting parameter. We find that exchanges up to 6th nearest-neighbor (nn) are important and obtain a very good level of agreement for both the dispersions and intensities as shown by comparing Figs. 1a) with b), and c) with d) (for more details see Supplementary Notes 5B, 5C, and 6C).
Quantum Order-by-Disorder -The presence of the finite magnon spectral gap ∆ is important as it indicates preferential moment orientations inside the easy plane. The magnetic ground state of Co 2+ (3d 7 ) ions in the local crystal field environment is a Kramers doublet with pseudospin-1/2, for which there is no local anisotropy, so any preferential orientation must be selected by interactions beyond the minimal H XXZ Hamiltonian. We focus our attention on bilinear couplings in the pseudospin as higher order two-site couplings project down to such couplings. As outlined in Supplementary Note 10, multi-site couplings will be suppressed by the large charge gap. As there is no detectable distortion of the crystal lattice following the onset of the magnetic order, we perform the analysis of bi-linear couplings between cobalt moments that are symmetry-allowed by the crystal structure space group. We find that whilst various bond-dependent exchange couplings can be present in principle, at the classical level, surprisingly, the ground state energy remains independent of the moment orientation in the ab planesee Supplementary Notes 7 and 8.
This degeneracy must however be an artefact of the mean-field approximation, as the real material Hamiltonian has only discrete, rather than continuous rotational symmetry around the c-axis. Such degeneracies would in general be expected to be lifted by quantum fluctuations via an order-by-disorder mechanism [19][20][21][22][23][24], when the ground state energy (per site) acquires a contribution from zero-point fluctuations of the form qu (φ) = 1 2 m ω m (k) , where φ defines the moments' orientation in the ab-plane relative to the a-axis and ω m (k) is the average energy of dispersive branch m = 1 to 4 over the Brillouin zone. The possibility that an order-by-disorder mechanism might be relevant for the ground state selection in CoTiO 3 was mentioned in [18], but no quantitative model was proposed. We show by direct calculations in Supplementary Note 8 that the semi-classical degeneracy is indeed lifted by zero-point fluctuations from bonddependent anisotropic couplings such as η ≡ J yy − J xx on the 1st neighbor bond where y defines the local bond direction and x is in-plane transverse to y, and we find an induced gap that scales as ∆ ∼ |η| 3/2 at leading order. At the level of the low energy pseudospin-1/2 moments this provides a natural qualitative mechanism for the observed gap. One can also place this finding in the context of a theory that operates within the full set of 12 single-ion spin and orbital states. In fact, working within the pseudospin-1/2 picture suggests an unphysically large coupling η calculated in Supplementary Note 8 compared to the coupling η fitted in Supplementary Note 7. Since the crystal field excitations are comparable to the exchange scale, an entirely natural mechanism for order-by-disorder to arise is through virtual crystal field fluctuations in a model that includes small spinorbital exchange. The virtual crystal field mechanism has been discussed in the context of Er 2 Ti 2 O 7 [26,27] − essentially the only other well-characterized example of order-by-disorder − where the linear spin wave mechanism and virtual crystal field mechanism are complementary. However, in CoTiO 3 virtual crystal field excitations are the leading cause of the discrete symmetry breaking.
Neutron Intensity Fingerprint of Magnon Isospin Winding -Having established the presence of bond-dependent exchange in this material, we now focus on the Dirac points in the magnon spectrum which provide an ideal setting to explore predicted intensity modulations associated with the isospin winding around nodal points. To explain this physics we use the simple example of a two-dimensional (2D) honeycomb Heisenberg ferromagnet H = −J i,j S i · S j (J > 0) taken from Ref. [30] to which we refer for further details, generalizations to band structures in 3D, and different types of touching points. The magnon band structure for this model computed within linear spin wave theory around the collinear ferromagnetic ground state has Dirac points at finite frequency at the corners (K-points) of the 2D hexagonal Brillouin zone (dashed outline in Fig. 2d). For a small momentum δk measured from a Dirac node, the effective spin wave Hamiltonian takes the famous form H eff = v δk · σ where v = 3JSa 0 /2 is the Dirac velocity (a 0 is the nearest-neighbor distance) and the isospin encoded in the Pauli matrices σ originates from the two sublattice honeycomb structure. By analogy with the Zeeman Hamiltonian, it follows that magnon wavefunctions carry an isospin polarization that is locked to the offset momentum δk thus winding around each Dirac point, see Fig. 2b). This feature is directly observable via INS because, in the vicinity of these points, the intensity is, up to a constant, the projection of the isospin polarization onto some directionn characteristic of each Dirac point [30], illustrated by the pink radial arrows in Fig. 2d). Explicitly, the intensity takes the form 1 ± cos(α − α 0 ) where α is the polar angle around the K point and α 0 defines the direction ofn, with the upper/lower sign for the top/bottom band, respectively. Therefore, the intensity winds smoothly around the Dirac point (as illustrated by the colour shading on the two conical bands in Fig. 2a).
Isospin of Dirac Magnons -CoTiO 3 provides a nearly ideal experimental platform to see the theoretically predicted winding of neutron intensity in the vicinity of the Dirac points. Fig. 2f) shows the INS data along the (1,1) in-plane direction through the nominal Dirac point at (2/3,2/3) where a clear near-linear band crossing is observed. In contrast, Fig. 2g) shows that the INS data through the same K point, but along the orthogonal (1,1) direction, has vanishingly small intensity in one of the two crossing bands. This strong intensity asymmetry in orthogonal scans is precisely what is expected based on the predicted isospin winding around a Dirac node in Fig. 2a). This can be seen more directly in Fig. 2c), which plots the intensity dependence as a function of angle α winding around the Dirac node in the top/bottom bands (filled/open symbols), the maxima and minima in each band are 180 • apart and in anti-phase between the two bands, the solid lines show fits to the generic form A ± ± B ± cos(α − α 0 ) with the upper/lower sign for top/bottom band. The fits give α 0 = −80(3) • , in good agreement with the XXZ∆ model for the same scan −81(1) • , the offset from −60 • is due to the buckling of the honeycomb layers, which rotate then vectors in plane upon varying L, for more details see Supplementary Note 5. The observed two-fold angular dependence is precisely the fingerprint of the predicted isospin winding for the near-nodal quasiparticles.
at the same energy and wind along L in a 'double-helix' [see Fig. 3b)], in opposite senses between adjacent Ktype points due to the3 point group symmetry of the crystal lattice. However, neither of those cases can explain the fine structure observed by the energy scan in Fig. 3d) centred at K-points, where two peaks are clearly resolved, 0.75(5) meV apart, instead of a single peak (XXZ∆ model, dashed red line, case ii) above; for case i) the single peak would be even sharper). This fine structure was not detected by earlier lower-resolution studies [18] and accounting for it requires anisotropic coupling terms beyond H XXZ . We have already argued that such terms must be present in order to account for the spectral gap. As shown in Supplementary Note 9, these terms all preserve the Dirac nodal lines while shifting their position in momentum space along in-plane directions related to the moment orientation in the ground state. To make quantitative contact with the experiment, we demonstrate that adding a finite nearest neighbor bonddependent exchange η leaves the dispersions largely unaffected relative to the η = 0 case in the magnetic Brillouin zone interior while leading to the observed double peak structure in Fig. 3d)(solid line).
Dirac Excitons -We now describe high-energy excitations, which we attribute to transitions to higher crystal field levels, where we also observe propagating excitations with linear band touching points and intensity winding around nodal points. The local spin-orbit coupled and trigonally distorted octahedral crystal field scheme for a Co 2+ (3d 7 ) ion (L = 3 and S = 3/2) is shown in Fig. 4a). Fig. 4b) shows INS measurements observing two peaks centred near 28 and 58 meV, which we identify with the (exciton) transitions to the two trigonally-split doublets of the j eff = 3/2 excited quadruplet (blue and red vertical thick arrow in Fig. 4a).
Fig. 4c) shows higher resolution INS measurements observing clear in-plane dispersions for the lower exciton modes near 28 meV, attributed to hopping due to spin and orbital exchange. Two modes are expected due to the two sublattices of the honeycomb structure and Fig. 4f) shows clear evidence for mode crossing at the two labelled nodal positions. Angular intensity maps around a nodal point in Fig. 4e) show a clear two-fold angular dependence, in anti-phase between the top/bottom bands (filled/open symbols), as expected from the intensity winding picture, again in complete analogy with the spectroscopic signature seen for the Dirac magnon wavefunctions in Fig. 2c). The observed dispersions and relative intensities of the two exciton modes can be well captured by a tight-binding model, detailed in Supplementary Note 4. The experimental and modelled exciton dispersions are compared in Figs. 4c) and d). We note that after this work was completed, Ref. [37] appeared, also reporting INS measurements of the exciton dispersion in CoTiO 3 .

Spin Orbit + Trigonal
Octahedral cubic [26,28] meV bottom band [28,30] meV top band d) c) Figure 4. Spin-orbit excitons: dispersions and Dirac node. a) Schematic level splitting for a Co 2+ ion in an octahedral crystal field of trigonal symmetry including spin-orbit coupling. b) INS energy scan observing transitions to the first two excited crystal field levels (the blue/red arrows above the peaks show the transitions indicated by matching colour vertical arrows in a), the solid line is a guide to the eye. c) INS data probing the dispersions of the first crystal level along high-symmetry directions, compared in d) with a tight-binding model (thick solid/dashed lines through both graphs show best fit dispersions). e) Angular intensity dependence around the nodal point (2/3,5/3) for the top/bottom exciton bands fitted to an A± ±B± cos(α−α0) form (solid lines, α0 = 155(3) • , calculated 153(1) • , in-plane radial wavevector range [0.075, 0.3]Å −1 ). The black squares/white circles denote the inelastic neutron scattering intensity for the top/bottom exciton bands, respectively, with error bars representing one standard deviation. Note the analogous behaviour to the intensity dependence in azimuthal scans for Dirac magnons in Fig. 2c). f) Exciton bands crossing at the two labelled nodal Dirac points, analogous to the magnon bands crossing in Fig. 2f). In e-f) intensities are averaged for L = [0, 3.5], in c) for a transverse wavevector range ±0.1Å −1 , and in f) for a transverse in-plane wavevector range ±0.025Å −1 . Data were collected at 8 K with Ei = 83 meV in b) and 45 meV in c,e,f). The colour bar in f) also applies to c) and d), and indicate scattering intensity in arbitrary units on a linear scale.

DISCUSSION
To summarise, we have reported INS measurements of the magnon dispersions in the stacked honeycomb CoTiO 3 , which reveal the presence of a spectral gap and Dirac nodal lines. We have shown that the gap implies the presence of significant bond-dependent anisotropic exchange originating from spin-orbit coupling and we have proposed a minimal model compatible with the experimental data to explain the discrete symmetry breaking via a quantum order-by-disorder mechanism. We have also observed key signatures of proximity to Dirac magnon physics through near-linear band touching and characteristic two-fold intensity periodicity in azimuthal scans attributed to the isospin winding around the Dirac node. The similar features seen also at the nodal band crossing in the spin-orbit excitons show that neutron scattering provides a window into the universal properties of highly constrained wavefunctions around linear band-touching points in bosonic systems in the solid state.

DATA AVAILABILITY
The experimental data in this study is available from Ref. [38]. support from a Royal Society University Research Fellowship. RC acknowledges support from the National Science Foundation under Grant No. NSF PHY-1748958 and hospitality from KITP where part of this work was completed. The neutron scattering measurements at the ISIS Facility were supported by a beamtime allocation [39] from the Science and Technology Facilities Council.

AUTHOR CONTRIBUTIONS
RC and PAM conceived research, DP synthesized the crystal and powder samples, ME prepared the multicrystal mount, ME, RC, PAM and HCW performed the INS experiments, ME, RC and PAM performed the analysis, PAM developed relevant theoretical models, PAM, ME and RC performed theoretical calculations, PM and RDJ performed powder neutron diffraction measurements and RDJ analysed this data, PAM, RC, ME and RDJ wrote the paper and the supplementary information with input from all co-authors.

COMPETING INTERESTS
The authors declare no competing interests.

Supplementary Information
Here we provide additional technical details on 1) the refinement of the crystal and magnetic structures from powder neutron diffraction, 2-3) calculation of the single-ion levels in the presence of spin-orbit coupling, trigonal crystal field, and exchange mean field to determine the spin and orbital contributions to the ordered moment in the ground state, 4) a tight binding model to capture the exciton dispersion, 5) the spinwave calculations for the minimal effective S = 1/2 XXZ model used to parametrize the magnon dispersions, 6) details of the INS experiments and quantitative fit of the observed dispersions, 7) symmetry-allowed bonddependent anisotropic exchanges, 8) quantum order-bydisorder from those terms as the origin of the spectral gap, 9) topology of the nodal lines of Dirac magnons and Hamiltonian symmetries, 10) flavor-wave theory based on a model with spin-orbital exchange that captures the discrete symmetry-breaking, the magnons and their spectral gap, as well as the dispersive excitons in a single model.

Supplementary Note 1. REFINEMENT OF CRYSTAL AND MAGNETIC STRUCTURES
Here we present neutron powder diffraction (NPD) measurements to determine the magnitude of the ordered  (2) moment in the ground state, which is an important ingredient in the parametrization of spin and orbital character of the cobalt magnetic moments. The experiments were performed using the WISH time-of-flight diffractometer [40] at ISIS, the UK Neutron and Muon Source. A high quality, single phase powder sample of CoTiO 3 (mass 3.125 g) was loaded into a 6 mm diameter vanadium can and mounted within an Oxford Instruments 4 He cryostat. High counting statistics data were collected at 1.5 and 150 K, representative of the magnetically ordered and paramagnetic phases, respectively (N.B. the paramagnetic data were collected well above T N = 38(3) K as magnetic diffuse scattering was found to persist above the transition). Additional lower counting statistics data were also collected on warming in 2.5 K steps between 1.5 and 50 K to obtain the order parameter. In the following analysis, Rietveld refinements of nuclear and magnetic structural models were performed using Fullprof [41], simultaneously against data measured in detector banks 2 and 9 (medium resolution, large d-spacing range) and banks 5 and 6 (high resolution, short d-spacing) of the WISH instrument. A small absorption correction was included in the refinements to account for moderate neutron absorption by cobalt.
The published ilmenite crystal structure of CoTiO 3 [42] (space group R3, herein defined using hexagonal axes in the obverse setting) was refined against the paramagnetic data (Supplementary Figures 1a-c). Excellent agreement between model and data was achieved and the crystal structure parameters are summarised in Supplementary Table I.
Below T N , more than 10 new diffraction peaks appeared (labelled "M" in Supplementary Figure 1d), which could be indexed using the T-point propagation vector Q = (0, 0, 3/2). Symmetry analysis performed using isodistort [43,44], showed that the full T-point magnetic representation for the cobalt Wyckoff positions decomposed into two 1D irreducible representations, T + , highlighting magnetic diffraction intensities labelled "M". The weak diffraction peak labelled by an asterisk (*) in (c) likely originates from a small CoTi2O5 impurity. Its presence did not affect the quantitative analysis of the diffraction pattern. Fits to the data are shown as solid black lines, and the difference I obs − I calc is given as a blue line at the bottom of the panes. In (a) and (c) nuclear peak positions are denoted by black tick marks, and in (b) and (d) nuclear and magnetic peak positions are denoted by top and bottom black tick marks, respectively. The temperature dependence of the cobalt magnetic moment evaluated by fitting variable temperature neutron powder diffraction data is shown in the inset to panel a).
tations, T + 2 ⊕T + 3 and T − 2 ⊕T − 3 . There exist four, symmetry distinct magnetic structures that transform by these four representations, respectively: (a) Ferromagnetic (FM) honeycomb layers stacked via antiferromagnetic (AFM) bonds with magnetic moments parallel to the c-axis (magnetic space group R I3 ), (b) AFM honeycomb layers (Néel-type) stacked via FM bonds with magnetic moments parallel to the c-axis (magnetic space group R I3 ), (c) FM honeycomb layers stacked via AFM bonds with magnetic moments perpendicular to the c-axis (magnetic space group P S1 ), and (d) AFM honeycomb layers (Néel-type) stacked via FM bonds with magnetic moments perpendicular to the c-axis (magnetic space group P S1 ).
We note that for structures (c) and (d) all in-plane moment directions are indistinguishable by symmetry. Furthermore, the T + 2 ⊕ T + 3 (T − 2 ⊕ T − 3 ) symmetry allows the T + 1 (T − 1 ) mode to appear via a secondary order parameter, which describes a global rotation of all moments out of the ab plane towards the hexagonal c axis whilst maintaining a collinear magnetic structure.
The largest magnetic diffraction intensity occurs for the magnetic Bragg peak indexed by the propagation vector Q = (0, 0, 3/2). Given that the magnetic neutron diffraction intensity is proportional to the component of the magnetic moments perpendicular to the scattering vector, this observation alone conclusively rules out structures (a) and (b) that have moments strictly parallel to the c axis. Furthermore, one can show that in the case of perfectly flat cobalt honeycomb planes (z Co = 1/3 in Supplementary Table I) the magnetic structure factor at Q is maximal for FM honeycomb planes and exactly zero for AFM honeycomb planes. The honeycomb planes of the true crystal structure are not perfectly flat, but the small buckling of these planes leads to only a few percent change in the predicted diffraction intensities. Hence, case (c) (illustrated in Supplementary Figure 5) is uniquely identified as the primary magnetic structure of CoTiO 3 by the observation of the largest intensity at the propagation vector alone, in agreement with earlier neutron powder diffraction results [42].
A magnetic structure model based on (c) was refined against the neutron powder diffraction data collected at 1.5 K (Supplementary Figures 1b and d). Excellent agreement between model and data was achieved (R p = 4.9%, R wp = 4.3%, R Mag = 3.1%). The in-plane direction of the magnetic moments cannot be determined from powder averaged diffraction data, and symmetry allowed out-of-plane tilting of the magnetic moments was found to be statistically insignificant. At 1.5 K the cobalt magnetic moment refined to 3.08(1) µ B . The temperature dependence of the magnetic moment was extracted from fits to data collected on warming and is shown in the inset to Supplementary Figure 1a).
The above magnetic structure has lower symmetry (P S1 ) than the paramagnetic crystal structure (R3). In this case the crystal symmetry can be lowered via magnetostriction. However we found that a hexagonal unit cell metric could be used to achieve excellent fits to our data at all measured temperatures and no peak splitting or significant peak broadening could be observed within the experimental resolution upon cooling below T N . We therefore estimate that any symmetry lowering of the hexagonal metric by the magnetic ordering involves changes in the lattice parameters below a conservative threshold of 0.02%.

Supplementary Note 2. SINGLE-ION PHYSICS
Here we discuss the ground state and higher-energy excited states of the Co 2+ (3d 7 ) ions given their local, octahedrally-coordinated crystal field environment and spin-orbit interaction, fitted to inter-level transitions observed in INS data. Hund's rules -appropriate to the case where the Coulomb interaction is greater than the crystal field -give a bare d 7 shell orbital triplet L = 3 and high spin S = 3/2. For an ideal octahedron, the crystal field acting on those levels has Hamiltonian H CF = −B O 0 4 + 5O 4 4 with B > 0, leading to a ground state triplet (Γ 4 ), and excited triplet and singlet levels, above energy gaps of 480 B and 1080 B, respectively. Those level splittings are of the order of 1 eV. Viewed another way, the crystal field levels are populated in the high spin t 5 2g e 2 g configuration, which is the aforementioned S = 3/2 orbital triplet. Cobalt(II) ions in octahedral environments may also occur in a low spin configuration with spin-1/2 degree of freedom and a two-fold orbital degeneracy, however CoTiO 3 is consistent with the high spin single-ion configuration because, as we show below, this offers a natural explanation for i) the observed transitions to higher single-ion levels and ii) the experimentally determined magnitude of the ordered moment in the ground state determined in Supplementary Note 1.
Empirically, the exchange scale and spin-orbit coupling in CoTiO 3 are both of order 10 meV. Since the octahedral crystal field splitting is larger than any other relevant magnetic scales, we may focus on the ground state orbital triplet as an effective l = 1 orbital angular momentum state with wavefunctions [45] in terms of the |L z states of the full L = 3 Hilbert space. The full angular momentum operator when projected onto the restricted l = 1 Hilbert space is expressed as L ≡ (−3/2)l. The spin-orbit coupling H SO = (3/2)λl · S with λ > 0 acts on the l = 1 and S = 3/2 states numbering 12 in all. It is convenient to define an effective angular momentum J eff = l + S as J eff is a good quantum number for the eigenstates of H SO . The spectrum is a spin-orbital ground state doublet with J eff = 1/2 at energy −15λ/4, a quartet (J eff = 3/2) at −3λ/2 and a 6-fold degenerate set (J eff = 5/2) at 9λ/4. The lowest doublet wavefunctions take the form in the |l z , S z basis.
In the actual crystal structure the oxygen octahedra around the Co ions are slightly trigonally distorted (the local point group symmetry at the Co sites is 3 instead of43m for a cubic octahedron) and this distortion can be parameterized by the term H trig = δ[l 2 z − (2/3)], where z denotes the c-axis. The level scheme in the presence of spin-orbit and trigonal distortion is summarized in Supplementary Figure 2(a). The trigonal distortion splits the levels into six Kramers doublets, with J z eff remaining a good quantum number. We now use the available experimental data to constrain the single-ion parameters. Inelastic neutron scattering has revealed the existence of single-ion levels at 28 and 58 meV (see Fig. 4a). A best fit to those levels gives δ = 45 meV and λ = 18 meV, consistent with earlier reports [18]. In the above analysis we have assumed δ > 0, as this gives larger magnetic moment in the ab-plane compared to along the c-axis, in agreement with single-crystal magnetic susceptibility measurements [31].
For the trigonally distorted case we examine the anisotropy of the magnetic moment in the ground state. This means that we compute matrix elements of the Zeeman coupled moment g l l + g S S within the ground state Kramers doublet described by an effective spin S = 1/2. Throughout we use Serif symbols S and S to refer to the effective spin-1/2 and SansSerif symbols S and S to refer to the real spin-3/2. Here g l = −3/2 and g S ≈ 2. Supplementary Figure 2(c) shows how the g-factors along and perpendicular to the trigonal axis vary with the reduced parameter δ/λ. The moments are isotropic for no trigonal distortion, but develop a strong easy-plane/axis anisotropy in the g-factor for +/−ve δ/λ. For the fitted 4 < l a t e x i t s h a 1 _ b a s e 6 4 = " f D 5 9 9 e g P J G R a + w 1 O h X U o Q J s c R i k = " > A A A B 7 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 B I v g q S R S 0 G P R g x 4 r 2 A 9 o Q 5 l s N + 3 S 3 U 3 c 3 Q g l 9 E 9 4 8 a C I V / + O N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Y c K Z N p 7 3 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p I r R J Y h 6 r T o i a c i Z p 0 z D D a S d R F E X I a T s c 3 8 z 8 9 h N V m s X y w U w S G g g c S h Y x g s Z K n d 4 t C o H 9 W r 9 c 8 a r e H O 4 q 8 X N S g R y N f v m r N 4 h J K q g 0 h K P W X d 9 L T J C h M o x w O i 3 1 U k 0 T J G M c 0 q 6 l E g X V Q T a / d + q e W W X g R r G y J Y 0 7 V 3 9 P Z C i 0 n o j Q d g o 0 I 7 3 s z c T / v G 5 q o q s g Y z J J D Z V k s S h K u W t i d / a 8 O 2 C K E s M n l i B R z N 7 q k h E q J M Z G V L I h + M s v r 5 L W R d X 3 q v 5 9 r V K / z u M o w g m c w j n 4 c A l 1 u I M G N I E A h 2 d 4 h T f n 0 X l x 3 p 2 P R W v B y W e O 4 Q + c z x + D g Y + c < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " f D 5 9 9 e g P J G R a + w 1 O h X U o Q J s c R i k = " > A A A B 7 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 B I v g q S R S 0 G P R g x 4 r 2 A 9 o Q 5 l s N + 3 S 3 U 3 c 3 Q g l 9 E 9 4 8 a C I V / + O N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Y c K Z N p 7 3 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p I r R J Y h 6 r T o i a c i Z p 0 z D D a S d R F E X I a T s c 3 8 z 8 9 h N V m s X y w U w S G g g c S h Y x g s Z K n d 4 t C o H 9 W r 9 c 8 a r e H O 4 q 8 X N S g R y N f v m r N 4 h J K q g 0 h K P W X d 9 L T J C h M o x w O i 3 1 U k 0 T J G M c 0 q 6 l E g X V Q T a / d + q e W W X g R r G y J Y 0 7 V 3 9 P Z C i 0 n o j Q d g o 0 I 7 3 s z c T / v G 5 q o q s g Y z J J D Z V k s S h K u W t i d / a 8 O 2 C K E s M n l i B R z N 7 q k h E q J M Z G V L I h + M s v r 5 L W R d X 3 q v 5 9 r V K / z u M o w g m c w j n 4 c A l 1 u I M G N I E A h 2 d 4 h T f n 0 X l x 3 p 2 P R W v B y W e O 4 Q + c z x + D g Y + c < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " f D 5 9 9 e g P J G R a + w 1 O h X U o Q J s c R i k = " > A A A B 7 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 B I v g q S R S 0 G P R g x 4 r 2 A 9 o Q 5 l s N + 3 S 3 U 3 c 3 Q g l 9 E 9 4 8 a C I V / + O N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Y c K Z N p 7 3 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p I r R J Y h 6 r T o i a c i Z p 0 z D D a S d R F E X I a T s c 3 8 z 8 9 h N V m s X y w U w S G g g c S h Y x g s Z K n d 4 t C o H 9 W r 9 c 8 a r e H O 4 q 8 X N S g R y N f v m r N 4 h J K q g 0 h K P W X d 9 L T J C h M o x w O i 3 1 U k 0 T J G M c 0 q 6 l E g X V Q T a / d + q e W W X g R r G y J Y 0 7 V 3 9 P Z C i 0 n o j Q d g o 0 I 7 3 s z c T / v G 5 q o q s g Y z J J D Z V k s S h K u W t i d / a 8 O 2 C K E s M n l i B R z N 7 q k h E q J M Z G V L I h + M s v r 5 L W R d X 3 q v 5 9 r V K / z u M o w g m c w j n 4 c A l 1 u I M G N I E A h 2 d 4 h T f n 0 X l x 3 p 2 P R W v B y W e O 4 Q + c z x + D g Y + c < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " f D 5 9 9 e g P J G R a + w 1 O h X U o Q J s c R i k = " > A A A B 7 3 i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 B I v g q S R S 0 G P R g x 4 r 2 A 9 o Q 5 l s N + 3 S 3 U 3 c 3 Q g l 9 E 9 4 8 a C I V / + O N / + N 2 z Y H b X 0 w 8 H h v h p l 5 Y c K Z N p 7 3 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p I r R J Y h 6 r T o i a c i Z p 0 z D D a S d R F E X I a T s c 3 8 z 8 9 h N V m s X y w U w S G g g c S h Y x g s Z K n d 4 t C o H 9 W r 9 c 8 a r e H O 4 q 8 X N S g R y N f v m r N 4 h J K q g 0 h K P W X d 9 L T J C h M o x w O i 3 1 U k 0 T J G M c 0 q 6 l E g X V Q T a / d + q e W W X g R r G y J Y 0 7 V 3 9 P Z C i 0 n o j Q d g o 0 I 7 3 s z c T / v G 5 q o q s g Y z J J D Z V k s S h K u W t i d / 4 v 5 B X c e p Y l h j s Y h V M 6 A a B Y + w Z r g R 2 E w U U h k I b A S D 6 4 n f G K L S P I 7 u z S h B X 9 J e x E P O q L G S f 9 v J 2 k o S D M P x w 1 O n W H L L 7 h R k m X h z U o I 5 q p 3 i V 7 s b s 1 R i Z J i g W r c 8 N z F + R p X h T O C 4 0 E 4 1 J p Q N a A 9 b l k Z U o v a z 6 d F j c m K V L g l j Z S s y Z K r + n s i o 1 H o k A 9 s p q e n r R W 8 i / u e 1 U h N e + h m P k t R g x G a L w l Q Q E 5 N J A q T L F T I j R p Z Q p r i 9 l b A + V Z Q Z m 1 P B h u A t v r x M 6 m d l z y 1 7 d + e l y t U 8 j j w c w T G c g g c X U I E b q E I N G D z C M 7 z C m z N 0 X p x 3 5 2 P W m n P m M 4 f w B 8 7 n D 6 S l k g E = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " Q 2 A 6 8 5 s E z m G 2 k 6 z y B F 1 b 9 n 6 v M T I = " > A A A B 9 H i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e g F / E U w T w g W c P s p D c Z M r O 7 z s w G 4 p L v 8 O J B E a 9 + j D f / x s n j o I k F D U V V N 9 1 d Q S K 4 N q 7 7 7 e R W V t f W N / K b h a 3 t n d 2 9 4 v 5 B X c e p Y l h j s Y h V M 6 A a B Y + w Z r g R 2 E w U U h k I b A S D 6 4 n f G K L S P I 7 u z S h B X 9 J e x E P O q L G S f 9 v J 2 k o S D M P x w 1 O n W H L L 7 h R k m X h z U o I 5 q p 3 i V 7 s b s 1 R i Z J i g W r c 8 N z F + R p X h T O C 4 0 E 4 1 J p Q N a A 9 b l k Z U o v a z 6 d F j c m K V L g l j Z S s y Z K r + n s i o 1 H o k A 9 s p q e n r R W 8 i / u e 1 U h N e + h m P k t R g x G a L w l Q Q E 5 N J A q T L F T I j R p Z Q p r i 9 l b A + V Z Q Z m 1 P B h u A t v r x M 6 m d l z y 1 7 d + e l y t U 8 j j w c w T G c g g c X U I E b q E I N G D z C M 7 z C m z N 0 X p x 3 5 2 P W m n P m M 4 f w B 8 7 n D 6 S l k g E = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " Q 2 A 6 8 5 s E z m G 2 k 6 z y B F 1 b 9 n 6 v M T I = " > A A A B 9 H i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e g F / E U w T w g W c P s p D c Z M r O 7 z s w G 4 p L v 8 O J B E a 9 + j D f / x s n j o I k F D U V V N 9 1 d Q S K 4 N q 7 7 7 e R W V t f W N / K b h a 3 t n d 2 9 4 v 5 B X c e p Y l h j s Y h V M 6 A a B Y + w Z r g R 2 E w U U h k I b A S D 6 4 n f G K L S P I 7 u z S h B X 9 J e x E P O q L G S f 9 v J 2 k o S D M P x w 1 O n W H L L 7 h R k m X h z U o I 5 q p 3 i V 7 s b s 1 R i Z J i g W r c 8 N z F + R p X h T O C 4 0 E 4 1 J p Q N a A 9 b l k Z U o v a z 6 d F j c m K V L g l j Z S s y Z K r + n s i o 1 H o k A 9 s p q e n r R W 8 i / u e 1 U h N e + h m P k t R g x G a L w l Q Q E 5 N J A q T L F T I j R p Z Q p r i 9 l b A + V Z Q Z m 1 P B h u A t v r x M 6 m d l z y 1 7 d + e l y t U 8 j j w c w T G c g g c X U I E b q E I N G D z C M 7 z C m z N 0 X p x 3 5 2 P W m n P m M 4 f w B 8 7 n D 6 S l k g E = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " Q 2 A 6 8 5 s E z m G 2 k 6 z y B F 1 b 9 n 6 v M T I = " > A A A B 9 H i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e g F / E U w T w g W c P s p D c Z M r O 7 z s w G 4 p L v 8 O J B E a 9 + j D f / x s n j o I k F D U V V N 9 1 d Q S K 4 N q 7 7 7 e R W V t f W N / K b h a 3 t n d 2 9 4 v 5 B X c e p Y l h j s Y h V M 6 A a B Y + w Z r g R 2 E w U U h k I b A S D 6 4 n f G K L S P I 7 u z S h B X 9 J e x E P O q L G S f 9 v J 2 k o S D M P x w 1 O n W H L L 7 h R k m X h z U o I 5 q p 3 i V 7 s b s 1 R i Z J i g W r c 8 N z F + R p X h T O C 4 0 E 4 1 J p Q N a A 9 b l k Z U o v a z 6 d F j c m K V L g l j Z S s y Z K r + n s i o 1 H o k A 9 s p q e n r R W 8 i / u e 1 U h N e + h m P k t R g x G a L w l Q Q E 5 N J A q T L F T I j R p Z Q p r i 9 l b A + V Z Q Z m 1 P B h u A t v r x M 6 m d l z y 1 7 d + e l y t U 8 j j w c w T G c g g c X U I E b q E I N G D z C M 7 z C m z N 0 X p x 3 5 2 P W m n P m M 4 f w B 8 7 n D 6 S l k g E = < / l a t e x i t > ± 1 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " Q R K s r K 9 P r l d f X a 7 N b n c N w b L r 4 h 8 = " < l a t e x i t s h a 1 _ b a s e 6 4 = " Q R K s r K 9 P r l d f X a 7 N b n c N w b L r 4 h 8 = "

l G 4 K 3 + P I y a d d r n l v z 7 s 6 r j e s i j h I c w w m c g Q e X 0 I B b a E I L C C h 4 h l d 4 c 5 6 c F + f d + Z i 3 r j j F z B H 8 g f P 5 A 5 d r k o 4 = < / l a t e x i t >
± 1 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " Q R K s r K 9 P r l d f X a 7 N b n c N w b L r 4 h 8 = "     crystal field scheme, g = 2.9 and g ⊥ = 4.95 so the ratio g ⊥ /g ≈ 1.7, and the expected in-plane moment would be g ⊥ µ B S = 2.5 µ B , which falls short of the 3 µ B found experimentally.  Here we show that the shortfall between the calculated magnetic moment in the single-ion picture and the experimentally-determined ordered moment can be explained naturally if exchange mean-field effects are incorporated into the single-ion picture. These can be parameterized by the Zeeman Hamiltonian H MF = −h ⊥ · S ⊥ where the ⊥ subscript indicates that the mean field is oriented in-plane, along the ordered moment direction in the magnetic structure.

Supplementary Note 3. MAGNITUDE OF THE ORDERED MOMENT
We treat h ⊥ as a variable parameter and solve for the magnetic moment in the ground state with the result shown in Supplementary Figure 3(a). We also compute the single-ion spectrum as a function of h ⊥ in Supplementary Figure 3(b) showing that the splitting of the first exciton around 28 meV is very small, consistent with the experimental finding.
We can estimate the magnitude of the actual mean field experienced by the magnetic moments from the XXZ exchange parameters that provide a quantitative description of the observed magnon dispersions in Supplementary Note 6. Those exchange parameters refer to an effective S = 1/2 spin model and, for this onsite moment, the magnitude of the mean field is h ⊥ = |S 6 n=1 (−) n+1 z n J ⊥ n | where z n is the number of n'th nearest neighbors and J ⊥ n is the n'th neighbor in-plane exchange. The sign (−) n+1 takes care of the antiferromagnetic arrangement between layers. For the exchange couplings given in Supplementary Equation (21), we find h ⊥ = 9.1 meV. We must then adjust for the bare moment by matching the Zeeman splittings for E MF = −h ⊥ · S ⊥ giving h ⊥ = 5.0 meV from which we read off from Supplementary Figure 3(a) a renormalized moment of about 3.0 µ B , as deduced experimentally. We also note that the spin-orbital mean field of Supplementary Note 10, with parameters chosen principally to match the spin wave bandwidth, gives a ground state ordered moment of about 2.8 µ B , again predicting an enhancement of the ordered moment compared to that of isolated ions and towards the value seen experimentally. In summary, we conclude that the enhanced magnetic moment seen experimentally is due to mean-field exchange effects.

Supplementary Note 4. TIGHT-BINDING MODEL FOR THE EXCITON DISPERSION
Here we outline the tight-binding model used to describe the observed dispersion [in Fig. 4c)] of the lowest crystal field excited level near 28 meV, attributed to hopping due to spin and orbital exchange. In a first approximation we neglect the effect of magnetic ordering on the crystal-field excitations and consider hopping of crystal-field excitations only between sites in the same honeycomb layer, so between sites of the A and B sublattices indicated in Supplementary Figure 5. Because the Kramers degeneracy of the crystal field modes is preserved and there are two sublattices in the paramagnetic regime, two dispersive bands are expected, analogous to the two bands of mobile electrons in graphene, which touch at the corners (K points) of the two-dimensional Brillouin zone [46]. A tight-binding description including 1st, 2nd and 3rd nearest neighbor in-plane (1st, 3rd and 5th in the full crystal structure) with hopping integrals t 1,3,5 on the same paths as J 1,3,5 in Supplementary Figure 5 gives the dispersions relations where and the cobalt positions r 1,2 and geometric factors γ nk are defined (later) in Supplementary Equation. (9) and (15). The above equations can capture well the observed dispersions of the exciton modes, see solid and dashed lines in Figure 4c). To find the model parameters experimental dispersion points were extracted from fitting Gaussian peaks to constant-energy and -wavevector scans through the high-energy INS data. From a best-fit to the experimental dispersion points we obtain The uncertainties in the fit parameters were obtained by adding Gaussian noise with a representative standard deviation σ = 0.3 meV to the energies of the experimentally extracted exciton dispersion points and fitting the model parameters for many such data sets. This resulted in a distribution of values for each of the model parameters, the quoted uncertainties are the standard deviations of those distributions. The hopping terms obtained above are of the same order of magnitude as the fitted exchange couplings presented in Supplementary Note 6 C. Fig. 4d) shows the intensity dependence assuming it is determined solely by interference scattering from the A and B sublattices, which takes the form [47] I ± (k) ∼ 1 ± cos ϕ k with the upper/lower sign for the top/bottom band and the phase angle ϕ k defined above. For wavevectors with in-plane projection in the vicinity of a K point, the phase ϕ k is directly related to the polar (azimuthal) angle α δk of the in-plane wavevector displacement δk away from K. In the limit |δk| → 0 and t 5 /t 1 → 0, the following relations are obtained for representative K-points where the azimuthal angle α δk is measured with reference to the ( 1 21 0) direction (horizontal axis in Fig. 2d). The relations near other K-points are obtained by3m symmetry. Here = 2(z Co − 1/3) characterises the buckling of the cobalt honeycomb layer ( = 0 for flat planes). The above intensity form captures well the overall intensity distribution and explains why only one exciton mode carries significant weight for most wavevectors in Fig. 4c) except the last panel where both modes are visible. Note that by simultaneously changing the sign of both t 1 and t 5 leaves the dispersion relations in Supplementary Equation (2) unchanged, however the intensities of the two modes then become completely the opposite way round to what is seen experimentally in Fig. 4c), so the parameter signs as listed in Supplementary Equation (4) are uniquely determined by combining dispersions and intensities constraints.
The above intensity form I(k) ± also explains the angular dependence of the intensity in the azimuthal scans in Fig. 4e) with maximum intensity in the top band occurring near α δ k = 155(5) • , compared to the predicted value of 153(1) • based on Supplementary Equation (5) [π − 2π L averaged for the appropriate L-integration range of the scan]. The tight-binding model in this Section provides a good empirical fit to the observed exciton data. In Supplementary Note 10 we treat the excitons and magnons in a unified way showing that the antiferromagnetic order should lead to further splitting of the exciton modes, however such splittings are expected to be small, beyond the resolution of the present measurements.

A. Structural and magnetic Brillouin zones
Here we describe the Brillouin zone relevant for the reciprocal space periodicity of the magnetic dispersion relations. We introduce the hexagonal primitive vectors a, b and c indicated in Supplementary Figure 5 and the primitive rhombohedral unit cell with vectors such that c = A 1 + A 2 + A 3 . The Brillouin zone corresponding to this primitive structural cell is illustrated in Supplementary Figure 4a) and belongs to the elongated (c > 3/2a) rhombohedral case [48]. It has top and bottom regular hexagonal faces with midpoints at ± 00 3 2 and twelve side faces that alternate between rectangular and slightly-distorted hexagonal with midpoints at 1 2 01 and 1 2 0 1 2 , respectively, with other faces obtained by3 symmetry.
The magnetic structure is illustrated in Supplementary  Figure 5 and has moments parallel in each layer and antiparallel between layers. This magnetic periodicity can be captured by a doubled-volume rhombohedral primitive cell shown by the dashed outline in Supplementary Figure 5, with basis vectors rotated by 60 • and elongated a factor of 2 along c compared to the rhombohedral primitive structural cell in Supplementary Equation (8), with the magnetic primitive unit cell vectors given by The Brillouin zone corresponding to this primitive magnetic unit cell is shown in Supplementary Figure 4c are hexagonal with midpoints at ± 00 3 4 and the twelve side faces alternate between rectangular and strongly distorted hexagonal with midpoints at 1 2 0 1 2 and 1 2 01 4 , respectively, with other faces obtained by3 symmetry.
The projection of the magnetic Brillouin zone in the hk plane is illustrated in Supplementary Figure 4d), where the inner hexagon corresponds to the top face at L = 3/4. In projection, this is located inside the 2D hexagonal Brillouin zone (red dashed outline) of a single honeycomb layer. Upon decreasing L, the corners of the magnetic Brillouin zone move initially outwards, along the set of small black segments, followed by two other small segments, such that in projection they describe small equilateral triangles centred at the nominal K-points of the 2D hexagonal Brillouin zone. Viewed in 3D, the magnetic Brillouin zone edges wrap around the straight lines that project onto K-points, as illustrated in Supplementary Figure 4b). This has the consequence that points along those straight lines have no special symmetry, they act like general points in the magnetic Brillouin zone. Therefore, there are no symmetry-imposed constraints for touching points in the magnetic dispersion bands to be pinned at those positions. Indeed, as illustrated in Fig. 3b) and detailed in Supplementary Note 7 A and Supplementary Note 9, we find that in the general case nodal lines wind along L and precess in-plane around positions that can be displaced away from the nominal K-points.

B. The XXZ Model and Further Neighbor Couplings
Here we give details of the analytical spin-wave calculations for the minimal easy-plane exchange model that captures the principal features of the observed magnon dispersion relations. To describe the full structural arrangement of the Co ions we use the following primitive unit cell with vectors and vectors defining the positions of the two cobalt ions in this primitive cell where = 2(z Co − 1/3) characterises the buckling of the cobalt honeycombs, with the Co z-coordinate z Co given in Supplementary Table I. The above primitive cell was chosen to emphasize the ab planes as a 'natural' building block of the Co structural arrangement. The full structure of cobalt ions is then generated by the set of positional vectors where the integers n ia select the primitive unit cell and m = 1, 2 is the cobalt sublattice index. The minimal model that we find describes all but the fine structure of the spin wave spectrum is an XXZ model on all these couplings: where J z n is the (Ising) coupling for spin components along z c and J ⊥ n is the coupling for spin components in the ab plane, the summation is over all interacting pairs i, j n of n'th nearest neighbors counted once, and we include n=1 to 6.

C. Spin Wave Calculations
The magnetic structure shown in Supplementary  primitive magnetic cell. For the spin-wave calculation it is convenient to define a global Cartesian xyz frame with x a, z c and y = z × x, as illustrated in Supplementary Figure 18d). We also define a local frame denoted xỹz wherez lies along the direction of the local ordered moment andx is in-plane, where the two frames are related by where φ is the in-plane angle of the ordered spins of sublattices A and B, measured relative to x in the positive sense if rotating around z and n = 0 for sublattices A and B, and n = 1 for sublattices C and D.
In the local frame all moments are parallel and the magnetic primitive cell is reduced to the structural primitive cell, i.e. in this frame there are only two magnetic sublattices as opposed to four in the original frame. For the XXZ model in Supplementary Equation (11) the spin Hamiltonian expressed in the local frame has the same periodicity as that of the structural cell and in this case the problem is reduced to obtaining the dispersion relations for a Hamiltonian of the generic form , with a † k and b † k being magnon creation operators for the two magnetic sublattices. Here the sum extends over all wavevectors k in the structural Brillouin zone in Supplementary Figure 4a).
Introducing A, B, C and D as implicit functions of k, the dynamical matrix D(k) has the form Including couplings up to 6th nearest neighbor we find Here with coefficients given in Supplementary Table II for In order to compute the neutron scattering intensities, we require the right eigenvectors of GD. The components are The eigenvectors are then the columns of  where the bar meansW (ω) ≡ W (ω)/ N (ω) and so on.
Let us define where implicitly the functions W , X, Y , Z, N are evaluated at wavevector k, and energy ω that satisfies the delta functions δ(ω − ω ± (k)). The in-and out-of-plane dynamical correlations in the local frame are obtained as Sxx 0 (k, ω) and Sỹỹ(k, ω), respectively, so both magnon modes ω ± (k) occur in both polarizations with different intensities. Note that the magnon dispersion relations are independent of the buckling parameter , which only has the effect of modulating the intensities through the exponential phase factor in Supplementary Equation (18). Upon transformation to the global frame the out-ofplane correlations are unchanged. The in-plane polarized dispersions are momentum shifted by the propagation vector Q such that the total INS intensity is proportional to where we have included the neutron polarization factors p ζ = 1 − (k ζ /k) 2 and the anisotropic g-factor components for in-and out-of-plane moment directions, g ⊥ and g , respectively. Here k ζ = k · ζ is the wavevector transfer component along the ζ =x,ỹ direction. In the global frame at a given wavevector k there are four magnon branches, ω ± (k) (polarised out-of-plane along z =ỹ) and ω ± (k + Q) (polarized in-plane), i.e. one recovers the expected number of modes for the four magnetic sublattices in the original problem. In the above we have used the fact that 2Q is a vector of the reciprocal lattice of the structural cell, so ω(k + Q) and ω(k − Q) are identical by reciprocal space translational symmetry. The above analytical expressions for the dispersion relations and intensities have been checked against a numerical spin-wave code for the full four-sublattice spin-wave Hamiltonian, and also against SpinW [49]. The spectrum is gapless at the origin and at the magnetic propagation vector Q, as expected given the U (1) symmetry of the XXZ Hamiltonian, i.e. there is no energy cost in rotating the spins in the ab plane. The dispersion relations as well as the functional form of the dynamical structure factors are independent of the inplane angle φ, the only dependence of the INS intensity on φ comes through the neutron polarization factor for inplane fluctuations, px. This is the case for a single magnetic domain, however assuming six equally-populated magnetic domains with moment directions at φ + nπ/3 (n = 0 to 5) as expected due to the3 lattice point group symmetry of the crystal structure, the domain average px is independent of φ and the total INS intensity in this case is proportional to where k z = k · z is the wavevector transfer component along the c-axis.
For the purpose of comparison with data it is helpful to discuss the overall reciprocal space periodicity and symmetry of the dispersions and dynamical correlations. The Brillouin zone folding that occurs upon going from the local to the global frame has the consequence that each of the dispersive modes in the global frame when labelled as ω m (k) with m = 1 to 4 in order of increasing energy, has the translational periodicity of the magnetic Brillouin zone. The magnetic structure in Supplementary Figure 5 breaks the 3-fold rotational symmetry of the crystal structure, however the XXZ Hamiltonian has a higher symmetry, U (1), than required by the crystal structure, with the consequence that the magnon dispersions ω m (k) are independent of the in-plane moment's angle φ and have the rotational symmetry of the cobalt structural arrangement, which is3m. This rotational symmetry implies also that all magnetic domains will have identical dispersion relations, which justifies Supplementary Equation (20) for a multi magnetic domain sample. Each of the two intensity terms in that equation, separately has the same rotational point group symmetrȳ 3m. In the case of flat plane honeycombs ( = 0), each of those two intensity terms also has the translational periodicity of the structural Brillouin zone, however the buckling of the layers breaks the translational periodicity of the intensity along the L-direction as it introduces an intensity modulation factor due to interference scattering from the two cobalt sites in the same honeycomb layer being offset along z, this intensity modulation term has a long period, 1/ along L.
Finally, we note that following the general arguments presented in [30], the dynamical structure factor for small in-plane wavevector displacements δk away from the nodal points will have the azimuthal angular dependence 1 ± cos ϕ k , where the phase angle ϕ k is related to the azimuthal angle α δk of δk via Supplementary Equations (5)- (7). This leads to a two-fold intensity modulation in azimuthal scans, in anti-phase between the two touching bands, as observed by the data in Fig. 2c). The relation between ϕ k and α δk varies between neighboring K-points following a3m symmetry. This is illustrated in Fig. 2d) by the radial thick magenta arrows which show the directionsn away from the nearby K-points along which the intensity in the the top band is maximal in azimuthal scans at L = 0. At finite L, due to scattering interference from the two cobalt sites offset along z, then vectors rotate in-plane by an angle 2π L, in opposite senses for adjacent K-points, following a3m symmetry. This L-dependence provides a natural explanation for the observed angular intensity dependence in the azimuthal scan in Fig. 2c) around the (2/3,2/3) Dirac node, with maximum intensity in the top band observed near α δ k = −80(3) • , compared to −81(1) • calculated based on Supplementary Equation (6) averaged for the appropriate L-integration range of the scan. The L-dependence of the azimuthal scans is illustrated in Supplementary  Figure 6.

A. Experimental Details
Here we provide details of the INS experiments [39] to probe the spin dynamics, performed using the MERLIN direct-geometry time-of-flight spectrometer [50] at ISIS. The sample consisted of two co-aligned single crystals of CoTiO 3 (total mass 5.8 g) grown via the floating zone method, mounted with the (hk0) scattering plane horizontal. Full Horace maps of the inelastic scattering were collected at a base temperature of 8 K (cooling was provided by a closed-cycle refrigerator) by rotating the sample around the vertical axis in steps of 0.5 • over an angular range of 120 • , with each step counted for 9 mins at an average proton current of 170 µA. The temperaturedependence of the inelastic scattering up to 300 K was measured for one representative sample orientation. The spectrometer was operated in repetition rate multiplication (RRM) mode to collect the inelastic scattering simultaneously for monochromatic incident neutrons with energies E i =9.6, 18 and 45 meV, with energy resolutions on the elastic line of 0.36(2), 0.72(2) and 2.7(1) (full width at half maximum, FWHM), respectively. Additional measurements were collected with E i =83 meV to probe transitions to higher crystal field levels. The elastic line in all runs was centred on zero energy transfer to better than 0.25% of E i . The time-of-flight neutron data were processed using the mantid [51] and horace [52]  where interference scattering from Co, Ti and O leads to a strong intensity for domain A, but near cancellation for domain B, and viceversa for reflection (231). The observed diffraction pattern showed almost equal intensities for those two reference reflections, so we conclude that the sample contained equal amounts of the A and B domains, which would imply a3m point group symmetry for the (diffraction and inelastic) signal. Indeed the inelastic intensity showed to a very good degree3m symmetry with mirrors at (h0l) and to enhance the counting statistics the wavevector transfers k of the pixels in the four-dimensional Horace scans were remapped using symmetry operations of the above point group to a minimal 60 • sector in the hk plane and l > 0.
We note that magnetic ordering with moments in plane breaks the 3-fold rotation, so the dispersion relations and dynamical structure factor for a single magnetic domain would in principle have point group1 (a minimal model that exhibits this lower symmetry of its magnetic spectrum is the XXZ Hamiltonian augmented by finite diagonal exchange anisotropy η = 0 discussed in Supplementary Note 7 A. However, averaging over three equal-weight magnetic domains with moments rotated by ±120 • as expected in a macroscopic sample, would restore the 3-fold symmetry for the intensity pattern. This combined with the A and B structural domains then would restore the higher point group symmetry3m for the intensity pattern, justifying the pixel averaging used.  Fig. 1c). The colour bar indicates the scattering intensity in arbitrary units on a linear scale.

C. Parameterization of magnon dispersions by an XXZ∆ model
The XXZ Hamiltonian discussed in the previous section has a U (1) symmetry, however the crystal structure has only discrete rotational symmetry and moreover the observed magnon spectrum is clearly gapped, as shown in Fig. 1c) and Supplementary Figure 8, indicating not a continuous, but a discrete set of allowed φ values. To account for the presence of a spectral gap at this stage we introduce a phenomenological gap parameter ∆ and assume that the effect of the symmetrybreaking interactions that generate this gap can be accounted for, in a first approximation, by simply adding this gap in quadrature to the analytical XXZ dispersions, i.e. the experimental dispersion points are compared withω m (k) = ∆ 2 + ω 2 m (k), where m = 1 to 4 labels the four magnon modes at a given wavevector k in order of increasing energy. Empirical (h, k, l, E, m) dispersion points (where E is energy) were extracted from fitting Gaussian peaks to constant-wavevector and/orenergy scans through the four-dimensional inelastic neutron scattering data. Supplementary Figure 9 illustrates the level of agreement that can be obtained when comparing nearly 2, 000 empirical dispersion points with dispersions of the XXZ∆ model for a representative set of exchange parameters below, all in meV, and∆ = 1.23(7) meV, where −/+ve signs for the exchanges mean FM/AFM coupling. In the above we have used the symbol∆ to indicate that the gap is overestimated through this analysis. There is a net shift of the scattering weight towards higher energies originating from the finite wavevector integration around the lowest energy mode because the integration range captures intensity from the mode away from the minimum, thus shifting the average upwards. We account for this effect in a first approximation by assuming all exchange parameters fixed as per Supplementary Equation (21) and calculating the expected scattering for the slice in Fig. 1c), which is most sensitive to the gap, allowing for a variable ∆ in the fit. We include the full wavevector averaging in the transverse (highly-dispersive, in-plane) direction as in the data slice, and optimise ∆ to get the best agreement between scans through the data and simulation, as shown in Supplementary Figure 10. This gave ∆ = 1.0(1) meV, renormalized down from∆. Fig. 1c)  dispersions (up to the fine structure around the K points to be discussed later) and qualitatively of the intensities as well, compare Fig. 1a) and b). However, the number of Hamiltonian parameters considered is large and we have found that some parameters are strongly correlated in their effects on the dispersions. Therefore, the dispersions alone are not constraining enough to uniquely determine the values of all the individual exchange parameters. For example, moving in parameter space away from the set in Supplementary Equation (21) by varying the value of J z 1 , fixing J ⊥ 1 = −J z 1 − 4.4 meV and optimising all the remaining exchange parameters results in a very small relative change in χ 2 = i |ω obs (i) − ω calc (i)| 2 , of only a few percent when J z 1 is reduced all the way to 0, such a small variation in χ 2 is at the level that the change in the agreement with the data is almost indistinguishable.
Here ω obs (i)/ω calc (i) is the observed/calculated energy for the ith dispersion point. With the exception of J ⊥ 1 , any one of the other 11 exchange parameters can be set equal to 0 and optimising the rest of the parameters gives a comparable agreement to that in Supplementary Figure 9. Therefore, more constraints are needed to uniquely identify the values of the individual exchange parameters and we therefore regard the set in Supplementary Equation (21) as representative of the best agreement that can be obtained with the measured dispersions, and in the following we focus on the key features of the measured dispersions.
Whilst the overall dispersion trends are in general in agreement with the minimal model parametrization proposed in [18], our higher-resolution INS data reveal additional dispersion modulations and fine structure (splitting of modes) that require additional couplings and anisotropies. For example Supplementary Figure 11 shows a clear splitting between the two lower modes (gray and magenta solid dots) in the region in-between the two labelled K-points, those two lower modes would be almost degenerate in this region in the parametrization used in [18]. The model in Supplementary Equation (21) (lines) predicts a substantial splitting between those modes, although it still underestimates the magnitude of the splitting seen experimentally. The agreement with the data can be improved by adding a bonddependent anisotropic exchange η, which, we argue, is physically responsible also for generating the finite spectral gap above the magnetic Bragg peaks, to be discussed in the following two sections.  (21) with anisotropic XXZ inter-layer couplings. Heisenberg inter-layer couplings as in [18] would give almost degenerate lower modes in this region (not shown). Solid dots are empirical dispersion points (colour indicates the mode index in order of increasing energy as per the legend in Supplementary Figure 9. The upper (yellow) and lower (magenta) pair of modes do not touch along this wavevector path. The nodal lines are present for these parameters, but are displaced away from the two K points along directions that make a finite angle with the plotted wavevector path. The colour bar indicates scattering intensity in arbitrary units on a linear scale.
For completeness, we note that the dipolar couplings are negligible compared to the scale of the above exchanges, i.e. the dipolar energy scale is µ 0 µ 2 /4πa 3 0 = 0.018 meV, where µ is the ordered magnetic moment per site in the ground state (3 µ B ) and a 0 ≈ 3Å is the nearest-neighbor Co-Co distance.
For the above J 1 − J 6 XXZ model the reduction of the ordered moment due to zero point fluctuations ∆S within linear spin wave theory is merely 0.021 compared to about 0.1 for the nearest-neighbor XXZ model (with only J ⊥ 1 and J z 1 nonzero).

Supplementary Note 7. SYMMETRY AND ANISOTROPIC EXCHANGE
In Supplementary Note 6 we showed that the spin waves computed from an XXZ model capture most of the features of the experimental neutron scattering data very well. However, it was necessary to include a phenomenological gap, which is absent in the XXZ model. We now begin to address the microscopic origin of the gap and we do this in two parts. The first is to recognize that spin-orbit coupling can lead to anisotropic exchanges beyond the XXZ model which arises from a projection of a Heisenberg model down to the spin-orbital ground state Kramers doublet. In this section we discuss the possible anisotropic exchange from a phenomenological point of view. These additional anisotropies in the effective spin one-half description generally lead to quantum order-bydisorder as outlined in Supplementary Note 8 in which we also compute the order-by-disorder gap exactly to leading order in 1/S for particular anisotropic exchange couplings in order to estimate the required magnitude of the anisotropies. In the second part, Supplementary Note 10, we consider the microscopic origin of the anisotropic exchange and directly compute the spectral gap through a flavor wave mean field theory.
We start with the nearest neighbor bonds. The inversion symmetry at the midpoint of the bond forces this exchange to be symmetric. In other words, six independent coupling terms are allowed, which, of course, includes the XXZ form. Once the anisotropic exchange is defined for one of the bonds chosen as "reference", the exchange on all other bonds of the same type in the full crystal lattice is obtained using crystal lattice symmetry operations such as the 3-fold rotation at the cobalt sites or primitive lattice translations. It is insightful to consider two (bond-dependent) reference frames to define the anisotropic exchange. The global xyz frame used so far is a natural frame for the vertical diamond-shaded bond in Supplementary Figure 18d) as the y axis is along the bond direction and z is along c. In this frame, the exchange matrix is It is also convenient to write the exchange in a Cartesian 111 frame (also illustrated in Supplementary Figure 18d) and denoted by SansSerif symbols xyz), where the axes have the property that the hexagonal c-axis is along the symmetric combination (x +ŷ +ẑ)/ √ 3. This frame transforms simply under the lattice point group symmetry 3 and is natural from the point of view of the underlying superexchange mechanism. In this frame the exchange matrix for the same bond has the form where three of the six independent terms have a natural interpretation, i.e. J is the isotropic (Heisenberg) exchange, K is a bond-Ising, or Kitaev-like, coupling, and Γ is an off-diagonal symmetric exchange for components in the plane orthogonal to the Kitaev axis. It is possible to justify microscopically the origin of those three couplings using a minimal superexchange model for the pair of 90 • Co-O-Co bonds (Supplementary Note 10). The transformation that converts between the exchange terms between the two coordinate frames is In the following we shall consider anisotropic exchange beyond nearest neighbor in particular to address the effect of these couplings on the fine structure of the spin wave spectrum at the Dirac nodes. For 2nd nearest neighbors, symmetry is highly constraining and restricts the exchange to the XXZ form expressed in the xyz frame as The 3rd nearest neighbor bonds are all constrained by lattice symmetries once a single such bond is fixed. However, a single bond has no symmetry on its own and nine independent exchange couplings are allowed including three diagonal couplings, three symmetric off-diagonal couplings and three Dzyaloshinskii-Moriya couplings. The same is true for 4th nearest neighbor couplings, while 5th and 6th nearest neighbor exchange couplings are constrained not to have antisymmetric exchange (as is the case for nearest neighbors) so they have six couplings each.

A. Bond-dependent Anisotropic Exchange
In the next section, we shall find it useful to consider the anisotropic exchange component on the nearest neighbor bond η ≡ J yy − J xx . One can compute the spin wave spectrum in the presence of this coupling within the local frame of Supplementary Note 5 C with the B and D parameters in Supplementary Equation (14) acquiring the additive terms B and D , respectively, with where β i = k·R 1i with the vectors R 1i , i = 1−3 given in Supplementary Table II. Note that for finite η the spectrum depends on the spin orientation angle φ in the easy plane and we will show in the subsequent section that this feature is responsible for selecting discrete orientations via quantum order-by-disorder, φ = nπ/3 for η > 0, and φ = π/2 + nπ/3 for η < 0, n integer.
For small η, significant changes to the spectrum occur only near the magnetic Brillouin zone boundary where the dispersion surfaces are distorted/shifted along a direction that correlates with the moments' orientation in the ground state, leading to an overall1 point group symmetry, the same as the magnetic structure. Magnetic domains obtained by ±120 • rotation therefore can have distinct dispersion relations and this can provide a physical mechanism to account for several features in the data, in particular the splitting between the two lower modes in Supplementary Figure 12, where solid/dashed lines correspond to magnetic domains rotated by 120 • . This path was chosen because the observed splitting cannot be captured within an XXZ model, where the combined rotational (3m) and translational symmetry of the spectrum requires the two lower modes to be degenerate. This can be seen as follows: when mapped to the magnetic Brillouin zone in Supplementary Figure 4b), this path is equivalent to one of the main diagonals (blue dotted line) of the hexagonal top face (the parallel blue dotted segment on one of the side faces is simply an extension of the top segment into the next zone, mapped back into the first magnetic Brillouin zone). The top diagonal is mapped by rotational3m symmetry operations, i.e. 2 (1 2 10) followed by 3 (001), onto the same diagonal on the bottom face, translated by −Q. This means that for any wavevector k on the top diagonal ω − (k) = ω − (k − Q), implying that in the global frame ω 1 (k) = ω 2 (k), i.e. the lower two modes are degenerate along this path. Furthermore, the same rotational symmetry implies that all magnetic domains have the same dispersion relations, so the observed splitting cannot be explained within an XXZ model. η = 0 breaks this symmetry requirement allowing magnetic domains with spins rotated by 120 • to have non-overlapping dispersions, thus providing a natural mechanism to explain the observed splitting.
The η interaction can also account for the fine structure in the energy scan centred at K-points in Fig. 3d), which cannot be explained by the XXZ model (dashed red line). To illustrate this, it is helpful to consider first the effects of adding a finite η to a single, isolated honeycomb with only the dominant nearest-neighbor exchange J ⊥ 1 < 0 with all other exchanges set to zero. A contour plot of the upper magnon band in this case is shown in Supplementary Figure 13a), the Dirac cones are centred at the K-point corners of the 2D hexagonal Brillouin zone with the magnon band displaying a3m point group symmetry around the zone centre. Supplementary Figure 13b) shows the corresponding plot for small η > 0 and the moment's orientation in the ground state along the φ = 0 direction, indicated by the dashed white arrow. The magnon dispersion surface has now1 symmetry, being distorted along the ordered in-plane spin direction in the ground state, with the Dirac nodes moving along this direction (two nodes move in and four nodes move out of the Brillouin zone hexagon). The corresponding dispersion plots for the above two cases are shown in Supplementary Figure 13c), dashed red lines for case a), and solid blue lines for case b), with a zoomed-in version near the K point shown in Supplementary Figure 13d), note the Dirac nodal points have moved away from K-points, so a scan at K would show a double-peak structure, as in Fig. 3d).
For the full XXZ model the nodal points occur in the form of double helix nodal lines that precess around Kpoints, adding a finite η has the effect of shifting in-plane the centre of precession away from K-points along a direction parallel to the moments' direction in the ground state if η > 0 and φ = nπ/3, and transverse to the moments' direction if η < 0 and φ = π/2 + nπ/3. In both cases, the precession centres of the double-helix nodal lines in Figure 3b) are displaced away from the K-points by an in-plane wavevector of magnitude δκ, with the consequence that at K-points there is an energy gap 2 v δκ between the mean of the top two bands and the mean of the bottom two bands, where v is the Dirac velocity, similar to the simplified 2D case illustrated in Supplementary  Figure 13d). The scan in Fig. 3d) includes wavevectors in a narrow cylindrical region centred at K and extending out to just touch the double helix nodal lines, so for all of those wavevectors there will be a finite energy gap between the top and bottom sets of bands, so a two-peak structure would be expected in the energy scan, as indeed seen experimentally. For a quantitative comparison with the data, we use the cross-section model in Supplementary Equation (19) averaged over all magnetic domains. In the present case, it is sufficient to consider only three magnetic domains (of the A structural domain), so for η > 0 we take φ = 0, ±2π/3, as for those φ values each B magnetic domain has the same response as the A magnetic domain with spins along the same direction, and domains obtained by reversing the spins on each site also have identical signatures. We consider two scenarios related by the transformation η → −η and φ → φ + π/2, as this leaves B and D in Supplementary Equation (26) invariant, and subsequently the dispersions and dynamical structure factor in Supplementary Equations. (17,18) are invariant as well. The above two scenarios indeed produce a similar, but not identical intensity profile in Supplementary Figure 14 (dashed gray and solid black lines), the small differences are due to the polarization factor px in Supplementary Equation (19), which changes upon rotating the spins by 90 • between the two scenarios. The first peak in the energy scan is identified with crossing the two lower bands (almost overlapping ω 1,2 modes) and the higher peak with crossing the top two bands (almost overlapping ω 3,4 modes), with the observed energy separation well accounted for in Supplementary Figure 14 using |η| = 1.7 meV. The calculated lineshape for the case η < 0 (black solid line, lower peak more intense) is in better agreement with the data than for η > 0 (gray dashed line), from which we conclude that the former is the more likely of the two scenarios. For the parameters with best agreement we show in Supplementary Figure 15c-d) the calculated momentum intensity maps, the agreement with the data in panels a-b) is excellent. In particular, the observed dramatic change when moving from energies below (bottom panels) to energies above the nodal energy (top panels) -the intensity shift from inside to outside of the central hexagonal Brillouin zone (dashed hexagonal outline) -is well captured, even the subtle local rotations of the intensity patterns around the zone corners -most visible in the lower-right corners -are well reproduced (those intensity modulations arise from the finite buckling of the cobalt honeycombs).
The coupling η is one of four anisotropic nearest neighbor couplings beyond the XXZ model. Establishing in detail the magnitudes and signs of all such anisotropic exchanges is beyond the scope of the present work. However, we have demonstrated in this Section that a finite η provides a natural explanation for the double peak structure in Supplementary Figure 14  Left/right panels are data/calculation, top/bottom panels are above/below the Dirac node energy. The data panels are as in Fig. 2d-e) and the calculation parameters are as in Supplementary Figure 14  pling can also provide a mechanism to account for the appearance of the spectral gap.

A. Mean-field ground-state degeneracy
Upon including exchange anisotropy terms as discussed in the previous Section, the spin Hamiltonian symmetry is reduced down to discrete rotations, however at the mean-field level the ground state energy remains independent of the in-plane moment angle φ and consequently the linear spin-wave spectrum remains gapless. To see this, we parameterize the easy-plane spin configuration by m i = m 0 Re e iφ x i + iẑ i where thex andẑ axes define the local frame, which rotates 180 • around thex-axis between adjacent layers. Here, m 0 is the size of the ordered moment. We then note that the mean field free energy for quadratic spin interactions can be written as F MFT [Ψ] = α|Ψ| 2 + βΨ 2 + β (Ψ * ) 2 where Ψ = m 0 e iφ and invariance under 3-fold symmetry forces β = 0 as Ψ → e 2πi/3 Ψ under this symmetry operation. This means that the mean field free energy cannot depend on φ, even when symmetry-allowed two-spin exchange anisotropy terms are included. A similar argument forbids symmetry breaking from four-body couplings that may arise from spin-lattice coupling. Indeed, the lowest order terms that lead to discrete symmetry breaking are six-body terms that may be attributed to fluctuations. In short, we expect order-by-disorder to arise in CoTiO 3 purely on the basis of the spectral gap.

B. Quantum order-by-disorder in the presence of bond-dependent anisotropic exchange η
A more direct way to see this is that anisotropic exchange couplings break the U (1) Hamiltonian symmetry down to the lattice symmetries so one would expect the classical ground-state degeneracy as a function of φ to be lifted by zero-point quantum fluctuations. To leading order the zero point contribution to the ground state energy (per magnetic unit cell) is where N c is the number of magnetic unit cells, each containing N s = 4 magnetic sublattices. The sum runs over all wavevectors k in the magnetic Brillouin zone and all magnon branches indexed by m = 1 to N s . To illustrate the phenomenon we add to the XXZ model in Supplementary Note 6 a symmetry-allowed nearest-neighbor anisotropic diagonal coupling η defined in the xyz frame by J xx = J ⊥ 1 −η/2 and J yy = J ⊥ 1 +η/2, i.e. +/−ve η favors spins pointing orthogonal/parallel to the bond direction (as J ⊥ 1 < 0). At finite η the (numerically calculated) spin-wave dispersions do depend on the in-plane moment angle φ and the leading zero-point correction to the ground state energy is six-fold modulated in φ with period π/3 as shown in Supplementary Figure 16a) and can be fitted (solid line) to a cosinusoidal form where Λ is an odd polynomial in the anisotropic exchange η [see Supplementary Figure 16b)]. In other words, there is a quantum order-by-disorder mechanism that lifts the U (1) classical ground state degeneracy leading to a sixfold symmetric set of ground states, i.e. +/−ve η select the family of orientations φ = 0 or π/6 (modulo π/3), respectively.

C. Order-by-disorder spectral gap
To determine the anisotropy value that gives a gap comparable to what is seen experimentally we use the expression for the pseudo-Goldstone gap to leading order in 1/S given by [53] where the energy densities are defined through the total ground state energy where N = N c N s is the total number of spins, θ is the uniform tilt of the moments out of the ab plane towards the c-axis. The derivatives in Supplementary Equation (29) are evaluated about the quantum selected ground state configuration. We determine and so Note that the classical energy cl is independent of the anisotropy parameter η. For completeness we note that in the 111 frame, more closely tied to the underlying microscopic superexchange mechanism, the nearest neighbor part of ∂ 2 cl ∂θ 2 0 For η = J ⊥ 1 we find Equation (21) we found ∆S = 0.021, this increases to ∆S = 0.0237 for η = J ⊥ 1 . Such a large value of η comparable to the largest exchange J ⊥ 1 , is however not realistic, as in this case the dispersion relations would not be compatible with the experimental data, based on the analysis of the previous section a value of |η| 1.7 meV would be more in line with the observed dispersions. This suggests that either (i) the order-by-disorder gap formula Supplementary Equation (29) underestimates the gap, being derived in the large-S limit and applied here for S = 1/2, or ii) there are other effects present in the actual material that also contribute to the gap, and we consider such a possible gap generation mechanism via spin-orbital exchange later in Supplementary Note 10. The phase in the lower left part of the phase diagram is parameterized by moment angles φ = 0 and the other phase has φ = π/6. The contours and color scheme show the order-bydisorder pseudo-Goldstone gap computed from Supplementary Equation (29). The phase boundary is the locus of points (solid white dots) where the gap ∆ closes, vertical/horizontal error bars indicate the estimated uncertainty of the gap closing location in vertical/horizontal scans.
To complete this section, in the following we discuss the effects of the other anisotropic exchange terms in Supplementary Equation (22). The mixed in-plane-outof-plane terms J xz and J yz on their own make no contribution to the classical and leading order quantum ground state energy. The effect of an in-plane mixed term J xy can most easily be described by working in a reference frame x y z rotated with respect to the xyz frame around z-axis by an angle φ 0 such that we rotate into the principal axes of and select the axis that corresponds to the minimal eigenvalue. One finds that the angle φ 0 is half the polar angle of the vector (η, −2J xy ). For example, this establishes that for η < 0 and J xy = 0, φ 0 = π/6. Since the model has a 3 point group symmetry and is time reversal invariant the solution is determined by the above equation up to integer multiples of π/3. The zero-point quantum energy contribution has the same form as in Supplementary Equation (28), but with an origin angle φ 0 that depends on the value of J xy as above. In other words, the minimum energy angle is not symmetry-restricted to take only the discrete values 0 or π/6 (up to integer multiplets of π/3), but can take any real value in-between those reference ones for finite J xy . We note that this exchange anisotropy term is symmetry-allowed by the absence at the midpoint of the bond of a mirror plane normal to the nearest-neighbor bond due to (i) the layer stacking and (ii) the buckling of the cobalt honeycombs and (iii) distortions of the CoO 6 octahedra away from regular octahedra.

D. Order-by-disorder in a generalized Kitaev-Γ model
For completeness, we also consider the effects of exchange anisotropy terms described in the 111 frame. Starting from the minimal J 1 −J 6 XXZ model and adding Kitaev K and Γ terms in Supplementary Equation (23) we find that both terms can select the φ = 0 or π/6 (modulo π/3) family of ground states via order-by-disorder according to the phase diagram plotted in Supplementary  Figure 17, where the color represents the zero-point gap ∆ and the white dots indicate the boundary between the two families of ground states.

E. Single-ion anisotropy
Before moving on to discuss the nodal lines, we briefly dwell on the question of whether discrete symmetry breaking can arise directly from single ion anisotropy in CoTiO 3 . In other words, can the observed excitation gap, which implies discrete moment orientations in the ab plane, occur in the absence of bond-dependent anisotropic exchange? The short answer is no because any high order single ion anisotropies have to be filtered through the octahedral crystal field splitting, the spin-orbit coupling and trigonal distortion, with mixing of states provided by exchange interactions, by which point it is more fruitful to view them as effective bonddependent exchange between moments in the lowest Kramers doublet. The longer answer is as follows. The local site symmetry alone (point group 3) constrains the possible single ion crystal-field Hamiltonian to have the form [54] where O m l are Stevens operators acting on the orbital sector, expressed in a reference Cartesian XY Z frame with Z along c and X in the ab plane along some reference direction, unconstrained by symmetry for point group 3 (previously in Supplementary Note 2 we expressed the crystal field Hamiltonian H CF for an ideal octahedron in a reference frame with respect to the four-fold octahedron axes). In terms of the angular momentum operators O 0 2 = 3L 2 z − L(L + 1), which arises from the trigonal distortion of the ideal CoO 6 octahedron. The only operators that have a nontrivial in-plane anisotropy are O 3 4 , O 3 6 and O 6 6 . Semiclassically, for ϑ the polar angle from the Z axis and φ the azimuthal angle from the X-axis in the XY plane, O 3 4 ∼ sin 3 ϑ cos ϑ cos 3φ which vanishes for in-plane moments (ϑ = π/2) as does O 3 6 .
The operator O 6 6 ∼ L 6 − +L 6 + or, semiclassically, cos 6φ, has a 6-fold periodic angular dependence in the ab plane and can in principle lead to in-plane discrete symmetry breaking. However, within the low energy effective l = 1 Γ 4 orbital sector (which is the orbital ground state for an ideal octahedron), this term is inoperative. To obtain an effect from this operator we are forced to extend the model to the full free-ion L = 3 Hilbert space. We may estimate the size of the spectral gap realizing that it can only originate through perturbative mixing of excited orbital levels (of energy G ∼ 1 eV above the l = 1 Γ 4 orbital ground state) via the exchange J (with energy scale ∼ 10 meV). In other words, the mechanism, as previously argued, is order-by-disorder through virtual crystal field fluctuations [26,27]. A strong coupling calculation [27] reveals that the spectral gap will scale as J 3 /G 2 ∼ 10 −3 meV. In order to account for the magnitude of the observed spectral gap (∆ = 1 meV), multiplicative factors including those coming from perturbation theory combinatorics and matrix elements would have to boost this by three orders of magnitude. At the level of the full d 7 Hilbert space, both the single ion anisotropy and the spin-orbital exchange will contribute to the observed gap. To the extent that these mechanisms can be disentangled (since the O 6 6 anisotropy will also affect the exchange) the latter mechanism will be the primary one as this appears straightforwardly from mixing of an order 1 meV spin-orbital exchange across a 10 meV crystal field gap. This mechanism can also be viewed within the effective l = 1, S = 3/2 subspace or within the pseudo-spin one-half picture where the virtual crystal field excitations would manifest themselves as bond-dependent multiple spin exchange anisotropies.

Supplementary Note 9. NODAL LINES TOPOLOGY AND SYMMETRY
We now consider the implications of various discrete symmetries of the spin wave Hamiltonian for the magnon band structure and its topological properties. To consider the full generality of the problem, in the following we work with the (full) four-sublattice magnetic unit cell, for which the spin wave Hamiltonian has the form operators create magnons on the A, B, C and D magnetic sublattices, respectively, and the dynamical matrix takes the block form Without giving the explicit form of the A and B matrices for the J 1 to J 6 XXZ Hamiltonian, we note that the dynamical matrix D 8×8 (k) has time reversal symmetry defined through T −1 D 8×8 (k)T = D 8×8 (−k), where T is antiunitary. The time reversal operator is, explicitly, the complex conjugation operator times the unit operator acting on the sublattice indices. The same Hamiltonian also has spatial inversion symmetry − a unitary transformation: P −1 D 8×8 (k)P = D 8×8 (−k) with P = Γ 1 , the anti-diagonal matrix with ones along the anti-diagonal.
It follows that the model has A ≡ PT symmetry: where A is anti-unitary. It is this symmetry that is responsible for protecting the Dirac nodal lines. This is because the PT symmetry imposes a reality condition on the Hamiltonian. Consider now a closed loop in the 3D Brillouin zone. The Hamiltonian defined along this loop is real and along this loop there is a Z 2 topological classification meaning that there is a winding number that assumes values 0 (topological trivial) and π (topologically nontrivial). If the winding number on the loop is π, it must enclose a singular point -a nodal pointand since the winding number cannot change continuously, say by deforming the loop in 3D, there must be a nodal line in 3D such that the winding number is π on any loop winding around the nodal line. We now examine the robustness of the PT symmetry under perturbations. From the foregoing explicit form of the operators, it is straightforward to see that is the most general spin wave Hamiltonian that preserves PT . In other words, in order to break PT , we must find an exchange coupling that breaks these weak constraints on the Hamiltonian.
The J 1 to J 6 XXZ model has both time reversal symmetry and parity symmetry and is therefore compatible with Supplementary Equations (35) and (36) with a much more restricted form -for example the diagonal elements are identical. If we go beyond the XXZ model to the full set of symmetry-allowed exchange couplings, we have checked that the resulting spin wave Hamiltonian preserves the PT symmetry out to and including the sixth nearest neighbor couplings which includes 38 independent exchange terms. We have additionally verified explicitly (by numerical solution of the corresponding 8×8 spin-wave Hamiltonian) that none of these couplings gap out the nodal lines at least when the magnetic ground state is preserved.
We now consider the structure of the nodal lines. When the model is fine-tuned to have J 2 of Heisenberg form (J ⊥ 2 = J z 2 ) and the third neighbor exchanges are set to zero (J ⊥ 3 = J z 3 = 0), one finds that nodal lines are degenerate, occur at the K points, and are straight lines along L. For anisotropic J 2 (J ⊥ 2 = J z 2 ), the nodal lines around each K-point split into a pair of helical nodal lines that each wind around the K-points as illustrated by the red and blue lines in Fig. 3b) with periodicity L = 3. The chirality of the double helix nodal lines, i.e. the sense in which they wind along the L direction, is opposite at neighboring corner K-points due to the3 symmetry of the lattice. The minimal model that produces these "double helix" nodal lines is the XXZ J 1 -J 2 model. In this model, the anisotropy of the J 2 coupling is responsible for the splitting of the nodal lines and the winding along L, i.e. − the radius of the helix in momentum space in planes perpendicular to the L direction is Note that δk = 0 (the pair of nodal lines are fused and straight) for isotropic J 2 couplings. In this case adding anisotropic J 3 couplings can also produce helical nodal lines, but the effect comes in at higher order with |δk| For completeness, we note that the magnitude of δk is not affected by the buckling parameter , or by adding further neighbor XXZ couplings on the J 4 , J 5 or J 6 bonds.
We also note that the two nodal lines in a pair that wind around each K-point are translated versions of each other along L by the propagation vector Q (translation by 2Q leaves each nodal line invariant). This is most easily seen by first working in the local frame, where a single nodal line appears near each K point (the red helical lines in Fig. 3b) corresponding to the nodal points κ where ω + (κ) = ω − (κ). In the global frame, each of the modes ω ± acquires a 'pair' shifted in momentum by Q (as per Supplementary Equation (19)), so the original nodal points κ are translated into a new set of nodal points κ + Q, giving rise to the blue helical lines in Fig. 3b).
One can understand the origin of the double helix nodal lines in terms of the underlying exchange couplings by first observing that the γ nk functions for n = 1, 4, 5, 6 connect groups of three sites in the same layer and those always appear in the spin-wave Hamiltonian multiplying a function of the form j exp(ik · δ j ) where δ j vectors connect the central site to the three neighbors related by three-fold symmetry. These functions vanish at the hexagonal zone corner (independent of L) implying that all such terms must preserve a nodal line running along the L direction. In our J 1 − J 6 model, the relevant couplings for the existence of a helical nodal line are J 2 and J 3 . If the third neighbor couplings vanish, the helical line appears only if J 2 is anisotropic. Examining the 4 × 4 spin wave Hamiltonian, we observe that the only term that depends on the anisotropic part of J 2 is in the element B. We also note that such a term is identical to the one that appears in the related rhombohedral J 1 − J 2 Heisenberg ferromagnet (J 1,2 < 0), already studied in [55]. In the present cobalt lattice arrangement, the spin-wave Hamiltonian for this ferromagnetic model is where γ jk = γ jk e ik·(r2−r1) with j = 1, 2 and γ jk defined in Supplementary Equation (15). The presence of nodal points is determined by the condition |J 1 γ 1k +J 2 γ 2k | = 0. For J 2 = 0 this gives straight Dirac nodal lines at Kpoints. For finite J 2 the nodal lines become helical and precess around the K-points upon varying L. In detail, considering for concreteness the nodal line near the Kpoint (5/3,2/3), the above condition gives the in-plane wavevector offset δk from K, in Cartesian coordinates as to first order in J 2 /J 1 . Here we have used the Cartesian xyz frame defined in Supplementary Figure 18d), where k x is along (11 2 0) and k y along (010). The above equations describe a helical nodal line winding clockwise in the +L direction with period 3 along L and with radius proportional to J 2 /J 1 .

Supplementary Note 10. SPIN-ORBITAL FLAVOR WAVE THEORY
We have presented detailed, quantitative parametrization of the experimental data within the context of an effective XXZ spin one-half model for the cobalt moment.
The validity of this model can be argued on the basis of single-ion spectrum in the presence of spin-orbit coupling and trigonal distortion of the oxygen octahedra that leave a doublet on each magnetic site separated from the first excited state by about 28 meV. The quality of the agreement with data provide an ex post facto justification for this effective model. However, the exchange scale is some significant fraction of the low-lying single-ion splitting so one is naturally led to consider a more microscopic route to describing CoTiO 3 that incorporates the spin S = 3/2 and effective orbital l = 1 degrees of freedom in full. One important motivation for building on the effective model is to provide a direct calculation of the clearly resolved spectral gap of about 1 meV above the magnetic Bragg peaks. We have shown that the effective spin one-half model including exchange anisotropies cannot open up a gap at the mean field level, while the leading fluctuation contribution to the energy does select a discrete set of ground states. One then infers that a gap will arise in the spin wave spectrum and we have calculated this gap to leading order in 1/S. However, by enlarging the local Hilbert space and considering spin-orbital exchange there is hope that mixing of states does lift this U (1) degeneracy. In this section we develop a flavor wave or multi-boson expansion [56,57] for the excitations in CoTiO 3 that also allows us to include arbitrary couplings in the Hamiltonian and indeed find that the U (1) degeneracy is lifted by spin-orbital exchange terms, which select discrete in-plane angular orientations for the magnetic moments in the ground state and open a gap in the magnon spectrum.
A further motivation for investigating a more microscopic model is to deepen our understanding of the orderby-disorder mechanism. We have investigated in detail the effect of bilinear anisotropic exchange terms acting within the pseudospin one-half subspace. From the point of view of this effective model, any microscopic exchange of higher order acting on two sites projects down to an such a pseudospin one-half exchange coupling to leading order in an expansion in the inverse crystal field gap. However, virtual crystal field fluctuations will generate multi-site effective exchange terms that will break the accidental degeneracy of the spin bilinear model down to a discrete set of ground states. These effects can be more efficiently captured by enlarging the Hilbert space to consider the low-lying spin and orbital coupled crystal field as we describe below.
In this enlarged Hilbert space we shall consider the leading order spin-orbital exchange that produces orderby-disorder through virtual crystal field fluctuations. Higher order terms -such as four spin terms or two-site multipolar exchange -are possible in principle but they are suppressed by powers of the large charge gap. In cases where higher order exchange couplings are significant, such as in cuprates, the reason for the size of these couplings is the combination of proximity to a Mott transition and the largeness of the exchange [58,59]. CoTiO 3 , in contrast, is deep in the Mott insulating regime with an exchange about 5% of that in the cuprates leading to an estimate of the ratio of the biquadratic to the Heisenberg exchange of O(10 −3 ). As we discuss next, there are spinorbital exchange terms that are both much larger than this and that lead to discrete symmetry breaking.

A. Flavor Wave Expansion
The computation of the excitation spectrum proceeds as follows. We consider a general one-and two-body Hamiltonian acting on the spin S = 3/2, orbital l = 1 subspace where it is understood that the one-body terms include the spin-orbit coupling H SO and the trigonal distortion H trig . A local mean field theory yields a spectrum |i, a, p = l z ,S z c (p) l z ,S z |l, l z , S, S z on each distinct mag-netic sublattice a where p = 0, . . . , d − 1 and d = 12 is the local Hilbert space dimension. The ground state can then be written as |Ψ MF = i a |i, a, 0 .
Each operator can be written aŝ The leading order term is simply the mean field energy, the first order terms in A vanish when the mean field energy is minimal and the quadratic contribution from the single-ion physics is: The quadratic terms coming from the interactions are In Fourier space this can be written as is a 2(d − 1)N s × 2(d − 1)N s matrix (N s = 4 being the number of magnetic sublattices). The component matrices are The spectrum is determined by carrying out a Bogoliubov transformation which amounts to diagonalizing the matrix σ 3 H FW (k) where σ 3 = diag(1, . . . , 1, −1 . . . , −1) -the matrix with (d−1)N s ones and (d−1)N s minus ones along the diagonal. The transformation T k that brings the Hamiltonian to diagonal form (T † k ) −1 H FW (k)T −1 k = Λ k is generally non-unitary satisfying instead the constraint T k σ 3 T † k = σ 3 . From this, we compute the inelastic neutron cross section which is proportional to with where µ = g l l + g S S and the g-factor for the effective orbital moment is g l = −3/2 and g S ≈ 2.

B. Spin-Orbital Exchange
In the context of the effective spin one-half model used to fit the data, we considered various types of anisotropic exchange. To nearest neighbor, for example, we saw that six couplings are allowed by symmetry including Kitaev and Γ couplings. In this section, we briefly review how the effective spin one-half model can be understood through a superexchange calculation mediated by cobaltoxygen-cobalt bonds and we use these results to carry out a mean field calculation of the ground states of coupled spin-3/2 effective orbital moments with spin-orbital exchange couplings.
First we consider the geometry of the principal exchange pathways in CoTiO 3 The crystal structure consists of edge-sharing CoO 6 octahedra forming a honeycomb arrangement in the ab plane as illustrated in Supplementary Figure 7 left panel. The cobalt honeycomb network is buckled like cyclohexane in its chair conformation (± signs on the blue balls indicate alternating cobalt positions above/below the plane) so the only rotational symmetry at the cobalt sites is 3-fold and the oxygen geometry breaks the cobalt sublattice out-of-plane mirror symmetries. Exchange interactions between neighboring cobalt ions is primarily mediated by a pair of oxygen ions forming a planar unit with two near 90 • Co-O-Co bonds (at least within the errors of the crystal structure refinement).
For the purpose of identifying the most relevant exchange mechanisms, in the following we consider an idealized version of the actual crystal structure, where the cobalt honeycomb is planar (no buckling) and the oxygen octahedra are regular, as illustrated schematically in Supplementary Figure 18d). The spin-orbital exchange interactions for this idealized Co-O-Co bonding geometry have been considered by Liu and Khaliullin [15] and Sano, Kato and Motome [16]. The calculation of the exchange in this paper was carried out within the set of effective orbital l = 1 and spin 3/2 degrees of freedom on the cobalt ions. The geometry of the low-lying triplet of d orbitals is central to this calculation as hopping through the mediating oxygens with right-angle bonds allows the orbital character to change. For example, in the natural Cartesian frame xyz shown in Supplementary Figure 18d) the d yz orbital on one site is connected via the oxygen to d zx on the other site.
There are several possible exchange processes and hence several distinct spin-orbital exchange couplings. These processes can be grouped into a pair of classes. There are then, naively, nine possible exchange processes obtained from pairs with one from each class. However, processes IIIA and IIIC vanish by symmetry leaving seven routes. The resulting interactions are the product of isotropic exchange in spin space and angular-momentum violating exchange in orbital space. We consider the following four distinct spin-orbital couplings: on a single honeycomb nearest neighbor bond in the frame illustrated in Supplementary Figure 18d) where, following Liu and Khaliullin, we use the notation a = d yz , b = d zx , c = d xy (where xyz are the axes of the Cartesian 111 frame) and n a = a † a etc. Coupling t 1 originates from IA and IC, t 2 from IA, t 3 from IA and t 4 from IIB, IB and IA. While the couplings are fixed by microscopic terms, in the ensuing calculations we treat them as free parameters. Pure Heisenberg spin exchange may also arise microscopically, through process IIIB, and we also include these J n couplings for n = 1 to 6 in H Heisenberg = i,j n J n S i · S j . The above calculations assume an idealized cobaltoxide structural bonding and below we find that the spinorbital exchange terms derived in this case in Supplementary Equation (46) are sufficient to account for the main features of the ground state selection and excitation spectrum. The actual crystal structure has additional distortions, in particular the cobalt network is buckled, which amounts to a roughly 12 • tilt of the cobalt-oxygen-cobalt unit about an axis through the pair of oxygens mediating that bond. Further refinement of the ground state selection and fine structure of the spin wave spectrum may necessitate including the local rotation of the exchange coming from this buckling or perhaps including higher order contributions to the exchange coupling beyond the terms in Supplementary Equation (46). So far we have not discussed the spin-orbit coupling and trigonal distortion. As the hopping and Coulomb scales are the dominant energy scales the superexchange calculation is carried out without them and they are then included on an equal footing with the large U or Hund coupling spin-orbital exchange. To obtain the effective spin one-half exchange model of previous sections of this paper, one may project the spin-orbital exchange onto the spin-orbit coupled trigonally distorted doublet; the anisotropy in the effective spin model is inherited from the angular momentum violating orbital couplings.

C. Mean Field -Flavor Wave Results
As described in the previous section, the microscopic exchange to nearest neighbor acting on the spin 3/2 and effective orbital angular momentum 1 states will couple these degrees of freedom leading to an effective anisotropic exchange within the effective spin one-half model obtained by projecting the exchange onto the single-ion ground state doublet.
We consider a mean field theory including the microscopic single-ion physics, the spin-orbital exchange t n (n = 1, 2, 3, 4) to nearest neighbor and pure spin isotropic exchange coupling the nth neighbor J n for n = 1, . . . , 6, together with single-ion spin-orbit coupling λ and trigonal distortion parameter δ as given in Supplementary Note 2. The collinear, easy plane ordered state of CoTiO 3 with antiferromagnetically coupled layers is obtained straightforwardly by setting J 1 < 0, J 2 > 0. In Supplementary Figure 18a-b) we illustrate the magnon spectrum obtained via the flavor-wave approach for a representative set of exchange parameters chosen such as to approximately reproduce the in-plane and out-of-plane spin wave bandwidths seen in experiments. In the calculation of the dynamical structure factor, as discussed above, we use idealized spin-orbital exchange that omits the effects of buckling of the cobalt honeycombs, but we do include the effects of the buckling on the spin wave intensities using this simplified exchange, i.e. we use the actual cobalt positions in the crystal structure in the calculation of dynamical structure factor. Panel a) shows the case for spin-only exchange, when the ground state energy is independent of the in-plane moments' orientation angle φ and consequently the magnon spectrum has a gapless Goldstone mode, emerging out of the magnetic Bragg peak position (1, 1, 3/2) and the magnon spectrum has double-helix nodal lines as illustrated for the XXZ model in Fig. 3b). Supplementary Figure 18b) shows the case when a finite spin-orbital exchange perturbation t 3 is switched on. This selects the family of φ = 0 modulo π/3 ground state moment orientations and consequently opens a gap in the magnon spectrum. In addition to capturing the discrete ground state selection and spectral gap, a further advantage of the flavor wave picture is that it also gives the spectrum of exciton modes. Supplementary Figure 18c) shows the obtained dynamical structure factor for the lowest-energy exciton modes, the calculated spectrum bears strong resemblance to the data in Fig. 4c). For the parameters used in the above calculations, the largest effect of the finite spinorbital exchange is in opening of a magnon spectral gap, the magnon spectrum still displays nodal lines, and the effect on the exciton modes is relatively small.
We note that in the present treatment of the spinorbital exchange, each of the four t 1−4 terms in Supple-mentary Equation (46), irrespective of their sign, selects the family of φ = 0 modulo π/3 ground states at least when their magnitude is compatible with the experimentally observed spin wave gap. We leave it as subject for future research to investigate whether this is true in general for these couplings and, if so, how to understand this perturbatively in the spin-orbital exchange coupling over the crystal field gap. We also leave for the future the question of whether other symmetry-allowed spin-orbital exchange terms not explicitly listed in Supplementary Equation (6), could select the alternative set of φ = π/2 modulo π/3 ground states.
We established in Supplementary Note 8 that the effective spin one-half model at the mean field level has an accidental U (1) degeneracy that is lifted by quantum fluctuations. In contrast, the spin-orbital exchange model discussed in this section can be viewed as an example of order arising from virtual crystal field fluctuations first discussed in the context of Er 2 Ti 2 O 7 [26,27]. The theory developed in this section is based around a mean field theory that omits the order-by-disorder corrections discussed in Supplementary Note 8 that act within the effective spin one-half set of states. Yet the accidental degeneracy, that is present in the spin-orbital model when projected down to the low-energy doublets on each site, is lifted within the full mean field theory leading to a discrete set of ground states. The discrete symmetry breaking in this case originates from the enlarged Hilbert space and the admixing of excited crystal field levels into the ground state and is therefore suppressed in powers of the inverse crystal field gap. In CoTiO 3 one expects that both order-by-disorder mechanisms are operative. While disentangling the relative contributions of the two effects is non-trivial, the fact that the exchange scale is a significant proportion of the crystal field gap in the material strongly suggests that order by virtual crystal field fluctuations is an important factor in the ground state selection in the system. o S w 6 e W Y K K Y z Y r I G C t M j C 2 q Y k v w l r + 8 S j o X d c + t e 3 e X t e Z 1 U U c Z T u A U z s G D B j T h F l r Q B g I K n u E V 3 p w n 5 8 V 5 d z 4 W o y W n 2 D m G P 3 A + f w A w W Z L y < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " t k N Q X I x Z Y o a D 4 / g g 3 O h 4 R Z F r z 6 E = " > A A A B 9 X i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I U J d F N y 4 r 2 A e 0 Y 8 l k M m 1 o J h m S j F K G / o c b F 4 q 4 9 V / c + T d m 2 l l o 6 4 G Q w z n 3 k p M T J J x p 4 7 r f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 5 1 t E w V o W 0 i u V S 9 A G v K m a B t w w y n v U R R H A e c d o P J T e 5 3 H 6 n S T I p 7 M 0 2 o H + O R Y B E j 2 F j p Y R B I H u p p b K 8 M z 4 b V m l t 3 5 0 C r x C t I D Q q 0 h t W v Q S h J G l N h C M d a 9 z 0 3 M X 6 G l W G E 0 1 l l k G q a Y D L B I 9 q 3 V O C Y a j + b p 5 6 h M 6 u E K J L K H m H Q X P 2 9 k e F Y 5 9 H s Z I z N W C 9 7 u f i f 1 0 9 N d O V n T C S p o Y I s H o p S j o x E e Q U o Z I o S w 6 e W Y K K Y z Y r I G C t M j C 2 q Y k v w l r + 8 S j o X d c + t e 3 e X t e Z 1 U U c Z T u A U z s G D B j T h F l r Q B g I K n u E V 3 p w n 5 8 V 5 d z 4 W o y W n 2 D m G P 3 A + f w A w W Z L y < / l a t e x i t > c < l a t e x i t s h a 1 _ b a s e 6 4 = " U C K F p 3 f n w 4 Q E s e T v D N q u 6 n Z y X 6 0 = " > A A A B 9 X i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I U J d F N y 4 r 2 A e 0 Y 8 l k M m 1 o J h m S j F K G / o c b F 4 q 4 9 V / c + T d m 2 l l o 6 4 G Q w z n 3 k p M T J J x p 4 7 r f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 5 1 t E w V o W 0 i u V S 9 A G v K m a B t w w y n v U R R H A e c d o P J T e 5 3 H 6 n S T I p 7 M 0 2 o H + O R Y B E j 2 F j p Y R B I H u p p b K + M z I b V m l t 3 5 0 C r x C t I D Q q 0 h t W v Q S h J G l N h C M d a 9 z 0 3 M X 6 G l W G E 0 1 l l k G q a Y D L B I 9 q 3 V O C Y a j + b p 5 6 h M 6 u E K J L K H m H Q X P 2 9 k e F Y 5 9 H s Z I z N W C 9 7 u f i f 1 0 9 N d O V n T C S p o Y I s H o p S j o x E e Q U o Z I o S w 6 e W Y K K Y z Y r I G C t M j C 2 q Y k v w l r + 8 S j o X d c + t e 3 e X t e Z 1 U U c Z T u A U z s G D B j T h F l r Q B g I K n u E V 3 p w n 5 8 V 5 d z 4 W o y W n 2 D m G P 3 A + f w A z Y 5 L 0 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " U C K F p 3 f n w 4 Q E s e T v D N q u 6 n Z y X 6 0 = " > A A A B 9 X i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I U J d F N y 4 r 2 A e 0 Y 8 l k M m 1 o J h m S j F K G / o c b F 4 q 4 9 V / c + T d m 2 l l o 6 4 G Q w z n 3 k p M T J J x p 4 7 r f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 5 1 t E w V o W 0 i u V S 9 A G v K m a B t w w y n v U R R H A e c d o P J T e 5 3 H 6 n S T I p 7 M 0 2 o H + O R Y B E j 2 F j p Y R B I H u p p b K + M z I b V m l t 3 5 0 C r x C t I D Q q 0 h t W v Q S h J G l N h C M d a 9 z 0 3 M X 6 G l W G E 0 1 l l k G q a Y D L B I 9 q 3 V O C Y a j + b p 5 6 h M 6 u E K J L K H m H Q X P 2 9 k e F Y 5 9 H s Z I z N W C 9 7 u f i f 1 0 9 N d O V n T C S p o Y I s H o p S j o x E e Q U o Z I o S w 6 e W Y K K Y z Y r I G C t M j C 2 q Y k v w l r + 8 S j o X d c + t e 3 e X t e Z 1 U U c Z T u A U z s G D B j T h F l r Q B g I K n u E V 3 p w n 5 8 V 5 d z 4 W o y W n 2 D m G P 3 A + f w A z Y 5 L 0 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " U C K F p 3 f n w 4 Q E s e T v D N q u 6 n Z y X 6 0 = " > A A A B 9 X i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I U J d F N y 4 r 2 A e 0 Y 8 l k M m 1 o J h m S j F K G / o c b F 4 q 4 9 V / c + T d m 2 l l o 6 4 G Q w z n 3 k p M T J J x p 4 7 r f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 5 1 t E w V o W 0 i u V S 9 A G v K m a B t w w y n v U R R H A e c d o P J T e 5 3 H 6 n S T I p 7 M 0 2 o H + O R Y B E j 2 F j p Y R B I H u p p b K + M z I b V m l t 3 5 0 C r x C t I D Q q 0 h t W v Q S h J G l N h C M d a 9 z 0 3 M X 6 G l W G E 0 1 l l k G q a Y D L B I 9 q 3 V O C Y a j + b p 5 6 h M 6 u E K J L K H m H Q X P 2 9 k e F Y 5 9 H s Z I z N W C 9 7 u f i f 1 0 9 N d O V n T C S p o Y I s H o p S j o x E e Q U o Z I o S w 6 e W Y K K Y z Y r I G C t M j C 2 q Y k v w l r + 8 S j o X d c + t e 3 e X t e Z 1 U U c Z T u A U z s G D B j T h F l r Q B g I K n u E V 3 p w n 5 8 V 5 d z 4 W o y W n 2 D m G P 3 A + f w A z Y 5 L 0 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " U C K F p 3 f n w 4 Q E s e T v D N q u 6 n Z y X 6 0 = " > A A A B 9 X i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I U J d F N y 4 r 2 A e 0 Y 8 l k M m 1 o J h m S j F K G / o c b F 4 q 4 9 V / c + T d m 2 l l o 6 4 G Q w z n 3 k p M T J J x p 4 7 r f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 5 1 t E w V o W 0 i u V S 9 A G v K m a B t w w y n v U R R H A e c d o P J T e 5 3 H 6 n S T I p 7 M 0 2 o H + O R Y B E j 2 F j p Y R B I H u p p b K + M z I b V m l t 3 5 0 C r x C t I D Q q 0 h t W v Q S h J G l N h C M d a 9 z 0 3 M X 6 G l W G E 0 1 l l k G q a Y D L B I 9 q 3 V O C Y a j + b p 5 6 h M 6 u E K J L K H m H Q X P 2 9 k e F Y 5 9 H s Z I z N W C 9 7 u f i f 1 0 9 N d O V n T C S p o Y I s H o p S j o x E e Q U o Z I o S w 6 e W Y K K Y z Y r I G C t M j C 2 q Y k v w l r + 8 S j o X d c + t e 3 e X t e Z 1 U U c Z T u A U z s G D B j T h F l r Q B g I K n u E V 3 p w n 5 8 V 5 d z 4 W o y W n 2 D m G P 3 A + f w A z Y 5 L 0 < / l a t e x i t >