Charged elastic rings: deformation and dynamics

The mechanical relaxation of the electrically charged elastic ring system is dominated by the competing electrostatic and elastic forces. Under the long-range electrostatic force, the point charges tend to be uniformly distributed over an expanding sphere in three-dimensional space, which is suppressed by the elastic force. To investigate the ground state shape of the charged elastic ring, we first analyze the interaction between elementary charge pairs along the ring.

A charge pair consists of two point charges of the same sign separated by $a$. For two charge pairs at distance $R$ in the x-y plane, as shown in the upper panel in figure 1(b), the total electrostatic energy is

$\displaystyle U(\theta_{1},\theta_{2})=\frac{2\beta}{R}\left[2+\tilde{a}^{2}(4-3(\sin^{2}\theta_{1}+\sin^{2}\theta_{2}))\right],$

(3)

where $\tilde{a}=a/R$. From equation (3), we see that the two charge pairs tend to be parallel to each other and perpendicular to the connecting line for minimizing the interaction energy. For the three-dimensional case, where one charge pair is fixed and the other charge pair can freely rotate in the x-z plane [see the lower panel in figure 1(b)], the total electrostatic energy is

$\displaystyle U(\phi)=\frac{2\beta}{R}\left[2(1-\tilde{a}^{2})+3\tilde{a}^{4}(1+\cos^{2}\phi)\right].$

(4)

Equation (4) shows that the charge pair at $y=R$ tends to be perpendicular to the other charge pair by taking $\phi=\pi/2$ for minimizing the interaction energy. The analysis of equations (3) and (4) suggests that a circular charged ring might be unstable.

Now, we discuss the electrostatic force on a closed charged ring in the continuum limit. Variational calculation shows that the electrostatic force is always along the principal normal direction goldstein1995nonlinear ; santiago2013elastic . Specifically, the energy functional associated with the electrostatics of a closed, uniformly charged space curve $\vec{r}(s)$ is

$\displaystyle\mathcal{E}=\frac{1}{2}\oint ds\oint ds^{\prime}\Phi(|\vec{r}(s)-\vec{r}(s^{\prime})|).$

(5)

The corresponding electrostatic force

$\displaystyle F_{e}=-(\kappa(s)+\hat{n}\cdot\nabla)\oint ds^{\prime}\Phi\left(r(s,s^{\prime})\right),$

(6)

where $\kappa(s)$ is the curvature, $\hat{n}$ is the principal normal vector, and $r(s,s^{\prime})=|\vec{r}(s)-\vec{r}(s^{\prime})|$ goldstein1995nonlinear ; santiago2013elastic . The electrostatic force acting on the length element at $\vec{r}(s)$ is along the normal direction $\hat{n}(s)$. Note that equation (6) is valid for any pair interaction potential $\Phi(r)$. While the purely normal electrostatic force could support the planar circular shape with constant curvature, it is unknown if the circular ring is stable, especially by the interplay of electrostatics and elasticity.

To explore the lowest-energy configuration of charged elastic rings under the competing electrostatic and elastic forces in three-dimensional space, we resort to numerical simulations at high precision based on the steepest descent method. Numerical simulation has proven to be a suitable tool to investigate the mechanically equilibrium properties and dynamics of long-range interacting particle systems yao2016electrostatics ; yao2021epl .

Preliminary numerical simulations show that the deformed rings in three-dimensional space are ultimately flattened to circles under the combined electrostatic and elastic forces. However, closer examination shows the existence of slight, but persistent sinusoidal deformations on the rings. Notably, in the long-time relaxation process, wave structure is developed on the originally circular ring. The wave number keeps reducing as the system evolves towards the lowest-energy shape. Simulation details are provided in SI. Note that the wave structure on the ring may be stabilized by introducing the attribute of spontaneous curvature to the ring system hossein2020spontaneous . The undulated shapes also provide clues for the ring geometry due to both thermal fluctuations and combination of the electrostatic and elastic interactions. In figure 2(a), we plot the distance $r_{i}$ from each particle to the center of mass for the numerically obtained lowest-energy shape under the precision of step size $s=10^{-8}$. The $r_{i}$ curve could be precisely fitted by the shape of the $n=2$ mode: $r(\theta)=r_{1}+h\cos(2\theta+\theta_{0})$. Both shapes of the $n=1$ and $n=2$ modes are shown in figure 2(b); the reference circles (dashed, red) are also shown. The relative deviation of the fitting function and the data, which is measured by $\sum_{i=1}^{N}|r_{i}-r(\theta_{i})|/\langle r_{i}\rangle$, is as small as $10^{-10}$.

