A renormalization approach to the universality of scaling in phyllotaxis

Physica D 298–299 (2015) 68–86

Contents lists available at ScienceDirect

Physica D journal homepage: www.elsevier.com/locate/physd

A renormalization approach to the universality of scaling in phyllotaxis Christian H. Reick Max Planck Institute for Meteorology, Bundesstraße 53, D-20146, Germany

highlights • Self similarity in Fourier transformed dynamics of Thornley’s model of phyllotaxis. • Construction of a renormalization theory for Thornley’s model of phyllotaxis. • Universality of scaling and divergence selection in Thornley’s model of phyllotaxis.

article

info

Article history: Received 11 August 2014 Received in revised form 7 January 2015 Accepted 9 February 2015 Available online 17 February 2015 Communicated by A. Mikhailov Keywords: Phyllotaxis Divergence spectrum Renormalization Universality Mathematical biology Critical phenomena

abstract Phyllotaxis, i.e. the arrangement of plant organs like leaves, florets, scales, bracts etc. around a shoot, stem, or cone, is often highly regular. Across the plant kingdom phyllotaxis shows not only qualitatively, but also quantitatively identical features, like the occurrence of divergence angles close to noble irrationals. In a previous study (Reick, 2012) a mechanism has been identified that explains the selection of these particular divergence angles on the basis of self-similarity and scaling, numerically found in the bifurcation diagrams of simple dynamical models of phyllataxis. In the present paper, by constructing a renormalization theory, the universality of this scaling is proved for a whole class of models, prototypically represented by Thornley’s model of phyllotaxis (Thornley, 1975). The renormalization is constructed from another self-similarity found numerically for the Fourier transform of the abstract potential governing the mutual inhibition of primordia. Surprisingly, the resulting renormalization transformation is already known from the treatment of the quasiperiodic transition to chaos but operates here on a different function space. It turns out that the fixed points of the renormalization transformation are characterized by divergences of the form Θ (κ) = 1/τ (κ) , where, written as continued fraction, τ (κ) = [κ; κ, κ, . . .], κ ∈ N+ . To show the universality of the scaling, it is demonstrated that the fixed points are unstable and that the associated scaling factors α (κ) = −(τ (κ) )2 and β (κ) = τ (κ) are exactly those that were numerically found in (Reick, 2012) to rule the selfsimilarity of the bifurcation structure. Thereby, the present paper puts forward an explanation for the universal appearance of certain phyllotactic patterns that is independent of physiological detail of plant growth. © 2015 Elsevier B.V. All rights reserved.

1. Introduction During the last years impressive advances in realistic modeling of the growth of phyllotactic patterns have been made [1,2]. By simulating explicitly the cell growth close to the tip of shoots these models reveal qualitatively and quantitatively the phyllotactic patterns observed in nature. But although these models are based on considerable physiological and biochemical insight, the mere occurrence of the correct patterns does not explain why only particular phyllotactic patterns occur in these simulations and

E-mail address: [email protected]. http://dx.doi.org/10.1016/j.physd.2015.02.003 0167-2789/© 2015 Elsevier B.V. All rights reserved.

throughout the plant kingdom in general [3]. For approaching this question the basic structural elements of phyllotactic growth that are responsible for the selection of the very few empirically observed patterns have to be identified. One such structural element is related to the observation Hofmeister made more than hundred years ago that new primordia – these are the early buds that on further growth form leaves, florets, bracts, etc. – occur where already existing primordia left the widest gap [4]. This could be easily understood if existing primordia would somehow inhibit the initiation of new primordia in their immediate neighborhood. This idea, published early in the 20th century by Schoute [5], is one such key structural element. Other key structural elements are circular symmetry and a superposition principle for the inhibitory ‘‘forces’’, as simple dynamical

C.H. Reick / Physica D 298–299 (2015) 68–86

models that implement these principles show by their ability to reproduce observed phyllotactic patterns [6–13]. The nature of these models is quite diverse, e.g. with respect to the representation of time (discrete vs. continuous), space (1D vs. 2D), or the conditions for the appearance of new primordia. It is interesting to note that also the above mentioned realistic tissue growth models [1,2] implement these basic elements: In the circular geometry of the shoot apex, they follow the idea of Reinhardt et al. [14] that the growth hormone auxin – at sufficiently high concentration responsible for the initiation of new primordia – is actively pumped against its gradient away from existing primordia (see also [15]). In this way, in their immediate neighborhood the growth of new primordia is inhibited. Another type of models has emerged from the assumption that primordia arise from a buckling instability of the surface tissues [16]. The analysis of the full nonlinear elastic problem close to the bucklung threshold reveals the full spectrum of observed phyllotactic patterns [17–19]. Basic to this mode coupling approach is besides a circular topology the scalar nature of the buckling field that leads to a generic form of the mode coupling equations close to the threshold. In this type of models, the appearance of an intrinsic length scale for the buckling prevents the appearance of new primordia close to existing primordia. It is especially interesting that by a mean field approximation of the realistic tissue growth model of Jönsson et al. [2] the mode coupling equations of the buckling instability approach are obtained [20]. Most striking in phyllotaxis is the re-appearance of certain number theoretic features in the spiral arrangements of plant organs. For example the numbers of parallel right and left spirals around e.g. a pine cone are typically successive members of a Fibonacci series. Another regularity refers to the divergence angle, which denotes the polar angle between two successive leaves or other plant organs measured around the central axis of e.g. the shoot. At large density of plant organs these angles are found to have values close to√ some noble irrational numbers, like the inverse golden mean ( 5 − 1)/2 = 0.61803 . . . ≈ 222.5°. During the last centuries this mystery has been the subject of numerous studies, and lots of explanations have been proposed [21], often concentrating on this golden divergence, although divergence angles close to other irrational numbers are well documented [3]. One particular challenge is to explain why the same numbers are observed in many different plant species and plant parts (shoots, flowers, fruits) despite a large diversity in growth form and structure across the plant kingdom. From a modeling point of view one thus has to show that the selection of certain numbers is largely independent of model details. It can thus be expected that simple models are much more promising to understand this universality than realistic ones that are typically too complicated to distinguish the relevant from the irrelevant structural elements. Several studies following this strategy have been published: Marzec and Kappraff [7] considered a variant of Thornley’s model [6]. Here each primordium is equipped with a function decreasing with distance that inhibits the initialization of primordia closeby. In their study Marzec and Kappraff could relate the predominance of the golden divergence to a sharp decrease of the Fourier components of this inhibitor function, other details are irrelevant. Using a model introduced by Douady and Couder [22], Kunz has proposed another explanation [23,24]. He conjectures that in the parameter region of sufficiently fast growth the bifurcation structure is largely independent of the type of functions used for inhibition of primordia growth, and exemplifies this for two different inhibitor functions. This bifurcation structure, closely following a Farey tree, selects divergence values close to the observed ones. In a similar direction goes the explanation by Guerreiro [11] (see also [25]): Using a Thornley type model he shows for a large class of inhibitor functions that solutions avoid the vicinity of rational divergence values with small denominators whereby noble

69

divergences are favored. Phyllotactic patterns also appear in infinite 2D lattices selected by minimizing a mutual repulsion energy between lattice points. For such systems Levitov [26,27] could relate the solutions of this minimization problem to the inherent self-similarity of the hyperbolic plane. The selected divergences are shown to be largely independent of the chosen repulsion forces. The cylindric geometry underlying phyllotaxis in plants is different from that of an infinite plane, so that Levitov’s theory is not directly applicable. Nevertheless, Atela et al. [28] formulated a model for plant phyllotaxis where the solution structure can as well be understood with help of the modular group. In the present study another approach to the problem of the universal re-appearance of certain phyllotactic patterns is pursued, starting from self-similarity and scaling. In a foregoing study [29] – called henceforth ‘‘study I’’ – scaling properties have been derived for two different types of models, namely Douady–Couder and Thornley type models (as they were called in that study). The following is a short summary of results from study I as far as they are relevant here. It focuses on Thornley type models because only these are considered in the present study; a more detailed description of this class of models follows in Section 2. As for all ‘‘field models’’ models of phyllotaxis [29] – some have already been mentioned above [7,11,22,28] – also in Thornley type models solutions are determined from a potential V (x), where x is an angular location around the shoot displaying phyllotaxis. This potential is constructed by superposing the inhibitory influences from all existing primordia. Concentrating on solutions forming a regular lattice of primordia, the potential depends on two parameters: the divergence angle Θ , and a control parameter ϵ that represents the intensity by which existing primordia inhibit the birth of new primordia. Although each pair of these parameters (Θ , ϵ) defines a lattice, only certain combinations of these parameters form solutions of the model, namely those for which the inhibitory forces are smallest, i.e. where the potential assumes a minimum: dVΘ ,ϵ (x)  dx

  

= 0.

(1.1)

x =0

As an example in Fig. 1 part of the divergence spectrum of one such model of Thornley type is plotted. Besides solutions to the latter equation (‘‘proper solutions’’) also parameter pairs are shown (‘‘improper solutions’’) where upon parameter variation the left hand side of (1.1) crosses zero discontinuously. These improper solutions form a kind of backbone structure, limiting possible behavior of proper solutions [23,29]. In study I the solution structure of equation (1.1) was analyzed for Thornley type models in the limit of vanishing growth ϵ → 0, which is the parameter region of high plant organ density. In this limit, called criticality in the following, the divergence spectrum shows self-similarity and scaling. This gets obvious by considering the neighborhood of points (Θ∞ , ϵ = 0) at criticality. Writing neighboring points as (Θ , ϵ) = (Θ∞ + ϑ, ϵ) it was numerically demonstrated in study I that the solution structure in such a neighborhood is to high accuracy self-similar under the rescaling

(ϑ, ϵ) → (αϑ, βϵ),

(1.2)

where α and β are the scaling factors, further investigated in the present study. Analyzed were some cases with irrational divergences having a continued fraction expansion1 of type Θ∞ = [a1 , a2 , . . . , ak , κ, κ, . . .]. Introducing

τ (κ) = [κ; κ, κ, κ, . . .]

(1.3)

1 Notation used for continued fractions: [a ; a , a , . . .] := a + 1/a + 1/a +· · ·, 0 1 2 0 1 2 with a0 , a1 , a2 , . . . ∈ N. For a0 = 0 the shorthand [a1 , a2 , . . .] is used.

70

C.H. Reick / Physica D 298–299 (2015) 68–86

Fig. 1. Selfsimilarity and scaling of the divergence spectrum in the vicinity of the divergence Θ∞ = 1/(τ (1) )2 = [2, 1, 1, 1, . . .] = 0.381966 . . . . The figures show solutions for a generalized Thornley model (see Section 2) using a bicubic inhibitor function (see Appendix A) and a non-monotonic sequence of aging coefficients an = n(1 − ϵ)n−1 . Subfigure (b) is a magnification of region b in subfigure (a) under the re-scaling (1.2) with scaling factors (1.4) for κ = 1. Similarly, using the same re-scaling, subfigure (c) is a magnification of region c in subfigure (b), etc. The small arrows in the middle of the Θ -axis indicate the location of Θ∞ , which is for this system for ϵ → 0 the limit value of the so called fundamental solution born at large ϵ in a symmetry breaking bifurcation and entering the graphs from the top. The solutions plotted in the graph are obtained by determining the parameter pairs (Θ , ϵ) for which dVΘ ,ϵ /dx|x=0 crosses zero (compare (1.1)). Continuous crossings reveal ‘‘proper solutions’’ of the model, whereas for ‘‘improper solutions’’ zero is crossed discontinuously. Proper and improper solutions are born as pairs in bifurcations seen as cusps in the divergence spectrum. In the graphs, improper solutions are visible as straight vertical lines at constant rational divergences. Accordingly, bifurcations can be identified with the rational divergence of the respective improper solution, as indicated by the fractional values seen in the graphs. More details on such divergence spectra can be found in study I.

it turned out that within numerical accuracy for the cases studied (κ = 1 and κ = 2)

 2 α = − τ (κ)

(1.4)

β = τ (κ) .

This self-similarity under re-scaling with α and β is demonstrated here once more in the sequence of plots shown in Fig. 1 for one particular Thornley type model. The value of α can be understood from the Farey tree structure of the divergence spectrum (see study I), whereas evidence for the value of β is purely numerical.

2

Surprisingly, for Douady–Couder models β = τ (κ) was found; why this value differs from the above one obtained for Thornley type models remains a puzzle, but is not topic of the present study. Instead, the central topic of the present study is the derivation of the scaling factors α and β of (1.4) for a whole class of Thornley type models thereby proving universality of scaling for this type of models. As a consequence, the present study provides further theoretical support for the universal selection of asymptotic divergence values by the ‘‘early determination mechanism’’ described in study I. The universality of scaling is demonstrated by employing a renormalization approach that is encouraged by the observed self-similarity in the divergence spectra. The development of



renormalization theory has been a major breakthrough in understanding universality of second order phase transitions [30–32] and also of certain transitions to chaos [33–36]. Basic to this theory is a renormalization transformation that expresses a self-similarity present in the considered class of systems. Close to criticality the scaling observed in these systems is related to the stability properties of the fixed points of the renormalization transformation. The existence of a renormalization implies that all systems close to the stable manifold of a particular fixed point show the same scaling behavior, i.e. they belong to the same universality class. Accordingly, to prove universality by means of a renormalization theory one has (i) to identify the renormalization transformation, (ii) determine its fixed points, and (iii) analyze the stability properties of the fixed points. In the present study this program is applied to Thornley type models of phyllotaxis. After introducing this class of models in Section 2, the derivation of the renormalization transformation follows in Section 3. This derivation is based on a self-similarity different from that found in study I, namely a self-similarity in the Fourier transform of the inhibitory potential. This self-similarity is the same as found in the spectra of quasiperiodic systems at the edge of the transition to chaos. The renormalization transformation turns out to factorize into three separate operators, so that the analysis of its fixed points and their stability reduces to a fixed

