Predicting the Characteristics of Defect Transitions on Curved Surfaces

We follow the covariant formalism of Bowick et al. [21], Giomi and Bowick [24] with stress-free boundary conditions. The great accuracy of this approach for the present type of problem, and its relation to other formalisms detailed in Li et al. [37], have been discussed in Agarwal and Hilgenfeldt [31]. The elastic energy for an arbitrary disclination position ${\bf x}_{D}$ on a manifold $\mathbb{P}$ reads

$\displaystyle F_{\text{el}}(\mathbf{x}_{D})=\frac{1}{2Y}\int\Gamma^{2}(\mathbf{x},\mathbf{x}_{D})\mathrm{d}\mathbf{x},$

(2)

where the isotropic stress $\Gamma$ can be decomposed as

$\displaystyle\Gamma(\mathbf{x},\mathbf{x}_{D})=-\Gamma_{D}(\mathbf{x},\mathbf{x}_{D})-\Gamma_{s}(\mathbf{x})+U(\mathbf{x},\mathbf{x}_{D}),$

(3)

and

	$\displaystyle\Gamma_{D}(\mathbf{x},\mathbf{x}_{D})=-\frac{\pi}{3}YG_{L}(\mathbf{x},\mathbf{x}_{D}),$		(4a)
	$\displaystyle\Gamma_{S}(\mathbf{x})=Y\int K_{G}(\mathbf{y})G_{L}(\mathbf{x},\mathbf{y})\mathrm{d}\mathbf{y},$		(4b)

where

$\displaystyle G_{L}(\mathbf{x};\mathbf{x}_{D})=\frac{1}{2\pi}\log\bigg{|}\frac{z(\mathbf{x})-z(\mathbf{x}_{D})}{1-z(\mathbf{x})\overline{z(\mathbf{x}_{D})}}\bigg{|}$

(5)

is the Green’s function of the covariant Laplace operator on $\mathbb{P}$ and $z(\mathbf{x})=\varrho(r)e^{i\phi}$ is a point on the unit disk on the complex plane while $\varrho(r)$ is the conformal radius of a map from the surface onto the unit disk. Here, (4a) is the contribution due to a disclination positioned arbitrarily at ${\bf x}_{D}$ while (4b) captures the screening effect of Gaussian curvature $K_{G}(\mathbf{x})$. The harmonic term $U({\bf x},{\bf x}_{D})$ is determined by the boundary conditions at the rim of the surface; the sum of these three contributions is a fine balance of different effects. Balancing $\Gamma_{D}$ and $\Gamma_{S}$ represents a form of local curvature argument, where local Gaussian curvature compensates for a defect charge, while $U({\bf x},{\bf x}_{D})$ introduces non-trivial boundary-dependent modifications. Going beyond Giomi and Bowick [24], in the present work we perform the explicit computation of the harmonic function $U$ on the manifold $\mathbb{P}$ for arbitrary positioning of a disclination; one finds

$\displaystyle U({\bf x},{\bf x}_{D})=-Y\int\mathrm{d}\mathbf{y}H(\mathbf{x},\mathbf{y})\left(\frac{\pi}{3}\delta(\mathbf{y},\mathbf{x}_{D})-K_{G}(\mathbf{y})\right),$

(6)

where the harmonic kernel $H(\mathbf{x})$ is given by [38]:

$$H(\mathbf{x},\mathbf{y})=-\frac{1}{2\pi}f_{0}(\varrho_{y})-\sum_{n\geq 1}\frac{1}{\pi}\varrho_{x}^{n}\varrho_{y}^{n}\cos n(\phi_{y}-\phi_{x})f_{n}(\varrho_{y})$$

(7)

This formalism relies on being able to explicitly find a conformal mapping from $\mathbb{P}$ onto the unit disk where points ${\bf x},{\bf y}$ on the surface are mapped to complex numbers with moduli $\varrho_{x}$, $\varrho_{y}$ and arguments $\phi_{x},\phi_{y}$, respectively (see Supplementary Information for details and the functional forms of $\varrho,f_{n}$). The first term of (7) is a radially symmetric contribution due to the boundary conditions while the infinite sum acknowledges asymmetry introduced by, e.g., an intermediate position of the disclination. Accordingly, (6) can now be split into two terms: $U(\mathbf{x},\mathbf{x}_{D})=U_{K}({\bf x})+U_{D}(\mathbf{x},\mathbf{x}_{D})$, corresponding to the boundary contributions due to Gaussian curvature and the non-trivial disclination singularity, respectively, viz.

	$\displaystyle U_{K}(\mathbf{x})=$	$\displaystyle Y\int\mathrm{d}\mathbf{y}H(\mathbf{x},\mathbf{y})K_{G}(\mathbf{y}),$		(8a)
	$\displaystyle U_{D}(\mathbf{x},\mathbf{x}_{D})=$	$\displaystyle\frac{f_{0}(\varrho_{D})}{6}+\sum_{n\geq 1}\frac{(\varrho\varrho_{D})^{n}}{3}f_{n}(\varrho_{D})\cos n(\phi-\phi_{D}).$		(8b)