The shape of the $n=1$ mode is absent in numerical simulations, because this mode essentially describes the translation of the entire ring in the perturbation regime. For the $n=1$ mode, the translation $\delta x$ along the x-axis satisfies the relation $(r_{1}+h_{0}\cos\theta-r_{0})/\delta x=\cos\theta$, which leads to $\delta x=h_{0}+O(h_{0}^{2})$. It is only for the $n=1$ mode that the linear term of $h_{0}$ in the expression for $\delta x$ is independent of $\theta$. As such, the $n=2$ mode represents the first deformation mode in the perturbation regime that breaks the continuous rotational symmetry of the circular ring.

We further quantify the degree of the out-of-plane deformation by the mean squared distance $\Delta$ of the particles with respect to the reference plane as spanned by the three particles $i_{1}=1$, $i_{2}={\textrm{floor}}(N/3)$, and $i_{3}={\textrm{floor}}(2N/3)$.

$\displaystyle\Delta=\frac{1}{\ell_{ext}}\sqrt{\frac{\sum_{i=1}^{N}d_{i}^{2}}{N}}.$

(7)

where $d_{i}$ is the distance from the particle $i$ to the reference plane. $\ell_{ext}$ is the lateral extension of the deformed ring, which is equal to the mean side length of the reference triangle. Simulations show that the magnitude of the out-of-plane deformation is much smaller than that of the sinusoidal deformation by two orders of magnitude; the detailed information is provided in SI. In other words, the sinusoidal deformation occurs approximately in the plane. It is further observed that the standard deviation of the bond length distribution is much smaller than the amplitude of the sinusoidal deformation by two orders of magnitude.

Systematic numerical simulations show that the persistence of the quasi-planar sinusoidal deformation in the lowest-energy configurations is a common phenomenon for varying $N$ and $\beta$. We examine typical cases in the parameter space of $N\in[100,4000]$ and $\beta\in[10^{-5},10^{6}]$ by high-precision simulations, and observe that the lowest-energy configurations of the charged elastic rings uniformly converge to the shape of the $n=2$ mode; the maximum residual force on the particles is typically at the order of $10^{-7}f_{0}$. We shall emphasize that the amplitude of the sinusoidal deformation is much smaller than the mean radius of the ring. For example, analysis of the lowest-energy configurations as obtained for the step size $s$ at the order of $10^{-4}$ shows that the ratio of the amplitude of the sinusoidal deformation and the mean radius is at the order of $10^{-5}$ when the value of $\beta$ is changed by three orders of magnitude from unity to 1000. As such, it is a challenge to experimentally observe the sinusoidal deformation in the lowest-energy charged elastic ring. The exclusion of thermal fluctuations and acoustic waves in numerical simulations allows us to reveal the sinusoidal deformation.

In the electrostatics-driven deformation of the charged elastic ring, the electrostatic and elastic energies have opposite dependence on the system size. As such, the ring size ultimately reaches an equilibrium value in the relaxation process. The variations of energy and ring size in the relaxation process are plotted in figure 3(a) and 3(b), respectively. In the following, we perform scaling analysis of the equilibrium size of the ring. The asymptotic expression for the electrostatic energy of $N$ point charges evenly distributed on a circle of radius $R$ in the large $N$ limit is Saff_asymptoticsfor ; Borodachov2012

$\displaystyle E_{{\textrm{e}}}\propto\frac{1}{R}N^{2}\ln N.$

(8)

The elastic stretching energy

$\displaystyle E_{{\textrm{s}}}\propto\frac{1}{N}(R-R_{0})^{2},$