C.H. Reick / Physica D 298–299 (2015) 68–86

point analysis for these sub-operators; these sub-operators are analyzed separately in Sections 4–6. Main result of this stability analysis is the recovery of α and β from Eq. (1.4) as scaling factors controlling the dynamics of the renormalization into the unstable directions of the fixed points. In the final Section 7, the applicability of the foregoing results to phyllotaxis in general is discussed. – For better reading, mathematical detail has been shifted to the various appendices and specific mathematical notation is explained in footnotes. 2. Generalized Thornley’s model The class of models for which a renormalization is constructed in the present paper, generalizes Thornley’s model of phyllotaxis [6] with respect to the functional form of the aging term. Thornley’s model describes the growth of phyllotactic patterns in a highly idealized way. In plants new primordia (leaf buds) arise at the meristem that encircles the shoot apex close to its tip. This meristem is represented in Thornley’s model as a circle; positions on the circle are measured by an angular coordinate x ∈ [−1/2, 1/2); for technical reasons these coordinates are cyclically extended to x ∈ R. New primordia appear at discrete time steps. If the first primordium arises at time step n0 at angular position xn0 , subsequent primordia arise at positions xn1 , xn2 , . . . at time steps nk = n0 + k, k = 1, 2, . . . . To account for Hofmeister’s observation that the distance to existing primordia is decisive for the location of newly emerging primordia, in Thornley’s model every primordium is equipped with an inhibitory ‘‘potential’’ that suppresses the appearance of new primordia in its neighborhood. This potential decreases with angular distance and with the age of the primordium (number of time steps since it was born). The angular contribution to this potential is represented by an inhibitor function f (x) identical for all primordia, decreasing with angular distance x from the primordium. The age aspect of inhibition is modeled by prescribing a sequence of aging factors (an )n=1,2,... that reduce the inhibitory strength according to their age n: The older a primordium the less it is inhibiting. Finally, at time N, a new primordium is assumed to arise at angular position x, where the superposed inhibition from the existing primordia is least, or, more precisely, where the ‘‘potential’’ VN ( x ) =

N −1 

aN −n f (x − xn )

(2.1)

n =n 0

assumes its minimum. Most studies on phyllotaxis employing this type of model [6,7, 25,29] used aging factors of the form an = λn , λ ∈ [0, 1), only Guerreiro [11] presented some results for a more general class of aging factors. In the present study it is convenient to start out from sequences of aging factors that are weakly restricted by assuming an at least exponential decrease, i.e. for the sequences considered here there must exist a κ ∈ R and a γ with 0 < γ < 1 such that

|an | ≤ κγ n ,

n = 1, 2, 3, . . . .

(2.2)

This restriction is chosen for technical reasons to become clear in Section 5. The class of systems for which the results of the present study apply is likely larger, covering e.g. also sequences of the type an = n−(1+ϵ) , ϵ > 0. Condition (2.2) ∞implies in particular that the series an is summable, i.e. that | n=1 an | < ∞. Nevertheless, as will be discussed below, criticality, which is of particular interest here, is reached exactly where this summability breaks down. The class of inhibitor functions is made precise as follows: f (x) is twice continuously differentiable over (0, 1) with bounded second derivative, and f (x) = f (x + 1)

periodicity

(2.3)

f (x) = f (−x) d f (x) 2

dx2

71

symmetry

(2.4)

> 0 for x ∈ (0, 1) convexity.

(2.5)

Periodicity is assumed here only for technical reasons, namely to extend the domain of definition from [−1/2, +1/2) to R. In contrast, the symmetry of f (x) arises from the observation that there is no preferred angular direction around the meristem. The reasons for the assumed convexity is more subtle and will be discussed below. – Some immediate consequences of these assumptions are: the assumed convexity together with the symmetry implies that f (x) has a minimum at x = ±1/2, assumes its maximum at x = 0, is monotonous in between, and shows a cusp at x = 0 (compare Fig. 4 in Appendix E). In addition, a normalization is needed for f (x). Numerically convenient is f (0) = 1

⟨f (x)⟩ =

(2.6)



1 2

− 12

dx′ f (x′ ) = 0

(2.7)

but later on also other normalizations will be used. Note that the location of new primordia as determined from the minimum of (2.1) is independent from the choice of normalization. The striking regularities observed in spiral phyllotaxis show up at growth stages when already many primordia exist. These tend to be located at angular positions that differ by a constant angle, the divergence introduced above. Mathematically, these regular arrangements are obtained by assuming that the primordia are located at predefined positions xn =  x + nΘ , where Θ denotes the divergence and  x is an arbitrary angular offset that can be chosen at convenience. By letting n0 → −∞ and choosing  x such that the new primordium at time step N = 0 is located at x = 0, i.e. x0 = 0 solves Eq. (1.1), the potential becomes VΘ (x) =

∞ 

an f (x + nΘ ).

(2.8)

n =1

The whole study presented in the following is based on this ‘‘regular’’ potential. Note that the time index N of the potential seen in (2.1) has been replaced by Θ in (2.8) because the assumed regularity in the arrangement of primordia makes the regular potential independent of time but dependent on the assumed divergence angle. In study I the parameter λ in the exponential-type aging factors an = λn−1 played the role of a control parameter, since by increasing λ from zero to values close to one, the model can be driven into the critical region where it shows self-similarity and scaling. Analogously, here, families of aging factors (an,ϵ )n=1,2,... depending on a parameter ϵ > 0 will be considered, such that criticality is reached for ϵ → 0. For an = λn−1 one could e.g. choose ϵ = − ln(λ) or ϵ = 1 − λ. For such families, solutions of the model are determined by Eq. (1.1), which can be efficiently solved numerically (see e.g. study I). It is important to note that study I gives a hint on how to define criticality without referring to self-similarity: For the two rather different models of phyllotaxis considered there, in the undercritical situation primordia that appeared in the infinite past contribute only insignificantly to the potential (2.8), their contribution is effectively suppressed by a sufficiently fast decrease of the sequence of an . Criticality is reached when also primordia from the infinite past determine the location of new primordia. In the present context this means that criticality is related to an infinitely long memory in the sum (2.8) (i.e. λ → 1 for exponential aging factors). Mathematically, criticality is thus reached when the family (an,ϵ )n=1,2,... looses for ϵ → 0 its summability. Accordingly, the

72

C.H. Reick / Physica D 298–299 (2015) 68–86

restriction to a class of summable series of aging factors – like those introduced in (2.2) – is sufficient to study the approach to criticality numerically. But in the course of the theoretical analysis below, this class has to be extended to non-summable series that appear as limit values of families of such summable series in order to consider criticality explicitly. It remains to explain the reason for adding the convexity condition (2.5) on the inhibitory function. Considering the series of aging factors an = λn−1 Koch et al. [37] showed that the convexity condition guarantees that in the limit λ → 1 the proper solutions of Thornley’s model converge to regular arrangements of primordia associated with irrational divergences, whereas rational divergences, including a finite neighborhood of divergence values, are forbidden. This is in agreement with botanical observations where at high densities of plant organs, like e.g. in the capitulum of sunflowers, divergence values close to simple irrationals are found (see e.g. [3]). In the present study it will turn out, that the convexity condition is also essential for the much wider class of aging factors considered here (see Eq. (6.11) derived in Appendix B). 3. Renormalization of the potential The renormalization introduced below is motivated by a selfsimilarity found in the Fourier spectra of the regular potentials (2.8). This self-similarity is obtained as follows. The Fourier representation of the inhibitor functions reads because of symmetry (2.4) and normalization (2.7) f (x) =

∞ 

fn e2π inx = 2

n=−∞

∞ 

fn cos(2π nx)

(3.1)

n=1

with 1 2

 fn =

− 12

dxf (x) cos(2π nx).

(3.2)

∞ 

fn a(nΘ )e2π inx

(3.3)

n=−∞

where the aging function a(x) was introduced as a one-sided Fourier transform of the sequence of aging factors by a(x) =

In the following it is therefore sufficient to consider exclusively the odd part VΘ′′ (x) of the potential, and for shortness this part will itself be called ‘‘potential’’. Similarly, since with (3.8) also only the odd part of the aging function is relevant here, a′′ (x) will be called ‘‘aging function’’ without further mention of the symmetry. For the derivation of a renormalization transformation it is interesting to consider V ′′ (x) at criticality. Let Vn′′ = −2fn a′′ (nΘ ) denote the coefficients of the Fourier sine expansion (3.8). At criticality one finds for exponential aging factors an = λn−1 Vn′′ = −fn cot(π nΘ ).

Plugging (3.1) into (2.8) gives VΘ ( x ) =

Fig. 2. Coefficients Vn′′ of the cosine expansion of V ′′ (x) (see Eq. (3.10)) for golden divergence ΘG = [1, 1, 1, . . .] = 0.61803398 . . . at criticality. The coefficients are shown for three different inhibitor functions (biquadratic, bicubic and sine; see Appendix A). They are scaled by n/af to make the regularity and the asymptotic independence from the particular inhibitor function more visible. For the values af see Appendix A.

∞ 

an e2π inx .

(3.4)

n =1

It is convenient to split a(x) into its real and imaginary part: a(x) = a′ (x) + ia′′ (x),

(3.5)

with a′ (x) and a′′ (x) real for x ∈ R. Because of the definition (3.4), a′ (x) is an even and a′′ (x) an odd function of x. Using a similar notation for the even and odd parts of the potential, i.e. VΘ (x) = VΘ′ (x) + VΘ′′ (x),

(3.6)

one finds VΘ′ (x) = 2

∞ 

fn a′ (nΘ ) cos(2π nx)

(3.7)

n=1

VΘ′′ (x) = −2

∞ 

fn a′′ (nΘ ) sin(2π nx).

(3.8)

n=1

The solutions of Thornley’s model are completely independent from the even part of the potential because with (1.1) they are solutions of

  dVΘ′′    . = 0= dx x=0 dx x=0 dVΘ 

(3.9)

(3.10)

This result has been obtained by letting an → 1 ∀n ∈ N+ , which technically is achieved by performing the limit λ → 1 (compare Table 1). Fig. 2 shows (up to a scaling) these coefficients for three different inhibitor functions using as an example the golden divergence ΘG = [1, 1, 1, . . .]. First, it should be noted that for large n the Fourier coefficients get (up to factor) independent from the particular inhibitor function. This can be explained by the Dirichlet theory of Fourier sums (see e.g. [39, Vol. 2]): Because f (x) is symmetric and has a discontinuous derivative in x = 0 but is otherwise continuously differentiable, the Fourier coefficients fn decrease asymptotically as 1/n2 (see also Appendix B). Second, the Fourier coefficients displayed in this figure are obtained from the particular case of exponential aging coefficients in the limit λ → 1. But the result of this limit is largely independent of how it is reached, one could have taken also another set of aging factors reaching an = 1 in the critical limit. Therefore it is expected that the spectrum of Fourier coefficients at criticality is largely independent from the considered family of aging factors. Third, the spectrum shows a high degree of regularity, that improves for large n. E.g. there is a regular sequence of peaks at n = 3, 5, 8, 13, 21, . . . , which are successive Fibonacci numbers defined by Fn+1 = Fn + Fn−1 , F0 = F1 = 1. Because ΘG = limn→∞ Fn /Fn+1 , the structure of the peaks obviously reflects the chosen divergence. And indeed, although for other divergences the peak structure of critical spectra of V ′′ (x) is different, it shows similar number theoretic regularities (not shown). The regularity of the spectrum gets even more evident when plotting the asymptotic Fourier coefficients (3.10) in a somewhat different manner. For circular remapping of values x ∈ R to the interval [−0.5, 0.5) the abbreviation

(x)c := x − [x + 0.5]

(3.11)

C.H. Reick / Physica D 298–299 (2015) 68–86

73

Table 1 Some sequences of aging factors an , and the associated aging functionsa′′ (x). 3rdcolumn: behavior of the aging functions at criticality, n i.e. for λ → 1. 4th column: the associated fixed point  a′′ϵ (z ) = limn→∞ Raext a′′ϵ (z ) of Raext ; for shortness the normalization factor has been left away. Last column: type of singularity of limλ→1 a′′ (x) at x = 0, i.e. at criticality; note that this singularity type is identical to that of  a′′ϵ=0 (x) because limλ→1 a′′ (x) is member of the stable manifold of the fixed point  a′′ϵ=0 (x) of Ra . The singularity type reveals that only the first three series of aging factors lead to meromorphic aging functions at criticality. – Last line: ℑ() denotes the imaginary part of the argument, and Φ () denotes the Lerch function; the dash means that no simple analytic expression could be found. The entries in the last line have been obtained with help from [38, formulas 9.550 and 9.556]. an

a′′ (x)

limλ→1 a′′ (x)

λ

sin 2π x 1−2λ cos 2π x+λ2

1 2

n−1 n 1λ−λ2

(1−λ)2 sin 2π x (1−2λ cos 2π x+λ2 )2

cos π x 8 sin3 π x

n2 λn−1

1−6λ2 +λ4 +2λ(1+λ2 ) cos 2π x

n−1

λn−1 n

n−(2−λ)

(1−2λ cos 2π x+λ2 )3

1

λ

sin 2π x

λ sin 2π x arctan 1−λ cos 2π x

ℑ(yΦ (y, 2 − λ, 1)), y = e−2π ix

 a′′ϵ (x)

cot π x

x

Singularity

1

x−1

ϵ 1+( x )2 ϵ x

1

x−3

x

1

x−3

ϵ (1+( x )2 )2 ϵ

πx − 4cos sin3 π x

ϵ (1+( x )2 )2 ϵ

arctan(cot π x) π

sign(x) − π x 2

arctan(2π ϵx )

sign(x)

–

sign(x)

