Stability of a twisted Plateau border with line tension and bending stiffness

The stability of a Plateau border between three soap films is considered, taking into account the effects of line tension and bending stiffness in the border. A simple geometry is considered, in which the border initially lies in equilibrium along the axis of a circular cylinder, with three equally spaced films radiating outwards to meet the inside wall of the cylinder. The films are pinned at the two ends of the cylinder with a fixed relative twist, so the initial film surfaces are helicoids. The stability of this system to small perturbations, involving both the films and the border, is investigated as a function of the cylinder aspect ratio, twist angle, film surface tension, border line tension and border bending stiffness. Analytically, the stability problem is reduced to finding the first occurrence of a zero eigenvalue of an infinite matrix, which is then estimated numerically. The results from this calculation are in good agreement with full numerical simulations.


Introduction
Aqueous foams have been widely used in both domestic and industrial applications for many years, particularly due to their high interfacial area, surface activity and yield stress.Recently, renewed interest in foams, particularly their stability, has been driven by architecture (e.g.ARUP Beijing Water Cube, which also highlights their visual attractivenes; Carfrae, 2006) and their use as well-defined proxies This is because deformations of a border under negative tension with amplitude ε and wavenumber k will lead to an O(ε 2 k 2 ) energy gain from increasing the border length, whereas the energy cost from increasing the area of the adjacent surfaces is only O(ε 2 k) since the surface perturbation is confined to an O(k −1 ) distance away from the border.
So, in order to obtain a well-posed problem, as well as a negative line tension in the border, an additional effect is needed to stabilize the border at larger wavenumbers.The mathematically convenient and physically reasonable thing to do is to also introduce a (positive) bending stiffness to the border, with an associated energy cost of O(ε 2 k 4 ).At smaller wavenumbers k, the negative tension will dominate, leading to a destabilizing effect.At larger wavenumbers, the bending stiffness will dominate, leading to a net stabilizing effect.
For this initial study, we adopt the simplest well-posed border-mechanics model in a simple model system.The precise details of the border mechanics are likely to be important, not just for quantitative determination of the stability boundary but also when analysing foam configurations with curved Plateau borders, so this deserves further attention.(The latter can be caused by pressure differences between cells and/or the sagging of borders under their own weight.See, e.g.Embley & Grassia, 2007;Weaire et al., 2004.)In this paper we shall consider theoretically the effects of border tension and bending stiffness on the stability of a single Plateau border in a simple model system.Previously, Cox & Jones (2014) and Whittaker & Cox (2015) considered the stability of an ideal (zero-volume) Plateau border occupying the centreline of a circular cylinder of radius R, with equally spaced films radiating out to the cylinder wall, as shown in Fig. 1(b).The ends of the film were pinned to fixed equally spaced radii, while the sides were free to move on the internal curved wall of the cylinder.Different geometries were considered by varying the length L of the cylinder and the relative twist θ between the radii at the two ends.
The films and border were assumed to have zero thickness.The mechanics of the system was then governed by a uniform surface tension in the films, with the border being taken as entirely passive.In all geometries, there is an equilibrium configuration in which each of the films occupies part of a helicoid whose axis lies along the centreline of the cylinder.This configuration was found to be stable in shorter and less twisted geometries, and unstable in longer and more twisted geometries.The stability boundary was computed as a critical twist angle θ c as a function of the aspect ratio = L/R.For a given , the surface is stable for θ < θ c and unstable for θ > θ c .The instabilities manifested themselves by the Plateau border moving off the centreline and one or more of the films re-connecting on the curved surface of the cylinder.
Here we shall consider the same problem as in Whittaker & Cox (2015), but include the effects of both a line tension and bending stiffness in the Plateau border.Kern & Weaire (2003) showed that a liquid-filled Plateau border is well approximated by an idealized Plateau border with negative line tension, which depends on the liquid volume; we extend that idea here.Our aim is to determine the stability of a twisted Plateau border inside a circular cylinder, in terms of the twist, cylinder aspect ratio and the new border mechanics.We employ the same method as Whittaker & Cox (2015) of looking for a 'minimal-energy O(ε) perturbation that results in zero net energy change at O(ε 2 )' in order to find the stability boundary.

Set-up
We adopt the same set-up as in Whittaker & Cox (2015), with N 2 equally spaced radial films inside a circular cylinder of length L and radius R, as shown in Fig. 1.The relevant case for Plateau borders in foams is N = 3.However, the other cases N = 2 and N 4 may be relevant in artificial situations, where, e.g. one might wish to construct a frame and membrane structure that is energetically stable.The case N 4 turns out to be equivalent to N = 3, and so we use a general development below to cover N 3.The case N = 2 is slightly different, and is dealt with in Appendix D.