(9)

where $R_{0}=N\ell_{0}$. $\ell_{0}$ is the balance length of the spring. From $\delta E_{{\textrm{tot}}}/\delta R|_{{\textrm{equil.}}}=0$, we obtain the scaling law:

$\displaystyle R_{{\textrm{equil.}}}\propto N(\ln N)^{1/3}.$

(10)

Equation (10) is confirmed in our numerical simulations (see SI).

Here, it is of interest to discuss the division of a charged ring into two identical rings based on equation (8); the perimeter of the ring is conserved in the division process. Without considering the interaction of the two rings, the change of the electrostatic energy is

$\displaystyle\frac{\Delta E}{E_{0}}=-\frac{\ln 2}{\ln N}<0,$

(11)

which indicates the reduction of the electrostatic energy as the ring is equally divided. This result is the two-dimensional version of the electrostatics-driven instability of charged spheres in the seminal work by Rayleigh rayleigh1882xx .

In this subsection, we perform analytical perturbation calculation to further substantiate the main conclusion of the preceding numerical analysis that the ground state shape of the charged elastic ring is not a perfect circle as one may expect intuitively.

According to the numerical results, both the out-of-plane deformation and the standard deviation of the bond length distribution in the lowest-energy configurations are vanishingly small. As such, we consider the configurations of uniformly distributed particles on a perfect circle and on a slightly deformed ring in the plane. Here, for the sake of analytical tractability as well as considering the numerical results, the bond length on both configurations is set to be identical, so their elastic energies are the same. The variation of the energy due to the deformation of the ring is completely determined by the electrostatic part. Could the electrostatic energy be lowered by deforming the uniformly charged circular ring? This question is also closely related to the stability of a uniformly charged circular ring.

In the following, we perform perturbation analysis in the continuum limit to address this question. According to the preceding discussion, the linear charge density is invariant in the deformation, and the deformation is confined in the plane. Now, consider a slight deformation of the uniformly charged ring in the plane:

$\displaystyle\vec{r}(\theta)=\vec{r}_{0}+h(\theta)\hat{n}(\theta),$

(12)

where $\hat{n}(\theta)$ is the outward unit normal vector on the original circular ring of radius $r_{0}$. $\vec{r}_{0}=r_{0}\hat{n}(\theta)$. The resulting change of the total electrostatic energy is

$\displaystyle\Delta E=\lambda^{2}\oint d\ell\oint d\ell^{\prime}V(r_{pp^{\prime}})-\lambda^{2}\oint d\ell_{0}\oint d\ell^{\prime}_{0}V(r_{p_{0}p_{0}^{\prime}}),$

(13)

where the first and the second terms are the integrals over the perturbed and circularly distributed charges. $\lambda$ is the linear charge density. $V(r)$ is the interaction potential of two unit charges separated by distance $r$. $d\ell$ and $d\ell^{\prime}$ are the line elements at any two points $p$ and $p^{\prime}$ on the perturbed ring. And $d\ell_{0}$ and $d\ell^{\prime}_{0}$ are the line elements at any two points $p_{0}$ and $p^{\prime}_{0}$ on the circular ring.

In the following, we expand the expression of $\Delta E$ in terms of the perturbation $h(\theta)$. We first expand $|\vec{r}_{pp^{\prime}}|$, which appears in the first term in equation (13), in terms of $h$. The position vector connecting the points $p$ and $p^{\prime}$:

$\displaystyle\vec{r}_{pp^{\prime}}=\vec{r}_{p_{0}p_{0}^{\prime}}+\delta\vec{r},$

(14)

where $\delta\vec{r}=h(\theta^{\prime})\hat{n}(\theta^{\prime})-h(\theta)\hat{n}(\theta)$. By using the following relation:

$\displaystyle\left|\vec{f}+\delta\vec{f}\,\right|=f+\Delta f,$

(15)

where $f=|\vec{f}\,|$ and

