Elasto-capillary necking, bulging and Maxwell states in soft compressible cylinders

Localized pattern formations and"two-phase"deformations are studied theoretically in soft compressible cylinders subject to surface tension and axial loading through several force-controlled loading scenarios. By drawing upon known results for separate yet mathematically similar elastic localization problems, a concise family of analytical bifurcation conditions for localized bulging or necking are derived in terms of a general compressible strain energy function. The effect of material compressibility and strain-stiffening behaviour on the bifurcation point is analysed, and comparisons between our theoretical bifurcation conditions and the corresponding numerical simulation results of Dortdivanlioglu and Javili (Extreme Mech. Lett. 55, 2022) are made. It is then explained how the fully developed"two-phase"(Maxwell) state which evolves from the initial localized bifurcation solution can be comprehensively understood using the simple analytical expressions for the force parameters corresponding to the primary axial tension deformation. The power of this simple analytical approach in validating numerical simulation results for elastic localization and phase-separation-like problems of this nature is highlighted


Introduction
The treatment of localized pattern formation in solid cylinders and hollow tubes as a bifurcation problem has become increasingly prevalent over the last 15 years. A problem which serves as the foundation for this area of research is the localized bulging of a hollow tube subject to the combined effects of axial loading and internal inflation. Despite many experimental observations in the past (Mallock, 1891;Kyriakides and Yu-Chung, 1990), this was only recognized as a bifurcation phenomenon with zero wavenumber under the membrane assumption relatively recently by Fu et al. (2008). Fu et al. (2016) demonstrated that, for a tube of arbitrary thickness, the bifurcation condition for localized bulging is that the Jacobian determinant of the inflation pressure P and the resultant axial force N as functions of the axial stretch and the circumferential stretch on the inner surface must vanish. This is equivalently the condition for an axi-symmetric bifurcation mode with zero wavenumber to exist (Yu and Fu, 2022). Previous linear bifurcation analyses (Haughton and Ogden, 1979a,b) had focussed on periodic axi-symmetric modes, and the zero wavenumber mode was incorrectly thought to correspond to an alternate uniformally inflated state. Since the revelation of Fu et al. (2016), many additional effects such as rotation (Wang et al., 2017), double fibre-reinforcement (Wang and Fu, 2018), bi-layering (Liu et al., 2019) and torsion (Althobaiti, 2022) have been incorporated into the analysis.
The inflation problem has become prototypical in the sense that it often has a very similar mathematical structure to other more complicated elastic localization problems. For instance, through a reformulation of the Jacobian determinant bifurcation condition for the inflation problem, Fu et al. (2018) demonstrated that the bifurcation condition for localized necking in a dielectric membrane under in-plane mechanical stretching and an electric field is that the Hessian of the total free-energy function vanishes. The problem of localized bulging or necking in soft incompressible cylinders and tubes under axial loading and surface tension has also become very well understood as a result of the inflation problem.
Surface tension operates in solids at the elasto-capillary length scale γ/µ, where γ is the surface tension and µ is the ground-state shear modulus (Bico et al., 2018). For extremely soft materials such as gels, creams and biological tissue, surface tension effects can be of the same order of magnitude as the bulk elastic modulus at length scales ranging from tens of nanometres to millimetres. Thus, when modelling the large deformations of these small-scale materials, surface tension effects cannot be neglected. The localized axi-symmetric "beading" of soft slender cylinders, such as axons, has been widely observed experimentally (Matsuo and Tanaka, 1992;Bar-Ziv and Moses, 1994;Fong et al., 1999). The implication of this phenomenon in nerve damage due to traumatic brain injuries (Kilinc et al., 2009) and neurodegenerative conditions such as Alzheimer's and Parkinson's diseases (Datar et al., 2019) has motivated many recent studies surrounding the bifurcation behaviour of an incompressible solid cylinder under a resultant axial force N and a surface tension γ.
The aforementioned elasto-capillary problem was initially studied using non-linear elasticity theory by Taffetani and Ciarletta (2015a,b) and Xuan and Biggins (2016). However, through a weakly non-linear analysis, Fu et al. (2021) was the first to demonstrate that the initial bifurcation is a sub-critical localized necking or bulging solution (depending on the nature of the loading), and this is again associated with zero axial wavenumber. The corresponding bifurcation condition was also shown to take a simple analytical form in this case; for fixed N (fixed γ) and increasing γ (varying N ), localized necking or bulging oc-curs when dγ/dλ z = 0 (dN /dλ z = 0), where λ z is the axial stretch. Xuan and Biggins (2017) and Giudici and Biggins (2020) highlighted that the complete bifurcation process is a phase-transition-like phenomenon which culminates in a "two-phase" state consisting of two regions of distinct but uniform axial stretch λ z connected by a smooth transition zone. The connection between the initial localized bifurcation solution and the final "two-phase" state was explained both theoretically and through numerical simulations by Fu et al. (2021). The beading instability in solid cylinders has since been studied dynamically (Pandey et al., 2021) and through the active strain approach (Riccobelli, 2021). The observations of Fu et al. (2021) for a solid cylinder have also been fully extended to the case of an incompressible hollow tube by Emery and Fu (2021b,c). Circumferential buckling instabilities in cylinders and tubes under surface tension and axial loading (Emery and Fu, 2021a), growth (Bevilacqua et al., 2020), and uniform pressure and geometric everting (Wang et al., 2021) have also been extensively studied in recent years.
The case of a compressible solid cylinder has received very little attention in the literature. This is surprising since, whilst the incompressibility assumption typically makes the bifurcation analysis far easier, soft hydrogels can often possess a large degree of compressibility. Furthermore, in the case of soft biological tissue, whilst incompressibility is often assumed due to the high water content of the material, there is very little supporting experimental evidence for this assumption. Only Carew et al. (1968) has provided evidence that incompressibility is a suitable assumption in modelling arterial tissue. Very recently, however, Dortdivanlioglu and Javili (2022) (hereafter abbreviated as "DJ") analysed the effect of material compressibility via numerical simulations, extending what is already known for the incompressible case. Numerical simulation predictions for the initial bifurcation points are presented for two separate compressible neo-Hookean strain energy functions, and an extensive post-bifurcation analysis tracking the axial propagation of the localized solution is performed. The work offers a different numerical perspective on the existing literature for the incompressible case, but does not present any analytical results to compare with the numerical predictions for the compressible case. The aim of this paper is two-fold. Firstly, we look to extend the known analytical bifurcation conditions of Fu et al. (2021) to the compressible case for a general strain energy function, and make comparisons with the numerical bifurcation conditions in DJ. Secondly, we wish to highlight how the entire post-bifurcation process can be captured through analytical means, and it is hoped that this will encourage future studies to use analytical approaches to guide and corroborate numerical simulation results.
The remainder of this paper is organised as follows. In the next section, we present the primary axial tension deformation and derive corresponding analytical expressions for the resultant axial force N and the surface tension γ. In section 3, analytical bifurcation conditions for localized bulging or necking are derived in terms of a general strain energy function for several force-controlled loading scenarios. The effect of compressibility and strain-stiffening behaviour on the bifurcation point is considered, and comparisons are made with the numerical results of DJ. In section 4, we highlight a powerful analytical approach to describing the full post-bifurcation process of the cylinder, and make comparisons with the corresponding results in DJ. Concluding remarks are finally given in section 5.

