text:
# Criticality in confined ionic fluids

## I Introduction

Over the last decade the need to understand coulombic criticality has provided a new impetus to the study of electrolyte solutions. The current wave of exploration can be traced to the pioneering experiments of K. S. Pitzer et. al., who have first reported a surprising finding that Coulomb interactions might belong the mean-field universality class . This suggestion has not gone unchallenged and, in fact, later experiments are consistent with a crossover from mean-field to Ising universality class very close to the critical point . The crossover, if indeed it exists, is much closer to the critical point than for any other known fluid.

From the theoretical point of view it is very hard to justify anything but Ising criticality . The goal for the theorists must then lie in a seemingly simpler task of finding why the crossover region for electrolytes is so narrow. Unfortunately even this programme has failed to produce any satisfactory explanation . Most calculations suggest that the critical region for electrolytes should be comparable to that of a Lennard-Jones fluid. To further confound the mystery, Monte Carlo simulations are once again pointing in the direction of mean-field criticality .

A seemingly unrelated problem concerns the disappearance of the anticipated liquid-gas transition in a system of dipolar hard spheres (DHS). Since the DHS is the simplest realization of a polar fluid, for a long time it has been believed that it must exhibit a liquid-gas phase separation. It came, therefore, as quite a surprise when the Monte Carlo simulations failed to located this singularity . Instead what they found was that as the temperature was lowered, the dipolar particles aligned, forming polymer-like chains. Since these chains interact weakly , it has been argued that the liquid-gas transition must be driven entirely by the free unassociated dipoles . In fact a critical density of free dipoles is needed for the phase transition to occur. At low temperatures, where most of the theories localize the transition, formation of dipolar chains depletes the density of free dipoles bellow the critical threshold necessary for the transition to take place .

In this paper we shall study a third model which we hope might span the gap between the three dimensional Coulomb gas and the dipolar hard spheres, and thus shed some additional light on the criticality in these interesting and important systems. Our new model consists of a neutral electrolyte confined to a two dimensional plane . This can be visualized as oppositely charged surfactant molecules adsorbed to a water-oil interface. An example of such a system is cetyltrimethyl amonium-hydroxy naphthalene carboxylate $`(CTAHNC)`$, which is composed of two surface active parts, $`CTA^+`$ and $`HNC^{}`$ .

We shall argue that unlike the $`3d`$ electrolyte , the confined $`2d`$ plasma might not phase separate. Instead as the temperature is lowered, chains composed of alternating positive and negative ions will begin to form $`(\mathrm{}++\mathrm{})`$ . Just as dipolar chains, these new clusters interact weakly between themselves. However, they diminish the concentration of free ions bellow the critical value necessary for the transition to take place.

## II The Model

Our system consists of an ionic fluid of $`N_+`$ positive and $`N_{}`$ negative ions confined to a plane of area $`A`$ located at $`z=0`$, separating two different solvents occupying the half-spaces at $`z>0`$ and $`z<0`$. We shall restrict our attention to the neutral electrolytes, for which $`N_+=N_{}=N/2`$. The solvents are treated as uniform mediums with dielectric constants $`ϵ_2`$ and $`ϵ_1`$ for $`z>0`$ and $`z<0`$, respectively. The ions of both species are modeled as hard spheres of diameter $`a`$ and charge $`\pm q`$ located at the center. It will be convenient to define the dimensionless particle density as $`\rho _\pm ^{}=\rho _\pm a^2`$, where $`\rho _i=N_i/A`$. The total density of ions is $`\rho =\rho _++\rho _{}`$.

All the relevant thermodynamic information is contained in the free energy density, $`f=F/A`$. Unfortunately due to the complexity of interactions, no exact expression for $`f`$ can be found. We shall, therefore, attempt to construct the approximate free energy using the most relevant contributions. These can be divided into entropic and electrostatic,

$$f(T,\rho _+,\rho _{})=f_{ent}^{(1)}(T,\rho _+,\rho _{})+f_{el}(T,\rho _+,\rho _{}).$$
(1)
The entropic (mixing) free energy is given by

$`\beta f_{ent}^{(1)}(T,\rho _+,\rho _{})`$ $`=`$ $`\left[\rho _+\mathrm{ln}(\rho _+\mathrm{\Lambda }^2)\rho _++\rho _{}\mathrm{ln}(\rho _{}\mathrm{\Lambda }^2)\rho _{}\right]`$ (2)
where $`\mathrm{\Lambda }=\sqrt{2\pi \mathrm{}^2/mk_BT}`$ is the thermal de Broglie wavelength of the ions. The second term in Eq. (1) is due to electrostatic interactions. It is important to note that the electrostatic free energy is purely correlational, since the mean-field contribution is zero. To calculate $`f_{el}`$, let us fix one ion at the origin. Adopting the cylindrical coordinate system $`(\varrho ,\phi ,z)`$, the central ion is located at $`\varrho =0,z=0`$. Due to the electrostatic interactions the other particles will arrange themselves within the plane in accordance with the Boltzmann distribution. Since no charge is present in the regions $`z>0`$ and $`z<0`$, the electrostatic potential there must satisfy the Laplace equation. Appealing to the azimuthal symmetry and taking into account the fact that the potential should vanish at infinity we find ,

$$\varphi _1(\varrho ,z)=_0^{\mathrm{}}A_1(k)J_0(k\varrho )\text{e}^{kz}𝑑kforz<0$$
(3)
and

$$\varphi _2(\varrho ,z)=_0^{\mathrm{}}A_2(k)J_0(k\varrho )\text{e}^{kz}𝑑kforz>0,$$
(4)
where $`J_0(x)`$ is the Bessel function of order zero.

The functions $`A_1(k)`$ and $`A_2(k)`$ are determined by the boundary conditions: continuity of electrostatic potential,

$$\varphi _2(\varrho ,0)=\varphi _1(\varrho ,0),$$
(5)
and discontinuity of displacement field across the $`z=0`$ plane,

$$\left(ϵ_2𝐄_2(\varrho ,z)ϵ_1𝐄_1(\varrho ,z)\right)\widehat{n}=4\pi \sigma _{eff}(\varrho ),$$
(6)
where $`\sigma _{eff}(\varrho )`$ is the surface charge density and $`\widehat{n}`$ is a unit vector normal to the interface, pointing from region 1 to 2. The continuity of electrostatic potential results in $`A_1(k)=A_2(k)=A(k)`$, while Eq. (6) requires that

$$2_0^{\mathrm{}}kA(k)J_0(k\varrho )𝑑k=\frac{4\pi \sigma _{eff}(\varrho )}{D},$$
(7)
where $`D=(ϵ_1+ϵ_2)/2`$. The surface charge density  is given by

$$\sigma _{eff}(\varrho )=\sigma _s(\varrho )+\frac{q\delta (\varrho )}{2\pi \varrho }.$$
(8)
The term $`\sigma _s`$ is the charge density of the ”ionic cloud” around the central ion,

$$\sigma _s(\varrho )=q\rho _+e^{\beta q\varphi }q\rho _{}e^{+\beta q\varphi }.$$
(9)
Eqs. (7) and (8) are only valid in the limit $`a0`$. Unfortunately it is highly nontrivial to take into account the boundary condition associate with the finite ionic size. To circumvent this difficulty, we shall first solve all the equations in the point particle limit. Then to account for the finite particle size, we shall replace the bare charge of the central ion $`q`$, in Eq. (8), by an effective charge $`Q`$, $`qQ`$. The effective charge $`Q`$ will be determined by the condition of an overall charge neutrality,

$$2\pi _a^{\mathrm{}}\sigma _s(\varrho )\varrho 𝑑\varrho =q.$$
(10)
For the $`3d`$ electrolyte this procedure leads to an electrostatic potential identical with the one found from the exact solution of the Debye-Hückel equations with the appropriate hard core boundary conditions . In the present geometry the procedure outlined above will only be an approximation, although we believe a rather good one.

In the spirit of Debye-Hückel theory  we shall now linearize the Boltzmann factor in Eq. (9). The surface charge density then reduces to

$$\sigma _s(\varrho )=\frac{D\kappa _s\varphi (\varrho )}{2\pi },$$
(11)
where $`\kappa _s=2\pi \rho ^{}/T^{}a`$ is the inverse Gouy-Chapman length and $`T^{}=k_BTDa/q^2`$ is the reduced temperature. Eq. (7) can now be solved yielding the expression for electrostatic potential,

$$\varphi (\varrho ,z)=\frac{Q}{D}_0^{\mathrm{}}\frac{k}{k+\kappa _s}J_0(k\varrho )\text{e}^{k|z|}𝑑k.$$
(12)
For $`z=0`$, the in plane potential agrees with the one obtained by Velazquez and Blum  and can be conveniently rewritten as

$$\varphi _>(\varrho )=\frac{Q\tau _0(\kappa _s\varrho )}{D\varrho },$$
(13)
The subscript ”$`>`$” is included to stress that for rigid particles this form will be appropriate only outside the hard core exclusion region, $`\varrho >a`$. The function $`\tau _\nu (x)`$ is defined as

$$\tau _\nu (x)=1\frac{\pi x^{1\nu }}{2}[H_\nu (x)N_\nu (x)]$$
(14)
where $`H_\nu (x)`$ and $`N_\nu (x)`$, are the Struve and the Bessel functions of order $`\nu `$, respectively . The charge neutrality, Eq. (10), together with Eqs. (11) and (13), determines the effective charge $`Q`$,

$$Q=\frac{q}{\kappa _sa\tau _1(\kappa _sa)}.$$
(15)
In the limit of $`a0`$, the renormalized charge reduces to the bare charge $`q`$.

To obtain the electrostatic free energy, we require the potential inside the excluded volume region, $`\varphi _<(\varrho )`$. Unfortunately the procedure presented above leaves this undetermined. The only statement that we can make is that $`\varphi _<(\varrho )`$ must be of the form

$$\varphi _<(\varrho )=\frac{q}{D\varrho }+C(\varrho ).$$
(16)
For a three dimensional unconstrained electrolyte, $`C`$ is a constant. This, however, is not the case in the present geometry and $`C(\varrho )`$ is a function of position. In particular the value $`C(0)`$ is the potential that the central ion feels due to the presence of other particles. In their earlier study, Velasquez and Blum  approximated $`C(\varrho )`$ by a constant which they then determined by requiring continuity of electrostatic potential across the exclusion boundary $`\varphi _<(a)=\varphi _>(a)`$. This, however, is a very rough approximation since there is nothing to prevent $`C(\varrho )`$ from being a very strongly varying function of position. To avoid this difficulty we shall use an alternate method of obtaining the value of $`C(0)`$. To this end we note that the potential at the center of a circularly symmetric charge distribution is

$$C(0;\rho )=\frac{2\pi }{D}_a^{\mathrm{}}\sigma _s(\varrho )𝑑\varrho ,$$
(17)
where in order to emphasize that the potential depends on the density of ions, we have explicitly included $`\rho `$ in its definition. Using Eq. (11) with the potential given by Eq. (13) we find

$$C(0;\rho )=\frac{q}{Da\tau _1(k_sa)}_{\kappa _sa}^{\mathrm{}}\frac{\tau _0(z)}{z}𝑑z.$$
(18)
The excess chemical potential can now be calculated straight forwardly by appealing to the Güntelberg charging process . We find $`\mu _\pm ^{ex}=qC(0;\rho )/2`$. The chemical potential of the positive and negative ions is

