Designing three-dimensional flat bands in nodal-line semimetals

Electrons with large kinetic energy have a superconducting instability for infinitesimal attractive interactions. Quenching the kinetic energy and creating a flat band renders an infinitesimal repulsive interaction the relevant perturbation. Thus, flat band systems are an ideal platform to study the competition of superconductivity and magnetism and their possible coexistence. Recent advances in the field of twisted bilayer graphene highlight this in the context of two-dimensional materials. Two dimensions, however, put severe restrictions on the stability of the low-temperature phases due to enhanced fluctuations. Only three-dimensional flat bands can solve the conundrum of combining the exotic flat-band phases with stable order existing at high temperatures. Here, we present a way to generate such flat bands through strain engineering in topological nodal-line semimetals. We present analytical and numerical evidence for this scenario and study the competition of the arising superconducting and magnetic orders as a function of externally controlled parameters. We show that the order parameter is rigid because the quantum geometry of the Bloch wave functions leads to a large superfluid stiffness. Using density-functional theory and numerical tight-binding calculations we further apply our theory to strained rhombohedral graphite and CaAgP materials.


I. INTRODUCTION
The study of correlated many-particle states in flatband systems goes back to the consideration of fewparticle nuclear-physics systems in the 1960s when S. T. Belyaev demonstrated that, in the presence of degenerate single-particle states, interactions can lead to a pairing gap increasing linearly with the interaction strength 1 . Since then, there has been a fruitful exchange of ideas between the nuclear-physics and the condensed-matter communities exploring analogies between nuclear-physics systems and ultrasmall superconducting grains 2-4 . Recently, two-dimensional (2D) flat bands have provided new ground for exotic states of condensed matter 5-9 . In particular, the advances in the fabrication of flat bands in twisted bilayer graphene have attracted a lot of attention due to the novel exotic phases becoming available experimentally [10][11][12][13][14][15][16] . Other setups and materials realizing 2D flat bands have been studied [17][18][19] but, to our knowledge, their three-dimensional (3D) counterparts 20,21 have not been explored in realistic materials. In the present manuscript, we carry this development to its logical end point by proposing a viable approach for realizing 3D flat-band systems.
We combine two actively-studied ingredients to manifest a 3D flat band: nodal-line semimetals (NLSMs) 22 and strain engineering in topological semimetals 23 . NLSMs are materials where topologically protected band crossings form a line (nodal line) in the Brillouin zone, provided that certain symmetries are satisfied. Additionally, topologically protected drumhead surface states ap-pear inside the region bounded by the projection of the nodal line onto the 2D surface Brillouin zone. The drumhead surface states exhibit a nearly flat dispersion if the material has approximate chiral symmetry, for example due to a sublattice structure 7 .
Strain engineering produces opportunities to generate Landau level-like flat bands in the absence of external magnetic fields 24-30 . This is because the action of the strain near the Fermi surface resembles that of a magnetic field in the local area of the Brillouin zone. In the case of Weyl/Dirac semimetals, there are simply two Fermi pockets near the Fermi energy and the action of strain near the pockets can be described as electromagnetic fields 31-33 . Here, we show that such properties of strain are more generic and can be applied to materials with nodal lines as well. For stationary strain, the resulting pseudo-magnetic field then depends on the position along the nodal line. There is, however, always one commonality in nodal systems insensitive to the strength and direction of the pseudo-magnetic field -the zeroth pseudo-Landau level (PLL). Therefore, one would expect that the zeroth PLL forms a 3D flat band, while the higher PLLs are flat only in two dimensions and have a dispersion along the nodal line.
We confirm this intuitive argument above using numerical and analytical arguments. We show that the zeroth PLL indeed forms a 3D flat band, which evolves from the drumhead surface states of the NLSM, and obtain wavefunctions thereof. Using these wavefunctions and assuming competition between magnetism and superconductivity, we obtain the phase diagram of the system as a function of the filling factor and the interaction strengths.

arXiv:2008.09122v2 [cond-mat.mes-hall] 26 Mar 2021
We show this system is a promising platform for studying intertwined phases 9,34 as it can be tuned in-situ from the magnetic to the superconducting phase by controlling the magnitude of the strain.
Moreover, we go beyond the mean-field picture by analyzing the properties of the collective modes and argue that the emerging order is more stable in the 3D case than in the previously studied 2D cases. Importantly, we show that the system supports a hitherto unexplored 3D quantum geometry of the Bloch wave functions that leads to a large superfluid stiffness in all directions despite the flatness of the bands. The 3D quantum geometry is a nontrivial consequence of the momentum dependence of the direction of the pseudo-magnetic field along the nodal line. Finally, we study the properties of 3D flat bands for the material examples of NLSMs belonging to CaAgP material class and of rhombohedral graphite.

A. Two-band NLSM Hamiltonian and drumhead surface states
A NLSM is a material for which the energy gap between two bands near the Fermi level closes along a line in the Brillouin zone [ Fig. 1(a)]. The minimal model for a NLSM, believed to be applicable, for instance, to the CaAgP materials class (see below) and Ca 3 P 2 35-37 , reads: i=x,y,z cos k i + 2t 2 σ x sin k z . (1) Here σ x,y are Pauli matrices, and t, t 1 , and t 2 are parameters determining the size of the nodal line and the Fermi velocity. Throughout this paper, the length scales are given in units of the lattice constant and the momentum is measured in units of the inverse lattice constant. For t 1 t, this Hamiltonian has a gap closing along a circular path in the Brillouin zone -a nodal loop -given by k z = 0 and (k x , k y ) = Q = Q(cos θ, sin θ), where Q = t 1 /t and θ is the azimuthal angle. The crossing between the two bands is protected by x − y plane reflection symmetry. The bands are σ z eigenstates in the k z = 0 plane corresponding to reflection eigenvalues ±1. Taking a different viewpoint, the nodal line can also be considered to be protected by the chiral symmetry σ y or the combination of time-reversal and inversion symmetries. For this Hamiltonian, the drumhead surface states, appearing at momenta k < Q, have zero energy due to chiral symmetry. In the limit of small nodal-line radii Q, the corresponding surface wave functions at momentum k = Q + q can be written as Ψ S,± Q, q = 1 2L x L y l S sinh L l S 1 ±i e i(Q+q)·r e ± z l S , where q = q(cos θ, sin θ) and q < 0 describes the deviation from the nodal line. Here, ± correspond to top and bottom surfaces, −L/2 < z < L/2, l S = −t 2 /(tQq) is the localization length of the drumhead surface states, and L x (L y ) is the length of the system in the x (y) direction along which we have applied periodic boundary conditions.