the mappings defining their dynamics, it is applied here to another object, namely the inhibitory potential underlying the phyllotactic dynamics. Therefore, despite their functional similarity, the renormalizations differ by the function space on which they operate. Mathematically, the renormalizations described by Ostlund et al. [35] and Shraiman [40] differ slightly in the definition of frequency bands. In the following the formulation by Shraiman, called removal of heads, is used. In this ‘removal of heads’ formulation step (i) from above is implemented by removing from (Vn′′ )n=1,2,... all Fourier components with frequencies larger than Θ /2, where Θ is an irrational divergence chosen for the particular construction, like ΘG in Fig. 3. The associated Fourier indices are collected in the set of ‘‘heads’’3



Fig. 3. Here the data for Vn′′ from the biquadratic case of Fig. 2 are plotted in a different way to bring out more clearly their regularity: abscissa values are the absolute value of the frequencies νn = (nΘG )c and the ordinate shows |Vn′′ | divided by νn . The arrows at the top indicate the abscissa values at midpoints between the major peaks.

HΘ

= hn ∈ N |n = 1, 2, . . . ; hn < hn+1 ; |(hn Θ )c | > +

(i) removal of the frequency band close to νn = 0.5, (ii) renumbering of all frequency components so that frequencies once more cover the whole range up to νn = 0.5, and (iii) rescaling of the amplitudes. Thereby the bands are shifted to higher frequencies and selfsimilarity means that the shifted bands approximately match the original bands. This kind of renormalization (as it is called) has previously been described in a different context, namely for the quasiperiodic transition to chaos, by Ostlund et al. [35] and Shraiman [40]. That the same renormalization applies to phyllotaxis is a consequence of identical band structures of the spectra (compare Fig. 3 of the present study with Fig. 11 in [35]). But whereas for quasiperiodic systems the renormalization is applied directly to

2 Here [x] denotes the Gauss bracket, i.e. the largest integer smaller or equal than x; confusion with the notation for continued fraction expansions is not expected.

2



.

(3.12) Consider its complement with respect to N+ :

Θ

 



is used.2 Using this circular remapping, Fig. 3 shows the absolute values of Vn′′ /(nΘG )c mapped against the modulus of the ‘‘frequency’’ νn = (nΘG )c . (Note that the mapping n → (nΘG )c is invertible because ΘG is irrational.) The similarity of the peak structure in the ‘‘bands’’ between major peaks that improves at small frequencies is quite striking. This kind of self-similarity is found also for phyllotactic patterns at other irrational divergences with periodic continuous fraction expansion (not shown). Intuitively, the self-similarity seen in Fig. 3 comes about by performing the following sequence of operations on (Vn′′ )n=1,2,... :

Θ

HΘ = N+ \ HΘ = n ∈ N+ |(nΘ )c | ≤ 2



.

(3.13)

Let the elements of HΘ be denoted by η1 , η2 , η3 , . . . such that ηn < ηn+1 , n ∈ N+ . Then one can show (see Appendix C) that



 



HΘ = ηn ηn =

n

Θ

+

1 2



 , n = 1, 2, 3, . . . .

(3.14)

Thereby the three transformation steps can collectively be summarized in the transformation

(Vn )n=1,2,... → ′′



1

N

Vη′′n

 n=1,2,...

,

ηn ∈ HΘ ,

(3.15)

where ηn assures the appropriate renumbering (step (ii)), and N is the necessary normalization factor requested by step (iii) (to be determined). Alternatively, by a glance at (3.8), one can define the renormalization as acting directly on V ′′ (x):

RVΘ′′ (x) := −





∞ 2 

N

fηn a′′ (ηn Θ ) sin(2π nx).

(3.16)

n=1

Note that the renormalization operator R introduced here acts independently of the system being critical or non-critical. This definition thus goes beyond the numerical evidence for self-similarity expressed in Fig. 3 for a system at criticality. But this more general definition is plausible insofar as it is expected that like for other

3 N+ = {1, 2, 3, . . .}.

74

C.H. Reick / Physica D 298–299 (2015) 68–86

critical phenomena self-similarity builds up gradually when approaching criticality. In the following it is shown that the renormalized potential (3.16) can be rewritten in a form structurally identical to the original potential (3.8), except that the inhibitor coefficients and the aging function are replaced by renormalized terms. The final expression is

RVΘ′′ (x) = −2





∞  [Rf f ]n [Ra a′′ ](nΘ ′ ) sin(2π nx)

(3.17)

n =1

where Θ ′ is a renormalized divergence 1

Θ ′ = G (Θ ) =

Θ

 −

1

Θ



.

(3.18)

Here G is the so called Gauss map (see e.g. [41]) and Rf and Ra are operators that will be specified below. The renormalized inhibitor coefficients showing up in (3.17) are obtained by comparison with (3.16):



Rf f

 n

=

fηn

Nf

,

ηn ∈ HΘ , n ∈ N+ ,

(3.19)

with a suitably chosen normalization Nf , e.g. such that the renormalized inhibitor function obeys (2.6). To obtain an expression for the renormalized aging function one can use the relation (proved in Appendix C)

(ηn Θ )c = −Θ (nΘ ′ )c .

(3.20)

With this relationship one finds a (ηn Θ ) = a (−Θ (nΘ )c ) = −a′′ (Θ (nΘ ′ )c ) (because a′′ (x) is odd). Taking this into account one can identify from (3.17) and (3.16) the renormalized aging function as ′′

Ra a′′ (x) =





a′′ (Θ x)

Na

,

′′

′

4. The renormalization G of the divergence In this section the divergence renormalization (3.18) (‘‘Gauss map’’) is analyzed for its fixed points and their stability. The Gauss map shows also up as part of the renormalization of the quasiperiodic transition to chaos [35,40]. There it describes the renormalization of the winding number, which is as well a possible interpretation of the divergence Θ to which the Gauss map applies here. The fixed points of the Gauss map are the divergences

Θ (κ) = [κ, κ, κ, . . .],

κ ∈ N+ ,

(4.1)

where [. . .] denotes a continued fraction expansion. This gets obvious by rewriting the Gauss map as

G ([a1 , a2 , a3 , . . .]) = [a2 , a3 , . . .].

(4.2)

(κ)

, G is expanded in their vicinity.

To analyze the stability of the Θ This gives

G(Θ (κ) + ϑ) = Θ (κ) + α (κ) ϑ + O (ϑ 2 )

(4.3)

where by noting that

τ (κ) = [κ; κ, κ, . . .] = 1/Θ (κ)

(4.4)

one obtains from (3.18) (3.21)

where (nΘ )c has been replaced by x. The normalization Na is related to the others by N = −Nf Na . Overall, by these considerations the renormalization R factorizes into the renormalizations Ra , Rf , and G, so that R can formally be written as (Ra , Rf , G). These three operators act on different spaces, the space of aging functions, the space of inhibitor coefficients, and the real space of divergences. Having determined the renormalization from the observed selfsimilarity, its fixed points and their stability have to be determined in order to prove the universality of scaling. This will take most of the remainder of this study. Here advantage will be taken of the factorization of the full renormalization into three separate operations: Obviously, the fixed points of the full renormalization can be determined by identifying the fixed points of the individual operators over their respective subspaces. Although not so obvious, also the stability analysis of the fixed points can be broken down to separate problems over the three subspaces: Ra and Rf are mutually independent, but each of them depends on Θ – Ra explicitly and Rf via the ηn – and therefore on the dynamics of the divergence under the action of G. Hence, the stability problem is not fully separating, so that it looks as if for each pair (Ra , G) and (Rf , G) their combined action has to be analyzed. But fortunately, the action of G itself is completely independent from the action of the other operators. Thereby, once the stability problem for the fixed points of G is solved, for the analysis of the stability of the fixed points of the combined problems (Ra , G) and (Rf , G), it will be sufficient to consider only the special case where Θ is a fixed point of G; this indeed reduces the combined problems to separate problems for Ra and Rf . This observation structures the next steps of the present investigation: in the following section ′

the fixed points of G are determined together with their stability. Then in the two subsequent sections the fixed point and stability problems for Ra and Rf are tackled, considering in accordance with the remarks from above only the special case where the divergence is a fixed point of G. Especially these latter two sections are rather technical. A reader interested only in the results may jump immediately to the final Section 7 where the findings are summarized and combined to draw the final conclusions.

α (κ) =

  dΘ  dG 

Θ (κ)

 (κ) 2 1 = − . 2 = − τ Θ (κ)

(4.5)

This value for α (κ) is exactly the value found for scaling of the divergence in study I (compare (1.2) and (1.4)). Moreover, since τ (κ) > 1, all fixed points of G are unstable. From the viewpoint of the full renormalization, α (κ) is an eigenvalue of the linearization of (Ra , Rf , G) at the respective fixed point, describing the unstable dynamics into the direction of divergence changes. Accordingly, all fixed points of the full renormalization are unstable with respect to perturbations of the divergence. This recovery of the scaling factors α (κ) for the divergences from the renormalization is the first major result of this study. It remains to recover also the scaling factor β for the control parameter rescaling (compare (1.2)). This will be the result of the next section. 5. The renormalization Ra of the aging function In this section the renormalization Ra of the aging function (3.21) is analyzed. The divergence Θ explicitly shows up in this renormalization, and is mapped as part of the full renormalization (3.16) to a new divergence value via the Gauss map G each time the renormalization is applied (see (3.18)). But as discussed at the end of Section 3, for the following analysis of the fixed points of Ra and their stability it is sufficient to let Θ be one of the fixed points Θ (κ) of G, so that this remapping of the divergence can be ignored here. To simplify the notation, these fixed point divergences Θ (κ) are simply denoted by Θ throughout this section. The analysis proceeds in several steps. First the necessity for a non-standard renormalization approach is discussed. Next, the special class of meromorphic aging functions is considered and an

C.H. Reick / Physica D 298–299 (2015) 68–86

‘‘extended’’ renormalization operator acting on families of aging functions is introduced. Then, by application of this extended renormalization, the control parameter scaling found numerically for Thornley’s model is recovered for the class of meromorphic aging functions. Since also non-meromorphic aging functions are relevant to phyllotaxis, the analysis proceeds with an analysis of the stability of fixed points belonging to other function classes. Finally, a simple plausibility argument is given, why despite many fixed points for every Θ only a single universality class is found. 5.1. Failure of the standard approach and a criterion for criticality According to the standard renormalization approach the natural way to proceed would be to determine the fixed points of Ra , linearize Ra around the fixed points, and determine the eigenvalues associated with the unstable eigendirections of this linearized operator. But for the case at hand, this standard procedure fails so that the problem has to be tackled in a different way. To understand why this standard approach fails, it is useful to have a glance at some of the fixed points of Ra . Choosing the normalization Na in (3.21) such that for a given value x = x0 the aging functions obey a′′ (x0 ) = 1, the renormalization becomes

  a′′ (Θ x) . a′′ (x) = Ra a′′ (x) = ′′ a ( Θ x0 )

(5.1)

This equation is a variant of Schröder’s equation f (g (x)) = af (x) that shows up also in other renormalization problems [42, Section 23], [43]. One easily sees that the functions

 a′′k (x) =

 x 2k+1 0

x

,

k∈Z

(5.2)

are fixed points of Ra , where the powers cover only odd numbers because the aging functions a′′ (x) must have an odd symmetry. Actually, there exist more odd fixed points of Ra , these will be considered in a separate subsection below. From the fixed points (5.2) only those representing critical systems are relevant here. As discussed at the end of Section 2, criticality is signaled by a loss of summability of the sequence of aging factors. By invoking the inverse of the Fourier expansion (3.4) this infinite sum can be represented in terms of a′′ (x) as ∞  n=1

1 2

 an = 2

dx a′′ (x) cot π x.

(5.3)

0

in the present context the behavior of Ra into the direction of subcritical systems is of particular interest. But into such directions a linear expansion of Ra has no meaning because a manifold of analytic aging functions (the sub-critical systems) cannot be continuously expanded at a singular aging function (representing a critical system). Accordingly, the standard renormalization procedure is not applicable here. And it will turn out that this conclusion remains also valid for other fixed points of Ra than those from the class (5.2). In conclusion, the stability of the fixed points has to be determined by other means. This is achieved in the following by extending the renormalization operator Ra from operating on functions to operating on families of functions. In this way whole stable and unstable manifolds of the fixed points of Ra will get fixed points of the extended renormalization. This will then allow to study explicitly the stability behavior of Ra along these manifolds. 5.2. Meromorphic aging functions and an extended renormalization operator To proceed in this way, the space of aging functions must be extended to include besides subcritical systems also critical systems, namely the non-analytic fixed points (5.2). Since the singularities of those fixed points are poles, it is natural to assume the aging functions to be member of the space of meromorphic functions; these functions are analytic in the finite complex plane except for a finite number of poles. Nevertheless, by the example of the special family of aging functions λn−1 /n it is clear that this class is too narrow in the present context: although all subcritical aging functions are analytic, this family gets non-meromorphic at criticality (compare Table 1). But fortunately, the analysis of meromorphic aging functions will also help to investigate other function types. Meromorphic functions can be written as the quotient of two entire functions, i.e. functions analytic over the finite complex plane. Employing the Weierstrass factor theorem for entire functions, and remembering that a′′ (x) must be odd, aging functions can be written as4

       2 2 2 z z z , , , . . . p1 p2 p3       , a′′ (z ) = z N 2 2 2 z z z Qb⃗ , q2 , q3 , . . . q1 Pa⃗

z ∈ C, N ∈ Z odd.

Because of the singular behavior of the co-tangens in x = 0, for representing a system at criticality a′′ (x) must obey

 ′′   a (x)     x  → ∞ for x → 0.

75

(5.4)