$$\beta \mu _\pm =\frac{f}{\rho _\pm }=\mathrm{ln}(\frac{\rho \mathrm{\Lambda }^2}{2})+\frac{\beta qC(0;\rho )}{2}.$$
(19)
The critical point is determined by the conditions, $`\mu _\pm /\rho =0`$ and $`^2\mu _\pm /\rho ^2=0`$, which reduce to

$`2T^{}\tau _1^2+I\tau _1I\tau _0\tau _1\tau _0=0`$ (20)
(21)
$`I\tau _1^24I\tau _0\tau _13\tau _0\tau _1^2+I\tau _0^2+2\tau _0^2\tau _1+x^2I\tau _1^2+I\tau _1+x^2\tau _1^3+\tau _1^2=0`$ (22)
with

$`I(x)={\displaystyle _x^{\mathrm{}}}{\displaystyle \frac{\tau _0(z)}{z}}𝑑z,`$ (23)
where $`x=\kappa _sa`$. Solving Eqs. (20), the critical point is located at $`T_c^{}=0.0517386`$ and $`\rho _{c+}^{}=\rho _c^{}=0.00121988`$.

## III Debye-Hückel-Bjerrum approximation.

Clearly the low temperature at which the critical point is located should make us worry about the approximations which have been adopted. Certainly at such extreme conditions the linearization of the Boltzmann factor, Eq. (11), is no longer valid. Fortunately, all is not lost. Evidently, linearization of the Boltzmann factor in Eq. (9) diminishes the weight of configurations in which the oppositely charged ions are in a close proximity. It is possible, therefore, to correct for the omitted non-linearities by explicitly allowing for the formation of ”clusters”. These clusters are assumed to be in a chemical equilibrium with the free unassociated ions, and their density is determined by the law of mass action . The most basic such cluster is a dipole formed by a $`(+)`$ pair. Within the simplest version of this theory the dipoles are treated as non-interacting ideal specie. For $`3d`$ electrolyte this Debye-Hückel-Bjerrum approximation (DHBj)  has proven extremely successful, predicting the location of the critical point in close agreement with the Monte Carlo simulations .

The total density of ions $`\rho `$ can then be subdivided into that of free unassociated monopoles $`\rho _1=\rho _++\rho _{}`$, and of dipolar pairs $`\rho _2`$, with $`\rho =\rho _1+2\rho _2`$. In the spirit of DHBj theory we shall first treat the dipoles as ideal non-interacting entities whose concentration is determined by the law of mass-action, $`\mu _++\mu _{}=\mu _2`$. The Helmholtz free energy density is then given by

$$\beta f_{DHBj}=f_{ent}^{(2)}(\rho _1,\rho _2,T)+f_{el}(\rho _1,T)$$
(24)
where the entropic contribution is

$$\beta f_{ent}^{(2)}(\rho _1,\rho _2,T)=\left[\rho _2\mathrm{ln}\left(\frac{\rho _2\mathrm{\Lambda }^4}{\xi _2}\right)\rho _2\right]+\left[\rho _1\mathrm{ln}\left(\frac{\rho _1\mathrm{\Lambda }^2}{2}\right)\rho _1\right].$$
(25)
The $`\xi _2`$ is the internal partition function of a dipolar pair,

$$\xi _2(T;R)=2\pi _a^R\text{e}^{\beta U_2(\varrho )}\varrho 𝑑\varrho ,$$
(26)
where $`\beta U_2=a/T^{}\varrho `$ is the electrostatic potential between the associated ions. In order to evaluate $`\xi _2(T;R)`$ we must define the distance $`R`$ at which two ions can be considered to be associated. Following Bjerrum  we choose the value of $`R_{Bj}`$ at which $`\xi _2(T;R)`$, as a function of $`R`$, has an inflection point, $`R_{Bj}=a/T^{}`$ . With this, the integral in Eq. (26) can be evaluated to yield

$$\xi _2=\frac{\pi a^2\text{e}^b}{b}\left[b^3\text{e}^b\left(Ei(b)Ei(1)+2\text{e}\right)b(1+b)\right],$$
(27)
where $`b=1/T^{}`$ and $`Ei(x)`$ is the exponential integral function.

Since within the DHBj approximation the dipoles are treated as ideal, the electrostatic free energy only depends on the density of free monopoles, $`\rho _1`$. Substituting the free energy into the law of mass-action, the density of dipoles is given by

$$\rho _2^{}=\xi _2\rho _+\rho _{}e^{\beta (\mu _+^{ex}+\mu _{}^{ex})}.$$
(28)
The charge neutrality requires that $`\rho _+=\rho _{}=\rho _1/2`$. The expression for excess chemical potential was calculated in Section 2, $`\mu _\pm ^{ex}=\beta qC(0;\rho _1)/2`$, and the inverse Gouy-Chapman length is now $`\kappa _s=2\pi \rho _1^{}/T^{}a`$.

Within the DHBj approximation the dipoles are treated as ideal specie, therefore, they cannot influence the location of the critical point. Thus, the critical temperature must still be $`T_c^{}=0.0517386`$ while the critical density of monopoles must remain $`\rho _{1c}^{}=0.00243975`$. Substituting these values into Eq. (28), we find that the density of dipoles is $`\rho _{2c}^{}=1.08847`$, which is extremely high. If there are so many dipoles is it also not possible that there will be higher order clusters as well?

## IV Linear Ionic Chains

Unfortunately as soon as we get to clusters of three ions the calculations get extremely difficult. The basic problem is the internal partition function of the higher order clusters, which can no longer be evaluated exactly. Furthermore, while it is evident that for a cluster of three ions the low temperature configurations are chain-like, $`(++)`$ or $`(+)`$, this is far from obvious for a neutral cluster of two positive and two negative ions. The entropy favors a chain-like configuration $`(++)`$, while the energy favors a compact square cluster. Which will gain in the critical region is hard to say. We note, however, that exactly the same situation was encountered for dipolar hard spheres . In that case, in the vicinity of the critical point, the chain like configurations dominated. Since it is much easier to study the linear clusters they, therefore, will provide a starting point for our analysis.

We begin by supposing that at low temperatures our system will be composed of monopoles of density $`\rho _1`$ and chains of $`n`$ monomers with densities $`\rho _n`$. Once again in the spirit of Bjerrum we shall first treat the clusters as ideal non-interacting species. The particle conservation requires that

$$\rho =\underset{n=1}{\overset{\mathrm{}}{}}n\rho _n.$$
(29)
Consider an alternating chain composed of $`t`$ positive and $`s=nt`$ negative ions. The partition function for such a cluster is,

$$\xi _n(T)=\frac{1}{s!t!}\frac{1}{A}_{\mathrm{\Omega }_n}\underset{i=1}{\overset{n}{}}d^2𝐫_i\text{e}^{\beta U_n}.$$
(30)
where $`\mathrm{\Omega }_n`$ is the configurational volume and $`U_n`$ is the total energy of interaction between the ions forming a chain,

$$U_n(𝐫_1,𝐫_2,\mathrm{},𝐫_n)=\underset{i<j}{\overset{n}{}}\phi _{ij}(r_{ij}).$$
(31)
At zero temperature chains are rigid and the particles are in contact with one another. The displacement vector between two ions $`i`$ and $`j`$ is $`\stackrel{}{r}_{i,j}^0a|ij|\widehat{x}`$, where $`\widehat{x}`$ is the unit vector along the chain. The electrostatic energy of interaction between the ions of the chain can be evaluated exactly yielding

$$U_n^0=\frac{q^2}{Da}S(n),$$
(32)
where

$$S(n)=\underset{k=1}{\overset{n1}{}}\frac{(1)^k(nk)}{k}.$$
(33)
For nonzero but small temperatures, deviations from $`\stackrel{}{r}_{i,j}^0`$ occur. These fluctuations can be taken into account by making a Taylor expansion of $`U_n`$ around the ground state up to quadrupolar order

$$U_n=U_n^0\frac{q^2}{2D}\underset{ij}{\overset{n}{}}(1)^{i+j}\{\frac{(\stackrel{}{r}_{ij}\stackrel{}{r}_{ij}^0)_x}{(\stackrel{}{r}_{ij}^0)^2}+\frac{(\stackrel{}{r}_{ij}\stackrel{}{r}_{ij}^0)_y^2}{2|\stackrel{}{r}_{ij}^0|^3}\frac{(\stackrel{}{r}_{ij}\stackrel{}{r}_{ij}^0)_x^2}{|\stackrel{}{r}_{ij}^0|^3}\}.$$
(34)
Here the subscripts $`x`$ and $`y`$ represent the components along the chain’s direction and perpendicular to it, respectively. Since we are assuming small deviations from the ground state, for each ion $`i`$ we shall consider only fluctuation contributions to $`\xi _n`$ arising from the interactions between the nearest and the second nearest neighbors. Choosing as the unit bases $`\widehat{x}`$ and $`\widehat{y}`$, vectors parallel and perpendicular to the direction of the chain in the plane $`z=0`$, the nearest and the second nearest displacement vectors can be written as

$$\stackrel{}{R}_1^{(i)}=a(1+\lambda _i)\left(\begin{array}{cc}\hfill \mathrm{cos}\varphi _i& \\ \hfill \mathrm{sin}\varphi _i& \end{array}\right)$$
(35)
and

$$\stackrel{}{R}_2^{(i)}=a\left[\begin{array}{cc}\hfill (1+\lambda _i)\mathrm{cos}\varphi _i+(1+\lambda _{i+1})\mathrm{cos}\varphi _{i+1}& \\ \hfill (1+\lambda _i)\mathrm{sin}\varphi _i+(1+\lambda _{i+1})\mathrm{sin}\varphi _{i+1}& \end{array}\right],$$
(36)
respectively. Here $`\lambda `$ is the radial and $`\varphi `$ is the angular deviation from the relative positions in the ground state.

Substituting Eqs. (35) and (36) into Eq. (34), for small fluctuations the electrostatic energy becomes

$$\frac{U_n}{k_BT}=\frac{1}{T^{}}\left\{S(n)+\underset{k=1}{\overset{n1}{}}\lambda _k\frac{1}{4}\underset{k=1}{\overset{n2}{}}\left[\lambda _k+\lambda _{k+1}\frac{1}{4}(\varphi _k\varphi _{k+1})^2\right]\right\}.$$
(37)
In the low temperature limit, configurational integral, Eq. (30), can be performed explicitly yielding,

$$\xi _n(T^{})=\left(\frac{a^{2n2}}{9}\right)2^{2n}\pi ^{\frac{n}{2}}T_{}^{}{}_{}{}^{\frac{3n}{2}2}\text{e}^{S(n)/T^{}}forn3.$$
(38)
The prefactor is the result of thermal fluctuations while the exponential is due to the ground state energy.

The condition of chemical equilibrium between the monopoles and the n-chains is expressed through the law of mass action $`\mu _n=t\mu _++s\mu _{}`$. The chemical potential for monopoles is given by Eq. (19) with $`\rho \rho _1`$, while for n-chains the chemical potential is,

$$\mu _n=k_BT\mathrm{ln}\left[\frac{\rho _n\mathrm{\Lambda }^{2n}}{\xi _n(T)}\right]$$
(39)
The law of mass-action reduces to

$$\rho _n=\left(\frac{\rho _1}{2}\right)^n\xi _n(T^{})\text{e}^{n\beta \mu _\pm ^{ex}}$$
(40)
At the level of approximation that we have adopted, the chains are treated as non-interacting ideal species. This means that just like in the case of dipoles in Section 3, they cannot affect location of the critical point. Therefore $`T_c^{}=0.0517386`$ and $`\rho _{1c}^{}=0.00243975`$. Substituting these values into Eq. (40) we find that the sum in Eq. (29) diverges. In fact according to the Cauchy criterion, the sum in Eq. (29) will converge absolutely if and only if,