Primary deformation
Consider a compressible, isotropic, hyperelastic solid cylinder with a reference configuration B 0 defined in terms of the cylindrical polar coordinates (R, Θ, Z), where (2.1) The finitely deformed configuration B e is in terms of the cylindrical polar coordinates (r, θ, z), and we assume that the solid cylinder undergoes a primary homogeneous deformation of the form where λ θ and λ z are the constant circumferential and axial stretches, respectively. Thus, we have that The position vectors of a representative material particle in B 0 and B e are given respectively by where (E R , E Θ , E Z ) and (e r , e θ , e z ) are the corresponding orthonormal bases. The primary deformation gradient F is then defined through dx = F dX and may be written as Given the one-to-one correspondence between B 0 and B e , we must have that J ≡ det F > 0. The associated left Cauchy-Green strain tensor B = F F T takes the form B = λ 2 θ (e r ⊗ e r + e θ ⊗ e θ ) + λ 2 z e z ⊗ e z , (2.6) and its three principal invariants are expressed through (2.7) We assume that the constitutive behaviour of the material is governed by a strain energy function of the form W = W (I 1 , I 3 ). (2.8) In the computation of our results, we will predominantly specify W to the following compressible Gent material model in order to account for strain-stiffening behaviour: The parameter J m is the material extensibility limit, µ is the ground-state shear modulus and λ = 2ν/(1 − 2ν), where ν ∈ [0, 1/2] is Poisson's ratio. For a fully compressible material, we have that ν → 0, whilst the incompressible limit corresponds to ν → 1/2. In order to facilitate a comparison between our theory and the numerical results of DJ, we will also consider the quadratic and logarithmic compressible neo-Hookean material models given respectively by: (2.11) Note that, in the limit J m → ∞, the compressible Gent material model (2.9) reduces to the quadratic neo-Hookean model (2.10). All of the algebraic manipulations and computations associated with the following work have been performed in Mathematica (Wolfram Research Inc., 2021).

