A computer-assisted proof of dynamo growth in the stretch-fold-shear map

The Stretch-Fold-Shear (SFS) operator is a functional linear operator acting on complex-valued functions of a real variable x on some domain containing in It arises from a stylized model in kinematic dynamo theory where magnetic field growth corresponds to an eigenvalue of modulus greater than 1. When the shear parameter α is zero, the spectrum of can be determined exactly, and the eigenfunctions corresponding to non-zero eigenvalues are related to the Bernoulli polynomials. The spectrum for has not been rigorously determined although the spectrum has been approximated numerically. In this paper, a computer-assisted proof is presented to provide rigorous bounds on the leading eigenvalue for , showing inter alia that has an eigenvalue of modulus greater than 1 for all α satisfying , thereby partially confirming an outstanding conjecture on the SFS operator.


Introduction
In kinematic dynamo theory, the induction equation for the growth of a magnetic field transported by a fluid flow is: [e.g. 3,10]. Here B is the magnetic field, u the fluid flow, ε the magnetic diffusivity and ε −1 the magnetic Reynolds number. The fluid flow u is assumed to be given, and could be steady or unsteady; it is usually taken to be incompressible, with ∇ · u = 0. For example, a seminal study of Arnold and Korkina [4] considers magnetic field evolution in the steady, three-dimensional ABC flow u = (cos y + sin z, cos z + sin x, cos x + sin y).
For any such given u as well as a general initial condition, the magnetic growth rate is γ (ε) = lim t→∞ t −1 log B(t) and the flow u is a dynamo if γ (ε) > 0 for some diffusivity ε > 0. For example, the above ABC flow is a dynamo as growth occurs for 8.9 ≤ ε −1 ≤ 17.5 and for ε −1 ≥ 27 up to the limits of present numerical simulations [15]. In astrophysical contexts, the magnetic Reynolds number ε −1 , a dimensionless measure of the electrical conductivity, is very large, and this motivates the study of kinematic dynamos in the delicate limit of zero diffusion ε → 0. If the dynamo growth rate γ (ε) remains positive and of order unity as ε tends to zero, the flow u is called a fast dynamo. In this case, we can define the fast dynamo growth rate γ 0 by γ 0 = lim inf ε→0 γ (ε). Delicate aspects of this limit become apparent if we look at the perfectly conducting situation ε = 0, when magnetic field evolution becomes intimately related to the streamline structure of the flow, as may be seen by using ∇ · u = ∇ · B = 0 to write (1) as Here the effect of the flow u on the magnetic field B is seen to be Lie-dragging: in other words, magnetic field vectors are advected, rotated and stretched with infinitesimal line elements in the flow. Another viewpoint is that magnetic field lines are materially transported in the flow, this being known as Alfvén's theorem. For a smooth flow u of any complexity, for example, one that is chaotic in the sense of having a positive topological entropy and so guaranteed exponential growth of line length in the flow, the magnetic energy increases exponentially, γ (0) > 0, in this case when the diffusion is strictly zero [21]. While it is easy enough to gain magnetic field growth in a chaotic flow with ε = 0, the growing field becomes indefinitely fine-scaled with frequent changes of sign, and the result is that the introduction of any magnetic diffusion ε > 0, no matter how small, can have a dramatic effect on the growth rates of field. For example, a chaotic flow in the periodic plane will have γ (0) > 0 but it is known through the Zeldovich anti-dynamo theorem [35] that γ (ε) < 0 for all ε and so γ 0 ≤ 0. The singular nature of the ε → 0 limit versus ε = 0 arises because the stretching is accompanied by the folding of field, and when these folds reach a diffusive scale, the effect of diffusion, no matter how weak, is dramatic and destructive.
Thus flows in the plane cannot be fast dynamos, or dynamos at all, and some threedimensional motions are needed. Perhaps the simplest class of fast dynamos, based largely on numerical evidence, consists of flows taking the convenient form u(x, y, t) = (u x , u y , u z ), independent of z. Field amplification operates via the so-called stretch-foldshear (SFS) mechanism [5], where magnetic field lines are stretched and folded through the action of u x and u y , and then shear in the z direction partially aligns the magnetic field vectors so that growth can occur. An example is the Galloway-Proctor 'circularly polarised' flow [16], Even in such flows, that in simulations show robust fast dynamo action, the limiting growth rate γ 0 is found numerically to be strictly less than γ (0); in other words, even very weak diffusion has a significant destructive effect on the highly convoluted magnetic fields that are generated, highlighting the delicate aspects of the ε → 0 limit. We remark that there is only one rigorous result of any generality, of Klapper and Young [21], that the topological entropy of a smooth flow field is an upper bound on the fast dynamo growth rate γ 0 , but there is no rigorous theory on lower bounds for smooth flows in (say) periodic three-dimensional space.
For this reason, a number of stylized models have been considered in fast dynamo theory. Here we consider one designed to illustrate the SFS mechanism, due originally to Childress and Bayly [5] and studied extensively in [17][18][19], which is where M SF is a folded baker's map: and M Sh is the shear in the z-direction: M Sh (x, y, z) = (x, y, z + αx), where α ≥ 0 is the shear parameter. Here a map is deployed as an idealized model of the motion of points and magnetic field vectors over a single period of a periodic flow such as (4).
is a complex-valued function of a real variable x on [−1, 1], the maps M SF and M Sh induce maps of b by the operators: This leads to the SFS operator T = T Sh T SF on b(x) : The phase shifts from the shear in Equation (8) generate an ever more complicated field T n b as n increases. To check whether this is a model for a fast dynamo, strictly speaking, some weak diffusive process should be introduced and growth rates studied in the limit ε → 0, of weak diffusion. However, an alternative route is to replace the averaging effect of diffusive smoothing by integration of T n b against a smooth test function c(x), namely to look at the growth of If (c, T n b) grows exponentially, with growth rate > 0 for general functions b and c, the map T is called a perfect dynamo. It is conjectured, with fair numerical evidence, that = γ 0 ; see [10] and references therein. In this way, via this so-called flux conjecture, one aims to link the growth for the fast dynamo limit of weak, non-zero diffusion, with the everevolving field for zero diffusion by making use of an appropriate measure of growth, namely integration against a smooth test function. With this in mind, the present paper will focus only on perfect dynamo action in the SFS mapping [17]. Once we look for exponential growth in (c, T n b), it makes sense to introduce the adjoint operator T * (with respect to the above inner product) such that (c, T n b) = (T * n c, b) for all b, c [5]. This adjoint operator T * is defined by Note that T * is the composition of two operators T * SF and T * Sh , where The key advantage in passing from consideration of the growth of T n b to growth in T * n c is that while iterating T on an initial condition creates ever finer oscillations in a complex field T n b, applying T * repeatedly giving T * n c stretches out and smooths structure, this being evident from numerical simulations [5].
Thus for the SFS model to be a perfect dynamo (and so by the flux conjecture a fast dynamo) we seek an eigenvalue of T and T * of modulus greater than unity; moreover, the smoothness of the eigenfunctions of the operator T or T * is the primary consideration. Since the operators T and T * are not compact in L 2 and their eigenfunctions are not generally smooth, it is convenient to restrict to a smaller function space which is a subset of L 2 .
Because of its analyticity-improving property, it is most convenient to study the operator T * (rather than T) restricted to a smaller space F for which the functions are not only smooth but also analytic in a disc on the complex plane. The restricted operator is denoted by S. This restricted SFS operator S is compact and growing smooth eigenfunctions of this operator can be obtained to establish perfect dynamo action in the SFS model. For further details on the SFS dynamo mechanism see [17] and for a discussion of passive scalar decay in baker's maps, where the effects of diffusion and boundary conditions are more delicate, see [19].
The above considerations lead to the following definition of the SFS operator which is the object of study in this paper.
where α ≥ 0 is the shear parameter. For α=0, S takes on the simpler form Following the above discussion, we seek magnetic flux growth in the form of a perfect dynamo, which corresponds to the existence of an eigenvalue λ of S, with |λ| > 1. In [17], it is conjectured that for all α > π 2 , S has at least one eigenvalue λ with |λ| > 1. In fact, numerical evidence suggests that for α > π 2 the principal eigenvalue satisfies |λ| > 1 and that, for all families of eigenvalues λ, lim α→∞ |λ| = √ 2. Figure 1 (reproduced from the paper [17]) shows the absolute values |λ| against α.
The purpose of this paper is to study the spectrum of S and, specifically, to prove the following theorem:  (1) c(x) ∈ F r , for some r > 1.
Note that in 2(a) it is not claimed that all eigenvalues have modulus less than one for α ∈ [0, π/2) (although numerical results suggest that this is the case), but rather that there is at least one eigenvalue with this property. Theorem 1.1 provides a partial proof of a conjecture contained in [17,18]. Its proof uses computer-assisted estimates to obtain mathematically rigorous results. In the context of kinematic dynamo theory, the theory demonstrates that the SFS model supports the growth of a perfect dynamo for the shear parameter α > π/2.