$$\mathrm{\Delta }\underset{n\mathrm{}}{lim}\rho _{n}^{}{}_{}{}^{1/n}<1.$$
(41)
Using Eq. (40),

$$\mathrm{\Delta }=2\rho _1^{}\sqrt{\pi T^3}\mathrm{exp}\left[\frac{\alpha +I(x)/\tau _1(x)}{2T^{}}\right]$$
(42)
where $`\alpha =2\mathrm{ln}2`$ and $`x=\kappa _s(\rho _1)a`$. Inserting the critical parameters into the expression above, we find that at criticality $`\mathrm{\Delta }_c=5.82706`$, and the Cauchy criterion is strongly violated. Therefore, the critical density of monopoles lies outside the radius of convergence of Eq. (29). This means that for any finite density $`\rho `$, the density of monopoles never reaches the threshold necessary for the phase separation to occur, $`\rho _1<\rho _{1c}`$.

## V Results and discussion

We have presented an argument which suggests that a confined $`2d`$ electrolyte should not phase separate. Instead, just as for the case of dipolar hard spheres, as the temperature is lowered, the ions will associate forming polymer-like chains of alternating positive and negative monomers. Clearly our argument is based on a number of assumptions. First, we have supposed that in the critical region the linear chains predominate over the compact clusters. This is far from obvious. If the compact clusters have lower free energy than the chains, they can provide nuclei for the condensation and the gas-liquid phase separation. A second assumption implicit in our calculations is that the chains interact only weakly. This is somewhat easier to justify. Consider an infinitely long rigid line of alternating charges $`(\mathrm{}+++\mathrm{})`$. Suppose that the center to center distance between the nearest neighbors is $`a`$. It is then possible to show that the potential produced by such a line of charge decays exponentially, $`\psi (r)\mathrm{exp}(\pi r/a)`$, where $`r`$ is the distance perpendicular to the chain. Thus, the interactions between long polymer-like clusters should, indeed, be quite weak. However, the shorter chains can still interact sufficiently strongly to drive phase transition . Finally, even if the formation of chains prevents the liquid-gas phase separation, it does not forestall other kinds of phase transitions from taking place . At sufficiently high densities, the two component plasma will crystallize exhibiting a pseudo long range order. At low temperatures and densities, where the polymer-like chains predominate, the Y-like defect formation  can lead to a coexistence between two phases, one with high and the other with low concentration of defects .

After this work was completed and submitted for publication, the referee drew our attention to a recently published simulation by Weis, Levesque, and Caillol (WLC)  of a $`2d`$ ionic fluid. Indeed, these authors found coexistence between high and low density phases. WLC estimated the critical temperature to be $`T_c0.04`$, which should be compared with our Debye-Hückel prediction of $`T_c=0.052`$. Furthermore, in the low density phase, WLC found predominance of chain and ring-like clusters. The high density phase resembles a percolating gel-like cluster . Although there is a phase coexistence, it is difficult to associate it with a traditional liquid-vapor transition. Instead the coexistence resembles more a sol-gel transition in polymer systems. The task for theorists must now be to quantitatively understand the phase transition found by WLC. We hope that the current paper will provide a first step in this direction.

Acknowledgments

This work was supported by the Brazilian science agencies CNPq, FINEP, and Fapergs. J.E.F.M. acknowledge the postdoctoral fellowship from the CNPq and CLAF during the development of this work, and partial financial support from CONACyT. M. C. B. is grateful for the hospitality at the Center of Polymer Studies at Boston University, where part of this work was done.


format: text

text:
# Contents

## 1 Introduction

Our object of study in this paper is a deflated rigid loop (elastica hypoarealis). To explain this, consider a closed loop in a plane. The loop is made of some elastic material and, if taken thin enough, the contribution of its longitudinal deformations to its elastic energy can be ignored. This gives rise to a hamiltonian density proportional to the square of its curvature — hence elastica. It is clear that fixing only the length $`L`$ of the loop will not lead to any surprises, the only possible equilibrium being a circle. What will happen though if we, in addition, fix its area $`A`$ to be less than that of the circle? (hence hypoarealis).

This particular hamiltonian (without the constraints) has a distinguished history, dating back at least to the Bernoullis. It makes a reappearance in Euler’s inspired analysis of a mundane but subtle problem in mechanical engineering: the buckling of a loaded beam . More recently, variants of this problem have attracted the attention of mathematicians. In particular, the problem of determining the curves of constant length which minimize the bending energy, on surfaces of constant curvature, is considered in Refs. , while interconnections with knot theory are explored in . Adding the constraint on the enclosed area, as we propose here, provides a particularly fruitful generalization of this work.

Hamiltonians depending on extrinsic curvature have also been studied in statistical physics . For space curves, an effective hamiltonian which depends on twist as well as bending provides a phenomenological description of stiff polymers, and in particular DNA (see e.g. Refs. ). Moving up one dimension, the leading term in the Helfrich hamiltonian, which describes the equilibrium configurations of lipid membranes, is proportional to the square of the traced extrinsic curvature $`K`$ integrated over the surface, namely, the conformally invariant Willmore functional . In the variant of the Helfrich model proposed by Svetina and Zeks, the so-called bilayer couple model, constraints are placed on the surface, such as constant area, constant enclosed volume, and constant integrated mean curvature (the latter breaking the $`KK`$ symmetry of the problem) .222 We note that there is no genuine one-dimensional analog of its refinement, the so-called area difference model . However, the equations which determine the equilibrium are highly non-trivial higher-order PDEs. One motivation for introducing our variant of the elastica, as a planar analogue of the above, is its potential as a toy model for understanding these higher dimensional membranes. At the simplest level, if a closed membrane possesses a symmetry along the $`z`$-axis, its profile at any $`z`$ will be a fixed loop. The dimensionally reduced Helfrich hamiltonian is then the bending energy of this loop. It is possible to examine the loop hamiltonian exactly. While a lot is known about axially symmetrical closed configurations , there is very little non-perturbative knowledge of the equilibria which exist when axial symmetry is broken: the loop configurations provide an analytic point of entry.

An additional application of the planar problem was pointed out recently by Willmore: any closed loop can be exploited to generate an axially symmetric toroidal geometry. The equilibria of the loop hamiltonian on a surface of constant negative curvature can be mapped into equilibria of the Helfrich functional .

The feature of our model we will exploit is its integrability, very well disguised when the problem is cast with respect to the variables embedding the loop in the plane. Indeed, the determination of the curvature at equilibrium can be reduced to the study of the motion of a particle in a one-dimensional quartic potential. The curvature can be solved as a quadrature in terms of elliptic integrals. Remarkably, we discover that the problem possesses a second, far less-obvious, level of integrability: the embedding can be expressed as a local function of the curvature without any need of the two integrations one would have expected. As a consequence, the equilibrium configuration is given by a geometrical construction of pythagorean simplicity.

The condition of closure of the loop results in a discrete spectrum for a given constrained perimeter and area. We show that scaling can be exploited to fix one of these, the perimeter say. One branch of the spectrum then consists of configurations with an $`n`$-fold symmetry ($`n=2,3,\mathrm{}`$) which inflate into a circle as the enclosed area is increased. If the area is evacuated these configurations will eventually self-intersect, tending to a limiting configuration of an inverted (i.e., negative area) circle, winding $`n1`$ times and decorated with $`n`$ infinitesimally small circles. The latter dominate, in this limit, the energy and give rise to a pole in the energy versus area diagram, the residue of which we evaluate exactly for each $`n`$.

The remaining branches of the spectrum correspond to rather complicated self-intersecting pairs which bifurcate from a limiting configuration. While self-intersecting configurations are undoubtedly of interest mathematically, they are not of primary relevance to the physical problem we pose — we determine sufficient conditions for avoiding them (nevertheless, we describe a physical situation where they may arise).

In general, one would expect equilibria connecting perturbations to be governed by a fourth order differential equation. We show, however, that this equation can be reduced (by three orders) to a first order one. The latter is singular at the circular configuration, where bifurcation occurs. Analytical approximations for the energy as a function of the area are derived, in the limit of sinusoidal perturbations of the circle of frequency $`n`$. Combined with the pole mentioned above, these expressions provide a reliable sketch of the energy – area relationship for each $`n`$.

The paper is organized as follows: Sect. 2 introduces the model and the shape equation that determines equilibria. We show how it can be reduced to the motion of a fictitious classical particle in a quartic potential. We address the behaviour under scaling of the shape equation and we show that the problem has essentially one free parameter. Moreover we analyze the loop statics; this provides a physical interpretation for the shape equation itself. In Sect. 3 we study the configuration space for this model. In particular, we describe the angle $`\mathrm{\Theta }_0`$ by which the normal to the loop gets rotated in a full oscillation of the fictitious particle in the potential. Equilibria connecting perturbations is the subject of Sect. 4. These are determined by a complicated fourth order differential equation, which is used to obtain the purely geometrical construction mentioned above. In turn, the latter permits the reduction of the order of the original equation from four to one, as well as the derivation of a sufficient condition for non-self-intersections. Finally, we obtain analytical expressions for the energy as a function of the area, in some limiting regimes, which allow for a reliable sketch of its behaviour. In the appendices we collect various expressions, useful in the calculation of variations, we derive a recursion relation for the average of the powers of $`K`$ and comment on an interesting Legendre transform in the space of parameters of the model.

## 2 First Considerations

### 2.1 Energy functional

A closed loop, parametrized by $`s[0,1]`$, is described by the embedding in the plane,

$$\stackrel{}{x}=\stackrel{}{X}(s).$$

The arc length $`\mathrm{}`$ along the loop is given by

$$\mathrm{}(s)=_0^s𝑑s^{}\left(\frac{d\stackrel{}{X}}{ds^{}}\frac{d\stackrel{}{X}}{ds^{}}\right)^{1/2}.$$
(1)
The most general expression for the energy of the configuration $`X(s)`$ which i) does not depend on the parametrization, ii) involves no higher than two derivatives of $`\stackrel{}{X}`$ and iii) is quadratic in these derivatives, is given by

$$F[\stackrel{}{X}]=\alpha 𝑑\mathrm{}K^2,$$
(2)
where $`K`$ is the geodesic curvature, equal to the inverse of the radius of curvature at each point of the loop. We will take henceforth $`\alpha `$ to be equal to 1 — $`\alpha `$ is dimensionful, unlike its two-dimensional analogue, and this means we are measuring length in units of $`\alpha `$. Let $`\widehat{t}`$ be the unit tangent to the loop (transversed counterclockwise) and $`\widehat{n}`$ its outwards normal. Then we have the Frenet-Serret equations for a plane curve333Notice that the opposite sign convention for $`K`$ is quite common in the literature. (the prime denotes a derivative with respect to arc length $`\mathrm{}`$),

$$\widehat{n}^{}=K\widehat{t},\widehat{t}^{}=K\widehat{n}.$$
(3)
In terms of the angle $`\mathrm{\Theta }`$ that $`\widehat{n}`$ makes with, say, the $`x`$-axis, $`K=\mathrm{\Theta }^{}`$.

To implement the constraints of fixed length and enclosed area, we introduce the constrained functional:

$$F_c[X]=F[X]+\mu \left(𝑑\mathrm{}L\right)+\sigma \left(_{\mathrm{int}}d^2xA\right).$$
(4)
$`\mu `$ appears as a lagrange multiplier enforcing the constraint fixing the length of the loop to some value, $`L`$. In the same way, $`\sigma `$ is associated with the constraint fixing the enclosed area to the value $`A`$. We note that, in general,

