remarkRemark hypothesisHypothesis claimClaim Generation and motion of
interfacesP. W. Miller et al.

# Generation and motion of interfaces in a mass-conserving reaction-diffusion
system

Pearson W. Miller  1 Daniel Fortunato  2 Matteo Novaga3
Stanislav Y. Shvartsman1,4,5 Cyrill B. Muratov3,6,†
1Center for Computational Biology, Flatiron Institute, New York, NY 10010, USA
2Center for Computational Mathematics, Flatiron Institute, New York, NY 10010,
USA
3Dipartimento di Matematica, Università di Pisa, Largo B. Pontecorvo 5, 56127
Pisa, Italy
4Department of Molecular Biology, Princeton University, Princeton, NJ 08540,
USA
5Lewis-Sigler Institute for Integrative Genomics, Princeton University,
Princeton, NJ 08540, USA
6Department of Mathematical Sciences, New Jersey Institute of Technology,
Newark, NJ 07102, USA
†To whom the correspondence should be addressed
E-mail<EMAIL_ADDRESS>
###### keywords:

pattern formation, reaction-diffusion, singular perturbations, Laplace-
Beltrami operator, long-time behavior

> Reaction-diffusion models with nonlocal constraints naturally arise as
> limiting cases of coupled bulk-surface models of intracellular signalling.
> In this paper, a minimal, mass-conserving model of cell-polarization on a
> curved membrane is analyzed in the limit of slow surface diffusion. Using
> the tools of formal asymptotics and calculus of variations, we study the
> characteristic wave-pinning behavior of this system on three dynamical
> timescales. On the short timescale, generation of an interface separating
> high- and low-concentration domains is established under suitable
> conditions. Intermediate timescale dynamics is shown to lead to a uniform
> growth or shrinking of these domains to sizes which are fixed by global
> parameters. Finally, the long time dynamics reduces to area-preserving
> geodesic curvature flow that may lead to multi-interface steady state
> solutions. These results provide a foundation for studying cell polarization
> and related phenomena in biologically relevant geometries.

35Q92, 35K57, 92C37

##  1 Introduction

Reaction-diffusion models are essential tools for understanding the spatial
self-organization of chemical patterns inside the cell. Over the past decade,
research interest has coalesced around two key features which set
intracellular dynamics apart from other classes of models. First, in contrast
to traditional models which exist on a single domain, these models are often
bulk-surface models, in that they feature distinct diffusion processes within
the 3D cytosolic volume of the cell and on the 2D cell membrane coupled by a
nonlinear boundary condition [46]. Second, many recent studies have emphasized
systems subject to mass-conservation [42, 24]. These properties are
particularly prevalent in models of cell polarization, a crucial process by
which the spatial distribution of proteins within a cell becomes highly
localized to a region of the membrane as a result of spontaneous symmetry
breaking or an external guiding cue [16]. Polarization is essential for a
great variety of biological phenomena, including guiding developmental
outcomes, establishing axes for cell division and guiding locomotion in motile
cells, and so has motivated significant interest from the mathematical biology
community [16].

The attention of applied mathematicians is increasingly turning to
understanding how various forms of spatial heterogeneity influence dynamics
and steady state behavior. Research that has emphasized the effect of
spatially varying kinetic parameters are the most obvious examples of this
trend [41, 21]. Over the last few years, the question of what role cell
geometry plays in guiding polarization has drawn considerable interest. While
initial polarization models were reduced to 1D systems, computational advances
have enabled numerical studies on 2D and fully 3D domains [15, 6]. Early steps
in addressing this problem have produced simulations that are highly
suggestive that localization is closely tied to surface curvature via
minimization of interfacial length, but this principle has yet to be
demonstrated by formal analysis [14, 20]. Very recently, a novel numerical
framework was introduced allowing for efficient simulation of cell
polarization models on surfaces of revolution, which in particular motivates
this particular study [37].

In this paper, we perform a formal asymptotic analysis of the surface-bound
version of the wave-pinning model first proposed in Ref. [39], which features
both mass-conservation and bulk-surface coupling. This model has attracted
considerable interest as a minimal theory for polarization, and has been
successfully applied to model systems such as the Rho-GTP pathway [56] or
Ezrin polarization in embryonic mouse cells [58]. The emergence of spontaneous
symmetry breaking in this system has been intensely studied, and the
dependence of the initial instability of the spatially uniform solution on
various parameters is well-characterized [56, 33, 48, 43]. While previous
studies considered the asymptotics of wave-pinning as surface diffusion
becomes small, they have generally been concerned with 1D domains or lower
orders of perturbation than we discover are needed to capture the full effect
of 3D-embedded domain geometry [40, 14]. Other works which have considered
fully 3D domains have limited analysis to linear stability studies [47, 20].
Our research here uses the timescale separation techniques to probe the
generation and propagation of interfaces in two-species reaction-diffusion
systems to the specific case of the wave-pinning model (for related rigorous
studies, see [27, 51, 12]). We note that for single bistable reaction-
diffusion equations in the Euclidean setting the generation and propagation of
interfaces is by now well understood mathematically [11, 5, 1]. In the context
closely related to our problem, results on propagation of interfaces for
bistable reaction-diffusion equations on Riemannian manifolds were recently
obtained in Refs. [45, 44].

In this work, we demonstrate that understanding the patterning outcome of
polar domains on a surface of arbitrary shape requires a thorough analysis of
the long-timescale behavior of the associated mass-conserving reaction-
diffusion equation. Working in the limit of slow surface diffusion, we
establish a separation of the dynamics into three distinct time scales: first,
the initial generation of the interface, then, uniform growth or shrinking of
the domains until their areas converge to the steady state values dependent
only on the global parameters, and finally, evolution of the interface by
area-preserving geodesic curvature flow that results in a finite union of
geodesic disks as time goes to infinity. In doing so, we demonstrate that
stable steady states with multiple disjoint interfaces exist on biologically
plausible domain shapes, something impossible without geometric effects.

This paper is organized as follows. In Sec. 2, we formulate our model of
polarization on a closed surface and define three sub-problems whose limiting
behavior should describe the effective dynamics on different asymptotic
timescales. In Sec. 3, we carry out a preliminary analysis of existence and
stability of the uniform states and, in particular, identify the parameter
regimes in which only nonuniform steady states can be stable. In Secs. 4
through 6, we asymptotically derive the limiting sub-problems as the surface
diffusion coefficient tends to zero and examine their dynamical behavior.
Finally, in Sec. 7 we perform a number of numerical tests to corroborate the
predictions of the asymptotic theory, and in Sec. 8 we make our concluding
remarks.

##  2 Model Summary

Let $\Omega\subset\mathbb{R}^{3}$ be an open, bounded, connected set with a
sufficiently regular boundary $\partial\Omega$. A basic model of cell
polarization can be written as (see [15] for more details):

(1) $\displaystyle\partial_{t}B$
$\displaystyle=D_{B}\nabla^{2}_{\partial\Omega}B+k_{b}\left(\beta+\frac{B^{\nu}}{G^{\nu}+B^{\nu}}\right)C-k_{d}B$
$\displaystyle\text{in }\partial\Omega\times(0,T),$ (2)
$\displaystyle\partial_{t}C$ $\displaystyle=D_{C}\nabla^{2}_{\Omega}C,$
$\displaystyle\text{in }\Omega\times(0,T),$ (3) $\displaystyle D_{C}(\nabla
C\cdot\hat{n})|_{\partial\Omega}$
$\displaystyle=-k_{b}\left(\beta+\frac{B^{\nu}}{G^{\nu}+B^{\nu}}\right)C+k_{d}B.$
$\displaystyle\text{in }\partial\Omega\times(0,T),$ (4) $\displaystyle
B(\cdot,0)$ $\displaystyle=B_{0}$ $\displaystyle\text{ in }\partial\Omega,$
(5) $\displaystyle\ C(\cdot,0)$ $\displaystyle=C_{0}$ $\displaystyle\text{ in
}\Omega.$

Here, (1) describes a reaction-diffusion process on the surface
$\partial\Omega$ of the surface-bound protein concentration $B$, while (2)
gives bulk diffusion of the concentration $C$ of the same protein in the
bounded volume $\Omega$, and (3) is the boundary condition coupling the two.
The operator $\nabla^{2}_{\partial\Omega}$ is the Laplace–Beltrami operator on
the surface, while $\nabla^{2}_{\Omega}$ refers to the standard 3D Laplacian.
While we start with a model where membrane-bound dynamics is already purely
two-dimensional, the validity of this class of model as a limiting case of a
membrane of finite thickness was demonstrated in Ref. [31].

We are principally interested in the regime where $D_{C}\gg D_{B}$, and in
this regime it is reasonable to treat the bulk concentration as spatially
uniform [15]. As mass is conserved globally, we define the quantity
$C_{\text{tot}}$ as the total number of protein molecules divided by the bulk
volume $V$ of $\Omega$, so that the bulk concentration can be expressed as

(6) $C=C_{\text{tot}}-\frac{1}{V}\int_{\partial\Omega}B\,dS.$

Substituting this back into (1) yields a non-local equation for surface-bound
species

(7)
$\displaystyle\partial_{t}B=D_{B}\nabla^{2}_{\partial\Omega}B-k_{d}B+k_{b}\left(\beta+\frac{B^{\nu}}{G^{\nu}+B^{\nu}}\right)\left(C_{\text{tot}}-\frac{1}{V}\int_{\partial\Omega}B\,dS\right)$
$\displaystyle\text{ on }\partial\Omega\times(0,T).$

This equation is rendered dimensionless by a rescaling $B(x,t)\rightarrow
u(\tilde{x},\tilde{t})$, with $u=k_{d}B/(k_{b}C_{\text{tot}})$,
$\tilde{t}=k_{d}t$ and $\tilde{x}=x/\sqrt{A}\in\partial\tilde{\Omega}$, where
$A=|\partial\Omega|$ is the surface area of $\partial\Omega$, and
$\tilde{\Omega}=\Omega/\sqrt{A}$ is a rescaling of $\Omega$ that ensures that
the rescaled domain $\tilde{\Omega}$ has boundary of unit area. Dropping the
tildes and the subscript $\partial\Omega$ for simplicity of notation from now
on, we arrive at the dimensionless form

(8) $\displaystyle\partial_{t}u$
$\displaystyle=\delta^{2}\nabla^{2}u-u+f(u)\left(1-\alpha U\right)$
$\displaystyle\text{ on }\partial\Omega\times(0,T),$ (9) $\displaystyle U(t)$
$\displaystyle=\int_{\partial\Omega}u(x,t)dS$ $\displaystyle\text{for
}t\in(0,T),$ (10) $\displaystyle u(x,0)$ $\displaystyle=g(x)$
$\displaystyle\text{for }x\in\partial\Omega,$

with the dimensionless parameters111We use a slightly different convention
from [15], in which the definitions of $\alpha$ and $\delta$ differ from the
present ones by constant factors.

(11) $\displaystyle\delta=\sqrt{D_{B}\over k_{d}A},\qquad\alpha={k_{b}A\over
k_{d}V},\qquad\gamma={k_{d}G\over k_{b}C_{\text{tot}}},$

and where we defined $f(u)=\beta+\frac{u^{\nu}}{\gamma^{\nu}+u^{\nu}}$ for
further notational convenience. Together with the already dimensionless
parameters $\beta$ and $\nu$, the parameters $\alpha$, $\gamma$ and $\delta$
define the parameter space for our problem. As was already noted, we have
$|\partial\Omega|=1$ now.

The long timescale dynamics of this problem for a purely 1D model was
previously studied via asymptotic expansions in [40], and later on a disc in
[14]. In Ref. [15], exact solutions were constructed for the special case of a
spherical domain and $\nu=\infty$. Here we confine ourselves to the specific
case of $\nu=2$ for the sake of analytical tractability, but treat general
spatial domains. A few useful remarks on this system can be made based upon
existing results. First, (1)–(5) represents a special case of a more general
class described in Ref. [52], and per the results therein, there exists a
unique classical solution $(B,C)$ such that the functions $B$ and $C$ are
smooth and uniformly bounded. Further, Theorem $2.6$ of Ref. [26] demonstrates
that our shadow system (8)–(10) likewise has a unique weak solution, and that
solution is the limit of the solutions to (1)–(5) as $D_{B}\rightarrow\infty$,
justifying the use of this reduced approach.

In the asymptotic analysis described in the following sections, we find it
convenient to frame our problem as three equivalent problems parametrized by
$\delta\ll 1$:

$\displaystyle(P_{0}^{\delta})$
$\displaystyle\qquad\partial_{t}u_{0}^{\delta}=\delta^{2}\nabla^{2}u_{0}^{\delta}-u_{0}^{\delta}+f(u_{0}^{\delta})\left(1-\alpha\int_{\partial\Omega}u_{0}^{\delta}dS\right),\qquad
u_{0}^{\delta}(x,0)=g_{0}^{\delta}(x),$ $\displaystyle(P_{1}^{\delta})$
$\displaystyle\qquad\partial_{t}u_{1}^{\delta}=\delta\nabla^{2}u_{1}^{\delta}+\delta^{-1}\left[-u_{1}^{\delta}+f(u_{1}^{\delta})\left(1-\alpha\int_{\partial\Omega}u_{1}^{\delta}dS\right)\right],\qquad
u_{1}^{\delta}(x,0)=g_{1}^{\delta}(x),$ $\displaystyle(P_{2}^{\delta})$
$\displaystyle\qquad\partial_{t}u_{2}^{\delta}=\nabla^{2}u_{2}^{\delta}+\delta^{-2}\left[-u_{2}^{\delta}+f(u_{2}^{\delta})\left(1-\alpha\int_{\partial\Omega}u_{2}^{\delta}dS\right)\right],\qquad
u_{2}^{\delta}(x,0)=g_{2}^{\delta}(x).$