B. Strain and pseudo-Landau levels
A simple way to realize strain with a constant gradient of the strain field is bending 23,38 , as shown in Fig. 1b. The bend can be achieved by a proper choice of the growth substrate or by mechanical bending of the device. Such strain creates a displacement field u = (xz/R, 0, 0) with corresponding strain tensor components u 11 = z/R and u 13 = u 31 = x/(2R). In the limit of small nodalline radii, the strain induces a pseudo-magnetic field (see Appendix below) where Q = Q(cos θ, sin θ) is a momentum along the nodal line and l B = √ RQ is the pseudo-magnetic length. In this case, the analytical solution for the zeroth PLL wave function is Here ξ q = (z + z q )/l z , z q = ql 2 B and the pseudo-magnetic localization length in the z-direction is given by Alternatively, one can perform a gauge transformation such that the wave function is localized in the x−y plane and a plane wave in the z-direction. In that case, the localization length is The prefactors in the localization lengths arise due to the elongations of the elliptical cyclotron orbits. Going beyond the lowest-order expansion in Q would lead to a θ dependence of the strength of the pseudomagnetic field because the strain shrinks the nodal line anisotropically. Nevertheless, the energy of the zeroth PLL is independent of the strength of the pseudomagnetic field and, hence, in all cases we robustly obtain the desired 3D flat band. Since the pseudo-magnetic field originates from the spatial variation of the radius of the nodal loop, it can also be created by varying the relative amount of P and As contents in the CaAgP 1−x As x alloys along the z-direction, for example, by changing the crystal growth conditions with time.
L a n d a u le v e ls drumhead states Figure 1. (a) Schematic representation of a nodal-line semimetal bandstructure. Here, Q is a vector pointing to a momentum on the nodal line such that |Q| = Q is the radius of a circular nodal loop. A Dirac cone is formed with respect to the perpendicular momentum components q and kz. (b) Example of a strain profile leading to the formation of a 3D flat band. The strain can be created e.g. by bending the sample or growing it on a cylindrical surface. R is the radius of the cylinder in the middle of the sample. (c) Bulk-boundary correspondence of the zeroth PLL of the NLSM. Due to the strain, the radius of the nodal circle varies as a function of z, such that the radius of the momentum-space area of the drumhead surface states on the top surface Qtop is different from the bottom surface Q bot . Zeroth PLL bulk states appear in the momentum space region between Qtop and Q bot . (d) Spectrum of the model Hamiltonian H tb [Eq. (1)] in the presence of strain [Eq.
(2)] with parameters t1 = 0.25t, t2 = 0.8t, L = 1000, and R = 8000. For comparison, we also show the energies of the effective Hamiltonian H eff (orange bands) constructed from the analytical low-energy solutions (see Appendix below). (e) Numerical (blue) and analytical (orange) wavefunctions of the zeroth PLL at different momenta kx = ki (i = 1, 2, 3, 4) indicated in (d). The drumhead surface states appear at both surfaces for |kx| < Qtop, Q bot . At |kx| = Qtop the drumhead surface states from the top surface are deformed into the zeroth PLL bulk states and move towards the bottom surface with increasing Qtop < |kx| < Q bot . At |kx| = Q bot the PLL bulk state hybridizes with the drumhead surface state at the bottom surface. (f) Degeneracy of the flat band as a function of the height L of the sample. Numerically, the degeneracy is calculated by integrating the density of states over an energy window (|E| < 10 −3 t). The model parameters t1,2 are the same as in (d). Bold lines are the numerical results and dashed lines correspond to the analytical formulas. In the presence of strain, the degeneracy grows linearly with L demonstrating the existence of a 3D flat band with degeneracy N flat proportional to the volume of the sample (green dashed line). In the absence of strain, the degeneracy of the zero-energy drumhead surface states saturates with increasing L demonstrating that it is proportional to the area of the surface of the sample, N0.

C. Connection of drumhead surface states and pseudo-Landau level
In Fig. 1(c, d), we plot the spectrum of the lowest PLL states. We find that the zeroth PLL is two-fold degenerate: the bulk states discussed above coexist with drumhead states localized to the bottom surface. We further observe that the PLL bulk states evolve from drumhead states localized to the top surface, similar to the connection between Fermi arcs and PLLs in Weyl semimetals 39 . By using approximate analytic solutions for the drumhead surface states and for the PLL bulk wave functions, and by considering their coupling at the surface of the sample, we obtain an effective low-energy Hamiltonian (see Appendix below) that accurately describes the ex-act numerical results as shown in Fig. 1(d, e). We confirm the 3D nature of the flat band by plotting the scaling of the total number of states in the flat band for different thicknesses L with and without strain in Fig. 1(f). For sufficiently large L, the number of states for the drumhead surface states without strain saturates to the value N 0 = Q 2 L x L y /(2π). On the contrary, the number of states, including the PLL bulk states, in the presence of strain grows linearly with L as N flat = QL x L y L/(2πl 2 B )confirming our hypothesis of a flat band that is genuinely 3D.
The interaction effects depend on the k-space extent of the flat bands. Unlike in twisted bilayer graphene, the flat band in these 3D systems covers only a fraction of the Brillouin zone. In twisted bilayer graphene, however, the area of the moiré-Brillouin zone is inversely proportional to the square of the superlattice lattice constant and the latter is large around the magic angle. Thus, we expect that interaction effects in the considered 3D flat band systems are as important as in the case of twisted bilayer graphene (see below).

III. EFFECTS OF INTERACTIONS
In flat-band systems, the density of states diverges and, therefore, interaction effects are important. Moreover, these systems have instabilities both with respect to repulsive and to attractive interactions, so that competition between various symmetry-broken phases is expected to be a generic feature of flat-band systems. We now turn our attention to the effect of an attractive pairing interaction and Coulomb repulsion between electrons.

