2.1. The vertical branch

Solving for equilibria of (Equation 1.4), we find that $u\,\nabla\left(u -{\mu} (V*u)\right)=\nu\in\mathbb{R}^2$, a constant. If $u(x,y)=0$ at some point, then $\nu=0$. Otherwise, $\partial_x\left(u -{\mu} (V*u)\right)=\nu_1/u$ has a sign, which contradicts periodicity, unless $\nu_1=0$. With the same reasoning for the $y$-derivative, we conclude that $\nu=0$. For a solution without vacuum, that is, for $u \gt 0$, we then find

(2.1)\begin{equation} u - {\mu}(V * u) = \rho, \end{equation}

for $\rho$ a constant. Given our simple choices of repulsive potentials, one then has the explicit vertical branches

(2.2)\begin{align} \!\!\!\!\!\!\!\,\,\,\,\,\quad\qquad\quad{\mu}_* = \frac{1}{\pi^2}, \quad u = &\,\,A_0 + A_1\cos(x)\cos(y) + B_1\sin(x)\sin(y)\nonumber \\ &+ a_1\cos(x)\sin(y) + b_1\sin(x)\cos(y) \quad &&\quad\quad\quad\quad\quad\quad\textrm{(sq.)}, \end{align}

(2.3)\begin{align} \qquad\quad{\mu}_* = \frac{1}{\sqrt{3}\pi^2}, \quad u = &\,\,A_0 + A_1\cos\left(x+\frac{y}{\sqrt{3}}\right) + B_1\cos\left(x-\frac{y}{\sqrt{3}}\right) + C_1\cos\left(\frac{2y}{\sqrt{3}}\right)\nonumber\\ &+a_1\sin\left(x+\frac{y}{\sqrt{3}}\right) + b_1\sin\left(x-\frac{y}{\sqrt{3}}\right) + c_1\sin\left(\frac{2y}{\sqrt{3}}\right) \quad &&\textrm{(hex.)}. \end{align}

We shall focus on solutions with $a_1=b_1=0$ in (Equation 2.2) and $a_1=b_1=c_1=0$ in (Equation 2.3) in the remainder of this paper. These subspaces consist of functions invariant under point reflections $(x,y)\mapsto (-x,-y)$. We only comment briefly on solutions outside of this fixed-point subspace. Solutions with maximal isotropy all lie within these subspaces, possibly after a spatial translation; see Remarks 1 and 4, below.

Fixing the average density to, say, 1, we may set $A_0=1$ and are left with two- and three-dimensional subspaces of solutions, respectively. These subspaces contain one-dimensional subspaces with larger symmetry (maximal isotropy) as follows. In the square symmetry case, $A_1 = B_1$ corresponds to fissures, solutions depending on $x-y$, only, and $B_1 = 0$ (or $A_1=0$) corresponds to bubble solutions invariant under the fourfold rotation $R$. Solutions are strictly positive, possessing no vacuum regions, as long as $\max(|A_1|,|B_1|) \lt A_0$.

In the hexagonal symmetry case, $A_1 = B_1 = C_1$ corresponds to bubble solutions with sixfold rotational symmetry $R$, while fissures correspond to solutions where only one of $A_1, B_1 ,C_1$ is nonzero. Bubble solutions are strictly positive if $-\frac{1}{3}A_0 \lt A_1,B_1, C_1 \lt \frac{2}{3}A_0$, fissure solutions for $|A_1| \lt A_0=1$, $B_1=C_1=0$.

2.2. Vacuum formation: square lattice symmetry

We analyse solutions in the vicinity of the boundary of the vertical branch, when $\max(|A_1|,|B_1|) = A_0$. We look for non-negative solutions with support $\Omega_0\subset\Omega$ that have average density $1$. They therefore satisfy

(2.4)\begin{equation} u - {\mu}(V * u) - \rho = 0, \text{for } x\in\Omega_0 \qquad\textrm{(McKean--Vlasov condition),} \end{equation}

(2.5)\begin{equation} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\frac{1}{\text{Area}(\Omega)} \iint_{\Omega_0} u(x,y) dydx = 1\!\qquad\qquad\qquad\quad\,\,\,\,\textrm{(average density condition),} \end{equation}

(2.6)\begin{equation}u(x,y)=0\;\;\text{on}\;\partial\Omega_0\qquad\qquad\quad\text{(boundary condition). }\end{equation}

The continuity condition (Equation 2.6) follows from the fact that $u$ is a weak solution.

Theorem 1

(Vacuum region scaling – square lattice). Let $V(x,y) = 2\cos(x)\cos(y)$, and set ${\mu} = \frac{1}{\pi^2} + {\widetilde{\mu}}$, with ${\widetilde{\mu}}$ sufficiently small. Then solutions $u$ to (Equation 2.4)–(Equation 2.6), invariant under point reflections $u(x,y)=u(-x,-y)$, with vacuum regions, that is, $\Omega_0\neq \Omega$, are of the form

(2.7)\begin{equation} u(x,y) = \left( A_0 + A_1 \cos (x) \cos (y)+B_1 \sin(x)\sin(y)\right)_+, \end{equation}

where $f_+=\max(f,0)$. There are then two solution branches, both supercritical, locally unique up to translations in $x$ and $y$, parameterized by ${\widetilde{\mu}} \gt rsim 0$ as $A_0=A_0({\widetilde{\mu}}),\, A_1 = A_1({\widetilde{\mu}}),\, B_1 = B_1({\widetilde{\mu}})$:

1. (1) fissures: $A_1=-B_1$,

   \begin{align*} A_0 =1-\frac{\pi}{2}{\widetilde{\mu}}+\mathcal{O}({\widetilde{\mu}}^2),\qquad A_1 =1+\left(\frac{3}{4\sqrt{2}}\pi^3{\widetilde{\mu}}\right)^{2/3}+\mathcal{O}({\widetilde{\mu}}). \end{align*}

   The half-width $\ell$ of the fissures and area $\mathcal{A}$ of vacuum region inside $\Omega$ are

   \begin{equation*} \ell=\frac{\pi}{2^{1/2}}\left(\frac{3}{2} {\widetilde{\mu}}\right)^{1/3}+\mathcal{O}({\widetilde{\mu}}^{2/3}) ,\qquad \mathcal{A}=4\pi^2\left(\frac{3}{2} {\widetilde{\mu}}\right)^{1/3}+\mathcal{O}({\widetilde{\mu}}^{2/3}). \end{equation*}

2. (2) bubbles: $B_1=0$,

   \begin{align*} A_0 = 1 - \frac{\pi^2}{4} {\widetilde{\mu}} + \mathcal{O}({\widetilde{\mu}}^{3/2}),\qquad A_1 = 1 + \frac{\pi^{3/2}}{\sqrt{2}} {\widetilde{\mu}}^{1/2} - \frac{\pi^2}{4} {\widetilde{\mu}} + \mathcal{O}({\widetilde{\mu}}^{3/2}) . \end{align*}

   The inner radius $\ell$ and area $\mathcal{A}$ of the vacuum region inside $\Omega$ are

   \begin{equation*} \ell = (2\pi^3 {\widetilde{\mu}})^{1/4} + \mathcal{O}({\widetilde{\mu}}^{3/4} ),\qquad \mathcal{A}=(2\pi^5 {\widetilde{\mu}})^{1/2}+\mathcal{O}({\widetilde{\mu}}). \end{equation*}

Proof. First, notice that (Equation 2.4) implies that for all $(x,y)$ where $u(x,y) \gt 0$, we have

\begin{equation*} u(x,y)=\rho+2{\mu}\iint_{\Omega_0} \cos(x-\xi)\cos(y-\eta)u(\xi,\eta) d\xi d\eta, \end{equation*}

which, using trigonometric addition theorems gives the form (Equation 2.7). Therefore, all solutions with (or without) vacuum are found as solutions to a finite system of equations for $A_0,A_1$ and $B_1$. We find solutions using a perturbative approach, starting with bubbles, $B_1=0$, $A_1\sim 1$ and $A_0\sim 1$. Substituting the ansatz $\left(A_0 + A_1 \cos (x) \cos (y)\right)_+$ into (Equation 2.4)–(Equation 2.6), we find

(2.8)\begin{equation} A_0 + A_1 \cos (x) \cos (y) - \left(\frac{1}{\pi^2}+{\widetilde{\mu}}\right) \iint_{\Omega_0} 2\cos(x-\xi)\cos(y-\eta) \left(A_0 + A_1 \cos (\xi) \cos (\eta)\right)d\xi d\eta - \rho = 0 \end{equation}

(2.9)\begin{equation}\frac{1}{2\pi^2} \iint_{\Omega_0} \left(A_0 + A_1 \cos (x) \cos (y) \right) dydx = 1 \end{equation}

(2.10)\begin{equation} \qquad \ A_0 + A_1 \cos (x) \cos (y) = 0 \text{on}\ \partial \Omega_0, \end{equation}

where

\begin{equation*} \Omega_0 = \left\{(x,y) \in \Omega \ \big| \ \ A_0 + A_1 \cos (x) \cos (y) \gt 0 \ \right\}; \end{equation*}

see Figure 2. Assuming that there are vacuum regions, we have $A_1 \gt A_0$, so that we may scale

\begin{align*} A_0 &= 1+a_0, \qquad A_1 = 1+ a_0 + z_1^2, \qquad a_0 = a_0(z_1), \end{align*}

for a small parameter $z_1$. Our explicit ansatz then allows us to write the free-boundary condition (Equation 2.10) in polar coordinates about the centre $(\pi, 0)$ of the bubble, as follows:

(2.11)\begin{equation} \begin{aligned} 0&= A_0 - A_1 \left ( 1- \frac{r^2}{2} + \frac{r^4(3-\cos(4\theta))}{48} + \mathcal{O}(r^6) \right) \\ &=\left( 1+a_0 \right)- \left( 1 + a_0 + z_1^2 \right) \left( 1-\frac{r^2}{2} + r^4 Q(\theta) + \mathcal{O}(r^6) \right), \end{aligned} \end{equation}

where $Q(\theta) := \frac{3-\cos(4\theta)}{48}.$ The inner radius $\ell$ of the bubble can be identified at leading order as the value of $r$ that solves (Equation 2.11). We may now further scale $\ell= \ell_1z_1$, which yields, after some algebraic manipulation,

(2.12)\begin{equation} 0 = -\frac{1}{1+a_0+z_1^2} + \frac{\ell_1^2}{2} - \ell_1^4 z_1^2 Q(\theta) + \mathcal{O}(z_1^4) =: F(\ell_1,z_1, a_0, \theta). \end{equation}

Figure 2. An illustration of vacuum bubbles (left) and fissures (right) in a square of width $2\pi$, comprising two (!) fundamental domains $\Omega$ of the square lattice. Blue regions represent areas of vacuum; $u$ is positive in the green regions $\Omega_0$. Note that there is one bubble and one fissure in each fundamental domain (boundaries of fundamental domain indicated by dashed lines); fissures have width $2\ell$ and total length $\sqrt{2}\pi$ in a fundamental domain, bubbles have radius $\ell$.