Thus the fixed points (5.2) represent critical systems only for k ≥ 0. Obviously a′′ (x) is non-analytic for these critical systems. This is the first ingredient necessary to understand the failure of the standard renormalization approach. The second ingredient is related to the assumption (2.2) of an at least exponentially fast decrease of the series of aging factors. From the theory of Fourier series it is known that the associated aging functions a′′ (x) are analytic (see e.g. [44, Section 25]). Since by this assumption the series of aging factors are summable, they represent sub-critical systems. Accordingly, for sub-critical systems the associated aging functions a′′ (x) are analytic. These observations now lead to the following conclusion. Following the standard renormalization approach one would linearize Ra around those fixed points of (5.2) representing critical systems and determine the eigendirections and associated eigenvalues. Considering the neighborhood of these fixed points,

(5.5)

Here Pa⃗ () and Qb⃗ () are entire functions that depend on a count-

⃗ = (a1 , a2 , a3 , . . .), b⃗ = (b1 , b2 , b3 , . . .), able set of parameters a and ±p1 , ±p2 , ±p3 , . . . , and ±q1 , ±q2 , ±q3 , . . . are the locations ⃗ ⃗ and b, of the zeros of Pa⃗ () and Qb⃗ (), respectively, independent of a

4 The Weierstrass factor theorem [45] states that an entire function f (z ) can be  

written as f (z ) = f (0) exp(G(z ))

kn

1 k= 1 k

 k z fn

∞

n=1

(1 −

z fn

) exp(hn (z ))

mn

, where hn (z ) =

for some suitably chosen kn ∈ N, and mn is the multiplicity of the

zero fn of f (z ). Moreover, G(z ) =

∞

1 n=1 n

 n z fN

is an entire function, written

here for convenience such that z is normalized by a suitably chosen zero fN of f (z ). This representation is valid for entire functions whose infinite sequence of zeros 0 < |f1 | ≤ |f2 | ≤ |f3 | ≤ . . . escapes to infinity. Accounting for the additional possibility of a zero in z = 0 with multiplicity m0 one sees that f (z ) can formally be written as f (z ) = z m0 Fg1 ,g2 ,g3 ,... ( fz , fz , fz , . . .), where g1 , g2 , g3 , . . . are all the 1

2

3

parameters involved in the Weierstrass representation of f (z ), in particular its zeros. Accounting for possible symmetries as reflected in the pairwise appearance of zeros, this leads to the representation (5.5) of a′′ (z ). In case f (z ) has only a finite number of zeros this representation is still valid, in that case the infinite product in the above representation of f (z ) is replaced by a finite product.

76

C.H. Reick / Physica D 298–299 (2015) 68–86

with qk ̸= 0. Because of this latter condition the behavior of a′′ (z ) depends for z → 0 only on the exponent N. In view of the criticality condition (5.4) the aging functions (5.5) are thus sub-critical for N ≥ 1, which is thus assumed in the following. Next, families of aging functions a′′ϵ (x) are considered, where ϵ ≥ 0 is a control parameter such that for ϵ > 0 the aging functions a′′ϵ (x) represent sub-critical systems, while criticality develops in the limit ϵ → 0. For the representation (5.5) of aging functions thus all continuous parameters must depend continuously on ϵ , ⃗ = a⃗(ϵ), b⃗ = b⃗(ϵ), as well as pk = pk (ϵ), qk = qk (ϵ), i.e. a k = 1, 2, 3, . . . . Since for ϵ > 0 all family members shall by construction represent sub-critical systems, i.e. N ≥ 1, and since N as a discrete parameter cannot depend on ϵ , criticality must develop for ϵ → 0 by the convergence of some zeros of Qb⃗(ϵ) () towards z = 0. Because of the symmetry of a′′ϵ (z ) this must happen by at least one pair of zeros simultaneously, e.g. by the pair ±q1 (ϵ). Since for x real, a′′ϵ (x) must be real, also the complex conjugate pair ±q∗1 (ϵ) must form zeros of Qb⃗(ϵ) . So there is a quadruple of zeros approaching zero for ϵ → 0, except for the case when q1 (ϵ) is purely imaginary so that a pair of zeros situated on the imaginary axis approaches z = 0. Note that q1 (ϵ) cannot be real because then a′′ϵ (x) would be singular at x = q1 (ϵ) and thus not be sub-critical for ϵ > 0. – The development of a singularity of a′′ϵ (z ) in z = 0 via the convergence of some zeros of Qb⃗(ϵ) () towards z = 0 for ϵ → 0 will be the key for the following construction of the extended renormalization operator. For ease of presentation, the case of a quadruple of zeros causing criticality will be temporarily left aside. To construct the extended renormalization it is important to understand the behavior of the aging functions (5.5) under the repeated application of Ra : n applications give



Ran a′′ (z ) =





zN (n)

Na

 2  n 2  n 2 zΘ zΘ Pa⃗ , , , . . . p1 p2 p3  , (5.6)      2 2 2 z Θn z Θn z Θn , q2 , q3 ,... Qb⃗ q1 z Θn

where for convenience the whole normalization has been absorbed (n) into the factor Na . Two observations are important here. First, for n → ∞, Pa⃗ () and Qb⃗ () converge because of 0 < Θ < 1 to values Pa⃗ (0, 0, 0, . . .) and Qb⃗ (0, 0, 0. . . .), which are non-zero by   construction so that limn→∞ Ran a′′ (x) ∼ xN ; for N ≥ 1 these are non-critical fixed points of Ra . Accordingly, all non-critical systems converge under repeated application of Ra to non-critical n n p fixed points. Second, since z pΘ and z qΘ can also be written as z / Θkn k

k

and z / Θ n , the repeated application of Ra can be understood as a shift of zeros of Pa⃗ () and Qb⃗ () towards infinity. Therefore one can use the modulus of that zero of Qb⃗ () located closest to z = 0 as a measure for the ‘‘distance to criticality’’. This idea of a distance is consistent with the observation that in order to reach criticality at least one pair of zeros ±qk (ϵ) of Qb⃗ () must reach zero, i.e. the ‘‘distance’’ must vanish. Assume for a given family a′′ϵ (x) that it gets critical when the pair of zeros ±q1 (ϵ) touches z = 0, i.e. ±q1 (ϵ = 0) = 0. Expansion into ϵ gives q1 (ϵ) = q ϵ + O (ϵ 2 ) with expansion parameter q. Since for the case of a pair of zeros the value q1 (ϵ) is purely imaginary (see above), and because the scale of ϵ can be chosen at convenience this expansion can be written as qk

q1 (ϵ) = iϵ + O (ϵ 2 ).

(5.7)

The extended renormalization Raext operating on whole families of aging functions is now constructed such that under its application the distance to criticality remains unchanged for every member of the family. This is achieved by assuring that z /q1 remains

unchanged. In view of (5.6) this is guaranteed to first order in ϵ by defining in analogy to (5.1) a′′Θ ϵ (Θ z )

Raext a′′ϵ (z ) =





a′′Θ ϵ (Θ x0 )

.

(5.8)

The key here is the multiplication of not only z by Θ but also ϵ by Θ that transforms z /q1 (ϵ) = z /(iϵ + O (ϵ 2 )) into z Θ /(iΘ ϵ + O (ϵ 2 )) = z /(iϵ + O (ϵ 2 )). The above assumption that the zeros ±q1 (ϵ) are responsible for a′′ϵ (x) getting critical for ϵ → 0 implies that no other zeros of Pa⃗(ϵ) or Qb⃗(ϵ) converge to zero in this limit. Assuming that these other zeros converge for ϵ → 0 to a finite value in the complex plane it is easy to see that under repeated application of the combined transform (z , ϵ) → (Θ z , Θ ϵ) involved in the repeated application of Raext they move to infinity. Accordingly, if in the beginning q1 determined the distance to criticality, this remains true under repeated application of Raext , so that indeed this operator conserves the distance. It should be noted that to first order in ϵ there is no alternative way to assure this constancy in distance. 5.3. Unstable manifolds and scaling Next the fate of families of aging functions under infinitely many applications of Raext is studied. Upon convergence, such families are fixed points  a′′ϵ (x) of Raext . With the knowledge from above on the fate of zeros under the repeated application one finds from (5.6):

 a′′ϵ (z ) = lim



n→∞

Raext

n

a′′ϵ (z ) =



zN (∞)

Na

Pa⃗(0) (0, 0, 0, . . .) Qb⃗(0)

  z 2 iϵ

. , 0, 0, . . . (5.9)

Here, for better readability, the normalization involved in (5.8) has (∞) been condensed into the factor Na . Since ±iϵ are the only zeros ′′ of  aϵ (z ), one can write Qb⃗(0)

  2 z iϵ

   z 2 m1  z 2  , 0, 0, . . . = 1 − R , iϵ iϵ

(5.10)

where m1 is the multiplicity of the pair of zeros, and R(z ) is an entire function without zeros in the finite complex plane. Such a function R(z ) can always be written as R(z ) = c exp(−G(z )), with a constant c and an entire function G(z ) (for a proof see [45, Section 7.6]). Thereby (5.9) gives when accounting explicitly for the normalization introduced in (5.8)

 a′′ϵ (x)

 =

N 

x

1+

 x0 2 m1 ϵ

 2

x0

1 + ϵx

  2   2  x x0 × exp G −G . ϵ ϵ

(5.11)

As far as the considered aging function families stay finite in the limit ϵ → 0 the exponential factor must converge towards 1 in this limit so that

 a′′ϵ=0 (x) =



x

N −2m1

x0

(5.12)

whereby the fixed points (5.2) of Ra are recovered. In particular, because by construction N ≥ 1, all fixed points (5.12) possess poles in x = 0 and thus represent critical systems as intended by the construction. From (5.11) now the main conclusion of this section follows. Application of the standard renormalization operator Ra to the family  a′′ϵ (x) reveals that



Ra a′′ϵ (x) =  a′′ϵ/Θ (x).



(5.13)

C.H. Reick / Physica D 298–299 (2015) 68–86

This means that Ra maps a family member at parameter ϵ to another family member at ϵ/Θ > ϵ and hence away from its fixed point  a′′ϵ=0 (x). Accordingly, the fixed points (5.11) of Raext are the unstable manifolds of the fixed points (5.2) of Ra , and along these unstable manifolds the control parameter ϵ scales as ϵ → ϵ/Θ . Thereby with (4.4) the second scaling parameter

β (κ) =

1

Θ (κ)

= τ (κ)

(5.14)

found empirically in study I (compare (1.4)) is recovered; here the shorthand Θ used for the fixed points of G throughout this section has been reverted to Θ (κ) . This same result is obtained when a quadruple of zeros is assumed to survive the limit n → ∞ performed in (5.9). In this case, instead of (5.7), one can set q1 (ϵ) = exp(−iϕ)ϵ + O (ϵ 2 ), where the angle ϕ defines the directions in the complex plane by which the quadruple of zeros approaches z = 0. As a consequence the factor [1 − z 2 /(iϵ)2 ]m1 in (5.10) is replaced by the product [1 − z 2 / exp(2iϕϵ)]m1 [1 − z 2 / exp(−2iϕϵ)]m1 , leading in (5.11) to a quartic polynomial x4 + ϵ 4 − 2x2 ϵ 2 cos(2ϕ) instead of a quadratic polynomial. This leads as well to the scaling (5.14). 5.4. Other types of fixed points Table 1 lists two examples of series of aging factors leading upon criticality to non-meromorphic aging functions. Accordingly, the above discussion of meromorphic aging functions is incomplete. And indeed, there exist more fixed points of Ra than only those listed in (5.2). These additional fixed points are not covered by the above discussion of meromorphic aging functions. Therefore in the following it will be discussed whether with these additional fixed points other types of scaling might appear. The fixed point equation (5.1) is a version of Schröder’s equation whose solutions have been identified for positive x by Aczél and Kuzma [46]. Extending their results by symmetry of the aging functions to negative x the full set of solutions is

 a′′ (x) = σ± (x)  with

 x k p(ln |x|/ ln Θ ) 0 , x p(ln x0 / ln Θ )

k odd k even

for σ± = σ+ for σ± = σ− ,

(5.15)

where σ± (x) is either the even function σ+ (x) = 1, or the odd function σ− (x) = sign(x), and p(x) is an arbitrary function with period 1, i.e. p(x + 1) = p(x); in addition a normalization with x0 > 0 has been assumed. From (5.4) it is seen that these fixed points can represent critical systems only for k ≥ 0, presumed that p(x) is bounded because otherwise  a′′ (x) would have infinitely many singularities so that criticality would loose its uniqueness. Eq. (5.15) covers all the meromorphic fixed points (5.2). The other fixed points of (5.15) fall into three groups that are discussed separately in the following. The fixed point  a′′ (x) = sign(x) This fixed point represents the special case k = 0, p(x) = 1, in (5.15). From Table 1 it is seen that this type of step function singularity appears in the critical aging function of e.g. the family of aging factors λn−1 /n. Following the same strategy as for meromorphic aging functions, this family can be used to identify the unstable manifold of this fixed point by applying the extended renormalization repeatedly—indeed, although Raext has been derived by analyzing meromorphic aging functions, there is no reason, why this operator should loose its applicability beyond this function class. Following this approach, one can easily check that – ignoring normalization – the function arctan(2π ϵx ) listed as limit limn→∞



Raext

n

a′′ϵ (z ) in the table is for Ra an unstable manifold



77