$\displaystyle\Delta f=\frac{\vec{f}\cdot\delta\vec{f}}{f}+\frac{1}{2}\frac{\delta\vec{f}\cdot\delta\vec{f}}{f}-\frac{1}{2}\frac{(\vec{f}\cdot\delta\vec{f})^{2}}{f^{3}}+O(\delta f^{3}),$

we have

$\displaystyle r_{pp^{\prime}}=r_{p_{0}p_{0}^{\prime}}+\Delta r,$

(17)

where

$\displaystyle r_{p_{0}p_{0}^{\prime}}=2r_{0}\biggr{|}\sin\frac{\theta-\theta^{\prime}}{2}\biggr{|},$

(18)

and

$\displaystyle\Delta r=\frac{r_{p_{0}p_{0}^{\prime}}}{2r_{0}}\left(h(\theta)+h(\theta^{\prime})\right)+O(h^{2}).$

(19)

Now, we insert the following equations to equation (13):

$\displaystyle d\ell=r_{0}\left(1+\frac{h(\theta)}{r_{0}}+\frac{1}{2}\frac{h^{\prime}(\theta)^{2}}{r_{0}^{2}}+O(h^{3})\right)d\theta,$

(20)

and

	$\displaystyle V(r_{pp^{\prime}})=V(r_{p_{0}p_{0}^{\prime}})$	$\displaystyle+$	$\displaystyle\left.\frac{dV(x)}{dx}\right|_{x=r_{p_{0}p_{0}^{\prime}}}\Delta r$
		$\displaystyle+$	$\displaystyle\frac{1}{2}\left.\frac{d^{2}V(x)}{dx^{2}}\right|_{{x=r_{p_{0}p_{0}^{\prime}}}}\Delta r^{2},$		(21)

where the expression for $\Delta r$ is given in equation (19). We finally obtain

	$\displaystyle\Delta E=\frac{\lambda^{2}}{r_{0}}\oint d\ell_{0}\oint d\ell_{0}^{\prime}\Big{[}(h(\theta)+h(\theta^{\prime}))\Big{(}V(r_{p_{0}p_{0}^{\prime}})$
	$\displaystyle+\frac{r_{p_{0}p_{0}^{\prime}}}{2}\frac{dV(r_{p_{0}p_{0}^{\prime}})}{dr_{p_{0}p_{0}^{\prime}}}\Big{)}+O(h^{2})\Big{]}.$		(22)

Now, we discuss the general case of $V(r)=\beta/r^{\gamma}$ ($\gamma$ is a positive integer). While we focus on the Coulomb potential ($\gamma=1$) in this work, scrutiny of the general case of $V(r)$ reveals the speciality of the long-ranged Coulomb potential as will be shown. Up to the linear term, we have

$\displaystyle\Delta E=\frac{(2-\gamma)\beta\lambda^{2}}{2r_{0}}\oint d\ell_{0}\oint d\ell_{0}^{\prime}\frac{h(\theta)+h(\theta^{\prime})}{r_{p_{0}p_{0}^{\prime}}^{\gamma}}.$

(23)

Now, consider the perturbation in the form of

$\displaystyle r(\theta)=r_{0}-\delta r_{0}+h_{0}\cos n\theta,$

(24)

where the deformation mode is characterized by the nonzero integer $n$. The shapes of the $n=1$ and $n=2$ modes are presented in figure 2(b); the red dashed reference circles are also shown. The $\delta r_{0}$ term in equation (24) (the zero mode) is introduced for the conservation of the contour length of the ring, which leads to the relation between $\delta r_{0}$, $h_{0}$ and $n$: $n^{2}h_{0}^{2}=4\delta r_{0}(r_{0}-\delta r_{0})$. By inserting equation (24) into equation (23), we have

$\displaystyle\frac{\Delta E}{2-\gamma}$

$\displaystyle\propto$

$\displaystyle\int_{\Omega}d\omega\frac{-\delta r_{0}+\frac{1}{4}h_{0}\cos(\frac{1}{2}nx)\cos(\frac{1}{2}ny)}{|\sin\frac{x}{2}|^{\gamma}},$

where