Stress-based formulation
Under the assumption (2.8), the Cauchy stress tensor σ is expressible as where the notation W i = ∂W/∂I i and W ij = ∂ 2 W/∂I i ∂I j for i, j = 1, 2, 3 is employed here and hereafter, and I is the identity tensor. On substituting (2.6) into (2.12), the Cauchy stresses in the radial, circumferential and axial directions are found to take the form These components are constant, and so the equilibrium equations div σ = 0 are automatically satisfied. The solid cylinder is under the combined effect of a surface tension γ and a resultant axial force N . We scale all lengths by R o , all stresses by the ground state shear modulus µ and the surface tension γ by µR o . As such, R o and µ may be set equal to unity without loss of generality; we use the same symbols to denote scaled quantities. The surface tension enters the analysis through the boundary condition (2.14) see Fig. 1. On substituting (2.13) 1 into (2.14), we obtain the following expression for γ in terms of λ θ and λ z : (2.16)

Bifurcation conditions for localization
We are interested in the bifurcation behaviour of the finitely deformed configuration B e . More specifically, we wish to determine the critical value of the load parameter at which localized patterns such as necking or bulging become theoretically possible. We consider three distinct types of loading which may be summarised as follows.
(1) Fixed γ and varying N : A fixed surface tension γ > 0 is applied to the cylinder which is also initially under zero resultant axial force N or a large strictly positive resultant axial force (due to the application of a dead weight at an end of the tube, say). In the former case, fixed surface tension produces an axial compression such that λ z < 1 initially, and from this point we may increase N monotonically from zero to trigger bifurcation. In the latter case, we can instead decrease N monotonically from its large starting point determined by the dead weight in order to incite a bifurcation. We refer to these cases as axial force controlled "loading" and "unloading", respectively.
(2) Fixed λ z and increasing γ: The total length of the tube is fixed both before and after bifurcation into a localized necking or bulging solution has taken place. In other words, the averaged axial stretch λ z ≥ 1, defined generally as the deformed length of the cylinder divided by the undeformed length, is fixed. The surface tension γ is then increased monotonically from zero to prompt a bifurcation.
(3) Fixed N and increasing γ: The resultant force N ≥ 0 is fixed, and this induces an initial axial stretch λ z ≥ 1. The surface tension γ is then increased monotonically from zero to trigger bifurcation, and the axial stretch λ z will decrease monotonically from its aforementioned starting value in tandem to preserve the constant value of N .

