Controlled transitions between phyllotactic states of repulsive particles confined on the surface of a cylinder

Phyllotactic states are regular lattice-like structures on cylinders and are a botanical classification scheme. In this communication, we report a sequence of transitions between phyllotactic states for particles with a repulsive particle-particle interaction on a cylindrical geometry at zero temperature. We can infer the transition points as a function of density via Monte Carlo simulations, as well as the mathematical descriptions of the ground states. The lattices we generate are described as phyllotactic states that fit onto the cylindrical surface as a set of helical chains. Our analysis shows how all state energies lie on the same parabola which we exploit to find the transitions.

Phyllotactic states are regular lattice-like structures on cylinders and are a botanical classification scheme. In this communication, we report a sequence of transitions between phyllotactic states for particles with a repulsive particle-particle interaction on a cylindrical geometry at zero temperature. We can infer the transition points as a function of density via Monte Carlo simulations, as well as the mathematical descriptions of the ground states. The lattices we generate are described as phyllotactic states that fit onto the cylindrical surface as a set of helical chains. Our analysis shows how all state energies lie on the same parabola which we exploit to find the transitions. The same state is projected into two dimensions where we see the origin of the phyllotactic notation. The periodicity vector, V =ŷc = ma + nb, has a length equal to the circumference and is a combination of lattice vectors [10]. The thick horizontal lines represent the periodic boundary condition. In this example, V = 3a + 2b.

I. INTRODUCTION
In this communication, we propose an algorithm that enables us to analytically construct the infinite sequence of transitions between phyllotactic states of repulsive particles confined to the surface of a cylinder. We build upon research on the close-packing of spheres inside cylindrical tubes [1][2][3][4] and disks in narrow channels [5,6], minimum energy structures of nano-particles in carbon nanotubes [7][8][9], and the arrangement of some areoles on cacti. All of these systems share the common mathematical description of phyllotaxis.
Phyllotaxis is a biological classification used to described the pattern formation of leaves on plant stems [11]. The phyllotactic notation treats the cylindrical stem as an "unravelled" two-dimensional structure with a triangular lattice (formed by two lattice vectors) with node sites corresponding to petiole (stalk) locations [10]. This can then be described as a set of N chains formed by following m and n multiples of the two lattice vectors to wrap once round the cylinder to the original node. The notation is expressed as [m + n, m, n] or [m, n, m − n], where m, n ∈ N. We distinguish these two representations because m and n generate one of the two phyllotactic structures depending on the system geometry and density of particles.
Optimal packing of spheres in cylinders and hard disks in periodic geometries are purely geometric problems. In three dimensions, shell-like structures of surface or core spheres form and are well described by phyllotaxis [12]. Similarly, hard disks form phyllotactic structures because of the periodicity vector defined in Fig. 1. The ground states of these systems can include helical grain boundaries (line-slips) [13]. Line-slips can appear due to geometrical constraints which are not present in our softmatter-type potential.
Related work [14] on parabolic confinement demonstrates the number of rows of repulsive particles that form become dependent on the background potential and interparticle interactions. Transitions between the number of rows of parabolically confined particles follow similar trends to the helical row transitions that we report despite the differing boundary conditions.
Interactions may be a function of the threedimensional distance between particles [15]. Alternatively, enforcing periodic boundary conditions in one direction causes interactions to depend on the arc lengths on the surface of the cylinder. Our study is the latter choice, although our results show qualitative similarity for both interaction distances.
We begin by numerically determining the energetics of the ground state lattices close to the structural transition points. The analysis leads to the conjecture of a scale invariant, α = c 2 ρ 2 , where c is the circumference of the cylinder and ρ 2 is two-dimensional particle density. The ground state is exclusively determined by α -and this is verified by our numerical results.

