Equilibrium configurations of hard spheres in a cylindrical harmonic potential

TL;DR

This paper investigates the equilibrium structures of hard spheres in a cylindrical harmonic potential, combining experimental observations in a rotating liquid-filled tube with theoretical modeling to understand buckled configurations and bifurcations.

Contribution

It provides a combined experimental and theoretical analysis of buckled structures of hard spheres in a cylindrical potential, highlighting bifurcation phenomena and stable configurations.

Findings

Experimental observation of buckled structures in a rotating liquid-filled tube.

Theoretical model predicts a wide range of bifurcating structures.

Stable structures observed both experimentally and in theory.

Abstract

A line of hard spheres confined by a transverse harmonic potential, with hard walls at its ends, exhibits a variety of buckled structures as it is compressed longitudinally. Here we show that these may be conveniently observed in a rotating liquid-filled tube (originally introduced by Lee et al. [T. Lee, K. Gizynski, and B. Grzybowski, Adv. Mater. 29, 1704274 (2017)] to assemble ordered three dimensional structures at higher compressions). The corresponding theoretical model is transparent and easily investigated numerically, as well as by analytic approximations. Hence we explore a wide range of predicted structures occurring via bifurcation, of which the stable ones are also observed in our experiments. Qualitatively similar structures have previously been found in trapped ion systems.

Equations

$$\Delta =(Nd-L)/d=N-L/d,$$

$$F_n$$

$$=G_0(\tan {\theta }_n+\tan {\theta }_{n+1}).$$

$$r_n+r_{n+1}$$

$$F_n+F_{n+1}$$

$${\theta }_{n+1}$$

$$F_{n+1}$$

$$E=\frac{E_{\text{rot}}}{m{\omega }^2{d}^2}=\frac{1}{2}\sum _{n=1}^N{r}_n^2,$$

$$\Delta =N-\sum _{n=1}^N\cos {\theta }_n.$$

$$E_S$$

$$+\left[{\left(\frac{{\delta }_1}{d}\right)}^2+{\left(\frac{{\delta }_N}{d}\right)}^2\right]).$$

$$\left(\begin{array}{c}{F}_n\\ {\theta }_n\end{array}\right)=\left(\begin{array}{cc}\frac{1}{G_0}-1& -1\\ \frac{1}{G_0}& -1\end{array}\right)\left(\begin{array}{c}{F}_{n-1}\\ {\theta }_{n-1}\end{array}\right).$$

$$\left(\begin{array}{c}{F}_n\\ {\theta }_n\end{array}\right)={\left(\begin{array}{cc}3+ϵ& -1\\ 4+ϵ& -1\end{array}\right)}^{n-1}\left(\begin{array}{c}{F}_1\\ {\theta }_1\end{array}\right).$$

$$F_n$$

$${\theta }_n$$

Full text

Granular systems Complex fluids and colloidal systems Theory of crystal structure, crystal symmetry; calculations and modeling

\nobibliography*

1 Introduction

Particles confined in the vicinity of a straight line by a transverse potential have been found to exhibit a rich variety of structures, particular distorted linear chains [1]. In detail they depend on the interactions between the particles, the confining potential and any boundary conditions at the end of a finite sample. Previous observations have been made with ion traps [2, 3, 4, 5, 6, 7, 8], but also finite dust clusters [9], overdamped colloidal systems [10] and microfluidic crystals comprising of droplets [11]. The complex scenario of the appearance of such structures, induced by bifurcations upon increasing compression, has been sketched by Landa et al. [12] for the case of quadripolar confining potential and Coulomb interactions.

Here we introduce a much more elementary experimental system for such an investigation. Our theoretical model is amenable to an analysis using simple numerical and analytical methods.

The experimental system consists of $N$ equal-size hard polypropylene spheres in a horizontal liquid-filled tube, rotating in a lathe. The spheres are buoyant, so that a centripetal force drives them towards the central axis; the rotational speed is high enough to make gravity negligible. This set-up was first introduced by Lee et al. [13] in order to determine the equilibrium phases (ordered three dimensional structures) over a wide range of dimensionless compressions $\Delta$

[TABLE]