We find by the implicit function theorem that there exists a unique $\ell_1 = \ell_1(z_1,a_0,\theta)$ for $a_0, z_1$ sufficiently small, with $\ell_1(0,0,\theta) = \sqrt{2}$. Furthermore, expanding in $z_1,a_0$, we calculate

(2.13)\begin{equation} \ell_1(z_1,a_0,\theta) = \sqrt{2} - \frac{1+4Q(\theta) }{2} z_1^2 - \frac{1}{\sqrt{2} } a_0 + \mathcal{O}( z_1^4 + a_0^2). \end{equation}

We can now use this expression for the boundary of the vacuum bubble to solve the mass constraint (Equation 2.9) for $a_0$. Substituting the expression from (Equation 2.13) into (Equation 2.9), we obtain

\begin{align*} 0 & = 2\pi^2 - \iint_{\Omega} \left(A_0 + A_1 \cos (x) \cos (y) \right) dydx + \iint_{\Omega_0^c} \left(A_0 + A_1 \cos (x) \cos (y) \right) dydx \\ &= 2\pi^2-2\pi^2(1 + a_0) + \iint_{\Omega_0^c} \left( 1+a_0 + \left( 1+ a_0 + z_1^2 \right)\cos (x) \cos (y) \right)dy dx \\ & = -2\pi^2a_0 + \int_{0}^{2\pi} \int_0^{z_1\ell_1(z_1,a_0)} \left(- z_1^2 + \frac{r^2}{2} + \mathcal{O}(a_0r^2 + r^4 + r^2 z_1^2) \right) r \ drd\theta \\ &= -2\pi^2 a_0 + \pi \left( -\frac{1}{2} z_1^4 \right) + \mathcal{O} (a_0z_1^4 + z_1^6). \end{align*}

We can again solve for $a_0 = a_0(z_1)$ near $a_0 = z_1 = 0$ by the implicit function theorem, obtaining

\begin{equation*} a_0(z_1) = -\frac{1}{2\pi} z_1^4 + \mathcal{O}(z_1^6). \end{equation*}

Finally, we use (Equation 2.8) to relate these quantities to the parameter ${\widetilde{\mu}}$. Equating the coefficients of $\cos(x)\cos(y)$ in (Equation 2.8), we find

\begin{align*} 0 = &\,\, A_1 - \left(\frac{1}{\pi^2} + {\widetilde{\mu}}\right)(A_1\pi^2) + \left(\frac{1}{\pi^2} + {\widetilde{\mu}}\right)\iint_{\Omega_0^c} 2\cos(\xi) \cos(\eta)\left( A_0 + A_1 \cos(\xi) \cos(\eta) \right) d \eta d \xi \\ = & -{\widetilde{\mu}} A_1\pi^2 \\ &+ 2 \left(\frac{1}{\pi^2} + {\widetilde{\mu}}\right)\int_0^{2\pi} \int_0^{z_1\ell_1(a_0,z_1,\theta)} \left( -1+\frac{r^2}{2} + \mathcal{O}(r^4)\right) \left( A_0 + A_1 \left( -1+\frac{r^2}{2} + \mathcal{O}(r^4)\right) \right) r d r d \theta \\ = & -\pi^2 {\widetilde{\mu}} \left( 1+ z_1^2 - \frac{1}{2\pi} z_1^4 \right) + \left( \frac{1}{\pi^2} + {\widetilde{\mu}} \right) \left( 2\pi z_1^4 \right) + \mathcal{O}(z_1^6). \end{align*}

We thus see that ${\widetilde{\mu}}$ scales like $z_1^4$. We solve for $z_1$ as a function of ${\widetilde{\mu}}^{1/4}$ near $z_1 = {\widetilde{\mu}}^{1/4} = 0$ by the implicit function theorem, obtaining

\begin{equation*} z_1 = \left(\frac{\pi^3}{2}\right)^{1/4}{\widetilde{\mu}}^{1/4} + \mathcal{O}({\widetilde{\mu}}^{3/4}). \end{equation*}

Recalling that $\ell = \ell_1z_1 = \sqrt{2}z_1 + \mathcal{O}(z_1^3)$, this gives the expansion for $\ell$ as desired. The case of fissures is simpler and was studied in [Reference Scheel and Stevens15].

2.3. Vacuum formation: hexagonal lattice symmetry

We next consider the case of hexagonal symmetry. For convenience of notation, we define the scaled lattice vectors and dual lattice vectors

(2.14)\begin{equation} \mathbf{e}_1 = \left( 1,0 \right), \quad \mathbf{e}_2 = \left( -\frac{1}{2}, \frac{\sqrt{3}}{2} \right) , \quad \mathbf{e}_1^\star = \left( 1, \frac{1}{\sqrt{3} } \right), \quad \mathbf{e}_2^\star = \left( 0, \frac{2}{\sqrt{3} } \right), \end{equation}

such that for $i,j\in\{1,2\}$,

(2.15)\begin{equation} \langle \mathbf{e}_i, \mathbf{e}_j^\star \rangle = \left\{ \begin{array}{ll}0,& i\neq j\\1,& i=j\end{array}\right. \end{equation}

The potential from (Equation 1.8) can then be written in the short form

(2.16)\begin{equation} V(x,y) = \cos \left( \mathbf{e}_1^\star \cdot (x,y) \right) + \cos \left( \mathbf{e}_2^\star \cdot (x,y) \right) + \cos \left( (\mathbf{e}_1^\star - \mathbf{e}_2^\star ) \cdot (x,y) \right). \end{equation}

Reasoning equivalent to the case of squares will show that solutions can be found in a three-dimensional reduced equation for the amplitudes of Fourier modes. Among those, we shall again find fissures and, similar to the bubbles in the square case, bubbles with dihedral isotropy $D_6$. Those bubbles arise in the form $u(x,y)=(A_0+A_1V(x,y))_+$. Different from the case of square lattices, however, the cases of $A_1 \gt 0$ and $A_1 \lt 0$ are different, leading to bubbles located in the corners of hexagons, taking an asymptotic shape of triangular corrections to circles, and bubbles with hexagonal corrections to circular shapes, located at the centres of hexagons; see Figure 3. We refer to the former patterns as triangles and to the latter ones as hexagons.

Figure 3. Inner figures (blue and green): An illustration of the vacuum bubbles on the hexagonal lattice. Left, triangles; right, hexagons. Blue regions represent areas of vacuum. In the case of triangles, bubbles are initially approximately circular, as pictured here; the triangular corrections become more apparent as ${\widetilde{\mu}}$ increases. All bubbles have (inner) radius $\ell$. Outer figures (pink and blue): example boundaries of triangle and hexagon vacuum regions from numerics, for various $\widetilde{\mu}$ (size not to scale with inner diagram). Triangles: $\widetilde{\mu} =$ 1e $-$5, 5e $-$5, 1e $-$4, 2.5e $-$4, 5e $-$4, 1e $-$3, 2e $-$3. Hexagons: $\widetilde{\mu} =$ 5e $-$5, .001, .005, .01, .03, .07.

Theorem 2

(Vacuum region scaling – hexagonal lattice). Let $V(x,y)$ as in (Equation 2.16), and ${\mu} = \frac{1}{\sqrt{3}\pi^2}+ {\widetilde{\mu}}$, with ${\widetilde{\mu}}$ sufficiently small. Then solutions $u$ to (Equation 2.4)–(Equation 2.6), invariant under point reflections $u(x,y)=u(-x,-y)$ and with vacuum regions, are of the form

(2.17)\begin{equation} u(x,y) =\left( A_0 + A_1 \cos \left( \mathbf{e}_1^\star \cdot (x,y) \right) + B_1 \cos \left( \mathbf{e}_2^\star \cdot (x,y) \right) +C_1 \cos \left( (\mathbf{e}_1^\star - \mathbf{e}_2^\star ) \cdot (x,y) \right)\right)_+, \end{equation}

where $f_+=\max(f,0)$. There are three solution branches, all supercritical, and parameterized by ${\widetilde{\mu}} \gt rsim 0$ as $A_0=A_0({\widetilde{\mu}}),\, A_1 = A_1({\widetilde{\mu}}),\, B_1 = B_1({\widetilde{\mu}}),\, C_1 = C_1({\widetilde{\mu}})$, locally unique up to translations in $x$ and $y$:

1. (1) Fissures: $B_1=C_1=0$,

   \begin{align*} A_0 =1-\frac{\pi}{2}{\widetilde{\mu}}+\mathcal{O}({\widetilde{\mu}}^2),\qquad A_1 =1+\left(\frac{3}{4\sqrt{2}}\pi^3{\widetilde{\mu}}\right)^{2/3}+\mathcal{O}({\widetilde{\mu}}). \end{align*}

   The half-width of the fissures $\ell$ and area $\mathcal{A}$ of vacuum region are

   \begin{equation*} \ell= \frac{3^{5/6}}{2^{4/3}}\pi{\widetilde{\mu}}^{1/3}+\mathcal{O}({\widetilde{\mu}}^{2/3}),\qquad \mathcal{A}= \frac{3^{5/6}}{2^{1/3}}\pi^2{\widetilde{\mu}}^{1/3}+\mathcal{O}({\widetilde{\mu}}^{2/3}). \end{equation*}

2. (2) Triangles: $A_1=B_1=C_1 \gt 0$,

   \begin{align*} A_0 = 1 -\frac{2\pi^2}{\sqrt{3}} {\widetilde{\mu}} + \mathcal{O}({\widetilde{\mu}}^{{3}/{2}}),\qquad A_1 = \frac{2}{3} + \frac{4\pi^{3/2}}{3\sqrt{3}}{\widetilde{\mu}}^{{1}/{2}} -\frac{4\pi^2 }{3\sqrt{3}}{\widetilde{\mu}} + \mathcal{O}({\widetilde{\mu}}^{{3}/{2}}). \end{align*}

   The inner radius of bubbles $\ell$ and area $\mathcal{A}$ of vacuum region are

   \begin{equation*} \ell= (12\pi^3 {\widetilde{\mu}})^{1/4} + \mathcal{O}({\widetilde{\mu}}^{{1}/{2}}),\qquad \mathcal{A}=2\sqrt{3}\pi^{5/2}{\widetilde{\mu}}^{1/2}+\mathcal{O}({\widetilde{\mu}}). \end{equation*}

3. (3) Hexagons: $A_1=B_1=C_1 \lt 0$,

   \begin{align*} A_0 = 1 -\frac{\pi^2}{18\sqrt{3}}{\widetilde{\mu}}+ \mathcal{O}({\widetilde{\mu}}^{{3}/{2}}),\qquad A_1 =-\frac{1}{3} - \frac{\sqrt{2\pi^3}}{9}{\widetilde{\mu}}^{1/2} +\frac{\pi^2}{18\sqrt{3}}{\widetilde{\mu}} + \mathcal{O}({\widetilde{\mu}}^{{3}/{2}}). \end{align*}

   The inner radius of bubbles $\ell$ and area $\mathcal{A}$ of vacuum region are

   \begin{equation*} \ell= (6\pi^3 {\widetilde{\mu}})^{1/4} + \mathcal{O}({\widetilde{\mu}}^{{3}/{4}}),\qquad \mathcal{A}=\sqrt{6}\pi^{5/2}{\widetilde{\mu}}^{1/2}+\mathcal{O}({\widetilde{\mu}}). \end{equation*}

Remark 2

(Vacuum area, hexagon versus triangle). Accounting for the number of vacuum bubbles per unit hexagon (see Figure 3), one can compare the areas of vacuum in the triangles and hexagon cases, to find

\begin{align*} \textrm{(hexagons)}& &&A_\textrm{hex} = \sqrt{6}\pi^{{5}/{2}}{\widetilde{\mu}}^{1/2} , \\ \textrm{(triangles)}& &&A_\textrm{tri} = 4\sqrt{3}\pi^{{5}/{2}}{\widetilde{\mu}}^{1/2} =2\sqrt{2}A_\textrm{hex} .&& \end{align*}

The triangle solutions thus possess not only a larger maximum density than the hexagons, but a larger area of vacuum.

We shall prove existence and establish expansions for triangles and hexagons in the next two sections. Existence and expansion for fissures again follow from the one-dimensional case [Reference Scheel and Stevens15]. Similar to the square case, Remark 1, the branches found here have maximal isotropy.

2.3.1. Proof of Theorem 2(ii) – triangles

We begin with the case $A_1 \gt 0$ and proceed analogously to the proof of Theorem 1. We scale

\begin{align*} A_0 &= 1+a_0, \qquad A_1 = \frac{2}{3}(1+ a_0 + z_1^2), \qquad a_0 = a_0(z_1), \end{align*}

for a small parameter $z_1$. The free-boundary condition (Equation 2.6), expanded in polar coordinates about the centres of the bubbles, gives

(2.18)\begin{equation} 0 = 1+a_0 + \frac{2}{3}\left(1+a_0+z_1^2\right)\left( -\frac{3}{2} + \frac{r^2}{2} \pm \ \frac{1}{6} \sin(3\theta) r^3 - \frac{r^4}{24} + \mathcal{O}(r^5)\right), \end{equation}

where the minus sign is taken for bubbles labelled $G_1$ and the plus sign is used for bubbles labelled $G_2$ (see Figure 3). We identify the bubble radius $\ell$ as the solution $r$ to (Equation 2.18), and further scale $\ell = \ell_1z_1$. In this scaling, after some algebraic manipulation, (Equation 2.18) becomes

(2.19)\begin{equation} 0 = - \frac{3}{2(1+a_0+z_1^2)} + \frac{\ell_1^2}{2}+ z_1\ell_1^3Q(\theta) - \frac{z_1^2\ell_1^4}{24} + \mathcal{O}(z_1^5), \end{equation}

with $Q(\theta) = \pm\frac{1}{6} \sin(3\theta).$ The implicit function theorem then yields a unique $\ell_1 = \ell_1(z_1,a_0,\theta)$, for $z_1, a_0$ sufficiently small and for arbitrary $\theta$. Expanding in $z_1, a_0$, we find

\begin{align*} \ell_1(z_1,a_0,\theta)&= \sqrt{3} - 3Q(\theta)z_1 +\left(\frac{45Q(\theta)^2}{2\sqrt{3} }-\frac{9}{8\sqrt{3} }\right) z_1^2 - \frac{\sqrt{3} }{2} a_0 + \mathcal{O}\left( |z_1^3| + |a_0z_1| + |a_0^2| \right). \end{align*}

Again, we can use this expression for $\ell_1(z_1,a_0,\theta)$ to solve for $a_0$ using the average-density condition (Equation 2.5). Substituting, we obtain

\begin{align*} 0 &= 2\sqrt{3}\pi^2 - \left(\iint_{\Omega}\left(A_0 + A_1V(x,y)\right)dydx - \iint_{ \Omega_0^c}\left(A_0 + A_1V(x,y)\right)dydx\right)\\ &= -2\sqrt{3}\pi^2\cdot a_0 + \iint_{\Omega_0^c}\left(\left(1 + a_0\right) + \frac{2}{3}(1 + a_0 + z_1^2)\left(2\cos (x) \cos\left( \frac{y}{\sqrt{3} } \right) + \cos \left( \frac{2y}{\sqrt{3} }\right) \right) \right)dydx \\ &=-2\sqrt{3}\pi^2\cdot a_0 + 2\int_0^{2\pi}\int_0^{z_1\ell_1(z_1,a_0,\theta)} \left[(1+a_0)\frac{r^2}{3} + z_1^2\left( -1 + \frac{r^2}{3} \right)+ A_1Q(\theta)r^3 + \mathcal{O}(r^4)\right]r \ drd\theta \\ &= -2\sqrt{3}\pi^2\cdot a_0 - 3\pi(1-a_0) z_1^4 + \mathcal{O} \left( z_1^6 \right), \end{align*}

and solve by the implicit function theorem to find, for $z_0$ sufficiently small,

(2.20)\begin{equation} a_0 = -\frac{\sqrt{3} }{2\pi} z_1^4 + \mathcal{O} \left( z_1^6 \right). \end{equation}

Lastly, we use (Equation 2.4) to relate these quantities to ${\widetilde{\mu}}$:

\begin{align*} 0 =&\, V(x,y)A_1- \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)\iint_{\Omega -\Omega_0^c} V(x-\xi,y-\eta) \left( A_0 + A_1 \cdot V(\xi,\eta) \right) d\xi d\eta \\ =&\, V(x,y)\left(A_1- \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)(A_1\sqrt{3}\pi^2)\right) \\ &+ \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)\iint_{\Omega_0^c} V(x-\xi,y-\eta) \left( A_0 + A_1 \cdot V(\xi,\eta) \right) d\xi d\eta. \end{align*}