$$\mu =\frac{F_c(L,A)}{L},\sigma =\frac{F_c(L,A)}{A}.$$
(5)
In particular, if the area constraint is relaxed so that $`\sigma =0`$, then $`F_c(L,A)/A=0`$. Looking at (5), one might be led to identify $`\mu `$, $`\sigma `$ with the tension and differential pressure on the loop but this is only half true — we discuss the physical meaning of these parameters in Sec. 2.5. We have not included, in our energy functional, a term proportional to the integral of $`K`$. This is because such a term is simply the total angle by which $`\widehat{n}`$ gets rotated in transversing the loop once, equal to $`2\pi `$ for a non-self-intersecting loop.

We emphasize that the bending energy is not scale invariant (notice that neither is the dimensionally reduced Helfrich hamiltonian). An unconstrained circular loop will expand without limit. The natural scale invariant expression, $`𝑑\mathrm{}|K|`$, though not in itself a topological invariant, does have a vanishing Euler-Lagrange derivative almost everywhere. Taken as an energy functional, its constrained configurations will be arcs of a circle (solutions of $`\mu K=\sigma `$ ) joined together at curvature discontinuities in such a way as to mimimize the corresponding constrained energy function. This model will be considered elsewhere.

### 2.2 Shape equation

The Euler-Lagrange equations follow from extremizing $`F_c[X]`$,

$$\frac{\delta F_c}{\delta X^\mu }=0.$$
(6)
The reparametrization invariance of $`F_c[X]`$ implies that the tangential projection of these equations is vacuous. The normal projection gives

$`2K^{\prime \prime }+K^3\mu K\sigma =0`$ (7)
(see the appendix for some relevant formulae). We note that with the identification, $`\mu \mu +K_0^2`$, Eq. (7) coincides with Eq. (7) of  (see also ). Is is worth pointing out that a loop on a surface of constant Gaussian curvature $`G`$ satisfies the above equation with $`\mu \mu 2G`$.

A plane curve is determined, up to rigid motions, by its curvature . As is well known, our one-dimensional system is completely integrable. This is obvious because we have cast derivatives with respect to arclength. With respect to an arbitrary parametrization of the loop, first derivative terms associated with the one-dimensional Laplacian appear which conceal this fact. Writing (7) in the form

$$K^{\prime \prime }=\frac{d}{dK}V(K),V(K)\frac{1}{8}K^4\frac{\mu }{4}K^2\frac{\sigma }{2}K,$$
(8)
we map our problem onto the problem of determining the motion of a fictitious particle in a quartic potential, $`K`$ being the displacement of the particle and $`\mathrm{}`$ playing the role of time. The total energy $`E`$ of the particle is conserved,

$$\frac{dE}{d\mathrm{}}=0,E\frac{1}{2}K_{}^{}{}_{}{}^{2}+V(K),$$
(9)
a fact that permits the expression of $`K^{}`$ in terms of $`K`$ and hence, of the arclength along the loop as an integral over $`K`$ (via $`d\mathrm{}=dK/K^{}`$)

$$\mathrm{}=\frac{dK}{\sqrt{2(EV(K))}}.$$
(10)
We emphasize the difference between the energy $`E`$ of the fictitious particle, on the one hand, and the bending energy $`F`$ of the loop, on the other. In particular, the configuration of least ‘energy’ in the analogue is not the minimum of the hamiltonian.

The motion of the particle in the potential $`V`$ is periodic. Closed loop configurations consist of an integer number of identical segments, each one corresponding to one full oscillation of the particle, and hence made up itself of two symmetric halves. Then, the condition of closure of the loop can be expressed as

$$\mathrm{\Theta }_0=\frac{2\pi }{n},n=2,3,\mathrm{}$$
(11)
where $`\mathrm{\Theta }_0`$ is the angle by which the normal $`\widehat{n}`$ gets rotated in one full oscillation of the particle, given by

$$\mathrm{\Theta }_0=2_{K_{min}}^{K_{max}}\frac{KdK}{\sqrt{2(EV(K))}}$$
(12)
($`K_{min}`$, $`K_{max}`$ denote the turning points). The equilibrium configurations have $`n`$-fold symmetry and a well-defined center. The value $`n=1`$ is omitted in (11) because it is special – see below. We examine in detail the resulting configurations in Sect. 3.

### 2.3 Scaling

Let us examine the behavior of Eq. (7) under scaling of the position vector, $`\stackrel{}{X}\lambda \stackrel{}{X}`$. We find

$$\mathrm{}\lambda \mathrm{},\frac{d}{d\mathrm{}}\frac{1}{\lambda }\frac{d}{d\mathrm{}},K\frac{1}{\lambda }K.$$
(13)
Let $`\stackrel{}{X}(\mathrm{})`$ correspond to some given solution of Eq. (7), with parameter values $`(\sigma _1,\mu _1)`$. Then, a scaled solution with the same shape, $`\lambda \stackrel{}{X}(\lambda \mathrm{})`$ and with ($`L_1`$, $`A_1`$, $`F_1`$, $`E_1`$) $``$ ($`\lambda L_1`$, $`\lambda ^2A_1`$, $`\lambda ^1F_1`$, $`\lambda ^4E_1`$), is obtained by rescaling the multipliers as follows:

$$\mu _1\frac{1}{\lambda ^2}\mu _1\mu _\lambda ,\sigma _1\frac{1}{\lambda ^3}\sigma _1\sigma _\lambda .$$
(14)
Eliminating $`\lambda `$ we find the orbits of scaling in the $`(\sigma ,\mu )`$-plane

$$\mu _\lambda =\frac{\mu _1}{\sigma _1^{2/3}}\sigma _\lambda ^{2/3}.$$
(15)
Furthermore, an inversion in the origin, $`X^\mu X^\mu `$ can be identified with a rescaling by $`\lambda =1`$, which maps a solution with a given $`\sigma `$ into one with $`\sigma `$. It follows that, for the purpose of identifying distinct configurations at least, one can set, e.g. $`\sigma =1`$ and scan the essentially one-dimensional configuration space by varying $`\mu `$. For each value of $`\mu `$, which fixes the form of the potential $`V(K)`$, one has still to vary the energy $`E`$ of the particle in the well so as to satisfy (11), resulting, in general, in an infinite discrete spectrum “above” the point $`(1,\mu )`$ in the $`(\sigma ,\mu )`$-plane.

### 2.4 Windings

The closed configurations we wish to classify possess two topological invariants, both given as the winding numbers of maps $`S^1S^1`$. The first such map is the normal map of the loop, associating with each value of the parameter $`s`$ (which ranges over the first $`S^1`$ above) the corresponding value of $`\widehat{n}`$ (which ranges over the unit circle, the second $`S^1`$ above). For each rotation of $`s`$, $`\widehat{n}`$ will generally rotate $`m`$ times, with $`m`$ an integer. Configurations with distinct $`m`$’s are topologically isolated, the physical implementation being the infinite energy barrier associated with the move . Non-self-intersecting loops correspond to normal maps with index $`m=\pm 1`$ (the converse is not true).

The second map is from the $`s`$-circle to the closed orbit ($`S^1`$) traced by the phase point in the plane $`(K,K^{})`$. For a configuration corresponding to $`n`$ complete oscillations of the particle, the phase point goes around the closed curve $`\frac{1}{2}K_{}^{}{}_{}{}^{2}+V(K)=E`$ $`n`$ times, Fig. 1. There is a finite energy barrier separating configurations with distinct $`n`$’s which prevents transitions between them.

### 2.5 Forces and torques

What can we learn about the loop configurations by looking at the equilibrium of forces and torques on an infinitesimal loop segment? First we ask, why would a rigid loop want to have area less than the maximal allowed by its length? One setup that supplies an answer is to imagine that our loop is actually the cross-section of an infinite cylinder, the interior and exterior of which are filled with an incompressible fluid. We start by filling the cylinder to its maximum capacity, this gives a circular loop. We then take out some fluid from the interior while making sure that the external pressure is large enough so as to prevent the formation of bubbles in the interior. Notice that, because of its rigidity, the cylinder takes up some part of the exterior pressure and only transmits to its interior a fraction of it, it’s the differential pressure that crumbles the cylinder walls. The second of (5) then points to the identification of $`\sigma `$ with this differential pressure, taking into account that the latter points inwards. One has to be careful in applying the same argument to $`\mu `$. The first of (5) seems to suggest that $`\mu `$ is the tension of the loop but this would presuppose the possibility of tangential deformations — there is nothing in our energy functional that tells us how much energy these cost and, indeed, we have already assumed that the loop cannot be stretched or compressed tangentially. The derivative in the first of (5) is computed by comparing distinct loops with infinitesimally differing lengths, it does not refer to the deformation of a single loop. With this in mind, we now turn to loop statics.

We denote by $`\stackrel{}{T}(\mathrm{}_0)`$ the total force from the segment of the loop with $`\mathrm{}<\mathrm{}_0`$ to the one with $`\mathrm{}>\mathrm{}_0`$ and by $`\tau (\mathrm{}_0)`$ the corresponding torque. The latter is equal to $`2K`$ — this follows from our normalization of the bending energy $`F`$. Then, balancing the torques on a segment extending from $`\mathrm{}`$ to $`\mathrm{}+d\mathrm{}`$ we get

$$\tau (\mathrm{})\tau (\mathrm{}+d\mathrm{})+T_n(\mathrm{}+d\mathrm{})d\mathrm{}=0\tau ^{}T_n=0,$$
(16)
where $`T_n`$ is the normal force, the reference point was taken at $`\mathrm{}`$ and the sign conventions are shown in Fig. 2.

The equilibrium of normal forces on the segment gives

$$T_n(\mathrm{})T_n(\mathrm{}+d\mathrm{})+T_t(\mathrm{}+d\mathrm{})K(\mathrm{})d\mathrm{}\sigma d\mathrm{}=0T_n^{}+KT_t\sigma =0,$$
(17)
while the tangential components give

$$T_t(\mathrm{})T_t(\mathrm{}+d\mathrm{})T_n(\mathrm{}+d\mathrm{})K(\mathrm{})d\mathrm{}=0T_t^{}+KT_n=0.$$
(18)
Notice how, in (17) and (18), the curvature $`K`$ is responsible for the tangential force at $`\mathrm{}+d\mathrm{}`$ contributing a normal component at $`\mathrm{}`$ and vice-versa. Solving the above system of equations we find

$`\tau =2K,T_n=2K^{},T_t=K^2\mu ,`$ (19)
where, in the third relation, the integration constant was fixed to the value $`\mu `$ by the requirement that one recover the differential equation for $`K`$, Eq. (7). We see from (19) that $`\mu `$ is the tension of the loop at its inflection points, if any. Notice also how (7) is obtained upon substitution of the last two of (19) in (17), thereby identifying the physical origin of each of the terms in the former. In a forthcoming publication, where we extend our considerations to the case of a loop in space, it is shown how the above expressions for the forces and torque follow, in a model independent way, from an application of Noether’s theorem.

Our submerged cylinder model for the loop leaves no room for self-intersecting configurations, could there be any use for these? Imagine the loop made of superconducting material, in the presence of a uniform magnetic field perpendicular to its plane. In the limit where the magnetic flux due to self-inductance is negligible compared to the one due to the external field, the area of the loop has to be constant to keep the flux constant. One can then adjust the area by changing the magnetic field and, for small enough areas, self-intersecting configurations will appear.