Geometry and mechanics
The radii of each film are pinned at the ends of the cylinder, but the contact lines of the films on the interior wall are free to move.The pinned radii are equally spaced at each end, but a fixed twist of angle θ is imposed between the two ends.The vanes are labelled by j ∈ {0, 1, . . .N − 1}.The angular offset of the jth vane relative to the 0th vane at the pinned ends is then θ j = 2π j/N.
Spatial positions are described using Cartesian coordinates (x * , y * , z * ).The origin is taken at the centre of one end of the cylinder, the z * direction corresponds to the axis of the cylinder and the x * direction corresponds to the direction of the 0th vane at z * = 0.
The vanes are all subject to a uniform surface tension (energy per unit area) γ .The Plateau border is subject to a line tension (energy per unit length) τ and a bending stiffness (coefficient of energy due to curvature) β.To ensure we have a well-posed and physically meaningful problem, we must have stability at large wavenumbers.This requires that we have γ > 0, together with either β > 0 or β = 0 and τ 0. (If β > 0 then τ can be negative.) Lengths are non-dimensionalized on the cylinder radius R. We therefore work with dimensionless Cartesian coordinates (x, y, z) = (x * , y * , z * )/R.We define the aspect ratio = L/R, which is thus the dimensionless length of the cylinder.We also introduce the twist parameter k = θ/ .Finally, we define two dimensionless parameters that give the relative importances of the two energies in the border compared with the energy in the surrounding vanes (The factor of N/2 in these definitions reflects the fact that a perturbation of the border will result in a average resistance equivalent to stretching N/2 vanes.)As in Whittaker & Cox (2015) we still have an initial equilibrium state in which the border lies along the axis of the cylinder, and the films are sections of helicoids.

Perturbation description
We must introduce notation to describe the perturbation to both the surfaces of the vanes and the line of the central border.So that we can work with O(1) perturbation functions, we introduce a small fixed parameter ε 1, chosen so that the typical dimensionless displacements are O(ε).Following Whittaker & Cox (2015), the unperturbed vane surface is parameterized by two dimensionless coordinates: ξ measuring radial distance from the cylinder axis and η measuring distance from one end of the cylinder.The perturbation of the jth vane is described by the dimensionless normal displacement εζ (j) (ξ , η) of the point originally parameterized by (ξ , η).The perturbed surface of the jth vane is therefore given by Each vane is pinned to the radii at the cylinder ends at η = 0, , meets the central border at ξ = ξ ( j) 0 (η) and meets the curved cylinder wall at ξ = ξ ( j) + (η).With no perturbations to the undeformed equilibrium configuration, we would have ζ ( j) (ξ , η) = 0, together with ξ ( j) 0 (η) = 0 and ξ ( j) The perturbation to the central border is described using two functions q(z) and h(z), which give the border's displacement at each dimensionless axial position z.The displacement is measured in the plane perpendicular to the cylinder's axis, with (dimensionless) components (εq(z), εh(z)) taken to be in the directions parallel and perpendicular to the equilibrium position of the 0th vane at that axial position z.The location r of the perturbed border is therefore given by Finally, it is also convenient to introduce notation for the perturbation to the angle at which each vane meets the Plateau border.We define εφ ( j) (z) as the change in the angle of the jth vane (from its equilibrium position) in the plane perpendicular to the cylinder axis at axial position z.
Observe that there is redundancy in this representation with ζ ( j) (ξ , η), q(z), h(z) and φ ( j) (z), as the vane displacements at the Plateau border must be compatible (kinematically and dynamically) with each other and with the displacements and rotation of the Plateau border.

Azimuthal expansions in θ j
In order to relate the border displacements (q, h) to the displacement ζ j of the jth vane, it is convenient to express the displacement of the border at each axial position in terms of its components (q ( j) , h ( j) ) Downloaded from https://academic.oup.com/imamat/advance-article-abstract/doi/10.1093/imamat/hxy063/5288764 by University of East Anglia user on 13 February 2019 parallel and perpendicular to the original orientation of the jth vane at that position.Using simple geometry, we obtain q ( j) (z) = q(z) cos θ j + h(z) sin θ j , (2.5) It is also convenient to decompose the rotation angle for each vane using Fourier-like series in terms of the angular offset θ j of each vane (2.7) For N = 3, the expression above is exact.For N 4, trigonometric functions with arguments of 2θ j and higher will also be needed.(For N = 2, see Appendix D.) To compute the coefficients ( φ, c, s) in terms of φ ( j) , we multiply the expansion (2.7) by either unity, cos θ j or sin θ j and then sum over j.The orthogonality of the trigonometric functions means that only one term on the right-hand side survives.Using the results from Appendix F of Whittaker & Cox (2015) for the sums of trigonometric functions and their products, we find that (for N 3) (2.9) (2.10)

Summary
The displacement of the jth vane is described by The transverse displacement of the border and the rotations of the vanes there are described by h(z), q(z) and φ(z), s(z), c(z) for 0 < z < 1 . (2.12) These border functions are related to the vane-specific functions h ( j) (z), q ( j) (z) and φ j (z) by (2.5)-(2.7).
The perturbation and domain functions in (2.11) and (2.12) are related and constrained by displacement boundary and matching conditions at the border, on the curved cylinder wall and at the cylinder ends.When considering the dynamics of perturbations, equilibrium conditions also apply at the border, at the other boundaries of the vanes and at each point of the vane surfaces.These conditions will be derived in the following two sections.

Displacement conditions at the border
The displacement (q, h) of the border and the rotation φ ( j) of the vanes at the border affect the domain of ξ for each vane and yield boundary conditions on ζ .These conditions were derived by Whittaker & Cox (2015) as As expected, the parallel displacement q ( j) of the border is associated with an expansion or contraction of the domain in ξ , while the perpendicular displacement h ( j) and the rotation φ ( j) set the boundary conditions on the film's normal displacement and displacement gradient at the border.See Fig. 2. (Note that the boundary conditions at the border at ξ = ξ ( j) 0 have been linearized back to ξ = 0, within the errors stated in (3.2) and (3.3).)