$U_{K}$ is a constant on rotationally symmetric surfaces, while $U_{D}(\mathbf{x},\mathbf{x}_{D})$, which was not taken into account in earlier work [24], has to be computed with the series truncated at a large but finite $n$; in general both computations are executed numerically. Note that the infinite sum of (8b) vanishes for a disclination at the boundary (where $f_{n\geq 1}=0$) as well as for a disclination at the center (where $\varrho_{D}=0$). For intermediate positions, as considered here, this term is generally non-zero.

Taking (2) together with (3) we have derived a rigorous covariant expression for the total energy that can be evaluated for arbitrary positions of the disclination. This reads

$\displaystyle F_{\text{el}}(\mathbf{x}_{D})=\frac{1}{2Y}\int\left[(U_{D}-\Gamma_{D})+(U_{K}-\Gamma_{S})\right]^{2}\mathrm{d}\mathbf{x}\,,$

(9)

indicating that the energy penalty can be interpreted as a combination of mismatches between topological charges and harmonic contributions from disclinations and Gaussian curvature.

The difference of $F_{\text{el}}$ for a configuration with one defect at arbitrary $\mathbf{x}_{D}$ and $F_{\text{el}}$ for a configuration with only boundary defects (defect position $\mathbf{x}_{D}=\mathbf{x}_{b}$) is more explicitly

$\displaystyle\Delta F_{\text{el}}(\mathbf{x}_{D})=\frac{1}{2Y}\int(U_{D}-\Gamma_{D})\left[(U_{D}-\Gamma_{D})+2(U_{K}-\Gamma_{S})\right]\mathrm{d}\mathbf{x},$

(10)

because only $U_{D}$ and $\Gamma_{D}$ depend on $\mathbf{x}_{D}$. Equation (10) is the generalization of the results of Agarwal and Hilgenfeldt [31] to arbitrary disclination position $\mathbf{x}_{D}$. For the portion of the present work that considers radially symmetric surfaces, we can write $\Delta F_{\text{el}}(r_{D})$. In the following, we scale out the (constant) material modulus, setting $Y=1$.

The formalism above simplifies if the defect position is the apex of an axisymmetric cap surface ($\mathbf{x}_{D}=0$), which eliminates all angular dependences and makes $U_{D}$ and $U_{K}$ straight surface averages of $\Gamma_{D}$ and $\Gamma_{S}$, respectively. For a given surface shape, a critical cap radius $r_{b}=r_{c}$ can then be defined as the extent of the surface for which $\Delta F_{\text{el}}(0)=0$ according to (10).

In Agarwal and Hilgenfeldt [31] it was shown that, rather than numerically evaluating (10), a simple criterion predicts $r_{c}$ with great accuracy, namely,

$\displaystyle\Gamma_{S}(0)\equiv\int\displaylimits_{0}^{r_{c}}K_{G}(r)\log\varrho(r)\sqrt{g}\,\mathrm{d}r=-\Gamma_{S_{0}}\,,$

(11)

where $\Gamma_{S_{0}}\equiv 1/6$ is a universal constant. Thus, the Gaussian curvature weighted with the characteristic singularity $\log\varrho(r)$ of the defect stress governs the transition in parameter space. Eq. (11) yields $r_{b}=r_{c}$ as a function of shape parameters of the surface. In deriving this analytical criterion, an approximation was pursued that uses direct leading-order Taylor expansions of $\Gamma_{D}(r)$ and $g(r)$ around $r=0$, i.e., the small-slope approximations

$\displaystyle\Gamma_{D}^{ss}(r,0)$

$\displaystyle=-\frac{1}{6}\log\left(r/r_{b}\right),\quad\sqrt{g}^{ss}=r,$

(12)

but uses a non-local approximation for $\Gamma_{S}(r)$ that matches the values of this function at $r=0$ and the boundary position $r=r_{b}$, i.e.,

$$\Gamma_{S}^{nl}(r)=\Gamma_{S}(0)\left(1-r^{2}/r_{b}^{2}\right),$$

(13)

where $\Gamma_{S}(0)$ is explicitly given by the integral in (11), rather than using the small-slope expansion

$$\Gamma_{S}^{ss}(r)=-\frac{1}{4}K_{G}(0)r_{b}^{2}\left(1-r^{2}/r_{b}^{2}\right)\,.$$

(14)

While the use of (11) allows for accurate prediction of the transition in parameter space, superior to the small-slope approximation [31], in this work we shall show that in order to predict the characteristics of the transition a further improvement in analytical theory is needed.