## 3 Configurations

### 3.1 Qualitative remarks

For the purposes of this section, we may set, as explained above, $`\sigma =1`$ in the expression for the potential $`V(K)`$. Its critical points are given by the zeros of its derivative, i.e., by the roots $`K_i`$, $`i=1,2,3`$, of

$$K^3\mu K1=0.$$
(20)
We find

$$K_1=K_++i\sqrt{3}K_{},K_2=K_+i\sqrt{3}K_{},K_3=2K_+,$$
(21)
where

$$K_+\frac{m}{12}+\frac{\mu }{m},K_{}\frac{m}{12}\frac{\mu }{m},m^3108+108\sqrt{1\frac{\mu ^3}{\mu _0^3}},$$
(22)
where $`\mu _03/2^{2/3}`$. Of these, $`K_3`$ is always real while the first two are real only for $`\mu \mu _0`$. Notice that, for any $`\mu `$, the effect of the quadratic and the linear term in the potential is relatively important only in a neighborhood of the origin. For motions of the particle with sufficiently high energy $`E`$, the time spent by the particle in this region is negligible and the quartic term dominates the motion. As a result, we may conclude that $`\mathrm{\Theta }_0`$ approaches zero with increasing $`E`$, due to the approximate symmetry $`KK`$ of the motion. The term linear in $`K`$ spoils this symmetry and, for $`\sigma >0`$ (as taken in this section), makes $`\mathrm{\Theta }_0`$ slightly positive for large $`E`$. As $`\mathrm{\Theta }_0`$ descends, with increasing $`E`$, from this positive value to 0, it will cross all critical values $`2\pi /n`$ for $`n`$ greater than some $`n_0`$. We expect therefore to encounter configurations with arbitrarily high $`n`$, for all $`\mu `$ — see Fig. 4.

On the other extreme, when the energy $`E`$ is only slightly greater than a local minimum of the potential (at $`K_1`$ or $`K_3`$), and the particle oscillates around that minimum, we may approximate $`V(K)`$ by a quadratic expression in $`K`$ and find for $`\mathrm{\Theta }_0`$ the limiting value

$$\mathrm{\Theta }_{0i}^{(max)}=2\pi \frac{\sqrt{2}}{\sqrt{3\frac{\mu }{K_{i}^{}{}_{}{}^{2}}}},i=1,3,$$
(23)
where the superscript (max) is used because, as we shall see shortly, this is actually a maximum of $`\mathrm{\Theta }_{0i}(\mu ,E)`$, for fixed $`\mu `$. Starting from (21), one infers the limiting values

$$\mathrm{\Theta }_{01}^{(max)}(\mu \mu _0^+)=\mathrm{},\mathrm{\Theta }_{01}^{(max)}(\mu \mathrm{})=2\pi ,$$
(24)
as well as

$$\mathrm{\Theta }_{03}^{(max)}(\mu \mathrm{})=0,\mathrm{\Theta }_{03}^{(max)}(\mu =0)=2\pi \sqrt{\frac{2}{3}},\mathrm{\Theta }_{03}^{(max)}(\mu \mathrm{})=2\pi .$$
(25)
A plot of $`\mathrm{\Theta }_{0i}^{(max)}(\mu )`$ is given in Fig. 3. We notice that, for $`\mu >\mu _0`$, there are no non-self-intersecting configurations corresponding to oscillations in the left well.

We explain now why the value $`n=1`$ was omitted in (11). The reasoning behind that relation was that a closed configuration should correspond to $`n`$ complete oscillations of the particle. This is not necessarily so in the case of the circle. Any number of “oscillations” (of zero amplitude) will fit into a circle, including irrational numbers, and indeed, (23) shows that the circle corresponding to the particle resting at $`K_i`$, $`i=1,3`$, is made of $`\sqrt{3\frac{\mu }{K_{i}^{}{}_{}{}^{2}}}/\sqrt{2}`$ complete oscillations.

To examine what happens for values of $`E`$ that render the quadratic and linear terms important, we analyze the cases $`\mu <\mu _0`$ and $`\mu \mu _0`$ separately.

### 3.2  $`𝝁\mathbf{<}𝝁_\mathrm{𝟎}`$

The potential possesses only one minimum, at $`K_3>0`$ given by (21) (see Fig. 1). $`\mathrm{\Theta }_{03}`$ starts at $`\mathrm{\Theta }_{03}^{(max)}`$, for $`E=V(K_3)`$, and for $`\mu `$ far from $`\mu _0`$, decreases monotonically to zero with increasing $`E`$. The configurations that do not appear are those for which $`2n<2\pi /\mathrm{\Theta }_0^{(max)}`$ as well as all with negative $`n`$. This condition determines a set $`\{\mu _n|n=2,3,\mathrm{}\}`$ of critical values of $`\mu `$, such that, for $`\mu <\mu _n`$, all configurations with $`n^{}=2,3,\mathrm{},n`$ are absent (the circle is, of course, always there). As $`\mu `$ approaches $`\mu _0`$, the left wall of the potential develops a plateau that tends to the horizontal, as $`\mu \mu _0`$. The particle spends a relatively long time in this (negative $`K`$) region, which results in a negative bump in the $`\mathrm{\Theta }_0`$ curve. As $`\mu \mu _0`$, the minimum of this bump tends to $`\mathrm{}`$ — negative $`n`$ configurations appear accordingly. Direct evaluation of the integral (12) gives

$$\mathrm{\Theta }_{03}(\mu ,E)=\frac{8g(aB+bA)}{AB}\left[\alpha _2𝒦(k)+\frac{\alpha \alpha _2}{1\alpha ^2}\mathrm{\Pi }(\frac{\alpha ^2}{\alpha ^21},k)\right],$$
(26)
where $`k^2=\frac{(ab)^2(AB)^2}{4AB}`$,

$$\begin{array}{ccccccccccc}\hfill g& =& \frac{1}{\sqrt{AB}}\hfill & & \hfill \alpha & =& \frac{AB}{A+B}\hfill & & \hfill \alpha _2& =& \frac{bAaB}{aB+bA}\hfill \\ \hfill A^2& =& (ab_1)^2+a_1^2\hfill & & \hfill B^2& =& (bb_1)^2+a_1^2\hfill & & \hfill b_1& =& \frac{1}{2}(a+b)\hfill \\ \hfill a_1^2& =& \frac{4}{a+b}+ab\frac{1}{4}(a+b)^2\hfill & & \hfill a& =& K_{max}(\mu ,E)\hfill & & \hfill b& =& K_{min}(\mu ,E)\hfill \end{array}$$
(27)
and $`𝒦`$, $`\mathrm{\Pi }`$ are the complete elliptic integrals of the first and third kind respectively, given by

$$𝒦(k)=_0^{\pi /2}\frac{d\theta }{\sqrt{1k^2\mathrm{sin}^2\theta }},\mathrm{\Pi }(\alpha ^2,k)=_0^{\pi /2}\frac{d\theta }{(1\alpha ^2\mathrm{sin}^2\theta )\sqrt{1k^2\mathrm{sin}^2\theta }}.$$
(28)
The result (26) holds for $`\mu >\mu _0`$ as well, with the appropriate choice of $`K_{min}`$, $`K_{max}`$. A plot of $`\mathrm{\Theta }_{03}(\mu ,E)`$, for various values of $`\mu `$, appears in Fig. 4.

A couple of remarks are in order at this point. Consider the point $`A`$ in Fig. 4, which is the intersection of the $`\overline{\mu }=0.5`$ curve with the $`n=2`$ line, corresponding to an $`n=2`$ configuration. Imagine now that $`\overline{\mu }`$ is diminished continuously — the corresponding curve will move more or less downwards forcing $`A`$ to move to the left. The corresponding $`\overline{E}`$ then diminishes which means that the fictitious particle oscillates in the potential well with smaller amplitude. When $`A`$ hits the $`\mathrm{\Theta }_{03}`$-axis, $`\overline{E}`$ is zero, the particle sits at the bottom of the well and the corresponding configuration becomes a circle. In other words, all configurations like the one corresponding to the point $`A`$ can be continuously deformed to a circle.

This is not the case though with the points $`B_1`$, $`B_2`$ in the same figure. Consider $`B_1`$ — it corresponds to an $`n=4`$ configuration for $`\overline{\mu }=1.01`$. Imagine now that $`\overline{\mu }`$ is diminished in a continuous way. As soon as it becomes smaller then 1, the infinite negative pole of the corresponding curve is softened to a negative minimum which, for $`\overline{\mu }`$ sufficiently close to 1, still intersects the $`n=4`$ line. Given that the $`\mu `$-curves in Fig. 4 become almost horizontal at high $`E`$, we give an exagerrated sketch of the situation in Fig. 5. As the value of $`\overline{\mu }`$ is lowered, this negative minimum rises and, for a critical value of $`\overline{\mu }`$, will just touch the $`n=4`$ line, i.e., $`B_1`$ and $`B_2`$ collapse to a single point $`B`$. The configuration that corresponds to this point is not a circle, since $`\overline{E}`$ is positive and the particle oscillates with a finite amplitude, i.e., $`K`$ is not constant. We conclude that configurations corresponding to points like $`B_1`$, $`B_2`$ cannot be continuously deformed to a circle. Moreover, there exist bifurcation points, like $`B`$ above, distinct from the circle. When $`\overline{\mu }`$ is lowered even more, the minimum of the $`\overline{\mu }`$-curve becomes positive and remarks similar to the above can be made about its points of intersection $`B_1^{}`$, $`B_2^{}`$ with sufficiently high $`n`$ lines (see Fig. 5).

### 3.3  $`𝝁\mathbf{}𝝁_\mathrm{𝟎}`$

The potential posseses two local minima, at $`K_1<0`$ and $`K_3>0`$ (with $`V(K_1)>V(K_3)`$) and one local maximum at $`K_2<0`$ (see Fig. 6). Depending on its energy and where it is started from, the particle is confined in the left well, the right well or visits both during every oscillation.

We note the following, regarding the asymptotic behaviour of $`\mathrm{\Theta }_{0i}`$ as $`\mu \mathrm{}`$. In this regime, the linear term becomes negligible in $`V(K)`$ and the integral giving $`\mathrm{\Theta }_{0i}`$ can be easily seen to reduce (with $`K^2x`$) to an integral giving (half) the period of a harmonic oscillator in a potential $`x^2/4\mu x/2`$, which is independent of its energy $`E`$, as well as the ‘constant force’ $`\mu `$. Another way to see the $`\mu `$-independence is by observing that, with $`\sigma =0`$ (which is equivalent to sending $`\mu `$ to infinity), one lies on the $`\mu `$-axis in the $`(\sigma ,\mu )`$-plane, which is an orbit under scaling, hence changes in $`\mu `$ in this regime leave $`\mathrm{\Theta }_{0i}`$ invariant. Another novelty in this case is that, as the energy approaches $`V(K_2)`$ (either from above or below), the period of the motion tends to infinity, most of which the particle spends at $`K_2`$. Assume we start the particle at the turning point on the right with zero velocity and with just the right energy to reach the central maximum of $`V`$, $`E=V(K_2)`$. The corresponding curve starts in a counterclockwise sense with decreasing curvature, passes through an inflection point and acquires infinite length spiralling forever clockwise as it approaches asymptotically a circle of radius $`1/K_2`$, see Fig. 6.