Equating coefficients of $V(x,y)$, using polar coordinates for the integral over the bubbles, we find

\begin{align*} 0 =& -\!\sqrt{3}{\widetilde{\mu}}\pi^2A_1 + \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)\iint_{\Omega_0^c}\frac{1}{3}V(\xi,\eta)(A_0+A_1 V(\xi,\eta))d\xi d\eta\\ =& -\!\sqrt{3}{\widetilde{\mu}}\pi^2A_1 + 2\left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)\int_0^{2\pi}\int_0^{z_1\ell_1}\left( \frac{1}{2}z_1^2 - \frac{1}{6}r^2 -\frac{Q(\theta)}{3}r^3 + \mathcal{O}(r^4,r^2z_1^2) \right)r \ dr d\theta \\ =& -\!\frac{2}{\sqrt{3}} {\widetilde{\mu}} \pi^2 \left(1 + z_1^2 -\frac{\sqrt{3} }{2\pi} z_1^4 \right) + \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)\left(\frac{3\pi}{2}z_1^4 + \mathcal{O}(z_1^6)\right). \end{align*}

As in the square symmetry case for bubbles, we find that $z_1$ scales like ${\widetilde{\mu}}^{{1}/{4}}$, and we obtain a unique $z_1 = z_1({\widetilde{\mu}}^{{1}/{4}})$, for ${\widetilde{\mu}}^{1/4}$ sufficiently small, by the implicit function theorem, with

(2.21)\begin{equation} z_1 = \left(\frac{4\pi^3}{3}{\widetilde{\mu}}\right)^{1/4} + \mathcal{O}({\widetilde{\mu}}^{1/2}). \end{equation}

Using the relation $\ell = \sqrt{3}z_1 +\mathcal{O}(z_1^2)$, this gives us the expression in Theorem 2 as desired.

2.3.2. Proof of Theorem 2(iii) – hexagons

We now turn to the case $A_1 \lt 0$. and proceed analogously. We scale

\begin{align*} A_0 &= 1+a_0, \qquad A_1 =- \frac{1}{3}\left(1 + a_0 + z_1^2\right), \qquad a_0 = a_0(z_1), \end{align*}

for a small parameter $z_1$. Setting the solution equal to 0 in the free-boundary Equation 2.6 and expanding in polar coordinates yields

(2.22)\begin{equation} 0= 1+a_0 - \frac{1}{3}\left(1+a_0 +z_1^2\right) \left( 3 - r^2 + \frac{r^4}{12} + \mathcal{O}(r^6) \right), \end{equation}

where $r$ in (Equation 2.22) will be the radius $\ell$ of the bubble. Scaling $\ell = \ell_1z_1$, and simplifying, we find

(2.23)\begin{equation} 0 = \frac{3}{1+a_0 + z_1^2} - \ell_1^2 + \frac{\ell_1^4 z_1^2}{12} + \mathcal{O}(\ell_1^6 z_1^4), \end{equation}

which we can solve by the implicit function theorem to obtain $\ell_1 = \ell_1(z_1,a_0,\theta)$, for $z_1, a_0$ small, and for arbitrary $\theta$. Expanding in $z_1, a_0$, we find

(2.24)\begin{equation} \ell_1 =\sqrt{3} -\frac{3\sqrt{3}}{8}z_0^2 -\frac{\sqrt{3}}{2}a_0 + \mathcal{O}(a_0^2 + z_1^4). \end{equation}

We substitute this expression for $\ell(z_1,a_0,\theta)$ to solve for $a_0$ in the average-density condition (Equation 2.5),

\begin{align*} 0 &= 2\sqrt{3}\pi^2 - \iint_{\Omega}\left(A_0 + A_1V(x,y)\right)dydx + \iint_{ \Omega_0^c}\left(A_0 + A_1V(x,y)\right)dydx\\ &= -2\sqrt{3}\pi^2\cdot a_0 + \iint_{\Omega_0^c}\left(\left(1 + a_0\right) -\frac{1}{3}(1 + a_0 + z_1^2)\left(2\cos (x) \cos\left( \frac{y}{\sqrt{3} } \right) + \cos \left( \frac{2y}{\sqrt{3} }\right) \right) \right)dydx \\ &=-2\sqrt{3}\pi^2\cdot a_0 + \int_0^{2\pi}\int_0^{z_1\ell_1(z_1,a_0,\theta)} \left[(1+a_0)\frac{r^2}{3} + z_1^2\left( -1 + \frac{r^2}{3} \right)+ \mathcal{O}(r^4)\right]r \ drd\theta \\ &= -2\sqrt{3}\pi^2\cdot a_0 - \frac{3\pi}{2}(1-a_0) z_1^4 + \mathcal{O} \left( z_1^6 \right), \end{align*}

and solve by the implicit function theorem to find, for $z_0$ sufficiently small,

(2.25)\begin{equation} a_0(z_1) = -\frac{\sqrt{3} }{4\pi} z_1^4 + \mathcal{O} \left( z_1^6 \right). \end{equation}

Lastly, we use (Equation 2.4) to relate previous quantities to ${\widetilde{\mu}}$:

\begin{align*} 0 =&\, V(x,y)A_1- \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)\iint_{\Omega -\Omega_0^c} V(x-\xi,y-\eta) \left( A_0 + A_1 \cdot V(\xi,\eta) \right) d\xi d\eta \\ =&\, V(x,y)\left(A_1 - \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)(A_1\sqrt{3}\pi^2) \right)\\ &+ \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)\iint_{\Omega_0^c} V(x-\xi,y-\eta) \left( A_0 - A_1 V(\xi,\eta) \right) dxdy. \end{align*}

Equating coefficients of $V(x,y)$, using polar coordinates for the integral over the bubble, we calculate