Fixed γ and varying N
For any fixed γ ≥ 0, we can define the circumferential stretch λ θ as an implicit function of the axial stretch λ z through (2.15). Based on the analysis of Fu et al. (2021) and Emery and Fu (2021b,c) for the incompressible solid cylinder and hollow tube case (respectively), we may then conjecture that the bifurcation condition for elasto-capillary localization is dN /dλ z = 0, where N is given in (2.16) and γ is fixed in the differentiation. This bifurcation condition takes the form where λ θd = dλ θ /dλ z , and the right-hand side of (3.1) is evaluated at the critical value of λ z , λ zcr . By differentiating equation (2.15) implicitly with respect to λ z , the following explicit expression for λ θd in terms of λ z and λ θ can be obtained: For a given fixed γ, the bifurcation values λ zcr can be determined numerically from (3.1). Then, with use of (2.16), the associated bifurcation values of the resultant axial force, N cr ≡ N (λ zcr ), may be computed. In Fig. 2, for the compressible Gent model (2.9) with J m = 100, we plot the resultant axial force N against λ z for several fixed γ ≥ 0 (blue curves), as well as the bifurcation criterion N cr = N (λ zcr ) (black curves), with (a) ν = 0.05 and (b) ν = 0.25. We observe that, as in the incompressible case, there exists a minimum value of γ, γ min , below which the N = N (λ z ) curve is monotonic increasing and localization cannot occur. Here, the value of γ min is dependent on the value of ν. For any fixed γ > γ min , there exists two bifurcation values of λ z , λ L zcr and λ R zcr > λ L zcr , which are located at the local maximum and minimum of the now non-monotonic N = N (λ z ) curve. In the limit γ → γ min , these two bifurcation values coalesce into a single value λ z = λ min , and the maximum and minimum of N coalesce into an inflection point (marked by the black cross). When "loading" from N = 0, the local maximum of N is the bifurcation point of interest, and we expect from the weakly non-linear analysis of Fu et al. (2021) and Emery and Fu (2021c) that localized necking will be triggered when this point is reached. In contrast, when "unloading" from some large strictly positive N , the local minimum of N is the relevant bifurcation point, and we expect that localized bulging will initiate here. We observe in Fig. 2 that a larger value of γ > γ min will delay the expected onset of localized necking when "loading" from N = 0, but incite the expected onset of localized bulging when "unloading" from some large strictly positive N . In Fig. 3 (a), we examine the variation of γ min with respect to ν for the Gent material model (2.9) with several fixed values of J m . We observe that γ min increases with both ν and J m . Thus, for materials with a greater level of compressibility (i.e. for values of ν closer to zero), or a lower level of extensibility, there is a greater range of values of γ for which localization can occur in this loading scenario. In (b), we plot this same relationship for the quadratic (solid dark blue curve) and logarithmic (solid light blue curve) neo-Hookean material models, and compare with the numerical simulation results presented in Fig. 6 (a) of DJ (squares). We note the exceptional agreement between both sets of results, and in the incompressible limit ν → 1/2, we recover the value γ min = 4 √ 2 which was reported in Fu et al. (2021).

Fixed λ z and increasing γ
When λ z ≥ 1 is fixed and γ is increased gradually from zero, the bifurcation condition for localization is also dN /dλ z = 0. Hence, it is equivalent to (3.1) except the left hand side is replaced with γ = γ cr here and the right hand side is evaluated at the fixed axial stretch rather than λ zcr . In Fig. 4, we plot γ cr against (a) ν for J m = 100 and several fixed values of λ z ≥ 1 and (b) λ z for ν = 0.4 and several fixed values of J m . In (a), we observe that highly compressible cylinders are the least susceptible to localization. We also observe that there exists a threshold value of ν below which a greater fixed axial stretch suppresses localization, and above which a greater fixed axial stretch encourages localization. In (b), we show that there exists a threshold value of λ z (λ z = 1.225 in the case presented) below which a greater extensibility limit suppresses localization, and above which a greater extensibility limit encourages localization. We observe that γ cr as a function of λ z also possesses a minimum. This behaviour was similarly observed in the incompressible case studied by Fu et al. (2021), and it was shown that the initial bifurcation solution was localized necking (localized bulging) for λ z < λ min (λ z > λ min ). Here, the value of λ min will vary with ν, and we plot this relationship in Fig. 5. We see that, for all values of J m considered, the value of λ min increases with ν. Thus, for cylinders with a greater degree of compressibility, there is a larger range of values of fixed λ z for which an initial localized bulging solution will emerge at γ = γ cr . In Fig. 5 (b), we demonstrate the exceptional agreement between our theoretical results and the numerical simulation results in Fig. 6 (b) of DJ for the quadratic and logarithmic neo-Hookean models. In Fig. 6, we plot the variation of γ cr with respect to ν for the quadratic and logarithmic neo-Hookean models with λ z = 1, and we note again that there is exceptional agreement between our theoretical results (solid curves) and the numerical simulation results given in Fig. 5 of DJ (squares).