Each problem describes the dynamics at the timescale of a different order in
$\delta$: $(P_{0}^{\delta})$ is the problem on the original $O(1)$ timescale,
$(P_{1}^{\delta})$ is the problem on the $O(\delta^{-1})$ long timescale and
$(P_{2}^{\delta})$ is the problem on the $O(\delta^{-2})$ longer timescale.
Crucially, we will show that as $\delta\to 0$ the initial condition of the
second and third problem may be taken to be the infinite time limit of the
solution at the previous stage, thus allowing to connect the solutions at
different timescales. In our analysis below, we will formally obtain the limit
behavior of each of these problems when $\delta\rightarrow 0$.

##  3 Preliminaries

In this section we introduce some basic facts about our system, and establish
some notation. We begin by examining the behavior of the spatially uniform
steady states of (8). These states satisfy

(12) $-u+f(u)(1-\alpha u)=0,$

and with our choice of the nonlinearity (recall that $\nu=2$) this reduces to
a cubic equation

(13) $\displaystyle-u(u^{2}+\gamma^{2})+(1-\alpha
u)(u^{2}(1+\beta)+\beta\gamma^{2})=0.$

For $\alpha,\beta,\gamma>0$, (13) will have at least one real positive root
$u_{0}$, with two additional positive roots $u_{\pm}$ emerging when its
discriminant

(14)
$\displaystyle\begin{split}\Delta(\alpha,\beta,\gamma)=\gamma^{4}(1&-4\beta(\alpha\beta+\alpha+1)(2\beta(\alpha\beta+\alpha+1)-5))\\\
&-4\gamma^{6}(\alpha\beta+1)^{3}(\alpha\beta+\alpha+1)-4\beta(\beta+1)^{3}\gamma^{2}>0.\end{split}$

We adopt the convention $u_{-}<u_{0}<u_{+}$, with $u_{\pm,0}$ denoting any one
of these roots for shorthand. Setting $\Delta=0$ allows one to derive the
conditions for the existence of three roots. We find the requirements for
three real roots to be

(15) $\displaystyle 0$ $\displaystyle<\beta<\frac{1}{8},$ (16) $\displaystyle
0$ $\displaystyle<\alpha<\frac{1-8\beta}{8\beta(1+\beta)},$ (17)
$\displaystyle\gamma_{-}(\alpha,\beta)$
$\displaystyle<\gamma<\gamma_{+}(\alpha,\beta),$

where

(18)
$\begin{split}\gamma_{\pm}(\alpha,\beta)=\sqrt{\frac{1-4\beta(\alpha\beta+\alpha+1)(2\beta(\alpha\beta+\alpha+1)-5)\pm\sqrt{(1-8\beta(\alpha\beta+\alpha+1))^{3}}}{8(\alpha\beta+1)^{3}(\alpha\beta+\alpha+1)}}.\end{split}$

Slices of the boundary between parameter space regions are plotted in the
$\gamma-\alpha$ plane in Fig. 1 for various values of $\beta$.

Figure 1: Boundaries between the bistable and monostable regimes plotted for
different values of $\beta$. Within the triangular-shaped regions, three
positive roots of (13) exist, while outside there is only one.

We next examine the linear stability of the uniformly stable states in the
usual manner, examining the behavior of perturbations of the form
$u_{\pm,0}(x,t)=u_{\pm,0}+\epsilon\eta_{\pm,0}(x,t)$ for $\epsilon\ll 1$.
Denoting the orthonormal eigenfunctions of the Laplace-Beltrami operator as
$v_{k}(x)$ such that $-\nabla^{2}v_{k}(x)=\lambda_{k}v_{k}(x)$ for
$k=0,1,...$, where $\lambda_{k}$ are the corresponding eigenvalues, we write
$\eta_{\pm,0}(x,t)=\sum_{k=0}^{\infty}a_{k}^{\pm,0}e^{\sigma_{k}^{\pm,0}t}v_{k}(x)$.
Note that for a compact surface, the eigenvalues $\lambda_{k}$ are real with
$0=\lambda_{0}<\lambda_{1}\leq\lambda_{2}\leq...$.

The linearized equation about the root $u^{0,\pm}$ takes the form

$\displaystyle\mathcal{L}\eta_{\pm,0}(x,t)$
$\displaystyle=\left[\partial_{t}-\delta^{2}\nabla^{2}+1-f^{\prime}(u_{\pm,0})(1-\alpha
u_{\pm,0})\right]\eta_{\pm,0}(x,t)+\alpha
f(u_{\pm,0})\int_{\partial\Omega}\eta_{\pm,0}(x,t)dS$ (19)
$\displaystyle=\sum_{k=0}^{\infty}\left(\sigma_{k}^{\pm,0}+\delta^{2}\lambda_{k}+1-f^{\prime}(u_{\pm,0})(1-\alpha
u_{\pm,0})+\alpha
f(u_{\pm,0})\delta_{k,0}\right)a_{k}^{\pm,0}e^{\sigma_{k}^{\pm,0}t}v_{k}(x)=0,$

where $\delta_{k,0}$ is the Kronecker delta symbol, and we took into account
that $v_{0}=1$ and $\int_{\partial\Omega}v_{k}dS=0$ for all $k>0$. Then the
following dispersion relation holds:

(20)
$\sigma_{k}^{\pm,0}=-\delta^{2}\lambda_{k}-1+f^{\prime}(u_{\pm,0})(1-\alpha
u_{\pm,0})-\alpha f(u_{\pm,0})\delta_{k,0}.$

In the three-root parameter region, $\sigma_{0}^{\pm}<0$ and
$\sigma_{0}^{0}>0$, while outside this region, $\sigma_{0}^{0}<0$, as can be
verified numerically. Furthermore, one can see that if $\delta$ is
sufficiently large then $\sigma_{k}^{\pm,0}<0$ for all $k>0$. So, as expected,
for sufficiently large diffusion the stability of each fixed point can be
determined entirely by its stability with respect to uniform perturbation
corresponding to $k=0$.

Conversely, for sufficiently small values of $\delta$ the steady states
$u_{\pm,0}$ may lose stability with respect to non-uniform perturbations
corresponding to $k>0$. In this case the fastest growing mode corresponds to
$k=1$, and for $f^{\prime}(u_{\pm,0})(1-\alpha u_{\pm,0})-1>0$ one can define

(21) $\delta_{0,\pm}=\sqrt{\frac{f^{\prime}(u_{\pm,0})(1-\alpha
u_{\pm,0})-1}{\lambda_{1}}}.$

Then for $\delta<\delta_{0,\pm}$ it follows that $\sigma_{1}^{0,\pm}>0$.

As our subsequent analysis emphasizes the limit where $\delta\rightarrow 0$,
it is useful to examine the boundaries of the region where
$f^{\prime}(u_{\pm,0})(1-\alpha u_{\pm,0})=1$ for each fixed point. Let
$b_{1}(\alpha,\gamma)<b_{2}(\alpha,\gamma)$ be the two real roots of the
depressed quartic

(22) $b^{4}+2\gamma^{2}(1+\alpha)b^{2}-2\gamma^{2}b+\gamma^{4}=0,$

which arises from rearranging $f^{\prime}(b)(1-\alpha b)=1$. A uniform fixed
point at $u_{\pm,0}$ is stable in the limit $\delta\rightarrow 0$ provided
$b_{1}(\alpha,\beta,\gamma)<u_{\pm,0}(\alpha,\beta,\gamma)<b_{2}(\alpha,\beta,\gamma)$.
We plot a slice of the parameter space at $\beta=0.025$ in Fig. 2 in order to
illustrate how the parameter space is partitioned. One can see that
$\sigma_{1}^{0}>0$ for all points within the three root region, whereas
$\sigma_{1}^{\pm}>0$ in distinct but overlapping regions. Importantly, there
is a large parameter region in which no uniform state is linearly stable and,
therefore, pattern formation must occur from a generic initial condition. Note
that for finite $\delta$ the diffusion exerts a stabilizing effect and shrinks
the regions with positive growth rates, but even for a relatively large value
of $\delta$ (e.g. $\delta\approx 0.1$) the stability regions look nearly
indistinguishable from those in Fig. 2. As a final remark, previous
examinations of wave-pinning models have observed that even in the absence of
linear instability, such systems tend to be easily destabilized by small
finite perturbations - it is hypothesized that many real systems exist
entirely within the stable regime but nonetheless polarize due to large scale
fluctuations [40].

Figure 2: Regions of the parameter space at $\beta=0.025$ colored according to
the stability of the uniform steady state in the limit $\delta\rightarrow 0$.
Regions $i$ and $ii$ denote where there is a unique uniform steady state
$u_{0}$. In the former, $\sigma_{k}^{0}<0$ for all $k$, while in the latter
$\sigma_{k}^{0}>0$ for some $k>0$. In $iii$, $iv$ and $v$, three uniform
steady states $u_{-}<u_{0}<u_{+}$ exist, with $\sigma_{0}^{\pm}<0$ and
$\sigma_{0}^{0}>0$. In $iv$, $\sigma_{k}^{+}>0$ for some $k>0$ while
$\sigma^{-}_{k}<0$ for all $k$, and in $iii$, $\sigma_{k}^{\pm}>0$ for some
$k>0$. Finally, in $v$, $\sigma^{\pm}_{k}<0$ for all $k$. A final region
exists on the boundary of $v$ and $i$ where $u_{-}$ becomes destabilized but
$u_{+}$ remains stable, but this region is too narrow to reasonably visualize
on this plot.

##  4 Generation of Interface

At the zeroth order in $\delta$, the diffusion term in $(P_{0}^{\delta})$ may
be dropped to the leading order to obtain a closed system of equations for
$u_{0}=\lim_{\delta\to 0}u_{0}^{\delta}$:

(23) $(P_{0}^{0})\qquad\begin{aligned}
\partial_{t}u_{0}&=-u_{0}+f(u_{0})(1-\alpha U_{0}),\\\ u_{0}(x,0)&=g(x),\\\
U_{0}(t)&=\int_{\partial\Omega}u_{0}(x,t)dS,\end{aligned}$

provided $g=\lim_{\delta\to 0}g_{0}^{\delta}$ is smooth and the limit is
uniform. All spatial coupling now occurs via the parameter $U_{0}$. In
essence, this problem reduces to an infinite system of ODEs connected via a
mean-field term. Our goal in this section is to demonstrate that at almost
every point $x\in\partial\Omega$, the value of $u(x,t)$ is expected to
converge to one of two specific steady state values as $t\to\infty$ for
generic choices of $g$. This problem has been considered before in previous
studies of wave-pinning, and conventionally the local bistability of (23) for
fixed $U_{0}$ was used to justify the assumption that $u_{0}$ indeed converges
[40, 14]. Nevertheless, for variable $U_{0}$, an additional non-trivial
argument is necessary to make such a conclusion. Here, we demonstrate that
problem $(P_{0}^{0})$ is in fact a gradient flow generated by the following
energy functional:

(24)
$H[u]=\int_{\partial\Omega}\Phi(u)dS+\frac{\alpha}{2}\left[\int_{\partial\Omega}udS\right]^{2},$

where

(25)
$\Phi(u)=\frac{u(1+\beta)(u-2(1+\beta))+\gamma^{2}\log\left(u^{2}(1+\beta)+\beta\gamma^{2})\right)}{2(1+\beta)^{2}}.$

This function is bounded from below for all $\beta>0$, and the second term in
(24) is nonnegative, so we have that $H[u]$ is bounded from below. Taking the
variational derivative of $H$ yields

(26) $\partial_{t}u_{0}=-f(u_{0})\frac{\delta H}{\delta u},\qquad\frac{\delta
H}{\delta
u}=-1+\frac{u(u^{2}+\gamma^{2})}{u^{2}(1+\beta)+\beta\gamma^{2}}+\alpha\int_{\partial\Omega}udS,$

which is a gradient flow generated by the energy functional $H[u]$ with the
mobility $f(u)$. As $f(u_{0})>0$ uniformly, we have

(27)
$\frac{dH(u_{0}(\cdot,t))}{dt}=-\int_{\partial\Omega}(f(u_{0}))^{-1}|\partial_{t}u_{0}|^{2}dS\leq
0.$