where $d$ is the sphere diameter and $L$ is the tube length. The present application is at very low compression, where the linear chain of contacting spheres is observed to buckle. In this range, the structures are all planar. This is assumed in the stepwise method described below, but not in the method of energy minimisation, which indeed finds planar configurations.

The rich scenario of structural transitions within this regime under increasing compression is predicted in detail by the analysis provided below. It uses two numerical methods: an iterative stepwise solution for force equilibrium positions and a simulation based on energy minimisation.

2 Theory and numerical analysis

2.1 Iterative stepwise method

In the experimental system, detailed below, each sphere of mass $m$ experiences a centripetal force $f_{c}=m\omega^{2}R$, where $R$ is the distance of its centre from the central axis of the tube and $\omega$ the rotational speed. For the case of no compression ($\Delta=0$) the spheres align in a linear chain along the central axis. At a finite compression ($\Delta>0$) the chain starts to buckle and the equilibrium structures take the form of modulated zig-zag structures as illustrated in Fig 1. We developed an elementary stepwise method to describe such structures for low energies.

In the following we will use the dimensionless distance from the central axis for each sphere $r=R/d$ and the dimensionless centripetal force $F=f_{c}/(m\omega^{2}d)=r$.

Our aim is to calculate the dimensionless forces $F_{n}=r_{n}$ and tilt angles $\theta_{n}$, as defined in Fig 1, for $n=1$ to $N$ spheres. Considerations of force equilibrium and geometrical equations yield iterative relations for $F_{n}$ (or $r_{n}$) and $\theta_{n}$ as follows.

The compressive forces $G_{n}$ between contacting spheres, are given by $G_{n}\cos\theta_{n}=G_{0}$ from the condition of force equilibrium in the $x$ direction, which is that of the central axis. $G_{0}$ is the magnitude of the compressive force at each end of the system.

The equilibrium of centripetal forces $F_{n}$ on the $n$th sphere gives

[TABLE]

The centres of contacting spheres are separated by their diameter. Hence the radial distances and forces are

[TABLE]

The above equations then relate $\theta_{n+1}$ and $F_{n+1}$ to $\theta_{n}$ and $F_{n}$, i.e.

[TABLE]

These equations may be used in a “shooting method” to find solutions for a specified value of $G_{0}$. The hard-wall boundary condition for sphere $n=1$ requires the first tilt-angle $\theta_{1}$ to be zero, with an arbitrary $F_{1}$. Using eqs. (4) we proceed iteratively to $(F_{N+1},\theta_{N+1})$. The angle $\theta_{N+1}$ corresponds to the contact of the $N$th sphere with the wall, as illustrated in Fig 1.

We search for values of $F_{1}$ (in general more than one) such that the angle $\theta_{N+1}$ is zero, satisfying the second hard-wall boundary conditions. This search is performed by coarse graining the initial force $F_{1}$ over a range of $0<F_{1}\leq 0.01$ in steps of $10^{-4}$. These values are then used as brackets in a bisection method.

The non-dimensional total energy $E$ of such a hard sphere structure can be calculated as

[TABLE]

where as in [14] we have omitted the (constant) energy contribution due to the moment of inertia of the spheres. The compression $\Delta$ from eq. (1) is given by

[TABLE]

By performing the root search at varying compressive force $G_{0}$, we can accumulate a data set, for which we can calculate the energies and compressions in this way.

In practice we encounter difficulties with the stepwise method beyond a compression of $\Delta\geq 0.9$ (for $N=20$). Above this point our implementation of the bisection search method has problems to find all solutions. We expect to successfully extend the application of the stepwise method to this regime in the future.

2.2 Simulations based on energy minimisation

To confirm and supplement the results of the stepwise method we also seek equilibrium configurations using energy minimisation starting from random configurations. These simulations are more general than the stepwise method since they are not restricted to be planar. Results for larger compressions can also be generated. We have used energy minimisation on a system of soft spheres, extrapolating to the limit of hard spheres, in order to corroborate the results of the stepwise calculations.

In the soft sphere model, the overlapping spheres repel each other according to Hooke’s law with a spring constant $k$. This crude formalism is often used in the context of foam structures and rheology [15].