Fixed N and increasing γ
By making appropriate rearrangements in (2.16), the surface tension γ can be expressed explicitly in terms of λ θ and λ z for any fixed N ≥ 0 as follows: Then, we may subtract (3.3) from (2.15) and define λ θ as an implicit function of λ z from the resulting equation. The bifurcation condition for localization is then conjectured to be dγ/dλ z = 0, where γ is given in (3.3) and N is fixed in the differentiation. Explicitly, this condition is expressible as

4)
and this equation is also evaluated at λ z = λ zcr . The expression for λ θd can again be obtained by differentiating (2.15) implicitly with respect to λ z . The left hand side of the resulting equation will vanish since dγ/dλ z = 0 is the bifurcation condition in this loading scenario, and we thus recover the relation presented in (3.2). Once the critical stretch λ zcr has been obtained from (3.4), we can substitute this value into the equation (3.3) to obtain the corresponding critical surface tension γ cr ≡ γ(λ zcr ). In Fig. 7, we plot the function γ = γ(λ z ) given in (3.3) for (a) ν = 0.25, J m = 100 and several fixed N ≥ 0, and (b) N = 8, J m = 100 and several fixed ν. As in the purely incompressible case, there is seen in (a) to be a minimum value of N , N min , below which the curve γ = γ(λ) is monotonic decreasing and localization is prohibited. This minimum value of N will depend on the value of ν. Above this minimum value of N , two bifurcation values of λ z , λ L zcr and λ R zcr > λ L zcr , emerge and correspond to the local minimum and maximum of the now non-monotonic γ = γ(λ z ) curve, respectively. When fixing N > N min with γ = 0 initially, an initial axial stretch λ z > λ R zcr will be produced. As we then increase γ monotonically from zero, we will approach λ z = λ R zcr from the right (recall that an increase in γ must be accompanied by a decrease in λ z to preserve the constant N ). Thus, the bifurcation point of interest in this loading scenario is situated at the local maximum of γ = γ(λ z ), and we expect that a localized bulge will initiate at this point given the findings of Fu et al. (2021) for the incompressible case.
In Fig. 7 (a), we observe that larger fixed N above N min correspond to larger values of γ cr (marked by the black dots), and so a greater fixed axial force will discourage localized bulging when the material is compressible. In (b), we find that there exists a maximum value of ν, ν max , above which the γ = γ(λ) curve becomes monotonic decreasing and localization becomes impossible. The value of ν max will vary with the value of the fixed N . We observe also that, for smaller values of ν (i.e. for materials with a greater level of compressibility), the associated value of γ cr is larger. Hence, increased compressibility discourages the initiation of a localized bulge in this loading scenario. In Fig. 8, we plot the variation of (a) N min against ν and (b) ν max against N . In (a), for any given Poisson ratio ν, localized bulging is only possible if the fixed axial force N is greater than the value of N min given by the appropriate blue curve. For all values of J m considered, the value of N min increases with ν. Thus, for materials with a greater degree of compressibility, there is a greater range of values of fixed N for which localized bulging can occur. We note also that, in the incompressible limit ν → 1/2, we recover the result N min = 9π/2 2/3 in the limit J m → ∞ which was originally given in Fig. 5 (b) of Fu et al. (2021). In (b), for any given fixed value of N , localized bulging is only possible provided the value of the Poisson's ratio is less than the value ν max on the appropriate blue curve shown. We observe that, for each value of J m , there exists a value of N below which localized bulging is prohibited in cylinders with any degree of compressibility. For instance, in the limit J m → ∞, localized bulging is impossible in any compressible cylinder if N < 5.825.