Here equality only occurs when $u_{0}$ is one of the steady state solutions to
(23), which due to the dissipative nature of the dynamics are the only
elements of the attractor of the long-time dynamics of the solutions of
$(P_{0}^{0})$ [25]. Away from the steady states, $H[u_{0}(\cdot,t)]$ is both
decreasing in time and bounded below, so as $t\rightarrow\infty$,
$H(u_{0}(\cdot,t))\rightarrow H^{\infty}$. We note that existence of solutions
for problem $(P_{0}^{0})$ with positive bounded initial data follows from the
standard theory of ODEs in Banach spaces, yielding a unique solution $u_{0}\in
C^{\infty}([0,\infty);L^{2}(\partial\Omega))$ [10]. In particular, this gives
$U_{0}\in C^{\infty}([0,\infty))$. Therefore, by the standard existence theory
for ODEs we have that $u(x,\cdot)\in C^{\infty}([0,\infty))$ for each
$x\in\partial\Omega$. In fact, it is not difficult to see that the solutions
of $(P_{0}^{0})$ remain positive and bounded independently of $t$, as they
should. Thus, the arguments leading to (27) are justified. Furthermore, the
boundedness of $u_{0}$ yields boundedness of $\partial_{t}^{2}u_{0}$ as well,
which together with (27) implies that $\partial_{t}u_{0}(\cdot,t)\rightarrow
0$ in $L^{2}(\partial\Omega)$ as $t\to\infty$.

Consider a sequence of $t_{n}>0$ such that $t_{n}\to\infty$ as $n\to\infty$.
Up to an extraction of a subsequence (not relabeled), we then have
$\partial_{t}u_{0}(\cdot,t_{n})\to 0$ a.e. in $\partial\Omega$. Since
$U_{0}(t_{n})$ is bounded, upon a further extraction of a subsequence we have
$U_{0}(t_{n})\to U_{0}^{\infty}$ as $n\to\infty$. Therefore, for a.e.
$x\in\partial\Omega$ we have that $u(x,t_{n})$ converges to one of the roots
of the equation

(28)
$-\Phi^{\prime}(h)=1-\frac{h(h^{2}+\gamma^{2})}{h^{2}(1+\beta)+\beta\gamma^{2}}=\alpha
U$

for $U=U_{0}^{\infty}$. This equation has at most three real roots, which we
label $h^{-}(U)<h^{0}(U)<h^{+}(U)$, and plot their dependence on $U$ in Fig.
3. The conditions for all three roots to be real are $0<\beta<1/8$ and
$U_{-}(\alpha,\beta,\gamma)<U<U_{+}(\alpha,\beta,\gamma)$, where
$U_{\pm}(\alpha,\beta,\gamma)$ correspond to saddle node bifurcations. We
calculate the explicit dependence of these on the kinetic parameters as

(29) $\displaystyle U_{\pm}(\alpha,\beta,\gamma)$
$\displaystyle=\alpha^{-1}\left[1-\frac{\left(3\pm\sqrt{1-8\beta}\right)\sqrt{1\pm\sqrt{1-8\beta}-2\beta)}}{\sqrt{2}(1\pm\sqrt{1-8\beta})(1+\beta)^{3/2}}\gamma\right].$

Notably, the form of $\Phi^{\prime}(h)$ implies that
$\frac{dh^{\pm}(U)}{dU}<0$. Therefore, the roots
$u_{0}=h^{\pm}(U_{0}^{\infty})$ represent the two linearly stable equilibria
of (23), while $u_{0}=h^{0}(U_{0}^{\infty})$ represents an unstable
equilibrium for $U_{0}=U_{0}^{\infty}$ fixed. Thus, generically we would
expect that for almost every $x\in\partial\Omega$ we have

(30)
$\lim_{n\to\infty}u_{0}(x,t_{n})=\begin{cases}h^{+}(U_{0}^{\infty})&\text{if}\
x\in\partial\Omega_{0}^{+}(U_{0}^{\infty}),\\\
h^{-}(U_{0}^{\infty})&\text{if}\
x\in\partial\Omega_{0}^{-}(U_{0}^{\infty}),\end{cases}$

for some $\partial\Omega_{0}^{\pm}\subset\partial\Omega$ (possibly depending
on the choice of the sequence) such that
$|\partial\Omega\backslash(\partial\Omega_{0}^{+}\cup\partial\Omega_{0}^{-})|=0$.

From the definitions of $\partial\Omega_{0}^{+}$,
$U_{0}^{\infty}=\lim_{n\to\infty}\int_{\partial\Omega}u_{0}(x,t_{n})dS$, and
$H^{\infty}=\lim_{n\to\infty}H[u_{0}(\cdot,t_{n})]$, the following system of
equations arises in the limit as $n\to\infty$:

(31) $\displaystyle|\partial\Omega_{0}^{+}|+|\partial\Omega_{0}^{-}|$
$\displaystyle=1,$ (32) $\displaystyle
h^{+}(U_{0}^{\infty})|\partial\Omega_{0}^{+}|+h^{-}(U_{0}^{\infty})|\partial\Omega_{0}^{-}|$
$\displaystyle=U_{0}^{\infty},$ (33)
$\displaystyle\Phi(h^{+}(U_{0}^{\infty}))|\partial\Omega_{0}^{+}|+\Phi(h^{-}(U_{0}^{\infty}))|\partial\Omega_{0}^{-}|+\frac{\alpha}{2}|U_{0}^{\infty}|^{2}$
$\displaystyle=H^{\infty}.$

Figure 3: Three solutions $h^{+}$, $h^{0}$, and $h^{-}$ of (28) as functions
of $U$. Here $U_{+}$ denotes the location of the saddle-node bifurcation of
$h^{+}$ and $h^{0}$. Parameters chosen for this plot are $\alpha=1$,
$\beta=0.01$ and $\gamma=0.25$.

We now consider the above equations as an algebraic system with
$U_{0}^{\infty}$ fixed and $|\partial\Omega_{0}^{\pm}|$ and $H^{\infty}$ as
variables. The conservation of area (31) ensures that
$\frac{d}{dU_{0}^{\infty}}|\partial\Omega_{0}^{+}|=-\frac{d}{dU_{0}^{\infty}}|\partial\Omega_{0}^{-}|$,
and differentiating (32) yields

(34)
$\frac{dh^{+}(U_{0}^{\infty})}{dU_{0}^{\infty}}|\partial\Omega_{0}^{+}|+\frac{dh^{-}(U_{0}^{\infty})}{dU_{0}^{\infty}}|\partial\Omega_{0}^{-}|+h^{+}(U_{0}^{\infty}){d\over
dU_{0}^{\infty}}|\partial\Omega_{0}^{+}|+h^{-}(U_{0}^{\infty}){d\over
dU_{0}^{\infty}}|\partial\Omega_{0}^{-}|=1.$

Rearranging the terms then establishes that $|\partial\Omega_{0}^{+}|$ is an
increasing function of $U_{0}^{\infty}$, since

(35) ${d\over
dU_{0}^{\infty}}|\partial\Omega_{0}^{+}|=\frac{1-|\partial\Omega_{0}^{+}|\frac{dh^{+}(U_{0}^{\infty})}{dU_{0}^{\infty}}-|\partial\Omega_{0}^{-}|\frac{dh^{-}(U_{0}^{\infty})}{dU_{0}^{\infty}}}{h^{+}(U_{0}^{\infty})-h^{-}(U_{0}^{\infty})}>0.$

As the final step, taking the derivative of (33) and using the identity
$\Phi^{\prime}(h^{\pm}(U_{0}^{\infty}))=-\alpha U_{0}^{\infty}$ yields the
relationship

(36) ${dH^{\infty}\over
dU_{0}^{\infty}}=\left(\frac{\Phi(h^{+}(U_{0}^{\infty}))-\Phi(h^{-}(U_{0}^{\infty}))}{h^{+}(U_{0}^{\infty})-h^{-}(U_{0}^{\infty})}-\Phi^{\prime}(h^{+}(U_{0}^{\infty}))\right)(h^{+}(U_{0}^{\infty})-h^{-}(U_{0}^{\infty})){d\over
dU_{0}^{\infty}}|\partial\Omega_{0}^{+}|.$

The term inside the brackets only equals zero at a single value of
$U_{0}^{\infty}$, corresponding to the double tangent construction (a straight
line that touches the graph of $\Phi(u)$ at two points). For all other values
of $U_{0}^{\infty}$, we thus have $dH^{\infty}/dU_{0}^{\infty}\neq 0$, as
shown in Fig. 4 for several choices of $\alpha$. It follows that there are at
most two isolated solutions $U_{0}^{\infty}$ to (31)–(33) for a given value of
$H^{\infty}$, i.e., the set of possible values of $U_{0}^{\infty}$ for
different choices of the sequence $(t_{n})$ is discrete. Thus, since the
$\omega$-limit set of $u_{0}$ is connected [25], convergence of
$H(u_{0}(\cdot,t_{n}))\to H^{\infty}$ as $n\to\infty$ fixes $U_{0}^{\infty}$
independently of the subsequence. So we have
$\lim_{t\to\infty}U_{0}(t)=U_{0}^{\infty}$, and by the ODE stability argument
we then have $\lim_{t\to\infty}u(x,t)=h^{\pm}(U_{0}^{\infty})$ as well for
each $x\in\partial\Omega_{0}^{\pm}$, respectively. This arguments also shows
that up to sets of measure zero the sets $\partial\Omega_{0}^{\pm}$ are
independent of the sequence, and so we finally obtain a full limit:

(37)
$\lim_{t\to\infty}U_{0}(t)=U_{0}^{\infty},\qquad\lim_{t\to\infty}u_{0}(x,t)=\begin{cases}h^{+}(U_{0}^{\infty})&x\in\partial\Omega_{0}^{+},\\\
h^{-}(U_{0}^{\infty})&x\in\partial\Omega_{0}^{-},\end{cases}$

for some $\partial\Omega_{0}^{\pm}\subset\partial\Omega$ depending only on the
parameters and the initial condition and satisfying
$|\partial\Omega_{0}^{-}|+|\partial\Omega_{0}^{+}|=1$. In other words, the
solution of problem $(P_{0}^{0})$ is expected to generically converge to a
piecewise-constant function corresponding to the two stable equilibrium
branches of the bistable nonlinearity for a certain long time limit value of
$U_{0}$. This result can be interpreted as follows: for small but finite value
of $\delta$ the solution of problem $(P_{0}^{\delta})$ will become nearly
piecewise constant on the timescale $1\ll t\ll\delta^{-1}$ almost everywhere
in $\partial\Omega$, except in a small transition region of width of order
$\delta$. Thus, an interface is expected to form on this timescale, which will
serve as an approximate initial condition for problem $(P_{1}^{\delta})$ at
the corresponding rescaled timescale $\delta\ll t\ll 1$ for the latter.

Figure 4: The function $H^{\infty}(U_{0}^{\infty})$ solving (31)–(33), plotted
at different values of $\alpha$, with $\beta=0.01$ and $\gamma=0.25$.

We note that while one should expect $u=h^{\pm}(U_{0}^{\infty})$ to be
generically selected as the limiting values of the solution away from the
interfaces, for arbitrary initial conditions the possibility of
$h^{0}(U_{0}^{\infty})$ on a set of positive measure may not be a priori
excluded. This prevents us to make (37) into a rigorous conclusion about the
long-time limit behavior for problem $(P_{0}^{0})$. This sort of issue is well
known in the studies of gradient flows and its further treatment would require
a more delicate analysis that goes beyond the scope of the present paper (for
a treatment of a closely related problem, see [3]).

Lastly, we point out that the arguments above do not a priori exclude the
possibility of $|\partial\Omega_{0}^{+}|=0$ or $|\partial\Omega_{0}^{-}|=0$,
depending on the initial conditions. In these cases the system ends up in a
uniform state and no further non-trivial dynamics governed by (8) – (10) is
expected to occur. This outcome of the dynamics may be interpreted as a
failure of cell polarization. However, such an outcome would be non-generic
for the kinetic parameters corresponding to the regions $ii$ and $iii$ in Fig.
2, in which all uniform states are linearly unstable, consistently with the
conclusion at the end of Sec. 3. In contrast, the case
$0<|\partial\Omega_{0}^{+}|<1$ corresponds to a patterned state as the outcome
of the dynamics of $(P_{0}^{0})$, in which the boundary between
$\partial\Omega_{0}^{+}$ and $\partial\Omega_{0}^{-}$ represents an interface
between the high- and low-concentration domains, respectively.

##  5 Convergence of surface concentration

Having established the conditions for generation of the interface bounding
uniform domains via non-local coupling, the next step is an examination of the
evolution of these interfaces. At $O(\delta^{-1})$ timescale (in the original
variables), the dynamics of our system is described by problem
$(P_{1}^{\delta})$. We claim that in the limit $\delta\to 0$ the solution
$u_{1}^{\delta}$ of this problem converges to that of

(38) $(P_{1}^{0})\qquad\begin{aligned}
u_{1}(x,t)&=\begin{cases}h^{+}(U_{1}(t))&x\in\partial\Omega_{1}^{+}(t),\\\
h^{-}(U_{1}(t))&x\in\partial\Omega_{1}^{-}(t),\end{cases}\qquad\partial\Omega_{1}^{+}(t)\cup\Gamma(t)\cup\partial\Omega_{1}^{-}(t)=\partial\Omega,\\\
U_{1}(t)&=h^{+}(U_{1}(t))|\partial\Omega_{1}^{+}(t)|+h^{-}(U_{1}(t))|\partial\Omega_{1}^{-}(t)|,\\\
{\partial\Gamma(x,t)\over\partial t}&=c(U_{1}(t)),\qquad
x\in\Gamma(t).\end{aligned}$