The example of spheroids illustrates both of these characteristics of continuous and discontinuous transitions: Fig. 2a plots the full covariant $\Delta F_{\text{el}}$ vs. defect position $r_{D}/r_{b}$ for spherical caps ($\kappa=1$), varying the cap extent $r_{b}$ through the critical value $r_{c}$. Here and in the following, we normalize energy differences by an elastic background energy scale

$\displaystyle F_{B}\equiv\pi\kappa^{4}r_{b}^{6}/384,$

(15)

whose value is obtained by inserting the small-slope (14) and $U_{K}=(1/A)\int\Gamma_{S}^{ss}dA=-(1/8)K_{G}(0)r_{b}^{2}$ into (2) for a defect free surface, i.e. $\Gamma_{D}=U_{D}=0$, and noting $K_{G}(0)=\kappa^{2}$.

The red line in Fig. 2a (for $r_{b}=r_{c}$) shows that, while $\Delta F_{\text{el}}(0)=0$, the energy is even lower for an intermediate defect position. Increasing $r_{b}$ smoothly moves this optimum disclination position from the boundary to the apex (crosses). This process begins at $r_{b}=r_{c}^{(1)}<r_{c}$ determined by $\Delta F_{\text{el}}^{\prime}(r_{D}=r_{b})=0$ (orange curve) and ends at $r_{b}=r_{c}^{(0)}>r_{c}$ determined by $\Delta F_{\text{el}}^{\prime}(r_{D}=0)=0$ (blue curve). Outside of this interval of $r_{b}$ values, $\Delta F_{\text{el}}(r_{D})$ is monotonic. As noted in Li et al. [28], spherical caps are thus an exponent of a continuous (”second-order”) disorder transition.

By contrast, Fig. 2b displays the energy landscapes for a prolate spheroid of $\kappa=5$. At $r_{b}=r_{c}$, the critical energy curve now shows a maximum, and this remains true for all $r_{b}$ in an interval $r_{c}^{(0)}<r_{b}<r_{c}^{(1)}$, where the $r_{c}^{(i)}$ are defined as above, but have switched numerical order. As a result, upon increasing $r_{b}$, an energy barrier prevents the disclination at the boundary from moving until $r_{b}>r_{c}^{(1)}>r_{c}$, when the defect abruptly jumps to the apex. This discontinuous transition thus requires a larger cap extent than $r_{c}$ or equivalently “overcharging” beyond the $q=1$ single disclination [27].

In Fig. 2c, we plot the $\kappa$-dependent values of $r_{c}^{(0)}$ (blue-dashed) and $r_{c}^{(1)}$ (orange-dashed) for the entire family of spheroidal caps. They flank the red line of the nominal transition criterion $\Delta F_{\text{el}}(0)=0$ closely and intersect it at a higher-order critical point where $r_{c}^{(0)}$ and $r_{c}^{(1)}$ switch order and, therefore, continuous transitions become discontinuous. The range of cap sizes over which the transition character manifests itself (for a given $\kappa$) is relatively small, as illustrated in the two insets of Fig. 2(c).

Can we observe such features more generally, going beyond spheroids? In Fig. 3, we take a closer look at the shape and scale of the energy landscape at critical extent $r_{b}=r_{c}$ for four different shape families, varying $\kappa=f^{\prime\prime}(0)$ for the first three and the fourth apical derivative $\lambda=f^{(4)}(0)$ (keeping $\kappa$ fixed) for the last one. Figure 3a again illustrates the transition from continuous to discontinuous in spheroids, and pinpoints the higher-order critical point at $\kappa\approx 1.3$. Hyperboloids (a shape family with potentially infinite cap extent, Fig. 3b) and a family of ”bell-shaped” surfaces whose Gaussian curvature changes sign on the surface (Fig. 3c), by contrast, never show continuous transitions. For large $\kappa$, however, the normalized energy difference $\Delta F_{\text{el}}/F_{B}$ for all shapes approaches a common energy barrier. Figure 3d shows that the transition character can be changed at constant $\kappa$ by changing the value of the fourth apical derivative $\lambda$.

This is an indication that higher apical derivatives of the surface shape play a central role here. Indeed, when solving this problem for the unique surface without higher derivatives (the paraboloid), the result accurately depicts the common asymptotic shape of the energy landscape (solid lines in Fig. 3).

The numerical integration of (10) yields these results, but gives little physical insight. When do we expect continuous vs. discontinuous transitions? Is the scale of the energy landscape (a few percent of the normalization value $F_{B}$) indeed universal, as suggested by these examples? Exactly which features of the surface are determinants of the transition character? To answer these questions, we turn to analytical approximations.