Post-bifurcation behaviour
Our analysis in the previous section provided a generalized analytical framework to support the numerical simulation predictions of DJ for the bifurcation points corresponding to elasto-capillary localization in compressible solid cylinders. In this section, we go one step further by showing that the post-bifurcation behaviour of the cylinder can be comprehensively described with the aid of our analytical expression (2.16) for the primary deformation. The importance of such an analytical approach in corroborating post-bifurcation results from numerical simulations is also highlighted. We focus primarily on the fixed λ z and increasing γ scenario, since this is arguably the least well understood case. To set the scene, we first review what is already known for the case of an incompressible hollow tube where the radial displacement of the inner surface is zero. This case was well studied in Emery and Fu (2021b,c), and we note that an incompressible solid cylinder is recovered when the fixed inner radius R i tends to zero. Finite Element Method simulations conducted in Abaqus (2013) are presented in Fig. 9 (solid blue curves) for the incompressible Gent material model with J m = 100, R i /R o = 0.4, (a) λ z = 1.5 > λ min and (b) λ z = λ min ≈ 1.16. The initial bifurcation occurs when dN /dλ z = 0, and we present this condition through the dashed blue curves. Through a weakly non-linear analysis, the initial bifurcation solution for λ z = 1.5 > λ min was determined in Emery and Fu (2021c) to emerge sub-critically and to take the explicit form of a localized bulge. The numerical simulations then showed that, if we attempt to increase γ beyond its bifurcation value γ cr , a snap-through to a "two-phase" state occurs. This configuration consists of an axially propagated bulge "phase" with constant stretch λ zL and a depressed "phase" with axial stretch λ zR < λ zL ; these "phases" are connected by a smooth yet sharp transition zone as shown in Fig. 9 (c). Note that the bulged "phase" could either be centred at z = 0 or situated as the two ends of the tube (depending on how the small amplitude imperfection in the simulations is introduced), and that the overall averaged axial stretch of the "two-phase" state remains fixed at 1.5. The stretches λ zL and λ zR are given by the left and right branches of numerical simulation curves in (a) and (b). In the case shown in (b) where λ z = λ min , an exceptionally super-critical bifurcation takes place when the bifurcation value γ = γ min is reached in which the tube evolves smoothly into the same "two-phase" state as in the λ z = 1.5 case without any initial localization. The difference between the two cases covered is that the overall averaged axial stretches are different, meaning that the proportion λ zL /λ z of the propagated bulge "phase" with respect to the overall fixed length of the tube will differ; see (c).   λ z = λ min λ z = 1.5 Figure 9: FEM simulation results (solid blue curve) for the case of a hollow tube with fixed inner radius R i = 0.4, J m = 100, (a) fixed λ = λ min and (b) fixed λ = 1.5. The dashed blue curves represent the theoretical bifurcation condition, and the black squares give the relationship between the surface tension γ and stretches λ L and λ R determined from (4.1). The black dots mark the bifurcation point in each case given by the simulations. In (c), we present the "two-phase" configuration of the tube for fixed λ = 1.5 and λ min when the surface tension has been increased beyond its bifurcation value to γ = 9. Both configurations consist of a bulged section with uniform axial stretch λ zL ≈ 0.59, in between two depressed sections with stretch λ zR ≈ 2.25. The proportion of the bulged "phases" differ in each case due to the different averaged axial stretches.
Numerical simulations aren't, however, necessary to determine λ zL and λ zR ; they can be determined with the aid of the analytical expression for N = N (λ z ). To elaborate, as was first elucidated by Clerk-Maxwell (1875), the stretches for each "phase" can be defined implicitly as functions of γ through the following equal area rule: (4.1) The resultant axial force N M W for the "two-phase" state is invariant across the bulged and depressed "phase", but will differ for each γ > γ cr . In the aforementioned hollow tube case, exceptional agreement has previously been shown between the values of the Maxwell stretches produced through FEM simulations and Maxwell's equal area rule; see Emery and Fu (2021b) and Fig. 9 (b).
We now turn our attention to the compressible cylinder case. In Fig. 10 (a) and (c), we present theoretical results for the bifurcation values of λ z (dashed curve) and the Maxwell stretches (solid curve) as functions of γ for the quadratic and logarithmic neo-Hookean models, respectively. The black stars are numerically simulated bifurcation points taken from Fig. 11 (d) of DJ, and we observe perfect agreement with our theory. However, the black dots, which are the numerically simulated Maxwell stretches for different values of γ from DJ, are at odds with our theoretical predictions. The basic principle of Maxwell's equal area rule in the present context is that the stretches λ zL and λ zR should be defined such that the magnitude of the areas between the horizontal line passing through the points (λ zL , N M W ) and (λ zR , N M W ), and the curve N = N (λ z ) above and below the horizontal line, should be equal. Only if this condition is satisfied can a "two-phase" state exist. We show in Fig. 10 (b) and (d) that this requirement is satisfied by our Maxwell stretches, but not by those predicted in DJ. Thus, it is clear that the analytical approach presented is a powerful tool in guiding numerical simulation studies of elastic phase-separation-like phenomena, and in validating the results of such studies. In Fig. 11, we plot the Maxwell stretches λ zL and λ zR as functions of γ obtained from the equal area rule for the compressible Gent material model (2.9). We set J m = 100 and fix ν at several different values. The right-most curve corresponds to ν → 1/2 (i.e. the incompressible limit), and to provide validation of our theoretical results, we compare with the corresponding FEM simulation results presented in Fig.  13 (a) of Fu et al. (2021) (black squares). We observe that there is exceptional agreement between the two sets of results, emphasising the point that the equal area rule should act as a consistency check for numerical simulation results (and vice versa).
We highlight the remarkable capability of the equal area rule in describing the postbifurcation behaviour of the compressible cylinder using only analytical expressions for force parameters corresponding to the primary deformation. Not only does the approach yield the Maxwell stretches for each portion of the fully developed "two-phase" state, it can also fully predict the nature of the axial propagation of the bulged or depressed phase. For example, in the fixed λ z and increasing γ scenario covered previously, say we fix λ z > λ min . Then, we understand that a localized bulge will emerge at γ cr , followed by a snap-through to a "two-phase" state comprising of an axially propagated bulged "phase" with stretch λ zL and a depressed "phase" with stretch λ zR . The longitudinal proportion of the bulged "phase" with respect to the overall "two-phase" state can easily be computed through λ zL /λ z . Now, from our results in Fig. 11, we observe that λ zL is generally a decreasing function of γ. Thus, as we increase γ beyond its bifurcation value, the proportion of the propagated bulge "phase" with respect to the overall length of the cylinder will also decrease. In other words, an increasing surface tension in the fully non-linear regime will cause a reduction in the axial length of the bulged "phase". We note, however, that dλ zL /dγ is typically small, and tends  Figure 10: In (a) and (c), the critical stretches λ zcr from our theoretical bifurcation condition (dashed blue curve) and the Maxwell stretches λ zL and λ zR computed using the equal area rule (solid blue curve) are presented in the (λ z , γ) plane for ν = 0.4, and the quadratic and logarithmic neo-Hookean material models, respectively. The black stars and dots give the numerical simulation results from DJ for the bifurcation points and the Maxwell stretches, respectively. In (b) and (d), we superpose our theoretical bifurcation points and Maxwell stretches as well as the numerically determined values of DJ on the N = N (λ z ) curve for γ = 6. This demonstrates that the theoretically determined Maxwell stretches satisfy the equal area rule, whereas the numerically simulated stretches don't.
to zero as γ gets sufficiently large. Hence, the reduction of "bulged" phase length as γ > γ cr is increased is very gradual, and when γ gets large enough, the length of the bulged "phase" will seemingly approach a non-zero limiting value.