The non-dimensional total energy $E_{S}$ for $N$ soft spheres (of diameter $d$), longitudinally confined between length $L$ is given by,

[TABLE]

The first term is the rotational energy of each sphere. The second term accounts for the overlap between any two spheres, where the overlap between spheres $n$ and $m$ is defined as $\delta_{nm}=|\bm{R}_{n}-\bm{R}_{m}|-d$, where $\bm{R}_{n}$ and $\bm{R}_{m}$ are the centre positions of two contacting spheres. The final term accounts for the overlaps $\delta_{1}$ and $\delta_{N}$ of the two end spheres with the two boundaries.

For any given values of compression $\Delta$ and $k/m\omega^{2}$ we find equilibrium solutions (stable or metastable) by varying the coordinates of the sphere centres. Finally, by performing a series of simulations with increasing values of $k/m\omega^{2}$ we can extrapolate to the hard sphere case (i.e. $k/m\omega^{2}\rightarrow\infty$) and compare directly with the stepwise method.

The solutions from the stepwise method are only in force equilibrium, i.e. they can be stable or unstable solutions. We have also used energy minimisation to check the stable/unstable character of the solutions.

3 Numerical results

We present results for a variety of low energy structures. Many other equilibrium structures could be found with higher energy, if we were to extend the range of our search parameters. These may be of limited physical significance and will not be pursued in the present paper.

3.1 Typical profiles

For low compressions our search yields only one structure that we will refer to as the symmetric structure S, since the profile for $F_{n}$ (or displacement $r_{n}$) is symmetric around the midpoint of the system. (Note that we have defined $F_{n}$ and $r_{n}$ to be positive.)

Examples of such a profile for $N=20$ are presented in Fig 2 for a low (green triangles) and high (blue stars) compression where we show the displacement $r_{n}$ from the central axis vs the (dimensionless) position $x_{n}~{}=~{}1/2~{}+~{}\sum_{i=2}^{n}~{}\cos(\theta_{i})$.

These results show perfect agreement with the symmetric structure generated by energy minimisation and extrapolated to the hard sphere limit. The structures obtained by energy minimisation are necessarily confined to stable cases.

For high compressions additional asymmetric structures are obtained from the stepwise method. An example for the displacement profile for such a structure is given by the red crosses in Fig 2.

3.2 Bifurcation diagrams

We have used the iterative stepwise method to search for low energy structures in the range of the compressive forces between $0.2\leq G_{0}\leq 0.25$ and initial forces between $0<F_{1}<0.01$. These structures correspond to relative compression below $\Delta<0.9$. They were computed for both an even ($N=20$) and an odd ($N=19$) number of spheres, for which the results differ qualitatively.

In these parameter ranges the total energy $E_{\text{Symm}}$ of the symmetric structure S increases from [math] to roughly $0.15$, whereas the difference between the energies of alternative structures is only of order $10^{-4}$. We therefore computed the energy $\Delta E=E-E_{\text{Symm}}$ relative to that of the symmetric structure at the same compressive force $G_{0}$ and plotted them against their compression for even (Fig 3(a)) and odd case (Fig 3(b)).

We present these relative energies in the vicinity of the compression range where the first asymmetric structures are created by bifurcation. While both cases of even and odd $N$ feature an increasing number of bifurcations as compression is increased, they are qualitatively different and will thus be discussed separately.

For the even case $N=20$ an increasing number of asymmetric structures (A-F) are introduced by bifurcation as compression increases. The unstable structures are marked with an asterisk. The first two additional branches A* and B emerge from an “out-of-the-blue” bifurcation at $\Delta=0.558$ without any preceding structure. Of these two branches, structures on branch B are stable, whereas structures from A* are unstable, as verified by energy minimisation.

Two further structures C* and D, appear via a pitch-fork bifurcation out of the previous stable structure B at $\Delta=0.588$. B and the additional branch of lower energy D are stable, whereas the upper one C* is not. A similar pitch-fork bifurcation of the D branch occurs for the next two structures at $\Delta=0.622$, from which the lower branch F is again stable and E* unstable.