As detailed in section II.2, the accurate prediction of the transition line in Fig. 2c was not possible with a small-slope Taylor expansion, but needed the non-local improvement (13) [31]. Using either $\Gamma_{S}^{ss}$ or $\Gamma_{S}^{nl}$, and the resulting changes in $U_{K}$, while maintaining the approximations (12), $\Delta F_{\text{el}}(r_{D})$ becomes analytically tractable. Figure 4a reproduces the continuous structure of the transition for $r_{b}$ around $r_{c}$ for spherical caps according to the rigorous numerical computation. Figures 4b and c show the analogous results for the small-slope and non-local approximations, respectively. Intriguingly, neither of these approaches produces any secondary characteristics at all – the energy difference at transition is flat, and monotonic at every other value of $r_{b}$. Thus, even though the non-local approximation reproduces the spread of energy values more accurately than the small-slope model, it does not probe the shape properties of the surface that are responsible for the order of the transition.

We now describe a minimal model that captures the dominant effects of the secondary transition character. Noting that a more accurate representation of the intrinsic isotropic stress $\Gamma_{S}(r)$ was crucial for predicting the transition location, we improve further on the approximation (13): in addition to matching the function value and first derivative at $r=0$ and function value at $r=r_{b}$, we now require a match of the second derivative at the apex. By symmetry, this requires a fourth-order polynomial in $r$, namely,

	$\displaystyle\Gamma_{S}^{nl}(r)=$	$\displaystyle-\frac{\Gamma_{S}^{\prime\prime}(0)r_{b}^{2}}{2}\left(1-\frac{r^{2}}{r_{b}^{2}}\right)$
		$\displaystyle+\left(\Gamma_{S}(0)+\frac{\Gamma_{S}^{\prime\prime}(0)r_{b}^{2}}{2}\right)\left(1-\frac{r^{4}}{r_{b}^{4}}\right),$		(16)

with the full rigorous expressions

$\displaystyle\Gamma_{S}(0)\equiv$

$\displaystyle\int\displaylimits_{0}^{r_{b}}K_{G}(r)\log\varrho(r)\sqrt{g}\,\mathrm{d}r,\,\Gamma_{S}^{\prime\prime}(0)=\frac{K_{G}(0)}{2}=\frac{\kappa^{2}}{2}\,.$

(17)

We use (4a) and (5) for arbitrary defect position $\mathbf{x}_{D}$ but employ the small-slope $\varrho^{ss}(r)=r/r_{b}$ and $\sqrt{g}^{ss}=r$. Inserting these into (10) and (8), the energy integral can be executed analytically (see Supplementary Information for details) and results in the following closed form expression:

	$\displaystyle\Delta F_{\text{el}}(r_{D})=\frac{\pi r_{b}^{2}}{864}\left(1-\frac{r_{D}^{2}}{r_{b}^{2}}\right)^{2}\bigg{[}$	$\displaystyle 3+8\left(2+\frac{r_{D}^{2}}{r_{b}^{2}}\right)\Gamma_{S}(0)$
		$\displaystyle-\left(1-4\frac{r_{D}^{2}}{r_{b}^{2}}\right)\frac{\kappa^{2}r_{b}^{2}}{2}\bigg{]}.$		(18)

The transition threshold is defined as the extent $r_{b}$ at which $\Delta F_{\text{el}}(0)=0$ and (18) results in

$\displaystyle\Gamma_{S}(0)=-3/16+\kappa^{2}r_{c}^{2}/32$

(19)

as the transition criterion within this level of approximation. The extremum of the energy landscape is obtained from

	$\displaystyle\Delta F_{\text{el}}^{\prime}(r_{D})=\frac{\pi r_{D}}{72}\left(\frac{r_{D}^{2}}{r_{b}^{2}}-1\right)$	$\displaystyle\bigg{[}1+4\left(1+\frac{r_{D}^{2}}{r_{b}^{2}}\right)\Gamma_{S}(0)$
		$\displaystyle-\left(1-2\frac{r_{D}^{2}}{r_{b}^{2}}\right)\frac{\kappa^{2}}{2}r_{b}^{2}\bigg{]}$		(20)

by setting $\Delta F_{\text{el}}^{\prime}(r_{D})|_{r_{b}=r_{c}}=0$. Inserting the transition criterion (19) leads to a straightforward non-trivial solution for the position of the minimum or maximum of the energy landscape $r_{D,m}=r_{c}/\sqrt{3}$. This result is indicated by dashed vertical lines in Fig. 3 and in good agreement for the vast majority of shapes. Inserting $r_{D,m}$ together with (19) into (18) yields the scale of the secondary structure. Normalizing by $F_{B}$, one obtains

$\displaystyle\Delta F_{\text{el}}^{\text{sec}}/F_{B}=\frac{4}{81\kappa^{4}r_{c}^{4}}\left(-2+3\kappa^{2}r_{c}^{2}\right)$

(21)