\begin{align*} 0 =& -\sqrt{3}{\widetilde{\mu}}\pi^2A_1 + \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)\iint_{\Omega_0^c}\frac{1}{3}V(\xi,\eta)(A_0+A_1 V(\xi,\eta))d\xi d\eta\\ =&\, \frac{{\widetilde{\mu}}\pi^2}{\sqrt{3}}\left(1 + z_1^2 -\frac{\sqrt{3} }{4\pi} z_1^4 \right) + \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)\int_0^{2\pi}\int_0^{z_1\ell_1}\left( -z_1^2 + \frac{r^3}{3} + \mathcal{O}(r^4, r^2z_1^2) \right)r \ dr d\theta \\ =&\, \frac{{\widetilde{\mu}}\pi^2}{\sqrt{3}}\left(1 + z_1^2 -\frac{\sqrt{3} }{4\pi} z_1^4 \right) + \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)\left(-\frac{3\pi}{2}z_1^4 + \mathcal{O}(z_1^6)\right). \end{align*}

Again, we find that $z_1$ scales like ${\widetilde{\mu}}^{1/4}$, and we obtain a unique $z_1 = z_1({\widetilde{\mu}}^{1/4})$, for ${\widetilde{\mu}}^{1/4}$ sufficiently small, by the implicit function theorem, with

(2.26)\begin{equation} z_1 = \left(\frac{2\pi^3}{3}{\widetilde{\mu}}\right)^{1/4} + \mathcal{O}({\widetilde{\mu}}^{3/4}). \end{equation}

Using the relation $\ell = \sqrt{3}z_1 + \mathcal{O}(z_1^3)$, this gives us the expression in Theorem 2 as desired.

3.1. Square lattice symmetry

We begin with the case of the square lattice.

Proposition 1

(Diffusive corrections – square bubbles). For $\varepsilon \gt rsim 0$, the McKean–Vlasov equation

(3.2)\begin{equation} u_t = \varepsilon \Delta u + \nabla \cdot \left(u \, \nabla\left( u - \left(\frac{1}{\pi^2}+{\widetilde{\mu}}\right)(V*u)\right)\right), \end{equation}

possesses an almost-vertical branch of stationary solutions $u_*(x,y;A,\varepsilon), 0 \lt A \lt 1$, $D_4$-symmetric and periodic, for ${\widetilde{\mu}} = {\widetilde{\mu}}_*(A,\varepsilon)$. Normalizing average density to 1, it has the expansion

(3.3)\begin{equation} {\widetilde{\mu}}_*(A,\varepsilon) = \varepsilon {\widetilde{\mu}}_1(A) + \mathcal{O}(\varepsilon^2), \qquad u_*(x,y;A,\varepsilon) = 1 + A\cos(x)\cos(y) + \mathcal{O}(\varepsilon), \end{equation}

where

(3.4)\begin{equation} A=\frac{1}{\pi^2}\int_{-\pi}^\pi \int_{-\pi}^\pi u_*(x,y)\cos(x)\cos(y)dxdy, \end{equation}

and, explicitly,

(3.5)\begin{equation} {\widetilde{\mu}}_1(A) =\frac{2}{\pi^3A^2} \left(2 \pi - 2 \left(\textrm{E}(A^2) + \sqrt{1 - A^2} \textrm{E}\left(A^2/(-1 + A^2)\right)\right)\right), \end{equation}

where $\textit{E}(k)$ is the complete elliptic integral, $E(k) = \int_{0}^{\frac{\pi}{2}}\sqrt{1-k^2 \sin^2\theta}d\theta$. The limits of ${\widetilde{\mu}}_1(A)$ are ${\widetilde{\mu}}_1(0) = \frac{1}{\pi^2}, {\widetilde{\mu}}_1(1) = \frac{4}{\pi^2}-\frac{8}{\pi^3}.$

Remark 5

(Diffusive corrections – square fissures). For fissures, ${\widetilde{\mu}} = {\widetilde{\mu}}_1\varepsilon + \mathcal{O}(\varepsilon^2)$,

\begin{equation*}{\widetilde{\mu}}_1(A) = \frac{4\pi^2\left(\frac{1-\sqrt{1-A^2}}{A^2}\right)}{2\pi^4} = \frac{2}{\pi^2}\left(\frac{1-\sqrt{1-A^2}}{A^2}\right), \text{with limits } {\widetilde{\mu}}_1(0) = \frac{1}{\pi^2}, {\widetilde{\mu}}_1(1)=\frac{2}{\pi^2}. \end{equation*}

The proof in this case is analogous to the proof in [Reference Scheel and Stevens15].

Proof. Noting that $\Delta u = \nabla \cdot \left(u\, \nabla (\log u)\right)$, we can write the fixed-point equation in the form

(3.6)\begin{equation} 0 = \nabla \cdot \left(u \,\nabla \left(\varepsilon \log u + u - \left(\frac{1}{\pi^2}+{\widetilde{\mu}}\right)(V*u)\right)\right). \end{equation}

Restricting to strictly positive solutions $u$, this is equivalent to

(3.7)\begin{equation} \varepsilon \log u + u - \left(\frac{1}{\pi^2}+{\widetilde{\mu}}\right)(V*u) = \rho, \end{equation}

for $\rho$ a constant. We append normalization conditions, solve for $\rho$ and define $F,F_0,F_1$ through

(3.8)\begin{equation} \begin{aligned} F(u,{\widetilde{\mu}},\varepsilon) &:= \tilde{F}(u,{\widetilde{\mu}},\varepsilon)-\rho,\quad \tilde{F}(u,{\widetilde{\mu}},\varepsilon)=\varepsilon \log u + u - \left(\frac{1}{\pi^2}+{\widetilde{\mu}}\right)(V*u), \quad \rho=\frac{1}{4\pi^2}\iint_\Omega \tilde{F} \\ F_0(u) &:= \int_{-\pi}^\pi \int_{-\pi}^\pi u(x,y)dxdy - 4\pi^2 \\ F_1(u) &:= \int_{-\pi}^\pi \int_{-\pi}^\pi u(x,y)\cos(x)\cos(y)dxdy - A\pi^2, \end{aligned} \end{equation}

where the parameter $A$ is fixed, $|A| \lt 1$. The system (Equation 3.8) defines a map

\begin{equation*} G(u,{\widetilde{\mu}},\varepsilon)=(F,F_0,F_1)(u,{\widetilde{\mu}},\varepsilon): \left( H^2_{\mathrm{e,p}} \right) \times \mathbb{R}^2 \to \overset{\circ}{L}{}^2_{\mathrm{e,p}} \times \mathbb{R}^2, \end{equation*}

where the subscripts $\{\mathrm{e},\mathrm{p}\}$ indicate the restriction of function spaces to $D_4$-symmetric, periodic functions, and $\overset{\circ}{L}{}^2$ indicates the restriction to zero average. One verifies that $G$ is well-defined and smooth, and

\begin{equation*} G(u_*^0(\cdot,\cdot;A),0,0) = 0, \quad \textrm{for } \ \ u_*^0(x,y;A) = 1 + A\cos(x)\cos(y). \end{equation*}

The derivative at $(u_*^0,0,0)$ is given by

\begin{equation*} \mathcal{A} = \begin{pmatrix} \mathcal{L} & \partial_{\widetilde{\mu}} F_* & \partial_\varepsilon F_* \\ \langle 1, \cdot \rangle & 0 & 0 \\ \langle \cos(x)\cos(y), \cdot\rangle & 0 & 0 \end{pmatrix}, \end{equation*}

where

\begin{align*} \mathcal{L}u = u - \left(\frac{1}{\pi^2}+{\widetilde{\mu}}\right)(V*u),\qquad \partial_{\widetilde{\mu}} F_* = -\pi^2A\cos(x)\cos(y),\qquad \partial_\varepsilon F_* = \log(1 + A\cos(x)\cos(y)). \end{align*}

The operator $\mathcal{L}$ is strictly elliptic since $|A| \lt 1$, and self-adjoint with domain $H^2_{\textrm{e,p}}$. It is therefore Fredholm index 1 due to the restriction of the codomain to average 0 functions. Bordering lemmas for Fredholm operators then imply that $\mathcal{A}$ is also Fredholm of index 1, so that we expect a one-dimensional family of solutions.

We turn ourselves to showing that the map $D_{u,{\widetilde{\mu}}}G(u_*^0(\cdot,\cdot;A),0,0)$, that is, the first two columns of $\mathcal{A}$, is invertible. Clearly, the kernel of $\mathcal{L}$ is spanned by $\{1, \cos(x)\cos(y) \}$, and the cokernel, by restriction to average zero functions, is spanned by $\{\cos(x)\cos(y) \}$.

As a consequence, we find that $\partial_{\widetilde{\mu}} F_* \notin \textrm{Rg}(\mathcal{L})$, since

\begin{equation*} \langle\pi^2A \cos(x)\cos(y), \cos(x)\cos(y) \rangle = \pi^4A \neq 0. \end{equation*}

We can thus see that $\mathcal{A}(u,{\widetilde{\mu}},0) = 0$ if and only if $ u = 0, {\widetilde{\mu}} = 0$, since the second and third components of $\mathcal{A}$ imply $u \notin \ker(\mathcal{L})$, and therefore $\mathcal{L}u \neq \partial_{\widetilde{\mu}} F_* {\widetilde{\mu}}$, unless $u = 0, {\widetilde{\mu}} = 0$. We therefore have that $D_{u,{\widetilde{\mu}}}G(u_*^0,0,0)$ is invertible, and we can solve for $u, {\widetilde{\mu}}$ as functions of $\varepsilon$. We then obtain the leading order expansion of ${\widetilde{\mu}}$ by projecting the first equation onto the cokernel. We evaluate

\begin{align*} \langle \partial_\varepsilon F_*, \cos(x)\cos(y)\rangle &= \int_{-\pi}^\pi \int_{-\pi}^\pi \cos(x)\cos(y)\cdot \log\left(1 + A\cos(x)\cos(y)\right)dxdy \\ &= \frac{2\pi}{A} \left(2 \pi - 2 \left(\textrm{E}(A^2) + \sqrt{1 - A^2} \textrm{E}(A^2/(-1 + A^2))\right)\right) \neq 0, \end{align*}

and thus find

\begin{align*} {\widetilde{\mu}} &= {\widetilde{\mu}}_1\varepsilon + \mathcal{O}(\varepsilon^2), \\ {\widetilde{\mu}}_1 &= -\frac{\left\langle \partial_\varepsilon F_*,\cos(x)\cos(y) \right\rangle }{\left\langle \partial_{\widetilde{\mu}} F_*,\cos(x)\cos(y) \right\rangle } = \frac{2}{\pi^3A^2} \left(2 \pi - 2 \left(\textrm{E}\left(A^2\right) + \sqrt{1 - A^2} \textrm{E}\left(A^2/(-1 + A^2)\right)\right)\right). \end{align*}

Remark 6

(Diffusive corrections – no elliptical bubbles). If one were to attempt to repeat the argument with $u_*^0 = 1 + A\cos (x) \cos (y) + B\sin(x)\sin(y) , $ for $B \neq A, \ \ A,B \neq 0$, one would not find nontrivial solutions for $\varepsilon \gt 0$. The corresponding nonlinear function