A. Magnetism
For the study of magnetism, we self-consistently solve the mean-field equations for the magnetic order parameter m z and the chemical potential µ in the presence of the density constraint that the filling factor of the band be fixed to C (see Appendix below) Here, V 0 is the effective interaction strength and n F is the Fermi function. The filling factor C is restricted to 0 ≤ C ≤ 2, such that C = 0 corresponds to the situation where the flat-band states are completely empty. whereas for C = 2 both spin-up and spin-down flat-band states are fully occupied. The effective interaction strength V 0 is computed by projecting the Coulomb interaction to the zeroth PLL. Furthermore, it can be tuned by varying the strain (see Appendix below). The critical temperature for magnetism depends on C as For conservative model parameters we estimate T c,m = 3 K (see Appendix below). We point out that this theory describes several different types of magnetic order parameters as their projections to the PLL bulk wave functions are the same (see Appendix below). However, the surface effects break the symmetry explicitly and distinguish the magnetic order parameters from each other. We find that the lowest energy state corresponds to the situation where the magnetic order is staggered with respect to the PLL bulk states and the drumhead states at the bottom surface (see Appendix below). We note that the exact nature of the magnetic order is not important for the discussion of competing phases below since we focus on the bulk states here.

B. Superconductivity
Starting from the reduced BCS Hamiltonian with the density constraint, the mean field equations for the superconducting order parameter ∆ and the chemical potential µ are (see Appendix below) where G 0 is the effective interaction strength and β = 1/(k B T ). The critical temperature is given by For typical model parameters we estimate T c,sc = 1 K (see Appendix below). We note that this value is similar to the critical temperatures observed in twisted bilayer graphene 10,11 . In the case of 3D flat bands in NLSMs, however, the stronger suppression of order-parameter fluctuations and the greater variability of parameters with strain can potentially lead to even larger critical temperatures. The same applies to the magnetic phase described above.

C. Competing phases
As recent experiments in twisted bilayer graphene indicate, the competition between different correlated phases is a generic feature of flat-band systems 15,16 . Moreover, in the case of competition between flat-band magnetism and superconductivity doping the system away from halffilling is generically expected to favor the superconducting phase 8 . Thus, we expect also magnetic and superconducting phases in 3D flat-band systems to be tunable using controllable parameters such as doping and the effective interaction strengths. The latter can be controlled, for example, by strain or an electrostatic environment.  6)] as a function of the filling factor C in the typical situation where the effective interaction strength for superconductivity is smaller than for magnetism G0 = 0.75V0. At half-filling C = 1, the system realizes magnetic order but doping away from half-filling leads to superconductivity, resembling the phase diagrams experimentally observed in cuprates and twisted bilayer graphene. (b) Phase diagram as a function of C and G0/V0 assuming that the system realizes the phase with the larger critical temperature. The dashed horizontal line represents a path through the phase diagram corresponding to (a). This expectation is confirmed by our calculations shown in Fig. 2. In Fig. 2(a), we plot the critical temperatures for magnetism [Eq. (4)] and superconductivity [Eq. (6)] as a function of the filling factor C. At halffilling C = 1, the system realizes the phase with the larger effective interaction strength. Hence, for the typical situation where the Coulomb repulsion dominates the pairing interaction, V 0 > G 0 , the magnetic phase is realized. On the contrary, doping the system away from half-filling favors the superconducting phase. Therefore, flat-band systems are expected to have phase diagrams similar to those observed experimentally in cuprates 40 and twisted bilayer graphene 10,11 . In contrast to the complicated, strongly correlated Mott physics of cuprates 40 , here the mechanisms of superconductivity and magnetism are independent of each other 9 . This is in agreement with experiments on twisted bilayer graphene 15,16 . In addi-tion to the filling factor C, the ground state can also be controlled by the relative interaction strength G 0 /V 0 , as shown in Fig. 2(b).
We note that our proposal offers a degree of control over G 0 and V 0 in a single sample. By varying the strength of the strain one changes the parameters as derived in the Appendix below, thus allowing to tune in situ between magnetic and superconducting phases.
We note that, in reality, the phase diagram may be more complicated due to possible co-existence of the order parameters near the mean-field phase boundary 9 . This may be analyzed using numerical techniques or functional renormalization group and are beyond the scope of the current manuscript. Nevertheless, we underline that 3D flat bands can shed light on the competition and intertwining of order parameters while avoiding the complications caused by strongly correlated Mott physics and the fluctuations present in low-dimensional systems.