This rigorously shows why a small-slope formalism is incapable of predicting secondary transition structure (since $\kappa^{2}r_{c}^{2}=2/3$ in this approximation [27]). Further analytical progress can be made by determining $r_{c}(\kappa)$ at transition within the present approximation. By symmetry, the large $\kappa$ expansion of $r_{c}$ for rotationally symmetric surfaces is $r_{c}=a/\kappa+b/\kappa^{3}+\dots$ (cf. Agarwal and Hilgenfeldt [31]). Inserting this into (19) and Taylor-expanding both sides for $\kappa\to\infty$ allows us to solve for the coefficients $a$, $b$, $\dots,$ order-wise. With the definition of $\Gamma_{S}(0)$, this reads

	$\displaystyle\Gamma_{S}(0)=$	$\displaystyle 1-\sqrt{1+a^{2}}+\log\left[\frac{\left(1+\sqrt{1+a^{2}}\right)}{2}\right]+\frac{b(1-\sqrt{a^{2}+1})}{a\kappa^{2}}$
		$\displaystyle+\frac{\left(-a^{4}+a^{2}-2\sqrt{a^{2}+1}+2\right)(\lambda/\kappa)}{18\sqrt{a^{2}+1}\kappa^{2}}+\dots$		(22)

Finally, expanding the RHS of (19) and equating with (22) we obtain $a\approx 0.844,\,b\approx-0.042\lambda/\kappa$. Note that a generic surface (with one dominant radial length scale) will have $\lambda=\mathcal{O}(\kappa)$ and $b=\mathcal{O}(1)$. Inserting into a large $\kappa$ expansion of (21), we obtain:

	$\displaystyle\Delta F_{\text{el}}^{sec}/F_{B}$	$\displaystyle=\frac{4(-2+3a^{2})}{81a^{4}}+\frac{8(4-3a^{2})b}{81a^{5}\kappa^{2}}+\dots$
		$\displaystyle\approx 0.0135-0.018\frac{\lambda}{\kappa^{3}}+\dots$		(23)

This suggests a universal scale for the energy barrier at transition for surfaces with large central curvature, consistent with the observations of Fig. 3. For $\kappa\to\infty$, higher derivatives are negligible and the behavior mimics that of a paraboloid. Furthermore, the result predicts a change from an energy maximum to a minimum at transition (a higher-order critical point), only if the quantity $\lambda/\kappa$ is positive. Note that for spheroids $\lambda/\kappa=3>0$, whereas for the hyperboloids or “bell” shapes of Fig. 3bc $\lambda/\kappa=-3<0$. Likewise, increasing $\lambda$ from zero should lead to the development of an energy minimum and thus a continuous transition, as shown in Fig. 3d.

The formalism outlined here also settles the question of the range of the secondary structure effects described in Sec. III.1 i.e., the values of $r_{c}^{(0)}$ and $r_{c}^{(1)}$. Using the defining criteria $\Delta F^{\prime}_{\text{el}}(r_{D}=0)|_{r_{b}=r_{c}^{(0)}}=0$ and $\Delta F^{\prime}_{\text{el}}(r_{D}=r_{b})|_{r_{b}=r_{c}^{(1)}}=0$ in (20) results in the following implicit equations:

$\displaystyle\Gamma_{S}(0)=-\frac{1}{8}+\frac{(\kappa r_{c}^{(0)})^{2}}{16},\quad\Gamma_{S}(0)=-\frac{1}{4}+\frac{(\kappa r_{c}^{(1)})^{2}}{8}\,,$

(24)

respectively. As before, we insert the large $\kappa$ expansion $r_{c}^{(i)}=a/\kappa+b/\kappa^{3}+\dots$, Taylor expand and determine the respective coefficients $a,b,\dots$ to finally obtain

$\displaystyle(r_{c}^{(1)}-r_{c}^{(0)})/r_{c}\approx 0.038-0.07\lambda/\kappa^{3}$

(25)

This suggests that the range of the secondary features is about $4\%$ relative to the critical $r_{c}$ for surfaces with large central curvature – again, the large $\kappa$ limit universally asymptotes to that of a paraboloid. On the other hand, shapes with small $\kappa$ could have larger ranges, but are not accurately described by this approach. Empirically, we see for all shape families analyzed that the range of secondary characteristics remains of this order up to the smallest $\kappa$ consistent with unique $f(r)$ shapes.

All of these predictions are in complete qualitative agreement with the numerical computations, while quantitative discrepancies in $\Delta F_{\text{el}}^{\text{sec}}$ can be systematically improved upon employing higher order approximations of the type that also improve the predictive error of $r_{c}(\kappa)$ [31]. Matching the third derivative of the metric $\sqrt{g}$ at $r=0$ yields

$\displaystyle\sqrt{g}=r+\frac{\kappa^{2}}{2}r^{3}.$

(26)

