Pressure-driven wrinkling of soft inner-lined tubes

Wrinkling of an inextensible elastic lining of an inner-lined tube under imposed pressure is considered. A simple equation modeling the elastic properties of the lining, the pressure, and the soft-substrate forces is derived. This equation aims to capture the wrinkling response of arterial endothelium to blood pressure changes. Numerical continuation is used to construct a bifurcation diagram as a function of the imposed pressure for in-plane deformations, in excellent agreement with weakly nonlinear theory, which we also develop. Our approach explains how the wavelength and amplitude of the wrinkles are selected as a function of the parameters in compressed wrinkling systems and shows how localized folds and mixed-mode states form in secondary bifurcations from wrinkled states.

Submitted to: New J. Phys.

Introduction
Lateral compression of a finite thin floating elastic sheet generates periodic wrinkles whose wavelength is the result of a balance between elastic forces and the restoring weight of the entrained liquid. On further compression, the sheet undergoes a transition from the wrinkled state to one characterised by a single fold [1]. However, wrinkling is not exclusive to floating elastica: the weight of the liquid can be replaced by other forces and used to generate wrinkling in both two-dimensional circular and three-dimensional spherical and curved geometries. Examples are provided by laterally compressed [2] or curved bilayer materials [3], as well as vertically loaded floating circular sheets [4,5,6] and spring-loaded interfaces [7]. In contrast, compressed or deflated spherical shells [8,9] exhibit buckling with no preferred length scale, as do elastic rings supporting a soap film [10,11,12,13]. Constrained buckling of elastic rings exhibits similar properties [14].
Understanding how surfaces wrinkle and then fold in different geometries under specific forces usually requires solving complicated systems of partial differential equations. The thin floating sheet in one dimension (1D) provides an exception. This system is not only modelled by a simple equation for in-plane deformations, but also turns out to be completely integrable in the limit of infinite extent [15,16,17,18,19]. As a result the remarkable shapes of both wrinkles and folds on thin floating sheets can be described using stunningly simple mathematical expressions [16,17], which naturally implies closed formulas for the wrinkling/folding thresholds in parameter space.
In this article we study the competition between in-plane wrinkling and buckling in a circular geometry within a similar framework. The results lead to greater understanding of a number of different systems where such competition is present. These include in-plane wrinkling of the elastic lining of an artery where wrinkled-tounwrinkled cycles driven by diastolic-to-systolic blood pressure changes may prevent clogging and adhesion of platelets via large changes in the local curvature of its endothelium [20,21]. Such cycling is likely to prove useful in other applications. A similar wrinkling instability is present in a rotating Hele-Shaw cell when a higher density fluid in the center is separated from a lower density fluid on the outside by an elastic membrane [22,23,24].
We construct an idealised two-dimensional model for this class of systems and compute strongly deformed states up to the point of self-contact, analyse their stability, and organise the results in the form of bifurcation diagrams. These diagrams describe the response of the system (compression, tension, maximum curvature) as a function of a control parameter, for example, the imposed pressure difference. We use the results to identify a transition from unwrinkled to periodic wrinkled states and then to folded states similar to what is observed in spring-loaded linings or tubular chitosan hydrogel surfaces [7,25]. Fold states arise via secondary bifurcations from the wrinkled state as in the one-dimensional case. Two cases are considered in detail, one with prominent wrinkling and a second one in which wrinkling is absent and only buckling remains.

