Construction of exact minimal parking garages: nonlinear helical motifs in optimally packed lamellar structures

Minimal surfaces arise as energy minimizers for fluid membranes and are thus found in a variety of biological systems. The tight lamellar structures of the endoplasmic reticulum and plant thylakoids are comprised of such minimal surfaces in which right and left handed helical motifs are embedded in stoichiometry suggesting global pitch balance. So far, the analytical treatment of helical motifs in minimal surfaces was limited to the small-slope approximation where motifs are represented by the graph of harmonic functions. However, in most biologically and physically relevant regimes the inter-motif separation is comparable with its pitch, and thus this approximation fails. Here, we present a recipe for constructing exact minimal surfaces with an arbitrary distribution of helical motifs, showing that any harmonic graph can be deformed into a minimal surface by exploiting lateral displacements only. We analyze in detail pairs of motifs of the similar and of opposite handedness and also an infinite chain of identical motifs with similar or alternating handedness. Last, we study the second variation of the area functional for collections of helical motifs with asymptotic helicoidal structure and show that in this subclass of minimal surfaces stability requires that the collection of motifs is pitch balanced.

Minimal surfaces arise as energy minimizers for fluid membranes and are thus found in a variety of biological systems. The tight lamellar structures of the endoplasmic reticulum and plant thylakoids are comprised of such minimal surfaces in which right and left handed helical motifs are embedded in stoichiometry suggesting global pitch balance. So far, the analytical treatment of helical motifs in minimal surfaces was limited to the small-slope approximation where motifs are represented by the graph of harmonic functions. However, in most biologically and physically relevant regimes the intermotif separation is comparable with its pitch, and thus this approximation fails. Here, we present a recipe for constructing exact minimal surfaces with an arbitrary distribution of helical motifs, showing that any harmonic graph can be deformed into a minimal surface by exploiting lateral displacements only. We analyze in detail pairs of motifs of the similar and of opposite handedness and also an infinite chain of identical motifs with similar or alternating handedness. Last, we study the second variation of the area functional for collections of helical motifs with asymptotic helicoidal structure and show that in this subclass of minimal surfaces stability requires that the collection of motifs is pitch balanced.

Introduction
From soap films and liquid crystals to plant thylakoids and the endoplasmic reticulum, minimal surfaces arise as the ground state of a variety of manmade and naturally occurring membranes and lamellar structures. Minimal surfaces minimize the area of a surface that passes through a given boundary [1], and thus arise naturally in cases where surface tension is dominant. However, minimal surfaces also arise in systems dominated by membrane bending energy as they constitute trivial critical points of the Helfrich free energy [2].
Recently, helical motifs were discovered in the lamellar structures of the endoplasmic reticulum [3] and plant and cyanobacteria thylakoids [4,5]. Both right and left handed motifs were observed in both systems, in stoichiometry that suggested global pitch balance. While in the endoplasmic reticulum the right and left handed helical motifs were mirror images of one another and appeared in equal amounts, in the thylakoid the right and left handed motifs differed in their pitch, core radius, and density yet preserved pitch balance on average [6].
To advance our understanding of these systems we require the ability to construct minimal surfaces in which the appropriate helical motifs are embedded. Twist grain boundaries, comprised of infinitely many helical motifs of a single handedness embedded along a line in a minimal surface, allow an exact formulation and are well understood [7][8][9]. In contrast, embedding finitely many motifs in a minimal surface, as well as combining motifs of different pitch values remains a challenge. Coming to address these structures more explicitly Guven et al [10] considered the small slope approximation for minimal surfaces. The main advantage of this approach is that minimal surfaces are obtained as harmonic graphs (i.e. surfaces given in a Monge-patch parameterization (x, y, h(x, y)), where h is harmonic), and thus this approach allows for the simple addition of helical motifs of arbitrary geometry and topology. However, the resulting surfaces are not exactly minimal. This is exceptionally apparent in the immediate vicinity of the helical structures, i.e. within a distance of a few pitch lengths from the core of each motif, as can be observed in Fig. 1.
In the biologically and physically relevant regimes, helical motifs are commonly separated by a distance comparable to the pitch, rendering the small slope approximation irrelevant for these structures [3,6]. With the purpose of bridging this gap, in this work, we provide a recipe for constructing exact minimal surfaces with an arbitrary distribution of helical motifs. We first show that any harmonic graph (approximately minimal surface) can be deformed into an exact minimal surface through an explicit, yet non-local operation. We moreover show that every minimal surface could be obtained through this recipe. We conclude by surveying key examples: a pair of motifs of the same and opposite pitch, and an infinite chain of identical motifs distributed along a line with similar or alternating handedness.