Concluding remarks
The complete bifurcation behaviour of an incompressible solid cylinder under axial loading and surface tension is fully understood. However, with the exception of the numerical study of DJ, analogous studies when the cylinder is compressible are scarce. In this paper, we have provided greater theoretical insights into the bifurcation behaviour of compressible solid cylinders and a source of comparison for existing and future numerical studies. By drawing upon known results for the mathematically similar problem of localized bulging in hollow tubes under axial loading and internal inflation, we derived analytical bifurcation conditions for elasto-capillary localized bulging or necking in compressible solid cylinders for three distinct loading scenarios. For the quadratic and logarithmic neo-Hookean material models, we found perfect agreement between our theoretical bifurcation conditions and the numerical simulation conditions presented in DJ. Based on the findings for an incompressible solid cylinder and hollow tube, we then conjectured that, when the axial stretch λ z is fixed and the surface tension γ is increased beyond its bifurcation value, the tube undergoes a phase-separation-like process where the initial localized solution evolves into a "two-phase" state. We explained how, with the use of the simple analytical expression (2.16) for N , the Maxwell stretches λ zL and λ zR associated with this "two-phase" state can be determined as functions of γ through the equal area rule. On comparing the results of the equal area rule approach with the corresponding numerical simulation results in Fu et al. (2021) and DJ, we found perfect agreement in the former case but disagreement in the latter. This highlighted that numerical studies of phase-separation-like phenomena don't always use the equal area rule as a consistency check on their results, and it is hoped that this paper will invoke change in this regard in future studies.