(3.9)\begin{equation} \begin{aligned} \widetilde{F}(u,{\widetilde{\mu}},\varepsilon) &:= \varepsilon \log u + u - \left(\frac{1}{\pi^2}+{\widetilde{\mu}}\right)(V*u) - \rho \\ \widetilde{F}_0(u) &:= \int_{-\pi}^\pi \int_{-\pi}^\pi u(x,y)dxdy - 4\pi^2 \\ \widetilde{F}_1(u) &:= \int_{-\pi}^\pi \int_{-\pi}^\pi u(x,y)\cos(x)\cos(y)dxdy - A\pi^2 \\ \widetilde{F}_2(u) &:=\int_{-\pi}^\pi \int_{-\pi}^\pi u(x,y)\sin(x)\sin(y)dxdy - B\pi^2 , \end{aligned} \end{equation}

where $\rho$ is the average of the other terms, defines a map $\widetilde{G}(u,{\widetilde{\mu}},\varepsilon): H^2_\mathrm{s,p}\times \mathbb{R}^2 \to \overset{\circ}{L}{}^2_{\mathrm{s,p}} \times \mathbb{R}^2$, where we can now restrict function spaces to $D_2$-symmetric, periodic functions, only, generated by reflections $x\leftrightarrow y$ and rotations by $\pi$, denoted by subscripts $\mathrm{s,p}$. It has linearization

\begin{equation*} \mathcal{\widetilde{A}} = \begin{pmatrix} \mathcal{L} & \partial_{\widetilde{\mu}} F_* & \partial_\varepsilon F_* \\ \langle 1, \cdot \rangle & 0 & 0 \\ \langle \cos(x)\cos(y), \cdot\rangle & 0 & 0 \\ \langle \sin(x)\sin(y), \cdot\rangle & 0 & 0 \end{pmatrix}, \end{equation*}

which is no longer Fredholm index 1, but Fredholm index 0 due to the additional normalization condition.

In order to determine whether $\mathcal{\widetilde{A}} $ possesses a kernel, we need to check whether the two vectors $\partial_{\widetilde{\mu}} F_* = A\cos(x)\cos(y)+B\sin(x)\sin(y)$ and $\partial_\varepsilon F_* = \log(1 + A\cos(x)\cos(y)+B\sin(x)\sin(y))$ span a complement of the range, that is, if their orthogonal projections on the basis $\{\cos(x)\cos(y), \sin(x)\sin(y)\}$ are linearly independent.

Evaluating the scalar products and changing integration variables leads to the condition

\begin{equation*} \int_{\eta=0}^{2\pi} \int_{\xi=0}^{2\pi} (A\cos\xi - B \cos\eta)\log(1+A\cos\eta+B\cos\xi)d\xi d\eta\neq 0, \end{equation*}

which after partial integration simplifies to

\begin{equation*} AB\int_{\eta=0}^{2\pi} \int_{\xi=0}^{2\pi} \frac{\cos^2\xi-\cos^2\eta}{1+A\cos\eta+B\cos\xi} d\xi d\eta\neq 0. \end{equation*}

This integral clearly vanishes when $A=B$. We evaluated the integral numerically and found that it is strictly monotone on lines $A+B=const$, hence nonzero for $A\neq B$. As a consequence, $\mathcal{\widetilde{A}} $ is invertible and the solution branch at $\varepsilon=0$ is unique, that is, there are no solutions with submaximal isotropy for $\varepsilon \gt rsim 0$. This is a notable difference from the $\varepsilon = 0$ case, where the vertical branch does contain mixed-mode solutions.

3.2. Hexagonal lattice symmetry

We now turn to the case of the hexagonal lattice.

Proposition 2

(Diffusive corrections – triangles and hexagons). For $\varepsilon \gt rsim 0$, the McKean–Vlasov equation

(3.10)\begin{equation} u_t = \varepsilon \Delta u + \nabla \cdot \left(u\,\nabla \left( u - \left(\frac{1}{\sqrt{3}\pi^2}+{\widetilde{\mu}}\right)(V*u)\right)\right) \end{equation}

possesses an almost-vertical branch of stationary solutions $u_*(x,y;A,\varepsilon)$ with $D_6$ isotropy and average 1, parameterized by $-\frac{1}{3} \lt A \lt \frac{2}{3}$, $A \neq 0$, for ${\widetilde{\mu}} = {\widetilde{\mu}}_*(A,\varepsilon)$ with expansion

(3.11)\begin{equation} {\widetilde{\mu}}_*(A,\varepsilon) = \varepsilon {\widetilde{\mu}}_1(A) + \mathcal{O}(\varepsilon^2), \qquad u_*(x,y;A,\varepsilon) = 1 + AV(x,y) + \mathcal{O}(\varepsilon), \end{equation}

where

(3.12)\begin{equation} A= \frac{1}{3\sqrt{3}\pi^2} \iint_{\Omega} u_*(x,y)V(x,y)dxdy , \end{equation}

and, explicitly,

(3.13)\begin{equation} {\widetilde{\mu}}_1(A)=\frac{1}{3\pi^3 A^2}\int_{-\frac{\sqrt{3}\pi}{2}}^{\frac{\sqrt{3}\pi}{2}}\left(1 - \sqrt{\left(1+A\cos\left(\frac{2y}{\sqrt{3}}\right)\right)^2 -4A^2\cos^2\left(\frac{y}{\sqrt{3}}\right) }\right)dy. \end{equation}

The limits of ${\widetilde{\mu}}_1(A)$ are ${\widetilde{\mu}}_1(0) = \frac{1}{\sqrt{3}\pi^2}, {\widetilde{\mu}}_1(\frac{2}{3}) = \frac{2}{\sqrt{3}\pi^2}-\frac{3}{2\pi^3}$, ${\widetilde{\mu}}_1(-\frac{1}{3}) = \frac{5}{\sqrt{3}\pi^2}-\frac{6}{\pi^3}$.

Remark 7

(Diffusive corrections – hexagonal fissures). For fissures in this setting, we find

\begin{equation*}{\widetilde{\mu}}_1 = \frac{2 (1 - \sqrt{1 - A^2})}{{\sqrt{3}A^2\pi^2}}, \quad \text{with limits } {\widetilde{\mu}}_1(0) = \frac{1}{\sqrt{3}\pi^2},\ {\widetilde{\mu}}_1(1) = {\widetilde{\mu}}_1(-1) = \frac{2}{\sqrt{3}\pi^2}, \end{equation*}

simply adapting the results in [Reference Scheel and Stevens15].

Proof. Using $\Delta u = \nabla \cdot \left(u \,\nabla (\log u)\right)$, we can again write the fixed-point equation in the form

(3.14)\begin{equation} 0 = \nabla \cdot \left(u \, \nabla \left(\varepsilon \log u + u - \left(\frac{1}{\sqrt{3}\pi^2}+{\widetilde{\mu}}\right)(V*u)\right)\right). \end{equation}

Restricting to strictly positive solutions $u$, this is equivalent to

(3.15)\begin{equation} \varepsilon \log u + u - \left(\frac{1}{\sqrt{3}\pi^2}+{\widetilde{\mu}}\right)(V*u) = \rho, \end{equation}

for $\rho$ a constant, which we necessarily choose as the average of the left-hand side. We consider the nonlinear functions

(3.16)\begin{equation} \begin{aligned} F(u,{\widetilde{\mu}},\varepsilon) &:= \varepsilon \log u + u - \left(\frac{1}{\sqrt{3}\pi^2}+{\widetilde{\mu}}\right)(V*u) - \rho \\ F_0(u) &:= \iint_{\Omega} u(x,y)dxdy - 2\sqrt{3}\pi^2 \\ F_1(u) &:= \iint_{\Omega} u(x,y)V(x,y) - 3\sqrt{3}A\pi^2, \end{aligned} \end{equation}

where $A$ is fixed, $A \neq 0, -\frac{1}{3} \lt A \lt \frac{2}{3}$. These functions (Equation 3.16) define a map

\begin{equation*}G(u,{\widetilde{\mu}},\varepsilon): H^2_{\textrm{e,p}}\times \mathbb{R}^2 \to \overset{\circ}{L}{}^2_{\textrm{e,p}} \times \mathbb{R}^2 ,\end{equation*}

where the subscript $\{\textrm{e,p}\}$ refers to the fact that functions are lattice-periodic and invariant under $D_6$, the symmetry of the fundamental hexagon and $\overset{\circ}{L}{}^2_{\textrm{e,p}} $ refers to the space of symmetric, mean-zero functions. One verifies that $G$ is well-defined and smooth, and

\begin{equation*} G(u_*^0(\cdot,\cdot;A),0,0) = 0, \quad \textrm{for } \ \ u_*^0(x,y;A) = 1 + AV(x,y). \end{equation*}

The derivative at $(u_*^0,0,0)$ is given by

\begin{equation*} \mathcal{A} = \begin{pmatrix} \mathcal{L} & \partial_{\widetilde{\mu}} F_* & \partial_\varepsilon F_* \\ \langle 1, \cdot \rangle & 0 & 0 \\ \langle V(x,y), \cdot\rangle & 0 & 0 \end{pmatrix}, \end{equation*}

where

\begin{align*} \mathcal{L}u = u - \left(\frac{1}{\sqrt{3}\pi^2}+{\widetilde{\mu}}\right)(V*u),\qquad \partial_{\widetilde{\mu}} F_* = -\sqrt{3} \pi^2AV(x),\qquad \partial_\varepsilon F_* = \log(1 + AV(x)). \end{align*}

The operator $\mathcal{L}$ is strictly elliptic since $\frac{1}{3} \lt A \lt \frac{2}{3}$, and self-adjoint with domain $H^2_{\textrm{e,p}}$. It is therefore Fredholm index 1 due to the restriction to average 0 functions in the codomain. Bordering lemmas for Fredholm operators then imply that $\mathcal{A}$ is Fredholm of index 1, so that we expect a one-dimensional family of solutions.

We now turn to showing that the map $D_{u,{\widetilde{\mu}}}G(u_*^0(\cdot,\cdot;A),0,0)$ consisting of the first two columns of $\mathcal{A}$, is invertible. We first notice that the kernel of $\mathcal{L}$ is explicitly the span of $\{1, V(x,y) \}$, and the cokernel is the span of $\{V(x,y) \}$. As a consequence, we find that $\partial_{\widetilde{\mu}} F_* \notin \textrm{Rg}(\mathcal{L})$, through

\begin{equation*} - \langle \sqrt{3}\pi^2A V(x,y), V(x,y) \rangle = -9\pi^4A \neq 0. \end{equation*}

We can thus see that $\mathcal{A}(u,{\widetilde{\mu}}) = 0$ if and only if $ u = 0, {\widetilde{\mu}} = 0$, since the second and third components of $\mathcal{A}$ imply $u \notin \ker(\mathcal{L})$, and therefore $\mathcal{L}u \neq \partial_{\widetilde{\mu}} F_* {\widetilde{\mu}}$, unless $u = 0, {\widetilde{\mu}} = 0$. We therefore have that $D_{u,{\widetilde{\mu}}}G(u_*^0,0,0)$ is invertible, and we can solve for $u, {\widetilde{\mu}}$ as functions of $\varepsilon$. We then obtain the leading order expansion of ${\widetilde{\mu}}$ by projecting the first equation onto the cokernel. We evaluate