Enneper immersions of minimal surfaces
We begin by presenting a somewhat underused representation of minimal surfaces termed Enneper immersions due to Andrade [11]. In essence, this technique allows us to "fix" an approximate minimal surface given as a graph over the plane such that it becomes an exact minimal surface by only exploiting lateral displacements and not varying the height data. Consider a surface given by a harmonic function h over a domain in the complex plane r(z) = (z, h(z)), z = x + iy ∈ Ω ⊂ C, ∂z∂zh = 0, (2.1) where for convenience we will from now on represent surfaces though complex (non necessarily analytic) functions, and also usez = x − iy, ∂z ≡ ∂ ∂z = 1 2 (∂x − i∂y), and ∂z ≡ ∂ ∂z = 1 2 (∂x + i∂y). The mean curvature H of such a surface reads H = (1 + h 2 x )hyy − 2hxhyhxy + (1 + h 2 y )hxx 2(1 + h 2 x + h 2 y ) 3 2 = −2 h 2 z hzz − (1 + 2hzhz)hzz + h 2 z hzz (1 + 4hzhz) 3 2 . As , such surfaces are considered minimal within the small slope approximation (where terms of order O( ∇h 2 ) are omitted). A great advantage of this method is that due to the additivity of harmonic functions one can simply "add" helical motifs. The surface yields a helicoid of pitch p 0 centered around z 0 which is exactly minimal. One can explicitly construct the surface which includes N helical motifs located at the points {z k } N k=1 and of corresponding pitches The mean curvature of this surface, however, does not vanish identically. (The helicoids and planes are the only minimal surfaces that are also the graph of a harmonic function.) In particular, considering the case of a dipole of pitches +p and −p separated a distance R from each other, say with axes located at z 1 = R 2 and z 2 = − R 2 , the mean curvature reads For example, as x = R/2 and y → 0 then H → −8/R, or expressed in dimensionless terms |H/ √ −K| = 2p/R. We thus see that the deviation from minimality becomes significant whenever the inter-motif separation becomes comparable to their pitch (which is the case in both [3] and [6]). To amend this we turn to construct a specific type of an Enneper immersion. We begin by stating the general result due to Andrade [11]: Let h : Ω ⊂ C → R be a harmonic function and L, P : Ω ⊂ C → C be holomorphic. We may define The mapping X is a conformal minimal surface when L and P satisfy L (z)P (z) = ∂h ∂z 2 and |L (z)| + |P (z)| = 0. (2.5) The first condition guarantees that X is conformal, i.e., a 11 = a 22 and a 12 = 0, while the second implies that the surface is regular, i.e., det a ij = 0. Conversely, every minimal surface can be written as the Enneper graph of a harmonic function. The reader is referred to [11]