Displacement conditions at the curved cylinder wall
The domain boundary for ξ at the wall is set by the displacement of the vane there.The jth vane meets the circular wall x 2 + y 2 = 1 at ξ = ξ ( j) + (η).Thus, we have (3.4)This implies where ζ ( j) , η) is the perturbation at the wall, linearized back to ξ = 1.See Fig. 2.

Displacement conditions at the cylinder ends
At the two ends of the cylinder, the vanes are pinned to equally spaced radii, with an overall twist θ between the two ends.The vanes therefore cannot move from their original positions there, and so we have Similarly, the central border cannot move from its original position at the ends and the vanes cannot rotate there, and so we must have

Equilibrium conditions and neutral perturbation
For a perturbed configuration to have minimal energy it must be in equilibrium, or else there would be a nearby configuration with a lower total energy.Therefore, the area of the vanes must be extremal with regard to further internal perturbations, and also the overall energy must be extremal with respect to any further motion of the Plateau border and the contact lines with the cylinder walls.We shall now impose each of these conditions.

Vane area and cylinder wall contact line
Requiring that the surface of each vane is in equilibrium yields the same Euler-Lagrange equation as derived by Whittaker & Cox (2015), Equilibrium of the contact line with the wall yields the condition that each vane meets the wall normally.For a circular wall, Whittaker & Cox (2015) obtained the condition ∂ζ ( j)  ∂ξ