The Model
To represent the lining on the inside of a soft tube, we consider an inextensible, infinitely thin membrane of length L = 2πR attached to a soft substrate as shown in figure 1. We suppose that in equilibrium (P = 0) the unlined soft tube has an inner radius r 0 < R (figure 1(b)) and hence that, when lined, the lining is forced to wrinkle. We model this force by an inward normal force per unit area F s = 1 2 K r(s) 2 − r 2 0 n (s) (figure 1(c)). Here r = r(s) denotes the lining profile (r is the distance from the tube centre) and s is the arclength. The substrate force F s is the simplest nonlinear model that is differentiable at r = 0 and that behaves like the classical Winkler foundation [26] when expanded around r 0 , with constant stiffness k = Kr 0 . Moreover, the quadratic contribution to the force vanishes in the flat-foundation limit, i.e. as r 0 → ∞, again recovering a Winkler-type foundation response. Although higher-order models [27], and in particular models that include nonlocal contributions [2], may provide a more realistic representation of the substrate forces, the Winkler model has been used extensively in studies of substrate-supported elastica and has provided important insights into the instabilities responsible for both wrinkled and localised states [7,28,29].
For in-plane deformations the resulting system is then described by where κ ≡ ∂ s φ is the local curvature. Here φ is the angle between the tangent plane and the horizontal or x-axis ( figure 1(a)). In terms of Cartesian coordinates (x, y) with origin at the tube center, ∂ s x = cos φ, ∂ s y = sin φ and r 2 ≡ x 2 + y 2 . The constants in (1) are the bending modulus B and the (unknown) tension T required to maintain the length L of the lining (T < 0 implies tangential compression). A brief derivation of (1) similar to that in [12] can be found in Appendix A. In the following we absorb the constant term 1 2 Kr 2 0 in the pressure P . The resulting system is then similar to a rotating Hele-Shaw cell filled with two fluids separated by an elastic membrane, with a higher density interior [22,23,24].
We define the natural length scale and introduce a dimensionless parameter that measures the perimeter of the lining in terms of λ, ≡ R/λ. We scale (1) according to The area within the lining, scaled relative to the area of the circle, is conveniently written via Stokes theorem as and, accordingly, its compression is ∆ ≡ 1 − S. The total energy, also scaled relative to the circle, is given by
This solution requires a simple relationship between the imposed pressure and the resulting tension, and serves as the starting point (order zero) for linear and weakly nonlinear analysis. Introducing a small parameter measuring the amplitude of a perturbation of the circle solution, we expand φ, x, y, T and P as follows: The coefficients of odd powers of in P and T vanish owing to the invariance of the system under rotations by half a wavelength. Substituting these expansions into (3) and the equations for x and y leads, at O( ), to To eliminate x 1 and y 1 , we compute (∂ 2 s L + L)[φ 1 , x 1 , y 1 ]: This equation reduces to an algebraic equation for the wavenumber m on assuming that φ 1 (s) ∝ sin(ms + δ): Modes with m = 0 (axisymmetric expansion) and m = 1 (translations) are excluded by inextensibility and pinning, respectively. Thus m ≥ 2 and solutions with integers m correspond to periodic states we refer to as wrinkles (W m ); δ corresponds to a rigid rotation of the solution, and can be set to zero. Thus Equation (8) is an important expression as it can be used to determine the critical pressure P * 0 for the onset of the wrinkling instability as the pressure increases and the wavenumber m = m * of the resulting wrinkles for a given . Figure 2 depicts P 0 as a function of m for three different values. The figure shows how the circular tube becomes wrinkled as P 0 overcomes the threshold P * 0 ≡ (− 5 +4 5/2 )/2 and the interior depressurises. It also shows how the choice of determines the order of appearance of new unstable wavenumbers. A simple formula gives the critical wavenumber at P * 0 : m * = 1 + √ 5 . When 5 < 9 the onset wavenumber is m * = 2 since m = 1 corresponds to translations (figure 2).
In terms of physical parameters, providing a key formula relating the critical pressure P * 0 for the onset of wrinkling to the geometry of the tube and the physical properties of the substrate and the lining. Expression (11) also indicates that the critical pressure can be tuned by a proper choice of r 0 and R, for instance, to generate lining wrinkles at pressure where the first subscript indicates the primary wavenumber and the second the new wavenumber introduced at the secondary bifurcation. The subscripts ± on F refer to the folded states with an extrusion (+) or intrusion (−). The subscripts s and a indicate whether these protrusions occur on the axis or off it. The letter B labels the buckling mode m = 2. A scale bar of unit length is included on the right. The inset shows the same results but over a larger range of P . equilibrium (P 0 = 0). As mentioned, this requires r 0 < R, i.e. that the lining has an excess of length over the unlined tube inner perimeter. Likewise, for large , the critical wavelength of the wrinkles in terms of physical parameters simplifies to We extend the above approach to compute periodic states with wavenumber m to higher order in (see Appendix B). We display the O( 2 ) expressions for φ 2 , x 2 , y 2 , P 2 and T 2 below: We computed the expansion to O 7 using computer algebra. From these results, we can compute the slope ∂P/∂T of the primary wrinkle branches at the bifurcation points given by (8). This slope is always positive unless m = 2 and 5 ≥ 81. The mode m = 2 is special, because of its maximum wavelength; this mode is the first one to emerge in the absence of the intrinsic scale [30], and we therefore refer to it as the buckling mode (B). In the following we extend the above results using numerical continuation and consider two cases. In the first (Section 4) substrate forces are substantial and wrinkling is present. In the second (Section 5) these forces are much weaker, wrinkling is absent and only buckling remains.