With only this modification ($\Gamma_{D}$ remains unchanged from (12) to permit analytical evaluation of the integrals) we go through the same steps outlined above and, while the explicit expressions are more complicated, we obtain a slightly different transition criterion $r_{c}\approx 0.829/\kappa-0.039(\lambda/\kappa)/\kappa^{3}$ and an altered secondary energy,

$\displaystyle\Delta F_{\text{el}}^{\text{sec}}/F_{B}\approx 0.052-0.017\frac{\lambda}{\kappa^{3}}+\dots,$

(27)

which more closely approximates the empirical universal barrier height for $\kappa\to\infty$, which is $\Delta F_{\text{el}}^{\text{sec}}/F_{B}\approx 0.034$ in the full covariant computation. Eq. (27) predicts that the higher-order critical point is given by $\lambda_{h}\approx 3\kappa_{h}^{3}$. For spheroids, this locates the transition from continuous to discontinuous behavior at $\kappa_{h}\approx 1$, again in good agreement with the covariant computation ($\kappa_{h}\approx 1.3$). The remaining quantitative discrepancies can be systematically alleviated by including higher order terms in the non-local expansions, although not all integrals may be tractable analytically.

The results from the previous subsection would suggest that only local information about the surface (such as $\kappa$ or $\lambda$) is sufficient to predict the nature of the transition, even though we employed a non-local formalism with the full $\Gamma_{S}(0)$. This raises the question whether a local higher-order small-slope formalism would yield similar quantitative agreement with the full covariant formalism. In order to test this, we instead employ a local fourth-order expansion of $\Gamma_{S}(r)$ that reads:

$\displaystyle\Gamma_{S}^{l}(r)=-\frac{\kappa^{2}r_{b}^{2}}{4}\left(1-\frac{r^{2}}{r_{b}^{2}}\right)+\frac{\kappa^{4}r_{b}^{4}}{96}\left(3-4\frac{\lambda}{\kappa^{3}}\right)\left(1-\frac{r^{4}}{r_{b}^{4}}\right),$

(28)

where we now use the small-slope expansion for $\Gamma_{S}(0)=-\kappa^{2}r_{b}^{2}/4+(3-4\lambda/\kappa)\kappa^{4}r_{b}^{4}/96+\dots$, retaining terms up to $\mathcal{O}(r^{4})$. Together with (26) ($\Gamma_{D}$ remains unmodified from (12)), we go through the same steps as outlined in the previous subsection, i.e. we evaluate the energy integral analytically, determine the transition criterion at this level of approximation, compute the position of the extremum and finally insert all these into $\Delta F_{\text{el}}/F_{B}$ to obtain the scale of the normalized secondary structure, which reads

$\displaystyle\Delta F_{\text{el}}^{\text{sec}}/F_{B}\approx 0.057-0.021\frac{\lambda}{\kappa^{3}}+\dots.$

(29)

We plot in Fig. 5(a) the full covariant $\Delta F_{\text{el}}^{\text{sec}}/F_{B}$ obtained by numerically integrating (10) (indicated by red dots) and compare the local approximation with the non-local one for the shape family of spheroids as our prototypical example. The large central curvature ($\kappa\to\infty$) limit for both converges to approximately the same asymptote with similar deviations from the covariant prediction (cf. (29) and (27)), however differences between the two become significant for surfaces with smaller values of $\kappa$ that are arguably of more practical relevance [31]. Fig. 5(b) exemplifies the significant discrepancy for the particular case of $\kappa=1$ (sphere) – while the local approximation qualitatively captures the behavior at transition it greatly underestimates the magnitude of the energy minimum and the disagreement only gets worse for surfaces with $\kappa<1$. Therefore, analogous to the transition criterion from our previous work [31], quantitatively describing the scale of the secondary structure around transition requires non-local information from the entire surface for shapes with $\kappa\sim\mathcal{O}(1)$.

We now turn our attention to the question of how robust this secondary-structure behavior is. In particular, the small range of cap sizes over which it occurs suggests it may be qualitatively altered by even slight deviations from the radial symmetry of the surface shape.

In general, the boundary of open surfaces will not obey radial symmetry; the Drosophila eye, for example, is well approximated by an ellipsoidal cap with an elliptical boundary [39]. While a fully covariant formalism in such a general scenario might be feasible numerically, we are not aware of any previous work that studies this role of rotational symmetry breaking. In order to estimate the magnitude of the effect of breaking boundary shape symmetry on the energy landscape, we shall utilize the small-slope approximation; as we have shown above, this does not, by itself, lead to any secondary structure of the transition (cf. Fig. 4b). As we have furthermore shown that the energy scale of intrinsic secondary structures is bounded at least for large central curvature $\kappa$, we can give a quantitative estimate of whether symmetry breaking will change this structure.

We apply a leading-order shape perturbation to a general surface of revolution, generating an ellipse from the circular boundary by stretching/contracting perpendicular axes by an amount $\epsilon$. Thus,