of the fixed point considered here. Hence, as in the meromorphic case, Eq. (5.13) holds and thus the families getting critical at the stable manifold of this fixed point also scale with β (κ) from (5.14). – It is interesting to note that also the family of aging factors n−(2−λ) (see table) develops at criticality a step singularity, although in contrast to the aging factors λn−1 /n its members are not from the class of exponentially bounded sequences (2.2), from which the present study started out. Note also that every family arctank (2π ϵx ), k odd, is an unstable manifold of the considered fixed point so that this unstable manifold is infinitely dimensional. Fixed points with p(x) ̸= 1 Because of the logarithmic argument of p() in (5.15) such fixed points involve an essential singularity in z = 0. It must be doubted that such fixed points have any relevance for phyllotaxis; the argument is as follows. From (5.1) it is seen that the particular structure of Ra implies that if a fixed point has the form of a product  a′′ (x) = f (x)g (x) then each of the functions f (x) and g (x) must be themselves fixed points of Ra (up to a normalization ignored here). Since p(ln |x|/ ln Θ ) shows up in the fixed points (5.15) as such a factor, this function must thus be itself a fixed point of Ra . Consider now an aging function a′′ (x) that belongs to the stable manifold of this fixed point and write it as a′′ (x) = p(ln |x|/ ln Θ )q(x). Upon these assumption it follows limn→∞ [(Ra )n q] (x) = 1. Conversely, any q(x) with this property makes a′′ (x) a member of the stable manifold of p(ln |x|/ ln Θ ). Thereby it is seen that any such a′′ (x) depends on Θ because p(ln |x|/ ln Θ ) does so. Accordingly, also a family of aging factors getting critical at a member a′′ (x) of the stable manifold of p(ln |x|/ ln Θ ) must encrypt Θ to hit a′′ (x) at criticality. Remember that Θ is here a shorthand for any of the Θ (κ) defined in (4.1), so that Θ is an element of a denumerable set of numbers. Therefore it is extremely improbable that upon a random choice of families of aging factors the associated critical aging function would indeed hit an a′′ (x) for the appropriate Θ . This argument considers Θ as a given number. But even relaxing this assumption does not change the picture: All Θ (κ) are unstable fixed points of the renormalization G. Accordingly, the dynamics of the full renormalization does not help in selecting such divergence values. Fixed points of the form sign(x)(x0 /x)k , k even For this class of fixed points unstable manifolds are found in analogy to (5.11) as

 a′′ϵ (x) =

1  x M sign(x)

N

ϵ

 m

1 + ϵx

(5.16)

with normalization N , M ≥ 0 even, and m > 0. An exponential factor as in (5.11) could be included but would not change the situation. One can easily check that (5.13) holds and thus once more the scaling (5.14) is recovered. It cannot be ruled out that other unstable manifolds exist that might show a different scaling. But generally it is currently unclear whether this class of fixed points is at all relevant to phyllotaxis. In principle one could derive from (5.16) the associated aging factor sequences to check their relevance, but so far the involved technical obstacles could not be overcome. 5.5. Many critical fixed points but only a single universality class The foregoing analysis revealed that fixed points of Ra either go along with the control parameter scaling found numerically for Thornley’s model (compare β (κ) in (1.4) and (5.14)), or the fixed points must be considered as irrelevant in the present context (the case p(x) ̸= 1 in (5.15)). The usual expectation from other renormalization theories would have been that each fixed point

78

C.H. Reick / Physica D 298–299 (2015) 68–86

defines its own universality class. A universality class covers all systems getting critical at the stable manifolds of the same fixed point. This concept of a universality class reflects the expectation that each fixed point goes along with its particular scaling obtained from the behavior along its unstable manifolds. But surprisingly, in the present case scaling turned out to be identical for all relevant fixed points so that in the sense of scaling only a single universality class is found despite the multitude of fixed points. That this is a consequence of the particular structure of the renormalization operator is elucidated by the following plausibility argument. Ra as defined in (5.1) involves the rescaling of x by Θ . Along the unstable manifold of a fixed point, Ra must map each manifold member to another member. Consider a family in the unstable manifold parametrized by a single control parameter ϵ such that for ϵ = 0 the fixed point is obtained; this family is denoted here as a′′ (ϵ, x). Then the application of Ra to a member of this family must be equivalent to a mapping of ϵ from one value to another, and this must compensate for the scaling of x. This is most simply achieved by assuming that x and ϵ show up in a′′ (ϵ, x) only in the combination x/ϵ , i.e. a′′ (ϵ, x) can be written as a′′ (x/ϵ). This implies by arguments already used above the scaling factor β (κ) from (5.14), independently of the fixed point considered. These considerations reveal that the uniform control parameter scaling across all fixed points is an immediate consequence of the fact that scaling of x shows up explicitly in the definition of Ra . This is different for renormalization operators in other contexts where the scaling is encoded implicitly in the nonlinear structure of the renormalization operator. 6. The renormalization Rf of the inhibitor function In this section the renormalization Rf from Eq. (3.19) acting on the space of inhibitor coefficients (3.2) is analyzed. It will be shown that the fixed points of Rf have for the considered class of inhibitor functions (2.3)–(2.5) no unstable manifold. For the discussion of universality of phyllotaxis not the fixed points themselves, but only their stability is relevant. Therefore the rather technical construction of the fixed points is shifted to Appendix E. The main result from this construction is that for every fixed point of the Gauss map, i.e. for every divergence of type (4.1), a unique fixed point of Rf exists. Examples for such fixed points are shown in Fig. 4 of the Appendix. With this result at hand, the discussion in this section concentrates on the stability of these fixed points. The stability is discussed only with respect to perturbations in the space of the inhibitor coefficients, the additional instability arising from the dynamics of the divergences under the Gauss map, need not be considered, as discussed at the end of Section 3. Accordingly, the divergences in this section are always the fixed point divergences Θ (κ) of G (see (4.1)), for which once more collectively the shorthand Θ will be used. The renormalization (3.19) investigated here, contains the normalization Nf , which was so far chosen such that f (0) = 1 (compare (2.6)). But for the following considerations it is more convenient to use instead the normalization fh1 := 1

(6.1)

where the index h1 is the smallest number of the set of heads (3.12). This choice can always be made because, as shown in Appendix B, for the class of inhibitor functions considered here, their Fourier coefficients are strictly positive. With this normalization the renormalization (3.19) reads

It is convenient to translate the renormalization (6.2) into an equivalent linear operation. Let the series ( fj )j=1,2,... denote a fixed point of the renormalization (6.2). Such fixed points exist (see Appendix E). Because the inhibitor coefficients fj are positive one can define for every fj a related coefficient ϕj via fj

ϕj := ln .  fj

Hereby (6.2) induces a related renormalization Rϕ on the series of coefficients (ϕj )j=1,2,3... :

 

Rϕ ϕ



Rf f

 j

=

fηj fηh

,

j ∈ N+ ,

(6.2)

1

i.e. Nf = fηh , where η1 , η2 , η3 , . . . are the elements of HΘ (see 1 (3.14)).



Rf f

= ln

j

j

 fj

,

(6.4)

which by entering into (6.2) can be identified as



Rϕ ϕ

= ϕηj − ϕηh1 ,

 j

j ∈ N+ .

(6.5)

Note that the normalization (6.1) translates into ϕh1 = 0 and that by (6.5) and in agreement with (6.3) the null series ( ϕj = 0)j=1,2,... is a unique fixed point of Rϕ . Via (6.4) the dynamics under the action of Rf is tightly linked to the dynamics under Rϕ . Thereby in particular the stability of the respective fixed points is linked: because of (6.3) a fixed point of Rf is stable/unstable exactly when the associated fixed point of Rϕ is stable/unstable. And one can even show that the stability is determined by the same eigenvalues in both cases. Accordingly, the following analysis of the stability of the fixed points of Rϕ immediately informs about the stability of the fixed points of Rf . The advantage of using Rϕ is, that in contrast to Rf , this operator is linear. Let the series (Φj )j=1,2,... be an eigenseries of Rϕ with eigenvalue δ , i.e.



Rϕ Φ

 j

= δ Φj ,

j = 1, 2, . . . .

(6.6)

With (6.5) this gives

δ Φj = Φηj − Φηh1 .

(6.7)

To solve this equation it is convenient to switch to the ‘‘train’’ notation (described in more detail Appendix C): By the definition of the set of heads (3.12) each irrational divergence – in the present context these are the fixed point divergences of G – introduces a particular arrangement of all natural numbers into separate sequences of numbers called trains; an example is shown in Table 2 of the Ap(n) pendix. The n-th train, n = 1, 2, . . . , is denoted by (tk )k=0,1,2,... , (n)

where each tk

is a natural number. In particular the n-th head is (n)

(1)

the zeroth element of the n-th train, i.e. t0 = hn , and t1 = ηh1 , where ηn is the n-th element of the complementary head set (3.13). Since every natural number shows up in exactly one train, there is (n) for each index j in (6.7) a train element tk = j. Using this, while (n)

noting that ηt (n) = tk+1 (see Appendix C), Eq. (6.7) gives k

δ Φt (n) = Φt (n) − Φt (1) . k

k+1

(6.8)

1

This recursion can be easily solved: For δ ̸= 1

Φt (n) = k

1 − δk 1−δ

Φt (1) + δ k Φhn ,

(6.9)

1

whereas for δ = 1 one finds

Φt (n) = kΦt (1) + Φhn . k



(6.3)

(6.10)

1

For the considered class of inhibitor functions it is demonstrated in Appendix B that in particular as a consequence of the convexity condition (2.5) the limit lim j2 fj

j→∞

(6.11)

C.H. Reick / Physica D 298–299 (2015) 68–86

exists and is non-zero. As shown in the following, this has consequences for the possible values of δ . Via (6.3) each eigenseries (ϕj = Φj )j=1,2,... is associated with a series of inhibitor coefficients fj . Applying the limit (6.11) to this series, and noting that this non-zero limit exists also for the series of fixed point coefficients ( fj )j=1,2,... , one sees that also the limit lim Φj

(6.12)

j→∞

(n)

exists. And this holds also for subseries. Hence, replacing j by tk , also the limit lim Φt (n)

k→∞

(6.13)

k

must exist. This limit is now applied to Eqs. (6.9) and (6.10). The first of these equations was derived under the assumption δ ̸= 1. A glance at this equation shows that limk→∞ Φt (n) exists k

only if |δ| < 1, because otherwise the series Φt (n) would diverge k

(if |δ| > 1), or the series would have two different accumulation points (if δ = −1). Applying the limit (6.13) to the second equation (6.10), which was derived for δ = 1, it can be seen that the limit exist only for Φt (1) = 0, which means that limk→∞ Φt (n) = Φhn . Because k

1

limj→∞ Φj exists, all subseries (Φt (n) )k=0,1,2,... must converge to the k

same value. Hence, Φh1 = Φh2 = Φh3 = · · ·. But via (6.3) the normalization (6.1) implies Φh1 = 0 so that Φhn = 0 for all n. Therefore, with (6.10), Φt (n) = 0 for all k and n. Hence for the k

considered class of inhibitor functions the only solution to (6.10) is (Φj )j=1,2,... = 0, which is not a valid eigenseries of Rϕ . In conclusion, the properties of the considered inhibitor functions imply |δ| < 1 so that the fixed point of Rϕ and thus also the fixed points of Rf are stable. Note that this is in particular a result of the convexity (2.5) of the inhibitor functions because, as seen in Appendix B, this condition guaranties the existence of the limit (6.11), which is basic to the foregoing stability proof. But it can also be seen that convexity is only a sufficient, but not a necessary condition: the proof would also work if limj→∞ jq fj would be non-zero with q ̸= 2, because also in this case the existence of the limit (6.12) would be guaranteed. 7. Conclusions The renormalization of Thornley’s model of phyllotaxis, introduced in Section 3 of this study, could be factorized into three separate operations: the remapping G of the divergences via the Gauss map (3.18), the renormalization Ra of the aging function (3.21), and the renormalization Rf of the inhibitor function (3.19). Each of these three operations describes separate degrees of freedom of the full renormalization, which formally thus can be written as (Ra , Rf , G), acting over the product space of aging functions, inhibitor functions, and divergences. Accordingly, fixed points of the full renormalization combine from the fixed points of the separate renormalizations. The analysis from the previous three sections revealed for G infinitely many fixed point divergences Θ (κ) , κ = 1, 2, . . . , for every Θ (κ) infinitely many fixed points of Ra , and for each Θ (κ) an otherwise unique fixed point of Rf . Altogether, the full renormalization thus possesses infinitely many fixed points for each Θ (κ) . Crucial for the question of universality is the stability of the fixed points. As already discussed at the end of Section 3, the stability of the full renormalization follows from the stability of the individual renormalizations. In Section 4 it was shown that all fixed points of G are unstable. Hence, all fixed points of the full

79

renormalization are unstable against perturbations of the divergence. Additional unstable directions arise from the renormalization of the aging functions (Section 5), whereas the fixed points are stable against perturbations of the inhibitor functions (Section 6). The key result of the foregoing renormalization analysis is the recovery of the scaling factors α (κ) and β (κ) (see (4.5) and (5.14)), that in study I were empirically obtained from an analysis of the self-similarity of the divergence spectrum of Thornley’s model. In the present study numerically identical values were obtained from a stability analysis of the renormalization. Although the starting point of the present study was a selfsimilarity in the Fourier spectrum of the inhibitory potential, which is different from the self-similarity analyzed in study I, this recovery of the scaling factors is not a mere numerical coincidence: Based on the self-similarity found in the Fourier spectrum of V ′′ (), the renormalization was introduced in Section 3 such that close to criticality it maps V ′′ () approximately to itself. This entails a rescaling of the divergence (compare (4.3)), and a rescaling of the control parameter ϵ (compare (5.13)). Writing Θ = Θ (κ) + ϑ , where Θ (κ) denotes as before a fixed point divergence of G, the renormalization thus maps the family VΘ′′ (κ) +ϑ,ϵ approximately