Numerical continuation: 5 = 576
To compute strongly nonlinear solutions, we implemented (3) as a boundary value problem in AUTO [31] (see Appendix C for details) and numerically continued different wrinkle states for a given starting from the circle branch satisfying (7). Each increment in P requires the solution of a nonlinear eigenvalue problem for the response T . The results show that the weakly nonlinear theory is remarkably accurate, even when = O(1) (see Appendix C for a comparison up to O( 7 ) when 5 = 576). The continuation approach also allows the computation of secondary branches of mixed modes (M) and folds (F). Figure 3 shows the compression ∆ as a function of the imposed pressure difference P for primary wrinkle states W m with different wavenumbers m, starting with W 5 corresponding to the onset wavenumber m * = 5. The figure shows not only the pressure required to initiate collapse of the tube (corresponding to ∆ = 0) but also its subsequent response to quasistatic increase in P , i.e., the figure represents the tube law describing the mechanical response for different modes of instability for the chosen value 5 = 576. Figure 4 shows another measure of the response of the system, the tension T , also as a function of P . The (P, T ) formulation provides the natural framework for numerical continuation. Both figures also show a number of secondary branches (the mixed states M and the fold states F) that bifurcate from the W states at finite amplitude, together with sample solution profiles at the locations indicated in the figures. All our plots use the same convention (colours and symbols).
While the circle solution (∆ = 0 in figure 3, black line in figure 4) exists for any pressure P , we observe primary branches W m of wrinkle states with different integer wavenumbers m only above the critical pressure P * 0 . Wrinkle solutions with wavenumbers below m * are interspersed with those above m * ; the wavenumber of the former decreases as P increases until m = 2; thereafter only wrinkle solutions with wavenumbers above m * are present and m increases monotonically with the pressure P . When m * is not an integer, the primary instability corresponds to the integer m nearest to m * provided m * ≥ 2. Figure 3 shows that the compression ∆ is almost proportional to the applied pressure P for all the wrinkle modes, i.e., that the modulus Y ≡ ∂P/∂∆ is approximately constant. Each W m branch ultimately results in selfcontact and at this point the continuation is terminated. Self-contact forces can be included as in [30,33], see also [14,34], but this has not been done here.
Besides wrinkle modes, numerical continuation reveals two types of secondary branches. Most commonly, secondary branches connect a primary mode with m ≥ m * to another primary mode with m < m * . Figure 3 shows that all intermediate solutions along the mixed-mode branch connecting m = 11 and m = 2 primary branches, i.e. connecting the points M 11,2 to M 2,11 , exhibit modulation at both wavenumbers. In fact, most of these interconnecting branches also result in self-contact, although longer, fully realisable interconnecting branches become possible as (and hence m * ) increases and the number of connections between W branches above and below m * grows. Secondary bifurcations that do not connect different primary modes are also present. These correspond to localised folds and come in pairs. The first pair F s ± bifurcates from W 5 with F s + representing a localised protrusion while F s − represents localised invagination. Both branches reach self-contact at almost the same point (figures 3 and 4). A family F a of asymmetric folds is also expected, but these states cannot be computed by AUTO with the imposed boundary conditions. Arrays of folds with different symmetries, analogous to those of [19], have also been found, with consistently higher degeneracy (see the yellow branches, e.g. F s + s − in figures 3 and 4). Figure 3 also reveals that the modulus Y drops dramatically along the F branches, a well-known consequence of the appearance of folds. In the case of the M branches, the modulus Y can be negative as is the case for the buckling mode B.
We also examined the energy E of the different wrinkled, folded and buckled states as a function of the compression ∆. For small compression the lowest energy solution corresponds to m * = 5, the natural wavenumber of the system for 5 = 576, as shown in figure 5. However, as the compression increases, the localised states F s ± bifurcate from the m * = 5 state, and the lowest energy state becomes F s + , with F s − at a slightly higher energy. This secondary bifurcation thus defines the wrinkle-to-fold transition, with threshold at ∆ c ≈ 0.084 for the particular case 5 = 576. The direction of branching of F s + and F s − is consistent with that leading to spatially localised states in the bistable Swift-Hohenberg equation [35]. For higher compressions, the F s ± are no longer realisable and other localised states correspond to global energy minima ( figure 5).
Finally, in figure 6 we plot the maximum curvature of the different states we have studied. A rapid increase in maximum curvature can be observed along all wrinkle branches after their bifurcation from the constant curvature circle solution. Larger m values result in faster increase in κ max . Folds and some mixed states display even faster increase in curvature after they emerge from secondary bifurcations. The transition between the wrinkle state W 5 and the fold states F s ± , the first one to take place, occurs at P = −217.7 (figure 4) and corresponds to κ max ≈ 3.09.
5. Numerical continuation: 5 = 0.005 When = 0 our problem becomes a pure buckling problem with no intrinsic length scale [30,33]. In this case it is known that the first buckling mode corresponds to m = 2 with more complex buckling modes requiring larger and larger pressures as the wavenumber m increases. Moreover, in this regime the governing equation involves the curvature κ only and the problem is analytically solvable [36,37,38].
To confirm that our model possesses the correct limiting behaviour and thereby validate our numerical continuation approach we take 5 to be very small and compare our results with those for = 0 and 5 = 576. Specifically, we take 5 = 0.005 and document the corresponding nonlinear results in Figure 7 for comparison with figures 3-6.
As expected, the first solution to emerge from the circle when 5 = 0.005 is m = 2, i.e., the buckling mode B, and the wavenumber of the subsequent solutions that emerge increases monotonically with the pressure difference P . Moreover, the appearance of these states requires positive values of P and the corresponding branches all behave in a similar fashion. These solutions are thus the expected buckled states. In figure  7 these states are still labelled W but this is only because is not identically zero. For these small values of there are no mixed modes or folds prior to self-contact. In fact, such secondary bifurcations move farther and farther out along each primary branch and beyond the point of self-contact as 5 → 0, and conversely, down each primary branch and towards the circle solution when 5 increases. This process leads, for sufficiently large , to the appearance of secondary bifurcations prior to self-contact and for negative values of P , as in figure 4. Figure 7a shows the compression ∆ as a function of P for 5 = 0.005 for comparison with figure 3 while figure 7b shows the tension T , also as a function of P . Figure 7a reveals that for smaller values of 5 the compression increases much more rapidly with P than for larger values 5 , a consequence of the absence of the stiffening effect of the substrate. These results are corroborated in 7b. The results in both figures are in accord with the weakly nonlinear theory: the modulus Y = ∂P/∂∆ is now positive for all wavenumbers m (figure 7a) and likewise all primary branches have positive slopes ∂P/∂T (figure 7b), even for m = 2, as predicted by the theory. Figure 7c shows the energy E as a function of the compression ∆ for 5 = 0.005 for comparison with figure 5. In contrast to the latter, E is now a monotonically increasing function of ∆ and the wavenumber m: for small 5 the bending energy dominates the substrate energy and its contribution grows with increasing compression. Thus the lowest energy state at a given compression is that with the lowest overall curvature, i.e., the wavenumber m = 2 state is the minimum energy state and so is stable until self-contact (cross). After this point, stability is transferred to the next lowest wavenumber solution, m = 3, etc.
In figure 7d, we plot the maximum curvature as a function of P for comparison with figure 6. The figure shows that for small 5 maximum curvature is reached much earlier as P increases than for larger 5 . However, in each case, the maximum value of κ max necessarily coincides with the point of self-contact and is identical to the corresponding curvature when = 576, i.e., κ max is independent of (figure 7, bottom right panel).
All this is in substantial contrast to the behaviour identified at larger 5 described in figures 3-6 but confirms that the solutions of (1) converge to the correct pure buckling limit as → 0.
Finally, the two lowest panels in figure 7 compare the profiles of the m = 2 and m = 5 solutions for the two different values of considered in this work. The comparison is made at a point 10% from the critical pressure P 0 for m = 5 and 2% from the critical pressure for m = 2 and again at the point of self-contact for both (red × and * symbols, respectively). We see that when is large the amount of compression for given ∆P is substantially less than for smaller . Thus the wrinkling or buckling process occurs over a smaller interval of P as decreases. However, at the point of self-contact the profiles in the two cases are identical and independent of the parameter as suggested by the weakly nonlinear analysis.
Evidently, as 5 decreases and the influence of the substrate wanes the bifurcation diagrams simplify dramatically and in the absence of the second length scale the system approaches the corresponding result for the unsupported ring ( = 0). This simplification arises because the secondary branches leading to both mixed modes and the folded states move past the point of self-contact thereby ceasing to be realisable. In this case the first primary mode is the lowest wavenumber mode, m = 2. Subsequent primary modes now come in monotonically with increasing m and all behave in a similar fashion. However, despite these changes the primary branches continue to bifurcate subcritically, in the sense that the lining loosens (tension T becomes less negative), as P increases.
On the other hand when 5 increases the wavenumber m * of the mode that first sets in also increases (figure 2). This fact leads to repeated mode jumping. For example, m * = 4 for 5 = 320 while m * = 5 for 5 = 576. Thus the mode m * = 4 remains dominant only over a finite interval of 5 , and as 5 increases m * = 4 is replaced by a new dominant mode, m * = 5. This transition is associated with a so-called codimension-two point where the dispersion relation (8) is simultaneously solved by two adjacent values of m, here m * = 4 and m * = 5. A similar situation occurs in the planar case, as described in detail in [19, Figure 5]. In particular, when m * = 4 the folds F bifurcate from W 4 ; as 5 increases towards the codimension-two point 5 4,5 = 360 the secondary bifurcations leading to the folds move down along the W 4 branch and reach zero amplitude when 5 = 5 4,5 . For 5 > 5 4,5 the dominant mode is m * = 5 and the secondary bifurcation to the fold state now takes place on W 5 . As 5 increases this bifurcation moves up along W 5 to a maximum amplitude before moving down again as the next codimension-two point is approached. This process repeats as 5 continues to increase, and ∆ c , the threshold for the onset of the fold state, therefore both oscillates and jumps from branch to branch. This behaviour is shown in figure 8(a) and is similar to that found in the planar case [19, Figure 5]; we expect that the tube problem studied here approaches the planar case once 5 is sufficiently large (sufficiently large tube radius).