Transverse force equilibrium on the border
In Whittaker & Cox (2015), the only forces on the displaced border were the surface-tension forces from the N attached vanes.Since the vanes each exert the same force γ per unit length of border, the equilibrium of the border meant that the angles between the vanes remained the same, i.e. φ ( j) (z) was independent of j.But here the perturbed border also generates forces internally from its line tension and bending stiffness.Therefore, imposing equilibrium will mean different vanes must rotate by different amounts to achieve the right net surface tension force to balance the border's internal forces.
The equilibrium condition that balances the forces from the surface tension in the vanes with the forces from the line tension and curvature in the border at r(z) can be expressed as where the non-dimensional forces-per-unit-length-of-border are f γ from surface tension in the jth vane, f τ from the line tension in the border and f β from the bending stiffness of the border, all nondimensionalized using the scale 1 2 Nγ .From the angles of the vanes and the border, the surface-tension force from the jth vane is where εψ(z) is the angle between the border tangent dr/ds and the cylinder axis, and ϕ(z) is the angle of dr/ds projected on to the x-y plane relative to the original angle kz of the 0th vane.The components of (4.4) in the x and y directions can be deduced from Fig. 2(a), given that the border remains almost aligned with the cylinder axis.The O(ε) component of (4.4) in the axial direction is the leading term from a rotation of the border through an angle εψ(z) with an orientation ϕ − θ j − εφ ( j) relative to the jth vane.We expand the trigonometric functions in (4.4), both for small ε and using the multiple angle formulae.We obtain cos(θ j + kz + εφ ( j)  When summed over j, the terms in these expansions either evaluate to zero (sin θ j and cos θ j )1 or can be replaced by c(z) or s(z) using (2.9) or (2.10) (sin θ j φ ( j) and cos θ j φ ( j) ).We then obtain From the standard linear approximations to the line tension and bending forces in a beam (see, e.g.Howell et al., 2009, Section 4.4), we model the bending and tension forces in the border as ds2 , (4.10) (4.11)where r * = Rr is the dimensional position of the border, and s * = Rs is the dimensional arc length along it.The dimensionless arc-length s is related to the dimensionless axial coordinate z by ds = dr dz dz. (4.12) From (2.4) we find Differentiating the expression (2.4) for r with respect to z, and substituting it in to (4.10) and (4.11) using (4.14), we find Substituting the expressions (4.8), (4.15) and (4.16) into the equilibrium equation ( 4.3), and taking the two components in the directions b and n at O(ε), we have These equations express the condition of transverse force equilibrium on the border at axial position z.

Equilibrium of border attachment points at z = 0,
Since we are including a bending stiffness in the border, another boundary condition is needed at z = 0, in addition to (3.7).We assume that the attachment points of the border at the tube ends cannot support any torques.Thus, for the border to be in equilibrium at its ends, we require the border to have zero torque there.In the simple bending model, this corresponds to zero curvature of the border at the ends, i.e. d 2 r/ds 2 = 0. Using (2.4) and (4.14), this condition becomes correct to O(ε 2 ).

Equivalence to zero second-order energy variation
It is shown in Appendix B that the displacement and equilibrium conditions (3.1)-(3.7),(4.1), (4.2) and (4.17)-(4.19),when taken together, imply that the overall O(ε 2 ) energy change under the perturbation is zero.This means that these conditions are sufficient to determine the critical perturbation on the stability boundary.
(Since the original configuration is in equilibrium, the energy change under any perturbation is at most O(ε 2 ).If there is a minimal-energy perturbation that results in an O(ε 2 ) decrease in energy, then the initial configuration is unstable to that perturbation.If all minimal-energy perturbations result in O(ε 2 ) increases in energy, then the initial configuration is stable.The stability boundary therefore coincides with a change in sign of the O(ε 2 ) energy change of a critical mode.)
The additional functions (corresponding to border displacements and rotations) are written as follows: s n sin nπη . (5.5) The functions being decomposed here all vanish at η = 0, by virtue of (3.6) and (3.7), and so are extended to odd periodic functions on (0, 2 ).In the extended functions, the function and its first axial derivative will be continuous everywhere, but we expect discontinuities in the second axial derivative at η = 0 and η = .Therefore, the coefficients in the Fourier expansions above will, in general, decay only as n −3 as n → ∞.

Euler-Lagrange equation for surface equilibrium
With the above Fourier representation, the Euler-Lagrange equation (4.1) decouples for each Fourier mode, and we obtain for each n, where (5.7) From (3.2) and (3.3), the boundary conditions at the border are (5.8)For k > 0, we can use the changes of variables (5.9) to transform the system (5.6)-(5.8) to the hypergeometric equation (5.11) 2 .From Section 9.15 of Gradshteyn & Ryzhik (2000), the solution for f n (w) can be written in terms of two Gaussian hypergeometric functions F(a, b, c; w) and F(a − c + 1, b − c + 1, c; w).This solution is then transformed back to the original variables to obtain where (5.14) For k = 0, care is needed because λ n becomes infinite if is finite.However, by writing λ n = nπ/(k ), (5.6) can be solved directly when k = 0. We obtain a solution in the same form as (5.12), but with (5.15) These solutions can also obtained by carefully taking the limits of (5.13) and (5.14) as k → 0 and λ n → ∞ with kλ n = nπ/ fixed.

Curved-wall boundary condition
Using (5.1)-(5.5)and (5.12), the Fourier decomposition of the boundary conditions (4.2) at ξ = 1 implies (5.16) for each n and j.We use the axial Fourier decompositions of the expressions (2.6) and (2.7) to obtain expressions for h ( j) n and φ ( j) n in terms of h n , q n , c n and s n .Substituting these into (5.16),we obtain −q n sin θ j + h n cos θ j R n + φn + c n cos θ j + s n sin θ j + . . .= 0, (5.17) where we have introduced (5.18) (It is shown in Appendix A that the denominator of R n is strictly negative, so the division by it is always acceptable.) Equating the Fourier coefficients for θ j in (5.17), we find2 φn = 0, Substituting from the solutions (5.13)-(5.14)for S n and A n for k > 0, we find where (5.21) (5.22) (5.23) The details of the calculation are provided in Appendix A. For k = 0, inserting (5.15) into (5.18),we obtain the simpler expression This expression can also be obtained as the limit of (5.20), as k → 0 and λ n → ∞ with kλ n = nπ/ fixed.

Plateau border equilibrium conditions
We now decompose the border equilibrium conditions (4.17) and (4.18) into a half-wave Fourier sine series in z.We note that the odd derivatives introduce cosines, which cause coupling between the Fourier modes.(Since we are using half-wave Fourier series, the sine and cosine functions are not orthogonal.) The relevant half-wave Fourier sine series for a general cosine function is where3 (5.26) Decomposing (4.17) and (4.18) into their Fourier components, we find that (5.27) (5.28) These decompositions are quite subtle.There are two different ways in which the third derivatives in (4.17) and (4.18) could be represented, depending on the order in which the differentiation and the decomposition of the cosine are carried out.It turns out that a combination of the two is required in order to ensure that the expressions obtained for c n and s n decay as n → ∞.Full details of how the decompositions are carried out can be found in Appendix C.
Recall that as n → ∞, the coefficients c n , s n , h n , q n are O(n −3 ) and λ n = O(n).Hence, in the expressions (5.27) and (5.28) above, it can be seen that some of the individual terms on the righthand sides do not decay as n → ∞.This may appear to be problematic as we should have c n and s n being O(n −3 ) on the left-hand sides.However, the boundary condition (4.19) places some additional constraints on the q n and h n , which force the non-decaying components to cancel out between the terms in (5.27) and (5.28).The use of the specific decompositions derived in Appendix C ensures that this is the case.

Summary
Equations (5.19), (5.27) and (5.28) form a homogeneous linear algebraic system for the coefficients {c n , s n , q n , h n }.A non-zero solution of this system represents a non-zero equilibrium perturbation with zero net energy change at O(ε 2 ).The stability boundary coincides with a set of such solutions.

Stability boundary for N 3
The stability boundary for any number of vanes N 3 coincides since N is absent from the mathematical description at this point.The case N = 2 is slightly different, as already mentioned, and is discussed in Appendix D.
Eliminating c n and s n between (5.19) and (5.27)-(5.28), the stability boundary for N 3 is given by points in θ -space where there are non-zero solutions for {h n , q n } to the linear system We observe that the odd q n and even h n decouple from the even q n and odd h n , but the two sets of terms satisfy essentially the same sets of equations.(This corresponds to two deformations of the border related by a rotation of π/2 about the cylinder axis.)Both coupled sets of equations can be written in the form Mx = 0, (6.3)where x = (h 1 , q 2 , h 3 , q 4 . ..)T or x = (q 1 , −h 2 , q 3 , −h 4 , . ..)T , and the infinite matrix M has components for n, m ∈ {1, 2, 3, . ..}.The function R n is given in terms of hypergeometric functions in (5.20), and λ n = nπ/θ = nπ/(k ) as in (5.7).So, given values for T, B, k and , the matrix M is known.
Observe that the matrix M is real and symmetric.Solutions to the homogeneous equations (6.1) and (6.2) occur precisely when the matrix M has a zero eigenvalue.It is then a matter of finding the values of the parameters T, B, k and that makes this happen.Downloaded from https://academic.oup.com/imamat/advance-article-abstract/doi/10.1093/imamat/hxy063/5288764 by University of East Anglia user on 13 February 2019

Limit of a weak border (small T and B)
When T = B = 0, the matrix M becomes diagonal with diagonal elements R n .Zero eigenvalues therefore occur when R n = 0 for each n.These are precisely the solutions found by Whittaker & Cox (2015), and the n = 1 root corresponds to the stability boundary.
When T, B 1, the eigenvalues are still approximated by the diagonal elements, as the off-diagonal elements are small.With errors of O(T 2 , B 2 ), the zero determinant close to the T = B = 0 solution at R 1 = 0 occurs when For a fixed k, we then have where θ = θ 0 (k) is the solution of R 1 = 0.However, θ 0 (k) must still be found numerically.

Limit of small twist (small k)
In the limit k → 0, the matrix M also becomes diagonal, though care needs to be taken with the behaviour of λ n and R n .The former behaves as k −1 (but all the terms are regularized by sufficient explicit multiples of k) and the latter tends to the expression in (5.24).Hence, the limiting diagonal entries are given by where μ n = nπ/ , while the off-diagonal elements all vanish.Hence, the critical perturbation modes are precisely the individual axial Fourier modes.The possible solutions to M nn = 0 can be examined by considering the function The equation M nn = 0 is then equivalent to T = f (μ n , B).Furthermore, we observe that the nth axial Fourier mode is stable when T > f (μ n , B) and unstable when T < f (μ n , B).
The function f (μ n , B) is illustrated in Fig. 3.As can be seen there, when B > 0, the nth axial Fourier mode is always stable for sufficiently small and unstable for sufficiently large , whatever the value of T. As varies while B, T and n are held fixed, there is generally a single solution to T = f (μ n , B) dividing these two regimes.However, when 0 < B < B * = 0.02803 there is a finite range of negative T in which there are three solutions rather than just one.(See Fig. 3 inset.)For such values of T and B, the mode is stable for small , then as increases it passes through a finite region of instability, then a finite region of stability, before becoming unstable again for large .
When T = B = 0, the nth mode is stable for /n 2.6187 and unstable for /n > 2.6187.The additional finite unstable region arises physically because the combined effect of the border mechanics is to disproportionally destabilize wavelengths of O( √ −B/T).When the destabilized wavelengths are sufficiently shorter than 2.6187 and the destabilization is large enough but not too large, this can cause the additional finite window of instability to appear below 2.6187.This leads to the possibility that for Downloaded from https://academic.oup.com/imamat/advance-article-abstract/doi/10.1093/imamat/hxy063/5288764 by University of East Anglia user on 13 February 2019 Fig. 3. Main graph: plot of f (μ n , B) (as defined in (6.8)) as a function of for various values of B. When θ = 0, solutions to det(M) = 0 occur precisely when T = f (μ n , B).For given B and T there are either 1, 2 or 3 solutions for /n.A cusp catastrophe occurs at (B, T) = (0.02803, −0.56141)where the curve of f (μ n , B) has a stationary point of inflection (marked by a filled circle on the graph).For given parameter values, the nth mode is stable if the point ( /n, T) lies on or above the appropriate curve on the graph.Inset: numerically computed bifurcation lines in the B-T plane, dividing the region of 1 solution (clear) and 3 solutions (shaded) for /n.For a given B, the boundaries are given by T = f when ∂ f /∂ = 0.The upper boundary is almost indistinguishable from its small-B asymptotic form T = −(3/2 2/3 )B 1/3 (shown by the dashed line).
given parameter values, a point will be in the finite stable region for the n = 1 mode, but in the finite unstable region for the n = 2 mode.Examining Fig. 3, an example of this occurring can be seen when B = 0.001, T = −0.3 and = 1.4.With these values and n = 1, we have /n = 1.4,which is stable.But for n = 2, we have /n = 0.7, which is unstable.

Full numerical solution
For the full problem (for general k, T and B), the elements on the leading diagonal of M increase at least as fast as O(1, Tn 2 , Bn 4 ) (owing to λ n ).The off-diagonal elements increase at most as fast as O(Tn, Bn 3 ) (when m = n ± 1).So on later rows of the matrix the diagonal elements dominate.We should therefore be able to obtain approximate solutions by looking for zeros of the determinant of a truncated version of the matrix.
It is then a matter of finding the curves in θ -space where the determinant vanishes, and identifying which corresponds to the critical case of the first mode becoming unstable.Results obtained using an interval bisection method to find zeros for a fixed k and varying are shown in Figs 4-6.The method was implemented in Matlab.To avoid numerical overflow it is convenient to scale the (n, m)th entry of M by Bn 2 m 2 π 4 / 4 est (for B = 0) or Tnmπ 2 / 2 est (for B = 0), where est is the centre of the initial interval for .This ensures that the dominating diagonal entries M nn tend to something close to unity as n becomes large.(A fixed est is used in preference to the current value of , in order to reduce the number of matrix entries that have to be re-computed at each iteration.)Rather than computing the determinant, we Downloaded from https://academic.oup.com/imamat/advance-article-abstract/doi/10.1093/imamat/hxy063/5288764 by University of East Anglia user on 13 February 2019 only need to know its sign for the interval bisection method.This can be more conveniently determined from a singular value decomposition (see, e.g.Trefethen, 1997), rather than explicitly computing the determinant.Also, to speed up the calculations, the asymptotic forms for hypergeometric functions given by Watson (1918) can be used to approximate R n when λ n is suitably large.

Simulations using the Surface Evolver
The geometry of soap films is, to a large extent, driven by their minimization of surface energy, or surface area (at constant surface tension γ ).In seeking to determine the shape of the twisted vanes in a cylinder in the case N = 3, we choose to compare our solutions with the output of a full numerical simulation performed in the Surface Evolver (Brakke, 1992), a piece of software designed expressly to find the equilibrium shape of energy-minimizing surfaces.The Surface Evolver also provides energy methods to accommodate line tensions and bending energies.
Our simulations are of a rather standard type.The initial geometry has three planar films, with surface tension set to γ = 1, constrained by the sides and ends of a fixed cylinder of radius R = 1.The films meet along a line, the Plateau border, which initially lies straight along the axis of the cylinder, and has a fixed line tension.To the border we also apply a curvature, or bending, energy.
To be precise, we construct the energy functional given in (B1)-(B6).The line tension is constructed as an additional energy by summing the lengths of the line segments making up the central Plateau border (Evolver's edge_length method) and weighting the result with the required line tension τ = 3 2 T. In this way negative line tensions, which reduce the total energy, can be accommodated.Similarly, the bending energy of the Plateau border is also an energy quantity, involving the integral of the squared curvature along the Plateau border.This is constructed using Evolver's sqcurve_string_marked method (effectively the sum of the angle deficits at each vertex of the discretized Plateau border, normalized by the adjacent edge lengths), with weight 1 2 β = 3 4 B. (The additional factor of one-half here arises from the factor multiplying the square-curvature integral in (B6).) At the start of each simulation, the films are subdivided (discretized) into small triangles and the border into short edges.In cases in which the theoretical prediction is that there is a unique critical length for given twist angle θ , an initial length is chosen, slightly smaller than the predicted critical one, and the films are twisted to a fixed θ .The length is then increased in steps of Δ = 0.01, with all vertices of the tessellation shifted up affinely, followed by 1400 area-minimizing iterations and any necessary adjustments to the tessellation as described below.The critical length beyond which instability occurs is detected by checking the eigenvalues of the Hessian of energy (see Brakke, 1996): when there is one or more negative eigenvalue, we stop the calculation and record the last stable value of .If the stability boundary is not a single-valued function of θ (e.g.Fig. 6 with T = −0.5 and B = 0.02) we invert this procedure and (for small values of θ ) increase the twist angle towards the instability at fixed length , later merging this with the data obtained at fixed, large, θ by stretching.
Setting the bounds on the area and length of the triangles and sub-edges, respectively, is the main subtlety of the simulation.As Kern & Weaire (2003) point out, in the negative line tension model we should expect further instabilities associated with the discretization of the Plateau border into short edges, and care is required in choosing the bounds.The triangles of the film tessellation are subdivided whenever one of their sides exceeds a certain length (which depends upon radial distance from the Plateau border, we use 0.015 within a distance of 0.1R, where greater accuracy is required, and 0.25 beyond that) and removed whenever their area shrinks below a critical value of 5 × 10 −5 .We are also careful to further refine at high twist angles to remain within these bounds, having found such high Downloaded from https://academic.oup.com/imamat/advance-article-abstract/doi/10.1093/imamat/hxy063/5288764 by University of East Anglia user on 13 February 2019  refinement necessary to achieve the desired accuracy, despite the high computational cost.(Each line of data shown in Fig. 4, e.g. can require up to 2 months to generate on a desktop PC, although this can be reduced if the theoretically predicted critical parameters are used to inform the initial condition.)

Results and discussion
The results of our theoretical solutions and the Surface Evolver simulations are shown in Figs 4-6, where we see excellent agreement between the two methods.The lines represent the stability boundaries for the different parameter values, with the system being stable to the left of each line and unstable to the right.As expected, positive line tension T and positive bending stiffness B are both stabilizing effects.A combination of positive bending stiffness and negative tension will stabilize short wavelength perturbations and destabilize longer ones.We can further understand the effect of the border on the stability of the system when T < 0 and B < 0 by considering the energy changes in the border under two specific perturbations.First we consider the fundamental (n = 1) planar axial mode q(z) = sin(π z/ ) cos(k ), h(z) = − sin(π z/ ) sin(k ).Using the expressions (B16) and (B19) we can see that this perturbation of the border is energetically neutral when = π √ −B/T.Any other perturbations at this length and any perturbations at a shorter length will result in an increase in energy in the border.Hence, for < π √ −B/T the system must be stabilized compared with the T = B = 0 case.Hence, in Fig. 6 the region to the left of = π √ −B/T and below the T = B = 0 curve is necessarily stable.
Secondly, we consider the n = 1 vane-aligned twisted axial mode q(z) = sin(π z/ ), h(z) = 0.This represents the most unstable mode for the vanes in the absence of the border mechanics.Again using the expressions (B16) and (B19) we see that this perturbation is energetically neutral for the border when For larger values of , the most unstable vane perturbation results in a decrease in the energy in the border.Hence, in Fig. 6, the region to the right of this line and above the T = B = 0 curve is necessarily unstable.
One striking aspect of the results is the strong stabilization given to highly twisted configurations at low aspect ratios (including really extreme twists).Even a very small positive tension or bending stiffness will stabilize a sufficiently short (small for any fixed θ = 0) border.This is partly because the vane-driven instability is relatively weak for short highly twisted borders.However, it is also because the unstable perturbations in the absence of any border mechanics are aligned with the vanes, and would therefore result in highly twisted borders with a very short axial pitch.These large wavenumber twists are strongly stabilized by the border mechanics.
The same applies for a combination of negative tension and positive bending stiffness (see Fig. 6).Even for the fundamental axial mode, the vane-driven instability at large k would lead to a highly twisted border, which has a short wavelength and is therefore stabilized by the border mechanics.
We also observe that for positive T and B in Figs 4 and 5 the neutral stability curves fold back on themselves as θ increases.In some regimes, increasing θ at fixed first reduces the stability and then increases it again.This can be explained by the competing mechanisms from the vanes and the border.Increasing θ destabilizes the vanes, but as θ increases further, the most unstable vane mode is aligned with the twisted vanes in the base state, and there results in a highly contorted border.The resulting Downloaded from https://academic.oup.com/imamat/advance-article-abstract/doi/10.1093/imamat/hxy063/5288764 by University of East Anglia user on 13 February 2019 Fig. 6.Stability boundary results for B > 0 and T < 0 with B/T = −0.04,so the border is stable to short wavelengths, but unstable to long wavelengths.Continuous lines: theoretical results obtained by looking for the first zero of det(M) using a truncation at n = 400 rows.Points: numerical results from the Surface Evolver.The vertical dot-dashed line is = π √ −B/T, where the fundamental planar axial mode is energetically neutral for the border in isolation.The dashed line is (7.1)where the fundamental vane-aligned twisted axial mode is energetically neutral for the border in isolation.
border deformation is therefore stabilized by the border mechanics at higher θ .This effect is not seen in Fig. 6 where T < 0 and B > 0, as the result of the border mechanics dominating at large θ is that the stability boundary there lies close to = π √ −B/T where the border perturbation is energetically neutral.
Another interesting phenomenon arises because of the preferential stabilizing and destabilizing of short wavelengths by the border mechanics.A negative tension T will preferentially destabilize short wavelengths, while a positive bending stiffness B will preferentially stabilize even shorter wavelengths.For certain parameter values, this can lead to an unstable 'island' of wavelengths between a short-wavelength stable region (stabilized by border bending) and a longer-wavelength stable region (stabilized by the vanes).(The usual long-wavelength unstable region due to the vanes occurs at even longer wavelengths.)With the right parameter values, this can result in an n = 2 (or a higher) axial mode being unstable, while lower-order axial modes are stable.This is contrary to the usual expectation that the modes become more stable as their axial wavenumber increases.For example, examining Fig. 3, we see that when θ = 0, = 1.4 B = 0.001 and T = −0.3, the n = 1 mode is stable, while the n = 2 mode is unstable.Similarly, it is also possible, as illustrated by the (T, B) = (−0.52,0.0208) curve in Fig. 6, for an unstable system to be stabilized by increasing .

Conclusions
In this paper, we have considered the effect of line tension and bending stiffness of a Plateau border on the stability of a simple system involving a border and three vanes with uniform surface tension.
The stability boundary is found in terms of an aspect ratio , the twist angle of the vanes θ , a dimensionless border tension T and a dimensionless border bending stiffness B. The method involves decomposing a general perturbation as an axial Fourier series, and then applying conditions for it to have minimal energy and be a critical perturbation in the sense that there is no second-order energy variation.Whittaker & Cox (2015) considered the same problem in the absence of border tension and bending stiffness.They showed that twisting the vanes or lengthening the border destabilizes the system.The most unstable mode is the fundamental mode that occupies the whole length of the border and displaces the border sideways in alignment with one of the vanes.
As would be expected, positive line tension and positive bending stiffness are both stabilizing effects, particularly at short wavelengths.A combination of positive bending stiffness and negative tension (as is anticipated in real-life foams) will stabilize short wavelength perturbations and destabilize longer ones.We have quantified these effects theoretically, as shown in Figs 4-6.The results of our theoretical calculations are in excellent agreement with simulations carried out using Surface Evolver.
This work represents only a first insight into the effects of the mechanics of a Plateau border on the stability of a foam.It introduces the possibility of a new mechanism to promote foam break-up, as in anti-foam applications, by constructing a flow in which the Plateau borders are lengthened and twisted.The precise effect of the border mechanics on this mechanism within a real foam is not yet clear, as the stabilizing and de-stabilizing effects are dependent on the values of T and B.
Further theoretical work is needed in two directions.First, modelling of the border itself should be used to derive the actual mechanical response to deformations (and hence appropriate values for B and T-assuming the linear model is appropriate).It is possible that in addition to bending and stretching energies, there may also be energy penalties associated with axial variation in the twist of the vanes about the border axis, and that the dihedral angles between pairs of vanes should also be taken into account.Secondly, the results for the stability of a single twisted border need to be up-scaled and applied to the network of borders in a foam.
It would also be interesting to perform experiments, both on bulk foam systems and on the single border set-up considered here.For the latter (which would extend the experiments of Cox & Jones, 2014) it would be necessary to ensure gravity did not play a role, and to find a suitable way of varying T and B independently.One variation can be achieved by varying the cylinder radius R, but for the second variation, the liquid content of the border or the surface tension of the fluid would need to be altered.
where F 1 and F 2 are as defined in (5.21) and (5.22).Using (5.14) we have that the denominator in R n is where F 3 is as defined in (5.23).Combining the results (A.2) and (A.3) gives us the expression (5.20) for R n .
We can further examine the behaviour of the denominator, by making use of the result from formula 9.131 of Gradshteyn & Ryzhik (2000).Applying this result to the hypergeometric function for F 3 in (5.23), we obtain ( A . 5 ) Now, we know λ n > 0, so for k > 0 all the arguments of the hypergeometric function are strictly positive.We also have k 2 /(1 + k 2 ) ∈ (0, 1), and so the standard series representation of F converges.Since each term in this series representation is positive, we then have that the hypergeometric function above is strictly positive.Hence, F 3 > 0, and thus, from (A.3), we have Hence, the denominator of R n is always strictly negative when k ≥ 0.

Appendix B. Proof that the equilibrium conditions imply zero second variation in energy
In this appendix, we shall verify the statement made in Section 4. The total (dimensional) energy E * in the vane and border system can be written as where A * is the total area of the vanes, L * is the total length of the border and K * is half the integral of the curvature squared along the border.We non-dimensionalize these quantities using the surface tension γ and the length scale R Recalling the definitions (2.1) of T and B, we obtain the non-dimensional energy equation where Here V is the two-dimensional region of space occupied by the vanes, dA is an area element, B is the line occupied by the border and s is the arc length along the border.The area integral (B.4) was derived by Whittaker & Cox (2015).The other two integrals (B.5) and (B.6) arise straightforwardly by using (4.12) to convert from the arc-length coordinate s to the axial coordinate z.

B.2 Energy in the vanes
The dimensionless energy in the vanes due to surface tension effects is given by the integral for A in (B.4).An expansion for A was derived by Whittaker & Cox (2015), using the surface geometry and the fact that the perturbed surface satisfies the Euler-Lagrange equations (4.1). 4 We have where Inserting the expression (2.5) for q ( j) into (B.9)reveals that since N−1 j=0 sin θ j = N−1 j=0 cos θ j = 0. Using the equilibrium condition (4.2) and the boundary conditions (3.2) and (3.3) in (B.10), we find that We now substitute for h ( j) using (2.6), and make use of the expressions (2.9) and (2.10) in order to evaluate the sums involving φ ( j)

Energy in the border
From (B.3), the dimensionless energy in the border from the tension and bending stiffness is given by (N/2)(TL + BK ).
To evaluate the integral (B.5) for L , we use the expression After expanding the squares, we can use integration by parts to obtain Evaluating the right-hand side using (2.4) and substituting the result into (B.6),we obtain which are curvatures of the border in the directions perpendicular and parallel to the 0th vane.By (4.19), H and Q vanish at η = 0, and so, by the arguments above, the coefficients H n and Q n in the Fourier sine series will decay at least as O(n −3 ).(Here we have again used the shorthand notation μ n = nπ/ .)But since the series (5.3) for h and q converge sufficiently fast to allow first and second derivatives to be computed we have, by decomposing (C.1) and making use of (5.25), nm q m n 2 − m 2 , ( C . 3 ) ( C . 4 ) While the individual terms on the right-hand sides of (C.3) and (C.4) may not decay as quickly as O(n −3 ), overall H n and Q n will be O(n −3 ) as n → ∞.We can therefore differentiate twice term by term to obtain series for the second derivatives n sin(μ n η), ( C . 5 ) where By differentiating (C.2) once and then using (5.25), we can derive series for the first derivatives where We are now in a position to evaluate the Fourier decomposition of (4.17) and (4.18).Starting with (4.17), we first re-write the right-hand side using H and Q, to obtain c(z) = B H (z) + 2kQ (z) − k 2 H(z) − T H(z).
(C.13) Now, using the Fourier series in (5.5) and (C.2)-(C.12),and equating coefficients, we obtain This simplifies to (5.27).The same procedure is used to obtain (5.28) from (4.18).While some of the individual terms do not converge (e.g.μ 4 n h n = O(n)) as n → ∞, the boundary conditions at z = 0, Labelling these matrices M (0) and M (1) , their coefficients are given by The stability boundary for N = 2 can be found as above, by seeking points in parameter space at which one of M (0) or M (1) has a zero eigenvalue.

Fig. 1 .
Fig. 1.The set-up used in Whittaker & Cox (2015), with N vanes inside a right circular cylinder, bounded by between equally spaced radii at the ends, the curved interior surface and a central border.(a) N = 2 and (b) N = 3. (Figure adapted from Whittaker & Cox (2015) by permission of Oxford University Press.)

Fig. 2 .
Fig. 2. The geometry of the film near the Plateau border, adapted from Whittaker & Cox (2015) by permission of Oxford University Press.(a) The central region at axial position z with N = 3 vanes.The Plateau border is originally at x = y = 0 with the vanes shown by thick dashed lines.It is then displaced by ε(q (0) , h (0) ) and the vanes each rotated by an additional angle εφ.The vanes in this displaced and rotated configuration are shown by thick solid lines.(b) The coordinates for describing the perturbation to the jth vane of the film, between the displaced Plateau border at ξ = ξ ( j) 0 (η) and the cylinder wall at ξ = ξ ( j) + (η).

Fig. 4 .
Fig. 4. Stability boundary results for with B = 0 and varying values of T 0. Lines: theoretical results results looking for the first zero of det(M) using a truncation at 400 rows.Points: numerical results from the Surface Evolver.The region of stability lies to the left of each curve.

Fig. 5 .
Fig. 5. Stability boundary results for with T = 0 and varying values of B 0. Lines: theoretical results obtained by looking for the first zero of det(M) using a truncation at 400 rows.Points: numerical results from the Surface Evolver.As in Fig. 4, the region of stability lies to the left of each curve.Downloaded from https://academic.oup.com/imamat/advance-article-abstract/doi/10.1093/imamat/hxy063/5288764 by University of East Anglia user on 13 February 2019

)
The boundary condition (3.7) means that the boundary terms all vanish.To evaluate the integral (B.6) for K , we first use (4