Transitions between phyllotactic lattice states in curved geometries

Phyllotaxis, the regular arrangement of leaves or other lateral organs in plants including pineapples, sunflowers and some cacti, has attracted scientific interest for centuries. More recently there has been interest in phyllotaxis within physical systems, especially for cylindrical geometry. In this letter, we expand from a cylindrical geometry and investigate transitions between phyllotactic states of soft vortex matter confined to a conical frustum. We show that the ground states of this system are consistent with previous results for cylindrical confinement and discuss the resulting defect structures at the transitions. We then eliminate these defects from the system by introducing a density gradient to create a configuration in a single state. The nature of the density gradient limits this approach to a small parameter range on the conical system. We therefore seek a new surface, the horn, for which a defect-free state can be maintained for a larger range of parameters.

A α values at transitions Table S1 states the transition values α T at which the lattice state transitions into the given phyllotactic state as the value of α is increased from α = 0. The transition values are calculated for the case of cylindrical confinement. This  Table S1. Transition values α T for which the system adopts the given phyllotactic state as the value of α is increased.

B Comparison between molecular dynamics and Monte Carlo results
In this section we give a direct comparison between results obtained for ground states using molecular dynamics (MD) and Monte Carlo (MC) methods. We use as our example the transition between the (4, 3, 1) and (4, 4, 0) states, which is expected to occur at α = 16.73. Figures S1(a) and S2(a) show snapshots of the ground states obtained using MD and MC respectively. In each snapshot there are 130 vortices and the system dimensions are L = 50, c 0 = 6.0 and ∆c = 0.5. In Figs. S1(b) and S2(b) we plot the density variation within each system, while in Figs. S1(c) and S2(c) we plot the resulting variation in the parameter α. The dark grey dashed lines in each graph correspond to the position of the transition within the system, which is marked by the presence of a single dislocation in the system snapshots. The light grey lines on each plot represent the expected α value of the transition and the corresponding position within the system at which the local value of α takes on that value, which is calculated from the red best-fit curve. In both cases, the difference between the expected and actual positions of the transition within the system corresponds to less than 4% of the system length. This difference is the length of the average lattice parameter in the snapshots.

3/6 C Molecular dynamics on the horn
We parameterise the horn surface using cylindrical coordinates (r, θ , z), with positions on the surface described by r = (r(z) cos θ , r(z) sin θ , z). The surface has a profile is the radius at z = 0 and z is the fractional distance along the length of the surface, such that 0 < z < 1, and η > 0. Throughout this section subscripts z or θ indicate a derivative with respect to that variable, while subscripts i and j are used to index particular vortices. It is useful to note the following relations In order to do molecular dynamics, we need to know the geodesics on the surface between pairs of points. The distance d i j between a pair of vortices is the length of the geodesic connecting the pair along the surface, a derivation for which is given in the following section. We numerically solve for the geodesics and calculate d i j for each pair within the cut-off radius.
The net force on a chosen vortex in a given time step is the sum of the individual forces acting upon it. We set a distance length scale a 0 as the length of the lattice parameter for the expected state at z = 0.5, i.e. f i j = − f 0 K 1 (d i j /a 0 ), where f i j is the magnitude of the force on vortex i due to vortex j, f 0 = 1 is a constant and K 1 is a modified Bessel function of the second kind. The distance δ d moved by a vortex is restricted to a small value through the choice of the value of time step δt. We therefore approximate that the vortex moves a distance δ d along the tangent vector to the net force, which is approximately parallel to the geodesic along the direction of the net force over the distance δ d.

C.1 Determining the distance between two particles -boundary value problem
The distance between two points is determined by numerically calculating the length of the geodesic connecting the pair. Depending on the state on the surface, there may be occasions where either θ z or z θ is singular. As such, geodesic equations for both θ (z) and z(θ ) must be known so that if one has a singularity, the other can be used.

4/6
Substituting in for r(z) then leads to the equation satisfied by the geodesic θ (z) 0 = 2(1 − ηz)(4(1 − ηz) 3 + η 2 r 2 0 )θ zz + 4ηr 2 0 (1 − ηz) 2 θ 3 z + η(8(1 − ηz) 3 − η 2 r 2 0 )θ z (C.12) Equation (C.12) can be solved numerically, with boundary conditions given by the positions of the two vortices, to find θ (z) and θ z (z). The distance d i j between vortices i and j is then the length of the geodesic and is determined by numerically integrating We repeat the process from the last section, this time solving for the geodesic z(θ ). In this case the line element ds is given by ds = dθ r(z(θ )) 2 + z 2 θ + r θ (z(θ )) 2 (C.14) Solving the Euler-Lagrange equation for z gives Simplifying this leads to the equation for the geodesic for z(θ ): Equation (C.20) can be solved numerically to find z(θ ) and z θ (θ ). The distance d i j between a vortex pair is once again the length of the geodesic and is determined by numerically integrating

D Determining the expected vortex locations on the horn
We derive the equation of the set of curves which can be used to construct the phyllotactic state with constant α on the curved surfaces. The vertices at which these curves intersect define the expected locations of vortex sites. The curves, known as loxodromes, always have their tangent vector at a fixed angle relative to the parallels and meridians of the surface. For the surface of revolution of the curve r(z) about the z axis, parameterised in cylindrical coordinates as r = (r(z) cos θ , r(z) sin θ , z), a set of orthonormal unit vectors along the surface can be defined aŝ where, in this section, subscript z denotes a derivative with respect to z. The parallels of the surface are the lines of constant v (or equivalently z) and the meridians are the lines of constant θ .