Conclusion
In this article, we provided a simple model of an inextensible elastic lining of an innerlined tube subjected to an imposed pressure difference, and described its buckled, wrinkled and folded solutions. We showed that wrinkling is statically generated by a competition between bending, soft-substrate forces and the applied pressure, and explored the limiting behaviour of our model as the strength of the substrate support is reduced eliminating the possibility of wrinkling. We showed that for sufficiently strong substrate support, increasing the applied pressure leads not only to a wrinkleto-fold transition, but also to mixed states. The energies of these states were calculated using weakly nonlinear theory and by numerical continuation for strongly nonlinear solutions. The wrinkle state with wavelength closest to natural is initially the state with the least energy and is thus stable until a single-fold state bifurcates from it. As 5 increases, additional mixed modes arise prior to self-contact, and states with an increasing number of localised folds become possible. The solution profiles match well with observations and resemble structures in growing composite rings [7]. Our approach explains how the wavelength and amplitude of the wrinkles are selected as a function of parameters in pressure-driven wrinkling systems. This is in turn key to understanding, for example, the artery self-cleaning process arising from wrinkledto-unwrinkled cycles triggered by blood pressure changes [20,21] and can be a good starting point for more refined models that include adhesion. A natural question that arises is how the bending modulus, the size of the system and the substrate properties may be optimised to maximise in-plane curvature, thereby optimising the self-cleaning properties for a given pressure jump, while avoiding the wrinkle-to-fold transition. For weaker substrate support the first primary bifurcation is to the m = 2 buckling mode, and the secondary bifurcations move to large amplitudes, beyond the point of self-contact. Thus all bifurcation diagrams simplify and wavenumber of the primary branches increases monotonically with increasing pressure.
Applications of this work to the time-dependent artery problem and to other systems exhibiting competition between buckling, wrinkling and folding will be described elsewhere. LG was funded by grant Conicyt Fondecyt Iniciación 11170700. We thank E. Cerda for valuable discussions.