Here $\partial\Omega_{1}^{\pm}(t)$ are the time-dependent subsets of
$\partial\Omega$ corresponding to the quasi-steady states $h^{\pm}(U_{1}(t))$
and $\Gamma(t)\subset\partial\Omega$ is a smooth closed curve or a collection
of curves representing the boundary between $\partial\Omega_{1}^{+}(t)$ and
$\partial\Omega_{1}^{-}(t)$, i.e., $\Gamma(t)$ is the interface. With some
abuse of notation, we denote by $\partial\Gamma(x,t)/\partial t$ the normal
velocity of $\Gamma(t)$ at each point $x\in\Gamma(t)$ in the direction of
$\partial\Omega_{1}^{-}(t)$, with the function $c(U_{1})$ to be specified. As
was already discussed at the end of Sec. 4, the sets
$\partial\Omega_{1}^{\pm}(t)$ at $t=0$ coincide with the sets
$\partial\Omega_{0}^{\pm}$ obtained in the long time limit of the solution of
problem $(P_{0}^{0})$.

To aid in the derivation of $(P_{1}^{0})$, it is convenient to define
$\Psi=\Psi(x,t)$ which is the signed distance function from $\Gamma(t)$, with
$\Psi>0$ corresponding to $u_{1}=h^{+}(U_{1}(t))$. Defining the stretched
distance $z=\Psi(x,t)/\delta$ to $\Gamma(t)$, we seek the solution of
$(P_{1}^{\delta})$ in the neighborhood of $\Gamma(t)$ in the form
$u_{1}^{\delta}(x,t)\simeq v(z,t)$, where $v$ is an unknown function [18].
Then to the leading order in $\delta$ we obtain in some tubular neighborhood
of $\Gamma(t)$:

(39) $\partial_{t}\Psi(x,t)=c(U_{1}(t)),\qquad|\nabla\Psi(x,t)|=1,$

where the function $c(U)$ is obtained by solving the traveling wave equation

(40) $\partial^{2}_{z}v-c(U)\partial_{z}v-v+f(v)(1-\alpha
U)=0,\qquad\lim_{z\to\pm\infty}v=h^{\pm}(U).$

The interface $\Gamma(t)$ is then reconstructed from the condition
$\Psi(x,t)=0$ for $x\in\Gamma(t)$. It is well known that for all values of $U$
for which the nonlinearity in (40) is of bistable type, this equation has a
unique solution (up to translations) for a unique value of $c(U)$ [19]. In
particular, there is a unique value $U=U_{1}^{\infty}$ such that
$c(U_{1}^{\infty})=0$, which is obtained by solving the equation

(41) $\int_{h^{-}(U)}^{h^{+}(U)}\left[(1-\alpha U)f(u)-u\right]du=0,$

and the sign of $c(U)$ coincides with the sign of the integral in (41) for
$U\not=U_{1}^{\infty}$. It is also possible to show that $dc(U)/dU<0$. Note
that (39) is equivalent to the last equation in (38), with the function $c(U)$
given implicitly by the solution of (40). The rest of the equations in (38)
are obtained by passing to the limit $\delta\to 0$ in the ansatz for
$u_{1}^{\delta}$.

The dependence $c(U)$ obtained from the numerical solution of (40) for a
particular choice of the parameters is illustrated in Fig. 5. This calculation
entailed time integration of

(42)
$\displaystyle\partial_{t}u(x,t)=\partial^{2}_{x}u(x,t)-u(x,t)+f(u(x,t))(1-\alpha
U)$

for fixed values of $U$, the initial condition

(43) $u(x,0)=\begin{cases}h^{+}(U)&x>0,\\\ h^{-}(U)&x\leq 0,\end{cases}$

and the boundary conditions $\partial_{x}u(\pm L,0)=0$. Integration was
performed using the Dedalus spectral package on a domain $[-L,L]$ with
$L=100$, discretized into $1024$ modes [8]. Each simulation was allowed to
relax until $t=100$ with $dt=0.1$, at which point the front velocity was
calculated as

(44) $\displaystyle c(U)={\int_{-L}^{L}\partial_{t}u(x,t)dx\over
h^{+}(U)-h^{-}(U)}.$

Figure 5: The dependence $c(U)$ obtained numerically for $\beta=0.01$ and
$\gamma=0.25$.

To study the long time behavior of the solutions of problem $(P_{1}^{0})$, we
argue as in [27] and note that the growth rate of
$|\partial\Omega_{1}^{\pm}(t)|$ is related to $c(U_{1})$ as

(45) ${d\over dt}|\partial\Omega_{+}(t)|=c(U_{1}(t))|\Gamma(t)|.$

From this, we have

(46) ${dU_{1}\over dt}=k(t)(U_{1}^{\infty}-U_{1}),$

where

(47) $\qquad
k(t)=\left(\frac{h^{+}(U_{1}(t))-h^{-}(U_{1}(t)))}{1-{dh^{+}(U_{1}(t)\over
dt}|\partial\Omega_{1}^{+}(t)|-{dh^{-}(U_{1}(t))\over
dt}|\partial\Omega_{1}^{-}(t)|}\right){c(U_{1}(t))|\Gamma(t)|\over
U_{1}^{\infty}-U_{1}(t)}.$

As $0<k(t)<\infty$ uniformly so long as the interface persists (recall that
$dc(U)/dU<0$), as $t\rightarrow\infty$ we find that $U_{1}(t)\to
U_{1}^{\infty}$ exponentially fast and, therefore,
$\Gamma(t)\to\Gamma_{1}^{\infty}$, with $\Gamma_{1}^{\infty}$ separating the
two limiting domains $\partial\Omega_{1,\infty}^{\pm}$ in $\partial\Omega$
such that $\lim_{t\to\infty}u_{1}(x,t)=h^{\pm}(U_{1}^{\infty})$ for all
$x\in\partial\Omega_{1,\infty}^{\pm}$. Again, this behavior can be interpreted
as follows: for times $1\ll t\ll\delta^{-1}$ the solution of problem
$(P_{1}^{\delta})$ will approach the long time limit of the solution of
problem $(P_{1}^{0})$, except in a transition layer of width $\delta$. This
solution corresponds to the times $\delta\ll t\ll 1$ for problem
$(P_{2}^{\delta})$ and will, therefore, serve as the approximate initial
condition for the latter.

An interesting feature of the long time behavior of problem $(P_{1}^{0})$ is
that as $t\rightarrow\infty$, the total amount $U$ of material on the surface
converges to a value that depends only on $\alpha$, $\beta$, and $\gamma$, and
is independent of the initial condition. Furthermore, the long time value of

(48)
$|\partial\Omega_{1,\infty}^{+}|={U_{1}^{\infty}-h^{-}(U_{1}^{\infty})\over
h^{+}(U_{1}^{\infty})-h^{-}(U_{1}^{\infty})}$

carries no explicit dependence of diffusion or domain geometry, other than
through the dependence of the dimensionless parameter $\alpha$ on the surface
to volume ratio $A/V$. The dependence of $U_{1}^{\infty}$ on $\alpha$ and
$\gamma$ can be expressed as $U_{1}^{\infty}=\frac{1-b(\beta)\gamma}{\alpha}$
for some function $b(\beta)$. We illustrate this fact in Fig. 7, which
demonstrates the convergence of $U(t)\rightarrow U_{\infty}$ from three
distinct initial conditions.

Some consideration is merited on the possible elimination of fronts: arbitrary
initial conditions imply an arbitrary number of closed interfaces might be
created during the initial dynamics, and some of these interfaces may vanish
under the dynamics described above. However, provided we restrict ourselves to
a parameter choice where
$h^{-}(U_{1}^{\infty}(\alpha,\beta,\gamma))<U_{1}^{\infty}(\alpha,\beta,\gamma)<h^{+}(U_{1}^{\infty}(\alpha,\beta,\gamma))$,
as well as $U_{-}<U_{1}^{\infty}<U_{+}$ as per the previous section, both
$h^{+}$ and $h^{-}$ will exist as real roots and conservation of mass will
require $|\partial\Omega_{1}^{+}(t)|,\ |\partial\Omega_{1}^{-}(t)|>0$ for all
$t$. Thus, in this regime, at least one interface exists as
$t\rightarrow\infty$.

##  6 Area-preserving geodesic curvature flow

In the previous section, we demonstrated that on the $O(\delta^{-1})$
timescale (in the original variables) the dynamics results in a stationary
interface with a prescribed quantity of matter on the surface and a fixed area
of the domains bounded by the interface. We now show that on longer timescales
there is no longer a net transport of mass between the boundary
$\partial\Omega$ and the enclosed volume $\Omega$, and the dynamics may be
reduced to a gradient flow driven by an interface length subject to domain
area conservation [34].

The dynamics on the $O(\delta^{-2})$ timescale (again, in the original
variables) is described by problem $(P_{2}^{\delta})$. We now claim that in
the limit $\delta\to 0$ the solution $u_{2}^{\delta}$ of this problem
converges to that of

(49) $(P_{2}^{0})\qquad\begin{aligned}
u_{2}(x,t)&=\begin{cases}h^{+}(U_{1}^{\infty})&x\in\partial\Omega_{2}^{+}(t),\\\
h^{-}(U_{1}^{\infty})&x\in\partial\Omega_{2}^{-}(t),\end{cases}\qquad\partial\Omega_{2}^{+}(t)\cup\Gamma(t)\cup\partial\Omega_{2}^{-}(t)=\partial\Omega,\\\
U_{1}^{\infty}&=h^{+}(U_{1}^{\infty})|\partial\Omega_{2}^{+}(t)|+h^{-}(U_{1}^{\infty})|\partial\Omega_{2}^{-}(t)|,\\\
{\partial\Gamma(x,t)\over\partial t}&=k_{g}(x,t)+\xi(t),\qquad
x\in\Gamma(t).\end{aligned}$

Here, as before, $\partial\Omega_{2}^{\pm}(t)$ are the time-dependent subsets
of $\partial\Omega$ corresponding to the quasi-steady states
$h^{\pm}(U_{1}^{\infty})$ and $\Gamma(t)\subset\partial\Omega$ is a smooth
closed curve or a collection of curves representing the boundary between
$\partial\Omega_{2}^{+}(t)$ and $\partial\Omega_{2}^{-}(t)$. The function
$k_{g}(x,t)$ refers to the geodesic curvature of $\Gamma(t)$ at
$x\in\Gamma(t)$, with the sign convention that $k_{g}<0$ if
$\partial\Omega_{2}^{+}$ is a small geodesic disk, and $\xi(t)$ is a Lagrange
multiplier ensuring conservation of $|\partial\Omega_{2}^{\pm}(t)|$. The sets
$\partial\Omega_{2}^{\pm}(t)$ at $t=0$ now coincide with the sets
$\partial\Omega_{1,\infty}^{\pm}$ obtained in the long time limit of the
solution of problem $(P_{1}^{0})$.

To derive this problem, we again denote by $z=\Psi(x,t)/\delta$ the stretched
distance to the interface and seek the solution in the form [18]

(50) $u_{2}^{\delta}(x,t)\simeq v(z-\zeta(t))+\delta w(z,x,t),$

where $v$, $w$ and $\zeta$ are to be found. Substituting this ansatz into
$(P_{2}^{\delta})$ yields

(51) $\displaystyle(\partial_{t}\Psi$
$\displaystyle-\nabla^{2}\Psi-\partial_{t}\zeta)\partial_{z}v=\delta^{-1}\left[\partial_{zz}^{2}v-v+f(v)\left(1-\alpha\int_{\partial\Omega}vdS\right)\right]$
$\displaystyle+\partial_{zz}^{2}w^{\delta}+\left[-1+f^{\prime}(v)\left(1-\alpha\int_{\partial\Omega}vdS\right)\right]w-\alpha
f(v)\int_{\partial\Omega}wdS+O(\delta).$

To eliminate the $O(\delta^{-1})$ term from the equation, we choose $v$ to be
a solution of (40) with $c(U)=0$, which exists if and only if
$\int_{\partial\Omega}vdS=U_{1}^{\infty}$. By continuity, the latter can
always be fixed by a suitable choice of $\zeta(t)=O(1)$ for a given
$\Gamma(t)$ for which the condition
$U_{1}^{\infty}=h^{+}(U_{1}^{\infty})|\partial\Omega_{2}^{+}(t)|+h^{-}(U_{1}^{\infty})|\partial\Omega_{2}^{-}(t)|$
holds.

Having now eliminated the $O(\delta^{-1})$ term from (51), we can write the
solvability condition for $w$, which is obtained by multiplying the equation
by $\partial_{z}v$ and integrating over $z$ [18]:

(52)
$\left(\partial_{t}\zeta(t)-\partial_{t}\Psi(x,t)+\nabla^{2}\Psi(y,t)\right)\int_{-\infty}^{\infty}|\partial_{z}v|^{2}dz-\alpha\int_{h^{-}(U_{1}^{\infty})}^{h^{+}(U_{1}^{\infty})}f(u)du\int_{\partial\Omega}wdS=O(\delta).$

We note that if $y\in\Gamma(t)$ is the projection of $x\in\partial\Omega$ from
a small tubular neighborhood of $\Gamma(t)$ onto $\Gamma(t)$, then
$\nabla^{2}\Psi(x,t)=\nabla^{2}\Psi(y,t)+o(1)$ for $|x-y|\ll 1$. Thus, after
dropping the $o(1)$ terms and rearranging the formula, we can rewrite (52) as