\begin{align*} \langle \partial_\varepsilon F_*, V(x,y)\rangle &= \int_{-\frac{\sqrt{3}\pi}{2}}^{\frac{\sqrt{3}\pi}{2}}\int_{-\pi}^\pi V(x,y)\cdot \log\left(1 + AV(x,y)\right)dxdy \\ &= \int_{-\frac{\sqrt{3}\pi}{2}}^{\frac{\sqrt{3}\pi}{2}}\frac{3\pi}{A}\left(1 - \sqrt{\left(1+A\cos\left(\frac{2y}{\sqrt{3}}\right)\right)^2 -4A^2\cos^2\left(\frac{y}{\sqrt{3}}\right) }\right)dy \\ &= \int_{-\frac{\sqrt{3}\pi}{2}}^{\frac{\sqrt{3}\pi}{2}}\frac{3\pi}{A}\left(1 - \sqrt{\left(1+A\cos\left(\frac{2y}{\sqrt{3}}\right)\right)^2 -4A^2\cos^2\left(\frac{y}{\sqrt{3}}\right) }+A\cos\left(\frac{2y}{\sqrt{3}}\right)\right)dy \\ & \gt 0, \end{align*}

since the integrand is nonnegative and not identically equal to 0. We thus find

\begin{align*} {\widetilde{\mu}} &= {\widetilde{\mu}}_1\varepsilon + \mathcal{O}(\varepsilon^2), \\ {\widetilde{\mu}}_1 &= -\frac{\left\langle \partial_\varepsilon F_*,\cos(x)\cos(y) \right\rangle }{\left\langle \partial_{\widetilde{\mu}} F_*,\cos(x)\cos(y) \right\rangle } \\ \ &= \int_{-\frac{\sqrt{3}\pi}{2}}^{\frac{\sqrt{3}\pi}{2}}\frac{1}{3\pi^3 A^2}\left(1 - \sqrt{\left(1+A\cos\left(\frac{2y}{\sqrt{3}}\right)\right)^2 -4A^2\cos^2\left(\frac{y}{\sqrt{3}}\right) }\right)dy, \end{align*}

with explicit limits ${\widetilde{\mu}}_1(0) = \frac{1}{\sqrt{3}\pi^2}, {\widetilde{\mu}}_1(\frac{2}{3}) = \frac{2}{\sqrt{3}\pi^2}-\frac{3}{2\pi^3}$, ${\widetilde{\mu}}_1(-\frac{1}{3}) = \frac{5}{\sqrt{3}\pi^2}-\frac{6}{\pi^3}$.

Numerical comparison and agreeement for the almost vertical branches is shown in Figure 4, for both square and hexagonal symmetry.

Figure 4. Numerical agreement and comparison for the almost-vertical branches. Dots are numerical; curves are the $\mathcal{O}(\varepsilon)$ predictions from (Equation 3.5) and (Equation 3.13); square symmetry (left) and hexagonal symmetry (right) shown. Bubble branches are in green; fissure branches are in blue.

3.3. Stability near onset: centre-manifold expansions

The characterization of the almost-vertical branches also encodes information about the stability of the solutions. The explicit dependence of the solutions on ${\widetilde{\mu}}$, along with the symmetries of the system, allow one to compute the coefficients in the centre-manifold expansion for small ${\widetilde{\mu}}$.

Square lattice. We take $\varepsilon \gt 0$ fixed, small, and let ${\widetilde{\mu}}_c = {\widetilde{\mu}}_*(0,\varepsilon)$ defined in Prop. 1. We can then perform a parameter-dependent centre-manifold reduction in the space of functions that are invariant under the reflections $x\leftrightarrow y$ and $x\leftrightarrow -y$. We decompose $u = u_c + u_h$, $u_c = A\cos(x-y)+B\cos(x+y)$, parameterizing the centre eigenspace, and $u_h\perp u_c$. The centre manifold is then given as a graph $u_h=h(u_c)$, with $ h(0) = 0$ and $h'(0)=0$, where we suppress parameter dependence. Within the space of even functions, we have the action of the dihedral group $D_4$ generated by the rotation of order 4 in the plane and the reflection $x\leftrightarrow -x$, as well as translations by $\pi$ in $x$ and $y$. On the centre eigenspace, and therefore on the reduced vector field, we find they induce a standard action of $D_2$, generated by $A\leftrightarrow B$ and $A\rightarrow -A$; see for instance [Reference Golubitsky, Stewart and Schaeffer12]. This enforces a Taylor expansion of the reduced vector field of the form,

(3.17)\begin{equation} \begin{aligned} A' &= ({\widetilde{\mu}}-{\widetilde{\mu}}_c)A + A(\alpha A^2 + \beta B^2) + \mathcal{O}(({\widetilde{\mu}}-{\widetilde{\mu}}_c)^2(|A|+|B|) +|A|^5+|B|^5)\\ B' &= ({\widetilde{\mu}}-{\widetilde{\mu}}_c)B +B(\alpha B^2 + \beta A^2) + \mathcal{O}(({\widetilde{\mu}}-{\widetilde{\mu}}_c)^2 (|A|+|B|)+|A|^5+|B|^5). \end{aligned} \end{equation}

Fissure solutions correspond to $B=0,$ and bubbles correspond to $A=B$, both one-dimensional invariant lines in this reduced phase plane description, enforced by symmetry. Thus steady-state solutions have, for small ${\widetilde{\mu}}-{\widetilde{\mu}}_c$, leading-order amplitudes

\begin{equation*}A = \pm\sqrt{\frac{-({\widetilde{\mu}}-{\widetilde{\mu}}_c)}{\alpha}} \quad \textrm{(fissures); \qquad and } \qquad B = \pm\sqrt{\frac{-({\widetilde{\mu}}-{\widetilde{\mu}}_c)}{\alpha+\beta}} \quad \textrm{(bubbles)}.\end{equation*}

We can compare these expressions to the asymptotics of ${\widetilde{\mu}}_1$ derived above in §3.1. Starting with fissure solutions, we expand ${\widetilde{\mu}}_1(A) = \frac{2}{\pi^2}\left(\frac{1-\sqrt{1-A^2}}{A^2}\right)$ about $A=0$, and find

\begin{equation*} {\widetilde{\mu}}_1 - {\widetilde{\mu}}_1(0) = \frac{1}{4\pi^2}A^2 + \mathcal{O}(A^4). \end{equation*}

This implies that $\alpha = \frac{-1}{4\pi^2}\varepsilon + \mathcal{O}(\varepsilon^2)$.

Similarly, we expand the expression for ${\widetilde{\mu}}_1(A)$ from §3.1 for the bubble solutions, noting that in the coordinates used there, $A$ corresponds to $2B$ here. We then obtain

\begin{equation*} {\widetilde{\mu}}_1 - {\widetilde{\mu}}_1(0) = \frac{3}{4\pi^2}B^2 + \mathcal{O}(B^4), \end{equation*}

from which we deduce $\alpha + \beta = \frac{-3}{4\pi^2}\varepsilon + \mathcal{O}(\varepsilon^2)$, so $\beta = \frac{-1}{2\pi^2}\varepsilon + \mathcal{O}(\varepsilon^2)$. Summarizing, we have

(3.18)\begin{equation} \alpha= \frac{-4}{2\pi^2}\varepsilon + \mathcal{O}(\varepsilon^2),\qquad \beta = \frac{-1}{2\pi^2}\varepsilon + \mathcal{O}(\varepsilon^2), \end{equation}

in (Equation 3.17). Notably, these dynamics for $A$ and $B$ near onset correspond to strong competition favouring the fissure branch. In other words, near onset, the fissure equilibria with $B=0$ (and the ones with $A=0$) are stable while the $A=B$ bubble equilibria are unstable.

Hexagonal lattice. We proceed in an analogous fashion. We restrict to a function space with reflection symmetry $y\mapsto -y$ and sixfold rotation in the plane. The centre subspace is three-dimensional, parameterized via $u_c = A\cos(x + \frac{y}{\sqrt{3}}) + B\cos(x - \frac{y}{\sqrt{3}}) + C\cos(\frac{2y}{\sqrt{3}})$. On the centre-subspace, symmetries are generated by rotations of $\pi/3$, $A\rightarrow B\rightarrow C\rightarrow A$, reflections in $y$, $A\leftrightarrow B$, and translations in $x$ by $\pi$, $(A,B)\mapsto -(A,B)$ and give the reduced equation, at leading order,

(3.19)\begin{equation} \begin{aligned} A' &= ({\widetilde{\mu}}-{\widetilde{\mu}}_c)A + \alpha BC + A(\beta A^2 + \gamma( B^2+ C^2)) \\ B' &= ({\widetilde{\mu}}-{\widetilde{\mu}}_c)B + \alpha AC +B(\beta B^2 + \gamma ( A^2+C^2)) \\ C' &= ({\widetilde{\mu}}-{\widetilde{\mu}}_c)B+ \alpha AB +C(\beta C^2 + \gamma ( A^2+B^2)). \end{aligned} \end{equation}

Fissure solutions correspond to $B = C = 0$, and bubbles to $|A| = |B| = |C|$, with amplitudes solving

\begin{align*} \ {\widetilde{\mu}}-{\widetilde{\mu}}_c &= -\beta A^2 + \mathcal{O}(A^4) \ \ \textrm{(fissures), } \\ \ {\widetilde{\mu}}-{\widetilde{\mu}}_c &= -\alpha A - (\beta + 2\gamma)A^2 + \mathcal{O}(A^3) \ \ \textrm{(bubbles). } \end{align*}

In order to determine $\beta$, we expand our analytical expression for the fissure coefficient ${\widetilde{\mu}}_1(A) = \frac{2}{\pi^2}\left(\frac{1-\sqrt{1-A^2}}{A^2}\right)$ at $A=0$, and find

\begin{equation*} {\widetilde{\mu}}_1 - {\widetilde{\mu}}_1(0) = \frac{1}{4\sqrt{3}\pi^2}A^2 + \mathcal{O}(A^4), \end{equation*}

thus $\beta = \frac{-1}{4\sqrt{3}\pi^2}\varepsilon + \mathcal{O}(\varepsilon^2)$. Next, to find $\alpha, \gamma$, we expand the analytical expression from (Equation 3.13) for the bubble coefficient ${\widetilde{\mu}}_1(A)$ at $A = 0$ to find

\begin{equation*} {\widetilde{\mu}}_1 - {\widetilde{\mu}}_1(0) = -\frac{1}{2\sqrt{3}\pi^2}A + \frac{5}{4\sqrt{3}\pi^2}A^2 + \mathcal{O}(A^3), \end{equation*}

which, when comparing, gives $\alpha = \frac{1}{2\sqrt{3}\pi^2}\varepsilon + \mathcal{O}(\varepsilon^2)$ and $\gamma = \frac{-1}{2\sqrt{3}\pi^2}\varepsilon + \mathcal{O}(\varepsilon^2)$.