Properties of the SFS operator
The Stretch-Fold-Shear (SFS) operator is a family S α of linear operators, depending on the real non-negative parameter α, which acts on complex-valued functions c(x) defined on a domain containing the real closed interval [−1, 1]. In general, we shall write S for S α whenever the value of α is clear from the context or of no particular relevance.
For the analysis below, we restrict to S acting on a space F r of complex-valued functions analytic on the disc D r about 0 of radius r. Here r > 1 and may depend on the specific value of the parameter α (for further details see Section 3) but, in any case, [−1, 1] ⊂ D r . The operator S may be considered on other classes of functions but analytic functions are particularly suited for our purposes.
Because historically the focus is on the real closed interval [−1, 1], we denote the function argument by x rather than z. Moreover, depending on the context, we refer to the function as c or c(x), not restricting c(x) to denote the value of the function at the point x. The function spaces F r in which we work will be discussed in detail in Section 3.
We note that this operator is well defined for functions on D r , since for all x ∈ D r , both (x − 1)/2 and (1 − x)/2 ∈ D r provided r > 1. Below we review some properties of this operator, but here we note, in passing, that on the Banach space F r , S is a compact operator and therefore its spectrum consists of discrete eigenvalues together with 0. We refer to an eigenvector c(x) corresponding to an eigenvalue λ of S as an eigenfunction and the pair (λ, c(x)) as an eigenvalue-eigenfunction pair.
The operator S may readily be rewritten twice the odd part of f. The form (13) illustrates that for any even function d e (x) the function c(x) = e −iαx d e (x) is mapped to 0 by S, so that 0 is an eigenvalue of infinite multiplicity. Indeed, functions of this type form the 0-eigenspace of S.
Below, when required, we shall choose a normalization c(0) = 1, although often formulae are simplified by using other normalizations. From the formula (13) it follows immediately that any eigenfunction with non-zero eigenvalue is an odd function of x−1 and, hence, Except for α = 0 (considered in detail below) there does not appear to be a straightforward closed-form for the eigenfunctions except in one case due to Bayly. Taking c(x) = e iα(x−1) − e iα(1−x) = 2i sin α(x − 1), a straightforward calculation shows that (14) is satisfied for λ = e −iα . We refer to this example, or its normalized form sin α(1 − x)/ sin α, as the Bayly eigenfunction, denoted by c B (x). Although we have not proved it, the Bayly eigenfunction appears to be the only eigenfunction with eigenvalue λ = e −iα .