$\displaystyle x=(1+\epsilon)r\cos\phi,\quad y=(1-\epsilon)r\sin\phi,\quad z=f(r)$

(30)

with $\epsilon\ll 1$, while the metric tensor in the small-slope limit is

$\displaystyle g_{ij}=\begin{bmatrix}1+\epsilon(2+\epsilon)\cos 2\phi&-2r\epsilon\sin 2\phi\\ -2r\epsilon\sin 2\phi&r^{2}\left(1+\epsilon(\epsilon-2)\cos 2\phi\right)\end{bmatrix}.$

(31)

Here, $\sqrt{g}=r(1-\epsilon^{2})$ and eccentricity $e=\sqrt{1-\frac{(1-\epsilon)^{2}}{(1+\epsilon)^{2}}}=2\sqrt{\epsilon}+\mathcal{O}(\epsilon^{3/2})$, while the small-slope Gaussian curvature is constant to $\mathcal{O}(\epsilon)$ and is given by

$\displaystyle K_{G}(r,\phi)=\kappa^{2}+\mathcal{O}\left(\epsilon^{2}\right)\,.$

(32)

The elliptical boundary of a section cut parallel to the $xy$-plane is given by

$\displaystyle r(\phi)=\frac{r_{b}(1-\epsilon^{2})}{\sqrt{1+\epsilon(\epsilon-2)\cos 2\phi}}\approx r_{b}(1+\epsilon\cos 2\phi).$

(33)

Breaking rotational symmetry of the surface leads to coupling of stresses in the radial and azimuthal directions. Therefore, the isotropic stresses are not as easily computed as in the previous subsections. Instead, we start by observing that the Airy stress function $\chi$ satisfies the following equation in polar coordinates [27]:

$\displaystyle\nabla^{4}_{\perp}\chi$

$\displaystyle=\frac{\pi}{3}\frac{\delta(r-r_{D})\delta(\phi-\phi_{D})}{\sqrt{g}}-K_{G}(r,\phi)$

(34)

subject to a zero normal stress boundary condition and regularity conditions at $r=0$. Here, $\nabla^{2}_{\perp}$ is the 2D Laplace-Beltrami operator in the small-slope limit and is given by $\nabla^{2}_{\perp}f=\frac{1}{\sqrt{g}}\partial_{i}\left(\sqrt{|g|}g^{ij}\partial_{j}f\right)$ using the small-slope limit (31). Here, $g^{ij}$ is the inverse of the metric tensor such that $g^{ij}g_{jk}=\delta^{i}_{k}$. The Airy function is related to the stress components in the usual way,

$$\sigma_{rr}=\frac{1}{r}\frac{\partial\chi}{\partial r}+\frac{1}{r^{2}}\frac{\partial^{2}\chi}{\partial\phi^{2}},\,\sigma_{r\phi}=-\frac{\partial}{\partial r}\left(\frac{1}{r}\frac{\partial\chi}{\partial\phi}\right),\,\sigma_{\phi\phi}=\frac{\partial^{2}\chi}{\partial r^{2}}.$$

(35)

Analogous to the case of a circular boundary, we impose vanishing normal stress, i.e. $\bm{\sigma}\cdot\hat{\bm{n}}=0$ on (34), where $\hat{\bm{n}}=n_{r}\hat{\bm{e}}_{r}+n_{\phi}\hat{\bm{e}}_{\phi}$ is the normal vector to the elliptical boundary. Therefore, one obtains two scalar equations:

$\displaystyle\bigg{[}\left(n_{r}\sigma_{rr}+n_{\phi}\sigma_{r\phi}\right)\hat{\bm{e}}_{r}+\left(n_{r}\sigma_{r\phi}+n_{\phi}\sigma_{\phi\phi}\right)\hat{\bm{e}}_{\phi}\bigg{]}_{r=r_{b}(\phi)}=0\,.$

(36)

The total in-plane elastic energy for a surface with stress-free boundary and metric $g$ in terms of the stress components (assuming a linear constitutive relation) is given by [24, 27]:

$\displaystyle F_{\text{el}}=\frac{1}{2}\int\Gamma(r,\phi)^{2}dA=\frac{1}{2}\int\left(\sigma_{rr}+\sigma_{\phi\phi}\right)^{2}dA,$

(37)

where $\Gamma=\sigma_{rr}+\sigma_{\phi\phi}$ is the trace of the stress tensor.

Writing $\chi$, and thus $\sigma_{ij}$, as a Fourier series expansion in $(\cos n\phi,\sin n\phi)$ and consistently expanding the stress components as well as $\nabla^{2}_{\perp}$ and the boundary condition in powers of $\epsilon$, we evaluate the elastic energy up to $\mathcal{O}(\epsilon)$ (see Supplemental Information for details). Inserting the expansion $\Gamma=\Gamma_{0}+\epsilon\Gamma_{1}$ into (37), the elastic energy can be cast explicitly as