onto itself, with (ϑ, ϵ) → (α (κ) ϑ, β (κ) ϵ). Hence, if (Θ (κ) + ϑ, ϵ) is the location of a solution of Thornley’s model (compare Eq. (3.9)), another solution is found approximately for (Θ (κ) + α (κ) ϑ, β (κ) ϵ). This is exactly the self-similarity observed in study I for the divergence spectrum, summarized in the introduction (compare (1.2)), where a slightly different notation was used: Θ∞ ≡ Θ (κ) . Accordingly, the two types of self-similarity and scaling are closely related. This recovery of the scaling factors α (κ) and β (κ) from the unstable manifolds of the renormalization proves the universality of the scaling, because by standard arguments of renormalization theory (see e.g. [30,32,47]) every family approaching the stable manifold of a fixed point scales at criticality according to the behavior along the unstable directions of that very fixed point. In addition, because a multitude of fixed points has been found, the present renormalization theory predicts several universality classes: there is one universality class for each divergence Θ (κ) , because there is one fixed point of Ra for each Θ (κ) . And for each divergence Θ (κ) there are separate scaling factors α (κ) and β (κ) . – All these renormalization predictions are in full agreement with the results from study I (compare (1.4)). Surprisingly, the ‘‘momentum’’ variant of the renormalization developed for systems showing a quasiperiodic transition to chaos [35,40] was re-derived here as an expression of selfsimilarity for phyllotactic systems close to criticality. Both types of systems are similar in that the underlying dynamics maps the unit circle onto itself. But this is not sufficient to explain why for both systems the spectra show similar band structures, which is the reason why the same renormalization applies. To understand this better, it might be useful to re-formulate Thornley’s model more clearly as a particular type of circle map, which is the framework used to describe the quasiperiodic transition to chaos (see e.g. [48]). The scaling and self-similarity discussed so far applies only to regions of the divergence spectrum around the fixed point divergences Θ (κ) , because these are the divergence values at the fixed points of the renormalization. By the symmetry of Thornley’s model under reflection around x = 0, which implies identical values of the potential for divergences symmetric around 0.5, it is obvious that the results obtained so far apply also to the vicinity  (1) = 1 − Θ (1) = [2, 1, 1, . . .], or of the mirrored divergence Θ  (κ) = 1 − Θ (κ) = [1, κ − 1, κ, κ, . . .] for more generally around Θ κ > 1. This opens the question at which other divergences scaling

80

C.H. Reick / Physica D 298–299 (2015) 68–86

can be found. Simulations (unpublished) indicate that divergences equivalent to the fixed point divergences (for the notion of ‘‘equivalence’’ compare the Appendix of study I), i.e. divergences of type (κ) ΘK := [a1 , a2 , . . . , aK , κ, κ, . . .], one finds the same scaling as for Θ (κ) . From the viewpoint of renormalization this can be (κ) understood as follows: Because GK (ΘK ) = Θ (κ) the family of po′′ tentials V (κ) is mapped by K applications of the full renormalΘK ,ϵ1

ization R to a family of potentials VΘ′′ (κ) ,ϵ , where ϵ1 and ϵ2 are the 2

respective control parameters. For this latter family the renormalization theory of the present study applies. Hence, assuming that close to criticality the mapping between the two families can be linearized, the universal scaling of the latter transfers to the former. – Universal scaling can also be expected at divergences with periodic continued fraction expansion. E.g. for a divergence with period-2 expansion Θ = [κ1 , κ2 , κ1 , κ2 , . . .] two applications of the renormalization map Θ always back to itself. Therefore around such a divergence universal scaling with the scaling factors α (κ1 ) and β (κ1 ) raised to the square is expected. The generalization to divergences with continued fraction expansions of other periods is obvious. Thereby these considerations corroborate the conjecture from study I that scaling of the divergence spectrum is to be expected in the neighborhood of every divergence of the form [a1 , . . . , aN , κ1 , . . . , κP ], N and P finite, where the overlining denotes infinite repetition. Another possible generalization of the theory presented in this study pertains to multijugacy. In Thornley’s model a single primordium is added in each iteration step. This generates so called ‘‘unijugate’’ phyllotaxis. But in many plants two or more primordia are born at the same time. Already 100 years ago, van Iterson [49] showed in detail how these so called ‘‘multijugate’’ phyllotactic lattices are constructed from unijugate lattices: by squeezing the unijugate lattice in angular direction by a factor 1/K (where K is the number of primordia born simultaneously), and then repeating this squeezed lattice K times around the shoot. Similarly, Thornley’s model can easily be extended to cover multijugate phyllotaxis by adding in (2.1) inhibitory terms for primordia born at the positions 1/K , 2/K , . . . , (K − 1)/K giving (K )

VΘ ( x ) =

∞  n =1

an

K −1  

f

k=0

x + nΘ +

k K



.

(7.1)

After a bit of algebra one obtains an expression identical to the potential (3.3), except that Θ and x are scaled by the van Iterson factor 1/K and that the Fourier coefficients fk of the inhibitor function are replaced by KfKk . Hence, because the structure of the potential is the same as in the unijugate case, the renormalization results of this study apply also to the multijugate case. The present study, based on the class of Thornley models, involves a number of idealizations. One may wonder to what extent these restrict the relevance of the obtained results for botanically observed phyllotaxis. The renormalization theory presented here was developed for the idealized case of regular solutions, i.e. solutions with a single divergence across the phyllotactic lattice (compare (2.8)). In practice this means that the considerations are applicable whenever for successive primordia the divergence angle is approximately constant, which is often realized at least for finite sections of a specimen (see e.g. the empirical papers in [3]). Another idealization concerns the assumed convexity of the inhibitor functions which implies a sharp cusp at x = 0 (compare Section 2). But primordia are not points, hence there is a natural limit for their density, so that the situation of criticality, where the pointlike primordia are infinitely close and dense, is never reached by plants. But for the mechanism underlying Thornley’s model to be applicable to the realistic situation of finite bud size it is only important

that the newly forming primordia experience the convex parts of the inhibitory potentials of the already existing primordia, the cusp itself is irrelevant. One may speculate that for specimen showing malformations of phyllotactic patterns, like those discussed in [50], this convexity condition is violated. Related to the impossibility of reaching the critical line in nature, is the problem of the finiteness of the critical region for the applicability of a theory like the one presented in this study. The extension of the critical region, defined by the validity of scaling laws, may differ from system to system. This latter point seems at a first glance to pose a serious limitation to the applicability of the present results. But the renormalization picture entails that even outside the critical region the system ‘‘feels’’ its presence, in the sense that self-similarity reaches in a universal way down into the subcritical region, although with distance from criticality the scaling gets more and more distorted. This distortion develops in a systematic and universal way, quantitatively determined by the irrelevant eigenvalues of the linearized renormalization operator (see e.g. [47,32,51]), so that the universal solution structure extends beyond the critical region. The universality thus introduces serious restrictions on the possible behavior of phyllotactic systems not only in the critical region (see also study I). What do the foregoing considerations mean for phyllotaxis in general? Strictly speaking, the renormalization approach to universality discussed here applies only to Thornley’s model of phyllotaxis. But self-similarity and scaling was found in study I also for the Douady–Couder model of phyllotaxis, with the same value for α but another value for β (compare (1.4)). The reason for this difference is currently not understood (see the discussion in study I). Nevertheless, the fact that scaling was found for both models, together with the successful construction of a renormalization in this study, hints to a much wider universality beyond the two types of models considered so far. This universality may not cover the value for parameter rescaling, but the biological evidence for universality is anyway about the behavior of the divergences, and not about this other dimension of primordial growth. Despite this unsatisfying situation, the universality found here thus adds confidence in the relevance of the early determination mechanism for a large class of systems—this mechanism, that crucially depends on the assumption of scaling, was discussed in study I to explain the preferred selection of certain divergence angles found in phyllotactic patterns. Nevertheless, to further assure universality of scaling and divergence selection beyond the class of models considered here, a renormalization of other models of phyllotaxis should be attempted. As a prerequisite these other models should as well show the self-similarity of the odd part of the Fourier spectrum seen in Fig. 3. Whether the subsequent technical difficulties in analyzing the resulting renormalization transformation could be tackled needs to be seen. In conclusion, a theory like the one presented here may contribute to a better understanding of the principles underlying the emergence of the regularities of phyllotaxis across different plant species and other realms of nature where phyllotactic patterns are observed [3]. It seems that for an explanation the fundamental mechanisms like e.g. the biochemical details of meristematic growth need not be known. Instead, certain high level properties like circular symmetry and the mere existence of an inhibitory mechanism together with a convexity property are sufficient. Acknowledgments I dedicate this study to my teacher Prof. Dr. Hartwig Schmidt who showed me the beauty and power of renormalization theory. – Thanks are due to Dr. Otto Böhringer for careful inspection of an early stage of the manuscript.

C.H. Reick / Physica D 298–299 (2015) 68–86

Appendix A. Fourier representation of some inhibitor functions Here the Fourier representations of two types of inhibitor functions obeying the defining equations (2.3)–(2.7) are presented together with their asymptotic behavior. Polynomial inhibitor functions: This class of inhibitor functions is defined by



1

f (n) (x) = a x −

2n

2

− b for x ∈ [0, 1], n = 1, 2, . . .

(A.1)

with n

a=4

and b =

2n

1 2n

.

(A.2)

By applying (2.3) they can be mapped cyclically to all x ∈ R and the basic interval [−1/2, 1/2] in particular. In addition f (n) (x = 0) = 1. The cases n = 2 and n = 3 are called here biquadratic and bicubic. Fig. 4 in Appendix E shows a plot of the biquadratic inhibitor function. For the Fourier coefficients (3.2) one finds (n)

fk

=

n (2n + 1)! 

2n

(−1)

j+1

 −1 (2(n − j) + 1)!(kπ )2j ,

(n)

and f0 (n)

(A.3)

2n + 1

k→∞

π2

.

(A.4)

Sine inhibitor function: This inhibitor function is defined here by f (x) = 1 −

1 2

|sin(π x)| ,

x ∈ R.

(A.5)

Its Fourier coefficients (3.2) are 1 4k2

−1

,

k ̸= 0

(A.6)

and f0 = 0. Their asymptotic behavior is found to be af = lim k2 fk = k→∞

1 4

.

(B.3)

ϵ→0

Let n > 0. Then using the definition (3.2) for the Fourier components of f one finds after two partial integrations while invoking (iv) 1

lim f ′ (−ϵ) −

2π 2 n2 ϵ→0

1 2



1 2π 2 n2

dx f ′′ (x) cos(2π nx).

(A.7)

Appendix B. Positivity and asymptotic behavior of the Fourier coefficients of inhibitor functions

Because f ′′ is bounded and (piecewise) monotonic (‘‘Dirichlet condition’’), a number M exists such that (see e.g. [52, Section 104])

 1   2  M   ′′ . dx f (x) cos(2π nx) <   0  n

(B.5)

Thereby (B.4) gives as claimed relation (B.2). And with cos(2π nx) = 1 − 2 sin2 (π nx) one can rewrite (B.4) to obtain (B.1): fn =

1 2



1

π 2 n2

dx f ′′ (x) sin2 (π nx) > 0,

The renormalization introduced in Section 3 is based on the successive removal of Fourier coefficients of the inhibitory potential. This is essentially an operation on the index set of the Fourier coefficients. More generally, this operation can be formulated as a mapping over the natural numbers N+ , that is implied by a particular irrational divergence. This view is made more precise in this appendix. The notation introduced here, as well as the results obtained, are particularly needed for the derivation of the fixed point inhibitor functions in Appendix E and the analysis of their stability in Section 6. For better readability, the proof of formulas is shifted to additional subsections of this appendix. Throughout this appendix the divergence Θ is assumed to be irrational. Remember the notation (x)c := x −[x + 0.5] from (3.11), where [x] denotes the largest integer less or equal x. Confusion with the bracket notation for continued fractions is not expected.

Proposition. Consider functions f : R → R with the following properties:

HΘ = hn ∈ N+ , n = 1, 2, . . . hn < hn+1 ; |(hn Θ )c | >

In Section 3 the renormalization introduced in this study was termed ‘‘removal of heads’’. Heads were defined in Eq. (3.12) as the ordered set

 





  1  hn = n + Θ n − , 2

Then the Fourier coefficients fn (except f0 ) obey:

2

lim n fn =

n→∞

1

lim f (−ϵ) > 0. 2π 2 ϵ→0 ′

Θ 2

 . (C.1)

Note that the definition of HΘ depends on the choice of Θ . Below in Appendix C.2 it is shown that all elements of HΘ are obtained by

f is twice continuously differentiable over (0, 1) f ′′ is bounded over (0, 1). f is periodic: f (x) = f (x + 1) f is symmetric: f (x) = f (−x) f is convex over (0, 1), i.e. f ′′ (x) > 0. n = 1, 2, . . . ,



Appendix C. Heads and trains

C.1. Removal of heads as an integer mapping

fn > 0,

(B.6)

0

For the class of inhibitor functions considered in this study (compare Eqs. (2.3)–(2.5)) the positivity of their Fourier coefficients is proved in this appendix, together with their asymptotic 1/n2 decrease. In this paragraph f ′ and f ′′ denote the first and second derivative of the inhibitor function f (x).

(i) (ii) (iii) (iv) (v)

(B.4)

0

where the last inequality follows from the convexity of f .

= 0. Their asymptotic behavior is found to be

= lim k2 fk(n) =

fk =

lim f ′ (ϵ) = − lim f ′ (−ϵ) > 0.

ϵ→0

j =1

k ̸= 0

af

Proof. In the following, without further mention, continuity and boundedness of f and f ′ over (0,1) will be invoked, that follow from (i) and (ii). This implies in particular the existence of the limits limϵ→0 f ′ (±ϵ). Moreover, (iii)–(v) imply that f ′ (1/2) = 0 and f ′ is discontinuous in x = 0, or, more precisely

fn = 2n + 1

81

(B.1) (B.2)

 := Θ

Θ , n = 1, 2, . . . . 1−Θ

(C.2)

As proved below in Appendix C.3 the natural numbers ηk , k = 1, 2, . . . of the complementary set HΘ = N+ \ HΘ (compare (3.13)) are obtained from

ηk =



k

Θ

+

1 2



,

k = 1, 2, 3, . . . ,

(C.3)

82

C.H. Reick / Physica D 298–299 (2015) 68–86

Table 2 Initial elements√of the first eleven trains for the golden divergence ΘG [1, 1, 1, . . .] = ( 5 − 1)/2 = 0.618033988 . . . .

=

