Annealing a Magnetic Cactus into Phyllotaxis
Cristiano Nisoli, Nathaniel M. Gabor, Paul E. Lammert, J. D. Maynard, Vincent H. Crespi
AAnnealing a Magnetic Cactus into Phyllotaxis
Cristiano Nisoli , Nathaniel M. Gabor , Paul E. Lammert , J. D. Maynard and Vincent H. Crespi . Theoretical Division and Center for Nonlinear Studies,Los Alamos National Laboratory, Los Alamos NM 87545 Department of Physics Cornell University, 109 Clark Hall, Ithaca, NY 14853-2501 Department of Physics The Pennsylvania State University, University Park, PA 16802-6300 (Dated: November 8, 2018)The appearance of mathematical regularities in the disposition of leaves on a stem, scales on apine-cone and spines on a cactus has puzzled scholars for millennia; similar so-called phyllotacticpatterns are seen in self-organized growth, polypeptides, convection, magnetic flux lattices and ionbeams. Levitov showed that a cylindrical lattice of repulsive particles can reproduce phyllotaxisunder the (unproved) assumption that minimum of energy would be achieved by 2-D Bravais lat-tices. Here we provide experimental and numerical evidence that the Phyllotactic lattice is actuallya ground state. When mechanically annealed, our experimental “magnetic cactus” precisely repro-duces botanical phyllotaxis, along with domain boundaries (called transitions in Botany) betweendifferent phyllotactic patterns. We employ a structural genetic algorithm to explore the more gen-eral axially unconstrained case, which reveals multijugate (multiple spirals) as well as monojugate(single spiral) phyllotaxis.
PACS numbers: 87.10.-e, 68.65.-k, 89.75.Fb
I. INTRODUCTION
Symmetrical morphologies and regular patterns in liv-ing organisms (Fig 1) have been credited with originatingthe idea of beauty, the notion of art as an imitation ofnature, and humanity’s first mathematical inquiries [1–4]. The fascinating symmetrical patterns of organs inplants, called phyllotaxis [1–3], were known to the Ro-mans (Pliny) and ancient Greeks (Theofrastus), whileearly recognitions are found in sources as ancient as theText of the Pyramids [1]. Leonardo da Vinci [5], AndreaCesalpino, and Charles Bonnet [6] studied phyllotaxis inthe modern era. Kepler proposed that the Fibonacci se-quence (1, 2, 3, 5, 8. . . ), where each term is the sum ofthe two preceeding ones [7], describes these phyllotacticpatterns.A discipline that thrived on multidisciplinary interac-tions [8], phyllotaxis found its standard mathematicaldescription when August and Louis Bravais [9] intro-duced the point lattice on a cylinder to represent thedispositions of leaves in 1837 (see Fig. 2), thirteen years before
August’s seminal work on crystallography [10].Unfortunately botanists neglected the work of the Bra-vais brothers, and it wasn’t until Church rediscovered iteighty years later that more progress was achieved in thefield [11].The geometrical description of cylindrical phyllotaxisrelies, in the simplest case, on the phyllotactic lattice in-troduced by the Bravais brothers [1–3, 9]. It consists ofa so-called generative spiral of divergence angle Ω. Wecan visually decompose the resulting lattice in crossingspirals that join nearest neighbors, as in Fig. 2, whichbotanists call parastichies. It is a fundamental observa-tion (made first by Kepler) that the numbers n , m ofcrossing parastichies needed to cover the lattice are con-secutive terms of the standard Fibonacci sequence, or less frequently the variants obtained by changing the secondterm, also called Lucas numbers: 1, 3, 4, 7, 11. . . and 1, 4,5, 9. . . often referred to as second and third phyllotaxis.From that, one can prove that the divergence angle of thegenerative spirals in plants assumes values close to [2, 12]Ω p = 360 ◦ ( τ + p ) , (1)where p = 1 , , τ = (cid:0) √ (cid:1) / n, m ) = ( kn (cid:48) , km (cid:48) ), k being the number of generative spirals [2, 12]. Not un-like domain boundaries in crystals, plants show kinks be-tween domains, called transitions by botanists [1, 13].In the last 50 years, phyllotactic patterns havebeen seen or predicted outside of botany: polypep-tide chains [14, 15], tubular packings of spheres [16],convection cells [17], layered superconductors [18], self-assembled microstructures [19], and cooled particlebeams [20, 21]. While it is still debated whether suchsystems might shed light on botanical phyllotaxis, theoccurrence of such mathematical regularities outside ofbotany is fascinating and leads to generalizations that –unlike quasistatic botany – allows for dynamics [22].In a groundbreaking work Levitov recognized phyl-lotaxis in vortices of layered superconductors [18]. Henext described how phyllotactic patterns represent statesof minimal energy of a cylindrical lattice (that is of a lat-tice with cylindrical boundary conditions) of mutuallyrepelling objects, the repulsion mimicking the interac-tions between spines, leaves, or seeds in plant morphol-ogy [23, 24]. Yet such a constraint to a lattice is ab-sent both in botany and in the physical systems to whichthis energetic model might apply, such as adatoms orlow-density electrons on nanotubes and ions or dipolar a r X i v : . [ c ond - m a t . s o f t ] F e b FIG. 1: Natural and Magnetic Cacti. A specimen of
Mam-millaria elongata displaying a helical morphology ubiquitousto nature, a magnetic cactus of dipoles on stacked bearings,and a schematic of a wrapped Bravais lattice showing the an-gular offset (divergence angle) Ω and the axial separation a between particles. molecules in cylindrical traps.Following up on earlier work that focused on the dy-namics of rotons and solitons in physical phyllotacticsystems [22], we provide here a detailed experimentaland numerical demonstration that Levitov’s constraintis not necessary, and that the lowest energy states ofrepulsive particles in cylindrical geometries are indeedphyllotactic lattices. In addition, we describe the exper-imental and numerical generation of multijugate phyl-lotaxis, static kink-like domain boundaries between dif-ferent phyllotactic lattices, and unusual disordered yetreflection-symmetric structures that may be a static relicof soliton propagation.We show that when a “magnetic cactus” of magnets(spines) equally spaced on co-axial bearings (stem) withsouth poles all pointing outward is annealed, it preciselyreproduces botanical phyllotaxis. When studied numeri-cally via a structural genetic algorithm, the fully uncon-strained case reveals both multijugate and mono jugatephyllotaxis. In addition to our macro-scale implementa-tion, such systems could also be created at the quantumlevel in nanotubes or cold atomic gases.In section II we describe the statics of repulsive par-ticles in cylindrical geometries. In section III we detailthe experiment on the magnetic cactus. In Section IV wediscuss the more general case of multijugate phyllotaxis. II. PHYLLOTAXIS OF REPULSIVE PARTICLESIN CYLINDRICAL GEOMETRIES
In this section we will recall Levitov’s model [23, 24]and some of our own findings [22]. Following Levitov,let us assume that the lowest energy configuration for aset of particles with repulsive interactions, confined to acylindrical shell of radius R , is a helix with a fixed angu-lar offset Ω between consecutive particles and a uniformaxial spacing a , as in Fig. 1 (this so far unproved ansatzwill be investigated later both numerically and experi- FIG. 2: The Bravais lattice with cylindrical boundary con-ditions that defines a phyllotactic spiral. The cylinder axisis vertical, while the horizontal direction contains three cir-cumferential repeats. The solid line is the generative spiral:this one-dimensional Bravais lattice generates the full struc-ture. The dashed lines are the so-called parastichies or visiblesecondary spirals: they connect nearest neighbors on the sur-face of the cylinder. Adapted from A. Bravais and L. Bravais,1837 [9]. mentally). For a generic pair-wise repulsive interaction v ij between particles i and j , the energy of the helix is V = (cid:80) i (cid:54) = j v i,j . Since the lattice structure is defined byΩ, we can write V (Ω).In Fig. 3 we plot V (Ω) for various values of the ratio a/R : for specificity we employed a dipole dipole interac-tion v i,j = p i · p j /r i,j − p i · r i,j )( p j · r i,j ) /r i,j , repulsiveat the densities considered here. However, the followingconsiderations only depend upon geometry and thereforeapply to a vast range of reasonably behaved, long rangerepulsive interactions.When a/R (cid:29)
1, the angle Ω = π maximizes distancebetween neighboring particles and therefore V (Ω) has aminimum in π . The angle between second nearest neigh-bors along the helix is 2 π , which means that they faceeach other. And thus, as the density increases, inter-action between the facing second nearest neighbors be-comes predominant, and Ω = π is not a minimum for V (Ω) anymore. If whe shift the helical angle from π , therepulsive interaction between second nearest neighbors isreduced, with minimal penalty from nearest neighbors.In terms of V (Ω), that means a local maximum Ω = π .This argument can be iterated for every commensu-rate winding that allows particles separated by j neigh-bors to face each other. As density increases further,the angles 2 π/ π/ πi/j with i, j relatively primecorresponds to a configuration where each particle faceseach jth neighbor. For every j there will be a value of a/R low enough such that Ω = 2 πi/j is a local maximum,which we call a peak of rank j .The proliferation of peaks for increasing linear densityis shown in Fig 3. We can see that at high density, peaksof equal rank are nearly degenerate; that is natural, sincetheir principal defining energetic contribution arises fromparticles facing each other at a distance ja . The minimaalso become more nearly degenerate as the density in-creases. That can be explained intuitively, since for an-gles incommensurate to π each particle “sees” the othersas incommensurately smeared around the cylinder, andis therefore embedded in a nearly uniform backgroundcharge from the other particles. The degenerate energyof the ground state can be well approximated by an uni-form continuum distribution (cid:15) , whereas the energy of apeak of order j will be V (2 πi/j ) (cid:39) v ( ja )+ (cid:15) , where v ( r )is the energy of two particles facing at a distance r : forour dipole interaction v ( ja ) = p /a j .The first step to calculate the degeneracy of our sys-tem at a given density, is to find the corresponding max-imum rank of the peaks. As all of the peaks of the samerank have the same energy, and appear in the spectrumtogether, we can focus on the emergence of 2 π/J . For a/R (cid:28)
1, this new peak will emerge when the distancebetween particles separated by a distance Ja equals thatof particles separated by 2 πR/J . Therefore one finds forthe maximum rank J = (cid:104)(cid:104)(cid:114) πRa (cid:105)(cid:105) , (2)which as expected only depends on purely geometricalparameters. A little more tricky is to compute the de-generacy, given J . The set of all the peaks has the car-dinality of the class of all the fractions i/j , with i, j coprime and j ≤ J . This can be considered as the dis-jointed union of other classes, called Farey classes of orderj, defined as follows: P j ≡ { Ω = 2 πi/j | for i, j coprimeand i ≤ j } , i.e. all fractions in lowest terms between0 and 1 whose denominators do not exceed j [25]. Theunion of all Farey classes up to a certain order J has thecardinality of the set of peaks for a spectrum of maxi-mum rank J . Now, the cardinality of P j is know fromnumber theory to be Euler’s totient function, φ ( j ) [26].Therefore, the degeneracy D of the energy minima for asystem with a maximum peak rank J is [26] D = J (cid:88) j =1 φ ( j ) = 3 π J + O ( J log J ) , (3)which, from Eq. 2, scales as D ∼ R/a .Finally, we recall [3, 18] that the order j , j of thepeaks bracketing a minimum relates to its structure ina straighforward way: the helix corresponds to a rhom-bic lattice where each particle has its nearest neighborsat axial displacements of ± aj , ± aj and second nearestneighbors at ± a ( j + j ) or ± a ( j − j ) [18]. Also, j and j give the number of crossing secondary spirals (paras-tichies) needed to cover the lattice by connecting nearestneighbors. !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! " ! Π !!!!!Ε !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! " ! Π V !!!!!Ε FIG. 3: Lattice energy V (Ω) versus divergence angle for suc-cessively halving values of a/R starting from 0 . (cid:15) = p /a , where p is the magnetic dipole).Notice the proliferation of peaks as a/R decreases. Repro-duced from [22]. For completeness, let us now follow Levitov [18, 23, 24],and consider the adiabatic evolution of our system asthe linear density is increased. As new sets of maximaand minima emerge, the true minimum goes througha series of quasi-bifurcations, the consequence of anelusive symmetry whose explanation goes beyond ourscope. Suffice it to say that the system evolves quasi-statically from one of these optimal Ω to another as
R/a increases, asymptoting to the golden angle Ω =2 π/ ( τ + 1) (cid:2) τ = (cid:0) √ (cid:1) / (cid:3) , ubiquitous in botany, aseach minimum is bracketed by peaks whose ranks, be-cause of the Farey tree structure described above, areconsecutive elements of the Fibonacci sequence. Occa-sional “wrong turns” at later stages, will not shift theconvergence too far from the golden angle, yet the Fi-bonacci structure would be lost. However if one or twoconsecutive wrong turns happen at the second or secondand third bifurcations the system will converge to the al-ternative angles of second or third phyllotaxis, given byEq. (1).We have only surveyed so far spiraling lattices gener-ated by a single helix. A straightforward generalizationgives multijugate phyllotaxis, when two or more elementsgrow at the same axial coordinate [1–3]. This case, whichLevitov does not explore, can be easily mathematicallyreduced to monojugate case, by considering two or morereplicas of the phyllotactic lattice as in Fig. 2. In our ex-perimental realization we restrict ourselves to the mono-jugate phyllotaxis, and explore multijugate only numeri-cally. III. A MAGNETIC CACTUS
There is a long history of experimental reproductionsof phyllotactic patterns. Recently, Doady et al. describedphyllotaxis in terms of dynamical systems and then veri- f d FIG. 4: Experimental apparatus. Top: Each unit of the mag-netic cactus consists of a magnet element and a unit ringsecured to a central axis. The ring diameter d is 2.2 cm.Bottom: a schematic representation of the mounted magneticcactus and surrounding measurement devices. The viewer’seye is restricted by the viewing slit and the reference wires.Measurements are taken directly from the dividing head. fied it experimentally by examining dynamical processesin droplets of ferro-fluid [27]. But even more than a cen-tury ago, Airy showed that phyllotaxis emerged in opti-mal packing of hard spheres connected by a rubber band,once the band was twisted to increase density [28].Here we expand on what was announced in a recentlypublished Letter [22]: we verify experimentally the as-sumptions of Levitov’s energetic model, by studying thelow energy configurations of interacting magnets stackedevenly-spaced and free to rotate around a common axis.We constructed a mechanical system that it is free toexplore the three angles of botanical phyllotaxis (Eq. 1). A. Experimental Apparatus
We built a magnetic cactus by mounting permanentmagnets (spines) on stacked co-axial bearings (a stem)which are free to rotate about a central axis, as in Fig. 4.All the magnets point outward, to produce a repulsiveinteraction between all magnet pairs. To avoid effectsof gravity, the apparatus rests in the vertical position,and is non-magnetic. We built two different versions, thesecond with magnets twice as long as the first, as to havea larger effective radius which gives three rather than twostable structures.We employed cylindrical permanent iron-neodymiummagnets, 1.2 cm long and 0.6 cm in diameter. They aremounted on fifty aluminum rings of 2.2 cm outer diam-eter, each affixed to a non-magnetic radial ball bearing(acetal/silicon, Nordex) as in Fig. 4. These unit rings areevenly spaced on an aluminum rod in a stacked structureof 39.9 cm axial length.At static equilibrium we measure separation angle be-tween each magnet element, by rotating the cactus un-til a magnet element aligns with the reference wires. Atelescope and a vertical viewing slit accompanied by twovertical reference wires assist in data acquisition.
B. Annealing
By measuring the dipole-dipole interaction between anindividual magnet pair, we can reconstruct the curve ofthe lattice energy V (Ω) as a function of the angular offsetΩ between magnets. We find that the first arrangement,with short magnets, admits two minima, given by theangles of Eq. 1 for p = 1 ,
2. The second arrangement,with long magnets, has three minima corresponding tothe angles of Eq. 1 p = 1 , ,
3, one of which ( p = 3) is aweak metastable minimum. These divergence angles ofstable helices are very close to those predicted by phyl-lotaxis, of Eq. 1, and are all accessible by experimentalprocedure described below.Before every data acquisition, the cactus is disorderedand then athermally annealed into a low-energy state.The protocol involves repeatedly winding the bottom-most magnet to generate an ever-tightening spiral, un-til an explosive release of energy disorders the lattice.Next, an independent external magnet is oscillated insmall circular motions near randomly chosen points whilethe cylinder as a whole is slowly rotated, to further ran-domize magnet orientations. After 10-30 second of me-chanical annealing through applied vibrations, the sys-tem consistently enters a robust ordered state which doesnot anneal further on experimental timescales. C. Results
Figure 5 reports the experimental results for both ar-rangements by plotting the measured angle between con-
Axial Index o o o o o ! (deg) ! ! ! " E(mJ) ! " ! ! " ! ! Axial Index o o o o o ! (deg) ! ! ! " ! E(mJ) !" ! ! " ! ! ! FIG. 5: Top: A 3-D rendering of the experimental data andthe corresponding Bravais lattice of the magnetic cactus an-nealed in a spiral configuration of divergence angle Ω andparastichies (2,3) (blue and red dashed lines), Fibonacci num-bers. The bottom two panels show the experimentally mea-sured angular offsets Ω between successive magnets for mag-netic cacti with short (middle) and long (bottom) magnets,plotted versus magnet index, which simply counts the num-ber of magnets along the axis. Flat regions are perfect spiralswhile steps are boundaries between different phyllotactic do-mains. The dotted lines give the phyllotactic angles Ω , Ω ,Ω and 2 π − Ω defined in the text, whereas the dashed linesare minima of the magnetic lattice energy (insets) calculatedby interpolating the measured pair-wise magnet-magnet in-teraction. Data reproduced from Ref. [22]. Axial Index o o o ! (deg) ! ! ! " ! FIG. 6: A kink between domains of first and second phyl-lotaxis. From top to bottom: a 3-D rendering, and its un-wound 2-D Bravais lattice from numerical simulations. Below,the same kink plotted as angle increments between successivemagnets for the experimental system (crosses) and numericalsimulation (circles). secutive magnets. The more narrow (short-magnet) cac-tus self-organizes into the spirals with divergence anglesprecisely reproducing those of first phyllotaxis, Ω = Ω ,and second phyllotaxis, Ω = Ω , as in Eq. (1). When theresults are represented in a 2-D lattice, as in the top ofFig. 5, parastichies can be drawn. As parastichial num-bers for Ω = Ω we find the Fibonacci numbers (2 , , the Lucas numbers (3 , and Lucas numbers (1 , V (Ω) of thelattice obtained by interpolating measured values for thepair-wise magnet-magnet interaction, plotted as a func-tion of divergence angle Ω. As we cans see, local minimacorrespond to phyllotactic angles.Figure 5 also shows that in many instances the systemfragments into two or three distinct domains whose do-main walls always share a common parastichy, as seen inbotany [1, 13], and as expected in physics for a quasi-one-dimensional degenerate system. We have computed nu-merically one such transition via dynamical simulationsin a velocity-Verlet algorithm, in the following way: westart from a crude static step-like kink as an initial condi-tion, and allow it to radiate energy in the form of phononswaves until it stabilizes in a kink with superimposed vi-brations; we then average this configuration over time,to remove these residual oscillations. When the result is !"""""""""" Π !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! " ! Π !!!!!Ε !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! " ! Π V !!!!!Ε FIG. 7: A numerically calculated kink in a system of dipolesof high degeneracy (seven minima, a/R = 0 . / π = 1 / , / / π = 1 / , / (cid:15) = p /a , p being the magnetic dipole. used as new initial conditions, it proves to be a statickink. Fig. 6 reports our numerical results for a kink in asystem whose size and interaction reproduces the physi-cal realization of the magnetic cactus, along with the ex-perimental data for such a kink. The match is essentiallyperfect, indicating that the dissipative (i.e. frictional)forces neglected in our model do not significantly affectthe static configurations. We apply the same numericalprocedure to calculate a kink in a system of larger degen-eracy, among domains which are absent in our physicalrealization. We use a smaller a/R ratio and a differ-ent interaction between particles (ideal dipole instead ofphysical dipole). The result shown in Fig. 7 reproducesthe same qualitative shape of the previous, lower degen-eracy case. Similar kinks are present also in a fully un-constrained cactus, one in which the particles can movealong the axis, and are found in early generations of ourstructural genetic algoritms (see below). Finally, thesekinks can travel as novel topological solitons, with a richphenomenology that is explained elsewhere [22, 29]. Axial Index o o o o o ! (deg) ! ! ! " ! E(mJ) !" ! ! " ! ! ! FIG. 8: Measured angular offsets between successive magnetsfor magnetic cacti with long magnets, plotted versus mag-net index, showing symmetric kink/anti-kink domain bound-aries. Dotted lines give the phyllotactic angles Ω , Ω , Ω and 2 π − Ω defined in the text. Dashed lines are minima ofthe magnetic lattice energy (inset) calculated by interpolatingthe measured pair-wise magnet-magnet interaction. Finally, in the system with longer magnets, we occa-sionally found intriguing yet hard-to-interpret configura-tions that contain two nearly reflection symmetric do-main boundaries. Fig. 8 reports two such configurations,measured in independent experimental runs. Althoughwe do not have a firm explanation for these structures,we speculate that they form as frozen-in soliton wavesthat initiated symmetrically at both ends of the struc-ture, upon release of the wound-up elastic energy duringinitial preparation. Indeed an analytical, continuum the-ory for phyllotactic solitons which we have developed re-cently, and which explains results of the dynamical symu-lations also supports the existence of similar frozen-inpulses [29].
IV. FULLY UNCONSTRAINED CACTUS:STRUCTURAL GENETIC ALGORITHM
Our experimental apparatus is not fully unconstrained:the axial coordinates of the dipoles are fixed, and soonly the azimuthal movement is allowed. While this isan huge improvement toward the original helical con-straint of Levitov, many (most) physical systems thatcould manifest phyllotactic patterns do not posses sucha lesser azimuthal constraint. To corroborate and extendour experimental results to a completely unconstrainedsystem, we seek the energy minimum in a set of repulsiveparticles that can move axially as well as angularly on acylindrical surface, via a non-local numerical optimiza-tion. To this purpose, we developed a structural geneticalgorithm.
A. Genetic Algorithm
A genetic algorithm is a method of optimization thatmimics evolution to find the absolute minimum in a func-tion which shows a large number of metastable minima.The coordinates of the energy functions are called genes,and a set of genes is a particular specification of valuefor those variables. The routine typically starts with aset of “parents”, or specific points in the domain of theenergy function. At each stage of the routine, parents“mate” to produce children via exchange of genes: a sub-set of the coordinates of each of the two configurationsare swapped, therefore generating new points in the en-ergy domain, called children. Each of those children isthen locally relaxed to a minimum via a local search.The new population of parents and children undergoesgenetic selection and only the fittest (the lowest energyones) form a new population.There are many different implementations of this gen-eral idea: care is taken not to lose genetic diversity duringselection, to avoid a population of almost identical repli-cas; that is usually achieved with a more or less skilledgenetic selection, which might retain less geneticaly fitindividuals, and often by introducing mutations in theform of random alteration of the gene sequence, whichwould opefully prevent the routine from gettting stuckaround a metastable region. Choice of parameterizationof the structures (genes) and mating (crossover) is crucialto the performance of the algorithm.About fifteen years ago, Deaven and Ho [30] introduceda so-called structural genetic algorithm, which provedparticularly efficient in minimazing the energy of physicalstructures, as it allows for physical intuition in definingthe genes and mating procedure. With it, they found theC fullerene structure as a ground state of 60 carbonatoms interacting with suitable atomic potentials [30]and solved the celebrated Thomson problem of repulsivecharges on a sphere [31], a task quite similar to ours. B. Our Algorithm
In our implementation we use a population of 10 mem-bers. Each member reppresents a configuration of 101particles on the cylinder: more esplicitly, the geneticstructure of each member P k , k = 1 , . . . ,
10, of the popu-lation is a set of variables, or P k = { θ i , z i } i =1 which spec-ifies, in cylindrical coordinates, the positions of the par-ticles composing its structure. The particles interact viaa pair-wise inverse quadratic repulsion V = V o ( r o /r ) ,where r is the three dimensional distance between parti-cles ; we introduce a confining potential in the form of FIG. 9: Relative energy of the fittest member of the popula-tion in every generation. Early generations return very highenergy configurations corresponding to disordered metastablestates. Interrmediate generations show populations of phyl-lotactic domains separated by kinks between. Finally, thealgorithm converges to a single crystalline domain in the bulk(deformations at the boundaries simply accomodate the sys-tem to the confining potential). an external axial square-well of width L , which sets thelength of the cylinder, and hence the density. The choiceof the pairwise interaction is not fundamental, as long asit is long ranged, repulsive, and well behaved [22]; ourparticular choice simply speeds up the computation.We generate the first population randomly. At eachstep, we randomly couple mates, and excahge their genesemploying the following mating procedure: we order thegenes by increasing axial coordinates z < z < · · · 101 genes, where n isa random number, between randomly selected parents.The children obtained in this way are then relaxed to astable structure via a standard conjugate gradient algo-rithm. We then prepare the new generation by selectingthe lowest energy individuals in the population of parentsand children, yet making sure that the energy differencebetween members does not fall below a certain threshold,to preserve genetic diversity: when new children cannotproduce a new population of 10 in accordance with theenergy threshold, we introduce mutations by randomlyaltering a certain number of members. C. Numerical Results During the structural evolution, the earliest popula-tions contains metastable disordered states. Members ofintermediate populations show kinks between domains ofdifferent divergence angle, configurations which are alsoseen experimentally. After fifteen to twenty generations,the algorithm typically converges to a single crystallinedomain.Figure 9 reports the energy of the fittest member ofthe population at each generation in a typical run, show-ing a punctuated-equilibrium evolution where the most- Axial Index o o o o o ! (deg) " z (L/N) ! FIG. 10: Numerical optimization via structural genetic al-gorithm for N = 101 repulsive particles [ V = V o ( r o /r ) ]constrained to a cylindrical surface of length L and radius R = 1 . L/N . The resulting 2-D Bravais lattice has a nearlyconstant axial separation ∆ z = z i +1 − z i (top) and angulardivergence Ω between successive particles (bottom), neglect-ing fringe effects at the border of the potential well. In thebulk, particles self-organize on a single spiral of divergenceΩ = Ω . Oscillations at the boundaries are due to the effectof the confining potential. fit structure progressively decreases in energy in intermit-tent steps separated by plateaux. The final converged re-sults form well-defined two-dimensional cylindrical crys-tals away from boundaries.Figure 10 shows the crystalline structure to which thealgorithm converges, for R = 1 . L/N : a single spi-ral with Ω = Ω , as defined in Eqn. (1), correspond-ing to first phyllotaxis with parastichies (1 , z = z i +1 − z i returns the value L/N in the bulk, whichimplies a single generative spiral. This choice of RN/L corresponds to a density close to that of our experimentalapparatus. V. MULTIJUGATE PHYLLOTAXIS For highly degenerate systems the genetic algorithmreturns configurations with more than one generativespiral, corresponding to what in botany is called mul-tijugate phyllotaxis [2]. We have seen before that he-lices make cylindrically symmetric lattices. On the otherhand, every cylindrically symmetric lattice can be de-composed into a suitable number of equispaced genera-tive spirals [2, 9]. That is accomplished by discretizingthe cylinder along its axis into equally spaced rings andthen assigning at each ring n sites, equally spaced andseparated by a 2 π/n angular shift. As before, each ringis shifted consecutively by a divergence angle Ω. The case a ! aa A x i s Circumference ! ! Ω π/ Ω !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! !!!!! " ! Π !!!!!Ε FIG. 11: Top: 3-D schematics and 2-D lattice for a 2-jugateconfiguration. Bottom: Lattice energy versus angular offsetfor n-jugate configurations in a system of repulsive dipoles( (cid:15) = p /a , where p is the magnetic dipole) for 1-jugate(black), 2-jugate (red), 3-jugate (green), and 4-jugate (blue). n = 2 is shown at the top of Fig. 11. The case n = 1 isshown in our experimental arrangement of Fig. 1.By decomposing the n-jugate cylindrical lattice into n lateral replicas of single-spiral lattice, as in Fig. 2, thereader is easily convinced that multijugate phyllotaxis re-duces to the previously described monojugate case. Allthe considerations above apply, provided that one nowtakes the periodicity to be 2 π/n , and the distance be-tween rings to be na (with, as before, a = L/N ). Itfollows that a n-jugate configuration will have local max-ima in 2 πi/n , i = 1 . . . n and, following the discussion ofsection II, one finds that there will be other local maximacorresponding to angles 2 π/n × i/j when j ≤ [[ J/n ]], and J is the maximum rank given by Eq. 2.Note now that if two multijugate lattices of jugation n , n (cid:48) have a peak in the commensurate angle 2 πi/j , thenthe energy of the peak is the same, as is shown in Fig. 11,bottom, which compares the plots of the energy of suchan arrangement for different values of n . In fact bothconfigurations correspond to particles facing each otherafter j/n and j/n (cid:48) rings, and therefore at the same dis-tance na × j/n = n (cid:48) a × j/n (cid:48) = ja , independent of n or n (cid:48) .For small n the minima in the energy of n-jugate config-urations essentially degenerate with the monojugate onepreviously explored. For large n they have higher en-ergy. If a/R is small enough, the threshold is n > J , asinteraction between particles on the same ring becomecomparable to those facing in the minimal monojugatepeak. VI. CONCLUSION We have studied the lowest energy configurations ofrepulsive particles on cylindrical surfaces, both experi-mentally and numerically. We have found that they cor-respond to the spiraling lattices seen in the phyllotaxisof living beings, both monojugate and multijugate. By establishing experimentally and numerically that phyl-lotactic point lattices are ground states in the very gen-eral geometric scenario of unconstrained repulsive parti-cles on cylinders, we have opened the study of phyllotaxisto a much wider range of annealable physical systemswhere the particles could be electrons, adatoms, ions,dipolar molecules, nanoparticles, etc. constrained by ex-ternal potentials.Unlike plants, these multifarious, non-biological Phyl-lotactic systems could access various degrees of dynamics,providing new phenomenology well beyond that availableto over-damped, adiabatic botany. We have reportedelsewhere [22] on the dynamical richness of this physicalphyllotaxis, including classical rotons and a large familyof novel, inter-converting topological solitons. [1] I. Adler, D. Barabe, R. V. Jean, Ann. Bot. Phyllotaxis: A Systemic Study in Plant Mor-phogenesis (Cambridge Univ. Press, Cambridge, 1994).[3] http://maven.smith.edu/ ∼ phyllo/ provides an overviewand useful applets.[4] N. Grew, The anatomy of plants (Johnson Reprint Corp.,New York, 1965).[5] E. MacCurdy, The notebooks of Leonardo da Vinci (Braziller, New York, 1955).[6] C. Bonnet, Recherches sur l’usage des feuilles dans lesplantes (E. Luzac, fils, G¨ottingen and Leiden, 1754).[7] L. E. Sigler, Fibonacci’s Liber Abaci (Springer-Verlag,New York, 2002).[8] V. Jean, B. Denis, Multidisciplinarity: a key to phyl-lotaxis in Symmetry in Plants On the Relation of Phyllotaxis to Mechan-ical Laws. (Williams and Norgate, London, 1904).[12] I. Adler, Journal of Theoretical Biology 596 (1954).[15] S. F. Albdurnur, K. A. Laki, J. Theor. Biol. , 186103 (2009).[23] L. S. Levitov, EuroPhys. Lett. 533 (1991).[24] H. W. Lee, L. S. Levitov, Universality in Phyllotaxis:a Mechanical Theory in Symmetry in Plants , 385,(1816).[26] J. H. Conway, and R. K. Guy, The Book of Numbers. (Springer-Verlag, New York, 1996); T. Nagell, Introduc-tion to Number Theory (Wiley, New York, 1951).[27] S. Douady, Y. Couder, Phys. Rev. Lett. 176 (1873).[29] C. Nisoli, Phys. Rev. E , 026110 (2009).[30] D. M. Deaven, and K. M. Ho, Phys. Rev. Lett.53