(53)
$\partial_{t}\Psi(x,t)=\nabla^{2}\Psi(y,t)+\xi(t),\qquad|\nabla\Psi(x,t)|=1,$

where

(54)
$\xi(t)=\partial_{t}\zeta(t)-\frac{\alpha\int_{h^{-}(U_{1}^{\infty})}^{h^{+}(U_{1}^{\infty})}f(u)du}{\int_{-\infty}^{\infty}|\partial_{z}v|^{2}dz}\int_{\partial\Omega}w(x,t)dS.$

The obtained equation (53) is nothing but the last equation in (49) written in
terms of the distance function [2], recalling the definition (with our sign
convention) of geodesic curvature $k_{g}(x,t)=\nabla\cdot n(x,t)$, where
$n(x,t)=\nabla\Psi(x,t)$ is the unit normal to the curve $\Gamma(t)$ at
$x\in\Gamma(t)$ that lies in the tangent plane to $\partial\Omega$ at $x$ and
points towards $\partial\Omega_{2}^{+}(t)$. Finally, conservation of
$|\partial\Omega_{2}^{\pm}(t)|$ follows by passing to the limit $\delta\to 0$
in our ansatz, and the initial condition is given by the long time limit of
the solution of problem $(P_{1}^{0})$, i.e., we have
$\Gamma(t)=\Gamma_{1}^{\infty}$ for $t=0$. Note that the value of
$|\partial\Omega_{2}^{+}|$ is given explicitly by the right-hand side of (48).

Since the area of $|\partial\Omega_{2}^{\pm}(t)|$ is constant for the
solutions of problem $(P_{2}^{0})$, we have
$\partial_{t}|\partial\Omega_{2}^{\pm}(t)|=0$. Furthermore, since the growth
rate of a region’s area is equal to the integral of its normal velocity, we
can eliminate $\xi(t)$ to obtain a purely geometric expression for the normal
velocity:

(55) ${\partial\Gamma(x,t)\over\partial t}=k_{g}(x,t)-\langle
k_{g}(\cdot,t)\rangle_{\Gamma(t)},\qquad x\in\Gamma(t),$

where $\langle
k_{g}(\cdot,t)\rangle_{\Gamma(t)}=\frac{1}{|\Gamma(t)|}\int_{\Gamma(t)}k_{g}(x,t)ds$
and $|\Gamma(t)|$ is the length of $\Gamma(t)$. This equation describes an
area-preserving geodesic curvature flow on a surface. A geometric identity
relates the evolution of interface length to curvature as
$\frac{d}{dt}|\Gamma(t)|=-\int_{\Gamma(t)}k_{g}(x,t){\partial\Gamma(x,t)\over\partial
t}ds$ [30]. From this, we obtain an inequality

(56) $\frac{1}{|\Gamma(t)|}\frac{d}{dt}|\Gamma(t)|=\langle
k_{g}(\cdot,t)\rangle^{2}_{\Gamma(t)}-\langle
k_{g}^{2}(\cdot,t)\rangle_{\Gamma(t)}\leq 0,$

which demonstrates that the interface length $|\Gamma(t)|$ is strictly
decreasing except for when $k_{g}$ is constant along $\Gamma(t)$. The
interpretation of this result is that for small $\delta$ the long time limits
of the solutions of problem $(P_{2}^{\delta})$ will be close to domains of
prescribed area and locally minimal interfaces of constant curvature (geodesic
disks).

The problem of identifying locally minimal interfaces on a given surface has
been studied at length over the years, and while no comprehensive solution
exists, there is a considerable literature of results which can inform our
understanding [50, 38, 28]. In particular, it is known that for every $0<a<1$
there exists a set $\partial\Omega^{+}\subset\partial\Omega$ whose boundary
$\Gamma$ globally minimizes the perimeter among all sets with
$|\partial\Omega^{+}|=a$ [50]. In this case $\Gamma$ is a smooth curve of
constant geodesic curvature. Framed in terms of our reaction-diffusion
problem, this means that for appropriate choices of $\alpha$, $\beta$, and
$\gamma$, a non-uniform steady state for $(P_{2}^{0})$ should always exist
regardless of the particular shape of $\Omega$. Whether the solution of
problem $(P_{2}^{0})$ converges as $t\to\infty$ to a global minimizer, or,
indeed, converges, however, is not known a priori. We note that due to the
area conservation the interface $\Gamma(t)$ associated with the solutions of
problem $(P_{2}^{0})$ cannot vanish. Furthermore, due to the gradient flow
nature of the dynamics the $\omega$-limit set of these solutions consists only
of steady states, and convergence to one of the steady states is guaranteed if
the set of the steady states is discrete. In particular, this would be the
case if the global minimizer of the perimeter were the unique constant
curvature solution with the prescribed area. We note that the latter would
represent the most robust and reproducible scenario for cell polarization.

We can also use existing variational results to construct surfaces which yield
some surprising solutions. It has been widely suggested that strictly mass-
conserving reaction-diffusion equations will exhibit uninterrupted coarsening
and admit steady states with only a single domain, but this phenomenon has
only been studied in the Euclidean setting [7, 55]. Multi-domain polarization
has in fact been observed in real cells, and failure to account for this is
seen as a failure point of minimal models such as the one in this paper [13].
However, we find that when geometric effects are properly considered, multi-
cap solutions are possible on some surfaces, which we now demonstrate by
constructing examples. Per the theorems in Ref. [49], if $\partial\Omega$ is a
surface of revolution about the $z$-axis that is symmetric with respect to
reflections around the $xy$-plane, and if $\partial\Omega$ has Gaussian
curvature $K(z)$ that is strictly increasing in $|z|$, then the locally
minimal interfaces bounding a fixed area can be fully classified [49], and are
of the types depicted in Fig. 6. All stable interfaces lie in the planes
$z=\text{const}$, which means $k_{g}$ can be written as a function of $z$, and
it is convenient to define a function on $\partial\Omega$ of the form
$L(z)=k_{g}^{2}(z)+K(z)$. In addition to single-interface solutions, there are
two possible solutions with disjoint interfaces. Case 1 features two
interfaces at $\pm z_{0}$ for some $z_{0}>0$, respectively, which are stable
provided they lie within the region where $L(z_{0})<0$. Case 2 has one
interface at $z_{1}$ (s.t. $L(z_{1})<0)$ and the other at some $z_{2}$ such
that $z_{1}z_{2}<0$ and $L(z_{2})>0$. As the positions of $\partial\Omega^{+}$
and $\partial\Omega^{-}$ are interchangeable under geodesic curvature flow,
each case can either represent a two cap solution, or a single high-
concentration band containing $z=0$. In essence, these example solutions show
the existence of a geometry-induced interruption to coarsening. As an aside,
results such as those in Refs. [49, 50, 38] do rule out multi-cap solutions on
symmetric convex surfaces such as spheroids.

##  7 Numerical Study

Figure 6: Numerical steady state solutions of (8) – (10) on a surface of
revolution (see the main text for details). The region of $L\leq 0$, bounded
by black dashed lines here, can support solutions with multiple stable
interfaces. Parameter choices are $\beta=0.05$, $\gamma=0.2$ and
$\delta=0.01$, and from left to right the values of $\alpha$ used were $10$,
$3$, $13$, $10$ and $3.2$.

We now verify numerically the existence of nontrivial steady state solutions
predicted by the asymptotic analysis of $(P_{2}^{\delta})$. We simulated (8) –
(10) on the surface shown in Fig. 6 for each of the possible configurations,
using a recently developed numerical scheme [36]. The values of $\alpha$ used
were (from left to right): $1$, $0.33$, $1.3$, $1$ and $0.35$. For all cases,
$\beta=0.05$, $\gamma=0.2$ and $\delta=0.01$. The surface geometry was
constructed by revolution of polar curve $r(\theta)=2a^{2}\left(\cos
2\theta+\sqrt{(b/a)^{4}-\sin^{2}2\theta}\right)$ (for $a=0.2728$ and
$b=0.3896$) about the vertical axis. The simulations were performed on a
spatial discretization of $128\times 128$ modes and integrated with $dt=0.1$
until $t=500$, at which point a numerical steady state was achieved. In each
case, the initial conditions were step functions approximating the desired
steady state. The steady states were verified using a GMRES-based Newton-
Krylov method [29]. This method further allowed us to compute approximate
values of the dominant eigenvalues, and in each case these were negative,
consistent with the proposed stability of these states. While the constraint
on $U_{\infty}$ for the single interface case is as described in the previous
section, both case 1 and 2 of the double interface solutions have a narrower
range of acceptable $U_{\infty}$, owing to the constraints on the locations of
the interfaces. As such, we chose a different $\alpha$ for each simulation to
ensure a viable solution.

To corroborate the multiple stages of the dynamics predicted by our asymptotic
analysis, we performed two additional numerical studies of (8) – (10), again
using the integration methods introduced in Ref. [36]. Recall that the results
of Sec. 4 indicate convergence of $U(t)$ to a constant value $U_{1}$ depending
on the initial condition on a large $O(1)$ time scale. Furthermore, by the
results of Sec. 5 the function $U(t)$ converges to the final value
$U_{1}^{\infty}(\alpha,\beta,\gamma)$ independent of initial conditions on a
large $O(\delta^{-1})$ time scale. As the first test of this, we examine the
convergence of $U(t)$ in three simulations. These were initialized at
$u_{0}+\gamma\sigma(x)/2$, where $u_{0}$ is as defined in Sec. 3, and
$\sigma(x)$ is a random function taking on values between $-1$ and $1$, with
the spatial correlation length $0.6$. Further details on how this distribution
was generated can be found in the supplementary materials of Ref. [36].
Simulations are performed on a sphere of unit area, discretized to $256\times
256$ modes, and integrated with a timestep of $dt=0.1$. The results of the
simulations for $\delta=10^{-3}$ are plotted in Fig. 7 (see the figure caption
for the values of other parameters). The trajectories of all three examples
can be split into two types of behavior. Dynamics from $t_{0}=0$ to
$t_{1}\approx 10^{1}$ exhibit a rapid change in $U(t)$ towards an
approximately constant value: this is the timescale of the interface
generation. On the other hand, on the interval from $t_{1}\approx 10^{1}$
through $t_{2}\approx 10^{4}$, the different instances of $U(t)$ converge to
the predicted steady state at $U_{1}^{\infty}$ from Sec. 5.

Figure 7: Numerical trajectories of $U(t)$ from three distinct initial
conditions. The three simulations have different initial conditions, but the
same model parameters, and converge to the value
$U_{1}^{\infty}(\alpha,\beta,\gamma)$ from the asymptotic analysis of Sec. 5
as $t\rightarrow\infty$. Simulations performed at $\alpha=2.25$,
$\beta=0.0001$, $\gamma=0.25$ and $\delta=10^{-3}$ on a sphere.

A second numerical experiment is to study the behavior of the interface length
over time. For $\delta\ll 1$, an effective proxy for the total interface
length is the quantity

(57) $\langle|\nabla u|\rangle(t)=\int_{\partial\Omega}|\nabla u(\cdot,t)|dS,$

as in the last stage of the evolution the profile $u(x,t)$ will be close to
two fixed values everywhere outside the thin transition layers approximating
the interfaces. We initialize a simulation at a small perturbation about the
steady state as $u_{0}+0.01\sigma(x)$ and plot the evolution of the interface
in Fig. 8. The simulation domain was a unit-area prolate spheroid with an
aspect ratio of $2$, discretized into $128\times 128$ modes. As before, our
timestep was $dt=0.1$.

For our choice of $\delta=10^{-2}$ here, the distinct regimes corresponding to
the limit behaviors of $(P_{0}^{\delta})$, $(P_{1}^{\delta})$ and
$(P_{2}^{\delta})$ as $\delta\to 0$ are very clearly expressed. Much as in the
previous example, the timescale of interface formation is on the order of
$t\approx 10^{1}$: the order of magnitude increase in $\delta$ between the two
figures has a minimal effect on the zeroth order behavior of the system. At
the end of this interval, surface concentration is fully partitioned into
approximately uniform domains $\partial\Omega^{\pm}$. Once the quantity
$\langle|\nabla u|\rangle$ peaks, it begins a period of slower evolution
consistent with the predicted front motion of $(P_{1}^{0})$, which then
plateaus once $U_{1}^{\infty}$ is reached. A slower evolution of
$\langle|\nabla u|\rangle$, which is now monotonically decreasing,
subsequently occurs under the approximate $(P_{2}^{0})$ dynamics. A brief note
should be made of the step-like shape of the plot in this final interval.
Contraction under $(P_{1}^{0})$ front motion reduces this simulation to two
spots, and the rapid drop in interface length which occurs around $t=2000$
represents the annihilation of one of these. After this, the remaining spot
persists for a long period near midway between the two poles before finally
migrating to one end of the spheroid. The final steady state is at the point
of maximum Gaussian curvature, consistent with the geometric results noted in
the previous section. Taken together, these numerical tests represent strong
evidence that the formal asymptotics developed here accurately reproduces the
behavior of (8) – (10) for small but finite choices of $\delta$.