$\displaystyle\int_{\Omega}d\omega\equiv\int_{0}^{2\pi}dy\int_{-y}^{y}dx+\int_{2\pi}^{4\pi}dy\int_{-(4\pi-y)}^{4\pi-y}dx.$

The second term associated with $h_{0}$ in the integrand in equation (LABEL:DE_n) is zero for nonzero integers $n$ and $\gamma$. This conclusion is also true by replacing $h_{0}\cos n\theta$ for $h_{0}\sum_{m}\cos m\theta$ in equation (24), where $m$ is a nonzero integer. It is important to notice that the sign of the remaining term in equation (LABEL:DE_n) depends on the value of $\gamma$. The marginal value for $\gamma$ is 2 for the ring system, which is in contrast to the value of 4 for the charged sphere system jadhao2015coulomb . This observation suggests the connection of the instability of long-range interacting systems and the dimension of space. For $\gamma>2$, $\Delta E>0$, indicating the stability of the circular ring under relatively short-ranged interactions. However, $\Delta E<0$ at $\gamma=1$. This result indicates that the ground state shape of the charged ring under Coulomb potential is not a perfect circle, which is consistent with numerical observation.

We further perform numerical simulations to check the above analytical perturbation results. Specifically, we place a collection of point charges evenly along the plane curve of the $n$-mode shape specified by equation (24), and compute the total electrostatic energy. By comparison with the case of the corresponding circular ring of the same contour length and charge density, we obtain the change of the total electrostatic energy $\Delta E$ in the deformation. In figure 4, we present the plots of $\Delta E$ versus $\gamma$ for the $n=1$ and $n=2$ modes. Figure 4 shows that $\Delta E<0$ for $\gamma=1$, and $\Delta E>0$ for $\gamma>1$. These numerical results are consistent with the analytical perturbation analysis. The small deviation of $\Delta E$ from zero for $\gamma=2$ and $n=2$ in figure 4(b) shall be attributed to the discrepancy of the continuous and the discrete charge distributions adopted in the perturbation calculation and numerical simulations, respectively. Examination of higher order terms ($n=3,4,5,6$) also shows that $\Delta E<0$ for $\gamma=1$, and $\Delta E>0$ for $\gamma>1$.

To conclude, the analytical perturbation calculation reveals the instability of the circular shape of uniformly charged rings in the continuum limit. The combination of the analytical perturbation analysis and the numerical results shows that the ground state shape of the charged elastic ring is indeed not a perfect circle. While the $n=2$ mode has been identified as the lowest-energy shape by the high-precision steepest descent method, we shall note that it is still an open question to rigorously prove if this elliptic shape is the ground state of the system. The subtlety of the long-range nature of the interaction imposes challenges and brings richness to this question.

The crucial role of the long-range electrostatic force in the static deformation of the charged elastic ring system has been discussed in the preceding subsections. In this subsection, we proceed to discuss the dynamical effect of the electrostatic force. To this end, we impose random disturbance by specifying a random velocity field along the ring in the lowest-energy state, and analyze the temporally-varying kinetic energy curve. The dynamical evolution of the ring upon the disturbance is governed by the Hamiltonian dynamics. It is found that the dynamical effect of the electrostatic force is reflected in the modulation of the frequency. The relevant results are summarized in figure 5.

Figure 5(a) shows the temporally-varying kinetic energy $E_{K}(t)$. The total energy of the system is well conserved in numerical simulations; the fluctuation of the total energy is at the order of $10^{-10}\%$. We perform fourier transformation of the irregular kinetic energy curve, and plot the $\tilde{E}_{K}(\omega)$ curve in figure 5(b). Figure 5(b) shows that the dynamical rhythm of the interacting point charges connected by linear springs is dominated by a pair of frequencies $\omega_{1}$ and $\omega_{2}$. For comparison, we present the value for the eigenfrequency $\omega_{0}$ of the elementary freestanding system consisting of two identical masses $m$ connected by a linear spring of stiffness $k_{0}$. $\omega_{0}=\sqrt{2k_{0}/m}$, and the corresponding frequency of the kinetic energy curve is doubled: $\omega_{E_{k}}=2\sqrt{2}\approx 2.8$.