Examples of structures from all of the seven branches for the even case in Fig 3(a) are given in Fig 4. The vertical black solid line in these plots represents the midpoint of the structure; the vertical red dashed line indicates the peak position of the sphere profile, as estimated by a quadratic fit to the sphere positions around the maximum. For unstable structures the peak position coincides roughly with the centre position of a sphere. Note the degeneracy: asymmetric structures may have the peak left or right of the centre.

The energy diagram for the odd case of 19 spheres (in Fig 3(b)) differs with respect to the first bifurcation. Here only a single new stable structure (branch B) emerges. From then on bifurcations follow the pattern of the even case, in which previous structures remain stable and new structures of lower energy are stable (i.e. D and F are stable, while C* and E* are unstable).

While the structures that we have identified here appear to be the only equilibria within the specified range in energy and compression, structures with a more complicated profile occur at higher energy, which we have not addressed here. The displacement profiles of these structures can contain two or more off-centred peaks.

3.3 Maximum angles

We have also computed the maximum angle $\theta_{\text{max}}$ of the symmetric structure with varying compression for the stepwise method and energy minimisation, see Fig 5. This is a quantity that can readily be extracted from experimental data, see below.

While our results for the stepwise method stop at a compression of $0.9$, the maximum angle $\theta_{\text{max}}$ from the energy minimisation was computed up to a compression of $\Delta\lesssim 1.3$. At this point the modulated zigzag structure acquires an additional contact with the next-nearest neighbour sphere.

At low compressions ($\Delta\lessapprox 0.1$), where the displacement profile is of the type shown by the green triangles in Fig 2, $\theta_{\text{max}}$ varies as $\sqrt{\Delta}$.

4 Linear approximation

In order to better understand the above results, we have developed an approximate, linear analytic description as follows. For small angles $\theta_{n}$ and forces $F_{n}$, linearisation of eq. (4) leads to

[TABLE]

Recursive substitution of $F_{n}$ and $\theta_{n}$ and setting $\frac{1}{G_{0}}~{}=~{}4~{}+~{}\epsilon$ for $\epsilon$ small and positive, results in:

[TABLE]

The largest possible value for the compressive force is $G_{0}=1/4$, which, for the case of an infinitely long chain, corresponds to the uniform zigzag structure ($r_{n}=-r_{n-1}$).

A solution for $(F_{n},\theta_{n})^{T}$ may be expressed in terms of eigenvalues $\lambda_{1,2}$ and eigenvectors $\bm{V}_{1,2}$ of the above matrix as $(F_{n},\theta_{n})^{T}=a\lambda_{1}^{n-1}\bm{V}_{1}+b\lambda_{2}^{n-1}\bm{V}_{2}$. To lowest order in $\epsilon$ the eigenvalues are given by $\lambda_{1,2}=1\pm\sqrt{\epsilon}$ with the corresponding eigenvectors $\bm{V}_{1,2}=\left(\frac{1}{2}(1\pm\frac{\sqrt{\epsilon}}{2}),1\right)^{T}$. The prefactors $a$ and $b$ are obtained from the initial conditions $F_{1}$ and $\theta_{1}$.

The solution in the linearised approximation for $\theta_{1}=0$ is then given by

[TABLE]

with the offset in the forces $\phi~{}=~{}\operatorname{arctanh}(\sqrt{\epsilon}/2)$. Note that $\theta_{n}$ does not have an offset, since $\theta_{1}=0$.

A comparison of angles $\theta_{n}$ using the approximated linearised equation (11) and the previously numerical exact stepwise method is shown in Fig 6 for a compressive force of $G_{0}=0.234$, resulting in a compression of $\Delta=0.500$. The starting value $F_{1}$ in the linearised scheme was taken from the corresponding value in the stepwise method.

We find excellent agreement between the linearised theory and the stepwise method up to about $n=8$. The linear theory produces a monotonically increasing function (Fig 6), whereas the accurate solution “rolls over” and decreases towards the second boundary. This can be understood in terms of the role of nonlinearity, and approximated in an ad hoc manner: we will leave this to a subsequent paper.

5 Comparison with experiment