Figure 8: Tracking the quantity $\langle|\nabla
u|\rangle=\int_{\partial\Omega}|\nabla u|dS$ shows us the growth and shrinkage
of the interface over time. Evolution was simulated at
$\alpha=0.58,\beta=0.001,\gamma=0.4$, and $\delta=0.01$. We observe differing
dynamical regimes each roughly two decades in timescales: initial growth as
the interface forms, followed by evolution under an eikonal-type flow, and
further monotonic interface length decrease under geodesic curvature flow. The
three inset plots show the simulation near the approximate end time of each
regime.

##  8 Discussion

In summary, we performed an extensive analysis of the asymptotic behavior of a
minimal model of cell polarization on a 2D manifold. Our results establish a
hierarchy of three timescales of the dynamics in the case of small surface
diffusivity. On the short timescale, we derived an evolution equation in the
form of a globally coupled system of ODEs. While our initial problem itself
appears to lack a variational structure, the reduced dynamics is shown to be a
gradient flow generated by a nonlocal energy functional. From this, we were
able to demonstrate generation of interfaces between high- and low-
concentration regions for generic initial data. Our argument relied on the
assumption that the solution for almost all points of the surface approached
one of the two stable branches of the bistable nonlinearity, which is expected
to be generically true but cannot be established for all initial data to make
the result fully rigorous. Let us mention that in Ref. [3], in a closely
related nonlocal Allen-Cahn problem, the authors were able to avoid this kind
of assumption by constructing a second Lyapunov functional (in fact, an
infinite family of Lyapunov functionals) to prove full convergence in time at
the expense of allowing unstable solutions in the limit as well. However, the
direct coupling between the nonlocal and nonlinear terms in problem
$(P_{0}^{0})$ greatly complicates this sort of argument, so it remains unclear
whether this approach could be adapted to our problem.

On the intermediate timescale, we demonstrated in the sharp-interface limit
that the interfacial velocity is uniform and depends only on the enclosed
area, and that this enclosed area approaches a final value determined only by
the kinetic rates. This area is determined strictly by the dimensionless
parameters of our model, and thus is independent of the surface geometry. This
fact raises a question as to how robust cell polarization can be to
perturbations of the cell shape. This question could be particularly
interesting in light of recent interest in coupled chemo-mechanical models of
deforming cell membranes [37, 35, 17, 22].

Finally, on the long timescale, we demonstrated that the dynamics reduce to
area-preserving interface motion by geodesic curvature. This result formally
establishes the intuitive notion that on curved surfaces wave-pinning
phenomena is driven by a gradient flow generated by the interface length. The
steady state of such a flow has a highly nontrivial dependence on surface
shape - the ability of some surface geometries to support several distinct
solutions suggests important implications for cellular decision-making and
addresses a perceived failure of minimal mass-conserving reaction-diffusion
models. One particularly interesting question for future work characterizing
mechanisms to select between different steady states. Very recently
theoretical work has shown that chemical gradients can enforce state
transitions in simple geometries [9]. Examining similar effects on geometries
such as those in Fig. 6 could provide useful insight into biological control
mechanisms.

There is merit in some further comparison between our results and some
deceptively similar results from flat surfaces. In the 1D version of this
system, multi-domain solutions can also form, but these are strictly
metastable and vanish on exponentially slow timescales [40, 53]. By contrast,
the multi-domain solutions discussed here can be truly stable. Previous works
such as Ref. [57] have also showed that domain geometry can establish multi-
domain solutions in wave-pinning models. One key difference is that these
studies examined flat 2D domains bounded by no-flux boundaries of various
shapes, rather than curved surfaces. Both versions can produce multi-domain
solutions, but a key difference lies in the time needed to reach steady
states. Specifically, a circular domain in the flat scenario will be
metastable if it is sufficiently far from the boundary, and like in the 1D
case discussed above it will only approach the steady state location
exponentially slowly. This slowness stems from the approximate translational
invariance of the interfaces away from the boundaries. Non-uniform local
curvature breaks this invariance, making the convergence to steady state
domain locations considerably faster, on the algebraic rather than
exponentially long timescale in terms of the inverse diffusivity.

The wave-pinning model studied here represents one of the simplest possible
characterizations of coupled bulk-surface dynamics in a single cell. A great
number of more complex models have been proposed over the past decade,
including ones featuring multiple chemical species [43] and advection [23, 4].
Additionally, recent experiments have demonstrated that bulk-surface systems
are capable of much more complex pattern formation than discussed here, such
as the spiral wave dynamics observed in starfish embryo [32, 54], establishing
a long-term need for further analysis of this class of system.

## References

  * [1] M. Alfaro and H. Matano, On the validity of formal asymptotic expansions in Allen-Cahn equation and FitzHugh-Nagumo system with generic initial data, Discrete Contin. Dyn. Syst. Ser. B, 17 (2012), pp. 1639–1649.
  * [2] L. Ambrosio, Geometric evolution problems, distance function and viscosity solutions, in Calculus of Variations and Partial Differential Equations: Topics on Geometrical Evolution Problems and Degree Theory, G. Buttazzo, A. Marino, and M. Murthy, eds., Springer, Berlin, Heidelberg, 2000, pp. 5–93.
  * [3] J. M. Ball and Y. Şengül, Quasistatic nonlinear viscoelasticity and gradient flows, J. Dyn. Differ. Equ., 27 (2015), pp. 405–442.
  * [4] L. Barberi and K. Kruse, Localized states in active fluids, arXiv:2209.02581, (2022).
  * [5] G. Bellettini and M. Paolini, Quasi-optimal error estimates for the mean curvature flow with a forcing term, Differ. Integral Equ., 8 (1995), pp. 735 – 752.
  * [6] S. Bialecki, B. Kazmierczak, and T. Lipniacki, Polarization of concave domains by traveling wave pinning, PloS one, 12 (2017), p. e0190372.
  * [7] F. Brauns, H. Weyer, J. Halatek, J. Yoon, and E. Frey, Wavelength selection by interrupted coarsening in reaction-diffusion systems, Phys. Rev. Lett., 126 (2021), p. 104101.
  * [8] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, Dedalus: A flexible framework for numerical simulations with spectral methods, Phys. Rev. Research, 2 (2020), p. 023068.
  * [9] A. Buttenschön and L. Edelstein-Keshet, Cell repolarization: A bifurcation study of spatio-temporal perturbations of polar cells, Bull. Math. Biol., 84 (2022), pp. 1–29.
  * [10] H. Cartan, Calcul Differentiel: I-Calcul differetiel dans les espaces de Banach; II-Equations differentielles, Cours de mathematiques II, Hermann et Cie, Editeurs, 1967.
  * [11] X. Chen, Generation and propagation of interfaces for reaction-diffusion equations, J. Differ. Equ., 96 (1992), pp. 116–141.
  * [12] X. Chen, Generation and propagation of interfaces in reaction-diffusion systems, Trans. Amer. Math. Soc., 334 (1992), pp. 877–913.
  * [13] J.-G. Chiou, K. D. Moran, and D. J. Lew, How cells determine the number of polarity sites, eLife, 10 (2021), p. e58768.
  * [14] D. Cusseddu, L. Edelstein-Keshet, J. A. Mackenzie, S. Portet, and A. Madzvamuse, A coupled bulk-surface model for cell polarisation, J. Theor. Biol., 481 (2019), pp. 119–135.
  * [15] R. Diegmiller, H. Montanelli, C. B. Muratov, and S. Y. Shvartsman, Spherical caps in cell polarization, Biophys. J., 115 (2018), pp. 26–30.
  * [16] L. Edelstein-Keshet, W. R. Holmes, M. Zajac, and M. Dutot, From simple to detailed models for cell polarization, Philos. Trans. R. Soc. Lond., B, 368 (2013), p. 20130003.
  * [17] C. M. Elliott and L. Hatcher, Domain formation via phase separation for spherical biomembranes with small deformations, Eur. J. Appl. Math., 32 (2021), pp. 1127–1152.
  * [18] P. C. Fife, Dynamics of Internal Layers and Diffusive Interfaces, Society for Industrial and Applied Mathematics, Philadelphia, 1988.
  * [19] P. C. Fife and J. B. McLeod, The approach of solutions of nonlinear diffusion equations to traveling front solutions, Arch. Rational Mech. Anal., 65 (1977), pp. 335–361.
  * [20] R. Geßele, J. Halatek, L. Würthner, and E. Frey, Geometric cues stabilise long-axis polarisation of PAR protein patterns in _C. elegans_ , Nat. Commun., 11 (2020), pp. 1–12.
  * [21] D. Ghose, T. Elston, and D. Lew, Orientation of cell polarity by chemical gradients, Annu. Rev. Biophys., 51 (2022), pp. 431–451.
  * [22] D. Gomez, S. Iyaniwura, F. Paquin-Lefebvre, and M. Ward, Pattern forming systems coupling linear bulk diffusion to dynamically active membranes or cells, Philos. Trans. R. Soc. A, 379 (2021), p. 20200276.
  * [23] P. Gross, K. V. Kumar, N. W. Goehring, J. S. Bois, C. Hoege, F. Jülicher, and S. W. Grill, Guiding self-organized pattern formation in cell polarity establishment, Nat. Phys., 15 (2019), pp. 293–300.
  * [24] J. Halatek and E. Frey, Rethinking pattern formation in reaction–diffusion systems, Nat. Phys., 14 (2018), pp. 507–514.
  * [25] J. K. Hale, Asymptotic behavior of dissipative systems, vol. 25 of Mathematical Surveys and Monographs, American Mathematical Society, Providence, RI, 1988.
  * [26] S. Hausberg and M. Röger, Well-posedness and fast-diffusion limit for a bulk–surface reaction–diffusion system, Nonlinear Differ. Equ. Appl., 25 (2018), pp. 1–32.
  * [27] M. Henry, D. Hilhorst, and C. B. Muratov, A multiple scale pattern formation cascade in reaction-diffusion systems of activator-inhibitor type, Interfaces Free Boundaries, 20 (2018), pp. 297–336.
  * [28] H. Howards, M. Hutchings, and F. Morgan, The isoperimetric problem on surfaces, Am. Math. Mon., 106 (1999), pp. 430–439.
  * [29] C. T. Kelley, I. Kevrekidis, and L. Qiao, Newton-krylov solvers for time-steppers, arXiv:math/0404374, (2004).
  * [30] M. KolÁŘ, M. Beneš, and D. ŠevČoviČ, Area preserving geodesic curvature driven flow of closed curves on a surface, Discrete Contin. Dyn. Syst. - B., 22 (2017), pp. 3671–3689.
  * [31] J. Li, L. Su, X. Wang, and Y. Wang, Bulk-surface coupling: Derivation of two models, J. Differ. Equ., 289 (2021), pp. 1–34.
  * [32] J. Liu, J. F. Totz, P. W. Miller, A. D. Hastewell, Y.-C. Chao, J. Dunkel, and N. Fakhri, Topological braiding and virtual particles on the cell membrane, Proc. Natl. Acad. Sci. U.S.A., 118 (2021), p. e2104191118.
  * [33] A. Madzvamuse, A. H. Chung, and C. Venkataraman, Stability analysis and simulations of coupled bulk-surface reaction–diffusion systems, Proc. R. Soc. A, 471 (2015), p. 20140546.
  * [34] C. Mantegazza, Lecture Notes on Mean Curvature Flow, vol. 290 of Progress in Mathematics, Birkhauser, Basel, 2011.
  * [35] P. W. Miller and J. Dunkel, Gait-optimized locomotion of wave-driven soft sheets, Soft Matter, 16 (2020), pp. 3991–3999.
  * [36] P. W. Miller, D. Fortunato, C. Muratov, L. Greengard, and S. Shvartsman, Forced and spontaneous symmetry breaking in cell polarization, Nat. Comput. Sci., 2 (2022), pp. 504–511.
  * [37] P. W. Miller, N. Stoop, and J. Dunkel, Geometry of wave propagation on active deformable surfaces, Phys. Rev. Lett., 120 (2018), p. 268001.
  * [38] F. Morgan and D. L. Johnson, Some sharp isoperimetric theorems for Riemannian manifolds, Indiana Univ. Math. J., (2000), pp. 1017–1041.
  * [39] Y. Mori, A. Jilkine, and L. Edelstein-Keshet, Wave-pinning and cell polarity from a bistable reaction-diffusion system, Biophys. J., 94 (2008), pp. 3684–3697.
  * [40] Y. Mori, A. Jilkine, and L. Edelstein-Keshet, Asymptotic and bifurcation analysis of wave-pinning in a reaction-diffusion model for cell polarization, SIAM J. Appl. Math., 71 (2011), pp. 1401–1427.
  * [41] B. Niethammer, M. Röger, and J. J. Velázquez, A bulk-surface reaction-diffusion system for cell polarization, Interfaces Free Boundaries, 22 (2020), pp. 85–117.
  * [42] M. Otsuji, S. Ishihara, C. Co, K. Kaibuchi, A. Mochizuki, and S. Kuroda, A mass conserved reaction–diffusion system captures properties of cell polarity, PloS Comput. Biol., 3 (2007), p. e108.
  * [43] F. Paquin-Lefebvre, W. Nagata, and M. J. Ward, Pattern formation and oscillatory dynamics in a two-dimensional coupled bulk-surface reaction-diffusion system, SIAM J. Appl. Dyn., 18 (2019), pp. 1334–1390.
  * [44] A. Pisante and F. Punzo, Allen–Cahn approximation of mean curvature flow in riemannian manifolds, II: Brakke’s flows, Commun. Contemp. Math., 17 (2015), p. 1450041.
  * [45] A. Pisante and F. Punzo, Allen-Cahn approximation of mean curvature flow in riemannian manifolds I, uniform estimates, Ann. Sc. Norm. Super. Pisa Cl. Sci. (5), 15 (2016), pp. 309–341.
  * [46] W.-J. Rappel and L. Edelstein-Keshet, Mechanisms of cell polarization, Curr. Opin. Syst. Biol., 3 (2017), pp. 43–53.
  * [47] A. Rätz, Turing-type instabilities in bulk–surface reaction–diffusion systems, J. Comput. Appl. Math., 289 (2015), pp. 142–152.
  * [48] A. Rätz and M. Röger, Symmetry breaking in a bulk–surface reaction–diffusion model for signalling networks, Nonlinearity, 27 (2014), pp. 1805–1827.
  * [49] M. Ritoré, Constant geodesic curvature curves and isoperimetric domains in rotationally symmetric surfaces, Commun. Anal. Geom., 9 (2001), pp. 1093–1138.
  * [50] A. Ros, The isoperimetric problem, in Global theory of minimal surfaces, D. A. Hoffman, ed., vol. 2 of Clay mathematics proceedings, American Mathematical Society, Providence, Rhode Island, 2001, pp. 175–209.
  * [51] K. Sakamoto, Spatial homogenization and internal layers in a reaction-diffusion system, Hiroshima Math. J., 30 (2000), pp. 377–402.
  * [52] V. Sharma and J. Morgan, Global existence of solutions to reaction-diffusion systems with mass transport type boundary conditions, SIAM J. Math. Anal., 48 (2016), pp. 4202–4240.
  * [53] X. Sun and M. J. Ward, Dynamics and coarsening of interfaces for the viscous cahn—hilliard equation in one spatial dimension, Stud. Appl. Math., 105 (2000), pp. 203–234.
  * [54] T. H. Tan, J. Liu, P. W. Miller, M. Tekant, J. Dunkel, and N. Fakhri, Topological turbulence in the membrane of a living cell, Nat. Phys., 16 (2020), pp. 657–662.
  * [55] M. Tateno and S. Ishihara, Interfacial-curvature-driven coarsening in mass-conserved reaction-diffusion systems, Phys. Rev. Research, 3 (2021), p. 023198.
  * [56] P. K. Trong, E. M. Nicola, N. W. Goehring, K. V. Kumar, and S. W. Grill, Parameter-space topology of models for cell polarity, New J. Phys., 16 (2014), p. 065009.
  * [57] B. Vanderlei, J. J. Feng, and L. Edelstein-Keshet, A computational model of cell polarization and motility coupling mechanics and biochemistry, Multiscale Model. Simul., 9 (2011), pp. 1420–1443.
  * [58] M. Zhu, J. Cornwall-Scoones, P. Wang, C. E. Handford, J. Na, M. Thomson, and M. Zernicka-Goetz, Developmental clock and mechanism of de novo polarization of the mouse embryo, Science, 370 (2020), p. eabd2703.


