Ab initio framework for systems with helical symmetry: theory, numerical implementation and applications to torsional deformations in nanostructures

We formulate and implement Helical DFT -- a self-consistent first principles simulation method for nanostructures with helical symmetries. Such materials are well represented in all of nanotechnology, chemistry and biology, and are expected to be associated with unprecedented material properties. We rigorously demonstrate the existence and completeness of special solutions to the single electron problem for helical nanostructures, called helical Bloch waves. We describe how the Kohn-Sham Density Functional Theory equations for a helical nanostructure can be reduced to a fundamental domain with the aid of these solutions. A key component in our mathematical treatment is the definition and use of a helical Bloch-Floquet transform to perform a block-diagonalization of the Hamiltonian in the sense of direct integrals. We develop a symmetry-adapted finite-difference strategy in helical coordinates to discretize the governing equations, and obtain a working realization of the proposed approach. We verify the accuracy and convergence properties of our numerical implementation through examples. Finally, we employ Helical DFT to study the properties of zigzag and chiral single wall black phosphorus (i.e., phosphorene) nanotubes. We use our simulations to evaluate the torsional stiffness of a zigzag nanotube ab initio. Additionally, we observe an insulator-to-metal-like transition in the electronic properties of this nanotube as it is subjected to twisting. We also find that a similar transition can be effected in chiral phosphorene nanotubes by means of axial strains. Notably, self-consistent ab initio simulations of this nature are unprecedented and well outside the scope of any other systematic first principles method in existence. We end with a discussion on various future avenues and applications.