To explore the physical origin of the dominant frequencies, we examine several global and local dynamical modes, compute the relevant frequencies, and finally identify the mode possessing the identical spectral structure as in the $\tilde{E}_{K}(\omega)$ curve. This mode corresponds to the relative motion of nearest neighboring particles: $\delta\vec{x}_{i}=(\vec{x}_{i-1}-\vec{x}_{i})+(\vec{x}_{i+1}-\vec{x}_{i})$ [see the inset in figure 5(c)]. By averaging over all of the particles and applying fourier transformation, we plot the $\delta\tilde{x}(\omega)$ curve in figure 5(c). It turns out that the frequencies at the two peaks in the $\delta\tilde{x}(\omega)$ curve in figure 5(c) are identical to those in the $\tilde{E}_{K}(\omega)$ curve in figure 5(b). In general, oscillation of particles leads to variation of energy. Here, we identify the long-range electrostatic force driven short-wavelength dynamical mode (as characterized by $\delta\tilde{x}(\omega)$) that could explain the origin of the pair of dominant frequencies in the energy curve in figure 5(b).

We further find that the dominant frequencies, especially the value of $\omega_{1}$, could be modulated by the strength $\beta$ of the electrostatic force. Figure 5(d) shows the variation of the pair of the dominant frequencies with the increase of $\beta$ for the cases of $N=500$ (filled marks) and $N=1000$ (empty marks). The insensitivity of the frequency curves on the total number of particles could be attributed to the short-wavelength nature of the dominant oscillation as presented in figure 5(c). Here, we shall also emphasize the existence of the pair of the dominant frequencies in a wide range of $\beta$ covering several orders of magnitude.

We finally briefly discuss the model of the charged elastic ring. The simplicity of this model allows us to perform analytical perturbation calculation (under certain limiting conditions). This model exhibits the deformation under the purely repulsive Coulomb potential among all of the cases of $\gamma\in Z^{+}$ in the interaction potential $V(r)\sim 1/r^{\gamma}$, demonstrating the subtlety of the long-range force in the organization of matter. Considering the significance of electrostatic force in regulating materials at the nanoscale, the sinusoidal deformation phenomenon may occur in small-scale structures in electrolyte solutions Walker2011 . The smallness of the intrinsic sinusoidal deformation imposes a challenge for direct observation in real systems that are generally subject to thermal fluctuations. In this aspect, the idealized charged elastic ring model shows its value in revealing the subtle effect of long-range interaction. This effect, as demonstrated in the simple ring system, represents a kind of complexity that is distinct from what is caused by many degrees of freedom.

The ring model also serves as an example of shape instability under stretching. In contrast to the well studied compression driven instabilities of various elastic systems, the stretching caused instability under long-range repulsion has received much less attention timoshenko1951theory ; Landau1986 ; audoly2010elasticity . Note that the electromechanical stiffening effect has been reported in both elastic membrane ANDELMAN1995603 and rod PhysRevE.67.061805 ; Rudnick2004 systems.

Here, we may mention that the exploration of the instability phenomenon of the simple ring model (the $S^{1}$ geometry) is inspired by our previously discovered instability of the charged sphere system (the $S^{2}$ geometry) jadhao2015coulomb . To the best of our knowledge, the main results of the charged elastic ring system (despite of its simplicity) presented in this work have not been reported in previous work. Literature search shows that, since the seminal work by Maxwell on the equilibrium distribution of electric charges on a long narrow cylinder maxwell1877electrical , the one-dimensional electrostatics problem has been extended from straight lines griffiths1996charge ; jackson2000charge ; jackson2002charge ; amore2019thomson to circular Ball1977 ; Zypman2006 ; Hafeez2009 ; Villa2019 and eccentric Denisova1998 ; amore2019thomson rings of frozen geometry in both continuum and discrete regimes. The electrostatic force on a closed, uniformly charged space curve has been derived based on variational principle goldstein1995nonlinear ; santiago2013elastic . Statics and dynamics of flexible ring systems have been studied for the uncharged case PhysRevE.64.011909 .