Investigating stability in this reduced ODE, we find that at small amplitudes all solutions are unstable, as is common in the generic bifurcation scenario on hexagonal lattices [Reference Buzano and Golubitsky4]. The triangles stabilize at finite amplitude and parameter value ${\widetilde{\mu}}_- = \frac{-1}{20\sqrt{3}\pi^2}\varepsilon + \mathcal{O}(\varepsilon^2)$ in a saddle-node bifurcation. We compared predictions from centre manifold with the perturbative analysis, also indicating stability of solutions with good agreement away from vacuum formation; see Figure 5.

Figure 5. Centre-manifold predictions vs. numerical values, for $\varepsilon = 10^{-4}$ (square symmetry, left) and $\varepsilon = 2.5 10^{-3}$ (hexagonal symmetry, right). The parameter is shifted relative to the bifurcation point at $\varepsilon=0$. Centre-manifold predictions are in lighter colours. Fissures are in blue; bubbles are in green and black. Black represents branches found to be stable in centre-manifold predictions and in numerical continuation, respectively.

4.1. Numerical continuation: full equation with small diffusion

We implement secant continuation for the potentials in (Equation 1.8) on both square and hexagonal lattices. In the case of the hexagonal lattice, we transform coordinates, mapping the basis lattice vectors $e_1$ and $e_2$ from the hexagonal lattice onto the canonical basis and then computing on a square domain, with suitably transformed differential operators and potential. We impose phase conditions in $x$ and $y$ with reference solution $u^{(0)}$ to exclude the family of translations of the solution (see for instance [Reference Roose and Szalai14]), and a mass condition to select a branch with fixed average density. Mass and positional constraints are compensated for by dummy variables, adding a mass loss term $\alpha u$ and drift terms $s^xu_x$ and $s^yu_y$ to the equation. Numerically, we found $\alpha=s^x=s^y=0$ to numerical precision as expected. When solving for fissure branches, only one of the phase conditions is used and the drift in the $y$-direction is omitted. Together, we solve

(4.3)\begin{equation} \varepsilon \Delta u + \nabla \cdot \left( u \,\nabla \left( u - {\widetilde{\mu}} (V*u)\right) \right) -\alpha u + s^x u_x + s^yu_y = 0 \quad \text{(evolution equation),} \end{equation}

(4.4)\begin{equation} \qquad\qquad\qquad\qquad\quad\quad\quad\!\frac{1}{4\pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi} u(x,y)dydx = 1\quad \text{(mass condition),} \end{equation}

(4.5)\begin{equation} \begin{cases} \int_{-\pi}^\pi\int_{-\pi}^\pi u_x^{(0)}(x,y)(u^{(0)}(x,y)-u(x,y))dx dy = 0 \\ \int_{-\pi}^\pi\int_{-\pi}^\pi u_y^{(0)}(x,y)(u^{(0)}(x,y)-u(x,y))dx dy = 0 \end{cases} \text{(phase conditions),} \end{equation}

in the case of square lattice periodicity, with suitably transformed differential operators in the case of hexagonal lattice periodicity. For continuation, we add the usual secant condition and set $u^{(0)}$ as the reference solution from the previous step. The solution $u$ is discretized to $\underline{u}$ using $N^2$ Fourier modes in the two-dimensional domain. A Newton method is implemented in Matlab using its gmres linear solver for the $N^2 + 4$ variables $(\underline{u},s^x,s^y,\alpha,{\widetilde{\mu}})$ (resp. $N^2 +3$ variables for fissures). The diffusion term $\varepsilon$ serves as a regularizer and preconditioner for the linear solver.

Results are displayed in Figures 6 and 7. Figure 6 gives a quantitative comparison to the theory, displaying convergence of numerically computed branches to the explicit prediction. As predicted, the square of the vacuum area for bubbles grows linearly in the parameter, with slope matching our explicit calculation for the growth of vacuum bubbles, in both square and hexagonal periodicity (here for the triangle branch). Convergence in $\varepsilon$ is linear as predicted. Different branches, for both the square lattice and the hexagonal lattice periodicity are shown in Figure 7. Calculations include stability information which agrees with the information obtained from the centre-manifold analysis.

Figure 6. Numerical continuation of branches forming vacuum bubbles in case of square lattice periodicity (left) and hexagonal lattice periodicity (right, triangle branch), for different values of $\varepsilon$. Plotted is the square of the size of the vacuum area, computed as the area where $u \lt \varepsilon$, together with the theoretical predictions from Theorems 1 and 2, respectively.

Figure 7. Bifurcation diagrams from secant continuation for the full equation with diffusion, $\varepsilon = 0.05$ (square symmetry, left), $\varepsilon = 0.01$ (hexagonal symmetry, right). The colour along the fissure and bubble branches represents the size of the largest numerically computed eigenvalue of the Jacobian at that equilibrium, so that lighter branches are more unstable. Black represents a stable branch, where the largest computed eigenvalue is 0 due to translational symmetry.

4.2. Numerical continuation: finite rank system, no diffusion

In order to continue solutions of (Equation 4.1), we use the somewhat explicit form of solutions as in §2, and numerically solve Equation 2.4–(Equation 2.5) for the coefficients of Fourier modes.

Square symmetry. In order to find bubble and fissure solutions, we set

\begin{equation*} u =\big(A_0 + A_1\cos(x)\cos(y) + B_1\sin( x)\sin(y)\big)_+, \end{equation*}

where $f_+$ denotes the positive part of $f$. Bubble branches correspond to $B_1 = 0$, while fissure branches will have $A_1 = B_1 \neq 0$. We then numerically solve the system

(4.6)\begin{equation} \begin{aligned} 0& = 2\pi^2 - \iint_{\Omega_0} \left( A_0 + A_1 \cos(x) \cos(y)+ B_1\sin(x)\sin(y) \right) d x d y \\ 0& = A_1 - \left(\frac{1}{\pi^2} + {\widetilde{\mu}}\right)\iint_{\Omega_0} 2\cos(x) \cos(y)\left( A_0 + A_1 \cos(x) \cos(y) \right) d x d y\\ 0& = B_1 - \left(\frac{1}{\pi^2} + {\widetilde{\mu}}\right)\iint_{\Omega_0} 2B_1\sin^2(x) \sin^2(y) d x dy, \end{aligned} \end{equation}

where the domain $\Omega_0 $ is the subset of the domain on which $u$ is positive.

Hexagonal symmetry. For the case of hexagonal symmetry, we make the ansatz

\begin{equation*} u = \big(A_0 + A_1\cos(x-\frac{y}{\sqrt{3}}) + B_1\cos(x+\frac{y}{\sqrt{3}})+C_1\cos(\frac{2y}{\sqrt{3}})\big)_+. \end{equation*}

Here, a bubble branch satisfies $A_1 = B_1 = C_1 \neq 0$, while a fissure branch has $A_1 \neq 0, B_1 = C_1 = 0$. We solve the system

(4.7)\begin{equation} \begin{aligned} 0& = 2\sqrt{3}\pi^2 - \iint_{\Omega_0} \left( A_0 + A_1\cos(x-\frac{y}{\sqrt{3}}) + B_1\cos\bigg(x+\frac{y}{\sqrt{3}}\bigg) +C_1\cos\bigg(\frac{2y}{\sqrt{3}}\bigg) \right) dxdy \\ 0& = A_1 - \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)\iint_{\Omega_0} \cos\bigg(x-\frac{y}{\sqrt{3}}\bigg)\left( A_0 + A_1\cos\bigg(x-\frac{y}{\sqrt{3}}\bigg)\right) dxdy\\ 0& = B_1 - \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)\iint_{\Omega_0}\cos\bigg(x+\frac{y}{\sqrt{3}}\bigg)\left( A_0 + B_1\cos\bigg(x+\frac{y}{\sqrt{3}}\bigg)\right)dxdy\\ 0& = C_1 - \left(\frac{1}{\sqrt{3}\pi^2} + {\widetilde{\mu}}\right)\iint_{\Omega_0}\cos\bigg(\frac{2y}{\sqrt{3}}\bigg)\left( A_0 + C_1\cos\bigg(\frac{2y}{\sqrt{3}}\bigg)\right)dxdy, \end{aligned} \end{equation}

numerically, where $\Omega_0 $ is the subset of the domain on which $u$ is positive. Integrals are evaluated numerically in case of both square and hexagonal lattice periodicity, using a grid of $N^2$ grid points, $N = 3000$.

4.3. Branches and stability

In the case of square lattice symmetry, numerical continuation, both with and without diffusion, reveals a secondary bifurcation where the solution branches exchange stability. Such a bifurcation is accompanied by a branch of mixed-mode solutions connecting the bubble and fissure branches. In the case of no diffusion, this occurs at exactly ${\mu} = 2{\widetilde{\mu}}_*$. In this case, the branch is again vertical, and its existence is explicit; see §4.4. With diffusion, similarly to the primary vertical branch, this branch is almost vertical and slightly shifted in parameter space. The mixed-mode branch consists of solutions with elliptical bubbles that interpolate between bubbles and fissures. Modulations of fissures or more generally striped patterns have been observed in many other pattern-forming scenarios, for instance in the Taylor–Couette problem, where they are known as wavy rolls [Reference Chossat and Iooss8]. Again, the situation here is particular as the exchange of stability between fissures and bubbles does not display hysteresis at $\varepsilon=0$.

In addition to being curiously vertical, this secondary bifurcation also occurs precisely at the parameter value where the bubble branch changes topology, that is, where vacuum bubbles merge and confine the support of the solution $\Omega_0$ to disconnected clusters. In the finite-rank case, this change of topology happens precisely when $A_0 = 0$. At this secondary bifurcation, stability at both fissures and bubbles changes from what emerged from the local bifurcation and the primary vertical branch. Prior to the transition, fissure solutions are stable; afterwards, clusters (the continuation of bubbles) are stable. For square symmetry, then, the only stable solutions are fissures, first, and clusters, later; solutions with true vacuum bubbles are unstable; see Figure 7 for the bifurcation diagram with numerically computed eigenvalues.

We did not find such a secondary bifurcation for hexagonal symmetry. In this case, the only stable solutions are triangles, both near the bifurcation point and for larger values of ${\widetilde{\mu}}$; fissure and hexagon solutions are always unstable. This was confirmed both in the numerical continuation of the PDE solution with small diffusion and in the finite-rank continuation: the computed Morse index remains constant along all three solution branches for ${{\widetilde{\mu}}} \gt 0$.

4.4. Vertical branch II: the connecting branch

For (Equation 4.1), without diffusion and with square symmetry, we find an explicit vertical mixed-mode branch at ${\mu} = 2{\widetilde{\mu}}_*$. These solutions all have $A_0 = 0$.

Proposition 3. Let ${\widetilde{\mu}} = \frac{2}{\pi^2}, V (x, y) = 2\cos(x) \cos(y)$. For each $\tau \in [0,1]$, there exists a solution $u$ to Equation 2.4–(Equation 2.6) of the form

\begin{equation*} u(x,y) = \left(A_1 \cos (x) \cos (y) + \tau A_1\sin (x)\sin (y) \right)_+, \end{equation*}

for a unique $A_1=A_1(\tau) \gt 0$.

Proof. Let $\Omega_0' = \left\{(x,y) \in [-\pi,\pi^2] \ \big| \ A_1 \cos(x) \cos(y)+ \tau \ A_1\sin(x)\sin(y) \gt 0\ \right\}$. Substituting the ansatz ${\widetilde{\mu}} = \frac{2}{\pi^2},\ \ u = \left(A_1 \cos (x) \cos (y) + \tau A_1\sin (x)\sin (y) \right)_+$ into (Equation 2.4)–(Equation 2.6), we obtain the system