acronymnohyperfirsttrue [acronym]long-short [a]Laurin Pannullo

# Flattening of the quantum effective potential in fermionic theories

Gergely Endrődi  Tamás G. Kovács  Gergely Markó

###### Abstract

We present methods to constrain fermionic condensates on the level of the path
integral, which grant access to the quantum effective potential in the
infinite volume limit. In the case of a spontaneously broken symmetry, this
potential possesses a manifestly flat region, which is inaccessible to the
standard approach on the lattice. However, by constraining the appropriate
order parameters such as the chiral condensate, one is then able to probe the
flat region. We demonstrate our method of constraining fermionic condensates
in the 2-dimensional Gross-Neveu model, which exhibits a spontaneously broken
chiral symmetry. We show how the potential flattens for increasing volume and
that the flat region is dominated by inhomogeneous field configurations.

##  1 Introduction

In this work, we consider fermionic theories that spontaneously break chiral
symmetry, resulting in a non-zero chiral condensate. The standard procedure to
investigate the breaking of such a symmetry on the lattice would be to
introduce an explicit breaking of the symmetry and then taking the double
limit of infinite volume and subsequently zero explicit breaking. Such a
procedure can be extremely expensive in theories like Quantum chromodynamics
(QCD). An alternative approach would be to investigate the quantum effective
potential $\Gamma$, which is the Legendre transform of the free energy with
respect to the expectation value of the order parameter – in our case the
chiral condensate. It then incorporates all quantum effects of the theory.
This manifestly convex potential features a flat region in the case of a
broken symmetry, with the edges of the flat region corresponding to the ground
state manifold. Thus, if one has access to this object, it is rather trivial
to investigate whether the symmetry is broken and what the value of the order
parameter is.

While one cannot access this quantity directly in lattice simulations, it is
possible to employ the constraint path integral formalism. Here, one
introduces a constraint of the form $\delta(\phi-\varphi)$ on a given order
parameter $\varphi$ to the path integral. The resulting constraint path
integral $\mathcal{Z}_{\phi}$ only integrates over the configurations that
fulfill the constraint, i.e., $\varphi=\phi$. The so-called constraint
potential $\Omega(\phi)=-\ln\mathcal{Z}_{\phi}/V$ can be shown to coincide
with the quantum effective potential $\Gamma(\phi)$ in the infinite volume
limit [1]. This gives an alternative approach to investigate a spontaneous
symmetry breaking, where one only has to take the infinite volume limit
instead of the double limit that we discussed before. While a constraint of
bosonic order parameters has been successfully applied on the lattice [2, 3],
it has not been used so far for fermionic order parameters such as the chiral
condensate. In this work, we present methods to constrain the chiral
condensate and obtain the constraint potential. We test these methods in the
$2$-dimensional Gross-Neveu (GN) model [4] in the $N_{f}\to\infty$ limit.
While this presents a rather simple setup, our method is applicable to other
fermionic theories such as QCD. This work presents preliminary investigations
of our developments that we will discuss in full detail in an upcoming
publication [5].

##  2 The Gross-Neveu model

The $2$-dimensional GN model [4] is a so-called four-Fermi theory with its
bosonized Euclidean action obtained from a Hubbard-Stratonovich (HS)
transformation given by

$\displaystyle
S_{\sigma}=\int\mathrm{d}^{2}x\left[-\bar{\psi}(x)\left(\not{\partial}+\sigma(x)\right)\psi(x)+\tfrac{N_{f}}{2g^{2}}\sigma^{2}(x)\right].$
(1)

where $\bar{\psi},\psi$ are fermion spinors with $N_{\gamma}=2$ components and
$N_{f}$ species, and the coupling $g^{2}$ controls the coupling strength of
the scalar four-Fermi interaction $\left(\bar{\psi}\psi\right)^{2}$ in the
original action before bosonization. The bosonic field $\sigma$ is an
auxiliary field that is introduced by the HS transformation and its
expectation value is proportional to the chiral condensate as given by the
Ward identity

$\displaystyle\langle\bar{\psi}\psi(x)\rangle=\tfrac{N_{f}}{g^{2}}\langle\sigma(x)\rangle.$
(2)

One can integrate over the fermionic degrees of freedom in the path integral
to obtain

$\displaystyle\mathcal{Z}=\int\mathcal{D}[\sigma]\,\mathrm{e}^{-S_{\mathrm{eff}}},\quad
S_{\mathrm{eff}}=\int\mathrm{d}^{2}x\,\tfrac{N_{f}}{2g^{2}}\sigma^{2}(x)-\ln\mathrm{Det}\,Q[\sigma],\quad
Q[\sigma]=\not{\partial}+\sigma.$ (3)

An often considered limit is the limit $N_{f}\to\infty$, which suppresses all
quantum fluctuations for the bosonic field $\sigma$. In this limit, the only
configurations $\Sigma$ that contribute to the path integral are the ones with
the largest weight. These are the ones that minimize the effective action
$S_{\mathrm{eff}}$. Thus, the calculation of the path integral reduces to a
minimization problem. We exclusively consider this limit in this work.

The theory exhibits a discrete chiral symmetry under which the fermionic and
bosonic fields transform as

$\displaystyle\psi\to\gamma_{5}\psi,\quad\bar{\psi}\to-\bar{\psi}\gamma_{5},\quad\sigma\to-\sigma.$
(4)

This symmetry is spontaneously broken in the $N_{f}\to\infty$ limit and one
finds a non-zero expectation value $\langle\bar{\psi}\psi\rangle$. We
renormalize the theory, by requiring that the expectation value of $\sigma$
assumes a homogeneous non-zero value, i.e., $\langle\sigma\rangle=\sigma_{0}$,
which is achieved by choosing an appropriate value of the coupling (see Refs.
[6] for further details on the renormalization).

In this work, we consider a naive lattice discretization of the GN model. We
refer to Ref. [7] for a detailed discussion of the GN model within this
discretization.

##  3 Constraining the Gross-Neveu model

We are interested in constraining the expectation value of the chiral
condensate and obtaining the effective potential as a function of this
expectation value. We consider two options to achieve this:

  1. 1.

a fermionic constraint, where we constrain the chiral condensate directly,

  2. 2.

a bosonic constraint, where we constrain the auxiliary field, which is linked
to the chiral condensate via the Ward identity 2.

###  3.1 Fermionic constraint

With the fermionic constraint, we modify the bosonized path integral by
introducing a constraint on the spatial average of the condensate to obtain

$\displaystyle\mathcal{Z}_{\phi^{(\mathrm{f})}}=\int\mathcal{D}[\bar{\psi},\psi,\sigma]\
\mathrm{e}^{-S_{\sigma}}\,\delta\\!\left(\tfrac{1}{V}\int\mathrm{d}^{2}x\,\bar{\psi}(x)\psi(x)-\phi^{(\mathrm{f})}N_{f}\right).$
(5)

This type of constraint involves a combination of Grassmann numbers in a Dirac
delta distribution. A priori, it is not clear whether this is a well-defined
object. However, if we choose a representation of the Dirac delta as
$\delta(x)=\int_{-\infty}^{\infty}\tfrac{\mathrm{d}\eta}{2\uppi}\,\mathrm{e}^{\mathrm{i}\eta
x},$ one recognizes that it merely requires that arbitrary powers of $x$ are
non-vanishing. This is realized if $\bar{\psi}$ and $\psi$ have infinitely
many degrees of freedom as, e.g., in the infinite volume, and thus arbitrary
powers of $\tfrac{1}{V}\int\mathrm{d}x\,\bar{\psi}(x)\psi(x)$ do not vanish.
This aspect will be discussed in greater detail in Ref. [5].

Applying this expansion to the path integral with the fermionic constraint and
a subsequent integration over the fermionic degrees of freedom results in the
reduced path integral

$\displaystyle\mathcal{Z}_{\phi^{(\mathrm{f})}}=\int\\!\\!\mathcal{D}\sigma\int_{-\infty}^{\infty}\tfrac{\mathrm{d}\\!\eta}{2\uppi}\
\mathrm{e}^{-S^{\prime}_{\mathrm{eff}}(\eta)[\sigma]-\mathrm{i}\eta\phi^{(\mathrm{f})}N_{f}}\,,\quad
S^{\prime}_{\mathrm{eff}}(\eta)[\sigma]=-\ln\mathrm{Det}\left(Q\big{[}\sigma\big{]}+\tfrac{\mathrm{i}\eta}{V}\right)+\int\mathrm{d}^{2}x\,\tfrac{N_{f}\sigma^{2}}{2g^{2}}.$
(6)

One cannot perform the $\eta$-integral analytically in this form. We can,
however, expand the determinant in powers of $1/V$ to obtain

$\displaystyle
S^{\prime}_{\mathrm{eff}}(\eta)[\sigma]=-\ln\mathrm{Det}\left(Q[\sigma]\right)-\tfrac{\mathrm{i}\eta}{V}\mathrm{Tr}\,Q^{-1}-\tfrac{\eta^{2}}{2V^{2}}\mathrm{Tr}\,Q^{-2}+\mathcal{O}\\!\left(\tfrac{\eta^{3}}{V^{3}}\right)+\int\mathrm{d}^{2}x\,\tfrac{N_{f}\sigma^{2}}{2g^{2}}.$
(7)

Neglecting terms of order $\mathcal{O}\left(\tfrac{1}{V^{3}}\right)$, allows
us to carry out the $\eta$-integral in a closed form to obtain the constraint
path integral

$\displaystyle\mathcal{Z}_{\phi^{(\mathrm{f})}}=\mathcal{N}\int\\!\\!\mathcal{D}\sigma\,\mathrm{e}^{-S_{\mathrm{eff}}[\sigma]-\tfrac{N_{f}\left(\phi^{(\mathrm{f})}-\mathcal{M}\right)^{2}}{2\chi}-\tfrac{1}{2}\ln\mathrm{Tr}\,Q^{-2}},\quad\mathcal{M}\equiv\tfrac{1}{VN_{f}}\mathrm{Tr}\,Q^{-1},\quad\chi\equiv-\tfrac{1}{VN_{f}}\mathrm{Tr}\,Q^{-2},$
(8)