Appendix A. Derivation from Kirchhoff equations
In equilibrium, the forces acting on an element of the lining can be expressed in terms of the static Kirchhoff equations [12,39]: where F and M are the force and moment acting on the centerline of the element. The extra pressure in square brackets in (A.1) is due to the force per unit of area exerted by the substrate, modelled by a Winkler foundation [26] with a nonlinear quadratic term.
The moment is related to the local curvature of the ring by the constitutive relation M = (B∂ s φ) k, where φ is the angle between a tangent to the ring and fixed horizontal axis, and k is normal to the plane. Thus ∂ s r = (cos φ, sin φ). The unit vectors t and n are given by t = ∂ s r = (cos φ, sin φ) and n = (− sin φ, cos φ); n points towards the interior of the enclosed region. Accordingly, the system can be written in terms of five differential equations: The problem is defined after imposing the closed-curve boundary condition φ (L, t) = φ (0, t) + 2π, L = 2πR, and periodic boundary conditions on ∂ s φ, x, y, F x and F y . We show that this system of equations is equivalent to (3) of the text. For this purpose, we rewrite (A.2) using the constitutive relation M = (B∂ s φ) k and the identity (∂ s φ) k = ∂ s r × ∂ ss r, ∂ s r × (B∂ sss r + F) = 0, which is solved by F = −B∂ sss r + λ∂ s r, where we have introduced the Lagrange multiplier λ(s) to incorporate inextensibility. The latter expression, substituted in (A.1), yields −B∂ sss r + λ∂ ss r + ∂ s λ∂ s r + P − 1 2 K r 2 0 − r 2 n = 0.
Since n and t form an orthonormal basis, the two terms in parentheses must both vanish. From the second, we obtain a differential equation for λ whose solution is where T is a constant. Replacing λ (s) in the first set of parentheses by the above expression, we finally obtain: leading to eq. (3).
The same equation can also be derived from a constrained Lagrangian as done for the planar elastic sheet in [16].