(4.8)\begin{equation} 0 = 4\pi^2 - \iint_{\Omega_0'} \left( A_1 \cos(x) \cos(y)+ \tau A_1\sin(x)\sin(y) \right) d x d y \end{equation}

(4.9)\begin{equation} 0 = A_1 - \frac{2}{\pi^2} \iint_{\Omega_0'} A_1 \cos^2(x) \cos^2(y) d x d y\qquad\qquad\qquad\quad \end{equation}

(4.10)\begin{equation} 0 = \tau A_1 - \frac{2}{\pi^2} \iint_{\Omega_0'} \tau A_1\sin^2(x) \sin^2(y) d x dy.\qquad\qquad\qquad \end{equation}

We note two key facts:

1. (1) the complement $(\Omega_0')^c$ is exactly equal to $\Omega_0'$, shifted by $\pi$ in the $x$-direction (or, equivalently, by $\pi$ in the $y$-direction). This follows immediately from the symmetries $\cos(x + \pi) = -\cos(x),$ and $ \sin(x+\pi) = -\sin(x)$;

2. (2) the functions $\cos^2(x) \cos^2(y)$ and $\sin^2(x) \sin^2(y)$ are both $\pi$-periodic in $x$.

These facts imply that

\begin{align*} \iint_{\Omega_0'} \cos^2(x) \cos^2(y) d x d y &= \frac{1}{2}\int_{-\pi}^\pi \int_{-\pi}^\pi \cos^2(x) \cos^2(y) d x d y, \\ \iint_{\Omega_0'} \sin^2(x) \sin^2(y) d x d y &= \frac{1}{2}\int_{-\pi}^\pi \int_{-\pi}^\pi \sin^2(x) \sin^2(y) d x d y. \end{align*}

Then (Equation 4.9) and (Equation 4.10) are automatically satisfied no matter the value of $A_1, \tau$. Considering (Equation 4.8), we note that for fixed $\tau$, the regions $\Omega_0, \Omega_0'$ do not depend on $A_1$. Then the quantity

\begin{equation*}\iint_{\Omega_0'} \left( \cos(x) \cos(y)+ \tau\sin(x)\sin(y) \right) d x d y \end{equation*}

is simply a number, independent of $A_1$, and we can solve as

\begin{equation*} A_1(\tau) = \frac{4\pi^2}{\iint_{\Omega_0'} \left( \cos(x) \cos(y)+ \tau\sin(x)\sin(y) \right) d x d y }. \end{equation*}

Along the bubble branch, we note that such a solution with $A_0 = 0$ corresponds exactly to the transition between vacuum bubbles and isolated clusters.

Figure 8. Bifurcation diagrams from finite rank continuation, with corresponding states along the branches, showing square symmetry (left) and hexagonal symmetry (right). In the left panel, one sees in purple the interpolating branch of mixed-mode solutions. For hexagonal symmetry, we find no such branch. The computed Morse index remains constant along all three solution branches, with no observed crossing which would indicate a bifurcation.

5.1. Numerical continuation

We implement secant continuation to find parameter values where non-trivial branches bifurcate from the crystalline state in (Equation 5.1) and track branches of steady states. The Bessel potential is of course not periodic and needs to be replaced by its periodization, that is, the sum over all lattice translates. Since this sum is not explicit, we approximate the periodized potential by the sum of nine translated copies of $K^\zeta$. The translates are centred at all combinations of the respective lattice’s basis vectors around and including the origin (i.e. $(0,0), \pm\mathbf{e}_1, \pm \mathbf{e}_2$, etc.).

Adding phase conditions in both $x$- and $y$-directions, similarly to §4, we solve

\begin{align*} 0 &= -\frac{1}{N^2} \sum_{\underline{j} \neq \underline{\ell}}^{N} \partial_{1} W(x_{\underline{\ell}}-x_{\underline{j}}, y_{\underline{\ell}}-y_{\underline{j}} ) + s^x\mathbf{1},\\ 0&= -\frac{1}{N^2} \sum_{\underline{j} \neq \underline{\ell}}^{N} \partial_{2} W(x_{\underline{\ell}}-x_{\underline{j}}, y_{\underline{\ell}}-y_{\underline{j}} )+ s^y\mathbf{1},&&\textrm{(evolution equations)} \\ 0 &= \frac{1}{N^2}\sum_{\underline{j}} x_{\underline{j}}-\overline{x}_{\underline{j}}, \\ 0 &= \frac{1}{N^2}\sum_{\underline{j}} y_{\underline{j}}-\overline{y}_{\underline{j}}, &&\textrm{(phase conditions)} \end{align*}

where $\mathbf{1} = (1, 1, ... 1)$, and $s^x,s^y$ are dummy variables associated with the phase conditions. There is clearly no mass constraint necessary here. We solve the system in $2N^2+3$ variables $(x,y,,s^x,s^y,{\widetilde{\mu}})$ using Newton’s method and a secant normalization. All results described in the following were obtained with a fixed $\zeta = 0.3$.

Figure 10 shows the numerically calculated onset of instability of the crystalline state depending on the particle number, compared to the prediction from the continuum limit. We find convergence proportional to the inverse of the total particle number in both square and hexagonal lattice configurations.

Figure 10. Onset of instability of the pure crystalline state found using secant continuation for fixed $\zeta=0.3$ and varying total numbers of particles $N^2$ on the square lattice (left) and the hexagonal lattice (right). The finite- $N$ corrections to the bifurcation point appear to be $\mathcal{O}(N^{-2})$.

Figure 11 shows a numerically computed bifurcation diagram in the case of the hexagonal lattice. We find triangle, hexagonal and fissure branches. As predicted in the continuum limit, triangle branches are stable after an initial saddle-node bifurcation, while hexagon branches are unstable. Fissure branches are stable only for subcritical parameter values. The discreteness notably causes the bifurcation of fissures to be weakly subcritical, as opposed to the weak supercriticality caused by small noise $\varepsilon$ in §3, a phenomenon also observed in the one-dimensional case [Reference Scheel and Stevens15]. We noted that during continuation, there are several secondary bifurcations, and we do not claim at all that this bifurcation diagram is complete. An intuitive reason for secondary bifurcations are possible rearrangements and emergence of defects in the particle configuration, an effect which is notably absent in the one-dimensional case [Reference Scheel and Stevens15].

Figure 11. Bifurcation diagram obtained using secant continuation on the hexagonal lattice for $N^2= 49$ particles, with sample particle configurations observed in direct simulations. Triangle branch and hexagon branch are in green and dark green, respectively, fissure branch is in blue. Hexagons are not observed in direct simulations. On the vertical axis, we show the inverse of minimal distances between particles as a proxy for the maximal density.

In fact, such rearrangements and secondary bifurcations prevented us from creating the analogous figure in the case of a square lattice using numerical continuation. We therefore now turn to direct simulations, which do, in the case of square lattice symmetry, show a striking similarity to our results in the continuum limit.

5.2. Direct simulations and observed states

We collect here some observations from direct simulations.

Square lattice. On the square lattice, as in the continuum, we find that for a range of $\mu$ values above onset, perturbations from the crystal state lead to particles arranging to form fissures. Past a secondary bifurcation point, particles organize only transiently into fissures, which give way into single clusters. Figure 12 shows the states observed after long times for various $\mu$ values. In all cases observed, although the crystal state has a square particle microstructure, the vacuum states exhibit a hexagonal lattice microstructure (with occasional defects as expected in finite-particle simulations). In other words, in the supercritical pattern-forming regime, the predictions from the macrostructure of the problem appear to hold independently of particle microstructure.

Figure 12. Bifurcation on the square lattice, $N^2 = 144$ particles. States observed at $t = 10000$ in direct simulations for various values of ${\mu}$, from left to right: ${\mu}$ = 0.075, 0.1, 0.125, 0.15, 0.16, 0.2. Initial conditions were given by perturbing the crystal state with random noise at $t=0$, in all cases except ${\mu} = 0.075$. In that case, the initial perturbation was in the direction of the unstable eigenvector $\cos(\overline{x})\cos(\overline{y})$, due to the fact that random perturbations favoured the instability of the square lattice toward a uniform state with a hexagonal lattice microstructure.

Square lattice: hysteresis near bifurcation points. Near both the initial and secondary bifurcation points, the discrete structure induces hysteresis, more pronounced for small numbers of particles. In the subcritical region where the branches bend backwards, both fissures and bubbles can be seen at the same values of ${\mu}$, depending on whether the system evolves from a fissure or cluster state; see Figure 13. In each case, the selected state appears stable to small perturbations in the direction of the other. We note that for square symmetry of the potential, neither of the continuum models (the original system (Equation 2.4)–(Equation 2.6) or the corresponding system with diffusion) has a backward-bending branch – this form of hysteresis is an intrinsically discrete phenomenon.

Figure 13. Hysteresis near initial bifurcation in discrete model (square lattice), $N^2 = 324$ particles. Bubble and fissure states can be observed at the same ${\mu}$ values in the subcritical region where the branches bend backwards (second and fourth panels). The states in the second and fourth panels were obtained by first increasing ${\mu}$ to produce the states seen in the first and third panels, and then decreasing ${\mu}$ to 0.078, which is in the subcritical region for $N^2 = 324$.

Similarly, for a range of ${\mu}$ near the secondary bifurcation point where fissures and bubbles exchange stability, both fissures and clusters can be seen at the same values of ${\mu}$, depending on initial conditions; see Figure 14.

Figure 14. Hysteresis near secondary bifurcation in discrete model (square lattice), $N^2 = 324$ particles. Both bubble and fissure states can be observed after long times for a range of ${\mu}$ values near the secondary bifurcation, depending on the initial condition. Left two panels: ${\mu} = 0.155$. Right two panels: ${\mu} = 0.175$. For $N^2 = 324$, the interval $(0.155,0.177)$ appears to be the approximate range of bistability.

Hexagonal lattice. On the hexagonal lattice, in the supercritical range of ${\mu}$, only clusters are observed. Fissures and hexagons are not seen in direct simulations. This corroborates findings from §3 and §4 concerning stability.

All states observed in the supercritical regime were single connected clusters; true bubble solutions were seen only in the subcritical backward-bending part of the branch. It is worthwhile to note, however, that in the continuum model, the topology change from bubbles to clusters happens at quite a small value ${\widetilde{\mu}} \approx 0.003$. It is therefore likely that a very high number of particles, larger than the maximal number $N^2\sim 400$ considered here, would be needed to observe true bubble solutions supercritically. In the subcritical region, similarly to the square lattice, both fissures and bubbles can be observed, depending on initial conditions.

Figure 15 shows a hysteresis loop on the hexagonal lattice, where ${\mu}$ is initially increased past the bifurcation point, then decreased, so that the solution travels back along the upper bifurcation branch. Once the saddle-node is reached, when ${\mu}$ is decreased just a bit more, the solution‘falls off’ this branch and returns to a crystal-like formation.

Figure 15. Hysteresis loop on the hexagonal lattice, $N^2 = 49$ particles.