where $\mathcal{N}$ is an irrelevant normalization constant and the additional
$\ln\mathrm{Tr}\,Q^{-2}$ term is suppressed in the limit $N_{f}\to\infty$. The
relevant configurations $\Sigma$ in this limit are the ones that minimize the
complete exponent of the weight instead of only $S_{\mathrm{eff}}$. One can
intuitively understand that the modifications to the weight of the path
integral in Eq. 8 punish configurations, where $\mathcal{M}$ deviates from
$\phi^{(\mathrm{f})}$. As the expectation value of $\mathcal{M}$ would be the
chiral condensate in the unconstrained case, this result appears plausible.

In the presence of the constraint, it is, however, important to distinguish
${\cal M}$ from the chiral condensate, which reads then

$\displaystyle\frac{\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{f})}}}{VN_{f}}=\phi^{(\mathrm{f})}+\bigg{\langle}\frac{(\phi^{(\mathrm{f})}-\mathcal{M})^{2}}{\chi}\frac{\gamma}{\chi}\bigg{\rangle}_{\phi^{(\mathrm{f})}}-\frac{1}{VN_{f}}\bigg{\langle}\frac{\gamma}{\chi}\bigg{\rangle}_{\phi^{(\mathrm{f})}},\quad\gamma\equiv\tfrac{1}{VN_{f}}\mathrm{Tr}\,Q^{-3}.$
(9)

This might be surprising, since from Eq. 5, one would have expected the chiral
condensate to be $\phi^{(\mathrm{f})}$. The volume expansion destroys this
exactness and one retains additional terms as modifications, which stem from
the exchange of the infinite volume limit and the $\eta$-integration. However,
we will see in Section 4 that these corrections are actually small in the
relevant regions of $\phi^{(\mathrm{f})}$. Moreover, one finds that the Ward
identity 2 is preserved even after the volume expansion.

###  3.2 Bosonic constraint

The bosonic constraint, introduces a constraint on the spatial average of the
field $\sigma$ to the path integral as

$\displaystyle\mathcal{Z}_{{\phi^{(\mathrm{b})}}}=\int\mathcal{D}\big{[}\bar{\psi},\psi,\sigma\big{]}\
\mathrm{e}^{-S_{\sigma}}\,\delta\\!\left(\tfrac{1}{V}\int\mathrm{d}x\,\sigma(x)-\phi^{(\mathrm{b})}\right).$
(10)

Even though, we do not constrain the chiral condensate directly, one might
expect that it is constrained nevertheless due to the Ward identity 2.
However, one finds that this Ward identity receives an extra contribution in
the constrained path integral with the modified relation given by

$\displaystyle\frac{\langle\bar{\psi}\psi\rangle_{\phi^{(\mathrm{b})}}}{V}=\frac{N_{f}}{g^{2}}\phi^{(\mathrm{b})}-\frac{\partial\Omega_{\mathrm{b}}\left(\phi^{(\mathrm{b})}\right)}{\partial\phi^{(\mathrm{b})}}\equiv\frac{N_{f}}{g^{2}}\phi^{(\mathrm{b})}-j_{\mathrm{b}}\left(\phi^{(\mathrm{b})}\right),$
(11)

with $\Omega_{\mathrm{b}}=-\ln\mathcal{Z}_{\phi^{(\mathrm{b})}}$ Thus, we
expect that this bosonic constraint constrains the chiral condensate up to
$j_{\mathrm{b}}$, which is the derivative of the constraint potential with
respect to the constraint parameter. Moreover, in contrast to the fermionic
constraint, one finds that
$\langle\bar{\psi}\psi\rangle_{\phi^{(\mathrm{b})}}=V\bigl{<}\mathcal{M}\bigr{>}_{\phi^{(\mathrm{b})}}$
without any corrections. Bosonic constraints have been considered in various
studies [2, 3] and their application has been well understood. We exploit this
type of constraint in the special setup of the GN model to benchmark the
ability of the fermionic constraint to constrain the chiral condensate, which
we focus on in this work. As the bosonic constraint works on the level of the
bosonic fields, which are the degrees of freedom in which the minimization of
the action is done, we can enforce it numerically in the minimization routine.
The same is true for a Monte-Carlo simulation, where the $\sigma$-field would
be the degree of freedom that is integrated over.

##  4 Results of the constrained Gross-Neveu model

As a first test of the constraining methods, we consider a single coupling
$g^{2}\approx 0.49203$ on two lattices $V=L^{2}\in\\{40^{2},60^{2}\\}$. For
this coupling, the bosonic vacuum expectation value with the unconstrained
path integral is $\langle\sigma\rangle=\sigma_{0}=0.5=1.01620\,g^{2}$. In
order to facilitate the numerical minimizations, we restrict the
$\sigma$-field to be homogeneous in one of the two directions.

We start our investigation by calculating
$\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{x})}}$ for both constraints as
shown in Fig. 1 to estimate the extent of the consistency of the constraints.
With the bosonic constraint, we find that the chiral condensate is given by
$\phi^{(\mathrm{b})}/g^{2}$ up to a maximum value in $\phi^{(\mathrm{b})}$,
where the additional term in Eq. 11 contributes, i.e., $j\neq 0$. The chiral
condensate being constrained that strongly in some range of
$\phi^{(\mathrm{b})}$ even though we constrained the bosonic field, is
certainly remarkable. One finds that this maximum value of
$\phi^{(\mathrm{b})}$ increases with increasing volume. We find a similar
situation with the fermionic constraint, where the chiral condensate is
constrained again up to a maximum value in $\phi^{(\mathrm{f})}$, which
increases also with volume. However, the deviation from $\phi^{(\mathrm{f})}$
is minimal beyond this point and the chiral condensate is still constrained
approximately. This shows that the chiral condensate is controlled
considerably better by the fermionic constraint compared to the bosonic
constraint. For even larger $\phi^{(\mathrm{f})}$ (beyond $\sigma_{0}/g^{2}$),
one finds that there is a large deviation of
$\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{f})}}$ from
$\phi^{(\mathrm{f})}$, which appears to be independent of the volume. This
signals that the additional terms in Eq. 9 become relevant.

Figure 1: The expectation value
$\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{x})}}$ for the two constraint
methods.

To understand how the constraints works and why the jumps at the intermediate
values of $\phi^{(\mathrm{x})}$ happen, we consider the field configurations
$\Sigma$ that we find in the numerical minimization for $L=60$ in Fig. 2. We
show the configurations for both constraint types each at similar values of
$\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{x})}}$. For small values of
$\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{x})}}$, we find inhomogeneous
configurations in both cases, which interpolate between the minima of the
classical double well potential at $\pm\sigma_{0}=\pm 0.5$. The bosonic field
$\Sigma$ mostly stays in one of the two minima and transitions through the
unfavored regions with a tunneling profile that is independent of
$\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{x})}}$. For increasing value
of $\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{x})}}$, the ratio of space
that the bosonic field spends in one of the minima changes to accommodate the
change in the constraint. At a certain value of
$\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{x})}}$, the required ratio
cannot be fulfilled anymore leading to the observed breakdown of the
constraint. The approximate nature of the fermionic constraint enables
configurations such as the green curve, where the field does not tunnel all
the way. These configurations correspond to a value of the chiral condensate
just before the jump at $\phi^{(\mathrm{f})}\approx 0.75$ of the $L=60$ curve
in the right plot of Fig. 1. For large
$\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{x})}}$, one finds homogeneous
field configurations for both constraining methods. We find that the tunneling
profile has a fixed physical width and, thus, for larger volumes more extreme
ratios of space in the two minima can be fulfilled.

Figure 2: The minimizing field configuration $\Sigma$ for several values of
$\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{x})}}$ with the bosonic and
fermionic constraint for $L=60$.

Now we turn to the constraint potential
$\Omega_{\mathrm{x}}=-\ln\mathcal{Z}_{\phi^{(\mathrm{x})}}/V$. In a full
Monte-Carlo simulation, one would measure $j_{x}$ – the derivative of
$\Omega_{\mathrm{x}}$ w.r.t. to the constraint parameter – and reconstruct
$\Omega_{\mathrm{x}}$ from this [3]. In the considered $N_{f}\to\infty$ limit,
however, we can also directly calculate $\Omega_{\mathrm{x}}$, because there
is only a single configuration that is contained in the path integral. The two
methods yield the same result within numerical accuracy. Fig. 3 depicts the
constraint potential for both constraining methods as a function of the
constraint parameter. Both methods feature a flat region for small constraint
parameters, which extends to larger constraint parameters for larger volumes.
This regime features the inhomogeneous field configurations shown in Fig. 2.
The value of the flat region is finite but approaches the value corresponding
to the minimum for increasing volume. The flat region will extend all the way
to the minimum in the infinite volume limit. This is the behavior that we
expect from the quantum effective potential, which agrees with the constraint
potential for $V\to\infty$. We notice that the deviation from the flat
behavior is non-analytic for the bosonic constraint and smooth for the
fermionic constraint, which reflects the approximate nature of this
constraint.

Figure 3: The constraint potential $\Omega_{\mathrm{x}}$ as a function of the
constraint parameter with the bosonic and fermionic constraint for
$L\in\\{40,60\\}$.  Figure 4: The constraint potential $\Omega_{\mathrm{x}}$
as a function of $\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{x})}}$ with
both constraint types for $L\in\\{40,60\\}$.

Next, we consider the constraint potential as a function of
$\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{x})}}$ instead of the
constraint parameter, as this enables a real comparison of the two methods.
This is done in Fig. 4, where the results in the flat region from different
constraint methods for the same volume agree within numerical accuracy. For
the fermionic constraint, one finds that $\Omega_{\mathrm{f}}$ is double-
valued for some values of
$\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{f})}}$. There is, however, no
stability issue as this is only due to the fact that two values of
$\phi^{(\mathrm{f})}$ can correspond to a single value of
$\bigl{<}\bar{\psi}\psi\bigr{>}_{\phi^{(\mathrm{f})}}$ in the region where it
jumps in Fig. 1. We conclude that the fermionic constraint is actually able to
constrain the chiral condensate in order to calculate the constraint
potential. As expected, this flattens for increasing volume and we expect it
to be equivalent to the quantum effective potential in the infinite volume
limit.

##  5 Conclusion

In this work, we presented two methods to constrain the chiral condensate on
the level of the path integral and used these to successfully calculate the
constraint potential as a function of the expectation value
$\langle\bar{\psi}\psi\rangle$. While the bosonic constraint was used in this
form in others contexts, the fermionic constraint is an entirely new approach,
which directly constrains the fermion bilinear $\bar{\psi}\psi$. Neither of
the methods are limited to the $N_{f}\to\infty$ limit in which we tested the
methods, but could also be used in Monte-Carlo simulations. It was shown that
such constrained simulations are an excellent tool to investigate spontaneous
symmetry breaking on the lattice with a lesser numerical effort compared to
the traditional method with an explicit breaking [3]. Moreover, the fermionic
constraint is not unique to the GN model and can be implemented similarly in
other fermionic theories. This includes QCD, where it could provide a valuable
alternative to the standard approach in full lattice simulations. This
discussion and results of the constrained chiral GN model featuring a
continuous chiral symmetry will be presented in Ref. [5].

## Acknowledgments

We acknowledge useful discussions with Bastian Brandt. We acknowledge the
contribution of Marc Winstel to the code framework that was used for the
numerical minimization of the constrained Gross-Neveu model. G.E., G.M. and
L.P. acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) through the collaborative research center trans-
regio CRC-TR 211 “Strong-interaction matter under extreme conditions”– project
number 315477589 – TRR 211. L.P. acknowledges the support of the Helmholtz
Graduate School for Hadron and Ion Research and of the Giersch Foundation.

## References

  * [1] L. O’Raifeartaigh, A. Wipf and H. Yoneyama, _The constraint effective potential_ , _Nucl. Phys. B_ 271 (1986) 653.
  * [2] Z. Fodor, K. Holland, J. Kuti, D. Nogradi and C. Schroeder, _New Higgs physics from the lattice_ , _PoS_ LATTICE2007 (2007) 056.
  * [3] G. Endrődi, T.G. Kovács and G. Markó, _Spontaneous Symmetry Breaking via Inhomogeneities and the Differential Surface Tension_ , _Phys. Rev. Lett._ 127 (2021) 232002.
  * [4] D.J. Gross and A. Neveu, _Dynamical symmetry breaking in asymptotically free field theories_ , _Phys. Rev. D_ 10 (1974) 3235.
  * [5] G. Endrődi, T.G. Kovács, G. Markó and L. Pannullo, _Constraining fermionic condensates_ , in preparation.
  * [6] A. Koenigstein, L. Pannullo, S. Rechenberger, M. Winstel and M.J. Steil, _Detecting inhomogeneous chiral condensation from the bosonic two-point function in the (1 + 1)-dimensional Gross-Neveu model in the mean-field approximation_ , _J. Phys. A: Math. Theor._ 55 (2022) 375402.
  * [7] J. Lenz, L. Pannullo, M. Wagner, B. Wellegehausen and A. Wipf, _Inhomogeneous phases in the Gross-Neveu model in 1+1 dimensions at finite number of flavors_ , _Phys. Rev. D_ 101 (2020) 094512.