Banach space formalism
Let us now turn to the Banach spaces in which we shall work. Let r > 1 be a real constant and let D r = {x ∈ C : |x| < r} be the disc of radius r centred at 0. Let F r be the complex Banach space of complex functions f (x) analytic in D r such that the 1 -norm Where the value of r is not significant, we shall write F and D for brevity. It is possible to use other domains in C, in particular a disc centred at a point other than 0, e.g. centred at 1, but in this paper, a disc centred at 0 is sufficient for our purposes. Note that, for consistency, we use x to denote a complex variable as well as a real variable. Whichever variable is meant may be determined from context. Although other norms are possible in principle, the 1 -norm has some properties that make it useful for numerical computation. In particular, the following hold for functions f, g ∈ F r and scalars μ, One immediate application of these results is that the operator S : F r → F r is well defined and bounded. Indeed, | ± (x − 1)/2| = | ± r/2(x/r) ∓ 1/2| < ± r/2(x/r) ∓ 1/2 = (r + 1)/2 < r, since r > 1 and so, from (12), One advantage of working in the Banach space F r is that the operator S has good spectral properties. Indeed, it follows by standard results on composition operators that, for all α ≥ 0, the operator S = S α is a well-defined compact bounded linear operator S : F r → F r , so that its non-zero spectrum consists of discrete eigenvalues. (See [26] for a discussion of the principles involved in the proof of these facts.) In Sections 6 and 7 extended Banach spaces will be used, on which we place 1 -norms. In particular, we shall use the spaces C × F r (for eigenvalue-eigenfunction pairs) and C × F r × C × F r (for eigenvalue-eigenfunction pairs and their derivatives with respect to α, see Section 7) on which we define norms, both denoted by · as follows, for (a, c) ∈ C × F r and (a 1 , c 1 , a 2 , In both sections, Newton maps are defined to find solutions to eigenvalue problems.
We now turn to a discussion of the case α = 0, for which there is a complete theory.
The eigenfunctions corresponding to odd-degree polynomials for α = 0 can be readily computed using computer algebra and the first few are given in Table 1, normalized so that each is a monic polynomial.
The eigenfunctions in Table 1 may be expressed in terms of a generating function, analogously to classical Appel sequences, the most famous of which is the sequence of Bernoulli polynomials [8]. Indeed, following the approach in [33] and [34], we may write A standard calculation gives the identity so that, expanding G(x, t) as a power series in t we obtain the eigenfunctions given in Table 1 for odd n = 2k + 1 while c n (x) = 0 for n even.
The generating-function representation not only facilitates the calculation of the eigenfunctions for α = 0 but also provides a tool to obtain many properties of these eigenfunctions, similar to the large number of properties of the Bernoulli polynomials as tabulated in [1] and online resources such as [12]. This is no coincidence; indeed there is a strong connection between the eigenfunctions of S 0 and the Bernoulli polynomials. Recall that the sequence of Bernoulli polynomials B n (x) is given by the generating-function expansion Then, for odd n, the eigenfunctions c n (x) can be expressed in terms of B n (x): This means that many properties of the Bernoulli polynomials translate into analogous properties of the eigenfunctions c n (x) and vice versa.  (21) is zero, so that c n (x) is zero for n even in the expansion (21). This corresponds to the fact that every even polynomial is an eigenfunction with eigenvalue 0 for α = 0. Moreover, for odd n, the function c n (x) in (21) is a polynomial of degree n. Indeed, for m a positive integer, differentiating (19)  Finally, it is worth noting that for the type-(2k + 1) normalized eigenfunctions c 2k+1 (x)/c 2k+1 (0) → cos(π x/2) as k → ∞ and there are strong results on the asymptotics of the zeros of c 2k+1 (x) for α = 0. (See [11,25] for the relevant results for the Bernoulli polynomials.) The properties 1-3 in Theorem 4.1 are valid only for α = 0. But each eigenfunction corresponding to a non-zero eigenvalue is a polynomial of odd degree having precisely 2k + 1 roots in the complex plane, and it is possible to follow each eigenfunction as α increases from 0. We refer to the family of eigenfunctions emanating from the degree-(2k + 1) eigenfunction at α = 0 as the type-(2k + 1) eigenfunctions.
In what follows our aim is to choose k ≥ 1 and to show that for some α max sufficiently large, there is a type-(2k + 1) eigenvalue-eigenfunction pair for α satisfying 0 ≤ α ≤ α max , and by bounding |λ| to show that for π 2 < α ≤ α max , there exists an eigenvalue of modulus greater than 1, thereby confirming the conjecture in [17,18].

Recasting the eigenvalue problem as nonlinear root-finding
As is evident, the eigenfunctions in Table 1 are monic polynomials so that the type-(2k + 1) eigenfunction c 2k+1 (x) satisfies the normalization c (2k+1) 2k+1 (0) = (2k + 1)!. In general let us consider a continuous linear functional L : F → C and the normalization L(c) = 1. Note that the type-(2k + 1) eigenfunction in Table 1 . This particular normalization is computationally expensive to implement since it involves a (2k + 1)th derivative, so in practice, it is often better to use a simpler normalization such as L(c) = c(0). Note that, for a given normalization L, not all functions in F may necessarily be normalized since the linear functional L has a non-trivial kernel. However, for a given choice of k, a suitable normalization may be identified. In the work described here, the choice L(c) = c(0) has been used.
We now assume that a normalization functional L has been specified and we define a map G : Then a normalized eigenvalue-eigenfunction pair (λ, c) satisfies G(λ, c) = (0, 0), so that standard root-finding methods may be used, in particular the Newton method. Indeed, dG, the derivative of G at (λ, c) acting on tangent vectors (δλ, δc) ∈ C × F, is given by so that, for fixed α ≥ 0, the eigenvalue-eigenfunction pair can be found by the standard iteration scheme with quadratic convergence to a simple eigenvalue-eigenfunction pair. Indeed, close to a fixed eigenvalue-eigenfunction pair, the quadratic convergence ensures that the map G is a contraction mapping. For convenience, an approximation to the Newton map, chosen so that it shares the same fixed points as the exact map and still gives a contraction close to the eigenvalue-eigenfunction pair, is used to avoid inverting the derivative dG.
The Newton method may be adapted to rigorous computer-assisted functional analysis using standard interval analysis and its extension to function balls in the function space F, as described in outline below.

Analysis for α ∈ [0, 5]
To prove Theorem 1.1, we require two subsidiary, related theorems. The first, which is discussed in this section, gives the results on whether |λ| is greater or less than 1 for all α ∈ [0, 5], excepting an interval around π/2, namely, (1.5705, 1.571), for which less information is provided. (1) c(x) ∈ F r , for some r > 1. The proof of Theorem 6.1 is based on standard rigorous computer-assisted numerical functional analysis as outlined in [13]. We describe the method of proof in the following proof algorithm. The computer estimates are described in further detail below.

Proof algorithm
The method of proof is to divide the interval [0, 5] into a (rather large) set of subintervals such that the errors resulting by restricting α to each subinterval are sufficiently small for rigorous results to be obtained. In particular, a complex interval 1 is calculated that is guaranteed to contain an eigenvalue of the SFS operator for each α in the subinterval. These intervals provide the justification of the statements in Theorem 6.1. In practice, for each subinterval, high-precision calculations of an approximate eigenvalue-eigenfunction pair are performed non-rigorously for α set to the midpoint of the subinterval, and then rigorous calculations are performed on the whole subinterval. The parameters used in the computations are α-dependent; in general the calculation becomes increasingly timeconsuming as α increases and this puts a severe limitation on the maximum value of α for which it is possible to obtain rigorous results. A more detailed description follows.
Suppose α 0 = 0 and for i = 1, 2, . . . , we define a sequence of contiguous, possibly overlapping, real closed intervals For each i, we choose parameters, largely by trial and error, of the length δα i of the interval I i , of the radius r i > 1, of the precision of floating-point number calculations and of the degree of truncation of the Taylor series expanded on the domain D r i ⊇ [−1, 1]. We letᾱ i ∈ I i withᾱ i ≈ (α i + α i−1 )/2 and calculate to a high precision an approximate eigenvalue-eigenfunction pair (λ i ,c i ), either directly or using previously calculated results for i−1. In practice the non-rigorous Newton method is often an efficient and convenient method to obtain such a pair. Indeed we may use (a finite-degree version of) Newton's method (28) to obtain the approximation to a required degree of accuracy δ i , say. This process will also provide a finite-degree approximationd G i to dG (λ i ,c i ) .
In principle, many of the parameters may be allowed to vary with i, although in practice it is often convenient to keep the parameters largely fixed for a given run of the computer program covering a range of α in [0, 5]. The choice of parameters determines the Banach spaces F r i and C × F r i as described in Section 3.
Once an accurate approximate eigenvalue-eigenfunction pair has been calculated, the next step is to use computer-assisted methods to prove rigorous bounds for the eigenvalue-eigenfunction pair (λ, c) for all α ∈ I i . The idea is to fix α ∈ I i and to take a (small) closed disc-function-ball B = B ε (λ,c) around (λ i ,c i ) in C × F r i (equipped with a norm . ) and use an approximation to the Newton map to define a contraction on B using the following idea. Let N i be the approximate Newton map where F i is a fixed linear operator which approximates the inverse of dG (λ i ,c i ) evaluated atᾱ i . Note that N i does not have the quadratic convergence of the Newton method, but nevertheless should be a contraction, provided the parameters are chosen suitably. In fact, the derivative dN i(λ,c) = I − F i dG (λ,c) should be close to zero and certainly satisfy provided δ ≤ (1 − η)ε, a condition that must be checked. This means that N i maps B → B and is clearly a contraction (since η < 1) and so has a unique fixed point, for a fixed α ∈ I i . This fixed point gives the eigenvalue-eigenfunction pair. That N i is a contraction, and therefore each eigenvalue-eigenfunction pair is unique within the ball B i , will be an important part of the argument for the results in the next section. In this way, using several independent runs of the computer program, it is possible to obtain results for all the I i and thereby to cover the whole of [0,5]. This provides the proof of Theorem 6.1.

Analysis around α = π/2
The computer-assisted results described in the previous sections do not cover in sufficient detail the values of α close to π/2 and, in this section, we describe additional results which cover an interval of α straddling α = π/2. Specifically, we describe the proof of the following theorem: (3) (a) For α < π/2, |λ| < 1.
Combining this theorem with Theorem 6.1 gives Theorem 1.1. Again, we note that 3(a) makes no claim about all eigenvalues, only about a single eigenvalue found through the computer-assisted proof.
The proof of Theorem 7.1 requires a number of steps as described below. In particular, because of the existence of the Bayly family c B (x) of eigenfunctions described in Section 2, which crosses the type-3 family of eigenfunctions at α = π/2, the function c(x) is required to have zeros at x = −1, 1 to distinguish it from the Bayly family, so that the representation c(x) = (x + 1)(x − 1)c(x) has been used to ensure these conditions are satisfied. We describe the approach to proving Theorem 7.1 in a series of steps. The essential idea is to prove that, for each α in the interval [1.5705, 1.571], there exists a unique eigenvalue-eigenfunction pair in a closed ball in the Banach space C × F 1.1 (r is set to 1.1 for this work) and to obtain bounds on dλ/dα for α in this interval, to show that d|λ|/dα > 0, while at the same time showing that the Bayly eigenfunction for α = π/2 (suitably normalized) lies in this closed ball. This latter fact shows that the Bayly eigenfunction is in fact the unique eigenfunction in the ball for α = π/2 and since it is already established that λ = −i for the Bayly eigenfunction, the precise statement on λ in Theorem 7.1 is proved.
Let us write c(x) = (x − 1)(1 + x)ĉ(x), whereĉ(x) ∈ F r , for r > 1. Then the operator S induces an operatorŜ onĉ(x) as follows Note that the term x−1 in the denominator of Equation (30) cancels with x−1 in the term Sĉ(x) because the latter is an odd function of x−1. Modulo any function-space technicalities, the operators S andŜ have the same spectrum when S is restricted to functions c(x) satisfying c(±1) = 0. Now let us consider the operatorŜ given by Equation (30) on F r . Noting that for c(x) ∈ F r , e iαxĉ (x) ∈ F r with norm bounded by e αr ĉ , and [e αixĉ (x)] 2o /x ∈ F r with norm bounded by 2e αr ĉ , it follows thatŜĉ(x) ∈ F r with norm We note that ifĉ(x) ∈ F r then c(x) = (x − 1)(x + 1)ĉ(x) ∈ F r so that an eigenfunction c(x) ofŜ in F r gives an eigenfunction c(x) of S also in F r . For general α the Bayly eigenfunction c B (x) cannot be written in this way except for certain values of α (specifically α = nπ/2, n ∈ N). Indeed, for α = π/2, using the Euler product formula for the sine function gives sin whereĉ The steps in detail are enumerated below.
(1) An accurate approximate eigenvalue-eigenvalue pair (λ 0 , c 0 (x)) ∈ C × F 1.1 is calculated for α close to π/2 and a ball B around (λ 0 , c 0 (x)) is defined in C × F 1.1 . (2) For each α ∈ [1.5705, 1.571], it is shown, as in Section 6, that an approximate Newton map N is a contraction on B so that there is a unique eigenvalue-eigenvector pair (λ * , c * (x)) (dependent on α). (3) Using the expansion (34), it is shown that the Bayly eigenfunction pair corresponding to (−i, sin π(1 − x)/2) lies in B. This is a non-trivial computation for which the expansion (33) is used. Indeed, by truncating (34) at a suitably high, but finite, number of terms in the infinite product, and estimating the error in so doing, it is possible to show that the Bayly eigenfunction lies in a function ball contained in B. Because of the uniqueness of (λ * , c * (x)) ∈ B for each α, it follows that, for α = π/2, the pair (λ * , c * (x)) corresponds to (−i, sin π(1 − x)/2). (4) From the implicit function theorem, it follows that (λ * , c * ) depends smoothly on α.
Writingλ = dλ/dα andċ(x) = dc(x)/dα, the next step is to find a ball V containing the pair (λ,ċ(x)). The method used is an adaptation of the approximate-Newton method used to obtain the eigenvalue-eigenvector pair (λ * , c * (x)).

Computer-assisted estimates
The results in this paper depend on computer-assisted analytic estimates and, in this section, we outline the methodology we have used to obtain these estimates. The starting point for computer-assisted analytic estimates is the classical theory of interval arithmetic [2,20,[27][28][29] with rigorous error bounding to take into account computer rounding error. For finite closed intervals of real numbers, one defines unary and binary operations corresponding to the usual algebraic operations with output a finite closed interval of real numbers which contains the set of all results of performing the operations pointwise on the intervals, in such a way that the output interval is close to optimal. The Julia programming language [7] with the interval-arithmetic package IntervalArithmetic.jl [6] is a convenient programming environment that implements fixed (but arbitrary) precision interval arithmetic with an extension to intervals of complex numbers [24,29].
To obtain the functional analytic estimates, we use an extension of interval arithmetic to balls in function space, an approach that has been used extensively over the last 40 years [9,13,14,22,23,32]. Analytic functions are completely determined by their Taylor series about the centre of their domain. For computations in a computer with a finite amount of storage space, the Taylor series can be truncated after a certain number of terms. Detailed information on the Taylor polynomial of degree N (say) can be retained (by specifying complex intervals in which the Taylor coefficients lie); however, we must retain information on the norm of any higher-order terms. In fact, because certain operations on higher-order terms result in lower-order terms it makes sense to include an additional 'error' function to account for such terms, for which a bound on the may be kept. These considerations lead to the following definition of a complex function ball in F r [13].
Fixing r > 1 and the function space F r , for a given (fixed) positive integer N (the degree of truncation) and a sequence I P,i , i = 0, . . . , N of N + 1 intervals of complex numbers and two non-negative real numbers V H and V E , we may define a function ball B ⊆ F r to be the set of all functions f ∈ F r which can be written in the form where f P (x), the polynomial part of lower order, is a function f H (x), the higher-order part, is a function such that and f E (x), the error part, is a function such that In this definition, we have detailed information about f P since its Taylor coefficients lie in complex intervals that are known explicitly. But f H and f E are only known roughly by the bounds on their 1 -norms. The error function f E is introduced so that when lower-order terms are produced from higher-order terms as a result of, say, differentiating a higherorder function, we will not have detailed information about the coefficients, so that we may allocate the result to the error function, even if the terms are indeed of low order. The ideas of interval arithmetic may be extended to function balls, so that binary and unary operations on functions in F r may be extended to operations on function balls, and implemented on a computer, by adapting the operations to take into account computer rounding error. To achieve this, a suite of Julia programs was written, based on the IntervalArithmetic.jl package, to implement the algorithms described in Sections 6.1 and 7.
For the code that implements the computer proof see [30] and [31]. In practice, we have found that taking r = 1.1, polynomials of degree N between 10 and 25, and Float64precision has been sufficient to prove the results given in this article. The required choice of interval length for the parameter α varies with the size of α. For the proofs contained herein, interval lengths varying between 10 −5 (for small values of α) and 10 −4 (for α around 5.0) have been sufficient. Optimisation of the algorithm and code would likely increase the interval lengths.

Conclusion and further work
The work described in this article constitutes a preliminary study of the Stretch-Fold-Shear (SFS) operator for α ∈ [0, 5]. For α = 0, the spectrum is readily determined and the polynomial eigenfunctions may be obtained using generating-function techniques; the eigenfunctions are strongly connected with the classical Bernoulli polynomials. Despite the apparent simplicity of the SFS operator, the determination of its spectrum for α > 0 is a non-trivial endeavour and, although the spectrum has been approximated numerically, the computer-assisted results presented in this article are the first rigorous results to be obtained.
In this article, we have described a computer-assisted proof to obtain rigorous bounds for α ∈ [0, 5] on the eigenvalue of the type-3 eigenfunction, which corresponds, for α = 0, to a degree-3 polynomial eigenfunction. The methods can be readily adapted to other odd types (5, 7, etc.), and for larger values of α.
A word of caution is required here. The computer-assisted methods require accurate computations that become increasingly time consuming for higher eigenfunction types and even for relatively small values of α. Hence, while powerful in themselves, they do not appear to be applicable to the study of the spectrum for large values of α, let alone the determination of the asymptotic properties of the spectrum as α → ∞.
The generating-function approach is an elegant technique that provides powerful results for α = 0. It is unclear to what extent these techniques can be used for α > 0 and the investigation of their generalization to α > 0 is an important area of future research in this field.