II. SYSTEM
We model the cylinder as a two-dimensional planar system of repulsive particles with a periodic boundary condition imposed in the circumferential (y) direction. This model has been used previously for studying particle interactions on cylinders [16]. In the numerical simulations, the length of the cylinder is set to be much greater than the circumference, c, to approximate an infinitely long cylinder by imposing an additional periodic boundary condition in the axial (x) direction, modelling the cylinder as a torus. We study the bulk behaviour of the numerical system to make phenomenological inferences in our analytical model which assumes an infinitely long cylinder. Structural ground states that form in the bulk when L 10c are well-explained in the analysis. This is a zero temperature system.
To establish that the behaviour is not special to a particular interaction potential, we have worked with both a modified Bessel-function of the second kind, K 0 (r/λ), and the Yukawa-potential, exp(−κr)/r, where r denotes particle separation. These potentials are representative of two-dimensional soft matter system namely vortices in superconductors and Wigner Crystals. Qualitative similarities imply a generality to our results. Numerical results in this letter set κ = 1 and λ = 1 (arbitrary units).

III. INITIAL NUMERICAL RESULTS
We use a Metropolis Prescription of the Monte Carlo algorithm to anneal the system to zero temperature. Circumference and linear density, ρ 1 (number of particles per horizontal length along the cylinder), are varied and an initial phase diagram is generated. The linear density is appropriate since the two-dimensional density can be scaled out using ρ 2 = ρ 1 /c. Figure 1 is a typical observed state. Initial results indicate that the bulk of the system forms an isosceles (or equilateral in special cases) triangular lattice structure. The lattice is regular and unit cell sizes remain constant.
Because the system is numerically finite, commensurability effects can cause highly localised grain boundaries or lattice defects to form. This is due to helically defined states being unable to align correctly with themselves, so a discontinuity forms. However, the bulk of the system remains homogeneous which is where we derive our phenomenology.
Simulating between 0.5 ≤ c ≤ 5 and 0.5 ≤ ρ 1 ≤ 5 with a resolution of 0.1 × 0.1 shows that transition lines between states take the form c = α T /ρ 1 , where α T is a number to be determined for each transition. Figure 5 shows the phase diagram generated by our subsequent [  analysis, which is qualitatively similar to these initial results.

IV. PHENOMENOLOGICAL MODEL
The three key observations from the numerical results we can build a phenomenological model with are: • constant unit cell size, • the periodicity of the lattice in the vertical direction, • and an isosceles lattice. We combine these features to create a phenomenological model in terms of the lattice vectors.
TABLE I. For a regular isosceles lattice under periodic confinement in the y direction, the lattice vectors a and b are constrained as below,

Constraint
Physical Reason z · (b × a) = c/ρ 1 Constant unit cell size ma + nb =ŷc Periodicity vector |a| = |b| Isosceles triangular lattice Without loss of generality, the z−components of vectors a and b are set to zero. The constant unit cell size is defined using the linear density (usually 1/ρ 2 ). The periodicity means an integer (m and n) vector sum of a and b must arrive at the periodicity vector,ŷc, defined in Fig. 1. Having two vector lengths equal forms an isosceles lattice. Solving for a and b: where p = m 2 − n 2 . Although there are multiple solutions, they generate the lattices which are reflected on the x axis which correspond to degenerate energies and are not considered further. Note the factor of 1/ρ 1 and cρ 1 = α, where α is a continuous variable. For any given value of α, the angles between lattice vectors are conserved, since a·b/|a||b| only depends on α, and the absolute scale of the lattice is thus set by 1/ρ 1 . Furthermore, the numerical transitions are consistent with α T = cρ 1 (specific values of α that match up to transition curves).
For example, at the (m = 2, n = 2) → (m = 2, n = 1) transition, we construct a nonlinear model of the form ρ 1 = α/c which fits with α = 6.197 with an R-squared value of 0.999998. Note for this transition, the exact value of α = 8 3/5 ≃ 6.19677. Fitting the curves as inverse functions, we search along the line c = ρ 1 and cross each transition point once (we set c = ρ 1 = √ α, with dimensional constants absorbed by λ or κ). At more extreme scalings, one would expect this method to break down as the details of the interactions will cause higher order effects. We then express the lattice vectors as: For lim m→n : Three primitive lattice vectors are needed for phyllotactic notation. Vectors a and b are primitive, the third is c − = a − b or c + = a + b. |c − | = |c + | corresponds to a square lattice. This is only ever the ground state when (m = 1, n = 1) for α = 2, or α = m 2 + n 2 . When α < α , c + is primitive and when α > α , c − is primitive. Given the primitive vectors, the accompanying phyllotactic description is [m+n, m, n] when α > α and [m, n, m − n] when α < α . We use this notation alongside (m, n) for clarity between our original constraints and the phyllotactic nature of the problem.
We calculate the energies generated by a set of m and n for a particular value of α to find the ground state. In Appendix A, we show the maximum value of m needed to capture the lowest energy behaviour is ⌊ √ 2α max / 4 √ 3⌋. Where α max is the maximum value of α we search with. Figure 2 shows the ground states found for 0 < α < 25.
By numerically searching for ground states and transitions, we are able to categorise transitions and locate degenerate points.

V. ANALYSIS OF THE PHENOMENOLOGICAL MODEL
To explore the qualitative generality of our results we compare the numerical results for the Bessel-function and the Yukawa-potential highlighted in Fig. 2. The largest difference occurs when the gradients of the transitioning energies are similar. This is due to the energy difference of each state being small over a larger region than the average. We will see that these transitions have a slight dependence on the interaction potential which can cause small discrepancies in transition values. Both potentials show instantaneous degenerate ground states where two states are equivalent for a critical value of α but no associated structural transition is observed.
Degenerate ground states form when, for a given value of α(= α △ ), two states both form an equilateral triangular lattice. States can be found exactly by setting |a| = |a ± b|. In the two-dimensional phase diagram, this corresponds to lines of c = α △ /ρ 1 where the energy is instantaneously degenerate. We find: The lower value of α △ only exists if m 2 −4mn+n 2 < 0. We find no transitions at α = α △ , only instantaneous degeneracies where multiple states can exist if multiple values of m and n yield the same value of α △ .
Two types of structural transition exist: • two states sharing the exact same structure up to a global rotation or • two states with different structures that exist on opposite sides of an energetic parabola.  Since our analysis is along the line c = ρ 1 , the unit cell size is 1 and |a| = |b| still holds. Using these constraints without the periodicity needing to be satisfied, we generate a simpler description of the lattice, without any loss of generality, which is arbitrarily rotated so the chains are parallel to the horizontal axis. We parametrise this simpler lattice with γ:

VI. CALCULATING TRANSITION POINTS
This can be directly mapped to represent any of the full set of states via rotation and scaling. Figure 4 illustrates the variety of lattices generated. By varying γ and numerically calculating ground state energies, we find a double-minima curve where the minima correspond to values of α △ . Alternatively, if the energy is plotted against the vector length, L = |a| = 4γ 4 + 1/(2γ), a parabola is found where the variation of γ traces out a trajectory that doubles back on itself (Supplementary Video 1 demonstrates this in detail). We can then classify the structure of each state based on its position on this parabola and lattice vector length. Figure 4 shows the energy using Bessel function interactions, which is parabolic to leading order. Translating the curve by plotting as a function of (L − L 0 )/L 0 , where L 0 = 2/ √ 3, causes it to be symmetric at the origin. L 0 is the vector length corresponding to the per- with respect to scaled vector length, where L0 = 2/ √ 3 is the vector length that yields the minimum energy. The numbers correspond to the region or point on the curve where certain types of lattice form. State types 2 and 6 are the perfect lattice, 4 is the square lattice, and 1, 3, 5, and 6 are variations that describe the relative proximity of the nextnearest neighbour to the circular radius equal to the lattice vector length. Any lattice can be categorised by checking which of the following inequalities holds for the numbered state types: 1: |a| < |a − b|, 2: |a| = |a − b|, 3: |a − b| < |a| < |a + b|, 4: |a + b| = |a − b|, 5: |a + b| < |a| < |a − b| 6: |a| = |a + b|, 7: |a| > |a + b|. The parabolic fit is accurate in the region near to the minimum of the curve. Supplementary video S1 shows the structure as a function of γ and L in enhanced clarity.
fect equilateral triangular lattice when γ △ = 1/ 2 √ 3 or γ △ = √ 3/2. The Yukawa potential similarly yields a parabolic curve but a different minimum energy value when L = L 0 , further indicating the qualitative consistency for different potentials.
When two equal energies are on opposite sides of the origin, we define a separation distance between points meaning the lengths of the vectors at a transition point differ by 2d. This leads to (8) being true at the transition point. When two equal energies are on the same side of the origin, both states are the same distance away from the minimum meaning the vector lengths are equal and are expressed as (9). Solving either equation to calculates α T between states with indices (k, l) and (m, n).
For (8), Fig. 4 indicates the signs that should be used to find α T . Using the incorrect signs will yield d with the opposite sign, but the correct vector length. These equations are transcendental and can be numerically solved. It should be noted that the weak potential dependence means that these values are inexact so it is appropriate to truncate the values to a few significant figures. Figure  2 shows a slight variation in the location of some of these points as the potential is changed.
We recover the statement of equal vector lengths with (9), which is a special case of the parabola model. We calculate all of the exact values of α T for this case in Appendix B.
The two-dimensional phase diagram can be numerically generated by performing the previous routine of sampling the energies of the lattice parametrised by c and ρ 1 . Results show that transition lines behave as expected: c = α T /ρ 1 . These results are then checked by finding the transition lines with the developed model. Figure 5 shows the phase diagram for 0 < c < 5 and 0 < ρ 1 < 5. Initial Monte Carlo simulations, sampling energies, and analytic methods all conform to the same picture of the phase diagram being symmetric under interchange of circumference and density. Since at any value of α △ , the corresponding state must be the ground state, we can link together these minima by searching for the transition point that must occur between them. Multiple transitions can be found between minima, so all relevant states near the minima must be considered. This can be done with confidence computationally.
We initially found ground states up to α = 25 and correctly predicted ground states and transitions up to α = 50.
Similar helical structures are observed [17,18] in cylindrically confined systems with a general trend of increasing row numbers as one moves along the phase diagram.
We note that row transitions are a general property of confined particles. Since this system is not thermal at zero temperature, we cannot state the true order of the structural transitions. By calculating the energetic derivatives and searching for discontinuities, we find that all the transitions we observe 'appear' first-order. This is typical of confined systems [14], with the only secondorder transition being the zig-zag transition between one and two chains, [19,20], which is absent in our system since the single chain is always unstable due to lack of global confinement.
In conclusion, we have developed a model which predicts the zero temperature ground states of identical repulsive particles confined to a cylindrical system as a function of geometry and density. Lead by a geometrical picture which emerges from the initial data, we infer an idealised description of the system and search for ground state transitions. The parabolic behaviour of the perparticle energy when measured in the reference frame of the lattice vector length allows us to write down a pair of simultaneous equations that solve for α T . We find the number of rows of particles generally increases with circumference or density. The occurrences of the perfect equilateral triangular lattice divides the phase diagram into sectors that we search between in order to find transitions.
The analysis relies on the scale invariant α which alone determines the ground state structure on the cylinder. This result is robust and successfully predicts the lattice structure for parameters originally out of scope of the initial phase space used to inform it. More generally, we can determine the entire phase space for ground states in this system.
Although our analysis takes α to be a global invariant, since we found these lattices forming on finite size cylinders, we predict that these structures can also form in localised regions. Taking a local value of α to depend on local density and circumference, where the system could have varying circumference and density, local ground state structures might form. This result would then have a much wider application into more general systems, such as conical geometries.