This is seen by noting that the sets HΘ and HΘ are complementary. C.2. Characterization of the set of heads HΘ by (C.2) and some related relationships In this section it is proved that the numbers hn given by (C.2) fully exhaust the set of heads HΘ defined in (C.1). To this end first a number of preparatory relationships are derived. Throughout this section the following abbreviations are used:

 := Θ

Θ , 1−Θ

(C.9)

which was already introduced in (C.2), and

  1 (x)r := x − 2

which are ordered as well. This latter formula was introduced by Shraiman without proof in [40]. The renormalization introduced in Section 3 works such that successive Fourier components of the inhibitor potential with a head index hn ∈ HΘ are removed. At the same time the divergence is mapped to a new value via (see (3.18))

G (Θ ) =

1

Θ

 −

1



Θ

.

(C.4)

Hence, in subsequent renormalization steps heads related to the new divergence G (Θ ) are removed, and so on. To handle such repeated renormalization steps it is convenient to modify the notation slightly, by writing hΘ (n) := hn

(C.5)

ηΘ (k) := ηk .

In addition, the successive application of functions is represented as f ◦ g (x) := f (g (x)). By employing this new notation, in a second renormalization step the heads hG(Θ ) (n), n = 1, 2, . . . are removed from the set of Fourier components, obtained after the first renormalization. Accordingly, with respect to the original Fourier coefficients (those before the first renormalization), the coefficients with index ηΘ ◦ hG(Θ ) (n), n = 1, 2, . . . , are removed. Considering a third renormalization, the components ηΘ ◦ ηG(Θ ) ◦ hG2 (Θ ) (n), n = 1, 2, . . . , are removed. And so on. It is convenient to arrange the sequence of indices removed in successive renormalization steps into trains: The n-th train is defined by (n)

tk (Θ ) := ηΘ ◦ ηGΘ ◦ ηG2 Θ . . . ◦ ηGk−1 Θ ◦ hGk Θ (n), with k = 0, 1, 2, . . . .

(C.6)

Note that the first element (‘‘head’’) of the n-th train is the head (n) number hn = hΘ (n) = t0 . As an example Table 2 contains the initial elements of the first eleven trains for the golden divergence. The first column of table entries is the sequence of heads removed in the first renormalization step; in Fig. 2 these are the Fourier components at the local minima. The second column gives the indices of the original Fourier components removed in the second renormalization step. And so on. The rows are the trains, starting with the heads in the first column. Two additional features of the trains are of interest here. First, it follows from (C.6) that successive train values can be obtained recursively from (n)

(n)

tk+1 (Θ ) = ηΘ ◦ tk (GΘ ).

(C.7)

And, second, it is interesting to note that all trains together fully exhaust the natural numbers, i.e. for every irrational Θ



(n)



N+ = tk ; k = 0, 1, 2, . . . ; n = 1, 2, 3, . . . .

(C.8)

1

(3.11)

= x − [x] − ,

(C.10)

2

c