(a) Obtaining exact minimal surfaces from harmonic approximations
Consider an approximately minimal surface in the form of a harmonic graph (2.1), then the surface is an exact minimal surface. Conversely, every minimal surface can be locally parameterized as in Eq. (2.6). Moreover, the complex variable z provides a conformal parameterization of the surface r, and the area element dA and Gaussian curvature K read: and It is important to state that the x, y coordinates above are not Cartesian coordinates on the plane, i.e., they do not provide a Monge patch for the surface r, rather they are the image of these coordinates under the transformation z → z − P (z) 1 .
We next provide a rigorous proof of these claims.
In order to find P for a given minimal graph r(ζ) = (ζ, h(ζ)), we must be able to write a differential equation for P in terms of ζ andζ. To do that, we use Pz = P ζ ζz + Pζζz = P ζ − PzPζ 1 Distinguishing the coordinate systems is exceptionally important in view of the functional similarity between the above result and the expression for the Gaussian curvature of an arbitrary surface given as a graph: K (x,y,h(x,y)) = (hxxhyy − h and 0 = Pz = P ζ ζz + Pζζz = Pζ −PzP ζ to find the system of differential equations Notice, B is the "plus" solution for Pz and, consequently, the denominator of the ODE's in Eq. (2.9) is where we used that h ζ hζ = 1 4 ∇h 2 , which implies that the numerator of |B| 2 − 1 is a sum of positive quantities.
In short, the equations in (2.9) provide necessary conditions for the existence of the coordinate change ζ → ζ(z). It remains to show they are also sufficient. From now on, let r(ζ) = (ζ, h(ζ)) be a minimal graph. Thus, it makes sense to write the system of ODE's in Eq. (2.9), whose solvability condition is Pζ ζ = P ζζ . Using the definition of B, we have where H is the mean curvature of r(ζ) = (ζ, h(ζ)) as given in Eq. (2.2). In conclusion, we can integrate Eq. (2.9) if and only if r(ζ) = (ζ, h(ζ)) is minimal. Thus, after integrating Eq. (2.9) we can find P (ζ) and define z as z = ζ + P (ζ). Since the Jacobian of the transformation ζ → z(ζ) is we can invert ζ → z(ζ) and write ζ(z) = z − P (z).
Using the ODE's for P ζ and Pζ , we can solve this system and find that Pz = B(ζ(z)) and Pz = 0, which in particular implies P (z) must be holomorphic. In addition, using that Pz = B and that hz = h ζ ζz + hζζz = h ζ −Pzhζ , we can straightforwardly deduce that Pz = (hz) 2 .
In conclusion, given a minimal surface Σ parameterized as a graph r(ζ) = (ζ, h(ζ)), we can perform a change of coordinates ζ → z(ζ) and reparameterize Σ as an Enneper immersion in the particular form r(z) = (z − P (z), h(z)), i.e., we are able to see Σ as the deformation of a harmonic graph.

(b) Basic example: helicoid
Since a single helicoid is exactly minimal, one may expect that applying Eq. (2.6) to Eq. (2.3) would result in a trivial map, i.e., it would give r = (z, p 0 arg(z − z 0 )) back. However, an Enneper immersion is necessarily conformal, which is not the case for r = (z, p 0 arg(z − z 0 )). The above intuition is, however, not entirely wrong; the Enneper graph of h(z) = p 0 arg(z − z 0 ) yields the same helicoid but parameterized via a conformal immersion.
To be more precise, the parameterization is a double covering of the helicoid since ρ ∈ (−∞, 0] for 2r ≤ |p 0 | and ρ ∈ [0, ∞) for 2r ≥ |p 0 |. Notice that the entire boundary of the disc r = |p 0 |/2 is singular and can be associated with the axis of the helicoid. The large distortion associated with the mapping and the singularities is a consequence of the conformality of the parameterization.
The two copies of the helicoid can be easily distinguished using the surface's unit normal, which can be expressed as N = (p 2 0 + ρ 2 ) − 1 2 (p 0 sin θ, p 0 cos θ, ρ) ∈ S 2 . It is immediate to see that in the exterior (interior) of the disc of radius (2.11) These results will allow us to compare the surfaces obtained for interacting motifs with that of a single helicoid.

(c) Finite collections of helical motifs and the multipole expansion
The construction of an Enneper graph is associated with a harmonic function, whose singularities may be interpreted as charges in analogy to point charges in electrostatic. This analogy suggests that we can resort to a multipole expansion by considering helical motifs as point charges. While in the 2d electrostatic the potential associated with a point charge has the logarithmic singularity V (r, θ) = V (r) = ln 1 r , in the context of helical motifs point charges should have an arg-singularity h(r, θ) = h(θ) = θ. (Up to constants, h = θ is the only purely angular 2d harmonic function while h = ln r is the only purely radial 2d harmonic function.) Given an electric charge density µ(x, y) entirely contained in a disc D R = {(x, y) : x 2 + y 2 ≤ R}, the electrostatic potential outside the disc is harmonic ∆V = 0 in R 2 \D R , but inside D R we have ∆V = 2πµ. (We are employing a distinct sign convention that will prove to be more useful when extending the approach to our context.) Outside the disc D R , V has the multipole expansion [12] where c k = a k + ib k ∈ C and p ∈ R. The coefficients p, a k , and b k are written as functions of the charge density µ according to Now, considering the imaginary part, the multipole expansion associated with a "helicoidal" charge density µ is where the coefficients are the same as in Eq. (2.13). For a single helicoid of pitch p 0 located at 0 ∈ C, we know that the corresponding harmonic function is h(r, θ) = p 0 θ = p 0 Im ln z. Thus, a single helicoidal charge p 0 is, as expected, associated with the distribution µ = p 0 δ(z), where δ is the Dirac delta function. (This distribution leads to c k = 0.) On the other hand, for a helicoidal charge p 0 located at z 0 ∈ C, we have This expression is in agreement with Eq. (2.14) for the distribution µ = p 0 δ(z − z 0 ). Indeed, In general, for a set of helicoidal charges {p j } located at {z j } N j=1 ⊂ C, the charge density is given by µ = N j=1 p j δ(z − z j ). Then, the coefficients in the multipole expansion are Finally, the corresponding harmonic function is Notice that this power series converges absolutely as long as r > r k for all k, i.e., for all z outside the disc of radius max k {r k }. This expansion shows that the far-field behavior of the fundamental layer of any finite collection of helical motifs asymptotes to that of a helicoid whose pitch is given by the sum of the pitches of the individual motifs. Several such intertwined simple helicoids may be required to fully cover all leaves of the structure. Note, however, that an infinite collection of helical motifs, as well as motifs that are not perpendicular to the parallel layers they pierce, are not described by this expansion.

Concatenating helical motifs
Let us investigate in some detail the concatenation of two helical motifs, see Fig. 2. Consider two helical motifs of pitches p 1 and p 2 located at z 1 and z 2 , separated a distance R from each other. The corresponding harmonic function is given by h(z) = p 1 arg(z − z 1 ) + p 2 arg(z − z 2 ), whose graph is only minimal up to leading order in ∇h . Without loss of generality, we may assume (3.1) The first two terms above are the sum of the functions P (z; {z 1 , p 1 }) and P (z; {z 2 , p 2 }) associated with helicoids around z 1 = R 2 and z 2 = − R 2 . The surface obtained as the algebraic sum of the Enneper immersions of h(z; is not a minimal surface. The last term in P (z), which can be written as p1p2 R times the scaleindependent function 1 2 ln[(ζ + 1 2 )/(ζ − 1 2 )], ζ = z/R, thus introduces a pair-interaction between the two helicoids required to make the corresponding parameterization a conformal minimal immersion.
In addition, the stereographic projection of the unit normal of a minimal 2-helicoidal parking garage is given by (see Appendix A for the relation between the normal N and the Enneper data) .

(3.2)
As the example of a single helicoid teaches us, we should pay attention in selecting the correct domain for a minimal 2-helicoid. The analysis of g will help us identify the proper domain for the Enneper graph of a parking garage and we shall consider as the domain of the Enneper graph the points z ∈ C where the unit normal N takes values on the North hemisphere: Ω N = {z ∈ C : |g(z)| ≥ 1}. The curves |g| = 1, i.e., the boundary of Ω N , (Fig. 3) will be associated with the helices of each helical motif, as explained in the following subsections.

(a) Concatenating two identical helical motifs
For simplicity, we start by considering the symmetric problem of gluing two identical motifs such that p 1 = p 2 = p. For large distances (as we prove in Subsection 2.(c)) the resulting surface is well approximated by two intertwined and identical copies of helicoids of pitch 2p displaced a distance p along the helicoids' axis with respect to each other, Fig. 2 (Right). The near field structure is, however, non-trivial. The Gauss map in this case reads , from which follows that for every w ∈ C, the equation g(z) = w will have in general two distinct solutions. This means that the unit sphere is covered twice by the unit normal and, therefore, the total Gaussian curvature should be −8π. However, selecting as the domain of definition the points z where the unit normal takes values on the North hemisphere, Ω N , the correct value of the total curvature of a minimal nondipole is −4π, i.e., −2π for each helical motif. The equator of the unit sphere is sent by the stereographic projection of the unit normal on the curve z → |g(z)| = 1. Notice that where C(ρ, θ) = ρ 4 − 1 2 ρ 2 cos 2θ + 1 16 is a function whose level sets are the well known Cassini ovals (see Fig. 3, Left plot). It follows that the level curves of |g| 2 are modified Cassini ovals that can be analytically described in terms of fourth-degree polynomial curves whose shape is depicted in Fig. 3 (Right plot). For small pitches, the region C − Ω N is made of two disconnected pieces around z 1 and z 2 . However, when the pitches increase, the level sets display a topological transition at pc ≡ R and the complement of the domain Ω N is no longer disconnected. The curves |g| = 1, whose shape depends on p/R, parameterize the axes of the two helical motifs. By increasing p/R, we notice that under the deformation in Eq. (2.6) the two axes move toward each other and, consequently, the distance between them is smaller than the initial separation R of the undistorted surface, see Fig. 4 (Center). For p < pc, we may compute the points of closest approach of the two axes from the points of closest approach on the distinct connected components of |g(z)| = 1. For the critical value p = pc the two axes finally intersect. In addition to getting closer along the line containing the motifs, the two axes also incline away from the line connecting the motifs but in opposite senses, see Fig. 4 (Left). Finally, the effective pitch remains unchanged during the deformation process. (The height data provided by the harmonic function h is preserved.) To describe the geometry of the helical core in the symmetric case, we may compute the gradient and Hessian of h(z) = p arg(z + R 2 ) + p arg(z − R 2 ). The gradient and Hessian of h are and (3.5) Therefore, from Eq. (2.8), the Gaussian curvature of a minimal pair of identical motifs is The infinitesimal area element of a minimal pair of identical motifs is It follows that the area of an annulus with large enough radii r 2 > r 1 R is given by r2 r1 dA ≈ π(r 2 2 − r 2 1 ) + 4πp 2 ln for each helical motif. In addition, noticing that the curves defined by the level sets of |g| 2 are the well known Cassini ovals. Thus, the shape of the curve associated with the equator of S 2 , i.e., z → |g(z)| = 1, is completely determined by the parameterp = |p|/R, as described in Fig. 3 (Left plot). Here the curves |g| = 1, whose shape depends on p/R, parameterize the axes of the two helical motifs. Under the deformation in Eq. (2.6) the central part of two axes of a dipole pair do not converge toward each other upon increasing p/R, but their extremities do. For p < pc ≡ R 2 , the points of closest approach of the two axes can be computed from the points of closest approach on the distinct connected components of |g(z)| = 1, see Fig. 5 (Center plot). The level curves of |g(z)| have a topological transition at pc = R 2 and it follows that for p > pc the two axes of the Enneper graph of two oppositely handed helical motifs of pitch p and −p merge into a single smooth curve. In addition, the two axes effectively incline in the direction orthogonal to the line containing the motifs (both in the same sense), see Fig. 5 (Left plot). Finally, the pitch remains unchanged during the deformation process.
As in the previous case, to describe the geometry of the helical core in the dipole case, we may compute the gradient and Hessian of h(z) = p arg(z + R 2 ) − p arg(z − R 2 ). The gradient and Hessian of h are and Therefore, from Eq. (2.8), the Gaussian curvature of a minimal dipole is Notice that K decays to zero faster than the Gaussian curvature of a minimal non-dipole. On the other hand, the infinitesimal area of a minimal dipole is It follows that the area of an annulus with large enough radii r 2 > r 1 R is given by (3.14) The non-trivial contribution is nothing but the area of a planar annulus. This is compatible with the fact that for large distances, a minimal pair of opposite helical motifs is well approximated by a single helicoid of pitch p 0 = p − p = 0, i.e., a minimal dipole is approximately a plane at large distances. (See the multipole expansion developed in Subsection 2.(c).)

(c) Gluing finitely many helical motifs
For an arbitrary number of helicoids of pitches p 1 , . . . , p N located at z 1 , . . . , z N , we may consider the harmonic function Now, since P = (∂zh) 2 , we have after integration Each term in the first sum above corresponds to the function P (z; {z k , p k }) associated with a single helicoid of pitch p k located at z k . The remaining terms can be seen as pair interactions between distinct helicoids required to assure that the corresponding parameterization is a conformal minimal immersion. The stereographic projection of the layer normal N is given by The coefficients a k and d k are computed from {p j } and {z j } as where s k denotes the k-th symmetric polynomial, s 0 (x 1 , . . . , xm) = 1 and s k (x 1 , . . . , xm) = 1≤i1<···<i k ≤m x i1 · · · x i k .

(d) Twist grain boundary -linear chain of infinitely many helical motifs
Twist grain boundaries (TGB) are structures comprised of an infinite collection of identical helical motifs evenly spaced along a straight line. These structures have been studied extensively in the context of smectic liquid crystals [13], where the individual helical motifs are termed screw dislocations after their crystalline cousins. Unlike individual motifs of finite charge which asymptote to parallel planes, TGB's asymptote to two different families of parallel planes on the two sides of the TGB line, Fig. 6 (Left). The two families of planes form a finite angle with respect to one another, whose magnitude is determined by the ratio of inter-motif distance and its pitch. The distinct layers in a smectic structure are often modeled within the small slope approximation as level sets of the phase field [8] φ Equivalently, individual layers may be identified with the graph of the harmonic function Im ln sin πz d .
The resulting surface is minimal only to order O( ∇h 2 ). Thus this approximation was primarily employed to address small angle variations in the surfaces normals, and TGB's comprised of helical motifs that are well separated compared with their pitch. An exact minimal TGB surface allows estimating the elastic energy in cases where the inter-motif distance is comparable to the pitch.
To achieve an exact formulation in the most transparent form, we start from considering the Enneper immersion of a linear stack of infinitely many helical motifs of pitch p, evenly spaced a distance d from each other: Differentiating we obtain where we used the identity π cot(πz) = z −1 + 2z n≥1 (z 2 − n 2 ) −1 . Squaring and integrating we obtain the correction term to produce the desired Enneper immersion Alternatively, we can rewrite P (z) as an infinite sum of helicoids of pitch p located at {n d : n ∈ Z} plus a linear correction: The stereographic projection of the layer unit normal is As in the concatenation of finitely many helical motifs, here the level sets of |g| help to find the proper domain to compute the Enneper graph of a minimal TGB (see Fig. 7, left). In addition, the image of the surface unit normal N over the fundamental domain R : |g(x + iy)| ≥ 1} covers the unit sphere exactly twice, which results in a total curvature of Since there are infinitely many helical motifs periodically distributed on a line, the total curvature of a minimal TBG is infinite. The Gaussian curvature K is computed from the first and second , from which follows that Finally, upon application of the deformation in Eq. (2.6) leading to the Enneper graph of a minimal TGB, the effective distance between neighboring axes diminishes under the increase of p/ d , as depicted in Fig. 8 (Left). It then follows that in order to obtain a minimal TGB with prescribed pitch and axes distance, the value of d in the expression of h(z) has to be appropriately tuned. The pitch remains unchanged during the deformation process.
Similarly, we may also consider a Untwisted Grain Boundary (UtGB), i.e., an infinite chain of evenly spaced helical motifs of alternating handedness, see Fig. 6 right. The corresponding harmonic function is then given by where p and −p are the pitches of the helical motifs in each dipole pair that are separated by a distance d . From we can compute the auxiliary function P (z) = (∂zh) 2 dz as . Upon application of the deformation in Eq. (2.6) leading to the Enneper graph of a minimal UtGB, any pair of neighboring axes deformed under the increase of p/ d , as depicted in Fig. 8 (Right). On average, the distance between neighboring axes remains constant. However, the top and bottom portions of the two axes in the same unit cell approach while for two neighboring axes from distinct cells this happens around the center of axes. For a fixed d , neighboring axes finally touch at the critical value pc = 2 d π , which corresponds to the topological transition observed in the level curves of |g(z)| as depicted in Fig. 7. The pitch remains unchanged during the deformation process.
, the Gaussian curvature is Finally, the stereographic projection of the unit normal of the layers of a minimal UtGB is The level sets of |g| help to find the proper domain to compute the Enneper graph of a minimal UtGB (see Fig. 7, right). In addition, the image of the surface normal over the fundamental domain

Second variation of the area functional
The above procedures allow us to construct explicitly exact minimal surfaces with any desired arrangement of helical motifs. The physical motivation for this construction is that minimal surfaces arise naturally as critical points of both the area of a surface and of the Helfrich free energy. However, while minimal surfaces are stationary points of these functionals, they are not necessarily (local) minima. To establish that a given configuration is indeed of minimal energy the second variation should also be examined. For simplicity, we restrict ourselves to analyzing the area functional alone. Moreover, the configurations we will consider will be restricted to harmonic functions of the form of Eq. (3.15), i.e., arrangements of helical motifs that asymptote to helicoidal surfaces. This is not the most general case. Nonetheless, as will be next demonstrated, this restricted class singles out global pitch balance as a necessary condition for stability. While claims of instability are valid only within this restrictive set of asymptotically helicoidal surfaces, surfaces proved to be locally stable will remain so even when considering all possible harmonic functions.
Considering normal variations of a minimal surface r : , the second derivative of the area functional is [14] Therefore, a minimal surface is said to be stable, i.e., it is a local minimum, if I(v) > 0 for all v ∈ H(D). (WhenD = D ∪ ∂D is not compact, the surface is said to be stable if it is stable for all Ω ⊂ D such thatΩ is compact.) We note that our variations are assumed to vanish on the boundaries (which are also assumed stationary). Thus, while negative values in the second variation imply instability, a positive definite second variation does not necessarily imply stability if the helical motifs are allowed to change their relative orientation or move in space. A stability criterion for fixed boundaries can be established based on the area over the unit sphere covered by the surface unit normal N. If the area of this spherical image is smaller than 2π, then the corresponding minimal surface has to be stable [14]. (It is worth mentioning that for this stability criterion to work N does not need to be a one-to-one map.) In our case, the investigation of g, the stereographic projection of the surface unit normal N, will play an important role in the study of the stability of minimal helical motifs since the corresponding Gauss map is just the quotient of polynomials. Finally, a criterion of instability that will be useful to our purposes is the following: if the surface normal N is a one-to-one map from Σ to S 2 and N(Σ) is a hemisphere, in particular, the spherical area is 2π, then the corresponding minimal surface is unstable [1,14]. (When the area of N(Σ) ⊂ S 2 is precisely 2π, but N is not one-to-one or does not cover a hemisphere, the minimal surface may or may not be stable [15].) (a) Two equally handed helical motifs are unstable The stereographic projection of the Gauss map of a pair of equally handed helical motifs is given by For any given w ∈ C, the equation g(z) = w will generically have two distinct solutions and it then follows that the unit normal N covers the north hemisphere of the unit sphere S 2 exactly twice. (The domain Ω N is chosen such that N takes values on the north hemisphere only.) For the sake of our study of stability, we do not need to count the spherical area with multiplicity and, therefore, the area of the spherical image of any piece r : D ⊂ Ω N → R 3 of a minimal pair of equally handed motifs is bounded by 2π. We are going to show that there exists a bounded domain Dc ⊂ Ω N whose spherical image is precisely the north hemisphere and the corresponding normal is a 1-1 map, from which it follows that pairs of equally handed helical motifs are not local minima of the area functional with respect to normal variations that leave the boundary fixed. In general, any bounded region D ⊂ Ω N is necessarily contained in a sufficiently large disc D(0, ρ) ⊂ C, ρ 1. Under the stereographic projection Π : S 2 → C, a sufficiently small neighborhood of the north pole is mapped on the complement of D(0, ρ), in particular, it is mapped outside the domain D. Now, since g behaves linearly for large values of |z/R|, any sufficiently small neighborhood of the north pole is mapped under g = Π • N outside D ⊂ D(0, ρ) only once! Noticing that g(z = 0) = ∞, i.e., N(0) is the north pole of S 2 , we conclude that there exists a neighborhood of the north hemisphere which is the image under g of some region inside D(0, ρ). Thus, it follows that there should exist a domain Dc ⊂ Ω N such that its image under g is precisely the north hemisphere, as we are going to show below in more detail.
Using that the inverse of the stereographic projection is Π −1 (w) = ( 1+|w| 2 ), we immediately see that the parallels of the unit sphere, i.e., lines of constant latitude, are associated with the level curves of |g|. We then conclude that the unit normal of a pair of equally handed helical motifs when restricted to Dc = {z : r ≤ 1 2 } ∩ Ω N is a one-to-one map on the north hemisphere of the unit sphere. (The second copy comes from considering the outside of the disc {z : r ≤ 1 2 }). In conclusion, the normal is a one-to-one map on the north hemisphere when restricted to Dc and the minimal surface corresponding to a pair of equally handed motifs is unstable for any domain D containing Dc. (If D were stable, any subdomain D ⊆ D would have to be stable.) The proof above is done for the symmetric case, p 1 = p 2 , but similar arguments can be devised to show that any pair of helical motifs such that p 1 = −p 2 has to be unstable. (In fact, neutral total pitch is a necessary condition for stability, as will be shown in Theorem 4.2 below.)

(b) Minimal dipoles are stable
Since the Gauss map of a minimal dipole is a quadratic function on the complex plane, for any w, g(z) = w generically have two solutions and, as in the previous case, the unit normal N covers the north hemisphere of S 2 exactly twice. In addition, the area of the spherical image of any piece r : D ⊂ Ω N → R 3 of a minimal dipole is bounded by 2π. We are going to show that this area is in fact smaller than 2π for any bounded domain D ⊂ Ω N , from which we will conclude that minimal dipoles are local minima of the area functional for normal variations that leave the boundary fixed. Indeed, first notice that any bounded region D of Ω N is necessarily contained in a sufficiently large disc D(0, ρ). Since the north pole is mapped under stereographic projection Π : S 2 → C on infinity, each point on the sphere sufficiently close to the north pole should be mapped on a complex number of sufficiently large modulus. Now, using that g ∝ z 2 for large values of |z/R|, the two copies associated with a sufficiently small neighborhood of the north pole is necessarily mapped under the Gauss map g outside D ⊆ D(0, ρ). In conclusion, it follows that the spherical image of any bounded D ⊂ Ω N is strictly contained in a hemisphere and, therefore, the corresponding area is smaller than 2π. This can be alternatively confirmed by the fact that the parallels of the unit sphere are associated with the level curves of the Cassini ovals. Indeed, the z-coordinate of N = Π −1 • g is given by 1+|g| , which is constant if and only if |g| is constant and, in addition, N 3 → 1 ⇔ |g| → ∞ (See Fig. 3, Left).

(c) Stability criterion for finitely many helical motifs
We have seen on the previous subsections that minimal pairs of unequal handed helical motifs, p 1 + p 2 = 0, are not local minima of the area functional, while minimal dipoles are. Therefore, we arrive at the important conclusion that for two given helical motifs, neutral total pitch is a necessary and sufficient condition for the corresponding minimal surface be a local minimum of the area functional. For more than two helical motifs, extra conditions must be imposed to guarantee stability.
Following , the stability of a minimal surface may be investigated by computing the area under the Gauss map. For n helical motifs glued together, we will proceed as before and study how the complex plane is covered by the stereographic projection of the unit normal, which is a quotient of polynomials of degree n and n − 1. We now formulate a stability criterion for minimal parking garages depending on the position and pitches of the helical motifs. A proof for the first theorem below is based on the observation that n helical motifs are stable if only if g, which in general is a rational function, is a polynomial of degree n. The second theorem is a reformulation of the first and will allow us to check the stability of minimal helical motifs with respect to normal variations that leave the boundary fixed more easily.  Proof. The stereographic projection of the unit normal N of the Enneper graph Σ of n helical motifs is where the coefficients a k and d k are computed from {p i } n k=1 and {z i } n k=1 in terms of the symmetric polynomials as a k = s k (z 1 , . . . , zn) and d k = n j=1 p j s k−1 (z 1 , . . . , z j−1 , z j+1 , . . . , zn). Given w ∈ C, the equation g(z) = w can be rewritten as a polynomial equation of degree n whose coefficients depend on w, {p j }, and {z j }: Thus, g(z) = w generically has n solutions and, consequently, the Gauss map g should cover the north hemisphere of S 2 exactly n times. (The domain of the Enneper graph of helical motifs was defined as Ω N = {z ∈ C : N 3 > 0, N(z) = (N 1 (z), N 2 (z), N 3 (z))}.) On the one hand, if the denominator of g is constant, i.e., d 1 = · · · = d n−1 = 0, then g behaves as ± 2i dn z n for large values of |z/R| and the n copies associated with a sufficiently small neighborhood of the north pole is necessarily mapped under the Gauss map g outside a disc D(0, ρ), for some ρ 1. It follows that any bounded subdomain of Ω N has an area on the unit sphere under N smaller than 2π and, consequently, Σ is stable.
On the other hand, if the denominator of g is not constant, i.e., there exists some k 0 < n such that d k0 = 0, then g behaves as − 2i d k 0 z n−k0 for large values of |z/R|, say for all z ∈ D(0, ρ), ρ 1. This means that a sufficiently small neighborhood of the north pole is covered under N only n − k 0 times. Therefore, the remaining k 0 copies of a neighborhood of the north pole should come from the inside of the disc D(0, ρ). It then follows that there must exist a bounded subdomain Dc of Ω N such that N is 1-1 over Dc and N(Dc) ⊂ S 2 is precisely the north hemisphere. Consequently, Σ is unstable.
The coefficients d k+1 can be rewritten in a more convenient form. For k = 0, we obtain the total pitch Proof. The coefficients d k+1 and b k are related by d k+1 = k j=0 (−1) j b j s k−j (z). In addition, this correspondence can be rewritten as a system of linear equations: Since the determinant of the coefficient matrix of the system above does not vanish, (−1) n = 0, it follows that {b k } can be uniquely computed for a given {d k } and vice-versa. In particular,

(i) Example of stable pitch balanced helical motifs
We demonstrated that two helical motifs forming a minimal dipole is an example of a stable Enneper graph. We now provide another example consisting of helical motif of pitch −np balanced by n helical motifs of pitch p each. It is worth mentioning that this type of configuration resembles the basic units composing the helical geometry in plant thylakoids [6], where the central right-handed helical motif is balanced by several smaller left-handed motifs. Let p 0 = −np, p 1 = p, . . . , p n−1 = p, and pn = p be the pitches of n + 1 helical motifs located at the center and vertices of a regular n-gon, respectively. For the remaining equations, we may use the following identity for the n-th roots of unity ζ 0 = 1, ζ 1 , . . . , ζ n−1 : n−1 j=0 ζ k j = 0 if k is not a multiple of n, but n−1 j=0 ζ k j = n if otherwise [16]. Now, writing z 0 = 0 and z j = R ζ j−1 , j = 1, . . . , n, it follows that (4.10) From Theorem 4.2, we conclude that the configuration of helical motifs {(p 0 = −np, z 0 = 0)} ∪ {(p j = p, z j = Rζ j−1 )} n j=1 , where p = 0 and R > 0, is a stable minimal Enneper immersion.

Discussion
Examining the non-linear geometric interactions between helical motifs in lamellar structures is necessary for understanding the inter-motif interaction at distances comparable to their pitch. Recent experiments revealed such closely packed helical motifs in the Endoplasmic reticulum [3] and plant thylakoids [6]. By providing a constructive method for producing exact minimal surfaces from harmonic functions we are able to examine the geometry and interactions between such closely spaced helical motifs, which are not amenable to a small-slope approximation.
Focusing on asymptotically helicoidal surfaces we proved that for the obtained minimal surfaces to be locally stable minima of the area functional the embedded helical motifs must be pitch balanced. Note, however, that when considering compact domains with boundaries (as is the case for the pitch balanced helical motifs in the endoplasmic reticulum and plant thylakoids) the surfaces may not be asymptotically helicoidal, and thus non-pitch balanced helical motifs may form a stable configuration. Moreover, the stability for pitch balanced arrangements of motifs is obtained by considering perturbations that vanish on the boundaries, and in particular fix the helical motifs in place. In many physical and biological cases the helical motifs are mobile and can vary their relative orientation and position to further lower their energy. The systematic study of these relative forces and torques between motifs is left to future work. Extending the results of stability also requires that we resolve the core of each motif. The Enneper immersion yields in many cases not only the desired minimal surface and helical motif, but also a non-physical section of minimal surface to accommodate "excess material". In the present work, we chose the Gauss Map to define the boundaries between the physical and non-physical portions of the Enneper immersion. Other choices may also be possible and are expected to vary the inter-motif interactions. The constructive recipe provided here allows us to examine arbitrary arrangements of helical motifs at any separation and at any scale. For example, resolving the near field for a finite chain It follows that r is conformal, i.e., a 11 = a 22 and a 12 = 0, if and only if φ 2 1 + φ 2 2 + φ 2 3 = 0. In addition, r is regular if and only if |φ 1 | 2 + |φ 2 | 2 + |φ 3 | 2 = 0. Now, noticing that b 11 = ∂ 2 x r · N = −∂ 2 y r · N = −b 22 , the mean curvature of any conformal immersion r = φ+φ 2 dz vanishes. For an Enneper immersion r = (L −P , h) = Re( (L − P )dz, −i(L + P )dz, 2hzdz), it follows that the metric of the surface given in (2.6) takes the form: a = ∂xr · ∂xr ∂xr · ∂yr ∂xr · ∂yr ∂yr · ∂yr = (1 + hzhz) 2 1 0 0 1 .
It is immediate to verify that the conditions L P = (hz) 2 and |L | + |P | = 0 guarantee that any Enneper immersion is a regular conformal minimal immersion. Alternatively, to prove the minimality, we may use the fact that for a conformal metric a ij = F 2 δ ij , the Laplacian operator corresponding to a ij is ∇ 2 a = 1 F 2 ∇ 2 = 4 F 2 ∂z∂z. Now, using that 2H N = ∇ 2 a r, the mean curvature