$`\mathrm{\Theta }_{03}`$ tends accordingly to $`\mathrm{}`$, this corresponds to the negative poles in Fig. 4. As mentioned already, for $`V(K_1)<E<V(K_2)`$, $`\mathrm{\Theta }_0`$ will have two branches, $`\mathrm{\Theta }_{01}`$ and $`\mathrm{\Theta }_{03}`$, corresponding to the motions confined to either well. Both of these branches will tend to $`\mathrm{}`$ as $`EV(K_2)`$. For $`E>V(K_2)`$ there is only one branch, which starts from $`\mathrm{}`$ and approaches asymptotically zero (as $`E\mathrm{}`$), after having reached a positive local maximum. A three-dimensional plot of $`\mathrm{\Theta }_{03}(\mu ,E)`$ is given in Fig. 7. We present a representative collection of configurations in Fig. 8. We also show, for reference purposes, some self-intersecting configurations in Fig. 9.

## 4 Perturbations

### 4.1 Equilibrium-connecting deformations

We study here the following problem: given an equilibrium configuration $`\stackrel{}{X}`$, find an infinitesimal deformation $`\delta _ϵ\stackrel{}{X}=ϵ\widehat{n}`$ such that $`\stackrel{}{X}+\delta _ϵ\stackrel{}{X}`$ describes a nearby equilibrium. Notice that $`\delta _ϵ\stackrel{}{X}`$ is not necessarily the displacement vector of some point of the loop — the latter will generally have a tangential component as well. Our motivation is to obtain new solutions from known ones, with the closure condition automatically satisfied.

Using some of the formulae listed in the appendix, the linearization of the equilibrium condition (7) gives that the deformation $`ϵ`$ must satisfy

$`DE_4(ϵ)`$ $``$ $`2ϵ^{\prime \prime \prime \prime }+(5K^2\mu )ϵ^{\prime \prime }+10KK^{}ϵ^{}+(12E{\displaystyle \frac{5}{2}}K^4+6\mu K^2+10\sigma K)ϵ`$ (29)
$`=`$ $`\delta _ϵ\mu K\delta _ϵ\sigma ,`$
where the deformed configuration satisfies (7), with $`\mu \mu +\delta _ϵ\mu `$, $`\sigma \sigma +\delta _ϵ\sigma `$. The energy will change as well, $`EE+\delta _ϵE`$. This appears explicitly in the linearization of (8), which gives

$`DE_3(ϵ)`$ $``$ $`2K^{}ϵ^{\prime \prime \prime }+(K^3\mu K\sigma )ϵ^{\prime \prime }+2K^2K^{}ϵ^{}+K(12E{\displaystyle \frac{1}{2}}K^4+2\mu K^2+5\sigma K)ϵ`$ (30)
$`=`$ $`{\displaystyle \frac{1}{2}}K^2\delta _ϵ\mu K\delta _ϵ\sigma \delta _ϵE.`$
### 4.2 The $`𝝈`$-identities

We introduce a number of functions on the loop, relevant in the study of the solutions of (29). Any vector field $`\stackrel{}{a}`$ defines the following two functions on the loop

$$a_n\stackrel{}{a}\widehat{n},a_t\stackrel{}{a}\widehat{t}.$$
(31)
The position vector $`\stackrel{}{X}`$, in particular, supplies

$$h\stackrel{}{X}\widehat{n},p\stackrel{}{X}\widehat{t},$$
(32)
both functions evidently dependent on the choice of origin. One finds

$$a_n^{}=Ka_t,a_t^{}=Ka_n,h^{}=Kp,p^{}=1Kh,$$
(33)
where $`\stackrel{}{a}`$ is henceforth assumed constant. We will say that a function $`f`$ on the loop is the generator of a transformation iff, under the latter, $`\delta \stackrel{}{X}\widehat{n}=f`$. With this definition, $`a_n`$, $`a_t`$ are the generators of translations along $`\stackrel{}{a}`$ and perpendicularly to it respectively. $`p`$ and $`h`$ generate rotations and dilations respectively (both w.r.t. the origin). One easily verifies that $`DE_4(a_n)=DE_4(a_t)=DE_4(p)=0`$ while $`DE_4(h)=2\mu K+3\sigma `$, in accordance with the scaling behavior of $`\mu `$, $`\sigma `$ found earlier, Eq. (14). Experimenting a little with (29) we find, not without some surprise, two more solutions

$$DE_4(K^{})=0,DE_4\left(\sigma ^1(K^2\mu )\right)=2\mu K+3\sigma .$$
(34)
Noting, additionally, that when the origin is at the center of the loop, the zeros of $`p`$ coincide with the extrema of $`K`$, as well as the coincidence of the extrema of $`h`$ and $`K^2`$, we make the ansatz

$$h=\sigma ^1(K^2\mu )+f,p=2\sigma ^1K^{}+g,$$
(35)
with $`f`$, $`g`$ to be determined. Substituting in (33), we find $`f^{}=Kg`$, $`g^{}=Kf`$, with solution

$$f=a\mathrm{cos}(\mathrm{\Theta }\varphi _0),g=a\mathrm{sin}(\mathrm{\Theta }\varphi _0),$$
(36)
where $`a`$, $`\varphi _0`$ are arbitrary constants. One recognizes $`f`$, $`g`$ to be, respectively, the functions $`a_n`$, $`a_t`$ corresponding to the position vector $`\stackrel{}{a}`$ of the center of the loop, then (35) state

$`h=\sigma ^1(K^2\mu )+a_n,p=2\sigma ^1K^{}+a_t,`$ (37)
a remarkable geometrical property, the implications of which will occupy us shortly. Notice that neither of (37) survives in the $`\sigma 0`$ limit. Unless otherwise stated, we will take $`\stackrel{}{a}`$ equal to zero in what follows.

One might wonder whether (29) admits other polynomial solutions in $`K`$, apart from the second of (34). To investigate this, we rewrite (29) in terms of derivatives w.r.t. $`K`$, using $`\frac{d}{d\mathrm{}}=K^{}\frac{d}{dK}`$ and find for $`DE(K^n)`$ the leading term $`(\frac{1}{8}n^4+\frac{3}{4}n^3+\frac{1}{8}n^23n\frac{5}{2})K^{n+4}`$, with only positive integer root $`n=2`$. We conclude that no other polynomial solutions of (29) exist.

### 4.3 Geometrics

We are now in a position to give a purely geometrical construction of the equilibria. Indeed, starting from $`X^2=h^2+p^2`$ and using (37), we find

$`X^2X_0^2=4\sigma ^1K,X_0\sigma ^1\sqrt{8E+\mu ^2},`$ (38)
where, in the derivation, use was made of the first integral, Eq. (9). This remarkable formula expresses the embedding completely in terms of $`K(\mathrm{})`$. We emphasize that the shape is obtained directly from $`K`$, given $`(\sigma ,\mu ,E)`$, without any integration. Note that both the relations (37) follow by taking derivatives with respect to arclength of (38). It follows from (38) that the moment of inertia $``$ of the loop around an axis perpendicular to its plane and passing through its center, is given by

$$=LX_0^2+8\pi \sigma ^1.$$
(39)
To construct geometrically the equilibrium curve (for given $`\sigma `$, $`\mu `$, $`E`$), we draw a circle with radius $`X_0`$ and from a point $`P`$ outside it, bring the tangent $`PA`$ to the circle (see Fig. 10). Taking as base the square of this tangent, we construct a parallellepiped with volume $`4\sigma ^1`$, the height of which is the radius of curvature of the equilibrium curve at $`P`$ — this follows from (38), written in the form $`(X^2X_0^2)\rho =4\sigma ^1`$, where $`\rho `$ is the radius of curvature at $`P`$. To achieve closure of the resulting curve, one has to start with a particular slope, given by either of (37). In this way, one constructs the part of the curve lying outside the circle (the latter intersects the curve at its inflection points, if any) — the interior part, present only if the ficticious particle reaches into the negative $`K`$ region, is constructed in a similar manner.

Referring back to our expression for the force $`\stackrel{}{T}(\mathrm{})`$, Eq. (19), we realize that (37) implies that $`\stackrel{}{T}(\mathrm{})`$ is orthogonal to $`\stackrel{}{X}(\mathrm{})`$ while $`T(\mathrm{})=\sigma X(\mathrm{})`$, i.e.,

$`\stackrel{}{T}(\mathrm{})=\sigma (\widehat{n}\times \widehat{t})\times \stackrel{}{X}(\mathrm{}),T_n=\sigma p,T_t=\sigma h,`$ (40)
with $`\stackrel{}{T}(\mathrm{})`$ the force on the part of the loop pointed to by $`\widehat{t}`$. Notice that $`\stackrel{}{T}`$ is generally compressing but changes nevertheless to true tension at the points were $`\stackrel{}{X}`$ is tangential to the loop. Also, $`p`$ has acquired a direct physical interpretation as a consequence of (40): $`\sigma pd\mathrm{}`$ is just the torque, w.r.t. the origin, of the force due to the pressure on a segment of the loop with length $`d\mathrm{}`$. Then the vanishing of the total torque on the loop is guaranteed by the fact that $`p`$ is a derivative. In fact, one may derive a compact and rather pleasing formula for the torque due to pressure on any segment of the loop, like the one defined by $`A`$, $`B`$ in the sketch of an $`n=4`$ configuration of Fig. 11.

Suppose we take $`A`$ as reference point, then (37) is valid with $`\stackrel{}{a}`$ connecting $`A`$ with the center $`𝒪`$ of the loop. We get

$$\tau _{AB}^{(A)}=\sigma _A^B𝑑\mathrm{}p=\sigma _A^B𝑑\mathrm{}(2\sigma ^1K^{}+a_t)=\mathrm{\hspace{0.17em}2}(K_BK_A)+\sigma \stackrel{}{a}_A^B\stackrel{}{d\mathrm{}}$$

i.e.,

$`\tau _{AB}^{(A)}=2(K_BK_A)+\sigma \stackrel{}{a}\stackrel{}{AB}.`$ (41)
Moving the reference point simply moves one endpoint of $`\stackrel{}{a}`$. For $`\stackrel{}{a}=0`$ (torque w.r.t. $`𝒪`$), $`\tau _{AB}^{(O)}=2(K_BK_A)`$.

### 4.4 Self-intersections

Given that some potential applications of our model exclude self-intersecting configurations, we look now for sufficient conditions for non-self-intersection. We restrict our attention to self-intersections that can be reached by a continuous deformation of non-self-intersecting configurations. In other words, we consider a one-parameter family of configurations $`\stackrel{}{X}(\mathrm{},t)`$, $`t[0,1]`$, continuous in $`t`$, such that, for every $`t`$, the corresponding curvature satisfies (7) and we take $`\stackrel{}{X}(\mathrm{},0)`$, $`\stackrel{}{X}(\mathrm{},1)`$ to be non-self-intersecting and self-intersecting respectively. Then, we observe that as $`t`$ varies from 0 to 1, one necessarily encounters a “kiss”:  The position vector $`\stackrel{}{X}_A`$ is along the axis of symmetry of the lobe and tangential to the loop at $`A`$. The third of (40) then shows that the force at $`A`$ is normal to the loop, while the second of (40) gives its magnitude as

$$T_A=\sigma X_A.$$
(42)
The tangency of $`\stackrel{}{X}`$ at $`A`$ implies that $`h_A=0`$ and hence, using (37) once more, we get (assuming that $`K`$ is negative at a kiss444This assumption is true for all configurations we have studied numerically.)