Appendix B. Weakly nonlinear analysis
At each order in the weakly nonlinear analysis, we obtain a linear problem of the form for j = 1, 2, ..., with the first three N j given by To eliminate x j and y j from L[φ j , x j , y j ], we compute (∂ 2 s L + L)[φ j , x j , y j ]: . . . 5 2(∂ s x 0 )(∂ s x j ) + 2(∂ s y 0 )(∂ s y j ) + x 0 ∂ 2 s x j + y 0 ∂ 2 s y j = (∂ 2 s + 1)N j . Expansion of the geometric identities ∂ s x = cos φ and ∂ s y = sin φ now results in where the first three G j are given by Solving (B.1) for j = 1, 2 subject to the requirement that the solution is periodic yields the expressions for φ 1 , A 1 and for φ 2 , A 2 , P 2 , T 2 given in the text. For j even, the solvability condition imposed on G j + (∂ 2 s + 1)N j generates P j (T j ), while for j odd, it generates T j−1 (m, 5 ). Higher order expressions were obtained through symbolic calculations using the software Maple.

Appendix C. Numerical continuation with AUTO
We implemented the problem (1) in AUTO [31] as a 5-dimensional boundary value problem on the domain s ∈ [0, π], representing one half of the lining, with the boundary conditions φ (0) = π/2, φ (π) = 3π/2, x (0) = x 0 , x (π) = x 1 and y (0) = y (π) = 0 together with the force-free conditions φ (0) = φ (π) = 0 [12]. The boundary conditions constrain the rotation symmetry in φ and eliminate translations in y, while the force-free boundary conditions permit reflection in y = 0 to generate solutions on the full circle. A 5-dimensional system with 8 boundary conditions requires 4 degrees of freedom in the parameters [40], so we perform our continuation in (P, T, x 0 , x 1 ). This procedure allows T to adjust to increments in P and the endpoints x 1 , x 2 to change in accordance with the zero-force condition. Figure 3 of the text shows the resulting full circle profiles.
The imposed boundary conditions prevent the computation of asymmetric states F a that are also expected to appear via secondary bifurcations from wrinkled states.

Appendix D. Comparisons
In figure B1, we compare the results from numerical continuation of (3) of the text with the above boundary conditions and the corresponding results obtained above from weakly nonlinear theory carried out to O( 7 ). The values corresponding to the maximum displayed extent of each branch are summarized in table D1. The results demonstrate excellent agreement between perturbation theory and the numerically exact solutions for 1. Equally good agreement is found for the solution profiles as shown in figure C1.