D. Superfluid stiffness
To further substantiate our claims concerning the stability of the mean-field solutions, we analyze the collective modes of the system. In the case of 3D flat bands, the quasiparticle spectrum is fully gapped which means that the amplitude mode is gapped. Thus, the order parameter is stable against amplitude-mode fluctuations. This is in strong contrast to the order parameter appearing in the case of flat-band superconductivity due to 2D drumhead surface states, where the system is gapless and susceptible to strong amplitude-mode fluctuations 41 .
The phase rigidity, on the other hand, needs to be analyzed more carefully as kinetic contribution to it is negligible due to the flatness of the band. It is determined by the superfluid stiffness tensor D s , which is related to the supercurrent j in the system as where ϕ is the phase of the superconducting order parameter. The question whether a supercurrent even exists in the systems considered here is particularly relevant as the Fermi velocity within a featureless flat band is zero. Previous studies have shown that in 2D systems there can still exist a nonzero contribution to the superfluid stiffness caused by the quantum geometry of the Bloch wave functions 6,42-48 . However, the possibility of an analogous 3D quantum geometry that gives rise to superfluid stiffness in all directions of a 3D system has not been explored. In the Appendix below, we combine methods developed in the context of quantum Hall systems 6,42 and flat-band superconductors 43,46 to calculate the superfluid stiffness for the model Hamiltonian [Eq.
(2)]. We find that the superfluid stiffness perpendicular D s,zz and parallel D s,xx = D s,yy to the layers have large geometric contributions Therefore, the system has rigidity against phase fluctuations and supports large bulk supercurrents thereby confirming the 3D nature of the superconducting state.
Here, n 0 = Q/(2πl 2 B ) is the degeneracy of the flat band per volume and the directional dependence of the stiffness arises because of the elongations of the semiclassical cyclotron orbits.
We emphasize that the non-zero superfluid stiffness along all three directions, as obtained in this section, is not a property that flat bands exhibit by default, but the necessary 3D quantum geometry arises due to the variation of the pseudomagnetic field direction along the nodal line [Eq.

IV. MATERIAL CONSIDERATIONS
In this section, we provide a materials perspective on our proposal. In particular, we identify CaAgP and rhombohedral graphite NLSMs as possible material platforms for 3D flat bands with negligible spin-orbit interaction. We note that NLSMs can also be realized in systems with sizable spin-orbit interactions. Theoretically proposed candidates are TlTaSe 2 49 or alloys of SnTe 50 .

A. Strained CaAgP
CaAgP crystallizes in an ZrNiAl-type hexagonal structure 51 , which is illustrated in Fig. 3(a). Standard DFT calculations predict a topological NLSM phase for the nodal-line Dirac semimetal CaAgP. However, standard DFT is known to overestimate the band inversion and, in fact, experimental results for CaAgP show a trivially gapped phase for this compound. Nonetheless, it has been shown both theoretically and experimentally that strain and As doping can turn the phase from trivial to topological 52 . In the following, we focus on the NLSM phase of the CaAgP materials class.
The low-energy band structure is formed by the 5s orbitals of the Ag atoms and the 3p orbitals of the light P atoms. Therefore, spin-orbit interaction in CaAgP is small. In the atomic limit, the p orbitals are occupied and the s orbitals are unoccupied. At the Γ and at the A point of the Brillouin zone, there is a band inversion between one s band and one p z band, while the band orderings at M and K are trivial. As a consequence, a line node is observed along the paths Γ-M and Γ-K, which forms a circle in the k x -k y plane centred at the Γ point [see Fig. 3(a) and Appendix below]. The band dispersion at the line node is linear along both the radial and the k z directions. It should be noted that this line node is generally not protected from spin-orbit interaction but, in this case, the weak spin-orbit coupling induces only a small gap on the order of 10 K 51 .
We construct an effective two-orbital tight-binding model (see Appendix below) based on one 3p z -orbital Wannier function (WF) centred at one P atom and one 5s-orbital WF centred within a triangle of Ag atoms as previously done for other triangular systems 53 . With an energy cut-off of 6 meV for the tight-binding parameters, this allows us to accurately capture the low energy physics with a relatively simple model. Most importantly, our two-band model correctly reproduces the NLSM phase of the CaAgP materials class [see Fig. 3 In the next step, we implement the strain terms by modifying the tight-binding parameters in a way similar to our minimal model above (see Appendix below). We choose the bending direction to be aligned with the x axis. We compute the energy spectrum for a (001) slab and find two nearly flat bands, as shown in Fig. 3(c). Similar to our minimal model, these bands consist of superpositions of bottom-surface states and PLL bulk states. In particular, the center of the PLL bulk states shifts as function of the flat-band momentum, which happens symmetrically around the center of the flat band at Γ. However, due to the broken electron-hole symmetry in this material, the two flat bands acquire a small dispersion. Furthermore, we observe a small energy splitting between the two bands due to chiral-symmetry breaking, a symmetry that was present in the minimal model above due to its simplicity. The corresponding gap is smaller in the graphite case we turn to now.

B. Strained rhombohedral graphite
Another possible candidate material is rhombohedral graphite. It consists of light carbon atoms and has, therefore, negligible spin-orbit interaction. Rhombohedral graphite is composed of ABC-stacked graphene sheets [see Fig. 3(d)]. Isolated graphene sheets feature Dirac cones at the K and K' points of their 2D Brillouin zone. Stacking several sheets in the z direction in an ABCtype fashion leads to a k z -dependent shift of the Dirac cones away from the K and K' points thereby forming a nodal line. In particular, the nodal lines in rhombohedral graphite spiral around the vertical hinges of its hexagonal Brillouin zone [see Fig. 3 The low-energy properties of rhombohedral graphite, similar to graphene, can be captured by a two-band model (see Appendix below). This model reproduces the nodal spirals as can be inferred from comparing the spectrum shown in Fig. 3(e) to the Brillouin zone sketch in Fig. 3(d). Note that the electron-hole symmetry is only weakly broken in this material.
As before, we implement strain terms corresponding to a cylindrical strain profile with the bending direction aligned with the x axis (see Appendix below). The energy spectrum of a (001) slab of this system features two flat bands at the Fermi level, as illustrated in Fig. 3 contrast to CaAgP, their splitting and bandwidth is negligibly small. This is due to the approximate electron-hole symmetry of the material. The flat bands are again composed of PLL bulk states and drumhead surface states. However, their behavior and composition as a function of the flat-band momentum (q x , q y ) differs from our previous observations. Depending on the direction, we observe all possible combinations: bottom-surface and PLL bulk state, PLL bulk and top-surface state, bottom-and topsurface state, and even two PLL bulk states with opposite shift behavior (see Appendix below). Despite these differences we emphasize that in all cases a 3D flat band appears robustly in the presence of strain.

V. DISCUSSION AND CONCLUSIONS
We have proposed a feasible way to create 3D flat bands by applying strain to a nodal-line semimetal. In the process, we have discovered an inherent connection between the arising 3D flat band and the drumhead surface states of the parent nodal-line semimetal. Moreover, we have investigated the effects of interactions on the flat bands, highlighting the competition of superconductivity and magnetism in analogy with twisted bilayer graphene. By computing the superfluid stiffness, we have confirmed that the flatness of the bands does not impede the phase rigidity of the system. The flat-band system thus supports supercurrents and true long-range order solely due to the 3D quantum geometry of the Bloch wave functions, which arises because of the variation of the pseudomagnetic field direction along the nodal line. Going beyond our general idea, we have applied our theory to the potential candidate materials CaAgP and rhombohedral graphite. We have shown that both materials give rise to sufficiently flat bands under experimentally accessible strain, therefore representing viable candidate materials for the experimental realization of our ideas. Our conservative estimates for the critical temperatures are on the order of few Kelvins. However, we emphasize that they depend strongly on the magnitude of the applied strain or the pseudomagnetic field. In particular, if the pseudomagnetic field is controlled with the chemical composition we expect that much stronger pseudofield strengths can be obtained, such that critical temperatures on the order of tens of Kelvins could potentially be achieved. The critical temperatures are proportional to the k-space extent of the flat bands. Therefore, it is important to note that, even though the portion of the Brillouin zone covered by the flat band is small under the assumptions Q 1 and R L, the flat bands can generally cover an arbitrary portion of the Brillouin zone. In particular, the first assumption Q 1 was made only for the sake of analytical transparency. Moreover, implementing the effective magnetic field by varying the chemical composition so that the radius of the nodal line changes along the z-direction, we may also violate the second condition R L. In the extreme case Q ∼ 1 and R ∼ L we find numerically that the 3D flat bands can indeed cover most of the Brillouin zone (see Appendix below).
We have considered a specific geometry, realizable by bending the sample, leading to a suitable strain profile. We have, however, explicitly checked that only the strain tensor component u 11 matters for the appearance of the 3D flat band (see Appendix below) and, therefore, the only important ingredient for realizing a 3D flat band in nodal-line semimetals is an inhomogeneous tensile strain profile across the sample. Thus, a suitable strain tensor can be realized also in various other ways, such as by utilizing the type of multilayer structure designed in Ref. 54 for the creation of a giant linear strain gradient.
We emphasize that if the strain is realized by bending the sample, the different phases and phase transitions can be studied in situ in this system. As discussed above, tuning the strain changes the relative interaction strengths of superconductivity and magnetism thereby allowing to study the transition between these phases. Moreover, it is known that flat bands are an interesting playground for Anderson localization 55 . The in-situ tunability of the flat-band degeneracy by strain controls also the effective disorder strength, so that it might be possible to study the Anderson transition in a single sample.
In condensed matter physics, there is the immediate interest for the experimental community to realize this 3D analogue of twisted bilayer graphene in order to study the interplay of quantum geometry, flatness of the dispersion, disorder and intertwining of different types of order. Moreover, our proposal of realizing 3D flat bands by strain engineering nodal-line semimetals opens a new frontier of research beyond condensed matter physics because metamaterials and cold atomic gas systems can be engineered to exhibit dispersion relations of topological semimetals to large precision 56,57 . Controllable interactions in the cold atom systems makes realization of superconductivity and magnetism per our predictions possible in such systems. Beyond that, there are many related conceptual questions: how general is the connection between the surface states and the lowest pseudo-Landau levels in the presence of a gauge field. Can the picture be expanded to other semimetals and probably even to topological insulators, superconductors, weak, and fragile phases? We also note that Dirac nodal lines in the presence of effective gauge fields can lead to peculiar effective electrodynamics 58 , and it might be interesting to study the connection to surface states also in this context. 73 Note that this comes at considerably greater computational cost: previously, for each momentum (kx, ky) we had to diagonalize a 2Nz × 2Nz matrix to obtain the energies and states of the system, where Nz is the number of lattice sites in the z direction. Now, we instead have to deal with a matrix of size 2NzNx × 2NzNx for each ky, with Nx the number of sites in the x direction.

Appendix A: Strain-induced pseudo-magnetic field
In the tight-binding model the strain changes the x-bond length, which is implemented by a factor (1 − u 11 ) in front of cos k x (see Sec. H). Depending on the orbital structure also other modifications may appear, but these non-universal contributions are neglected in our consideration of the two-band NLSM model in this section (For completeness, we consider the effects of the non-universal orbital-mixing strain terms in Sec. L and demonstrate that they do not significantly influence our results.) Thus, the Hamiltonian in the presence of strain is assumed to bẽ We assume that the radius of the nodal loop Q is small and expand the Hamiltonian near an arbitrary nodal point Q = Q(cos θ, sin θ) = QQ, where θ is the polar angle in momentum space. In this way, we obtain where q = q(cos θ, sin θ) = qQ describes the deviation of the momentum from the nodal line. Notice that both here and in the following sections q takes both positive and negative values, i.e., it is not the absolute value of q. The Hamiltonian h(q) is parameterized by three independent parameters: θ, q, and q z . The strain term u 11 can be divided between q x and q y in such a way that the pseudo-magnetic field B 5 is continuous around the nodal loop. To see this, we notice that the strain-induced gauge potential is leading to the pseudo-magnetic field given by Eq.
(2) in the main text This can be used to define the pseudo-magnetic length l B = /eB 5 = √ RQ. The structure of the pseudo-magnetic field in k-space is illustrated in Fig. 4.

Appendix B: Bulk wave functions for the lowest pseudo Landau level
All pseudo-Landau levels of the Hamiltonian (A2) can be obtained analytically. We skip the explicit derivation here as it is analogous to the derivation of the Landau-level solutions for graphene. We concentrate on the zeroth Landau-level wave functions, which can be written as where ξ q = (z + z q )/l z , z q = ql 2 B and the localization length in the z-direction is given by Figure 4. Sketch of the momentum-dependent pseudo-magnetic field B5 along the nodal loop (red dotted line).
Here, the factor t2 tQ describes the elongation of the elliptical semiclassical cyclotron orbit in the z-direction. One can also perform a gauge transformation so that the wave function is localized in the x − y plane and is a plane wave in the z-direction. In that case, the localization length is Here, the factor tQ t2 describes the elongation of the elliptical cyclotron orbit within the x − y plane. This transformation is in analogy with the gauge transformation that can be performed on Landau-level wavefunctions in two dimensions.

Appendix C: Effective Hamiltonian for the connection of drumhead surface states to lowest pseudo-Landau level
To simplify our considerations, we look at solutions along a specific cut through momentum space. Without loss of generality, we choose Q y = 0 for a circular nodal line. Along this cut, after performing a basis rotation in the space of Pauli matrices of the form σ z → σ x and σ x → σ y , the effective nodal-line Hamiltonian under strain is By denoting Ψ q (x, z) = e iqx Φ(z) = e iqx (φ 1 , φ 2 ) T and assuming E = 0, we obtain two uncoupled equations We require that the solutions Φ(z) = (φ 1 , φ 2 ) T satisfy boundary conditions φ 1 (z = −L/2) = φ 2 (z = L/2) ≈ 0 up to corrections exponentially small in 1/L.

Drumhead surface states in the unstrained model
We start by inspecting the solutions of the unstrained limit, i.e., in the limit l z → ∞. In this case, the solutions are the drumhead surface states discussed in the main text, which in the new basis can be written as These solutions are states localized to opposite surfaces (top and bottom) and decay exponentially into the interior of the system.

Drumhead surface states and pseudo-Landau levels in the presence of strain
Once we switch on the strain, we obtain the following general solutions We immediately see that the drumhead-state solutions are recovered in the limit l z → ∞. At finite strain, the drumhead state from the top surface, corresponding to φ 1 , evolves into a PLL bulk state where Erf(ξ) is the Gauss error function. Interestingly, this solution satisfies the boundary condition φ 1 (z = −L/2) ≈ 0 also for q > 0. Hence, the radius of the plateau of zero-energy states grows with the system width.
On the other hand, the drumhead state from the bottom surface, corresponding to φ 2 , does not evolve into a PLL bulk state. It is modulated with a function e z 2 and therefore grows faster than exponentially for z → ±∞. For sufficiently small strain, we nevertheless expect to recover a state localized to the bottom surface that decays into the interior. By completing the square and normalizing the wave function, we obtain where we have used the imaginary error function Erfi(ξ) = −iErf(iξ). As before, this solution is subject to the boundary condition φ 2 (z = L/2) ≈ 0. This condition will be satisfied as long as the turning point of φ 2 , where its slope changes from negative to positive, is beyond the boundary of the system at z = L/2. This turning point is at −ql 2 B . Hence, φ 2 is only a valid solution for momenta deep inside the nodal circle, not even at k x = ±Q (q x = 0), which would contradict our numerical findings. This problem can be resolved by defining the approximate analytical solution asφ

Effective model in the presence of strain
Using the approximate solutions Φ 1 = (φ 1 , 0) and Φ 2 = (0, φ 2 ), we will now write down an effective Hamiltonian for the flat-band states. For that, we make use of our insights from numerics: the degenerate states are at zero energy in the interior of the flat band, whereas they split and acquire a finite dispersion at the flat-band boundary. This is due to the Landau-level state shifting towards the bottom surface where the other state is localized. The two states hybridize due to the surface effects, and we model this by introducing an effective coupling where the strength of the coupling δ can be determined by fitting the energies of the effective model E(q) = ±H 12 (q) to the numerical results obtained from the full tight-binding Hamiltonian. The coupling H 12 (q) can be expressed as We start by discussing the PLL ferromagnetism in the bulk. We point out that this theory describes several different types of magnetic order parameters because their projections to PLL bulk wave functions are the same. However, the surface effects distinguish some of these magnetic order parameters from each other, and in the end of the section we identify the magnetic order parameters favoured by the surface effects.
The Coulomb interactions projected to the zeroth PLL wave functions [Eq. (B1)] are described by the Hamiltonian where σ describes the spin of the electron, and Here the summations are subject to the restrictions: • Q, Q , Q + K and Q − K are on the nodal line, • q and q are perpendicular to the nodal line at the corresponding Q and Q , respectively, • q + k and q − k are perpendicular to the nodal line at the points described by Q + K and Q − K, respectively.
The direction of the magnetization can be chosen arbitrarily due to the SU (2)-symmetry, and here we have chosen it to be along the z-direction. This way we obtain We now additionally assume that the density is a constant independent of the position where 0 ≤ C ≤ 2 is the filling factor of the flat bands. Then the mean field Hamiltonian simplifies to a form where C † Q,q = (c † Q,q,↑ , c † Q,q,↓ ), Figure 5. Phase diagrams based on the mean-field equations for (a) magnetic order and for (b) superconducting order. V0 and G0 are the interaction strengths, mz is the magnetization, ∆ is the superconducting gap, µ is the chemical potential, and C is the density. The green lines correspond to the analytical formulas of the critical temperatures Tc. and The magnetization m z is independent of Q and q due to the spatial homogeneity, and thus m z and the chemical potential µ should be solved self-consistently using Eqs. (D6) and (D9). By straightforward calculation we obtain the mean-field equations where β = 1/k B T and we have defined V 0 = Q q V P (q, q , Q − Q , q − q ). We solve these equations numerically by reformulating them in terms of a minimization problem, for which we then compute the minima using a stochastic algorithm based on Basin-hopping. The results are presented in Fig. 5(a). Furthermore, from the equations above it is easy to see that the critical temperature for magnetism depends on the filling factor as At zero temperature we obtain The interaction strength depends on the parameters of the model as  where V C = e 2 /(4π 0 d) with the lattice constant d. Figure 6 shows the numerically calculated V 0 as a function of the curvature radius R. We have used model parameters estimated from a fit to the nodal line in CaAgP (see caption of Fig. 6), with the parameter t 1 tuned from 0.85t to 0.25t to obtain a smaller nodal circle in agreement with our initial assumptions of sufficiently small Q. We find that V 0 scales like 1/ √ R. By using R = 0.5 µm, we therefore estimate that V 0 = 114.5 meV / . We further estimate that typical values of the dielectric constant for these materials are on the order of 100. Hence, critical temperatures can be on the order of T c,m = 3 K.
Notice that we have assumed a larger dielectric constant than typically observed in bulk semimetals 60 . The reason is that the enhanced density of states in our system is expected to lead to larger screening effects. We also point out that the relevant length scales in our problem are similar to those in twisted bilayer graphene and our estimate for the critical temperature agrees with the experimentally observed critical temperatures in that system. Thus, our estimate can be considered conservative guided by the current knowledge about its 2D analogue, twisted bilayer graphene, but it might also be possible to observe larger critical temperatures in nodal-line semimetals due to larger stability to fluctuations in 3D and larger variability of parameters with strain.
We point out that our calculation is compatible with various types of spatially uniform magnetic orders because they can lead to the same projected order parameter within the low-energy theory. For this purpose we now consider the order parameters of the form M = m j σ i s j , where Pauli matrices σ i (i = 0, x, y, z) and s j (j = x, y, z) correspond to the orbital and spin degrees of freedom, respectively. Due to the SU (2) spin symmetry of the nonmagnetic phase, we can restrict our considerations to the case j = z.
The zeroth PLL wavefunctions are eigenstates of σ y , and therefore for order parameters with σ x s z and σ z s z the bulk PLL states stay at zero energy such that the system remains gapless. Thus, these order parameters are not energetically favored. On the other hand, order parameters σ 0 s z and σ y s z lead to the same projected order parameter within the PLL states and open a gap in the bulk. Thus, both of these order parameters are compatible with the calculation given above and good candidates for the ground state. However, there is a further distinction between these order parameters when the surface effects are taken into account as shown in Fig. 7. The order parameter σ y s z gives rise to a full gap both in the bulk and at the surface, but the order parameter σ 0 s z just shifts the two spin blocks of the Hamiltonian oppositely in energy so that the two shifted sets of bands cross at the edge of the flat bands. Therefore, we expect that the order parameter σ y s z will be energetically favored.
To shed light on the structure of this order parameter, we perform a change of basis throughM = U † M U with In this basis the order parameterM = m zσz s z is diagonal, and the magnetic order is staggered with respect to PLL bulk states and drumhead states at the bottom surface. In particular, if we restrict ourselves to the subspace of PLL bulk states, this corresponds to a ferromagnetic order parameter. Therefore, the calculations for the magnetic phase of the PLL bulk states, presented at the beginning of this section, is fully compatible also with respect to this order parameter.
Appendix E: Mean field theory for flat-band superconductivity When considering an attractive interaction between electrons within the reduced BCS Hamiltonian approach, the pairing Hamiltonian takes the form where the projected two-particle interaction potential, corresponding to an effective on-site attraction V (r 1 , r 2 ) = −gδ(r 2 − r 1 ), is Here g > 0 corresponds to the attractive interaction. Therefore, the gap equation can be written as This equation should be solved self-consistently in the presence of the density constraint Assuming a homogeneous order parameter ∆(Q, q) = ∆, we obtain and where we have defined From these equations one obtains a general solution for µ (valid at all temperatures) given by and then ∆ 0 can be solved from Eq. (E5) as a function of temperature. At zero temperature we get and the critical temperature is given by Moreover, we have also solved the mean-field equations numerically [see Fig. 5 Similarly as in the case of magnetism, the largest order parameter ∆(T = 0) = G 0 /2 and critical temperature k B T c,sc = G 0 /4 for superconductivity are obtained at half filling, but the important qualitative difference is that when the system is doped away from half filling the critical temperature for superconductivity decreases more slowly than the critical temperature for magnetism. For typical model parameters (see caption of Fig. 6), R = 0.5 µm, and g = 1.0 eV nm 3 , the effective interaction strength is G 0 = 0.3 meV indicating that the critical temperature can be on the order of T c,sc = 1 K.
Appendix F: Out-of-plane superfluid stiffness In this section, we calculate the superfluid stiffness in the z-direction by studying the energy cost of creating a phase-gradient ∆(z) = ∆ 0 exp(ikz) similarly as in Refs. 6 and 42. In our formalism, we can conveniently study this by assuming that For simplicity, we assume T → 0, C → 1 and k → 0, so that f 0 = 1/2 and the energy is where E c = − 1 2 Q,q ∆ = −V n 0 ∆/2 is the condensation energy, V is the volume, D s,z is the superfluid weight in the z-direction, and n 0 = Q/(2πl 2 B ) is the density of particles within the flat band. In this way, we identify The in-plane superfluid stiffness D s also relates the supercurrent j in a superconductor to a gauge potential A (in Coulomb gauge) by Here we have introduced a prefactor 4e 2 / 2 so that we have a common convention with Sec. F, where the stiffness was determined from the energy cost. Although the stiffness is generally a tensor, in the case studied here the off-diagonal components vanish, so that we are only interested in the stiffness in different directions.
In general the stiffness is composed of a conventional and a geometrical contribution but due to the dispersionless flat band the conventional contribution vanishes The geometric contribution at half-filling C = 1 can be calculated as 43-48 where g ij (k) is the Fubini-Study metric of the spin-up Bloch wave functions u(k). It is related to the quantum geometric tensor through the relation The Berry curvature of the system is also related to this quantity by From the last two equations we see that the second term in g ij is zero as it contains only symmetric integrals over odd functions in z and z . Hence, the components of the Fubini-Study metric g ij simplify to To evaluate the k-space integrals of g ij = g ij (α) we consider a system extending from z = −L/2 to +L/2. The PLL states at k = Q are centered at z = 0 and reach the top (bottom) surface at k top = Q − L/2l 2 B (k bot = Q + L/2l 2 B ). This defines our domain of integration and we obtain The expressions for the out-of-plane [Eq. (F2)] and in-plane [Eq. (G24)] stiffness are related to each other by replacement of l 2 z with l 2 xy /2 which arises because of the different elongations of the semiclassical cyclotron orbits in the z-direction and within the (x, y)-plane and due to the angular average within the (x, y)-plane. Similar results for the superfluid stiffness are expected also beyond mean-field approximation 61 .
We have checked our results for D s,geom also numerically (see Fig. 8). In the numerical calculation we adopt the essence of a method for calculating the Berry curvature in a discretized Brillouin zone 62 to efficiently compute the quantum geometric tensor G ij and utilize Eq. (G6). We find that the numerical and analytical results agree well up to the addition of a constant independent of L. We attribute this constant to the drumhead surface states coexisting with the zeroth PLL bulk states.
where r 0 is the hopping vector between the two orbitals involved in equilibrium, δr is the deviation from the equilibrium value r 0 , and t(r) is the hopping function. In general, there can be other correction terms due to the change of the angle between unlike orbitals, but we assume here that these contributions are sufficiently small. By writing δr = U r 0 , with the strain tensor U = (u ij ), as well as r 0 = r 0 e 0 and t(r 0 ) = t 0 , the equation above becomes, t(r 0 + δr) = t 0 + (e 0 · U e 0 ) r 0 ∂t ∂r r0 . (H2) For simplicity, we further assume that the hopping function can be approximated by t(r) = t 0 r 0 /r, such that we end up with the following expression: In the numerical strain implementation for CaAgP and rhombohedral graphite, all hopping terms were modified according to this formula, assuming that the dominant strain contribution comes from the u 11 component for a cylindrical substrate bent in the x direction.
Appendix I: Tight-binding model for rhombohedral graphite Our calculations for rhombohedral graphite in the main text are based on the following model 64 , Here, b is the distance between adjacent layers, δ i are the intra-layer nearest-neighbor hopping vectors, and n i are the intra-layer next-nearest-neighbor hopping vectors 65 . We take the model parameters from literature 66 : γ 0 = 2.58 eV, γ 1 = 0.34 eV, γ 2 = 0, γ 3 = 0.17 eV, and γ 4 = 0.04 eV. For the implementation of the cylindrical strain profile, we assume that the bending direction is along the x direction. This modifies the hopping amplitudes of our tight-binding model according to Eq. (H3). Hence, and We start from the Hamiltonian defined in Eq. (I1). To proceed analytically, we set γ 2 = γ 3 = γ 4 = 0, such that Θ(k) = 0. By using in Φ(k) the intra-layer nearest-neighbor hopping vectors and after implementing the strain terms according to Eq. (I4), we obtain with u 11 = z/R.
We now expand the strain Hamiltonian in Eq. (J4) around the K (K ) point with k = K ( ) + Q. This leads to Note that we have only expanded in the k x and k y directions, since the nodal line extends over the whole BZ in the k z direction. Before we proceed, we first look into the structure of these nodal lines.

Nodal lines of the unstrained system
The nodal lines of the unstrained system (u 11 = 0) spiral around axes parallel to the k z axis that go through the K and K points of the hexagonal BZ, Close to these axes, the energies of the Hamiltonian (with u 11 = 0) at momenta k = K ( ) + Q are to first order in Q x and Q y which leads to the following equation, Hence, the zero-energy states lie on a spiral given by where we have set q y ≡ q for Q 0,y > 0, and q y ≡ −q for Q 0,y < 0. With this definition, momenta with q < 0 (q > 0) are inside (outside) the projected nodal circle. For sufficiently small nodal lines, we can assume that γ 1 < γ 0 . In this case, the sign in front of γ 0 determines the overall sign of the terms γ 1 ± 3γ 0 . We therefore define, σ ± = 4Rbγ 1 (3γ 0 ± γ 1 ) (J19) and λ = 3aγ 0 2bγ 1 q. (J20) With this, we obtain the following zero-energy solutions: • for Q 0,y > 0, ψ 1 (z) = A 1 e −z 2 /2σ 2 + e λz ∝ e −(z−λσ 2 + ) ψ 2 (z) = A 2 e +z 2 /2σ 2 + e −λz ∝ e +(z−λσ 2 + ) • for Q 0,y < 0, ψ 1 (z) = B 1 e +z 2 /2σ 2 − e λz ∝ e +(z+λσ 2 − ) Here, the solutions describe a PLL bulk state and a top-surface state for Q 0,y > 0, but a PLL bulk state and a bottom-surface state for Q 0,y < 0. Furthermore, the PLL bulk-state shifts along this direction are not symmetric with respect to the center of the flat band: starting from Q 0,y < 0, the bulk-state shifts towards the top surface as we move towards the center of the flat band, whereas it shifts towards the bottom surface, if we do the same starting from Q 0,y > 0. This is in qualitative agreement with the numerical results [see Fig. 9(b)] and we have also confirmed that analytical and numerical solutions agree quantitatively.
3. Analytical solutions in the kx direction with respect to K Let us now turn to the perpendicular subset in momentum space for which Q y = 0. Before we proceed, recall from Fig. 9(b) that the flat band does not grow in the Q x direction. In particular, we know from numerics that the two PLL bulk states hybridize at the nodal-line momenta Q x,0 . Therefore, an expansion around a nodal point as before is not a promising approach to find zero-energy solutions.
Recall that we have not yet expanded the Hamiltonian along k z . From numerics, we obtain that the Fouriertransformed PLL bulk solutions are centered at k z values that depend on q x , suggesting that the corresponding expansion point along the k z direction should be a function of q x . In particular, we find the following relation between the bulk-state centers K z and q x q x (K z ) = Q cos (bK z + π/6),   (GGA) 69 has been applied. We have used a 12 × 12 × 18 k-point grid centered at Γ. After computing the Bloch wave functions in DFT, we construct corresponding Wannier functions (WFs) 70,71 using the WANNIER90 code 72 . To extract the orbital character of the electronic bands at low energies, we use the Slater-Koster interpolation scheme based on the WFs. Furthermore, we neglect spin-orbit interactions, which are small in this material.
As a first step, we construct a 12-orbital model based on the 3p orbitals of the three P atoms and on the 5s orbitals of the three Ag atoms in the unit cell. Our band structure results are shown in Figs. 10(a) and (b) and are in agreement with the literature 51 . The match between the DFT band structure and the interpolated band structure obtained from the WFs is good around the Fermi level. To obtain a simpler model catching the essential physics only at low energies, we next construct an effective two-orbital tight-binding model based on one p z orbital centered at one P atom and on one s-orbital centered at the middle of a triangle of Ag atoms. The band interpolation of the two-orbital model is shown in Figs. 10(c) and (d). As we can see from the comparison between Figs. 10(b) and (d), we do not lose accuracy between -0.60 and 0.85 eV moving from the 12-orbital to the two-orbital model. This model, as obtained from the Wannier interpolation, still has a large number of parameters. To reduce this number to a managable value, we finally set an energy cut-off of 6 meV for the tight-binding parameters such that the dispersion close to the Fermi level is still captured correctly. This results in a two-band model with 37 different parameters, which is used for the strain implementation in the main text. Figure 11 shows the energy bands of the system in such a geometry, i.e., the system is finite in x and z with widths of L x and L z , respectively, in units of the lattice constant of the unstrained system. For comparison, in Figs. 11(a) and (b) we also plot the energies of the unstrained system and the strained system without the s−p terms. We observe a number of flat subbands close to E = 0 comprising the top-and bottom-surface drumhead states in the case of the unstrained system, and surface and bulk pseudo-Landau level (PLL) states in the case of the strained system. Most importantly, we find that the k y extent and the number of flat subbands are larger than in the unstrained system. This reflects the growing density of states as strain is applied.
Finally, in Fig. 11(c) we plot the energy subbands of the strained system including the s − p strain terms. We find that the spectrum is strikingly similar to Fig. 11(b). In particular, number and extent of the flat subbands close to E = 0 are larger than in the unstrained system and almost identical to Fig. 11(b). The only difference is the subband level spacing away from E = 0. Moreover, we have checked that the flat subbands indeed consist of surface states and bulk PLL states behaving in a way identical to what we observe in the absence of the additional s − p strain terms: the centers of the bulk PLL shift as we move along k y or along the subband index.
This confirms that orbital-mixing strain terms have only a minor effect on the details of the spectrum and the system still supports 3D flat bands.