Introduction
The discovery and characterization of novel nanomaterials and nanostructures constitutes one of the principal areas of scientific research today [1,2]. Such materials and structures hold the promise of unlocking remarkable and unprecedented material properties that are otherwise unavailable in the bulk phase (i.e., crystalline materials). In recent years, the discovery of novel nanostructures has garnered much attention and acclaim [3,4], and the unusual properties of these new materials have led to ground breaking applications in almost every branch of science and engineering [5,6].
Nanostructures appear in various morphologies (including fullerenes, nanotubes, nanoclusters and two-dimensional materials), and are usually associated with non-periodic symmetries. 1 The mathematical framework for classifying nanostructures [7,9,10] shows that a vast class of these materials can be described as being helical, i.e., their spatial atomic arrangement possesses helical symmetries. Helical structures include important technological materials such as nanotubes (of any chirality), nanoribbons, nanowires and nanosprings; miscellaneous chiral structures encountered in chemistry; and examples from biology, including tail sheaths of viruses and many common proteins [7,11]. Figure 1 shows instances of helical structures that have been actively investigated in the literature.
Helical structures have been conjectured to be a fertile source of novel materials with unusual and attractive properties [7]. This is due to the fact that atoms in such structures find themselves in locally similar environments [7]. Coupled with the quasione-dimensional nature of these systems, as well as the presence of symmetries in the underlying governing equations, this makes it likely that collective or correlated electronic effects (such as those leading to ferromagnetism, ferroelectricity and superconductivity) can emerge in these materials [15]. On the other hand, helical structures are also inherently chiral and can therefore serve as natural examples of materials systems in which certain forms of symmetry breaking in the governing equations can lead to unconventional transport phenomena [16][17][18][19].
Given the relative abundance of helical nanostructures in existing materials, their likelihood of being associated with hitherto undiscovered forms of matter displaying exotic materials properties, and their overall scientific and technological importance, there appears to be a pressing need for reliable and efficient computational tools for studying such systems. The broad goal of the present contribution is to take important foundational steps in addressing the above scientific issue. Specifically, we present here the mathematical formulation and numerical implementation of a novel computational method called Helical DFT, that can simulate helical structures ab initio. We also obtain a practical working realization of this (density functional theory based) self-consistent first principles technique, and illustrate some of its capabilities through the study of an emergent nanotube material with interesting properties.
To put our work into perspective, we remark that the use of first principles (i.e., quantum mechanical) techniques to design and study materials is a very active area   [12][13][14].
of scientific endeavor today, and it forms the bulk of computational materials science research [20][21][22][23][24]. Among the wide array of first principles methods available, Kohn-Sham Density Functional Theory (KS-DFT) [25] enjoys widespread usage since it offers a good balance between computational cost and physical accuracy as compared to other techniques [26]. The pseudopotential plane-wave method, also called Plane-wave DFT, is the most widely used implementation of Kohn-Sham theory [27][28][29][30], and it involves expanding the unknowns into linear combinations of plane-waves. Since plane-waves are naturally associated with periodic symmetries (they are in fact eigenfunctions of translational symmetry operators), Plane-wave DFT is ideally suited for studying bulk (i.e. periodic or crystalline) systems, and is often found to be fundamentally inadequate for studying systems with non-periodic symmetries. In particular, using a Plane-wave DFT code for studying a helical structure such as a chiral nanotube can require the use of large periodic unit cells often containing many hundreds (or even thousands) of atoms. 2 In contrast, a computational method which is faithful to the underlying helical symmetry of such a structure would require a small helical unit cell, containing far fewer number of atoms. Since ground state electronic structure calculations using density functional theory (DFT) scale as the cube of the number of atoms in the unit cell, while excited state calculations scale as the fourth power, the difference in simulation run times for such calculations, in these two scenarios (i.e., the correct use of helical symmetry vs. incorrect use of periodic symmetry) can be drastically different in practice.
The above considerations form our point of departure from a conventional formulation and implementation of KS-DFT, to one that is adapted for helical systems. In order to formulate the equations of KS-DFT for a helical unit cell, an appropriate version of the Bloch Theorem [38,39] is required. We establish this result rigorously in this work, and use it to set up an electronic band theory for helical structures. Subsequently, we develop the notion of helical Bloch states, and use their properties to derive of the equations of KS-DFT, as they apply to helical systems. A key component in our mathematical treatment is the definition and use of a helical Bloch-Floquet transform to perform a block-diagonalization of the Hamiltonian in the sense of direct integrals. Our use of rigorous mathematical arguments and appropriate mathematical tools 3 is one of the highlights of our framework, and it allows the governing equations to be obtained systematically, and without recourse to an excessive amount of intuition. 4 As far as we are aware, our work is the first in presenting such a derivation, and also in expressing the detailed form of the equations of Kohn-Sham theory for helical structures. The final form of the equations are such that they are readily suited for implementation within systematically convergent electronic structure methods such as those based on finite differences [31][32][33][34], finite elements [35][36][37] or spectral basis functions [8,45,46]. We choose a symmetry adapted finite difference method in helical coordinates for discretizing the governing equations in this work, and set up a computational framework for numerically solving the discretized equations in a self-consistent manner. This gives us a working realization of an ab initio computational tool -called Helical DFT -that can be used to perform predictive simulations of helical systems in a systematic and efficient manner. It can therefore aid in the discovery, synthesis and characterization of helical structures. Subsequently, the remainder of this work focuses on illustrating various numerical and application oriented aspects of this novel computational tool through examples based on nanotube systems. To the best of our knowledge, Helical DFT is the first computational 2 In contrast to plane-waves, the use of real space techniques based on finite differences [31][32][33][34] or finite elements [35][36][37] allow for non-periodic boundary conditions to be imposed in a straight-forward manner. However there does not appear to be any prior work on using these techniques for self-consistent first principles calculations of helical systems. 3 Due to the infinite nature of helical groups, the mathematical arguments presented here are of somewhat different and more subtle nature as compared to the ones that can be employed for cyclic groups [40]. However, they can be seen as being broadly connected in the sense that they both deal with Fourier analysis of the respective symmetry groups [41,42]. 4 Since a rigorous thermodynamic limit theory for the Kohn-Sham problem is unknown [43,44], a derivation of the equations of the theory, as it applies to condensed matter systems often makes use of physical intuition. This process is prone to conceptual errors however, and we are aware of literature that lists certain terms of the equations incorrectly. In any case, the final form of the equations appear to be well known in the electronic structure community at large, since DFT codes routinely make use of them for simulating the crystalline phase.
method for helical systems that is based on first principles, and one that also behaves systematically with respect to convergence properties. This, among other reasons, is made possible by our use of the aforementioned helical coordinate system. To the best of our knowledge, this has not been employed in electronic structure calculations before. While the study of helical structures has much scientific and technological merit in of itself, the development of a computational method for studying such systems also brings with it the added benefit of being able to simulate the behavior of nanomaterials under torsional deformations. As explained in [7,40], homogeneous deformation modes are commensurate with periodic symmetries (i.e., applying a homogeneous deformation to a periodic structure results in another periodic structure), while certain inhomogeneous deformation modes can be associated with non-periodic symmetries. An attempt to study such inhomogeneous deformations while using a periodic method is likely to involve various uncontrolled approximations, complications and computational inefficiencies [47,48]. This issue appears to have been recognized for some time in the nanomechanics and materials literature, leading to a considerable body of work centered around suggestions presented in [7], whereby pure bending deformations in atomistic systems are simulated using cyclic symmetries, while helical symmetries are used to simulate torsion [49][50][51][52][53][54][55][56][57][58]. A persistent issue with the simulations in these studies however, is that they have all been carried out using interatomic potentials or tight binding methods. Due to the well known deficiencies of these techniques in simulating real materials [48,[59][60][61][62], true first principles simulation methods that behave systematically, and also take into account cyclic and/or helical symmetries have been deemed highly desirable [7,8,58]. There has been recent progress on this very issue with regard to cyclic symmetries [40,63], and the resulting computational methods have been used to study the bending behavior of nanoribbons and sheets of two dimensional materials ab initio. In this sense, the current contribution follows up on this line of work by making a first principles simulation framework for torsional deformations available. Consequently, through the use of this framework, we are able to extract the behavior of nanotubes of black phosphorus (i.e., phosphorene nanotubes) and study their mechanical and electronic response as they are subjected to twisting. 5 The coupling of these responses leads to some interesting electronic transitions in this material that is likely to make it an attractive candidate for sensing, modulation and actuation applications.
The rest of this work is organized as follows. Section 2 establishes the mathematical framework for a systematic formulation of the governing equations, and also derives the relevant expressions explicitly. Section 3 discusses formulation of a numerical scheme based on finite differences in helical coordinates, and Section 4 presents simulation studies. Section 5 summarizes the work and suggests avenues for future research. The appendices 5 Exploitation of helical symmetries in ab initio calculations has also been considered in the chemistry literature in the context of Linear Combination of Atomic Orbitals (LCAO) methods [64][65][66][67][68][69]. However, these methods differ in their perspective from the current contribution in that they concentrate on using symmetry-adapted basis functions for reducing the computational cost of the multi-center integrals and the Hamiltonian matrix elements, whereas our focus is on the formulation of symmetry-adapted cell problems (in helical coordinates), and a systematically convergent numerical treatment of these cell problems. Due to basis incompleteness and superposition errors, it is often non-trivial to systematically improve the quality of the numerical solutions obtained via LCAO methods, in contrast to the techniques presented here. Finally, the connection of helical symmetries with torsional deformations, as well as the effect of such deformations on other material properties does not appear to have been considered in the chemistry literature.
contain additional information and discussions on mathematical tools and results that allow for this work to be self-contained.

Formulation
In this section, we describe the key aspects of Helical DFT. We begin with a formal discussion of helical groups and helical structures in Section 2.1, and then discuss Kohn-Sham DFT, as it applies to such systems in Section 2.2. The atomic unit system with m e = 1, e = 1, = 1, 1 4π 0 = 1, is chosen for the rest of the work, unless otherwise mentioned.

Helical symmetry groups, fundamental domains and helical structures
A helical structure (i.e. a structure with helical symmetries) can be defined through the action of a helical group on a set of non-degenerate points in space. This definition makes it necessary for us to make the notion of a helical group precise. Following standard practice in the literature [7-10, 49, 70, 71], we introduce helical groups as subgroups of the Euclidean group in three dimensions. This requires us to introduce some relevant notation and basic rules regarding operations with isometries, as we now do.

Helical symmetry groups
Let e 1 , e 2 , e 3 denote the standard orthonormal basis 6 of R 3 and let (x 1 , x 2 , x 3 ) denote the Cartesian coordinates of a generic point x ∈ R 3 . An isometry (or rigid body motion) in R 3 will be denoted using the notation Υ = (R|c), with R ∈ SO(3) denoting the rotation part of the rigid body motion, and c ∈ R 3 denoting the translation part. The action 7 of Υ : R 3 → R 3 on a point x ∈ R 3 is written as Υ • x = Rx + c. Given a collection of points S ⊂ R 3 , we will use the notation Υ • S to denote the action of the isometry on each of the points in S, i.e., There is a natural multiplicative operation associated with isometries (denoted as • here) that arises as a composition of their maps. Specifically, given isometries Υ 1 = (R 1 |c 1 ), Υ 2 = (R 2 |c 2 ), we may define a third isometry . It follows that Υ 3 = (R 1 R 2 |R 1 c 2 + c 1 ), and that in general the operation • is not commutative (due to non-commutativity of finite rotations about arbitrary axes). The • operation also allows the definition of whole number powers of Υ, 6 We will use the following notation in what follows: f (·) will be used to denote a function when we do not wish to highlight the dependence of the function on its arguments. · will be used to denote the norm of a function or vector and ·, · will be used to denote the inner product. Often, we will attach a subscript to these symbols to denote the specific space in which the norm or inner product is being considered. Vectors and matrices in R 2 or R 3 will be denoted in boldface, with lower case letters reserved for vectors and uppercase letters used for matrices. We will sometimes use the · symbol between vectors in R 2 or R 3 , to denote the inner product. If a function has dependence on multiple arguments, we may choose to separate the arguments using ';' to emphasize a parametrized dependence of the function on the arguments following ';'. 7 As the name suggests, isometries preserve distances (and hence, also angles), i.e.,∀ x, y ∈ R 3 , and a generic isometry Υ, it holds that Υ(x) − Υ(y) i.e., for n = 1, 2, . . ., we may define Υ n := Υ • Υ • Υ . . . (n times). It is then easy to check that Υ n admits the expression Υ n = R n ( n−1 j=0 R j ) c , where the notation R 0 is used to denote the identity matrix. The identity isometry leaves every x ∈ R 3 invariant and can be written as (I|0), with I denoting the identity matrix and 0 denoting the null vector in R 3 . Given the isometry Υ = (R|c), we can form the isometry Υ = (R −1 | − R −1 c), which satisfies Υ • Υ = Υ • Υ = (I|0). Hence, we will denote Υ as Υ −1 -i.e., the inverse isometry to Υ. The set of all isometries so defined, i.e., E = {Υ = (R|c) : R ∈ SO(3), c ∈ R 3 }, together with the operation • and the inverse element defined above, form a group [72]. 8 Let α and τ be real numbers 9 such that 0 ≤ α < 1 and τ > 0, and let R 2πα denote a rotation around axis e 3 by angle 2πα. Then, the rigid body motion Υ h = (R 2πα |τ e 3 ) will be called a helical isometry 10 about axis e 3 . The action of Υ h on a point x ∈ R 3 is to rotate it by angle 2πα about axis e 3 , while also translating it by τ along the same axis. 11 Furthermore, applying the formulae for the powers of isometries and their inverses shown above, we see that for m = 1, 2, . . ., . Combining these, we may define Υ m h for any m ∈ Z as Υ m h = (R 2πmα |mτ e 3 ), with the m = 0 case automatically resulting in the identity isometry (I|0). We may therefore state: Proposition 2.1 (Helical group generated by a single element). The set of isometries forms a discrete group under the operation •.
Additionally, let N ∈ N, let Θ = 2π N and for n = 0, 1, . . . , N − 1, let R nΘ denote a rotation around axis e 3 by angle nΘ. Then the set of isometries endowed with the operation • forms a cyclic group [40] of order N. Note that since the rotational parts of the isometries in group G 1 and C all share e 3 as the common axis of rotation, the elements of G 1 and C commute (i.e., for any Υ m h ∈ G 1 and Υ n . We may now consider the direct product of the groups G 1 and C defined above to obtain a new helical 8 Since only pure rotations are included, this is the so called Euclidean group of direct isometries in three dimensions [72]. The full Euclidean group also includes improper rotations. 9 Most of the discussion in this work naturally also extends to the case when −1 < α < 0. However, we will not be considering that case here. 10 Alternately referred to as a screw transformation in the crystallography literature [9]. 11 A simple way to see this is to resolve x along and perpendicular to e 3 , i.e., 7 group 12,13 : Proposition 2.2 (Helical group generated by two elements). The set of isometries forms a discrete group under the operation •.
Since G 1 and C are generated by single elements, they are Abelian groups. Furthermore, since G 2 is generated by two elements (i.e., the generators of G 1 and C) which commute among themselves, it is an Abelian group as well.
The action of the groups G 1 and G 2 on points in space are easily described using cylindrical coordinates: if x ∈ R 3 is point with cylindrical coordinates (r, ϑ, z), then the action of the group element Υ m h ∈ G 1 is to send it to a point with cylindrical coordinates (r, ϑ + 2πmα, z + mτ ), while the action of Υ m h • Υ n c ∈ G 2 is to send it to the point with cylindrical coordinates (r, ϑ + 2πmα + nΘ, z + mτ ). In what follows, we will use the notation Υ to denote a generic isometry from G 1 or G 2 .

Fundamental domains
Given a point x ∈ R 3 , and a group of isometries G (which could be the helical groups G 1 or G 2 described above, for instance), the orbit of x under the group is the set Given a collection of points S ⊂ R 3 and a group of isometries G, we will use the notation G • S to denote the orbits of each of the points in S under the group: Let O ⊂ R 3 be a domain with regular boundary that is invariant under a given helical 12 While the discussion presented here already makes it evident that both G 1 and G 2 are groups (and are in fact Abelian groups), see [10] for a more complete derivation of these groups, as well as other types of helical groups not considered in this work. 13 Note that the groups G 1 and G 2 contain a group of translations as a normal subgroup if α is a rational number. In certain terminology [9,10], such cases would be identified as rod groups and the term helical group would be reserved only for cases for which α is an irrational number (i.e., when the group is not equivalent to a periodic group generated by a single translation.) However, we will not make this distinction here. 14 With these hypotheses, the boundary of O, denoted ∂O, can be shown to be invariant under the group as well. 15 In practice, we will require the fundamental domain to have some regularity properties in addition to the conditions in eq. 7 and 8, e.g. it should be connected and have compact closure.
8 and for Υ 1 , Υ 2 ∈ G: To see concrete examples of the sets O and D, let D R denote an open disk of radius R on the e 1 , e 2 plane, i.e., and let C denote the infinite cylinder obtained by translating D R along e 3 , i.e.: Then, the cylinder C has all the properties required of the domain O. Furthermore, we observe that the finite cylinder D G 1 = D R × {x 3 e 3 : 0 ≤ x 3 < τ } serves as the fundamental domain of G 1 in C. Finally, the sector with slanted walls, described in cylindrical coordinates as: serves as the fundamental domain of G 2 in C.

Helical Structures
A helical structure i.e., an atomic/molecular structure with helical symmetries is simply the orbit of a set of non-degenerate points under the action of one of the helical groups G 1 or G 2 . More precisely, let P G 1 ⊂ D G 1 or P G 2 ⊂ D G 2 in case of G 2 be a finite collection of distinct points labeled x k These points are representative of atomic positions within the fundamental domain and we will refer to them as simulated points or simulated atoms. The (valence) nuclear charges corresponding to these atoms will be denoted as Z k in case of G 2 . A helical structure is simply a set of the form: Additionally, for any Υ ∈ G 1 , the atom at the location Υ • x k is taken to be of the same species as the atom at x k ∈ P G 1 (similarly also for Υ ∈ G 2 and x k ∈ P G 2 ), and so it is associated with the same (valence) nuclear charge Z k .

Kohn-Sham Problem for Helical Structures
The Kohn-Sham equations, as they apply to finite structures can be found in numerous references [26,33,63]. In order to formulate an appropriate version of the Kohn-Sham equations for helical structures however, we need to keep in mind a few typical features of such a structure. In what follows, for the sake of simplicity, we will consider in detail the case of a structure associated with a helical group generated by a single element (i.e., the group G 1 described above). We will comment on modifications to the above case that need to be considered while dealing with a structure associated with a helical group generated two elements (i.e., the group G 2 described above), and present the final expressions/equations for this case in Appendix C. A more detailed discussion of the modifications and the application of resulting equations is the scope of ongoing and future work [73].
Helical structures are essentially quasi-one-dimensional in nature. This implies that they have limited spatial extent in the e 1 , e 2 plane, while being infinitely extended along the e 3 direction. Consequently, it is appropriate to set up the Kohn-Sham equations for such a structure in a computational domain which is of limited spatial extent in the e 1 , e 2 plane, while being infinite in extent along e 3 . This, along with the requirement that a symmetry adapted formulation of the Kohn-Sham equations needs to be solved on a domain that is also invariant with respect to the symmetry operations of the helical structure, suggests the cylinder C as being a natural choice for the computational domain (for a helical structure generated by a single element). The radius of this cylinder has to be consistent with the requirements that the all the atoms of the helical structure should be located sufficiently away from the lateral surface of the cylinder so as to allow sufficient decay of various fields that appear in the Kohn-Sham problem.
The quasi-one-dimensional nature of the systems under study results in additional complications. Specifically, due to the infinite extent of the system along the e 3 direction, the system is associated with an infinite number of electronic states 16 as well as an infinite number of nuclei. This potentially poses divergence issues while computing the electrostatics terms in the Kohn-Sham problem [63,74] and it is dealt with in this work by solving an appropriate symmetry adapted Poisson problem involving a neutral charge distribution -such a charge distribution arises as a combination of the electron density and the nuclear pseudocharges associated with the structure. Additionally, the infinitely many electronic states have to be incorporated into the Kohn-Sham problem in a manner that is consistent with the Pauli exclusion principle and the Aufbau principle [26,75]. Taking cue from the solid state/condensed matter physics literature -specifically, ab initio calculations of crystalline solids [75,76] -we address this issue here by formulating a band theory of electronic states for helical structures. This allows the Kohn-Sham problem for the entire helical structure, as posed on the cylinder C, to be reduced to computations on the fundamental domain when augmented with appropriate boundary conditions. 17 A key ingredient of the band theory for helical structures is an appropriate version of the Bloch theorem [38,76,77] for such systems. The form of this mathematical result can be guessed by looking at the analogous case of the Bloch Theorem for one dimensional 16 In general, these states would be expected to be delocalized over the entire volume of the cylinder C. 17 Note that we are not attempting to solve the thermodynamic limit problem associated with the helical structure in this work. Instead, we are postulating the form of the governing equations at the thermodynamic limit (i.e., for the infinite helical structure which is under study) and expressing them in a mathematically rigorous manner. This is necessary so that we can then numerically solve these equations and extract physical properties of systems of interest. In contrast, the thermodynamic limit problem would involve the passage from a finite (truncated) helical structure to the infinite one keeping various energetic contributions in mind, and is far beyond the scope of the current contribution. periodic systems 18 and the result appears to have been made use of by earlier researchers in various contexts [8,53,66,67,69,[78][79][80][81][82][83][84]. However, a rigorous mathematical derivation of the result does not seem to appear anywhere in the literature -other than in [8], where a proof of the existence of Helical Bloch waves was sketched by using tools from the theory of linear elliptic partial differential equations. In what follows, we address this gap in the literature and follow up on [8], by establishing the existence and completeness of helical Bloch waves, and then use this to gain insight into the spectrum of the single electron Hamiltonian associated with helical systems (i.e., to set up an electronic band theory for such systems). This information is subsequently used to set up the governing equations of the system. Our mathematical treatment closely follows the techniques presented in references [8,[85][86][87][88].

Analysis of the single electron problem for helical structures -helical Bloch waves
As a starting point, we consider the single electron Hamiltonian: with the real valued continuous potential V (x) invariant under the helical group G 1 , i.e., This operator naturally arises during each self-consistent field iteration cycle in Kohn-Sham calculations 19 and in that scenario, the invariance of the potential automatically follows from the invariance of the electron density [8].
We are interested in functions ψ that satisfy the equation Hψ = λψ within the region C in an appropriate manner. Additionally, to model the decay of the eigenstates as one moves away from the axis of the cylinder to infinity [89,90], we will enforce Dirichlet 18 See e.g. equations 25, 26 in [63]. 19 Within the setting of the local density approximation and the use of local pseudopotentials for example, V (x) can be identified as the total effective potential appearing in the Kohn-Sham equations and can be written as the sum of electrostatic and exchange-correlation terms, i.e., V (x) = V es (ρ(x)) + V xc (ρ(x)).
boundary conditions on the lateral surfaces of the cylinder 20 , i.e., ψ(x) = 0 for x ∈ ∂C. 21,22 Helical Bloch waves (or helical Bloch states) are solutions to the above equation which have the ansatz: Here φ(x; η) group invariant i.e., 20 This "wire" boundary condition is commonly employed in the literature for studying quasi-1D systems [34,91]. This boundary condition allows the operator H to have some convenient properties without having to enforce any specific decay conditions on V (x) as one moves away from the axis of the cylinder. 21 In what follows, we will use the following notation: if A is a measure space with measure µ, then for 1 ≤ p < ∞, we will use L p (A, B, µ) to denote Lebesgue measurable functions f : A → B which satisfy A f p B dµ < ∞, and we will use L ∞ (A, B, µ) to denote functions for which ess. sup. x∈A f (x) B < ∞. In particular, if A is a domain in R 3 , we will use L 2 (A) to denote the usual Hilbert space of complex valued functions on A which are square integrable (using the Lebsegue measure). The inner product of two functions on this space will be expressed as: Furthermore, H k (A) will denote the Sobolev space of tempered distributions whose k th weak derivative lies in L 2 (A), while H 1 0 (A) will denote the subspace of functions in H 1 (A) which vanish at the boundary of A in the trace sense. Finally, the rank one operator created as the tensor product of two functions f 1 , f 2 ∈ L 2 (A), i.e., f 1 ⊗ f 2 will act on a generic function f ∈ L 2 (A) to yield f, f 2 L 2 (A) f 1 . 22 We may view H as an unbounded operator on L 2 (C) with the function space Dom.(H) = H 2 (C)∩H 1 0 (C) as the domain of the operator. The operator H is formally symmetric (or, in linear algebra terminology, Hermitian since the underlying function spaces are complex): if f 1 , f 2 are Schwartz functions in C which obey the boundary condition f 1 (x) = f 2 (x) = 0 for x ∈ ∂C, we have: On using integration by parts [92] and the decay of f 1 and f 2 as x 3 → ∞, we get: Here ds denotes the oriented surface measure. The second term on the right-hand side above vanishes due to the boundary conditions obeyed by f 1 , f 2 on ∂C and so, this leaves us with: In a similar manner, we get: as the potential V (x) is real. Since Schwartz functions are dense in the domain of H, the result follows. The direct integral decomposition of H (Appendix B) makes it easy to appreciate that H is in fact self-adjoint. and obeys the boundary condition: commensurate with the boundary condition on ψ. The parameter η serves a role that is analogous to k-points in periodic calculations and as shown later, it can be chosen such that η ∈ [− 1 2 , 1 2 ). In what follows, we first show the existence of such solutions and then demonstrate their completeness. In essence, these results together give us information that certain special electronic states (i.e., helical Bloch states) can be always found to be associated with the single electron Hamiltonian of a helical structure, and they further inform us that such special states can be used to characterize all of the possible electronic states of the system (within the single electron model). Therefore, it is sufficient for us to restrict our attention to these states while discussing the spectrum of the single electron Hamiltonian associated with a helical structure. Our derivation of these results follows techniques employed in classic references on the mathematical theory of Bloch waves in crystals [85,87] and builds the theory in a "bottom up" manner using standard tools from functional analysis and the theory of linear elliptic operators (see [92][93][94] for relevant background material). In subsequent sections (Section 2.2.2, Appendix B), we use techniques presented in [88] to use helical Bloch waves for "block-diagonalizing" the single electron Hamiltonian through the apparatus of direct integrals, and then use this formalism to derive governing equations.
First, to demonstrate the existence of these special solutions, we have: Proof. We fix η ∈ R and substitute the helical Bloch wave ansatz in the equation Hψ = λψ to find that φ(x; η) should obey the following auxiliary equation in the region C: Additionally, φ should be group invariant and obey the zero Dirichlet boundary condition on ∂C. Let D denote the interior of the fundamental domain D G 1 , i.e., it is the open set described in cylindrical coordinates as D = {(r, ϑ, z) : 0 ≤ r < R, 0 < z < τ }. The boundary of D includes the lateral surface ∂D r=R = {(r, ϑ, z) : r = R, 0 < z < τ } that is shared with C, as well as the discs ∂D z=0 = {(r, ϑ, z) : 0 ≤ r ≤ R, z = 0} and ∂D z=τ = {(r, ϑ, z) : 0 ≤ r ≤ R, z = τ }, which are both parallel to the e 1 − e 2 plane. The group operation Υ h (i.e., the generator of the group G 1 ) maps ∂D z=0 to ∂D z=τ and conversely, Υ −1 h maps ∂D z=τ to ∂D z=0 . We now restrict the auxiliary eigenvalue problem h aux η φ = λφ as outlined in eq. 25, to the region D by imposing the boundary conditions φ( and (as before), φ(x, η) = 0 for x ∈ ∂D r=R . The operator h aux η on L 2 (D) is uniformly elliptic, and as shown in Appendix A, it is also symmetric with the above boundary conditions. Since D is a bounded domain and V (x) ∈ L ∞ (D), the operator h aux η satisfies the conditions of Gårding's Inequality (Theorems 9.17, 13 9.18 in [93]; Section 6.2 in [92]). This guarantees that h aux η has a unique self-adjoint extension in L 2 (D), which we also denote as h aux η here. Furthermore, as a consequence of the Rellich-Kondrachov Compactness Theorem (Theorem 7.29 in [93]; Section 5.7 in [92]), h aux η can be shown to have a compact resolvent (Lemma 9.20 in [93]). Consequently, h aux η has a discrete set of eigenvalues λ j (η) and corresponding eigenfunctions φ j (x; η) (Theorem 6.29 in [95]; Theorem 9.22 in [93]). Each eigenvalue is of finite multiplicity and such that λ j (η) → ∞ as j → ∞. Results from elliptic regularity theory (Sections 9.5, 9.6 in [93]; Section 6.3 in [92]) imply that φ j (·; η) ∈ H 2 (D). We now use the boundary conditions on φ outlined above to extend the eigenfunctions φ j (x; η) to all of C, noting that these boundary conditions are meaningful in the trace sense since the eigenfunctions are in H 2 (D). Thereafter, defining ψ j (x, η) = e −i2πη x 3 τ φ j (x, η), for j ∈ N and x ∈ C, establishes the theorem.
We define Λ = {λ j (η) : η ∈ R, j ∈ N} and Ψ = {ψ j (·; η) : η ∈ R, j ∈ N} as the collection of generalized eigenvalues and generalized eigenfunctions 23 associated with H. The first observation we make is that the sets Λ and Ψ are unchanged upon restricting η ∈ [− 1 2 , 1 2 ). To see this, we recall that λ j (η) and ψ j (·; η) are obtained by computing the spectrum of H when subjected to the conditions 24 in eqs. 22, 23. However, these equations can be equivalently recast as the following condition on ψ: or more generally, for m ∈ Z: In other words, solving Hψ = λψ while imposing the condition ψ(Υ h • x) = e −i2πη ψ(x) on ψ also gives us the sets Λ and Ψ. Since e −i2πη = e −i2π(η+n) for any n ∈ Z, we see that the boundary conditions on ψ do not change upon translating the value η by an integer. Thus, it suffices to restrict η ∈ [− 1 2 , 1 2 ). In what follows, we will denote I = [− 1 2 , 1 2 ), and we will re-define Λ = {λ j (η) : η ∈ I, j ∈ N} and Ψ = {ψ j (·; η) : η ∈ I, j ∈ N}. In keeping with solid state physics terminology, we will refer to the set I as reciprocal space (or more specifically, the Brillouin zone of the reciprocal space). Consequently, the dependence of a quantity on η will be termed as reciprocal space dependence while its dependence on usual physical space will be termed as real space dependence.
For a given j ∈ N, we will refer to the set Λ j = {λ j (η) : η ∈ I)} as a helical band. Results from the theory of regular perturbations of self-adjoint problems [95,96] imply that (for a fixed j ∈ N) the map η → λ j (η) is analytic. Therefore, the set Λ j is connected 23 The real numbers λ j (η) are generalized eigenvalues of H since (as discussed later) they are part of the essential spectrum of H and not its discrete spectrum. On a similar note, the functions ψ j (·; η) do not belong in L 2 (C) and therefore, they are not eigenfunctions of H in the usual sense. However, as discussed above, they do satisfy an equation of the form H ψ j (·; η) = λ(η) ψ j (·; η), thus suggesting their similarity to conventional eigenvalues and eigenfunctions. 24 The Dirichlet boundary condition in eq. 24 is also obeyed equivalently by ψ and does not need to be further considered here. and compact. 25 We will refer to the set Ψ = {ψ j (·; η) : η ∈ I, j ∈ N} as the collection of helical Bloch states corresponding to the helical bands. If we fix η ∈ I, then the set Ψ η = {ψ j (·; η) : j ∈ N} has the property that it is orthonormal and complete in L 2 (D). This follows directly from the properties of the group invariant functions φ j (x; η) defined above. Specifically, for j, j ∈ N: Furthermore, if h ∈ L 2 (D) such that h(·), ψ j (·; η) L 2 (D) = 0 for every j ∈ N, then we must have: Due to the completeness of the functions φ j (·; η) it then follows that h(x) e i2πη x 3 Due to the completeness of the set Ψ η for each η ∈ I, it actually follows that the set of helical Bloch states (i.e., the set Ψ) is complete in L 2 (C). To prove this important result, we first need to establish a few preliminaries related to the so-called helical Bloch-Floquet transform, i.e., an analogue of the classical Bloch-Floquet transform [86,88], as extended to the case of helical symmetries. Specifically, we show that there is a one-to-one correspondence between functions in L 2 (C) and L 2 (D × I) (this is the content of Lemmas 2.4 and 2.5), and we then identify the helical Bloch-Floquet transform as an operator 25 In contrast to the rigorous proof presented above, a formal derivation of the Bloch theorem for a helical structure, inspired by the solid state physics literature [38,39] is as follows: We observe that since the Laplacian commutes with all isometry operations -including those that constitute the group G 1 , and further, since the potential V (x) is group invariant (eq. 15), the operator H must commute with the symmetry operations in the group G 1 . Specifically, for any continuous function f defined over C, we may define the operators: Then, for any function f in the domain of H, the relationship T Υ Hf = HT Υ f holds for any T Υ ∈ T . This commutation property can be used to infer that the unitary representations of G 1 and the operator H can be "simultaneously diagonalized" in a suitable basis of common "eigenstates". Since G 1 is an Abelian group, its irreducible representations are all one-dimensional [42,97]. Furthermore, these irreducible representations can be used to decompose any unitary representation of the group [41,42]. This suggests therefore that the eigenstates associated with H transform under the group in a manner similar to the irreducible representations of G 1 , which then implies the helical Bloch theorem. While the above argument is perhaps correct in spirit and variants of the argument appear often in the physics literature (in the context of periodic systems) it has a number of technical deficiencies owing to the fact that H is an unbounded operator and the group G 1 is infinite. These issues prevent heuristic arguments like the one above -which are more suited to representations of finite groups on finite dimensional spaces -from being applied in the current context. In particular e.g., Bloch states are not eigenfunctions in the usual sense since they are not square integrable.
which maps between these spaces 26 . Thereafter, the completeness of Ψ η in L 2 (D) for each η ∈ I, in conjunction with the use of the helical Bloch-Floquet transform can be used to demonstrate the completeness of helical Bloch states in L 2 (C) (this being the content of Theorem 2.6). The completeness result of the helical Bloch states is intimately connected to the direct integral decomposition of the Hamiltonian, which we use for deriving the governing equations in the next section. We have: , η ∈ I and m ∈ Z. We define: Then g is defined almost everywhere in D and further, g ∈ L 2 (D × I).
This establishes that the function since h(x) = ∞ on a set of non-zero measure would violate eq. 32. Thus, for almost every x ∈ D, the sequence f m (x) m∈Z is square summable, and the expression for g in eq. 31 can be interpreted as a Fourier expansion (in the η variable). We may now use Parseval's identity [94] and eq. 31 to obtain: Integrating both sides of this expression for x ∈ D and using the steps in eq. 32 establishes that g ∈ L 2 (D × I), as required.
The following result is the converse of Lemma 2.4 and is established using the same tools as above: 26 Throughout this paper, we will often write functions in L 2 (D × I) as f (x, η) as well as f (x; η).
The latter notation is meant to emphasize such functions as being η-parametrized members of L 2 (D). However, as pointed out by an anonymous reviewer, it is perhaps not always possible to make this distinction consistently. 27 We would like to thank the anonymous reviewers for their comments which helped clarify and fix certain technical aspects of this proof, including the suggestion that Plancharel's Theorem [94] can be used to make certain statements in the above proof more precise.
Lemma 2.5. Let g ∈ L 2 (D × I) and for x ∈ D, m ∈ Z, let: Furthermore, let the function f be an extension off from the domain D to the domain C in the sense that for x ∈ D, Then, f ∈ L 2 (C).
Proof. By Tonelli's theorem [94], since g ∈ L 2 (D × I), it holds that g(x, ·) ∈ L 2 (I) for almost every x ∈ D. Then, we may interpret eq. 34 as a Fourier transform in η. By Parseval's identity [94], we have: Integrating both sides over x ∈ D and using g ∈ L 2 (D × I), we get: This shows that f ∈ L 2 (C), as required. Note that the interchange of the summation and the integral in the calculations above can be justified using the Fubini-Tonelli Theorem [94].
Lemma 2.4 establishes the existence of an operator U : L 2 (C) → L 2 (D × I) defined as: while Lemma 2.5 establishes the existence of its inverse U −1 : L 2 (D × I) → L 2 (C) defined as: To verify that eq. 39 indeed defines the inverse of the operator in eq. 38, we consider f ∈ L 2 (C) and g ∈ L 2 (D × I) such that g = Uf , i.e., We now multiply the above by e −i2πm η for m ∈ Z and integrate over η, to arrive at: Thus f = U −1 g in accordance with eq. 39.
We also observe, based on the calculations in eq. 37 that: and therefore, the operator U is an isometric-isomorphism 28 between the spaces L 2 (C) and L 2 (D × I). In analogy to the Bloch-Floquet transform in the literature used for studying periodic problems [86,88], we will refer to the operator U as the helical Bloch-Floquet transform 29,30 This operator allows us to demonstrate the completeness of the helical Bloch waves in L 2 (C). As mentioned earlier, the basic idea behind this proof is to map a given f ∈ L 2 (C) to its counterpart in L 2 (D × I) and to then use the completeness of the set Ψ η for each η ∈ I.
Proof. Since Uf ∈ L 2 (D × I), it follows from Fubini's theorem that Uf (·; η) ∈ L 2 (D) for almost every η ∈ I. Therefore, it can be approximated using the functions in the set Ψ η . Consequently, if we define: then g (·; η) → (Uf )(·; η) in L 2 (D) as → ∞ for almost every η ∈ I. In other words, the residual: 28 Eq. 42 shows that U is an isometry and Lemma 2.5 shows that it has a well defined inverse. Therefore, it is a unitary operator [94]. 29 This operator is closely related to the so-called Zak transform [98,99] associated with the group. 30 By use of the definition in eq. 38, it is easy to observe that the operator behaves in the following manner with respect to the action of the group: for any n ∈ Z.
has the property that r (η) → 0 for almost every η ∈ I, as → ∞. Furthermore using , as well as Bessel's inequality [94] on g → Uf , we get: However, (Uf )(·; η) 2 L 2 (D) is in L 1 (I) based on the calculations in Lemma 31. Therefore, by the Dominated Convergence Theorem [94]: and consequently: Since U is an isometric isomorphism, this implies that U −1 g → f in L 2 (C) as → ∞. Now, using eq. 39, we see that: On the other hand, evaluating eq. 44 at x = Υ m h • y, and using eq. 27, we see that: Comparing eqs. 50 and 52, it follows that U −1 g = f since y and m are generic, and therefore, f → f in L 2 (C) when → ∞, as required.
As mentioned earlier, the above results imply in essence that the spectral properties of H can be described completely in terms of helical bands and helical Bloch states (refer to Appendix B for further discussion along these lines). 31 Additionally, since the behavior of any helical Bloch state over all of C is completely specified based on its behavior over D (once a value of η ∈ I is chosen), the single electron problem posed on all of C 31 An immediate consequence of the completeness theorem for Bloch states is that the spectrum of H is completely contained in the set of helical bands, i.e., more precisely, spec.(H) ⊆ clos.(Λ), with clos.(·) denoting the (topological) closure. This is because, if κ ∈ R is such that κ / ∈ clos.(Λ), then the action of (H − κ) −1 on f ∈ L 2 (C) can be computed formally using eq. 44 in Theorem 2.6 as: can be reduced to a set of problems (indexed by η) posed on the fundamental domain (illustrated in Figure 2). Consequently, by appropriate use of the helical Bloch states and the helical bands, quantities of interest in Kohn-Sham theory (which can be described using the solutions to the single electron problem), can be formulated entirely in terms of quantities specified on the fundamental domain. We now look at this procedure in more detail.

Formulation of governing equations
In what follows, we will consider the helical structure to be at finite electronic temperature T e Kelvin and we will ignore spin polarization effects. For the sake of clarity of presentation, we will itemize the formulation/derivation of the various terms and equations, as we go along.
Electron Density and Density Matrix: A quantity of key importance in Kohn-Sham theory is the electron density. For a finite structure, such as a molecule or a cluster, this can be expressed in a straightforward manner in terms of the associated Kohn-Sham eigenstates and electronic occupations [33,63]. For a helical structure however, care has to be taken to express this quantity due to the fact that there are effectively an infinite number of electrons associated with the structure. In what follows, motivated by rigorous mathematical results related to the description of electronic states in crystalline systems [101][102][103], we address this issue by defining the single particle density operator [104][105][106] in terms of the single electron Hamiltonian, and then expressing the electron density in terms of the diagonal of the density operator.
To clarify the above procedure, let us first consider a finite system (i.e., an isolated molecule or an atomic cluster, for example) in R 3 , and let the single electron Hamiltonian, the single particle density operator (or density matrix), and the electron density for the system be denoted as H, D and (x) respectively. Then, H and D are related as [104-Now, using H ψ s (·; η) = λ s (η) ψ s (·; η), we have: so that: Since κ / ∈ clos.(Λ), the term 1 λs(η)−κ remains bounded even as s → ∞. Therefore, the right-hand side of eq. 55 can be interpreted as a bounded operator on f , and so, κ must belong to the resolvent set of H. Conversely, based on the techniques presented in [85,86] it is also possible to directly demonstrate that Λ ⊆ spec.(H), by constructing a suitable singular sequence of the form and using Weyl's criterion [95,96,100]. Here ζ(·) is a carefully chosen smooth cutoff function. Since, spec.(H) is always a closed set [95], and by definition, clos.(Λ) is the smallest closed set containing Λ, it follows that clos.(Λ) = spec.(H). Furthermore, if λ ∈ spec.(H), it can be immediately seen to be part of the essential spectrum of H. This is because if it were part of the point spectrum, then λ would be associated with eigenfunctions of finite multiplicity. Due to the fact that H commutes with the operators in T = T Υ : , these eigenfunctions would be left invariant by the operators in T as well (also see footnote 25). However, this would contradict the requirement that these eigenfunctions belong to L 2 (C). Note that these above results also follow from the direct integral decomposition of the Hamiltonian discussed in Section 2.2.2 and Appendix B.

Helical Bloch Theorem
Single electron problem posed on infinite cylinder Single electron problem reduced to fundamental domain

106]
: with f Te (·) denoting the Fermi-Dirac distribution function at electronic temperature T e , i.e.: Here, λ F and k B denote the Fermi level and the Boltzmann constant respectively. Due to this definition, D turns out to be a trace-class operator on L 2 (R 3 ), even though H is (generically) an unbounded self-adjoint operator on the same space. Assuming H has a pure point spectrum, denoting the eigenvalues of H as ω j ∈ R, and the corresponding eigenvectors as υ j ∈ L 2 (R 3 ), we may express H using its spectral representation as: Using this form we see that the action of H on any f ∈ Dom.(H) is expressible as: and that D can be expressed by means of spectral mapping [95] as: Due to this definition, the action of D on any f ∈ L 2 (R 3 ) can be expressed as: and it makes sense to write D in coordinate form as: for x, y ∈ R 3 . In this setting, the electron density (x) can be identified in terms of the diagonal of D, i.e., which leads to the well-known expression from Kohn-Sham theory [25] (also see footnote 33): Now, coming back to the case of the helical structure, we would analogously like to connect the single electron Hamiltonian H, the single particle density operator Γ and the electron density ρ(x). Accordingly, we define: as an operator on L 2 (C). The issue however, is that H does not admit a representation similar to eq. 58, and so the above definition does not immediately lead to transparent expressions for the electron density or the density operator in coordinate representation.
To adress this, it is useful to first recast eq. 65 in terms of helical Bloch states. The apparatus of direct integrals [88,107], discussed in Appendix B, allows us to do this in a mathematically rigorous manner. The key result from the appendix is that the helical Bloch-Floquet transform allows the single electron Hamiltonian to be "block-diagonalized" into a set of problems associated with the helical Bloch states that are posed over the fundamental domain, i.e.: Here, as before, H represents the operator − 1 2 ∆ + V (x) over the cylinder C along with the boundary condition ψ(x) = 0 for x ∈ ∂C. The potential V (x) is group invariant (eq. 15) and the unitary operator U : L 2 (C) → L 2 (D × I) represents the helical Bloch-Floquet transform (eq. 38). The operators {H η } η∈I represent the fibers of H (in the sense 22 of direct integrals) and are closely related to the operators h aux η introduced in the proof 32 of Theorem 2.3 (eq. 25). Specifically, for each η ∈ I, the operator H η represents the operator − 1 2 ∆ + V (x) over the interior of the fundamental domain (i.e., the set D) along with the boundary conditions ψ( The eigenstates of the operators {H η } η∈I are precisely the helical bands Λ = {λ j (η) : η ∈ I, j ∈ N} and the helical Bloch states Ψ = {ψ j (·; η) : η ∈ I, j ∈ N} restricted to the region D.
Above, eq. 66 expresses that H is unitarily equivalent to a "block-diagonal" operator whose "blocks" H η are indexed by η ∈ I. Therefore, upon computing Γ using eq. 65, we can expect to obtain another operator which is unitarily equivalent to a block-diagonal operator with blocks Γ η = f Te (H η ) η∈I . Since the function f Te (·) is analytic [108], these statements can be made mathematically precise by making use of the properties of the direct integral representation [88,Theorem XIII.85]. Thus, we may write for the density matrix Γ (as an operator on L 2 (C)): Next, if we are able to express f Te (H η ) in a transparent form, we may be able to further simplify the expression for Γ. Accordingly, we write the operators H η using spectral representation [95,109,110] as: and obtain: Thus, as an operator on L 2 (C), the density matrix admits the representation: While describing quantities over the fundamental domain, it is more appropriate and easier to deal with the density matrix as expressed as an operator on L 2 (D × I). This can be written asΓ = U Γ U −1 (i.e., the right hand side of eq. 67), and it admits the following expression in coordinate representation (with x, y ∈ D): 32 The key difference is that the operators h aux η include η dependence in the operators themselves and have group invariant solutions, whereas the operators H η include η dependence in the boundary conditions and have helical Bloch solutions (i.e., solutions which are group invariant up to an η dependent phase).

23
Now writing: we see that the electron density can be expressed 33 as (for x ∈ D): It is easy to see from the above expression 34 that the electron density is group invariant and also obeys a zero-Dirichlet boundary condition on the lateral surface of D (i.e., for x ∈ ∂D r=R ). It is also apparent from the expression 35 for the density matrix (eq. 71) that the following invariance relationship holds for any m ∈ Z and x, y ∈ D: For notational simplicity, it is convenient to introduce the scalars g j (η) = f Te λ j (η) , i.e. the thermalized occupation numbers of the electronic states of the system, that appear in eqs. 73, 71 and 70. We will denote the collection of occupation numbers as The electron density for an extended system is expected to obey the constraint of having a fixed number of electrons per unit fundamental domain of the system, even though the electronic states themselves are delocalized over the entire structure [26,44]. Denoting the number of electrons per unit cell as N e , in our case, this leads to: from which, using the orthonormality of the Bloch states over D, follows the constraint: 33 As pointed out by an anonymous reviewer, eq. 72 or (eq. 63 for the finite system case) can be somewhat subtle to interpret because the diagonal set (x, x) is of measure zero in R 3 × R 3 . As far as we can tell, it is quite common in the electronic structure calculations literature to define the electron density as the diagonal of the density matrix (see e.g. equation 22 in [106]), though this particular issue is never really addressed. The validity of this definition actually follows from the properties of nuclear operators (see [111,112] and references therein) and therefore, the case of the density matrices discussed here is covered. Alternately, if one were to accept the expression for the electron density as given in eq. 73 (or eq. 64 for the finite system case), then eq. 72 (or correspondingly eq. 63 ) can be seen as meaningful. 34 We would like to thank Eric Cances (Ecole des Ponts ParisTech) and Carlos Garcia Cervera (Univ. of California, Santa Barbara) for email communication related to technicalities of the above derivation of eq. 73 and also for providing useful references. 35 Note that the action of the operatorΓ on functions in L 2 (D × I) follows from the definition of the direct integral, specifically, eq. B.5 in Appendix B. Specifically, for f (x, η) ∈ L 2 (D × I):

24
In practice, the above equation can be used to compute the Fermi-level (λ F ) of the system. It also follows from this discussion that the density matrix Γ on L 2 (C) is locally trace class. Electronic Free Energy (per unit fundamental domain): With the above expressions in place, we can now use an energy-minimization formalism to deduce the governing equations of Kohn-Sham theory for the helical structure. Since the structure is infinite, the quantity of primary importance in this regard is the electronic free energy per unit fundamental domain, denoted here as F(Λ, Ψ, P G 1 , D, G 1 ). This notation emphasizes the dependence of this quantity on the helical Bloch states, the helical Bloch bands, the positions of the simulated atoms P G 1 = x k M G 1 k=1 within the fundamental domain, (the interior of) the fundamental domain and the helical group itself. Following [63], we express this quantity within the pseudopotential [113,114] and Local Density Approximations [25] as: with the terms on the right-hand side representing the kinetic energy of the electrons, the exchange correlation energy, the nonlocal pseudopotential energy, the electrostatic energy and the electronic entropy contribution, respectively. We now discuss each of the terms in the above equation in detail.
Kinetic Energy Term: The first term on the right-hand side of the above expression represents the kinetic energy of the electrons per unit fundamental domain. To motivate this term, we recall that for a finite system (in R 3 ) with single particle density matrix D, the kinetic energy can be expressed as [101,105,106]: with Tr.
[·] denoting the operator trace (of a trace-class operator on L 2 (R 3 )). Analogously, it would make sense to consider the trace per unit fundamental domain in case of the helical structure. As described in Appendix B, for an operator A which is invariant under the group G 1 and which is locally trace-class, it is possible to assign meaning to the trace per unit fundamental domain by means of the direct integral decomposition. Specifically, if the helical Bloch-Floquet transform block-diagonalizes the operator into its fibers as: then the trace per unit cell (denoted Tr. [·] henceforth) can be expressed as: with the trace inside the integral signifying the usual operator trace 36 in L 2 (D). The 36 The operator trace for any trace-class operator O on L 2 (D) can be computed as: where {f j } ∞ j=1 can be any orthonormal basis of L 2 (D). Refer e.g. to [101,105] for broader discussions of trace-class and locally trace-class operators in the context of electronic structure models. expression for the kinetic energy per unit fundamental domain for the helical structure therefore boils down to: 37 We now write − 1 2 ∆Γ η as (− 1 2 ∆) η Γ η , observe that Γ η is already available from eq. 69. Next, based on the discussion in Appendix B, we note that the fibers of the Laplacian on L 2 (C) are simply the Laplacian operators on L 2 (D) with (η -dependent) helical Bloch boundary conditions. Since the states ψ j (x; η) already satisfy these boundary conditions, and they form a basis of L 2 (D), it follows that: Exchange-Correlation Term: The second term on the right-hand side of eq. 78 represents the exchange correlation energy of the electrons per unit fundamental domain. Within the Local Density Approximation (LDA) [25], it can be written as: Note that it is also possible to modify this expression to use more sophisticated exchange correlation functionals such as the Generalized Gradient Approximation [115] and this will have little bearing on our subsequent discussion. Nonlocal Pseudopotential Energy Term: The third term on the right-hand side of eq. 78 represents the energetic contribution from the nonlocal part of the pseudopotential and models the effect of electronic core states. For a finite system of N at atoms located at the points {p k ∈ R 3 } Nat k=1 , if the single particle electron density is denoted as D, then this term has the following form: The operator V nl in the above equation is expressible in Kleinman-Bylander form [116] as : 37 Since the helical Bloch states ψ j (·; η) belong to the domain of the operator H η , it follows that the traces in eq. 83 are finite, making the expressions in that equation well defined. See [101,105] for further mathematical details along these lines.
Here, N k denotes the collection of projectors associated with the atom at p k , χ k,p are the projection functions, and γ k,p are the corresponding normalization constants. The functions χ k,p are themselves expressible in terms of atomic orbitals and are usually supported in a small region of space by design [117]. To obtain the correct analog of this expression for the helical structure i.e., the nonlocal pseudopotential energy per unit fundamental domain, it is useful to recall that this contribution to the energy is tied to the atoms in the fundamental domain as well as the electronic states in the system. It is of a somewhat different nature as compared to the kinetic energy term for instance, in which case the electrons are the only contributing source. Since the electrons in the extended structure are delocalized, the trace per unit fundamental domain leads to the appropriate expression in that case. In case of the nonlocal pseudopotential energy term however, the contribution from the electrons is delocalized, while those from the atoms are not. With this in mind, 38 we now focus on the atoms in the fundamental domain, and denote the non-local pseudoptential operator associated with these atoms as: Then, in analogy with eq. 86, the nonlocal pseudopotential energy per unit fundamental domain in case of the helical structure can be written by considering the action of V nl D to the density matrix operator defined in eq. 70, i.e.: To simplify this expression 39 , we employ eqs. 67, 70, the unitarity of the operator U, as well as the invariance of the trace under unitary transformations to obtain: Next, we observe that the operator V nl D = U V nl D U −1 acts on L 2 (D × I), and it admits a direct integral representation (i.e., V nl D = ⊕ I ( V nl D ) η dη ). The fibers of V nl D can be written as : In what follows, for the sake of brevity, we will denote the helical Bloch-Floquet transform of the projection functions, i.e., Uχ k,p (·; η; x k ) asχ k,p (·; η; x k ) and note that they can be 38 We would like to thank Phanish Suryanarayana, Georgia Institute of Technology, for discussions which helped clarify some of the properties of the nonlocal pseudoptential operator for the case of extended/condensed matter systems. 39 Note that Γ is locally trace class, while V nl D is a finite rank (and hence bounded) operator with a limited spatial extent. This makes eq. 89 well defined.
represented 40 via eq. 38 (for x ∈ D and η ∈ I as): With this notation, using properties tensor products, as well as eqs. 93 and 71, we may rewrite eq. 90 as: Now using the definition of the trace and that the helical Bloch states are a basis 41 , this 40 In practice, since the projection functions often have small support (centered about atomic positions), it is possible to truncate the above summation to just a few terms. Under certain circumstances, a somewhat more computationally convenient form forχ k,p (x; η; x k ) may be obtained by making use of the specific form of the projection functions χ k;p . The functions χ k;p are related to atomic orbitals and are therefore expressible as the product of a spherical harmonic with a compactly supported radially symmetric function. In the particular case that the projection functions arise from s-orbitals, it in fact follows that χ k;p (x; x k ) = χ k;p ( x − x k R 3 ). Then, we have: Thus, under these circumstances, the group action in the formula forχ k,p (x; η; x k ) has been shifted from x, to x k , which is easier to deal with computationally. 41 Note that for any fixed η ∈ I, {ψ j (·, η)} ∞ j=1 are a basis of L 2 (D). Additionally, as the proof of Theorem 2.6 shows, the entire set of helical Bloch states, i.e., Ψ = {ψ j (·; η) : η ∈ I, j ∈ N} forms a basis of L 2 (D × I). Thus, the trace of an operator O on L 2 (D × I) may be computed as: 28 reduces to: Electrostatic Energy Term: We now discuss the contribution of the electrostatic interaction energy to the free energy per unit fundamental domain. This is the fourth term on the right-hand side of eq. 78. Often, it is computationally advantageous to express this term using a so-called local formulation [63,[118][119][120], as we now do 42 . For a finite system with atomic nuclei located at the points {p k ∈ R 3 } Nat k=1 and electron density ρ, this term takes the form of the following optimization problem in the total electrostatic potential: Here b finite represents the total nuclear pseudocharge for the finite set of nuclei, and can be expressed in terms of the individual nuclear pseudocharges b k (x; p k ) Nat k=1 as: and the term E finite sc (p 1 , p 2 , . . . , p k ) corrects for self-interactions and overlaps of the nuclear pseudocharges [119]. Note that by design, the individual nuclear pseudocharges are usually smooth, radially symmetric functions centered at the nuclear positions, they have compact support, and they integrate to the (valence) nuclear charge of the nucleus in question. The electrostatic potential Φ that solves the maximization problem in eq. 97 is the Newtonian potential associated with the net charge in the system: 43 This is essentially the calculation described in eqs. 94, 96 above, with the operator O = UV nl D U −1Γ on L 2 (D × I). Note that the above results also directly follow from the properties of direct integral decomposition discussed in Appendix B. 42 The term local formulation is associated with the fact that the electrostatic potential Φ can be solved through a Poisson equation, which avoids evaluation of the non-local integrals in eq. 99. 43 Note that we have made a minor abuse of notation and used Φ to denote the "trial" electrostatic potentials involved in the maximization problems in eqs. 97 and 102, as well as the actual potentials that achieve these maxima (i.e., the arg max of the functionals listed in eqs. 97 and 102.). The latter are expressible succinctly as the corresponding Newtonian potentials in eqs. 99 and 101.
To extend the above formulation to a helical structure, we first write the total nuclear pseudocharge at any point x ∈ C in terms of the pseudocharges of the atoms in the fundamental domain as: and observe that this quantity is group invariant owing to the aforementioned properties of the individual nuclear pseudocharges [40]. Since the electron density is group invariant as well, it follows that the net electrostatic potential Φ expressed as 44 : is also group invariant [8]. Therefore, we may use Φ, b and ρ as defined over the fundamental domain, to define the analog of eq. 97 for the helical structure as: For the sake of brevity, we omit the details of the form of the corrections due to self interactions and overlaps of the nuclear pseudocharges, as reduced to the fundamental domain (i.e., the term E sc (P G 1 , G 1 , D) above) and instead point to [40,119] for relevant information. Electronic Entropy Term: Finally, the last term on the right-hand side of 78 represents the electronic entropy contribution to the free energy. Following [121], this term can be expressed for a finite system with density matrix D as: with I denoting the identity operator. To obtain the analogous expression for the case of the helical structure, we work with the trace per unit fundamental domain instead. This gives us: By means of spectral mapping [95], the use of eq. 69, and by noting again that the helical Bloch states ψ j (x; η) are a basis of L 2 (D), we readily obtain: 44 Using a Fourier series expansion, it can be shown (see e.g. [44] for similar arguments) that eq. 101 is well defined whenever the system is charge neutral i.e., when With the above terms explicitly defined, we now turn to the variational problem for deducing the governing equations. Variational Problem and Kohn-Sham Equations: The variational problem for Kohn-Sham ground state of a given helical structure (i.e., the atomic coordinates in the fundamental domain, the fundamental domain and the helical group are held fixed) consists of minimizing the electronic free energy F(Λ, Ψ, P G 1 , D, G 1 ) with respect to the helical Bloch states and the helical bands, subject to the constraint in eq. 77. In the literature, this minimization is often stated in terms of the helical Bloch states and the electronic occupation numbers instead [34,63]. Along those lines, we may define F(G, Ψ, P G 1 , D, G 1 ) = F(Λ, Ψ, P G 1 , D, G 1 ) to write the variational problem as: subject to: and the requirement that the states in Ψ be helical Bloch states. This requires that for any ψ j (·; η) ∈ Ψ and m ∈ Z, we have ψ j (Υ m h • x; η) = e −i2πmη ψ j (x; η) , as well as the orthonormality condition between two helical Bloch states ψ i (·; η), ψ j (·; η) ∈ Ψ: We take variations of the above constrained minimization problem, and obtain the Euler-Lagrange equations as the following helical symmetry adapted Kohn-Sham equations over the fundamental domain: Here, the helical symmetry adapted Kohn-Sham Hamiltonian operator (with its dependence on the helical Bloch sates, the occupation numbers, etc., made explicit) 45 is: in which V xc = δE xc (ρ, D) δρ is the exchange correlation potential, Φ is the net electrostatic potential and satisfies the following symmetry adapted Poisson problem over the fundamental domain: and the operator V nl D = U V nl D U −1 is as defined 46 in eqs. 91 and 93. 45 Up to a notational change, this operator is essentially the same as H KS [Λ,Ψ,P G 1 ,D,G1] , when the dependence on helical bands (instead of the occupation numbers) is highlighted. 46 Due to the propertyχ k,p (Υ n h • x; η; x k ) = e −i2πnηχ k,p (x; η; x k ) as observed in footnote 30, it follows that V nl D commutes with the group action in an appropriate sense (when viewed as an operator on L 2 (D × I)). This implies that the Kohn-Sham operator H KS commutes with the group action as well.
Harris-Foulkes Functional: The above set of expressions represent a set of coupled nonlinear partial differential equations in the fields ψ j (·; η) and the scalars g j (η). Once they have been solved self-consistently, the ground state electronic free energy F 0 (P G 1 , D, G 1 ) per unit fundamental domain can be computed through eq. 78. In practical calculations, since self-consistency is never achieved perfectly, a better estimate of the ground state electronic free energy may be found using the so-called Harris-Foulkes functional [122,123]. This can be written in helical symmetry-adapted form, using quantities expressed over the fundamental domain as: All the quantities on the right-hand side of the above equation are easily interpreted based on earlier discussion, except the first one, i.e., E band (Λ), which represents the electronic band energy per unit fundamental domain. For a finite system with a single electron Hamiltonian H and single particle density matrix D, this quantity is expressed as [106]: Analogously, for the helical structure, we use the trace per unit fundamental domain to write: Using eqs. 68 and 69 and using the completeness of the helical Bloch waves, we see that this is be expressible as: Atomic Forces: The Hellmann-Feynman forces on the atoms in the fundamental domain are (in Cartesian coordinates): By directly differentiating the various terms involved (eqs. 94, 102) we arrive at the 32 following expression for f k using quantities specified over the fundamental domain 47 : Here, Re. denotes the real part of the quantity in braces. This completes a discussion of the derivation of the various physically relevant terms, as well as the form of the equations of Kohn-Sham theory, as applied to a helical structure associated with a helical group generated by a single element. Comments on modifications of the above equations while dealing with a structure associated with a helical group generated two elements, and a presentation of the final expressions/equations for that case in appear in Appendix C.

Numerical Implementation
The Kohn-Sham equations for a helical structure (i.e., eq. 109 or eq. C.15) are a set of non-linear eigenvalue problems indexed by η (as well as ν in case of the group G 2 ) that are coupled to each other through the electron density ρ. The standard procedure for solving the equations of Kohn-Sham theory is through self-consistent field (SCF) iterations [25]. This amounts to starting from a reasonable guess of the electron density in the fundamental domain (e.g. superpositions of individual atomic densities, as is used in our simulations) and an appropriate set of trial orthonormal wavefunctions (randomly chosen in our simulations), and then evaluating the eigenstates of the Kohn-Sham operator with these guesses. Thus, a set of linear eigenvalue problems (i.e., those associated with the linearized Kohn-Sham operator evaluated at the given electron density) indexed by η (as well as ν in case of the group G 2 ) have to be solved. From this, the (trial) Fermi-level of the system and the (trial) occupation numbers maybe computed. The eigenfunctions and the occupation numbers may be then combined (in accordance with eq. 73 or eq. C.3) to yield the trial electron density for the next step of the iterations. The above procedure has to be repeated till the difference in the electron density (or the effective potential) between successive iterations reaches below the desired convergence threshold. We will now discuss several features of this self-consistent solution process as implemented in the Helical DFT code. 47 Motivated by [63,124] we may use integration by parts to modify the last term on the right-hand side of eq. 118, so that the derivatives of the projectors with respect to atomic coordinates can be eliminated in favor of Cartesian gradients of the wavefunctions instead. This tends to improve the accuracy of the computed forces in practical calculations -the orbitals are more smoothly varying than the projectors and therefore they tend to behave better upon taking derivatives. With this change as well as making use of the discussion in Footnote 40, it is possible to rewrite: under specific circumstances.

Discretization of reciprocal space
Many quantities described in Section 2.2.2 and Appendix C involve integrals over η ∈ I = [− 1 2 , 1 2 ) (as well as normalized summations over ν ∈ {0, 1, 2, . . . , N − 1} for the group G 2 ). To evaluate such integrals numerically, we employ quadrature based on the Monkhorst-Pack scheme [125]. Specifically, we sample the interval I using a grid of N η points, and write: Here, w b and η b denote the integration weights and integration nodes respectively. Summations over ν are left unchanged. The total number of points used for discretizing the reciprocal space (i.e., set B = I × {0, 1, 2, . . . , N − 1}), therefore, is N K = N η × N (with N = 1 for the group G 1 ). Based on considerations of time-reversal symmetry (which apply as long as e.g. magnetic fields are absent) [63,126], it follows that for η ∈ I and ν = 1, 2, N − 1: while for ν = 0: Effectively, the above considerations reduce the number of quadrature points over reciprocal space by half (i.e., N K ≈ N η × N 2 ). With the above discretization choices, the self-consistent field iterations for the Kohn-Sham problem amount to solving a series of N K linear eigenvalue problems on every iteration step. Based on the mathematical treatment presented earlier as well as in [40], it follows that eigenvalue problems associated with distinct values of η (i.e. η b in discretized form) and/or ν are disjoint from each other. This implies that these (linear) eigenvalue problems can be solved independently of each other, in an embarrassingly parallel manner, regardless of how the discretization in real space is carried out. We make use of this feature of the equations to assign these distinct eigenvalue problems to different computational cores. This serves as a natural parallelization scheme and helps in drastically reducing the wall time associated with the most computationally intensive part of the SCF iterations.

Truncation of infinite sums
Primarily, there are two distinct sources of infinite sums in the equations presented in Section 2.2.2 and Appendix C. The first arises due to summing over an infinite number of helical bands (e.g., eqs. 73 and C.3). To truncate such sums we assume that the electronic occupation numbers reduce to zero beyond the lowest N s bands and therefore, only N s eigenstates for each value of η b (and also each ν for G 2 ) need to be computed during the self-consistent field iterations. In effect, this is also an enforcement of the Aufbau principle for the system [26,44]. Depending on the size of the discretized reciprocal space (i.e., the value of the number N K ) we have found that including just a few extra bands beyond the minimum number required for holding the N e electrons per unit fundamental domain, suffices. 48 The second source of infinite sums arises from considering terms associated with group orbits (e.g., eqs. 100 and C.12), since helical groups by definition are infinite. However, these sums are also always associated with functions that are supported in a small ball around an atom of the structure (e.g. the nuclear pseudocharge in eq. 100 and the nonlocal pseudopotential projection function in eq. 93). Therefore, the influence of such sums on points in the fundamental domain is only dependent on terms of the summation that result in a nonzero overlap between the function support and the fundamental domain. This allows such infinite sums to be truncated as well.

Helical coordinate system
To carry out a discretization of the governing equations in a manner that is naturally adapted to the underlying helical symmetries of the structures being studied, it is useful to employ helical coordinates, as introduced in [8]. In order to have this coordinate system be commensurate with the helical groups G 1 or G 2 , the coordinate transformation formulae are as follows. For x ∈ C, if the Cartesian coordinates are (x 1 , x 2 , x 3 ), then the corresponding helical coordinates (r, θ 1 , θ 2 ) are 49 : The coordinates (r, θ 1 , θ 2 ) are a natural generalization of cylindrical coordinates in the sense that they effectively reduce to cylindrical coordinates when the twist angle parameter α of the system is set to zero. We may verify that these relations are onto and globally invertible on C\{t e 3 : t ∈ R}. Furthermore, the inverse relations: (r, θ 1 , θ 2 ) → (x 1 , x 2 , x 3 ) = r cos(2π(αθ 1 + θ 2 )), r sin(2π(αθ 1 + θ 2 )), θ 1 τ The action of a helical group can be easily computed in these coordinates as follows. Let x = (x 1 , x 2 , x 3 ) ∈ C have helical coordinates (r, θ 1 , θ 2 ). The action of the isometry Υ h that generates G 1 (and also G 2 ) is to map x to the point x = Υ h • x = R 2πα x + τ e 3 . Denoting the Cartesian and helical coordinates of this new point as (x 1 , x 2 , x 3 ) and (r , θ 1 , θ 2 ), respectively, we see that x 1 = x 1 cos(2πα) − x 2 sin(2πα), x 2 = x 2 cos(2πα) + x 1 sin(2πα) and x 3 = x 3 + τ . Now, using eqs. 123 and 122, we get r = r, θ 1 = θ 1 + 1, θ 2 = 48 This is a well used approximation strategy in the literature (see e.g. [33,34,37]). For finite systems at electronic temperatures that are less than a few thousand Kelvin, it suffices to choose the number of states to be equal to a multiple of half the number of electrons, with the multiplication factor being between 1.05 and 1.10 [127]. For an extended system like a helical structure, often a just a few extra bands beyond half the number of electrons is sufficient since this actually amounts to these few extra states being available for every value of η b or ν. As a result of this, tens or even hundreds of extra states (with occupation numbers approaching zero) get effectively included in the calculations. 49 For computational purposes, it is more apt to use arctan2(x 2 , x 1 ), along with a suitable modification for the negative Y quadrants, instead of arctan ( x2 x1 ) in eq. 122, such that a value between 0 and 2π is obtained. θ 2 . By a similar calculation, we see that the action of the second generator Υ c of the group G 2 (i.e., the pure rotation by 2π N about axis e 3 ), is to map the point with helical coordinates (r, θ 1 , θ 2 ) to the point (r, θ 1 , θ 2 + 1 N ). These calculations imply in particular that if a function is invariant under the group G 1 , then it is periodic in θ 1 , with period 1, when expressed in helical coordinates. Similarly, invariance of a function under the group G 2 , implies periodicity in θ 1 (with period 1), as well as periodicity in θ 2 , (with period 1/N) when expressed in helical coordinates. These observations make it easy to enforce helical Bloch boundary conditions on the wavefunctions, as well as the group invariance of the electrostatic potential, in simulations. Derivation of the Cartesian gradient operator, the Laplacian operator and the volume integral in helical coordinates appears in Appendix D.

Real space discretization: Finite difference scheme
We employ a finite difference strategy for discretizing the governing equations in real space. We choose an annular cylindrical region Ω with axis along e 3 as the simulation domain. This allows us to avoid the singularity associated with the helical coordinate system along the axis e 3 and does not present any issues as long as the atoms of the structure are located well within the annular region [40,63]. This latter condition is well satisfied by the nanotube structures simulated in this work. The set Ω can also serve adequately as a suitable fundamental domain for either group G 1 or G 2 in simulations (i.e., it can replace D or D in the formulae presented in Section 2.2.2 and Appendix C). In cylindrical coordinates (r, ϑ, z) we have: The boundary of Ω can be expressed as: with ∂R in and ∂R out denoting the surfaces r = R in and r = R out respectively, ∂ϑ 0 and ∂ϑ Θ denoting the surfaces ϑ = 2παz τ and ϑ = 2παz τ + Θ respectively, and finally, ∂Z 0 and ∂Z τ denoting the surfaces z = 0 and z = τ respectively. Figure 3 illustrates the simulation cell as well as the boundaries of this domain for an untwisted helical structure.
We set up a finite difference grid over Ω using helical coordinates, with spacing h r , h θ 1 and h θ 2 along the r, θ 1 and θ 2 directions respectively. Since the simulation domain can be represented in helical coordinates with r ranging from R in to R out , θ 1 ranging from 0 to 1, and θ 2 ranging from 0 to 1 N , it follows that R out − R in = N r h r , 1 = N θ 1 h θ 1 , and 1 N = N θ 2 h θ 2 , for natural numbers N r , N θ 1 and N θ 2 . We index the finite difference nodes using a triplet of natural numbers (i, j, k), for i = 1, 2, . . . , N r , j = 1, 2, . . . , N θ 1 and k = 1, 2, . . . , N θ 2 . We denote the value of a function f at the grid point (i, j, k) as f (i,j,k) .
We will denote h = Max. h r , τ h θ 1 , R in +Rout 2 h θ 2 as the mesh spacing.
The formulae presented in Section 2.2.2 and Appendix C require Cartesian gradients, the Laplacian operator and integrals over Ω to be evaluated using the finite difference scheme. The formulae for these quantities, as expressed in helical coordinates appears in Appendix D. With these in hand, we approximate the first order partial derivatives using central differences as: ∂f ∂r We approximate the pure (i.e., non-mixed) second order partial derivatives as: , while the mixed second order partial derivative is obtained by applying the above first order derivative formula first in θ 1 and then in θ 2 , i.e., In the above formulae, n o denotes half the finite difference order and is set to 6 for all our simulations (i.e., 12 th order finite differences are used). This choice has also been employed elsewhere [33,34,40,63] and is found to be adequate for attaining chemical accuracy. Letting s denote r, θ 1 or θ 2 , the weights that appear in the above formulae can be expressed as [128]: We employ the following quadrature rule for approximating integrals over Ω: with r i denoting the radial coordinate of the finite difference node (i, j, k).

Boundary conditions
We need to specify the boundary conditions on the various fields that are being solved for in the governing equations. These are the helical Bloch states ψ j (x; η, ν) (for η ∈ I and ν = {0, 1, 2, . . . , N−1}) which satisfy the Kohn-Sham equations (eq. 109 or eq. C.15), and the total electrostatic potential Φ which satisfies Poisson's equation (eq. 111 or eq. C.20). In helical coordinates, we may interpret the helical Bloch boundary conditions (eq. C.1) as the following conditions: which applies to the surfaces ∂Z 0 ∂Z τ ; as well as the condition: which applies to the surfaces ∂ϑ 0 ∂ϑ Θ . We assume that the atoms within Ω are sufficiently far away from the boundary surfaces ∂R in and ∂R out , so that the electron density decays to zero at these surfaces and zero Dirichlet boundary conditions on the wavefunctions (Section 2.2.1) can be applied. Thus, the conditions to be applied on the surfaces ∂R in ∂R out are: As already discussed, the electrostatic potential Φ inherits the symmetry of the helical structure (i.e., it is invariant under G 1 or G 2 ). Thus, it follows that the boundary conditions on this quantity on the surfaces ∂Z 0 ∂Z τ and ∂ϑ 0 ∂ϑ Θ are respectively: To apply the right boundary conditions on Φ on the surfaces ∂R in ∂R out , we may evaluate eq. 99 or eq. C.13 directly using Ewald summation [74] or multipole expansion [91] techniques. 50

Numerical linear algebra issues
On every SCF iteration step, the lowest N s eigenstates of a set of N K Kohn-Sham operators (indexed by η b for the case of G 1 , and η b , ν for the case of G 2 ) have to be computed. We employ iterative diagonalization based on Chebyshev polynomial filtered subspace iterations (CheFSI) [129][130][131][132] for this purpose. Due to the fact that the helical coordinate system is curvilinear, the finite difference discretization of the Laplacian operator results in a discretized operator that is non-Hermitian (even though the Laplacian as an operator on an infinite dimensional Hilbert space is Hermitian). This results in the discretized Hamiltonian operator also being non-Hermitian, which can lead to non-real eigenvalues of the matrix. However, as the discretization is made finer (i.e., h r , h θ 1 , h θ 2 become smaller) the discretized Laplacian and Hamiltonian matrices approach Hermitian matrices and so, the eigenvalues resulting from these discretized operators tend to have vanishingly small imaginary parts [40,133]. In practice, for the mesh spacings that have been considered in this work, the imaginary parts are small enough that they can be safely ignored without affecting the quality of the simulations (also see [40] and [63] where a similar situation was encountered).
We use Chebyshev polynomial filter orders in the range 55 to 80 for our simulations. This is somewhat higher than what is commonly employed in finite difference DFT calculations in affine coordinate systems with comparable mesh spacing [33,34,130]. We adopt it here to mitigate the effect of the larger spectral width of the discretized Hamiltonian that arises due to crowding of grid points as one approaches the origin in helical coordinates. We have also observed that the time for computing matrix-vector products using the discretized Hamiltonian in helical coordinates is larger compared to the time required for the case of a finite difference Hamitonian (of the same size) arising from a cylindrical system (obtained by setting the twist angle parameter α to 0 in the helical case). This is certainly due to the presence of cross derivatives in the helical case (eq. 128), which makes the discretized Hamiltonian somewhat less sparse when compared with the cylindrical one. Within the CheFSI method, we use Arnoldi iterations for computing the extremal eigenvalues of the Hamiltonian and a direct diagonalization method for computing the projected subspace problem.
We use the Generalized Minimal Residual method (GMRES) [134] to solve the Poisson problem associated with the total electrostatic potential Φ. To accelerate convergence, we use an incomplete LU factorization based preconditioner [135].

Matlab implementation: The Helical DFT code
We have implemented the above computational strategies using the MATLAB [136] software package into a code called Helical DFT. 51 The code parallelizes computation over the different η and ν values using MATLAB's Parallel Computing Toolbox (the parfor function). For maximum efficiency of the MATLAB implementation, code vectorization has been used as much as possible. However, for computing certain quantities (such as the atomic forces and the net nuclear pseudocharge), multiple levels of nested loops were found unavoidable. These routines were converted into machine code by use of the MATLAB Coder framework, which helped alleviate performance issues. In order to reduce the memory footprint associated with the storage of the different Hamiltonian matrices arising from the N K different values of η b (and also ν for G 2 ), we avoid computing matrix vector products through MATLAB's internal sparse matrix framework since that requires these matrices to be available explicitly. Instead, we store only the Laplacian part of the Hamiltonian matrix in compressed sparse row (CSR) format and apply the helical Bloch boundary conditions associated with the different values of η b and/or ν on the fly, while computing matrix vector products. The actual task of computing these matrix vector products is carried out using a C language routine which has been compiled and interfaced with our MATLAB code.
Some other relevant details related to the implementation are as follows. We use the periodic variant of Pulay's scheme [137,138] in the total potential to accelerate the convergence of the SCF iterations. The Fermi energy is calculated using a nonlinear equation root finder (MATLAB's fzero function). When required, structural relaxation is achieved by an implementation of the Fast Intertial Relaxation Engine (FIRE) algorithm [139].

Simulation Results and Discussion
We now turn to a discussion of numerical simulations and results. All simulations were run using a single node of the Mesabi cluster at the Minnesota Supercomputing Institute, or a single node of the Hoffman2 cluster at UCLA's Institute for Digital Research and Education. Each compute node of Mesabi has 24 Intel Haswell E5-2680v3 processors operating at 2.50 GHz, and 64 GB to 1 TB of RAM. Each compute node of the Hoffman2 cluster has two 18-core Intel Xeon Gold 6140 processors (with 24.75 MB cache, and running at 2.3 GHz), and 192 GB of RAM.
All calculations presented here use Troullier-Martins norm conserving pseudopotentials [113]. The Local Density Approximation [25] was used for modeling the exchange correlation energy, and the Perdew-Wang parametrization [140] of the correlation energy was employed. An electronic temperature of T e = 315.77 Kelvin was used for Fermi-Dirac smearing to help accelerate SCF convergence.
The large majority of the simulations here have focused on the study of single wall nanotubes. Starting from the sheet of an elemental two-dimensional material, nanotubes of any chirality can be formed using the so-called "roll-up" construction [141]. This procedure also allows us to see [7,49] that such nanotubes can be adequately represented using helical groups generated by two elements with just 4 atoms in the fundamental domain [63]. In this representation, the twist angle parameter α becomes related to the chirality of the tubes (α = 0 for achiral tubes), while the parameter N, associated with the cyclic group order, is related to the tube radius [49]. In our nanotube simulations, we have ensured that the atoms within the fundamental domain are always located 10 to 12 Bohrs away from the boundary surfaces ∂R in and ∂R out so as to allow sufficient decay of the electron density and the wavefunctions in the radial direction.

Materials system: Single layer black phosphorus nanotubes
Single-layer black phosphorus, or phosporene, is a two-dimensional nanomaterial that has been the object of intense investigation in recent years due its association with a number of unusual and fascinating material properties [142][143][144][145][146][147][148]. Nanotubes of this material, as formed by the roll-up construction have also received attention in the literature [149][150][151][152][153][154][155][156][157][158][159][160][161][162], due to their interesting optical and electronic properties, and the coupling of these properties to mechanical strains. This motivates our choice in selecting this material for the simulations presented in this work.
As a starting point, we obtained the ground state structure of a single layer of phosphorene ( Figure 4) using the same pseudopotential, exchange correlation functional and electronic temperature as the Helical DFT simulations subsequently described. We used the plane-wave DFT code ABINIT [30,163] to perform the geometry relaxation calculation. The periodic unit cell for this simulation contained 4 atoms. An energy cutoff of 40 Ha, along with 30 × 30 × 1 k-points and a cell vacuum of 25 Bohr in the Z-direction was used. At the end of this very refined calculation the atomic forces were all less than 10 −5 Ha/Bohr, while the cell stresses were of the order of 10 −8 Ha/Bohr 3 or lower. Some of the structural parameters obtained from the calculation are shown in Table 1. There appears to be generally good agreement with the literature 52 thus giving us confidence in the reliability of the subsequent simulations with regard to materials physics.
To represent the zigzag phosphorene nanotubes in the Helical DFT simulations described next, we roll up the phosphorene sheet along the X-axis and place the atoms from the aforementioned periodic unit cell (the shaded region in Figure 4) into the Helical DFT simulation cell (Ω). In the absence of relaxation effects, the pitch τ associated with the two-generator helical group, is equal to the lattice constant along the Y axis. Furthermore, the angle Θ associated with the cyclic group order N is related to the radius of the     . Here a 1 denotes the lattice vector along the X-axis in the phosphorene sheet, and R avg. denotes the average radial coordinate of the atoms in the fundamental domain of the nanotube. For chiral nanotubes, we additionally include a no-zero twist angle parameter α. Figure 5 shows examples of two phosphorene nanotubes studied in this work using Helical DFT.

Convergence and accuracy
To study the convergence properties of our numerical implementation, we consider a chiral phosphorene nanotube of radius 1.76 nanometers. The cyclic group order is N = 32 and the twist angle parameter is α = 0.0025. For a given structure and a fixed simulation domain, the two convergence parameters under study are the mesh spacing h = Max. h r , τ h θ 1 , R in +Rout 2 h θ 2 , and the number of points N η used for discretizing the set I in reciprocal space. As a reference calculation, we computed the electronic structure of the system with the finest mesh (h = 0.23 Bohr) and the highest value of N η (= 23) that we could afford under computational resource constraints. First, for the mesh convergence study, we fix N η = 23 and carry out a series of calculations with h ≈ 0.30, 0.40, 0.50, 0.60, 0.65, 0.70 Bohr. Next, for studying convergence with respect to N η , we fix h = 0.23 Bohr and carry out a series of calculations with N η = 1, 3, 7, 11, 15, 19. We plot the errors in the energy per atom and the atomic forces in each of the above  cases in Figure 6.
From the figures, it is clear that the code converges to the reference calculations systematically. Using straight-line fits to the convergence data with respect to h, we find slopes of 5.8 and 7.6 for the energies and forces, respectively. These numbers are very nearly identical to the convergence rates observed in finite difference calculations involving cylindrical coordinates [63] and are also comparable to finite difference calculations in affine coordinate systems [33,34]. From the data, we are also able to estimate that the parameters h = 0.40 Bohr and N η = 11 are more than sufficient to reach chemically accurate energies and forces (the exact errors for these choices with respect to the reference calculation were 2 × 10 −4 Ha/atom and 8 × 10 −5 Ha/Bohr in the energies and forces, respectively). For computational efficiency therefore, we use this set of parameters for all relaxation calculations described subsequently, and switch to the reference calculation parameters (i.e., h = 0.23 Bohr, N η = 23) only when more accurate energies / band gaps are required at the end of a relaxation procedure. 53 Verifying the accuracy of the Helical DFT code with respect to standard plane-wave codes can be challenging since the plane-wave codes may be required to include an enormous number of atoms in the periodic unit cell in order to replicate the exact system being studied by Helical DFT. In case of the above chiral nanotube system for example, over 50, 000 atoms would be needed, making the planewave calculation unfeasible.
Therefore, we carry out this accuracy check in two steps. First, we set up a Helical DFT calculation for a zigzag nanotube (i.e., α = 0) of radius 0.94 nanometers. For this tube, the cyclic group order N = 16. Then, we simulate this tube using the ABINIT code by employing a 64 atom unit cell (periodicity was enforced along the Z axis, Dirichlet boundary conditions were enforced along X and Y axes by padding with a large amount of vacuum). We converged both codes to the extent allowed by computational resources, and observed that the energies (in Ha/atom) and the forces (in Ha/Bohr) from these two calculations agreed with each other to 1 × 10 −4 or better. Next to study a case for which α = 0, we set up an artificial system consisting of atoms along a single helix (similar to the configuration in Figure 2). This system was generated using a single generator helical group (N = 1) and used α = 0.01. The helical unit cell had 2 atoms, while the periodic unit cell in ABINIT involved 200. Once again, upon convergence with respect to their respective discretization parameters, the codes produced results that differed from each other by about 1 × 10 −4 Ha/atom or Ha/Bohr. This completes the accuracy tests. 54 Through the above examples, we were also able to observe that for realistic helical nanostructure simulations, the wall time for Helical DFT can be up to orders of magnitude smaller compared to a well optimized plane-wave code like ABINIT, making it a powerful first principles tool in the study of such systems.    54 It is likely that Helical DFT can be made to agree with ABINIT results to finer levels of accuracy.
However, the quasi-one-dimensional nature of the systems being studied, and the slow convergence of the electrostatics requires the use of large amounts of vacuum padding in the ABINIT supercell, and this tends to cause serious convergence issues as the energy cutoff is increased.

Simulation of twisting: Ab initio computation of torsional stiffness
By using the two-generator helical group G 2 , it is possible to describe torsional deformations in nanostructures [49,55,[166][167][168]. We use this procedure here 55 to illustrate the utility of Helical DFT in extracting mechanical properties of nanomaterials ab initio. Specifically, we investigate the behavior of a zigzag phosphorene nanotube (nanotube radius 3.4 nanometer, cyclic group order N = 64.) under twisting deformations. We begin with the untwisted structure (α = 0) and relax the atoms in the simulation cell till all force components on all atoms are below 10 −3 Ha/Bohr. 56 Starting from this relaxed structure, we apply twisting deformations to the nanotube by prescribing non-zero values of the twist angle parameter α in the two-generator helical group G 2 used for describing the nanotubes. We varied α in steps of 0.001 and relaxed the resulting nanotube structure in each case, till all components of the forces on all the atoms in the simulation cell dropped below 10 −3 Ha/Bohr once again. The largest twist we considered is about 5 degrees per nanometer length of the tube, which corresponds to α = 0.006. Anticipating a quadratic dependence of the nanotube twist energy for small values of alpha [49], we plot the results as shown in Figure 7. In that figure, we have also plotted the twisting energies obtained when the atomic relaxation effects are not considered.
We write the twist energy per unit length of the tube as U twist = 1 2 k twist β 2 , with β denoting the twist per unit length of the tube (β = 2πα τ ) and k twist denoting the torsional stiffness. From the figure, we use the straight line fits near zero twist to estimate the value of k twist as 13.66 eV nm (deg) −2 and 27.33 eV nm (deg) −2 , for the relaxed and unrelaxed cases respectively 57 . It is also evident from the figure that non-linear effects start to play a noticeable role in this nanotube at around 4 degrees of twist per nanometer. These simulations serve as an example of determining constitutive parameters directly from quantum mechanics using Helical DFT.

Electronic properties: Behavior of zigzag tubes under torsional deformations
Since Helical DFT is an ab initio simulation tool, it can give us insights into the electronic, optical and transport properties of materials under study. In particular, it may be used to shed light into the coupling of mechanical strains with these properties. To illustrate these points, we consider again the case of the zigzag phosphorene nanotube discussed above (nanotube radius 3.4 nanometer, cyclic group order N = 64.). Figure 8 displays examples of helical band structure diagrams, which like their periodic counterparts, can be used to illustrate the electronic levels in the system. A unique feature of these diagrams however, is that they can be used to display the variation of Kohn-Sham eigenvalues with respect to η and ν. This makes them significantly easier to interpret 55 Due to the use of helical symmetries, the structures being modeled in the simulations are infinite. In practice, nanotube structures have finite extent, though some can have lengths of the order of macroscopic sizes [169]. The edge effects in these materials are expected to decay as one moves towards the interior of the material, both in the continuum elasticity sense [170,171] and at the level of the electronic structure [172]. This offers a justification for the conceptual correctness of the simulations. 56 Due to the relatively large nanotube radius, the atoms in the untwisted nanotube are in an environment similar to that in the phosphorene sheet. Since the phosphorene sheet itself had been relaxed well, the forces on the atoms in the nanotube were relatively small to begin with, and structural relaxation was usually completed in just a few steps. 57 These values implicitly use β in degrees per nanometer. They need to be multiplied with a factor of (180/π) 2 to be stated in units conventionally used in other work (e.g. [49,73]) i.e., eV nm. than traditional periodic band diagrams for quasi-one-dimesional (nanotube-like) structures, which have only one index (i.e., the wave vector k z in the axial direction) labeling the Kohn-Sham states. 58 We anticipate that helical band diagrams like the ones shown are likely to emerge as powerful tools in understanding instabilities in optical and electronic materials. The helical band diagrams described above can be used to compute the  size of the band gap of the system and infer whether it is conducting, semiconducting or insulating. 59 We plot in Figure 9 the variation in the band gap of the above nanotube, as it undergoes torsional deformations. The variation with and without atomic 58 Similar band diagrams have been considered in [173,174] in the context of phonon calculations of carbon nanotubes. 59 For an extended system, such as a helical structure, the band gap is defined as the difference between the conduction band minimum (CBM) and the valence band maximum (VBM). Within Helical DFT, this can be computed as the difference between the smallest eigenvalue above the Fermi level and the largest eigenvalue below the Fermi level, as η is varied in I and ν in 0, 1, 2, . . . N − 1.
relaxation effects are both displayed. From the figure, we observe that the nanotube has a semiconducting behavior overall, and can be made to go from an insulating state at no twist (direct band gap of about 0.81 eV), to a practically conducting one, once the twist reaches about 5 degrees per nanometer. This is a rather significant change in the electronic properties of the material, although in a mechanical sense, its deviation from simple linear elastic behavior (Figure 7) is still fairly modest at this level of twist. To further illustrate the above electronic transition, we compute the (electronic) density of states of the nanotube without and with twist (β ≈ 5 degrees/nanometer). Following [44], we write this at a given energy level E and an electronic temperature of T e as: and evaluate it along a fine mesh of values of E in the range [−1, 1]. The results are shown in Figure 10. It is apparent from the figure that the electronic states in the system undergo significant change as the nanotube is subjected to twisting. In particular, the number of states at or near the Fermi level λ F , increases to values well above 0, indicating onset of metallic behavior when the tube is twisted.

Electronic properties: Behavior of chiral tubes under axial strains
Finally, we use Helical DFT to study a chiral phosphorene nanotube. We choose a tube with N = 64 and τ = 8.244 Bohr as before, and also set α = 0.005. We relax the positions of the atoms within the simulation cell and use the resulting nanotube 60 Although it is well known that LDA is often unable to predict quantitatively accurate band gaps, the qualitative trends observed here are likely to be representative of actual physical behavior [175][176][177][178][179] in these nanotube systems. In any case, the formulation presented here does not have issues with regard to the use of more quantitatively accurate hybrid exchange correlation functionals [180,181] whose implementation within our framework is a subject worthy of further investigation.  (observed to have an indirect band gap of about 0.24 eV) as the starting structure for subsequent simulations. We subject this tube to (both tensile and compressive) axial strains by varying τ , and relax the atomic positions in each case. Figure 11 shows the variation in the band gap of this tube at different values of axial strain. The effect of not including atomic relaxation after the tube is subjected to the strains is also shown. From the figure, we see that in the range of strains considered, compressive strains tend to reduce the band gap, while tensile strains appear to increase it. The tube appears to transition into a metallic state at about 4% compressive strain. Based on the nature of the plot and motivated by earlier work on chiral carbon nanotubes [182,183], we fitted a one-term Fourier series to this band gap data, and found that this produces a high quality fit for both the relaxed and unrelaxed cases. This suggests that like the case of carbon nanotubes, it might be possible to build (approximate) tight-binding type models of the band gap behavior [184] for phosphorene nanotubes. This is a topic that warrants further investigation. Notably, the parameters in such models can be provided through high-quality ab initio simulations based on Helical DFT. The above ab initio simulations of the phosphorene nanotubes (both zigzag and chiral) suggests that these materials have highly adjustable electronic states. Therefore, they might find future applications as nanomaterials with tunable electronic/optical/transport properties. The simulations also highlight the possibility of using inhomogeneous strain modes for altering these properties, instead of homogeneous strain modes which are considered in the strain-engineering literature [185][186][187] more commonly. Further investigations of this material and others, along these lines, is the scope of future work.
Finally, we find it worthwhile to point out that the phosphorene nanotube simulations described above (both electronic and mechanical) would be very challenging, or well-nigh impossible to carry out using conventional first principles techniques, even with the aid of massively parallel high-performance computing resources. In contrast, our MATLAB implementation of Helical DFT often allows such simulations to be carried out within a few hours of simulation wall time on a desktop workstation or a single node of a supercomputing cluster, thus serving to reinforce the novelty and practical utility of the approach.

Conclusions and Future Directions
In summary, we have presented Helical DFT -a novel, systematic first principles simulation framework for systems with helical symmetries. We have presented a derivation of the equations of Kohn-Sham theory, as they apply to the case of helical structures. Our derivation is systematic, self-contained and for the most part, mathematically rigorous. We have then solved these governing equations numerically by using a finite difference method based on helical coordinates. Using this working realization of the proposed approach, we have carried out simulations involving phosphorene nanotubes, extracted their mechanical behavior ab initio, and identified changes in the electronic properties of this material as it undergoes twisting.
Having laid out this foundational work on Helical DFT, we now discuss a number of avenues for future investigation: • Development of an efficient spectral solution scheme: While the current finite difference based implementation of Helical DFT allows us to investigate a number of materials systems of interest, it also suffers a number of computational deficiencies (Section 3). Following our earlier work on the ab initio simulations of cluster systems with arbitrary point group symmetries, we have already formulated, and are currently in the process of implementing a spectral scheme for solving the governing equations [8,46,188]. This is expected to completely resolve the issues with the finite difference formulation and provide a more efficient numerical implementation, thus opening the door to the simulation of more complex helical structures.
• Mechanistic simulations of quasi-one-dimensional systems: The materials science literature is rich with examples of quasi-one-dimensional structures that have been discovered and/or synthesized through experimental means. The simulation tools developed here can be used to characterize these materials computationally. In particular, the use of helical symmetry adapted ab initio molecular dynamics [189] can be used to study the mechanical behavior of these materials under axial and torsional strains, as well as their instabilities and defects [49,51,190].
• Helical Wannier states and the helical Berry phase: The helical Bloch Floquet transform U introduced in eqs. 38 and C.2 allows us to consider the notion of helical Wannier states. These, like their periodic counterparts, are defined through the action of U −1 on the helical Bloch states [85], and can be expected to be exponentially localized for insulating systems [191]. Additionally, like the case of periodic systems, a geometric phase (i.e., the Berry phase [192]) associated with the helical Bloch phase factor may be defined. These observations are likely to spur the development of novel computational analysis methods [193,194] for helical materials, as well as the discovery of novel topological materials [195,196].
• Search for exotic materials properties and study of multi-physics coupling: Since Helical DFT is a first principles simulation technique, it allows investigation of the effect of torsional deformations on a material's optical, electronic, magnetic and transport properties. Like the case of the phospherene nanotube considered in this work, it is possible that torsional or axial deformations in certain helical structures might induce a significant redistribution of electronic states in the material, leading to the appearance of interesting collective properties. In this regard, the investigation of the nanoscale flexoelectric effect [197][198][199][200][201][202] in helical systems would be of particular interest, since owing to the quasi-one-dimensional nature of these materials, as well as the appearance of strain gradients in connection with torsional deformations, a significant polarization may appear along the axis of a nanostructure when twisted, thus leading to a strong flexoelectric effect.
• Search for new phases of matter and coherent phase transformations: Finally, as remarked in [70,71], the first principles techniques developed here might be instrumental in the discovery of novel phases of matter. Additionally, the discovery of transformations between such phases [70,203] -particularly, coherent ones [204,205] -are likely to lead to new classes of active materials. The methods developed here are likely to be very useful in the characterization of energy barriers of such transformations and help in their design.
-Appendix A. Verifying that the operator h aux η in Eq. 25 is symmetric The boundary conditions associated with the operator h aux η introduced in Eq. 25 are quite non-conventional. Here, we work through the steps of verifying that the operator is symmetric when augmented with these boundary conditions. For this, we consider smooth functions f 1 , f 2 over D (the interior of the fundamental domain D G 1 ), obeying the same boundary conditions as φ in the proof of Theorem 2.3, i.e., f 1,2 (x) = f 1,2 (Υ h • x), ∇f 1,2 (x) = R −1 2πα ∇f 1,2 (Υ h • x) for x ∈ ∂D z=0 and f 1,2 (x) = 0 for x ∈ ∂D r=R . Then we have: On using integration by parts [92], the first term on the right-hand side becomes: where ds is the oriented surface measure. The boundary terms on the right-hand side of eq. A.2 can be split as: 3) The first term on the right-hand side of eq. A.3 goes to zero due to the Dirichlet boundary condition f 2 (x) = 0 for x ∈ ∂D r=R . To evaluate the second and the third terms, collectively denoted as I 1 henceforth, we write the surface measures in terms of the local unit normals to get: We now express y ∈ D z=τ as Υ h • x with x ∈ D z=0 , so that the first integral on the right-hand side of eq. A.4 above can be rewritten as (with the dependence on x and y shown explicitly): Here, the notation ds y is used to denote the surface measure centered at the point y, and similarly ds Υ h •x denotes the surface measure centered at the point Υ h • x. Since Υ h is an isometry, ds Υ h •x has the same magnitude as ds x . Furthermore, the boundary conditions imply that f 2 (Υ h • x) = f 2 (x), and: since R 2πα (and hence R T 2πα ) has axis e 3 . Combining the above results, we arrive at: and therefore, eq. A.2 reduces to: The second term on the right-hand side of eq. A.1 is: with ds denoting the oriented surface measure (as earlier). The surface integral can be split as: from which the first term on the right-hand side vanishes due to Dirichlet boundary conditions on ∂D r=R . For the second and the third terms, denoted collectively as I 2 henceforth, we write the oriented surface measures in terms of the local unit normals to get: As earlier, we express y ∈ D z=τ as Υ h • x with x ∈ D z=0 , so that the first integral on the right-hand side of eq. A.12 above can be rewritten by use of the boundary conditions as (the dependence on x and y has been shown explicitly): It follows that, and therefore, the second term on the right-hand side of eq. A.1 is: Combining eqs. A.1, A.8 and A.15, we get: On the other hand, we can express f 1 , h aux η f 2 L 2 (D) as: Integrating by parts the first term on the right-hand side of eq. A.17, we have: The second term on the right-hand side of eq. A.18 can be sent to zero by application of the boundary conditions (using a procedure similar to the one outlined in eqs. A.3 -A.7). This leaves us with: which implies that the operator h aux η is symmetric on smooth functions obeying the boundary conditions outlined above. Since such functions are dense in the domain of h aux η (the boundary conditions being interpreted in the trace sense in that case), the symmetry of the operator follows.
Vector addition and scalar multiplication are defined pointwise in this space, i.e., for s ∈ M , f 1 , f 2 ∈ H and z ∈ C, we have f 1 + f 2 (s) = f 1 (s) + f 2 (s) and z f (s) = z f (s), while the inner product is defined as: The holds. We then call A decomposable, we refer to the operators A(s) as the fibers of A, and we denote this relationship as: Conversely, it is also possible to "build" operators on the space H, starting from operators on the space H . Specifically, given a measurable function A(·) from M to the set of (bounded or unbounded) self-adjoint operators 63  As before, we will use the notation: to denote the above relationship between the operator A and its fibers A(s).
With these definitions in place, we now apply the above apparatus to carry out a suitable decomposition of the single-electron Hamiltonian associated with a helical structure. We will consider the case of a structure associated with a helical group G 1 that is generated by a single element Υ h . As discussed in 2.2, the single-electron Hamiltonian in this case is the operator H = − 1 2 ∆ + V (x) over the space H = L 2 (C). The relevant boundary condition is ψ(x) = 0 for x ∈ ∂C and the potential V (x) is group invariant, i.e., V (x) = V ( Υ • x) for every Υ ∈ G 1 .
We may choose M = [− 1 2 , 1 2 ) = I, H = L 2 (D), µ as the Lebesgue measure on I (denoted as dη henceforth). Then, is the space L 2 (D × I). The helical Bloch-Floquet transform U introduced in eq. 38 allows us to map functions and operators in between the spaces H = L 2 (C) and H = L 2 (D × I).
Specifically, for each η ∈ I, let H η be the restriction of H to D along with the boundary ∇ψ(x) for x ∈ ∂D z=0 and, ψ(x) = 0 for x ∈ ∂D r=R . Then, the operator H can be decomposed into the set of operators H η η∈I (referred to as the fibers of H) using the direct integral decomposition, in the sense: To demonstrate this result, we will establish direct integral decompositions of the Laplacian and the potential operator individually as operators on L 2 (C), i.e., we will show successively that: and, The sought result will then follow from Theorem XIII.85g of [88].
To perform the direct integral decomposition of the Laplacian on L 2 (C), subject to zero Dirichlet boundary condition for x ∈ C, we first define the fibers (−∆) η as the restriction of −∆ to D along with the boundary conditions p(Υ h • x) = e −i2πη p(x), ∇p(x) for x ∈ ∂D z=0 and, p(x) = 0 for x ∈ ∂D r=R . Now, let B be the operator on the right-hand side of eq. B.9 acting on the space H = L 2 (D × I), and let f be a Schwartz class function that obeys f (x) = 0 for x ∈ C. Note that since f is a Schwartz class function, so are all its derivatives. With this setup, it suffices to show that Uf ∈ Dom.(B), and that U(−∆f ) = B(Uf ), in order to establish eq. B.9 (see e.g. Theorem XIII.87 and the lemma following that theorem in [88]).
By definition, Uf is given as: 11) and due to the properties of f , it is a smooth function in x over D that obeys (Uf )(x, η) = 0 for x ∈ ∂D r=R for any η ∈ I. Furthermore, by evaluating the above expression at Υ h •x, it is clear that (also see Footnote 30),the relationship: holds. Computing the gradient on both sides, we get: Now, evaluating the above expressions at x ∈ ∂D z=0 , we see that (Uf )(x, η) obeys all the boundary conditions necessary for it to be in Dom.(B). Next, we recall that the Laplacian is invariant under isometries 64 . Consequently, it holds that: (B.14) and so, U(−∆f ) = B(Uf ), as required. This establishes 65 eq. B.9. Next, to establish B.10, we set the fibered operator V η on L 2 (D) as: Then for any Schwartz class function f over C, we have that: which, using the invariance of the potential under G 1 , gives: Thus B.10 is established. With B.8 established, some conclusions can be immediately drawn regarding the structure of the spectrum of H. Using Theorem XIII.85 in [88] for example, it follows that the collection of eigenvalues of the operators H η together form the spectrum of H, i.e., Λ = spec.(H). Furthermore, Theorem XIII.86 of [88] implies that H has a purely continuous spectrum.
Finally, a result made use of in Section 2.2.2 is that if A is an operator that is invariant under G 1 , i.e., it commutes with the unitary operators in the set: and it is is locally trace-class, we may assign meaning to the trace per unit fundamental domain as: Discussion and rigorous proofs of this result appear in [206,207] for the case of periodic symmetries. To see why this result also holds true for the case of a helical symmetry group, we denote I D as the indicator function of the fundamental domain. Let {f j : j ∈ N} be an orthonormal basis of L 2 (C). Due to the fact that A is invariant under the group G 1 , we may write: Now, the trace per unit fundamental domain is: Using eq. B.20, this becomes: which, upon using B.5 and B.2, can be written as: Now, recognizing that the indicator function I D acts as a projection operator on L 2 (C), we see that U I D f j are simply basis functions of L 2 (D) for every η ∈ I. Denotingf j (x; η) = (U I D f j )(x; η) and exchanging the summation and the integral, the above equation can be re-written as: as claimed.
Appendix C. Helical structure associated with a group generated by two elements: Expressions for important physical quantities and governing equations For a helical structure associated with a group G 2 that is generated by two group elements (i.e., a screw transformation of the form (R 2πmα |mτ e 3 ) and a pure rotation of the form (R nΘ |0) with Θ = 2π N , and both rotations R 2πmα , R nΘ having e 3 as the common axis) the form of the governing equations as well as the expressions for the various quantities of interest can be easily deduced by using the discussion in Section 2.2 as a starting point, and using the following rules of substitution. The helical Bloch phase factors are of the form e −i2π(mη+ nν N ) with n, ν ∈ {0, 1, . . . , N − 1} and η ∈ I, instead of being e −i2πmη only. Consequently, helical Bloch states and the helical bands are labeled as ψ j (x; η, ν) and λ j (η, ν), respectively. Integrals over η ∈ I = [− 1 2 , 1 2 ) need to 57 be replaced with integrals over η ∈ I and summations over ν of the form . Taking into account the effect of all the symmetry operations of the group (i.e., computing the group orbit of a point for example) amounts to summing over all group elements of the form Υ m h • Υ n c = (R 2πmα+nΘ |mτ e 3 ) with m ∈ Z and n ∈ {0, 1, . . . , N − 1}. Also, in the expression for the forces (eq. 118), the rotation matrix (R 2πmα ) −1 needs to be replaced with (R 2πmα+nΘ ) −1 . Finally, spatial integrals over the fundamental domain D G 1 (or its interior D) have to be replaced by integrals over the fundamental domain D G 2 (or its interior D). The reciprocal space (or more specifically, the Brillouin zone of the reciprocal space) for the symmetry group G 2 will be denoted by the set B = I × {0, 1, 2, . . . , N − 1}.
The mathematical reasons behind the above substitution rules are as follows: The structure of the helical group generated by two elements is such that it can be expressed as the direct product of the helical group generated by a single element and a cyclic group. Therefore, the Bloch theorem for a structure associated with such a symmetry group can be established by first using Theorem 2.3, and then using an appropriate version 66 of Theorem 2.6 from [40]. The characters (i.e., one-dimensional complex irreducible representations) of the helical group generated by two elements can be obtained by multiplying out the characters of the helical group generated by a single element and the cyclic group. This leads to the helical Bloch states in this case to obey the condition: The completeness of these states follows by a combination of Theorem 2.6 and the properties of the Peter-Weyl projectors [40][41][42] for the cyclic group. This completeness result paves the way for a suitable Helical Bloch-Floquet transform U : L 2 (C) → L 2 ( D × B): (Uf )(x, η, ν) = as well as a direct integral representation of the single electron Hamiltonian. From these, the key quantities of interest and the governing equations may be systematically deduced (as demonstrated for the case of the helical group generated by a single element in Section 2.2). Keeping the above discussion in mind, we now present the final mathematical expressions for key quantities of interest for a helical structure generated by G 2 . For the sake of definiteness, we assume that the helical structure, is embedded in the infinite cylinder C (eq. 10), and as in eq. 11, we let D G 2 = (r, ϑ, z) : 0 ≤ r < R, 0 ≤ ϑ < 2π N , 0 ≤ z < τ denote a cylindrical sector that acts as the fundamental domain of the helical structure. Let D denote the interior of this fundamental domain. Furthermore, let the positions of the atoms within the fundamental domain be P G 2 = x k M G 2 k=1 and let b k (x; x k ) M G 2 k=1 denote the corresponding nuclear pseudocharges. The helical Bloch states presented above in eq. C.1 allow the single electron problem for the entire helical structure to be reduced to the fundamental domain D G 2 (or its interior D).
Let g j (η, ν) = f Te λ j (η, ν) denote the electronic occupation numbers at electronic temperature T e . The electron density for x ∈ D can be expressed as: g j (η, ν) |ψ j (x; η, ν)| 2 dη . (C. 3) The electron density needs to obey the constraint of integrating to a fixed number of electrons in the fundamental domain, i.e., g j (η, ν) ∆ψ j (·; η, ν), ψ j (·; η, ν) L 2 ( D) dη . (C.7) The second term is the exchange correlation energy per unit fundamental domain and is expressible as: g j (η, ν) χ k,p (·; η, ν; x k ), ψ j (·; η, ν) L 2 ( D) 2 dη , (C.9) with the functionsχ k,p (x; η, ν; x k ) being related to the helical Bloch Floquet transform associated with the group G 2 (eq. C.2) i.e.,: Here, as in eq. 87, χ k;p (·; p k ) denote the atom centered nonlocal pseudopotential operator projection functions. The fourth term on the right-hand side of eq. C.6 is the net electrostatic interaction energy per unit fundamental domain, and can be written as: Here, b(x, , P G 2 , G 2 ) is the net nuclear pseudocharge: and E sc (P G 2 , G 2 , D) represents self interaction and correction terms. The net electrostatic potential Φ can be expressed as the Newtonian potential: 67 The last term on the right-hand side of eq. C.6 is the electronic entropy contribution: g j (η, ν) log g j (η, ν) + 1 − g j (η, ν) log 1 − g j (η, ν) dη (C.14) The symmetry adapted Kohn-Sham equations for the system, as posed over the fundamental domain are:

H KS
[Λ,Ψ,P G 2 , D,G 2 ] ψ j (·; η, ν) = λ j (η, ν) ψ j (·; η, ν) , is the Kohn-Sham operator, with its dependence on the helical Bloch sates, the helical bands, etc., made explicit 68 . As before, V xc represents the exchange correlation potential, while the nonlocal pseudopotential operator V nl D has fibers (in coordinate representation): V nl D (x, y; η, ν) = M G 2 k=1 p∈N k γ k,pχk,p (x; η, ν; x k ) ⊗χ k,p (y; η, ν; x k ) (C.18) 67 As in the main text (Section 2.2.2), we have, with a minor abuse of notation, used Φ to denote the "trial" electrostatic potentials involved in the maximization problem in eq. C.11 as well as the actual potential that achieves this maximum, i.e., eq. C.13. 68 As in the discussion in Section 2.2.2, we may change notation and express the free energy per unit cell in terms of the occupations as F(G, Ψ, P G2 , D, G 2 ) = F(Λ, Ψ, P G2 , D, G 2 ), instead of using the eigenvalues. With this change in notation, the problem of determination of the Kohn-Sham ground state for a given helical structure, can be expressed as: The helical Bloch states obey (for i, j ∈ N,): ψ i (·; η, ν), ψ j (·; η, ν) L 2 ( D) = δ i,j .

(C.22)
A more comprehensive account of the above equations as well as extensive numerical simulations involving them appear in a follow-up contribution [73].

Appendix D. Cartesian gradient operator, Laplacian operator and integrals in helical coordinates
We derive expressions in helical coordinates for the Cartesian gradient operator (useful in expressing forces on the atoms), the Laplacian operator (useful for expressing the Schrodinger operator), and integrals (useful for computing energies) in this Appendix. The calculations involved are straight forward applications of the chain rule and they originally appear in [8]. We include these here for the sake of completeness.