$`K_A=\sqrt{\mu },`$ (43)
for all kisses, regardless of the order of the configuration. We conclude that a sufficient condition for non-self-intersection (of the type defined above), is $`\mu <0`$. (43) gives $`K_{}^{}{}_{A}{}^{2}`$ in terms of $`\sigma `$, $`\mu `$, $`E`$

$$K_{}^{}{}_{A}{}^{2}=2\left(E+\frac{1}{8}\mu ^2\frac{1}{2}\sigma \sqrt{\mu }\right).$$
(44)
On the other hand, from (38) we find that

$$X_A^2=\sigma ^2(8E+\mu ^2)4\sigma ^1\sqrt{\mu }.$$
(45)
The tangency of $`X_A`$ though means that $`2\sigma ^1K_A^{}=p_A=X_A`$ and the above two equations then give (setting $`\sigma =1`$)

$$E^2+\frac{1}{4}\left(\mu ^24\sqrt{\mu }\frac{1}{2}\right)E+\frac{1}{64}(\mu ^24\sqrt{\mu })(\mu ^24\sqrt{\mu }1)=0,$$
(46)
with roots

$$E_1=\frac{1}{2}\sqrt{\mu }\frac{1}{8}\mu ^2+\frac{1}{8},E_2=\frac{1}{2}\sqrt{\mu }\frac{1}{8}\mu ^2.$$
(47)
For the particular case $`n=2`$, $`K_A^{}=0`$ and we get $`E=E_2`$. We have seen in Sec. 3 that the relation $`\mathrm{\Theta }_0=2\pi /n`$ defines a curve in the $`(\mu ,E)`$-plane, consisting of all parameter pairs giving rise to a configuration of order $`n`$. The intersection of that curve with the ones just written above, consists of the points $`(\mu ,E)`$ giving rise to a kissing configuration of order $`n`$.

Consider how (40), (41) guarantee equilibrium in some particular examples. First, look at an even-$`n`$ configuration, say, $`n=2`$:  The total force due to the pressure, pushing together the two halves in the sketch, is $`\sigma |\stackrel{}{AB}|`$. (40) says that the tension at $`A`$, $`B`$ (purely tangential, compressing) is $`\sigma |\stackrel{}{AB}|/2`$, thus leaving each half at rest. As a second example, consider the lobe defined by a self-intersection:  The total force due to the pressure on the lobe is zero, and (40) says that the forces from the rest of the loop, on the two ends of the lobe that meet at $`A`$, are opposite (with direction so as to keep the lobe closed). As a further check on our results, one can verify the balancing of torques using (41) on, say, the right half of a kiss.

### 4.5 Connecting equilibria

There is more to be derived from (38). Taking $`\delta _ϵ`$ on both sides, we find

$$DE_2^{(a)}(ϵ)2ϵ^{\prime \prime }+(3K^2\mu )ϵ=\frac{\sigma }{2}\delta _ϵX_0^22K\sigma ^1\delta _ϵ\sigma ,$$
(48)
a considerable improvement over (29)555Notice, however, that (48) is only valid when the center of the loop is at the origin.. Moreover, by varying the first of (37) we also find

$$DE_2^{(b)}(ϵ)2Kϵ^{\prime \prime }2K^{}ϵ^{}+(\sigma +2K^3)ϵ=\sigma ^1(K^2\mu )\delta _ϵ\sigma \delta _ϵ\mu .$$
(49)
Comparison with (48) leads to a first order equation for $`ϵ`$

$$DE_1(ϵ)4K^{}ϵ^{}4K^{\prime \prime }ϵ=c_2K^2+c_1K+c_0,$$
(50)
where the constants $`c_i`$ are given by

$$c_22\sigma ^1\delta _ϵ\sigma ,c_1\sigma \delta _ϵX_0^2,c_02\delta _ϵ\mu 2\mu \sigma ^1\delta _ϵ\sigma .$$
(51)
Noting that $`\frac{d}{dK}=(K^{})^1\frac{d}{d\mathrm{}}`$, (50) can be written in the form

$$8(EV)\dot{ϵ}+4\dot{V}ϵ=c_2K^2+c_1K+c_0,$$
(52)
where the dot denotes differentiation w.r.t. $`K`$ and $`V=V(K)`$. This latter equation is readily integrated to give

$`ϵ={\displaystyle \frac{1}{4}}K^{}{\displaystyle 𝑑K\frac{c_2K^2+c_1K+c_0}{\left(2(EV)\right)^{3/2}}}.`$ (53)
One may add an arbitrary amount of $`K^{}`$ (rotation) to this (but not $`a_n`$ or $`a_t`$, since these move the origin). To get particular solutions from (53), we need to specify the direction of the deformation in the $`(\sigma ,\mu )`$-plane, i.e., the ratio $`\delta _ϵ\mu /\delta _ϵ\sigma `$. It will prove convenient for our further analysis of (53), to reparameterize the $`(\sigma ,\mu )`$-plane introducing new coordinates $`(\lambda ,\xi )`$ via

$$\lambda (\sigma ,\mu )=\sigma ^{1/3},\xi (\sigma ,\mu )=\sigma ^{2/3}\mu .$$
(54)
A point $`P`$ with coordinates $`(\sigma ,\mu )(0,0)`$ lies on a unique scaling orbit, which can be specified by the $`\mu `$-coordinate of its point of intersection $`P^{}`$ with the $`\sigma =1`$ line, this is the value of $`\xi `$ for $`P`$. One can get now from $`P^{}`$ to $`P`$ by scaling by $`\lambda `$ (see Fig. 12). The obvious advantage of these new coordinates is that the scaling orbits are constant-$`\xi `$ lines. The dependence on $`\lambda `$ of an arbitrary quantity $`S(\lambda ,\xi )`$, with length dimension $`q`$, is $`S(\lambda ,\xi )=\lambda ^q\stackrel{~}{S}(\xi )`$, where we denote by a tilde the remaining function of $`\xi `$. It follows that $`\frac{S}{\lambda }=q\lambda ^1S`$ so that $`\delta _ϵS=q\lambda ^1S\delta _ϵ\lambda +\lambda ^q\stackrel{~}{S}^{}\delta _ϵ\xi `$, the prime denoting here differentiation w.r.t. $`\xi `$.

We verify that (53) gives $`ϵh`$ for scaling. Eq. (51), written in terms of $`(\lambda ,\xi )`$, gives

$`c_2`$ $`=`$ $`6\lambda ^1\delta _ϵ\lambda `$
$`c_1`$ $`=`$ $`\lambda ^2\left(16\stackrel{~}{E}+2\xi ^2\right)\delta _ϵ\lambda +\lambda ^1\left(8\stackrel{~}{E}^{}+2\xi \right)\delta _ϵ\xi `$ (55)
$`c_0`$ $`=`$ $`2\lambda ^3\xi \delta _ϵ\lambda +2\lambda ^2\delta _ϵ\xi ,`$
where $`E(\lambda ,\xi )=\lambda ^4\stackrel{~}{E}(\xi )`$ is the energy that guarantees closure of some particular configuration. An increment $`\delta _ϵ\lambda `$ corresponds to scaling by a factor $`1+\delta _ϵ\lambda /\lambda `$, hence the corresponding $`ϵ`$ should be $`ϵ=(\delta _ϵ\lambda /\lambda )h`$. Putting $`\delta _ϵ\xi =0`$ in (4.5) one determines the $`c_i`$ for pure scaling, then computing $`d/dK((K^2\mu )/K^{})`$ one finds that the $`K^5`$ and $`K^3`$ terms in the numerator cancel and one recovers the integrand in (53) with just the right $`c`$’s.

Integrating the first of (37) w.r.t. $`d\mathrm{}`$, and using

$$A=\frac{1}{2}𝑑\mathrm{}h,$$
(56)
we get

$$F=\mu L+2\sigma A,$$
(57)
from which

$$\delta _ϵF=L\delta _ϵ\mu +\mu \delta _ϵL+2A\delta _ϵ\sigma +2\sigma \delta _ϵA$$
(58)
follows. On the other hand, direct variation of $`F`$, Eq. (2), gives

$$\delta _ϵF=\mu \delta _ϵL\sigma \delta _ϵA.$$
(59)
Comparing with (58) we find

$$L\delta _ϵ\mu +2A\delta _ϵ\sigma +3\sigma \delta _ϵA+2\mu \delta _ϵL=0.$$
(60)
For a length-preserving deformation, $`\delta _ϵL=0`$, and (60) reduces to

$$L\delta _ϵ\mu +2A\delta _ϵ\sigma +3\sigma \delta _ϵA=0.$$
(61)
We note in passing that the various differential operators we have defined above are related in the following way

$`[DE_1(ϵ)]^{}`$ $`=`$ $`2K^{}DE_2{}_{}{}^{(a)}(ϵ),`$
$`K[DE_1(ϵ)]^{}`$ $`=`$ $`2K^{}DE_2{}_{}{}^{(b)}(ϵ)+K^{}DE_1(ϵ),`$
$`DE_3(ϵ)`$ $`=`$ $`K^{}[DE_2{}_{}{}^{(a)}(ϵ)]^{}K^{\prime \prime }DE_2{}_{}{}^{(a)}(ϵ){\displaystyle \frac{1}{4}}(K^2\mu )DE_1(ϵ),`$
$`K^{}DE_4(ϵ)`$ $`=`$ $`2[DE_3(ϵ)]^{}.`$
It follows that any perturbation that satisfies Eq. (50) will necessarily satisfy all the higher order ones.

### 4.6 Perturbing the circle

Referring to our cylinder model of the loop (see Sec. 2.5), we consider here the following thought experiment: we start with the loop “filled” to capacity, i.e., in a circular configuration, and start deflating it by removing fluid from its interior while keeping its length fixed. We would like to follow the evolution of its shape, for distinct $`n`$’s, and derive the limiting form of the bending energy $`F`$ as a function of the area $`A`$ near the circular extreme. Going beyond the cylinder model, we would also like to allow for self-intersections — what happens if one just keeps subtracting area?

Consider the following perturbation to a circle of radius $`R`$

$$\rho (\theta )=R_\alpha \left(1+\alpha \mathrm{sin}(n\theta )\right),\alpha 1.$$
(62)
The line element and curvature in these coordinates are

$$d\mathrm{}=\sqrt{\rho ^2+\dot{\rho }^2}d\theta ,K=\frac{\rho ^2+2\dot{\rho }^2\rho \ddot{\rho }}{(\rho ^2+\dot{\rho }^2)^{3/2}},$$
(63)
where the dot here denotes differentiation w.r.t. $`\theta `$. Substituting (62) in the first of (63), integrating over $`\theta `$ and requiring the total length to be equal to the initial value $`2\pi R`$, we find

$$R_\alpha =R(1\frac{n^2\alpha ^2}{4})+𝒪(\alpha ^3),$$
(64)
while for the curvature the second of (63) gives

$$K=\frac{1}{R}\left(1+\alpha (n^21)\mathrm{sin}(n\theta )\right)+𝒪(\alpha ^2).$$
(65)
We substitute the above expression in the differential equation for $`K`$, Eq. (7), and demand that it be a solution to $`𝒪(\alpha )`$. The constant and $`𝒪(\alpha )`$ terms respectively give

$$1R^2\mu R^3\sigma =0,2n^4+(R^2\mu 5)n^2+3R^2\mu =0.$$
(66)
Notice that the first relation is simply the statement $`\dot{V}(1/R)=0`$ while the second also follows from (29) with $`ϵ`$ as in (62). For this perturbation the differential equation cannot be satisfied to $`𝒪(\alpha ^2)`$. From (66) we get