$\displaystyle F_{\text{el}}=\frac{1}{2}\int\limits_{0}^{2\pi}\int\limits_{0}^{r_{b}}\left(\Gamma_{0}^{2}+2\epsilon\Gamma_{0}\Gamma_{1}\right)rdrd\phi\,+\mathcal{O}(\epsilon^{2}).$

(38)

In terms of Fourier components, it is evident that only the squares of the modes will contribute to the energy while the cross terms will integrate out to zero.

Systematically evaluating the $\mathcal{O}(1)$ and $\mathcal{O}(\epsilon)$ contributions to the stress tensor, we obtain an angular correction to $\Delta F_{\text{el}}$. The final expression reads

	$\displaystyle\frac{\Delta F_{\text{el}}(\overline{r}_{D})}{F_{B}}$	$\displaystyle\approx\frac{2\left(2-3\kappa^{2}r_{b}^{2}\right)}{3\kappa^{4}r_{b}^{4}}\left(1-\overline{r}_{D}^{2}\right)^{2}$
		$\displaystyle+\frac{4}{3}\epsilon\overline{r}_{D}^{2}\cos 2\phi_{D}\left(\frac{-1+2\overline{r}_{D}^{2}-3\overline{r}_{D}^{4}+2\overline{r}_{D}^{6}-4\log\overline{r}_{D}}{\kappa^{4}r_{b}^{4}}\right),$		(39)

where $\overline{r}_{D}=r_{D}/r_{b}$. The leading order is the small-slope energy expression for a symmetric cap obtained, e.g., by Azadi and Grason [27], while the next term represents the leading effect of the shape perturbation and depends, in addition to $r_{D}$, on the angular position $\phi_{D}$ of the defect.

We display in Fig. 6, Eq. (39) as a function of $\overline{r}_{D}$ for $\epsilon=0.05$ (eccentricity of $e\approx 0.45$) along (a) the major axis ($\phi_{D}=0$) and (b) the minor axis ($\phi_{D}=\pi/2$). As expected, breaking the circular cross-sectional symmetry forces the system to select a preferred direction of migration of the defect, which in this situation is the minor axis: this perturbed small-slope formalism predicts a continuous variation of the defect position along $\phi=\pi/2$ until the apex position is established for a certain energy difference (here, $\Delta F_{\text{el}}/F_{B}\approx-0.35$ as marked by the green curve in Fig. 6(b)).

To quantify the scale of the secondary structure, we set $\Delta F_{\text{el}}(0)=0$ in $\eqref{DelFel_ellipse}$ to obtain the transition criterion $r_{c}=\sqrt{2/3}/\kappa$ (unchanged from the rotationally symmetric situation), while the position of the maxima/minima obtained by setting $\Delta F_{\text{el}}^{\prime}(0)=0$ results in $r_{D}\approx 0.55r_{b}$, a slight but significant difference from the position for intrinsic secondary structures described in section III.2. Finally, inserting into $\eqref{DelFel_ellipse}$ we obtain the energy scale

$\displaystyle\frac{\Delta F^{\text{sec}}_{\text{el}}}{F_{B}}\approx 1.61\epsilon\cos 2\phi_{D}+\mathcal{O}\left(\epsilon^{2}\right)\,.$

(40)

At the critical cap extent $r_{b}=r_{c}$, the small-slope $\Delta F_{\text{el}}$ vanishes, so that the secondary energy structure is proportional to $\epsilon$ and quadrupolar in $\phi_{D}$. Within the small-slope approximation this correction does not depend on $\kappa$ when normalized by $F_{B}$. The magnitude of secondary energy maxima and minima (cf. Fig. 6) is well described by (40). Having identified a universal scale of intrinsic secondary structure maxima in section III.2 as $\Delta F_{\text{el}}^{\text{sec}}/F_{B}\approx 0.034$ for moderate to large $\kappa$, we can say that symmetry breaking is likely to transform the transition character from discontinuous to continuous along the minor axis beyond a strain of $\epsilon\gg 0.02$. Quite subtle symmetry breaking is thus capable of qualitatively changing the disorder transition.

We note that inclusion of the next order $\mathcal{O}(r_{b}^{4})$ terms introduces non-trivial couplings between the shape perturbation and the secondary structure discussed in the previous subsections, leading to modification of (40). We know from Sec. III.2 that the leading term of the transition criterion $r_{c}$ becomes approximately $a=0.844$. Defining the deviation of $a$ from the small-slope value as $\delta=a-\sqrt{2/3}$, a Taylor expansion of (39) in $\delta\ll 1$ alters the numerical prefactor of (40) to $(1.61-7.9\delta)\approx 1.39$. This is a small quantitative deviation that does not change the nature of the effects discussed above.