where ‘‘r’’ stands for ‘‘remainder’’.

 (n − 0.5)] (see (C.2)) can be 1. The relation n → hn = [n + Θ inverted:  1+



hn



(C.2)

= 1+ n−

 1+Θ

Θ



2

= n,

(C.11)

where in the last step Θ < 1 was used. 2. One finds:

   1  n− (1 − Θ ) Θ = 2

(C.9) (C.10)

r

   1  n− − Θ 2 2 2    1 1 (C.2)  = Θ hn − − [ Θ hn ] + [ Θ hn ] − Θ n − 



 n− = Θ n+Θ

1



−

1

2

2

(C.2) (C.10)



   1  n− = (Θ hn )r + Θ n + Θ Θ 2    1  n− − Θ 2      1 1  n− + = (Θ hn )r + (1 − Θ ) Θ 2

= (Θ hn )r ,

r

2

(C.12)

where the last step was obtained by observing that the Gauss bracket is zero, because its argument is element of the open interval (0, 1). 3. For all x ∈ R

|(x)c | >

Θ 2

⇔ |(x)r | <

1−Θ 2

.

(C.13)

This is proved as follows. Because (x)r and (x)c are cyclic with period 1, it is sufficient to consider x ∈ [− 12 , + 12 ). First the case

0 ≤ x < 21 is considered. In this case (x)r = (x)c − 12 , (x)c > 0, and (x)r < 0. Using this, one can see that (x)c > Θ /2 is equivalent to (x)r > −(1 − Θ )/2 from which the claim follows. The case of negative x can be proved in a similar way. Now, the main claim of this section is proved, namely that indeed the set of head values HΘ is given by the values hn defined in (C.2). The proof proceeds in two steps: First hn ∈ HΘ is shown, and then that HΘ contains besides these hn no other numbers.

 (n − 1 ) ∈ [− 1 , + 1 ) Eq. (C.12) shows that Step 1: Noting that Θ 2 r 2 2 |(Θ hn )r | < (1 − Θ )/2, so that with (C.13) the elements hn obey the condition |(Θ hn )c | < 1/2 that defines the members of HΘ . Step 2: Here it is shown that the numbers hn defined in (C.2) indeed exhaust HΘ . The proof proceeds by showing that a number ζ that 



C.H. Reick / Physica D 298–299 (2015) 68–86

is not member of the series of hn defined in (C.2), does not belong to HΘ . Proof: Because ζ is not member of the series of hn , there exists a n such that hn < ζ < hn+1 . Then (C.11) and (C.9) imply

Therefore, one can distinguish two cases: (C.15)

(ζ Θ )r = −ζ (1 − Θ ) −

1 2

+ hn − [hn Θ ].

(C.17)

Introducing ∆ := ζ − hn , so that ∆ ≥ 1, integer, into this equation gives

(ζ Θ )r = (hn Θ )r − ∆(1 − Θ ).

(C.18)

Since by step 1 of this proof hn ∈ HΘ , Eq. (C.13) gives (hn Θ )r < (1 − Θ )/2, so that

(ζ Θ )r < =−

1−Θ 2

1−Θ 2

− ∆(1 − Θ ) ≤

1−Θ 2

− (1 − Θ )

.

(ζ Θ )r = (hn+1 Θ )r + ∆(1 − Θ ).

(C.20)

Because of (C.13) and ∆ ≥ 1

(ζ Θ )r > −

1−Θ 2

+ ∆(1 − Θ ) ≥

1−Θ 2

.

(C.21)

Hence, also in this case ζ is not member of HΘ . C.3. Characterization of the complementary set HΘ and some related relationships This section contains the proof that the values of ηk defined in Eqs. (3.14)/(C.3) indeed exhaust the full set HΘ complementary to the set of heads HΘ (compare (3.13)). In addition the formula (3.20), basic to the derivation of the renormalization of the age functions (3.21), is proved here. For the followingproofs it is convenient to show first that the mapping k → ηk = Θk + 12 in (3.14) can be inverted:

      1 k 1 1 Θ ηk + = Θ + + 2 Θ 2 2      k k k 1 1 = Θ − + + + Θ Θ Θ 2 2     k 1 = k−Θ + = k, Θ c 2

  − Θ Θ ′k c

= (Θ ηk )c .



1

Θ

−

1 Θ

(C.22)

the proof of (3.20) follows:

  1 = −Θ Θ ′ k − Θ ′ k + 2    k 1 1 ′ = −ΘΘ k + Θ + −k Θ 2 Θ

(C.23)

    (Θ ηk )c  = Θ  Θ ′ k  < Θ , c

(C.24) 2 where the inequality was obtained because (. . .)c ∈ (−0.5, 0.5). This shows that according to the definition (3.13) the numbers ηk are members of HΘ . Step 2: Consider a number ζ not belonging to the set of ηk , i.e. there is a number k such that ηk < ζ < ηk+1 . Because of (C.22) one thus has



k = Θ ηk +

1



2

    1 1 ≤ Θζ + ≤ Θ ηk+1 + = k + 1. (C.25) 2

2

    Because Θ ζ + 12 is integer this implies either (i) Θ ζ + 12 = k   or (ii) Θ ζ + 21 = k + 1. In case (i) we set ∆ := ζ − ηk and find (ζ Θ )c = ζ Θ − [ζ Θ + 0.5] = 1Θ + ηk Θ − k =   1 Θ (C.22) = 1Θ + (ηk Θ )c > ∆ − Θ≥ 2

2

(C.26)

because according to step 1, ηk ∈ HΘ , i.e. (ηk Θ )c ∈ (−Θ /2, +Θ /2), and ∆ ≥ 1. Therefore ζ is not member of HΘ . – In case (ii) we set ∆ := ηk+1 − ζ and find analogously

(ζ Θ )c = ζ Θ − [ζ Θ + 0.5] = −1Θ + ηk+1 Θ − (k + 1) = (C.22)

= −1Θ + (ηk+1 Θ )c   1 Θ < −∆ Θ ≤− 2

2

(C.27)

so that also in this case ζ is not element of HΘ . Appendix D. Train elements for fixed point divergences Θ (κ) As discussed in Appendix C.1, each irrational divergence Θ (n) leads to its own train series (tk )k=0,1,2,..., n=1,2,... . This is in particular true for divergences of type

where  k  the last equality was obtained  k  by1 noting that because ∈ (− 0 . 5 , 0 . 5 ) one has − Θ + 2 ∈ (0, 1). Noting that Θ c Θ c with (3.18) Θ ′ = G(Θ ) =

1



Next it is shown that the complement of the heads HΘ defined in (3.13) is indeed given by (3.14), i.e. by the set of ηk . The proof proceeds in two steps: First it is shown that every ηk belongs to HΘ . Then, in the second step, the proof is completed by showing that every other number ζ ∈ N+ outside the set of ηk is not member of HΘ . Step 1: From (C.23) it follows

(C.19)

Taking (C.13) into account, ζ is thus not member of HΘ . Case (ii): The procedure is analogous to case (i). Invoking assumption (ii) and using here ∆ = hn+1 − ζ one finds



2

(C.16)

These two cases are handled separately. Case (i): Noting that for non-integer x the relation −[x] = [1 − x] holds, one finds from the definition of (x)r and by invoking assumption (i) after some algebra:

1

k



+ − Θk = −ΘΘ ′ k + Θ Θ 2 Θ    1 = −kΘ Θ ′ + + Θ ηk Θ = −k + Θ ηk   1 (C.22) = − Θ ηk + + Θ ηk

n − 1 = [hn (1 − Θ )] ≤ [ζ (1 − Θ )] ≤ [hn+1 (1 − Θ )] = n. (C.14)

(i) n − 1 = [ζ (1 − Θ )] = [hn (1 − Θ )] (ii) n = [ζ (1 − Θ )] = [hn+1 (1 − Θ )] .

83



Θ (κ) = [κ, κ, κ, . . .],

κ ∈ N+ ,

(D.1)

which are fixed points of the Gauss map (compare (4.1)). As shown in this appendix, the associated train elements have certain simple properties. These are used to determine the fixed points of the inhibitor function renormalization in Appendix E. Proposition. For the divergences (D.1) the associated train elements are related by the one-step recursion (n)



(n)

tk−1 = Θ (κ) tk

+

1 2



,

n, k ∈ N + .

(D.2)

84

C.H. Reick / Physica D 298–299 (2015) 68–86

Proof. Note that Θ (κ) = G(Θ (κ) ). Then

Proof. By induction.

  1 (κ) (n) Θ tk + = 2   1 (C.7) n) = Θ (κ) η(tk(− ) + 1 2   (n) 

Proposition. For the divergences (D.1) the train elements asymptotically behave as

t

Θ (κ)

(C.3)

(n)

1

1

c

(n)

tk

2

Proposition. For the divergences (D.1) the associated train elements are related by the two-step recursion (n)

= κ tk +

(n) tk−1 ,

n, k ∈ N . +

(D.3)

Proof. This is a consequence of the previous proposition: (n) (C.7) tk+1 =

=

tk

(C.3)

η(tk ) = 

(n)



(n)

Θ (κ)

(κ + Θ (κ) )tk(n) +

+ 1 2

1



2



n) = κ tk(n) + tk(− 1. 

(D.2)

Proposition. For the divergences (D.1) the associated train elements can be written as (n)

= t1(n) Fk(κ) + t0(n) Fk(κ) −1 ,

n, k ∈ N + ,

(D.4)

where the extended Fibonacci numbers are recursively defined as in study I by (κ)

(κ)

+ Fk(κ) −1 ,

Fk+1 = κ Fk

k = 1, 2, 3, . . . , F0 = 0, F1 = 1.

Proof. With help of (D.3) this can be shown by induction.

(D.5) 

Proposition. For the divergences (D.1) one finds (1)

tk

lim

k→∞

(n)

=

tk

ηh1 + h1 Θ (κ) , ηhn + hn Θ (κ)

n = 1, 2, 3, . . . .

(D.6)

Proof. Noting that by (D.1) the fixed point divergences obey

−1  Θ (κ) = κ + Θ (κ) , one finds from the definition (D.5) of extended Fibonacci numbers: (κ)

Fk

lim

k→∞ F (κ) k+1

(κ)

 = κ + lim

k→∞

Fk−1

 −1

(κ)

Fk

lim

tk

k→∞

(n)

tk

(D.4)

=

(1)

+ t0 Θ

(n)

+ t0(n) Θ (κ)

t1

t1

(1)

(κ)

Appendix E. Fixed points of Rf In this appendix the fixed points of the inhibitor function renormalization Rf from (3.19) are constructed explicitly. For simplicity the index f on R and also on the normalization N is omitted here. Please note that this section draws heavily on the notation introduced in Appendix C and the results obtained there and in Appendix D. As obvious from (3.19), the defining equation for the fixed points is

(κ)

j

N

,

j = 1, 2, . . . ,

(n) Θ ◦tk (GΘ )

,

k = 0, 1, 2, . . . , n = 1, 2, . . . .

(E.2)

By applying the recursion (C.7) this is equivalent to

, 

k = 1, 2, . . . .

(D.7)

(E.3)

k+1

k

 (κ) −k  1 = √ (Θ ) − (−Θ (κ) )k , κ2 + 4

(E.1)

where  fj is the jth Fourier component of the inhibitor function at the fixed point (compare (3.2)) and ηj is the jth element of HΘ (compare (3.13)). Remember that ηj depends on the considered divergence. Since the fixed point problem (E.1) is of interest here only as part of the overall fixed point problem for the renormalization of the inhibitory potential (3.17), the divergences to be considered here must be fixed points of the associated renormalization equation for the divergence (3.18), namely Θ (κ) = [κ, κ, . . .] (see (4.1)). Since throughout this appendix divergences are always fixed point divergences, they are denoted simply by Θ in the following. Now advantage is taken of the train notation introduced in Appendix C. There it was shown that every irrational number defines a set of trains that all together pile up to the set of natural numbers (see (C.8)). Let us choose GΘ := G(Θ ) as the defining irrational number—we know that here GΘ = Θ , because Θ is a fixed point of G, but let us keep this for a moment. Hence, for given GΘ each index j in (E.1) can be identified with a particular (n) train element tk (GΘ ) so that the fixed point equation (E.1) can be rewritten as k

Proposition. For the divergences (D.1) extended Fibonacci numbers can explicitly be written as Fk

 fη

N ft (n) =  ft (n) ,

so that by applying (C.6) the claim follows.

(D.8)

   1 =  (Θ (κ) )−k t1(n) + Θ (κ) t0(n) (κ 2 + 4)   − (−Θ (κ) )k−1 t0(n) + Θ (κ) t1(n) .

N ft (n) (GΘ ) =  fη

= Θ (κ) .

Using this, one obtains (1)

n = 0, 1, 2, . . . .

Because Θ (κ) < 1, for k → ∞ the second term in the bracket goes to zero. Noting this, insertion into (D.8) proves the claim. 

 fj =

tk

= − log(Θ (κ) ),

Proof. Inserting (D.7) into (D.4) gives

where in the last step because of (. . .)c ∈ [−0.5, 0.5) and Θ (κ) < 1 the Gauss bracket turns out to be zero. 

(n) tk+1

k

k→∞



k−1 + + Θ (κ) 2 2     1 (3.11) (n) n) n) = tk(− = tk−1 + −Θ (κ) Θ (κ) tk(− + 1, 1

=

log(tk )

lim



where from hereon the reference to the fixed point divergence Θ is omitted. This shows that the renormalization is such that the Fourier coefficients of different trains, i.e. for different n, are unrelated. This implies that the Fourier coefficient at head index hn of the n-th train determines all the other Fourier coefficients (n) for that train, because with hn = t0 (compare (C.6) and (C.5)) it follows from (E.3)

 fhn  . ft (n) = k k

N

(E.4)

C.H. Reick / Physica D 298–299 (2015) 68–86

85

where in the last step the fact was used that the set of train elements is equal N+ (see (C.8)) so that the sum over index j could be replaced by a double-sum over all trains and their elements. Plugging in (E.10), evaluation of the geometric series and solving for  fh1 gives 1 − Θ2

 fh1 = 2

∞  η +h Θ 2  h1 1 n =1

Fig. 4. Fixed point inhibitor functions for three fixed point divergences computed from (E.10) and (E.12). For comparison also the biquadratic inhibitor function used in study I is shown (see also Appendix A).

It is now shown that for the fixed points considered here,

N can be derived from the asymptotic behavior of their Fourier coefficients. Since the limit (B.2) exists and is identical for all subseries one finds: lim j2 fj = lim

j→∞

(n)



tk

k→∞

2

 2 (E.4) (n)  ft (n) = lim tk N −k fhn . k→∞

k

(E.5)

(n)

Hence, N k ∼ (tk )2 for k → ∞, or log(N ) = 2 lim

(n)

log(tk ) k

k→∞

= −2 log(Θ ),

(E.6)

where in the second step (D.8) has been used. Accordingly

N = Θ −2 .

(E.7)

In a next step the relation between the Fourier coefficients at head indices is derived. This is obtained by noting that asymptotically the Fourier coefficients of all inhibitor functions approach zero in the same way, namely fn ∼ 1/n2 for n → ∞ (see Eq. (B.2)). To obtain the desired relation it is convenient to get rid of the normalization by dividing (E.4) by its special case n = 1. Moreover one can take the limit n → ∞ because the resulting relation holds for all n:

 ft (n)  fhn = lim k .  k→∞  fh f (1) 1

(E.8)

tk

Subseries have the same limit as the original series. Therefore the just mentioned asymptotic behavior leads to

  (1) 2  fhn tk . = lim (n)  k→∞ fh t 1

(E.9)

k

Here we are only interested in divergences Θ that are fixed points of G (compare (4.1)). Using their special properties, the right hand side limit in (E.9) was determined in (D.6) so that by entering (E.9) into (E.4) and using (E.7) gives

 ft (n) = Θ 2k k



ηh1 + h1 Θ ηhn + hn Θ

2

 fh1 .

(E.10)

By this equation all train elements are fully determined by  fh1 . Using the normalization (2.6) the value of  fh1 can be derived as follows: (2.6) (3.1) 1 =  f (x = 0) = 2

∞  j =1

 fj = 2

∞  ∞  n=1 k=0

 f (n) , tk

(E.11)

.

(E.12)

ηhn +hn Θ

Hence, Eqs. (E.10) and (E.12) fully determine the fixed points of the inhibitor function renormalization Rf . They show in particular that every fixed point divergence Θ = Θ (κ) is related to another fixed point inhibitor function. Some fixed point inhibitor functions are plotted in Fig. 4. References [1] R.S. Smith, S. Guyomarc’h, T. Mandel, D. Reinhardt, C. Kuhlemeier, P. Prusinkiewicz, A plausible model of phyllotaxis, PNAS 103 (2006) 1301–1306. [2] H. Jönsson, M.G. Heisler, B.E. Shapiro, E.M. Meyerowitz, E. Mjolsness, An auxindriven polarized transport model for phyllotaxis, PNAS 103 (2006) 1633–1638. [3] R.V. Jean, Phyllotaxis, Cambridge University Press, Cambridge, 1994. [4] W. Hofmeister, Allgemeine Morphologie der Gewächse, in: W. Hofmeister (Ed.), Handbuch der Physiologischen Botanik, Vol. 1, Engelmann, Leipzig, 1868, pp. 401–664. [5] J.C. Schoute, Beiträge zur Blattstellungslehre: I. Die Theorie, Recl. Trav. Bot. Néerlandais 10 (1913) 153–339. [6] J.H.M. Thornley, Phyllotaxis. I. A mechanistic model, Ann. Bot. 39 (1975) 491–507. [7] C. Marzec, J. Kappraff, Properties of maximal spacing on a circle related to phyllotaxis and to the golden mean, J. Theoret. Biol. 103 (1983) 201–226. [8] H.J. Scholz, Phyllotactic iterations, Ber. Bunsenges. Phys. Chem. 89 (1985) 699–703. [9] S. Douady, Y. Couder, Phyllotaxis as a physical self-organized growth process, Phys. Rev. Lett. 68 (1992) 2098–2101. [10] F. Rothen, A.-J. Koch, Phyllotaxis, or the properties of spiral lattices. I. Shape invariance under compression, J. Physique I 50 (1989) 633–657. [11] J. Guerreiro, Phyllotaxis: an interdisciplinary phenomenon, Physica D 80 (1995) 356–384. [12] H. Meinhardt, A.-J. Koch, G. Bernasconi, Models of pattern formation applied to plant development, in: R.V. Jean, D. Barabé (Eds.), Symmetry in Plants, World Scientific, Singapore, 1998, pp. 723–758. [13] G.P. Bernasconi, Reaction–diffusion model for phyllotaxis, Physica D 70 (1994) 90–99. [14] D. Reinhardt, E.R. Pesce, P. Stieger, T. Mandel, K. Baltensperger, M. Bennett, J. Traas, J. Friml, C. Kuhlemeier, Regulation of phyllotaxis by polar auxin transport, Nature 426 (2003) 255–260. [15] R.S. Smith, E.M. Bayer, Auxin transport-feedback models of patterning in plants, Plant Cell Environ. 32 (2009) 1258–1271. [16] P.B. Green, C.S. Steele, S.C. Rennich, Phyllotactic patterns: a biophysical mechanism for their origin, Ann. Bot. (1996) 515–527. [17] P.D. Shipman, A.C. Newell, Phyllotactic patterns on plants, Phys. Rev. Lett. 92 (2004) 168102. [18] P.D. Shipman, A.C. Newell, Polygonal planforms and phyllotaxis on plants, J. Theoret. Biol. 236 (2005) 154–197. [19] A.C. Newell, P.D. Shipman, Plants and Fibonacci, J. Stat. Phys. 121 (2005) 937–968. [20] A.C. Newell, P.D. Shipman, Z. Sun, Phyllotaxis: cooperation and competition between mechanical and biochemical processes, J. Theoret. Biol. 251 (2008) 421–439. [21] I. Adler, D. Barabé, R. Jean, A history of the study of phyllotaxis, Ann. Bot. 80 (1997) 231–244. [22] S. Douady, Y. Couder, Phyllotaxis as a dynamical self organizing process. Part I: the spiral modes resulting from time-periodic iterations, J. Theoret. Biol. 178 (1996) 255–274. [23] M. Kunz, Some analytical results about two physical models of phyllotaxis, Comm. Math. Phys. 169 (1995) 261–295. [24] M. Kunz, Dynamical models of phyllotaxis, Physica D 157 (2001) 147–165. [25] J. Guerreiro, F. Rothen, A global approach to botanic patterns, J. Theoret. Biol. 176 (1995) 233–245. [26] L.S. Levitov, Phyllotaxis of flux lattices in layered superconductors, Phys. Rev. Lett. 66 (1991) 224–227. [27] L.S. Levitov, Energetic approach to phyllotaxis, Europhys. Lett. 14 (1991) 533–539. [28] P. Atela, C. Golé, S. Hotton, A dynamical system for plant pattern formation: a rigorous analysis, J. Nonlinear Sci. 12 (2002) 641–676. [29] C.H. Reick, Self-similarity and scaling in two models of phyllotaxis and the selection of asymptotic divergence angles, J. Theoret. Biol. 313 (2012) 181–200.

86

C.H. Reick / Physica D 298–299 (2015) 68–86

[30] K.G. Wilson, J. Kogut, The renormalization group and the ϵ -expansion, Phys. Rep. 12 (1974) 75–200. [31] C. Domb, M.S. Green (Eds.), Phase Transitions and Critical Phenomena, Vol. 6, Academic, London, 1976. [32] M.E. Fisher, Renormalization group theory: its basis and formulation in statistical physics, Rev. Modern Phys. 70 (1998) 653–681. [33] M.J. Feigenbaum, Quantitative universality for a class of nonlinear transformations, J. Stat. Phys. 19 (1978) 25–52. [34] M.J. Feigenbaum, L.P. Kadanoff, S.J. Shenker, Quasiperiodicity in dissipative systems: a renormalization group analysis, Physica D 5 (1982) 370–386. [35] S. Ostlund, D. Rand, J. Sethna, E. Siggia, Universal properties of the transition from quasi-periodicity to chaos in dissipative systems, Physica D 8 (1983) 303–342. [36] B. Hu, Introduction to real-space renormalization-group methods in critical and chaotic phenomena, Phys. Rep. 91 (1982) 233–295. [37] A.J. Koch, J. Guerreiro, G.P. Bernasconi, J. Sadik, An analytical model of phyllotaxis, J. Phys. I France 4 (1994) 187–207. [38] I.S. Gradshteyn, M. Ryzhik, Tables of Integrals, Series, and Products, sixth ed., Academic Press, Orlando, Florida, 2000. [39] W.I. Smirnow, Lehrgang der höheren Mathematik, tenth ed., VEB Verlag der Wissenschaften, Berlin, 1971. [40] B. Shraiman, Transition from quasiperiodicity to chaos: a perturbative renormalization-group approach, Phys. Rev. A 29 (1984). [41] R.M. Corless, Continued fractions and chaos, Amer. Math. Monthly 99 (1992) 203–215.

[42] R.S. MacKay, Renormalization in area preserving maps (Ph.D. thesis), Princeton University, University Microfilm International, Ann Arbor, MI order 8301411, 1982. [43] T.L. Curtright, C.K. Zachos, Renormalization group functional equations, Phys. Rev. D 83 (2011) 065019. [44] N.K. Bary, A Treatise on Trigonometric Series, Vol. 1, Pergamon, Oxford, 1964. [45] E.T. Whittaker, G.N. Watson, A Course of Modern Analysis, fourth ed., Cambridge University Press, London, 1969. [46] J. Aczél, M. Kuzma, Generalizations of a ‘‘folk-theorem’’ on simple functional equations in a single variable, Results Math. 19 (1991) 5–21. [47] F.J. Wegener, The critical state, general aspects, in: C. Domb, M. Green (Eds.), Phase Transitions and Critical Phenomena, Vol. 6, Academic, London, 1976, pp. 8–126. [48] O.E. Lanford III, Circle mappings, in: H. Mitterer, et al. (Eds.), Recent Developments in Mathematical Physics, Springer, Berlin, Heidelberg, 1987, pp. 1–17. [49] G. van Iterson, Mathematische und Mikroskopisch-Anatomische Studien über Blattstellungen, Gustav Fischer Verlag, Jena, 1907. [50] J. Vieth, Biastrepsis and allomery of stems of dipsacus sylvestris mill, in: R. Jean, D. Barabé (Eds.), Symmetry in Plants, World Scientific, Singapore, 1994, pp. 231–245. [51] C. Reick, Universal corrections to parameter scaling in period-doubling systems: multiple scaling and crossover, Phys. Rev. A 45 (1992) 777–792. [52] H.S. Carslaw, Introduction to the Theory of Fourier’s Series and Integrals, third ed., Dover Publications, New York, 1930.