$$\sigma =\frac{2}{R^3}(n^21),\mu =\frac{1}{R^2}(32n^2).$$
(67)
During the deflating process, both $`\sigma `$ and $`\mu `$ will vary as functions of $`A`$. The point in the $`\sigma `$-$`\mu `$ plane corresponding to the configuration will trace out an orbit, starting at the above points (for each $`n`$), all of which lie on the line $`R^3\sigma +R^2\mu =1`$. For $`F`$ and $`A`$ we find

$$F=\frac{2\pi }{R}\left(1+\frac{1}{2}(n^21)^2\alpha ^2\right)+𝒪(\alpha ^3),A=\pi R^2\left(1\frac{1}{2}(n^21)\alpha ^2\right)+𝒪(\alpha ^3).$$
(68)
Notice that, near the circle,

$$\frac{\delta F}{\delta A}=\frac{2}{R^3}(n^21)=\sigma ,$$

in agreement with (59) (for $`\delta L=0`$). As $`A`$ keeps diminishing, our numerical analysis shows that the configurations start self-intersecting, giving rise to regions of negative area — this scenario is sketched in Fig. 13. Notice how the winding numbers ($`+1`$ for the little circles, $`1`$ for the big one) add up correctly to match that of the circle to which this configuration can be continuously deformed. This observation points to a feature of the above $`m=1`$ limiting configurations already alluded to in the introduction: the big circle in the configuration of order $`n`$ actually winds around itself $`n1`$ times. For example, to move from one lobe to the next along the loop, in the $`n=3`$ configuration, one must travel an angle of $`4\pi /3`$, not $`2\pi /3`$ (a glance at Fig. 9.b reveals how this comes about). The limiting shape in this direction then is a circular one, with radius slightly less than $`R/(n1)`$ and $`n`$ little circles attached to it. In this extreme, $`F(A)`$ is dominated by the little circles and assumes the limiting form

$$F(A)\frac{4\pi n^2}{R}\frac{1}{1+(n1)\frac{A}{\pi R^2}},$$
(69)
where the area is counted with the appropriate multiplicity due to the winding (e.g., plus three times the area of the little circle, minus twice that of the big one, in the $`n=3`$ case). A sketch of $`F(A)`$ that interpolates between (68) and (69) is given in Fig. 14.

## 5 Conclusions

In this paper, we have studied the equilibrium configurations of elastic planar loops with constant area and constant length. The condition of closure of the loop gives a discrete spectrum of configurations, lying along several branches in parameter space. We focus on non self-intersecting loops that inflate to a circle when the enclosed area is increased, due to their relevance as toy model for two-dimensional membranes. For this branch, starting from analytical expressions in the relevant limit cases, we obtain a reliable sketch of the dependence of the energy on the area. The equation that determines equilibria connecting perturbations led us, rather unexpectedly, to the $`\sigma `$-identities (37), from which a novel geometrical construction of the (constrained) elastic loop followed. In turn, these identities permit the expression of an equilibrium connecting normal deformation in terms of the curvature and the appropriate perturbations of the parameters.

In future work, we plan to address the important issue of the stability of the equilibrium configurations. In light of the complicated structure of the fourth order differential operator appearing in (29), this appears to be a non-trivial task in the general case. It is indicative of the intrinsic complexity of the question that even the relatively simple case of the figure eight configuration requires the elaborate analysis of .

Another interesting issue is the analysis of thermal fluctuations in this model, which would provide an analytical counterpart to the study of Ref.  of two-dimensional vesicles, using Monte Carlo techniques. In particular, there is a similarity between some of the configurations of Fig. 2 of the above work and the configuration of Fig. 8.a in this paper. It is also interesting to note that the shapes of Fig. 8 closely resemble a top view of the starfish vesicles of Ref. , which are almost planar. Finally, configurations similar to ours appear in , where our problem has been approached from a functional analytic point of view and important existence results have been derived, as well as in .

After the completion of this paper, we set out to study its most natural extension: geometric models for loops in space. We became aware of a large body of literature that explores the interconnections among hierarchies of functionals of the geometry of a curve, their corresponding generators of rigid motion and integrable systems, such as the KdV equation. We have found Ref.  an excellent point of entry to the subject. We also realized that our $`\sigma `$-identities can be obtained by adapting the cylindrical coordinates used by Langer and Singer  in the analysis of buckled rings. We expect that these explorations, apart from their intrinsic interest, will contribute to a deeper understanding of two-dimensional membranes.

## Acknowledgements

G.A. wishes to thank Gilberto Tavares for technical assistance and CONACyT for a graduate fellowship. R.C. was supported by CONACyT grant 32187-E. G.A., C.C. and J.G. were supported by CONACyT grant 32307-E and DGAPA-UNAM grant IN119792.

## Appendix A Some Useful Formulae

Under an infinitesimal deformation of the loop along the normal, $`\stackrel{}{X}\stackrel{}{X}+\delta _ϵ\stackrel{}{X}\stackrel{}{X}+ϵ(\mathrm{})\widehat{n}`$, we find

$$\begin{array}{ccccccccccc}\hfill \delta _ϵd\mathrm{}& =& Kϵd\mathrm{},\hfill & & \hfill \delta _ϵ\widehat{t}& =& ϵ^{}\widehat{n},\hfill & & \hfill \delta _ϵK& =& ϵ^{\prime \prime }K^2ϵ\hfill \\ \hfill \delta _ϵ\frac{d}{d\mathrm{}}& =& Kϵ\frac{d}{d\mathrm{}},\hfill & & \hfill \delta _ϵ\widehat{n}& =& ϵ^{}\widehat{t},\hfill & & \hfill \delta _ϵK^{}& =& ϵ^{\prime \prime \prime }K^2ϵ^{}3KK^{}ϵ\hfill \\ \hfill \delta _ϵh& =& ϵpϵ^{},\hfill & & \hfill \delta _ϵp& =& ϵ^{}h,\hfill & & \hfill \delta _ϵK^{\prime \prime }& =& ϵ^{\prime \prime \prime \prime }K^2ϵ^{\prime \prime }5KK^{}ϵ^{}\hfill \\ & & & & & & & & & & (6E\frac{11}{4}K^4+\frac{7}{2}\mu K^2+5\sigma K)ϵ\hfill \end{array}$$

(we have used the e.o.m and the first integral, Eqs. (7) and (9), to express $`K^{\prime \prime }`$ and $`(K^{})^2`$ in terms of $`K`$).

For the derivatives of $`h`$, $`p`$, $`K`$ we find

$`h^{}`$ $`=`$ $`Kp`$
$`h^{\prime \prime }`$ $`=`$ $`K^2h+K^{}p+k`$
$`h^{\prime \prime \prime }`$ $`=`$ $`3KK^{}h+({\displaystyle \frac{3}{2}}K^3+{\displaystyle \frac{\mu }{2}}K+{\displaystyle \frac{\sigma }{2}})p+2K^{}`$
$`h^{\prime \prime \prime \prime }`$ $`=`$ $`\left(6E+{\displaystyle \frac{15}{4}}K^4{\displaystyle \frac{7}{2}}\mu K^25\sigma K\right)h+\left({\displaystyle \frac{15}{2}}K^2K^{}+{\displaystyle \frac{\sigma }{2}}K^{}\right)p{\displaystyle \frac{5}{2}}K^3+{\displaystyle \frac{3}{2}}\mu K+{\displaystyle \frac{3}{2}}\sigma `$
$`p^{}`$ $`=`$ $`1Kh`$
$`p^{\prime \prime }`$ $`=`$ $`K^{}hK^2p`$
$`p^{\prime \prime \prime }`$ $`=`$ $`({\displaystyle \frac{3}{2}}K^3{\displaystyle \frac{\mu }{2}}K{\displaystyle \frac{\sigma }{2}})h3KK^{}pK^2`$
$`p^{\prime \prime \prime \prime }`$ $`=`$ $`({\displaystyle \frac{15}{2}}K^2K^{}{\displaystyle \frac{\mu }{2}}K^{})h+({\displaystyle \frac{15}{4}}K^4{\displaystyle \frac{7}{2}}\mu K^25\sigma K6E)p5KK^{}`$
$`K^{}`$ $`=`$ $`\sqrt{2}(E{\displaystyle \frac{1}{8}}K^4+{\displaystyle \frac{\mu }{4}}K^2+{\displaystyle \frac{\sigma }{2}}K)^{1/2}`$
$`K^{\prime \prime }`$ $`=`$ $`{\displaystyle \frac{1}{2}}(K^3+\mu K+\sigma )`$
$`K^{\prime \prime \prime }`$ $`=`$ $`{\displaystyle \frac{3}{2}}K^2K^{}+{\displaystyle \frac{\mu }{2}}K^{}`$
$`K^{\prime \prime \prime \prime }`$ $`=`$ $`{\displaystyle \frac{3}{2}}K^5{\displaystyle \frac{10}{4}}\mu K^3{\displaystyle \frac{15}{4}}\sigma K^2+({\displaystyle \frac{\mu ^2}{4}}6E)K+{\displaystyle \frac{\mu \sigma }{4}}`$
## Appendix B Averages

Consider the quantities $`W_n`$, defined by

$$W_n=_0^{L_0}K^n𝑑\mathrm{}=2_{K_{min}}^{K_{max}}K^n\frac{dK}{K^{}}.$$
(70)
$`L_0`$ here is the length along the loop corresponding to one full oscillation of the particle (we will use $`F_0`$ later with a similar meaning). $`W_n`$ is then (proportional to) the average of $`K^n`$ along the loop. Starting from

$$_{K_{min}}^{K_{max}}\frac{d}{dK}(K^nK^{})=0,$$
(71)
one may derive the recursion relation

$$\frac{1}{8}(n+2)W_{n+3}+\frac{\mu }{4}(n+1)W_{n+1}+\frac{\sigma }{4}(2n+1)W_n+EnW_{n1}=0,n=0,1,2,\mathrm{},$$
(72)
which permits, in principle, the calculation of the average of any power series in $`K`$, for a given equilibrium configuration. In particular,

$$W_3=\mu W_1+\sigma W_0,W_4=\frac{4}{3}W_2+2\sigma W_1+\frac{8}{3}EW_0.$$
(73)
Referring back to Fig. 1, we notice that the area $`S`$ enclosed by the orbit of the phase point, for some given $`\sigma `$, $`\mu `$, $`E`$, is given by

$$S(\sigma ,\mu ,E)=2_{K_{min}}^{K_{max}}K^{}𝑑K.$$
(74)
Differentiating w.r.t. $`\sigma `$ and using the fact that the integrand vanishes at the endpoints, we find

$$\mathrm{\Theta }_0=2\frac{S}{\sigma }.$$
(75)
Length and bending energy per full oscillation also follow by simple differentiation

$$L_0=\frac{S}{E},F_0=4\frac{S}{\mu }.$$
(76)
Notice that these quantities are defined whether or not the closure condition is satisfied, i.e., they are functions of the three independent variables $`\sigma `$, $`\mu `$, $`E`$. Combining (74) and (73) above, we get

$$S=\frac{4}{3}EL+\frac{1}{6}\mu F+\frac{1}{2}\sigma \mathrm{\Theta }.$$
(77)
One may regard $`L_0`$, $`\mathrm{\Theta }_0`$ and $`F_0`$ as new coordinates in the space of (not necessarily closed) configurations — the change of coordinates is non-singular, except for some special points, and is given by a Legendre transform. Then the closure condition, Eq. (11), restricts to the configurations that lie on the $`L_0`$-$`F_0`$ plane, with $`\mathrm{\Theta }_0=2\pi /n`$.