Our experimental procedure is similar to that of Lee et al. [13]. We placed an even number of $N=34$ polypropylene beads of density \rho=0.900\,\text{g/cm{}^{3}} and diameter $d=3.000\pm 0.001\,\text{mm}$ [16] in a cylindrical tube (inner diameter $15.91\pm 0.01\,\text{mm}$; outer diameter $20.17\pm 0.01\,\text{mm}$; length $130.55\pm 0.01\,\text{mm}$) filled with water (density \rho_{w}=1\,\text{g/cm{}^{3}}).

The tube is sealed on both ends with stoppers, making sure that no air bubbles remain within the system. The extent to which the stoppers intrude into the tube can be varied, allowing us to adjust the compression $\Delta$.

The tube is then mounted onto a commercial lathe (Charnwood W824), for which we set the rotation frequency to $\omega=1800\pm 50\,\text{rpm}$. In order to record the structures we used a stroboscopic lamp, whose frequency is matched to that of the lathe. A slight off-set between both frequencies is used so that recorded structures appear to be slowly rotating (see example in the supplemental video LatheExperiment.mp4).

Fig 7 shows images of the structures that we identified with the branches S, B, D, and F by comparing them with to the numerical results from 4. The experimental structures can be identified by the distance of the peak position to the structure centre. This is independent of $N$ for large number of spheres because the wall effects can be neglected. The identified structures correspond to all the stable structures of Fig 3. Structure S, as well as structure B, were found at a compression of $\Delta=0.44\pm 0.02$, while the compressions for structure D was $\Delta=0.59\pm 0.02$ and for F, $\Delta=0.68\pm 0.02$.

However, in order to reconcile these experimental results with the theoretical predictions of previous sections, it is necessary to introduce an effective diameter for the spheres, about $1\,\%$ greater than the true value. This increases the effective compressions by a constant shift of roughly $0.35$. We attribute it to the effects of vibration of the lathe, and will explore strategies for its mitigation in future work. This shift also features in previous results from Lee et al. [13, 14].

We extracted the maximum angle $\theta_{\text{max}}$ for the symmetric structure with varying compression for the experiments (see Fig 5). It shows very clearly the necessity for the adjustment of sphere diameter. Due to the neglecting wall effects, these results only depend on $N$ being odd or even for a large enough number of spheres.

6 Conclusion

The compressed and confined sphere chain presents a variety of fascinating observations, previously described in terms of “kinks” or “solitons” [6]. We have succeeded in exploring many of its properties, using simple apparatus and theoretical methods.

Recently we have found a yet simpler experimental method which should be useful, at least for purposes of demonstration. It consists of a horizontal tube into which ball bearings are introduced. Slight agitation enables them to settle in modulated zigzag structures similar to those depicted above [17]. A further variation, which appears to be promising, uses bubbles in a liquid-filled tube.

While we have so far investigated only simple structures with single peaks in the displacement profile, more complicated structures exist at higher energies, and may also be found with the stepwise method. Among these are structures that can be created by concatenating one of the single-peak structures with its mirrored counterpart. Their compression and energy will be doubled.

Other extensions will include the case of soft (elastic) spheres, for which we have already observed similar effects, using hydrogel particles and bubbles. Also the observations can be extended to much higher compression, making contact with the work of Lee et al. [13] and Winkelmann et al. [14], for the 3d structures generated. It may also be possible to take advantage of a technique that uses photoelastic material to indicate the magnitude of the compressive forces [18, 19].

Considerable current interest in the compressed sphere chain focusses on motion of kinks and the corresponding Peierls–Nabarro potential. This may be estimated by making a smooth interpolation of the energy values for stable and unstable states as calculated here.

We hope that the results presented here will find direct comparison with previous work, particularly with regards to ions confined in traps [2, 3, 12, 4, 5, 1, 6, 7, 8] as well as assemblies of magnetic particles in a channel [20]. This work should also be relevant to other systems in which buckling is a key feature: for example localised buckling has recently been observed in experiments involving an expanding (growing) elastic beam pinned to a substrate [21].

Acknowledgements.

We wish to thank S. Burke for help with the experimental set-up. This work was supported by: EPSRC grant numbers EP/K032208/1 and EP/R014604/1, as well as an Irish Research Council Postgraduate Scholarship (project ID GOIPG/2015/1998). We also acknowledge the support of Science Foundation Ireland (SFI) under grant number 13/IA/1926.