# The Effect of Geometry on Survival and Extinction in a Moving-Boundary
Problem Motivated by the Fisher–KPP Equation

Alexander K. Y. Tam  Matthew J. Simpson  School of Mathematical Sciences,
Queensland University of Technology, Brisbane QLD 4000, Australia.

###### Abstract

The Fisher–Stefan model involves solving the Fisher–KPP equation on a domain
whose boundary evolves according to a Stefan-like condition. The Fisher–Stefan
model alleviates two practical limitations of the standard Fisher–KPP model
when applied to biological invasion. First, unlike the Fisher–KPP equation,
solutions to the Fisher–Stefan model have compact support, enabling one to
define the interface between occupied and unoccupied regions unambiguously.
Second, the Fisher–Stefan model admits solutions for which the population
becomes extinct, which is not possible in the Fisher–KPP equation. Previous
research showed that population survival or extinction in the Fisher–Stefan
model depends on a critical length in one-dimensional Cartesian or radially-
symmetric geometry. However, the survival and extinction behaviour for general
two-dimensional regions remains unexplored. We combine analysis and level-set
numerical simulations of the Fisher–Stefan model to investigate the
survival–extinction conditions for rectangular-shaped initial conditions. We
show that it is insufficient to generalise the critical length conditions to
critical area in two-dimensions. Instead, knowledge of the region geometry is
required to determine whether a population will survive or become extinct.

###### keywords:

reaction–diffusion; Fisher–Stefan; level-set method; spreading–vanishing
dichotomy

##  1 Introduction

The dimensionless Fisher–KPP equation,

$\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}u}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}t}=\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}^{2}u}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}x^{2}}+u\left(1-u\right),$
(1.1)

is a classical prototype model in mathematical biology [1, 2, 3, 4] that
describes the spatiotemporal evolution of a population density, $u(x,t)>0,$
that evolves due to linear Fickian diffusion combined with a logistic source
term. A key property of the Fisher–KPP equation 1.1 is that it admits
travelling-wave solutions, with long-time speed $c=2$ for compactly-supported
initial conditions, when solved on an infinite domain. By relating population
invasion to travelling wave solutions, the Fisher–KPP equation and extensions
of this model have been used to represent species invasion in ecology [5, 6,
7], front propagation in chemical reactions [8], and biological cell invasion
[9, 10, 11, 12, 13, 14, 15]. However, despite its widespread usage, the
Fisher–KPP equation has practical limitations. Firstly, any initial condition
gives rise to population growth and complete colonisation. The Fisher–KPP
equation is thus unsuitable for populations where extinction [16] or arrested
invasion [17] is of interest. A second shortcoming is that solutions to 1.1 do
not have compact support. Thus, we cannot identify the interface between
occupied and unoccupied regions without ambiguity. This leads to difficulties
in applying the Fisher–KPP equation to processes with well-defined invasion
fronts, for example tumour cell invasion [18, 19] or wound healing [10, 20].

To address the shortcoming of non-compactly-supported solutions, a common
approach is to modify 1.1 to incorporate degenerate nonlinear diffusion [4,
21]

$\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}u}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}t}=\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}x}\left(D(u)\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}u}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}x}\right)+u\left(1-u\right),$
(1.2)

such that $D(0)=0.$ A common choice is $D(u)=u^{k}$ [22, 21] for $k>0.$
Setting $k=1$ leads to the _Porous-Fisher’s_ equation [21], and for $k>1$ 1.2
is sometimes called the _generalised Porous-Fisher’s_ equation [20, 23]. Like
the Fisher–KPP equation, the generalised Porous-Fisher’s equation admits
travelling-wave solutions. For the minimum wave speed, these travelling waves
are sharp-fronted and have compact support [24, 25], thus enabling one to
define an unambiguous front. However, although 1.2 admits compactly-supported
travelling-waves, like the Fisher–KPP equation it guarantees population
survival, and cannot be used to model population extinction or receding fronts
with local decrease in density [17, 16].

Another alternative to the standard Fisher–KPP equation is to consider a
moving-boundary problem, such that 1.1 holds on a compactly-supported region
$x<L(t),$ and $u(L(t),t)=0.$ The moving-boundary is assumed to evolve
according to

$\frac{\mathrm{d}L}{\mathrm{d}t}=-\kappa\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}u}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}x}\quad\text{
on }\quad x=L(t),$ (1.3)

where the parameter $\kappa$ relates the density gradient at $L(t)$ to the
speed of the boundary. Moving-boundary problems of this type are traditionally
used to model physical and industrial processes [26, 27, 28, 29, 30]. Indeed,
the boundary condition 1.3 is analogous to the classical Stefan condition [31]
for a material undergoing phase change, where $\kappa$ is the inverse Stefan
number. More recently, moving-boundary problems have been used to study
biological phenomena, including cancer invasion, cell motility, and wound
healing [32, 33, 34, 35, 36]. Applying the moving-boundary problem framework
to 1.1, we obtain

$\displaystyle\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}u}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}t}=\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}^{2}u}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}x^{2}}+u\left(1-u\right)\quad\text{
on }\quad 0<x<L(t),$ (1.4a)
$\displaystyle\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}u}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}x}=0\quad\text{
on }\quad x=0,$ (1.4b) $\displaystyle
u=0,\quad\frac{\mathrm{d}L}{\mathrm{d}t}=-\kappa\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}u}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}x}\quad\text{
on }\quad x=L(t),$ (1.4c) $\displaystyle u(x,0)=u_{0}(x)\quad\text{ on }\quad
0<x<L(0).$ (1.4d)

This extension to 1.1, known as the _Fisher–Stefan_ model, was first proposed
by Du and Lin [37], and Du and Guo [38, 39]. The Fisher–Stefan model
alleviates the two practical disadvantages of Fisher–KPP, because the boundary
$L(t)$ defines the front position explicitly, and it admits solutions for
population extinction [37, 38, 40].

The Fisher–Stefan model has been studied extensively in one-dimensional
Cartesian geometry. Du and Lin [37] proved that 1.4 admits solutions whereby
the population density evolves to a travelling wave solution as $t\to\infty,$
with asymptotic speed that depends on $\kappa.$ This corresponds to population
survival and successful invasion, with $u(x,t)\to 1$ as $t\to\infty.$ However,
Du and Lin [37] also proved that the Fisher–Stefan model admits solutions
where the population fails to establish, such that $u\to 0^{+}$ as
$t\to\infty.$ This corresponds to population extinction, which cannot occur in
the Fisher–KPP or generalised Porous-Fisher models. The survival–extinction
behaviour of the Fisher–Stefan model depends on the critical length,
$L_{c}=\pi/2$ [37, 41]. If $L(0)>L_{c},$ Du and Lin [37] showed that the
population will always survive. However, if $L(0)<L_{c},$ then survival only
occurs if $L(t)$ evolves such that $L(t)>L_{c}$ at some time, and otherwise
the population becomes extinct. In this $L(0)<L_{c}$ scenario, whether the
solution ever attains $L(t)>L_{c}$ depends on the initial condition
$u_{0}(x),$ and the parameter $\kappa$ [37].

Simpson [42] extended the work of Du and Lin [37] by considering the
Fisher–Stefan model on an $n$-sphere. By replacing the second derivative term
in 1.4a with the $n$-dimensional radially-symmetric Laplacian operator,
Simpson [42] showed that a critical radius, $R_{c},$ governs survival and
extinction analogously to the critical length. This critical radius depends on
the dimension $n.$ In the two-dimensional spreading-disc problem, the critical
radius is $R_{c}=\alpha_{01}\approx 2.4048,$ where $\alpha_{01}$ is the first
zero of $J_{0}(x),$ the zeroth-order Bessel function of the first kind [42].
However, the survival and extinction behaviour of the Fisher–Stefan model
remains unexplored for non-radially-symmetric two-dimensional geometries.
Investigating this is the subject of our study.

##  2 Results and Discussion

###  2.1 Mathematical Model

We consider a multidimensional extension to the dimensionless Fisher–Stefan
model 1.4. As illustrated in Figure 2.1A, we denote the region on which we
solve the Fisher–KPP equation as $\Omega(t),$ and its boundary as
$\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega(t).$
The boundary
$\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega(t)$
defines the interface between regions of non-zero population density and
regions of zero density. Biologically, this might represent the position of a
cell invasion front, or the boundary of a tumour. The multidimensional
Fisher–Stefan model is then

$\displaystyle\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}u}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}t}=\nabla^{2}u+u\left(1-u\right),\quad\text{
on }\quad\Omega(t),$ (2.1a) $\displaystyle u=0,\quad\text{ on
}\quad\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega(t),$
(2.1b) $\displaystyle V=-\kappa\nabla u\cdot\mathvec{\hat{n}}\quad\text{ on
}\quad\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega(t),$
(2.1c) $\displaystyle u(\mathvec{x},0)=u_{0}(\mathvec{x})\quad\text{ on
}\quad\Omega(0),$ (2.1d)

where $V$ is the normal speed of the interface, and $\mathvec{\hat{n}}$ is the
unit outward normal to the interface
$\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega(t).$

###  2.2 Numerical Methods

We use the level-set method [43, 44, 45, 46] to solve 2.1 on a two-dimensional
computational domain, $\mathcal{D}.$ This involves embedding the interface as
the zero level-set of a scalar function $\phi(\mathvec{x},t).$ That is,

$\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega(t)=\\{\mathvec{x}\mid\phi\left(\mathvec{x},t\right)=0\\},$
(2.2)

where $\phi(\mathvec{x},t)$ is defined for all $\mathvec{x}\in\mathcal{D},$
and is such that $\phi<0$ for $\mathvec{x}\in\Omega(t)$ and $\phi\geq 0$ for
$\mathvec{x}\notin\Omega(t).$ An example for a disc is illustrated in Figure
2.1B. The level-set function $\phi(\mathvec{x},t)$ evolves according to the
level-set equation

$\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\phi}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}t}+F\left\lvert\nabla\phi\right\rvert=0,$
(2.3)

where $F(\mathvec{x},t)$ is the _extension velocity field_ , a scalar function
defined for $\mathvec{x}\in\mathcal{D},$ such that $F=V$ on
$\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega(t).$
An example extension velocity field for a spreading disc is shown in Figure
2.1D. As $\phi$ evolves with time, so too does its zero level-set, which
defines the region on which we solve 2.1a. We solve the Fisher–Stefan model
using explicit finite-difference methods [47, 48, 49, 50].

Figure 2.1: Illustration of the level-set method used to solve the two-
dimensional Fisher–Stefan model 2.1. (A) Schematic of the rectangular
computational domain $\mathcal{D},$ which contains the blue region $\Omega(t)$
on which we solve the Fisher–KPP equation, and the red curve, which is the
boundary
$\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega(t).$
(B) An example level-set function, $\phi(x,y,t),$ for a disc-shaped
$\Omega(t).$ The black curve is the zero level-set, and the red curve is a
projection onto the $(x,y)$ plane marked with the grid. (C) An example two-
dimensional density profile $u(x,y,t)$ for a circular region. (D) An example
extension velocity field $F(x,y,t).$

For each numerical solution in this work, we use $201\times 201$ spatial
resolution, with $\Delta t=0.002,$ and $\kappa=0.1.$ This is sufficient to
obtain good agreement with grid-independent spreading disc solutions of
Simpson [42]. Full details on our level-set method, an open-source Julia
implementation, and numerical tests to validate the method are available on
Github.

###  2.3 Survival and Extinction in Two-Dimensions

We begin our investigation of survival and extinction in two-dimensions by
considering solutions with the disc initial condition,

$u(x,y,0)=\begin{cases}u_{0}&\quad\text{ if
}\quad\sqrt{\left(x-\frac{W}{2}\right)^{2}+\left(y-\frac{W}{2}\right)^{2}}\leq
R(0)\\\ 0&\quad\text{ otherwise,}\end{cases}$ (2.4)

where $u_{0}$ is a constant satisfying $0<u_{0}\leq 1,$ $R(0)$ is the initial
disc radius, and $W$ is the width of the square computational domain
$\mathcal{D}.$ Simpson [42] showed that the critical radius,
$R_{c}=\alpha_{01},$ explains the survival and extinction behaviour of a disc.
With the initial condition 2.4, the region $\Omega(t)$ evolves as a disc of
radius $R(t).$ If $R(t)>R_{c}$ at any time, the population survives.
Alternatively, if $R(t)<R_{c}$ for all $t,$ then the population eventually
becomes extinct. We solve the Fisher–Stefan model numerically using
$u_{0}=0.5,$ and compare solutions with $R(0)=2.1$ and $R(0)=2.4.$ Although,
both values of $R(0)$ are less than $R_{c},$ the Stefan-like condition 2.1c
suggests that $R(t)$ will increase for $\kappa>0,$ because $\nabla
u\cdot\hat{\mathvec{n}}<0$ on
$\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega(t).$
However, it is unclear whether $R(t)$ will ever exceed $R_{c},$ because the
rate of spread $V$ in 2.1c depends on $\kappa$ and the shape of the density
profile.

Level-set numerical solutions with $R(0)=2.1$ and $R(0)=2.4$ are presented in
Figure 2.2.

Figure 2.2: Numerical solutions to the Fisher–Stefan model with the disc
initial condition 2.4, and $u_{0}=0.5.$ (A–D): Solution with $R(0)=2.1.$ (A):
Initial condition, $t=0.$ (B): Density profile $u(x,y,t)$ at $t=5.$ (C):
Density profile $u(x,y,t)$ at $t=10.$ (D): Evolution of population density,
viewed as a one-dimensional slice through $y=10.$ The solid curve represents
the initial condition, and dashed curves are the densities at
$t\in\\{2,4,6,8,10\\}.$ Arrow indicates direction of increasing $t.$ (E–H):
solution with $R(0)=2.4.$ (E): Initial condition, $t=0.$ (F): Density profile
$u(x,y,t)$ at $t=50.$ (G): Density profile $u(x,y,t)$ at $t=100.$ (H):
Evolution of population density, viewed as a one-dimensional slice through
$y=10.$ The solid curve represents the initial condition, and dashed curves
are the densities at $t\in\\{20,40,60,80,100\\}.$ Arrow indicates direction of
increasing $t.$ (D&H): Red vertical line represents $x=10\pm R_{c},$ where
$R_{c}=\alpha_{01}.$

Panels A–D show that $R(0)=2.1$ results in population extinction, with density
$u\to 0^{+}$ as $t$ increases. This is because $R(t)<R_{c},$ shown as coloured
vertical lines in panel D. Eventually, population density $u\to 0^{+},$ which
prevents further spread because $\nabla u\to 0,$ and thus $V\to 0$ according
to 2.1c. In contrast, with $R(0)=2.4$ the population survives and spreads,
such that both density $u$ and the size of $\Omega(t)$ eventually increase
with time. Although the density $u(10,10,t)$ at the centre of $\Omega(t)$
decreases at first (see panel H), we eventually have $R(t)>R_{c}.$ This
enables the density to recover, and the population to survive. As
$t\to\infty,$ the solution approaches a radially-symmetric travelling wave.
The ability to capture both survival and extinction is an advantage of the
Fisher–Stefan model over the Fisher–KPP equation.

Despite being previously studied in radially-symmetric geometry, the survival
and extinction behaviour of the Fisher–Stefan model in general two-dimensional
geometry remains unexplored. For example, it is not immediately clear how the
results for a disc might translate to rectangular shapes. Given a disc of
radius $R_{c},$ one can conceive many rectangles that share certain properties
with the disc. Some examples are shown in Figure 2.3.

Figure 2.3: Potential square and rectangular analogues to the critical disc
with radius $R_{c}.$ (a) A square whose corners lie on the perimeter of the
circle. (b) A square with side width equal to the circle diameter. (c) A
square with the same area as the circle. (d) A rectangle with aspect ratio 4,
and the same area as the circle.

For example, in Figure 2.3(a) we have the square with width $\sqrt{2}R_{c},$
for which all points inside the square are within the critical disc. In Figure
2.3(b) we have the square with width $2R_{c},$ such that every point of the
critical disc is contained within the square. Based on the radially-symmetric
analysis, we would expect the initial square in Figure 2.3(a) to give rise to
extinction, and the initial square in Figure 2.3(b) to give rise to survival.
However, the situation is unclear for a square with side width between
$\sqrt{2}R_{c}$ and $2R_{c}.$ One example is the square with width
$\sqrt{\pi}R_{c}$ (Figure 2.3(c)), which has the same area as the critical
disc. In this case, consider drawing a family of rays from the centre of the
square or disc to the perimeter of each shape. Each ray is defined by its
polar angle $\theta\in[0,2\pi].$ For some $\theta,$ for example
$\theta=0,\pm\pi/2,$ $R_{c}$ exceeds the distance to the perimeter of the
square. In contrast, for $\theta=\pm\pi/4,\pm 3\pi/4,$ the distance to the
perimeter of the square exceeds $R_{c}.$ Therefore, it is unclear whether a
population of this shape will survive or become extinct. Furthermore, varying
the aspect ratio enables us to define a family of rectangles with the same
area, one of which is shown in Figure 2.3(d). The survival and extinction
behaviour of these shapes is also unclear according to previous radially-
symmetric analysis.

To address this question, we consider numerical solutions on initially
rectangular domains. We characterise the geometry of $\Omega(t)$ by defining
$L_{x}(t)$ and $L_{y}(t)$ to be the distances in the $x$ and $y$-directions
respectively occupied by the population, measured through the centre of
$\Omega.$ That is,

$\displaystyle L_{x}(t)$ $\displaystyle=\max\\{x\mid
u(x,W/2,t)>0\\}-\min\\{x\mid u(x,W/2,t)>0\\},$ (2.5a) $\displaystyle L_{y}(t)$
$\displaystyle=\max\\{y\mid u(W/2,y,t)>0\\}-\min\\{y\mid u(W/2,y,t)>0\\}.$
(2.5b)

For initially-rectangular domains, we then have the initial condition

$u(x,y,0)=\begin{cases}u_{0}&\quad\text{ if }\quad\left\lvert
x-\frac{W}{2}\right\rvert\leq\frac{L_{x}(0)}{2}\quad\text{ and
}\quad\left\lvert y-\frac{W}{2}\right\rvert\leq\frac{L_{y}(0)}{2},\\\
0&\quad\text{ otherwise.}\end{cases}$ (2.6)

Numerical results are shown in Figure 2.4. Panels A–E show that a solution
with the square initial condition $L_{x}(0)=L_{y}(0)=3.5$ leads to extinction.
Increasing the size of the initial square to $L_{x}(0)=L_{y}(0)=4$ yields
survival (Figure 2.4, panels F–J). A rectangular initial condition with
$L_{x}(0)=8$ and $L_{y}(0)=2$ leads to extinction, as panels K–O of Figure 2.4
show. This occurs despite this rectangle having the same initial area as the
square that survived. Thus, the survival or extinction of initially-
rectangular regions does not depend solely on area, and more detail about the
region geometry is required.

Figure 2.4: Numerical solutions with rectangular initial conditions 2.6, with
$u_{0}=0.5.$ (A–E): Solution with $L_{x}(0)=L_{y}(0)=3.5.$ (F–J): Solution
with $L_{x}(0)=L_{y}(0)=4.$ (K–O): Solution with $L_{x}(0)=8,$ $L_{y}(0)=2.$
(A): Initial density $u(x,y,0).$ (B): Density profile at $t=2.$ (C): Density
profile at $t=4.$ (D–E) One-dimensional slice $u(x,10,t).$ $u(10,y,t).$ Solid
curve is the initial condition, dashed curves are solutions plotted at
$t\in\\{1,2,3,4,5\\}.$ Arrow indicates direction of increasing $t.$ (F):
Initial density $u(x,y,0).$ (G): Density profile at $t=50.$ (H): Density
profile at $t=100.$ (I–J) One-dimensional slice $u(x,10,t).$ $u(10,y,t).$
Solid curve is the initial condition, dashed curves are solutions plotted at
$t\in\\{20,40,60,80,100\\}.$ Arrow indicates direction of increasing $t.$ (K):
Initial density $u(x,y,0).$ (L): Density profile at $t=1.$ (M): Density
profile at $t=2.$ (N–O) One-dimensional slice $u(x,10,t).$ $u(10,y,t).$ Solid
curve is the initial condition, dashed curves are solutions plotted at
$t\in\\{1,2,3,4,5\\}.$ Arrow indicates direction of increasing $t.$

To provide insight into the results in Figure 2.4, we perform similar analysis
to El-Hachem _et al._ [41] and Simpson [42]. This involves considering the
limit $u\to 0^{+},$ and obtaining critical conditions that determine survival
or extinction. In this limit, the leading-order population density
$\hat{u}(x,y,t)$ satisfies

$\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\hat{u}}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}t}=\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}^{2}\hat{u}}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}x^{2}}+\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}^{2}\hat{u}}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}y^{2}}+\hat{u}.$
(2.7)

We consider solutions to 2.7 on a fixed rectangular domain $0\leq x\leq X$ and
$0\leq y\leq Y,$ subject to arbitrary compactly-supported initial conditions,
and homogeneous Dirichlet boundary conditions on all boundaries. This problem
has the general solution

$\hat{u}(x,y,t)=\sum_{n=1}^{\infty}\sum_{m=1}^{\infty}A_{n,m}\sin\left(\frac{n\pi
x}{X}\right)\sin\left(\frac{m\pi
y}{Y}\right)\mathrm{e}^{-\left(\frac{n^{2}\pi^{2}}{X^{2}}+\frac{m^{2}\pi^{2}}{Y^{2}}-1\right)t},$
(2.8)

with coefficients $A_{n,m}$ chosen such that $\hat{u}(x,y,0)$ matches the
initial condition. We next consider the long-time solution $t\to\infty,$ such
that the leading-eigenvalue $n=m=1$ approximates the general solution 2.8.
Then,

$\hat{u}(x,y,t)\thicksim A_{1,1}\sin\left(\frac{\pi
x}{X}\right)\sin\left(\frac{\pi
y}{Y}\right)\mathrm{e}^{-\left(\frac{\pi^{2}}{X^{2}}+\frac{\pi^{2}}{Y^{2}}-1\right)t}\quad\text{
as }\quad t\to\infty.$ (2.9)

To obtain the conditions under which the population survives or becomes
extinct as $t\to\infty,$ we use the approximation 2.9 and consider
conservation of the population inside $\Omega.$ For $u\ll 1,$ the conservation
statement is

$\frac{\mathrm{d}M}{\mathrm{d}t}=\int_{\Omega}\hat{u}(x,y,t)-\int_{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega}-\nabla\hat{u}\cdot\mathvec{\hat{n}},$
(2.10)

where $M(t)=\int_{\Omega}\hat{u}(x,y,t)$ is the total population in
$\Omega(t).$ The first term on the right-hand side of 2.10 is the net
accumulation of $\hat{u}$ inside $\Omega(t)$ due to the (linearised) source
term. The second term on the right-hand side of 2.10 is the rate of population
loss through the boundary
$\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega(t)$
due to diffusion. The population survives if $\mathrm{d}M/\mathrm{d}t>0$ as
$t\to\infty,$ that is accumulation exceeds loss. Alternatively, the population
becomes extinct if $\mathrm{d}M/\mathrm{d}t<0$ as $t\to\infty,$ _i.e._ the
rate of loss exceeds accumulation. Critical behaviour occurs when
$\mathrm{d}M/\mathrm{d}t=0$ as $t\to\infty,$ such that the rates of population
gain and loss balance. This requires

$\int_{\Omega}\hat{u}(x,y,t)=\int_{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega}-\nabla\hat{u}\cdot\mathvec{\hat{n}}.$
(2.11)

In fixed rectangular geometry, $\Omega$ and
$\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega$
are readily parameterised to give

$\begin{gathered}\int_{0}^{Y}\int_{0}^{X}\hat{u}(x,y,t)\,\mathrm{d}x\,\mathrm{d}y=-\int_{0}^{Y}\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\hat{u}(X,y,t)}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}x}\,\mathrm{d}y\\\
+\int_{0}^{Y}\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\hat{u}(0,y,t)}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}x}\,\mathrm{d}y-\int_{0}^{X}\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\hat{u}(x,Y,t)}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}y}\,\mathrm{d}x+\int_{0}^{X}\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\hat{u}(x,0,t)}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}y}\,\mathrm{d}x.\end{gathered}$
(2.12)

Substituting the leading-eigenvalue approximation 2.9 for $\hat{u}$ in 2.12,
evaluating integrals, and simplifying then yields

$XY=\pi\sqrt{Y^{2}+X^{2}}.$ (2.13)

The condition 2.13 implies that the critical area for a rectangle with
dimensions $L_{x}\times L_{y}$ is $A_{c}=\pi\sqrt{L_{x}^{2}+L_{y}^{2}}.$ The
analysis suggests that a rectangular population will survive if its lengths
$L_{x}(t)$ and $L_{y}(t)$ ever satisfy $L_{x}L_{y}>A_{c},$ and will become
extinct if $L_{x}L_{y}<A_{c}$ for all time. Another interpretation of the
condition 2.13 is that extinction will occur if the exponential term in 2.9
decays, and the population will survive if the exponential term grows. Setting
the exponent in the leading-eigenvalue approximation 2.9 to zero then yields
the same condition 2.13.

Based on our analysis, the inequalities

$L_{y}>\pi\sqrt{\frac{L_{x}^{2}}{L_{x}^{2}-\pi^{2}}},\quad L_{x}>\pi,$ (2.14)

must be satisfied at some $t$ for population survival. Interestingly, the
population cannot survive if either $L_{x}$ or $L_{y}$ remain less than $\pi$
for all time. Furthermore, for each $L_{x}>\pi,$ there exists a unique
critical $L_{y,c}$ given by the right-hand side of the first inequality in
2.14, above which the population will survive. This completely characterises
the effect of rectangular geometry and aspect ratio on survival and extinction
in the two-dimensional Fisher–Stefan model. The conditions 2.14 enable
comparison with the radially-symmetric results of Simpson [42]. The critical
disc with radius $R_{c}=2.4048$ has area $A=18.17.$ By comparison, the
critical square has widths $L_{x}=L_{y}=\sqrt{2}\pi\approx 4.443,$ and area
$A=19.74.$

Our analysis explains the numerical solutions in Figure 2.4. For the square
solution in Figure 2.4, panels A–E with $L_{x}(0)=L_{y}(0)=3.5,$ the widths
never reach $\sqrt{2}\pi,$ and thus the population eventually becomes extinct.
However, with $L_{x}(0)=L_{y}(0)=4,$ (Figure 2.4, panels F–J), the square
widths eventually exceed $\sqrt{2}\pi,$ leading to eventual survival. For a
rectangle with $L_{x}(0)=8$ and $L_{y}(0)=2,$ the population becomes extinct
because $L_{y}(t)<\pi$ for all $t,$ which guarantees extinction. These and
other results are summarised in Figure 2.5, which illustrates how numerical
solutions with rectangular initial conditions evolve. If the inequalities 2.14
holds for any $t,$ the population survives and $L_{x}$ and $L_{y}$ continue to
increase. This is shown by the trajectories above the red curve in Figure 2.5.
Conversely, the population becomes extinct if 2.14 never holds. This occurs
for the trajectories that remain below the red curve in Figure 2.5.

Figure 2.5: Evolution of numerical solutions to the Fisher–Stefan model with
rectangular initial conditions. Black dots represent the initial rectangle
dimensions $L_{x}(0),$ and $L_{y}(0)$ for each solution plotted. Dashed curves
represent the evolution of $L_{x}(t)$ and $L_{y}(t)$ from $t=0$ to $t=10$ for
extinction solutions, or $t=0$ to $t=100$ for survival solutions. The red
curve defines the critical threshold between survival and extinction, as per
2.14.

The critical conditions 2.14 hold provided the domain $\Omega(t)$ remains
approximately rectangular. Our choice of $\kappa=0.1$ in all numerical
solutions of Figures 2.4 and 2.5 ensures that the rectangular shapes persist
sufficiently long for 2.14 to determine survival or extinction. With larger
$\kappa,$ for example $\kappa=1,$ initially-rectangular shapes satisfying the
survival conditions rapidly evolve to an expanding disc. In these scenarios,
the long-term survival–extinction behaviour instead obeys the critical radius
condition of Simpson [42].

##  3 Conclusion and Future Work

The Fisher–Stefan model provides an alternative to standard the Fisher–KPP
model that both defines an unambiguous invasion front, and admits solutions
whereby the population becomes extinct. In this study, we extended one-
dimensional and radially-symmetric results for the Fisher–Stefan model to two-
dimensional rectangular domains. Using the level-set method, we computed
numerical solutions to the two-dimensional Fisher–Stefan model with
$\kappa=0.1,$ for radially-symmetric and rectangular domains. Our objective
was to understand the conditions under which square and rectangular shapes
give rise to survival or extinction.

To explain the numerical observations, we obtained conditions for the
Fisher–Stefan model that govern survival and extinction of a rectangular-
shaped population. Using a conservation argument in the small density limit,
we showed that survival requires both lengths of the rectangle to exceed
$\pi,$ and for the rectangle area to exceed a critical threshold that depends
on the side lengths. Extending the one-dimensional idea of critical length or
critical radius to a critical area in two dimensions is insufficient to
explain survival and extinction. Instead, for rectangular-shaped domains
information about both widths is necessary. Interestingly, these findings are
similar to the behaviour of a population subject to the strong Allee effect.
In radially-symmetric geometry, Lewis and Kareiva [51] showed that a critical
radius governs survival and extinction. However, recent work by Li _et al._
[52] showed that additional geometric information is required to determine
survival or extinction in asymmetric two-dimensional rectangular domains.

The analytical methods established in this work apply directly to any domain
$\Omega(0)$ such that solution to 2.7 with homogeneous Dirichlet conditions on
$\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega(0)$
can be obtained using separation of variables. One avenue for future work is
to investigate the conditions for survival and extinction in these separable
shapes, for example an ellipse [53]. Furthermore, our numerical method can be
used to study survival and extinction for general initial shapes $\Omega(0),$
regardless of whether these are separable. This includes investigation of non-
simply-connected initial shapes.

All numerical solutions presented in this work have $\kappa=0.1.$ When solving
the Fisher–Stefan model numerically with larger $\kappa,$ for example
$\kappa=1,$ we observed that initially square and rectangular domains that
satisfied the conditions for survival quickly evolved to an expanding circular
shape. This numerical evidence suggests that radially-symmetric survival
solutions to the Fisher–Stefan model are stable with respect to small-
amplitude azimuthal perturbations. However, this has not been verified
analytically. The stability of inward-moving fronts in a hole-closing geometry
[20] is also yet to be analysed. We plan to investigate this stability problem
using analytical and numerical methods in future work.

## Acknowledgements

Computational resources and services used in this work were provided by HPC
and Research Support Group, Queensland University of Technology, Brisbane,
Australia. M. J. S. also acknowledges funding from the Australian Research
Council (Grant No. DP200100177).

## Author Contributions

M. J. S. designed the research. A. K. Y. T. wrote the numerical code, obtained
the results, performed the analysis, and wrote the manuscript, with guidance
from M. J. S.

## Competing Interests

We have no competing interests to declare.

## References

  * [1] A. Kolmogorov, I. Petrovsky, N. Piskunov, A study of the equation of diffusion with increase in the quantity of matter, and its application to a biological problem, Bulletin of Moscow State University Series A: Mathematics and Mechanics 1 (1937) 1–25.
  * [2] J. Canosa, On a nonlinear diffusion equation describing population growth, IBM Journal of Research and Development 17 (4) (1973) 307–313.  doi:10.1147/rd.174.0307.
  * [3] P. Grindrod, Patterns and Waves: The Theory and Applications of Reaction–Diffusion Equations, Oxford University Press, 1991.
  * [4] J. D. Murray, Mathematical Biology I: An Introduction, 3rd Edition, Springer, 2002.  doi:10.1007/b98868.
  * [5] B. H. Bradshaw-Hajek, P. Broadbridge, A robust cubic reaction–diffusion system for gene propagation, Mathematical and Computer Modelling 39 (9) (2004) 1151–1163.  doi:10.1016/S0895-7177(04)90537-7.
  * [6] J. G. Skellam, Random dispersal in theoretical populations, Bulletin of Mathematical Biology 53 (1) (1991) 135–165.  doi:10.1016/S0092-8240(05)80044-8.
  * [7] N. Shigesada, K. Kawasaki, Y. Takeda, Modeling stratified diffusion in biological invasions, The American Naturalist 146 (2) (1995) 229–251.  doi:10.1086/285796.
  * [8] G. N. Mercer, R. O. Weber, Combustion wave speed, Proceedings of the Royal Society of London Series A: Mathematical and Physical Sciences 450 (1938) (1995) 193–198.  doi:10.1098/rspa.1995.0079.
  * [9] R. A. Gatenby, E. T. Gawlinski, A reaction–diffusion model of cancer invasion, Cancer Research 56 (24) (1996) 5745–5753.
  * [10] P. K. Maini, D. L. S. McElwain, D. I. Leavesley, Traveling wave model to interpret a wound-healing cell migration assay for human peritoneal mesothelial cells, Tissue Engineering 10 (3-4) (2004) 475–482.  doi:10.1089/107632704323061834.
  * [11] B. G. Sengers, C. P. Please, R. O. C. Oreffo, Experimental characterization and computational modelling of two-dimensional cell spreading for skeletal regeneration, Journal of The Royal Society Interface 4 (17) (2007) 1107–1117.  doi:10.1098/rsif.2007.0233.
  * [12] J. A. Sherratt, J. D. Murray, Models of epidermal wound healing, Proceedings of the Royal Society of London Series B: Biological Sciences 241 (1300) (1990) 29–36.  doi:10.1098/rspb.1990.0061.
  * [13] M. J. Simpson, K. K. Treloar, B. J. Binder, P. Haridas, K. J. Manton, D. I. Leavesley, D. L. S. McElwain, R. E. Baker, Quantifying the roles of cell motility and cell proliferation in a circular barrier assay, Journal of The Royal Society Interface 10 (82) (2013) 20130007.  doi:10.1098/rsif.2013.0007.
  * [14] K. K. Treloar, M. J. Simpson, D. L. S. McElwain, R. E. Baker, Are _in vitro_ estimates of cell diffusivity and cell proliferation rate sensitive to assay geometry?, Journal of Theoretical Biology 356 (2014) 71–84.  doi:10.1016/j.jtbi.2014.04.026.
  * [15] S. T. Johnston, E. T. Shah, L. K. Chopin, D. L. S. McElwain, M. J. Simpson, Estimating cell diffusivity and cell proliferation rate by interpreting IncuCyte ZOOM™ assay data using the Fisher–Kolmogorov model, BMC Systems Biology 9 (1) (2015) 38\.  doi:10.1186/s12918-015-0182-y.
  * [16] M. El-Hachem, S. W. McCue, M. J. Simpson, Invading and receding sharp-fronted travelling waves, Bulletin of Mathematical Biology 83 (4) (2021) 35.  doi:10.1007/s11538-021-00862-y.
  * [17] K. A. Landman, G. J. Pettet, D. F. Newgreen, Mathematical models of cell colonization of uniformly growing domains, Bulletin of Mathematical Biology 65 (2) (2003) 235–262.  doi:10.1016/S0092-8240(02)00098-8.
  * [18] K. R. Swanson, C. Bridge, J. D. Murray, E. C. Alvord Jr, Virtual and real brain tumors: Using mathematical modeling to quantify glioma growth and invasion, Journal of the Neurological Sciences 216 (1) (2003) 1–10.  doi:10.1016/j.jns.2003.06.001.
  * [19] V. M. Pérez-García, G. F. Calvo, J. Belmonte-Beitia, D. Diego, L. Pérez-Romasanta, Bright solitary waves in malignant gliomas, Physical Review E 84 (2) (2011) 021921.  doi:10.1103/PhysRevE.84.021921.
  * [20] S. W. McCue, J. Wang, T. J. Moroney, K. Lo, S. Chou, M. J. Simpson, Hole-closing model reveals exponents for nonlinear degenerate diffusivity functions in cell biology, Physica D: Nonlinear Phenomena 398 (2019) 130–140.  doi:10.1016/j.physd.2019.06.005.
  * [21] T. P. Witelski, Merging traveling waves for the porous-Fisher’s equation, Applied Mathematics Letters 8 (4) (1995) 57–62.  doi:10.1016/0893-9659(95)00047-T.
  * [22] W. S. C. Gurney, R. M. Nisbet, The regulation of inhomogeneous populations, Journal of Theoretical Biology 52 (2) (1975) 441–457.  doi:10.1016/0022-5193(75)90011-9.
  * [23] D. J. Warne, R. E. Baker, M. J. Simpson, Using experimental data and information criteria to guide model selection for reaction–diffusion problems in mathematical biology, Bulletin of Mathematical Biology 81 (6) (2019) 1760–1804.  doi:10.1007/s11538-019-00589-x.
  * [24] F. Sánchez-Garduño, P. K. Maini, Existence and uniqueness of a sharp travelling wave in degenerate non-linear diffusion Fisher-KPP equations, Journal of Mathematical Biology 33 (2) (1994) 163–192.  doi:10.1007/BF00160178.
  * [25] J. A. Sherratt, B. P. Marchant, Nonsharp travelling wave fronts in the Fisher equation with degenerate nonlinear diffusion, Applied Mathematics Letters 9 (5) (1996) 33–38.  doi:10.1016/0893-9659(96)00069-9.
  * [26] S. L. Mitchell, M. Vynnycky, Finite-difference methods with increased accuracy and correct initialization for one-dimensional Stefan problems, Applied Mathematics and Computation 215 (4) (2009) 1609–1621.  doi:10.1016/j.amc.2009.07.054.
  * [27] S. L. Mitchell, T. G. Myers, Improving the accuracy of heat balance integral methods applied to thermal problems with time dependent boundary conditions, International Journal of Heat and Mass Transfer 53 (17) (2010) 3540–3551.  doi:10.1016/j.ijheatmasstransfer.2010.04.015.
  * [28] M. P. Dalwadi, S. L. Waters, H. M. Byrne, I. J. Hewitt, A mathematical framework for developing freezing protocols in the cryopreservation of cells, SIAM Journal on Applied Mathematics 80 (2) (2020) 657–689.  doi:10.1137/19M1275875.
  * [29] F. Brosa Planella, C. P. Please, R. A. Van Gorder, Extended Stefan problem for solidification of binary alloys in a finite planar domain, SIAM Journal on Applied Mathematics 79 (3) (2019) 876–913.  doi:10.1137/18M118699X.
  * [30] F. Brosa Planella, C. P. Please, R. A. Gorder, Extended Stefan problem for the solidification of binary alloys in a sphere, European Journal of Applied Mathematics 32 (2) (2021) 242–279.  doi:10.1017/S095679252000011X.
  * [31] J. Crank, Free and Moving Boundary Problems, Oxford University Press, 1987.
  * [32] R. Shuttleworth, D. Trucu, Multiscale modelling of fibres dynamics and cell adhesion within moving boundary cancer invasion, Bulletin of Mathematical Biology 81 (7) (2019) 2176–2219.  doi:10.1007/s11538-019-00598-w.
  * [33] L. S. Kimpton, J. P. Whiteley, S. L. Waters, J. R. King, J. M. Oliver, Multiple travelling-wave solutions in a minimal model for cell motility, Mathematical Medicine and Biology 30 (3) (2013) 241–272.  doi:10.1093/imammb/dqs023.
  * [34] N. T. Fadai, M. J. Simpson, New travelling wave solutions of the Porous–Fisher model with a moving boundary, Journal of Physics A: Mathematical and Theoretical 53 (9) (2020) 095601.  doi:10.1088/1751-8121/ab6d3c.
  * [35] E. A. Gaffney, P. K. Maini, C. D. McCaig, M. Zhao, J. V. Forrester, Modelling corneal epithelial wound closure in the presence of physiological electric fields via a moving boundary formalism, Mathematical Medicine and Biology: A Journal of the IMA 16 (4) (1999) 369–393.  doi:10.1093/imammb/16.4.369.
  * [36] J. P. Ward, J. R. King, Mathematical modelling of avascular-tumour growth, Mathematical Medicine and Biology: A Journal of the IMA 14 (1) (1997) 39–69.  doi:10.1093/imammb/14.1.39.
  * [37] Y. Du, Z. Lin, Spreading–vanishing dichotomy in the diffusive logistic model with a free boundary, SIAM Journal on Mathematical Analysis 42 (1) (2010) 377–405.  doi:10.1137/090771089.
  * [38] Y. Du, Z. Guo, Spreading–vanishing dichotomy in a diffusive logistic model with a free boundary, II, Journal of Differential Equations 250 (12) (2011) 4336–4366.  doi:10.1016/j.jde.2011.02.011.
  * [39] Y. Du, Z. Guo, The Stefan problem for the Fisher–KPP equation, Journal of Differential Equations 253 (3) (2012) 996–1035.  doi:10.1016/j.jde.2012.04.014.
  * [40] S. W. McCue, M. El-Hachem, M. J. Simpson, Exact sharp-fronted travelling wave solutions of the Fisher–KPP equation, Applied Mathematics Letters 114 (2021) 106918.  doi:10.1016/j.aml.2020.106918.
  * [41] M. El-Hachem, S. W. McCue, J. Wang, Y. Du, M. J. Simpson, Revisiting the Fisher–Kolmogorov–Petrovsky–Piskunov equation to interpret the spreading–extinction dichotomy, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 475 (2229) (2019) 20190378.  doi:10.1098/rspa.2019.0378.
  * [42] M. J. Simpson, Critical length for the spreading–vanishing dichotomy in higher dimensions, The ANZIAM Journal 62 (1) (2020) 3–17.  doi:10.1017/S1446181120000103.
  * [43] S. Osher, J. A. Sethian, Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations, Journal of Computational Physics 79 (1) (1988) 12–49.  doi:10.1016/0021-9991(88)90002-2.
  * [44] J. A. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, 2nd Edition, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, 32 Avenue of the Americas, New York, NY 10013-2473, USA, 1999.
  * [45] S. Osher, R. P. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Applied Mathematical Sciences, Springer-Verlag, New York, 2003.
  * [46] T. D. Aslam, A partial differential equation approach to multidimensional extrapolation, Journal of Computational Physics 193 (1) (2004) 349–355.  doi:10.1016/j.jcp.2003.08.001.
  * [47] C. Tsitouras, Runge–Kutta pairs of order 5(4) satisfying only the first column simplifying assumption, Computers & Mathematics with Applications 62 (2) (2011) 770–775.  doi:10.1016/j.camwa.2011.06.002.
  * [48] C. V. Rackauckas, Q. Nie, DifferentialEquations.jl – A performant and feature-rich ecosystem for solving differential equations in Julia, Journal of Open Research Software 5 (1) (2017) 15.  doi:10.5334/jors.151.
  * [49] C. T. Lin, E. Tadmor, High-resolution nonoscillatory central schemes for Hamilton–Jacobi equations, SIAM Journal on Scientific Computing 21 (6) (2000) 2163–2186.  doi:10.1137/S1064827598344856.
  * [50] G. S. Jiang, D. Levy, C. T. Lin, S. Osher, E. Tadmor, High-resolution nonoscillatory central schemes with nonstaggered grids for hyperbolic conservation laws, SIAM Journal on Numerical Analysis 35 (6) (1998) 2147–2168.  doi:10.1137/S0036142997317560.
  * [51] M. A. Lewis, P. Kareiva, Allee dynamics and the spread of invading organisms, Theoretical Population Biology 43 (2) (1993) 141–158.  doi:10.1006/tpbi.1993.1007.
  * [52] Y. Li, S. T. Johnston, P. R. Buenzli, P. van Heijster, M. J. Simpson, Extinction of bistable populations is affected by the shape of their initial spatial distribution, Bulletin of Mathematical Biology 84 (1) (2021) 21.  doi:10.1007/s11538-021-00974-5.
  * [53] M. Abramowitz, I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, tenth Edition, Applied Mathematics Series, United States Department of Commerce, National Bureau of Standards; Dover Publications, Washington, D. C., 1964.


††thanks<EMAIL_ADDRESS>
# Dropleton-Soliton Crossover mediated via Trap Modulation

Argha Debnath  Ayan Khan  Department of Physics, School of Engineering and
Applied Sciences, Bennett University, Greater Noida, Uttar Pradesh, India
Saurabh Basu  Department of Physics, Indian Institute of Technology -
Guwahati, Guwahati, Assam, India

###### Abstract

We report a droplet to a soliton crossover by tuning the external confinement
potential in a dilute Bose-Eienstein condensate by numerically solving the
modified Gross-Pitaevskii equation. The testimony of such a crossover is
presented via studying the fractional density of the condensate which smoothly
migrates from being a flat-head curve at weak confinement to a bright soliton
at strong confinement. Such a transition occurs across a region of the
potential whose strength varies over an order of magnitude and thus should be
fit to be termed as a crossover. We supplement our studies via exploring the
size of the bound pairs and the ramifications of the particle density therein.
Eventually, all of these aid us in arriving at a phase diagram in a space
defined by the trap strength and the particle number that shows the formation
of two phases consisting of droplets and solitons, along with a regime of
coexistence of these two.

##  I Introduction

The three states of matter, namely, solid, liquid and vapor are integral part
of our daily experience. However, under some extreme conditions other states
of matter can also emerge, such as plasmas at high fields and temperatures,
and superfluids at extremely low temperatures. These different phases are
manifestation of distinct interactions between their constituent particles.
For example, in a regular liquid, we observe the competition between van der
Waals inter particle attraction and the short-range repulsion due to the Pauli
exclusion principle. The transition from a liquid to a solid phase is also
observed by lowering the temperature. Interestingly, quantum liquids may not
exhibit similar properties. For example, helium remains a liquid close to
about absolute zero temperature (under atmospheric pressures) due to the large
zero-point motion of the constituent atoms. At this low temperature quantum
mechanical effects dominate, leading to the emergence of superfluidity.
Nevertheless, they can still poses some characteristics similar to classical
liquid such as droplet formation Böttcher  _et al._ (2020).

The experimental observation of atomic Bose-Einstein condensate (BEC) Anderson
_et al._ (1995); Bradley  _et al._ (1995); Davis  _et al._ (1995) had paved
the way for the emergence of new experimental platforms Dalfovo  _et al._
(1999); Giorgini  _et al._ (2008); Bloch  _et al._ (2008). This has allowed us
to study quantum phenomena and exotic states of matter with precise
controllability and tunability Chin  _et al._ (2010). As a result, lately
liquid-like state in a two component BEC Cabrera  _et al._ (2018) has been
reported. This unique development influences us to investigate the exotic
phase more closely as the existing understanding of liquid formation has
considerable influence from the theory of Van der Waals. Yet, these newly
emerged droplets in dilute ultra-cold atomic gases do not explicitly follow
the familiar notion of classical liquids Ferrier-Barbut (2019). These are
truly quantum liquids where quantum fluctuation plays an important role
Ferrier-Barbut  _et al._ (2016); Kadau  _et al._ (2016). These droplets are
small clusters of atoms self-bound by the interplay of attractive and
repulsive forces emerging from the mean-field (MF) and the beyond mean-field
(BMF) interactions Cappellaro  _et al._ (2018). The BMF contribution is
nothing but the Lee-Huang-Yang’s (LHY) correction Lee  _et al._ (1957) of the
Gross-Pitaevskii (GP) equation Gross (1961); Pitaevskii (1961). In a two
component BEC, the effective MF interaction can be significantly suppressed
due to the competitions between the intra-species and inter-species
interactions. This may result in a situation when the effective MF interaction
strength is in the order of LHY correction and competing in nature. This leads
to an ideal scenario for droplet formation. The experimental verification was
based on the theoretical proposition where it was claimed that, if the square
of the inter-species coupling is greater than the product of the intra-species
coupling then the collapse of the two-component BEC can be averted and a
dilute liquid-like droplet state emerges Petrov (2015).

However it must be noted that, the first experimental observation of droplet
was on dipolar condensate of ${}^{164}Dy$ Ferrier-Barbut  _et al._ (2016);
Kadau  _et al._ (2016) with subsequent theoretical developments coming to
light Wächtler and Santos (2016). Later the droplets were observed for a
binary BEC comprised of two hyper-fine states in ${}^{39}K$ Cabrera  _et al._
(2018); Ferioli  _et al._ (2019). Later, soliton to dropleton transition was
also reported Cheiney  _et al._ (2018). Quantum droplets were also observed in
a heteronuclear mixture of bosonic gases D’Errico  _et al._ (2019). In a
recent self-consistent theoretical formulation LHY correction is incorporated
in the GP equation via quantum fluctuation Salasnich (2018). We have also
noted theoretical assertions on the collective modes in a droplet-soliton
crossover Cappellaro  _et al._ (2018), existence of vortices in droplets Li
_et al._ (2018), dynamics of purely one-dimensional droplets Astrakharchik and
Malomed (2018) and its collective excitations Tylutki  _et al._ (2020) etc. At
this juncture, we also recall a numerical investigation of quantum liquids for
the dipolar BEC in quasi one-dimensional (Q1D) geometry Edmonds  _et al._
(2020). Apart from that, we have also seen considerable interest in analyzing
the origin of the droplets in lower dimensions Shamriz  _et al._ (2020). A
contemporary study also notes possible connection between the droplets and
modulational instability in an one dimensional system Mithun  _et al._ (2020).

Very recently, we have proposed an exact analytical solution for the droplet-
like state in a Q1D quasi-homogeneous cigar shaped BEC where we have assumed
the longitudinal confinement is sufficiently weak so that it can be neglected
Debnath and Khan (2021). In that article we have also commented on the soliton
to droplet transition. However, from an experimental perspective, it is more
prudent to consider the cigar shaped harmonic confinement in Q1D as depicted
in Fig. 1. Tuning the confinement potential should yield additional features
that may be worth exploring Dey  _et al._ (2012). So, we start with an
expectation of finding some novel consequences by tuning the trapping
potential. Our systematic analysis reveals that, the trap frequency indeed
plays a significant role in understanding the droplet-soliton transition.

In this article, we consider a binary BEC of two hyper-fine states while both
components occupy the same spatial mode Cabrera  _et al._ (2018). As a result,
the two-component nonlinear Schrödinger equation (NLSE) is reduced to an
effective one component equation. Next, we employ a dimensional reduction
reduction scheme following the prescription of Ref.Atre  _et al._ (2006). This
allows us to reduce the three dimensional system to a quasi one dimension
(Q1D) without loss of much generality. Physically, a dimensional reduction
makes sense, when the transverse components of the anisotropic harmonic
confinement is much stronger compared to the longitudinal trapping frequency.
The Q1D system now consists of a cubic and a quartic nonlineairty, contrary to
the prevailing concept of NLSE where we generally deal with only a cubic
nonlinearity. The cubic term defines the effective mean-field (EMF) two-body
interaction and the quartic term is the signature of beyond mean-field
contribution. We have recently carried out a systematic analysis for cubic-
quartic nonlinear Schrödinger equation (CQNLSE) Debnath and Khan (2021) in
quasi one dimension and reported exact analytical solutions. Further, we have
analyzed and portrayed the existence of cnoidal functions Abramowitz and
Stegun (1972) in the CQNLSE when subjected to cnoidal confinements Debnath and
Khan (2020). Lately, we have investigated on the existence of periodic modes
in CQNLSE and commented on the dimensional crossover from Q1D to an one
dimensional system Debnath  _et al._ (2021). In this article, our main
objective is to focus on the role played by the trap potential in droplet
formation and possible transition from a droplet to a solitonic state via
tuning the confinement potential and studied numerically by solving the
inhomogeneous CQNLSE.

In this paper, we communicate our results in the following order, in Sec.II,
we review the necessary theory which maps the two component BEC to one
component extended Gross-Pitaevskii (GP) equation. Thereafter we explicate the
formalism to reduce the three dimensional system to Q1D geometry. We briefly
recap the scheme to derive the analytical solution and compare them with our
current numerical code in Sec.III. In Sec.IV, first we investigate the effect
of the longitudinal trapping potential for fixed number of particles at
different values of the BMF interaction strength. Later we extend the analysis
toward the role of particle number and draw phase diagram. This allows us to
comment on the nature of droplet to soliton transition. We draw our conclusion
in Sec.V.

##  II Theoretical Model

Figure 1: (Color online) Pictorial illustration of binary BEC in a Q1D
confinement. The two different hyper-fine states are described by the bluish
and reddish spheres. The intra-species and inter-species interactions are
considered to be repulsive and attractive in nature.

Here, we consider a mixture of two hyper-fine states of ${}^{39}K$ in a Q1D
geometry following the prescription of Ref.Cabrera  _et al._ (2018). Here, we
like to note that there exists considerable research on purely one dimensional
(1D) systems Astrakharchik and Malomed (2018) where the energy of a weakly
interacting Bose gas can be predicted using Bogoliubov theory Lieb and Liniger
(1963); Popov (1972); Petrov and Astrakharchik (2016). However, it is also
true that condensate formation is not possible in 1D. Hence, it is common to
study the condensate formation in a Q1D geometry where the Bose gas is allowed
to expand in an optical waveguide. This dimensional engineering has enabled us
to observe exotic structures BEC like the bright soliton trains Khaykovich
_et al._ (2002); Strecker  _et al._ (2002). Here, by stating a Q1D system, we
presume that $\sqrt{na^{3}}<1$ although $\frac{1}{|\sqrt{n_{1D}a_{1D}}|}>1$
Ilg  _et al._ (2018). Here $n$ and $a$ denotes the particle density and
$s$-wave scattering length while $n_{1D}$ and $a_{1D}$ are the density and
scattering length in 1D respectively. Further we note that
$a_{1D}=2\hbar^{2}a/a_{\perp}^{2}m$, where $a_{\perp}$ is the characteristic
length scale of the trap and $m$ is the mass of the particle Edler  _et al._
(2017). Nevertheless, our objective is to remain more towards the right side
of the dimensional crossover, whereas a strict 1D system is more towards the
left of the crossover with $\frac{1}{|\sqrt{n_{1D}a_{1D}}|}<1$ Ilg  _et al._
(2018). From mathematical perspective, the fundamental difference between Q1D
and a 1D systems lies in the fact that equation of motion of the prior leads
to a CQNLSE where as the later can be described via quadratic-cubic NLSE
(QCNLSE) Petrov and Astrakharchik (2016).

A schematic description of the physical system can be visualized in Fig. 1
where the bluish and the reddish spheres depicts the two different hyper-fine
states of atoms distributed in an effectively cigar shaped or Q1D trap. The
intra-species interactions ($a_{11}$ and $a_{22}$) are repulsive in nature and
the inter-species interaction ($a_{12}$ or $a_{21}$) are attractive. The
experimentally observed droplets are described as small clusters of atoms
self-bound through the competition between the attractive and the repulsive
forces. It can be shown that when $-\sqrt{a_{11}a_{22}}>a_{12}$ then the
condensate is on the verge of collapse. Now, one can model the effective mean-
field interaction strength as $\delta a$ where $\delta
a\propto(a_{12}+\sqrt{a_{11}a_{22}})$. If $\delta a\lesssim 0$, then the
beyond mean-field (BMF) contribution can not be neglected. We can denote the
BMF interaction strength as $\delta
a^{\prime}\propto(\sqrt{a_{11}a_{22}})^{5/2}$ Cheiney  _et al._ (2018). In
this situation one can describe the system with an effective single component
GP equation by assuming that the two components occupy the same spatial mode.
Thus the resulting one component equation of motion, with cubic and quartic
nonlinearities, can be defined as Cabrera  _et al._ (2018); Debnath and Khan
(2021),

$\displaystyle i\hbar\frac{\partial\Psi}{\partial t}$ $\displaystyle=$
$\displaystyle\left[\left(-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{trap}\right)+U|\Psi|^{2}+U^{\prime}|\Psi|^{3}\right]\Psi,$

where, $U=\frac{4\pi\hbar}{m}\delta a$,
$U^{\prime}=\frac{256\sqrt{\pi}\hbar^{2}\delta a^{\prime}}{15m}$, $m$ being
the mass of the atoms and $V_{trap}$ describes the external potential. It is
possible to describe $V_{trap}$ as combination of the transverse and
longitudinal components of the harmonic confinement. The transverse trap can
be defined as, $V_{T}(y,z)=\frac{1}{2}m\omega_{\perp}^{2}(y^{2}+z^{2})$ and
the longitudinal component is noted as,
$V_{L}(x)=\frac{1}{2}m\omega_{0}^{2}x^{2}$. Here, $\omega_{\perp}$ is the
transverse trap frequency and $\omega_{0}$ signifies the longitudinal trap
frequency. In cigar-shaped BEC or in Q1D geometry transverse trap frequency is
typically $10$ times larger than the longitudinal frequency. Since, the
characteristic length scale happens to be
$a_{\perp}=\sqrt{\frac{\hbar}{m\omega_{\perp}}}$, therefore in a Q1D geometry,
one can assume $a_{x}/a_{\perp}\sim\sqrt{10}$. Here, $a_{x}$ is given by
$a_{x}=\sqrt{\frac{\hbar}{m\omega_{0}}}$ Salasnich  _et al._ (2002). These
assertions ensure that the interaction energy of the atoms remains much lesser
compared to the kinetic energy in the transverse direction.

Consequently, it is possible to reduce Eq.(LABEL:3dbgp) to an effective one-
dimensional equation. The dimensional reduction is performed by employing the
following ansatz,

$\displaystyle\Psi(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\pi
a_{B}}a_{\perp}}\psi\left(\frac{x}{a_{\perp}},\omega_{\perp}t\right)e^{\left(-i\omega_{\perp}t-\frac{y^{2}+z^{2}}{2a_{\perp}^{2}}\right)},$

where, $a_{B}$ is the Bohr radius.

Applying the ansatz from Eq.(LABEL:ansatz1) in Eq.(LABEL:3dbgp) we obtain the
extended GP equation as noted below,

$\displaystyle i\frac{\partial\psi(x,t)}{\partial t}=$
$\displaystyle\left[-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}+\frac{1}{2}Kx^{2}+\tilde{g}|\psi(x,t)|^{2}\right.$
$\displaystyle\qquad\left.+\frac{}{}\tilde{g}^{\prime}|\psi(x,t)|^{3}\right]\psi(x,t),$
(3)

where, $\tilde{g}=2\delta a/a_{B}$,
$\tilde{g}^{\prime}=(64\sqrt{2}/15\pi)\delta
a^{\prime}/(a_{B}^{3/2}a_{\perp})$ and $K=\omega^{2}_{0}/\omega^{2}_{\perp}$.
Here $x$ and $t$ are dimensionless as $x\equiv x/a_{\perp}$ and
$t\equiv\omega_{\perp}t$. From here onward, we shall follow this dimensionless
notations of $x$ and $t$.

In this article, our main focus is to explicate the role of $K$ in dropleton-
soliton transition. For this purpose, we solve Eq.(3) numerically for
different values of $K$. The numerical scheme involved a split-step Crank-
Nicolson (CN) method with imaginary time propagation Muruganandam and Adhikari
(2009). However, we prefer to briefly revisit our analytical solution for
better readability.

##  III Solution

To obtain an analytical solution, we first assume that the wave function can
be composed of a slow moving amplitude part and a fast moving phase part such
that ,
$\psi(x,t)=\rho(x,t)\exp{\left[i\left(\chi(x,t)+\mu_{0}t\right)\right]}$ where
$\rho(x,t)$ leads to the amplitude contribution and $\chi(x,t)$ is the non-
trivial phase, $\mu_{0}$ being the chemical potential. After some algebra, we
arrive at the stable analytical solution in the center of mass frame as
Debnath and Khan (2021),

$\displaystyle\rho(\zeta)$ $\displaystyle=$
$\displaystyle\frac{1+12\mu_{g}}{1+\sqrt{12\mu_{g}}\cosh{\left(\sqrt{g}\zeta\right)}}\,\textrm{for}\,|g|=3g^{\prime}.$
(4)

Here, $\zeta=x-ut$ described the comoving frame with $u$ being the velocity of
the condensate. In the analytical derivation we have also assumed that
$g=-2\tilde{g}$ and $g^{\prime}=2\tilde{g}^{\prime}$ implying a two-body
effective mean-field interaction is attractive and LHY contribution is
repulsive and $\mu_{g}=\mu_{0}/g$.

The relationship between the normalization, $N$ and chemical potential,
$\mu_{0}$ can be noted as,

$\displaystyle N$ $\displaystyle=$
$\displaystyle\frac{(1+\mu_{I})^{2}}{\sqrt{g}(1-\mu_{I})}\left[\ln\left[1-\frac{2}{\mu_{I}}\left(\sqrt{1-\mu_{I}}+1\right)\right]-2\right]$
(5)

Here, for the convenience of calculation we have denoted, $12\mu_{g}=\mu_{I}$.
Physically, $N$ denotes particle number associated with the formation of
localized wave and is scaled by $N_{0}$ where $N_{0}$ defines the number of
particle obtained from the constant background density solution such that,
$N_{0}=2a_{B}\left(\frac{15\pi a_{\perp}}{64}\right)^{2}\left(\frac{\delta
a}{\delta a^{\prime}}\right)^{2}$.

We test our numerical scheme, by comparing with the analytical result as
depicted in Fig. 2. In the figure, the solid blue line is our analytical
solution from Eq.4 and the red solid squares are the numerically obtained
solution of Eq.(3) for $K=0$. Fig. 2 corresponds to the droplet like state.
The computation appears in good agreement with the analytical result and thus
the scenario encourages us to proceed towards our actual objective of solving
Eq.(3) with $K\neq 0$.

Figure 2: (Color online) The figure depicts the comparison between the
analytical solution using Eq.(4) and numerical solution of Eq.(3) (with
$K=0$). The analytical solution is plotted with solid blue line whereas the
red squares represents numerically obtained solution. The density is
normalized by $n_{0}$ where $n_{0}$ is $n(x)|_{x=0}$.

##  IV Droplet to Soliton Transition

The satisfactory comparison of our numerical result with the analytical
solution motivates us to solve Eq.(3) with $K\neq 0$. We employ the imaginary
time propagation method over small time steps ($5\times 10^{-6}$) for total
$20000$ steps and solve the inhomogeneous CQNLSE for several trap frequencies
and different BMF interaction strengths. We introduce the interaction strength
in terms of $\kappa=g^{\prime}/|g|$ and we perform numerical calculations for
three arbitrary values of $\kappa$, such that $\kappa=0.1,0.5$, and $0.9$.
Here, it must be noted that in Ref.Debnath and Khan (2021), the analytical
solution was obtained for $\kappa=1/3$. In Fig. 3 we depict the density
distribution for $K=0.00001$ (purple dashed-dotted line), $K=0.001$ (green
dashed line) and $K=0.1$ (blue solid line) respectively. We can clearly see
that the density profile indicate towards a transition from a droplet to a
soliton as the trap frequency is increased. At very low trap frequency, that
is $K=0.00001$, the density plateau is distinctly visible indicating the
droplet state. However, with gradual increase in trap frequency (say at
$K=0.001$) we observe a bell like structure. A tighter confinement leads to
the generation of bright soliton as described by the solid blue line. The
figure was obtained for $\kappa=0.1$ however, we note that the profile is not
significantly altered for larger values of $\kappa$ till $\kappa=0.9$.

Figure 3: (Color online) The density profile of the numerical solution for
different $K$ value is depicted here. The purple dashed-dotted line described
the density for $K=0.00001$, the green dashed line describes $K=0.001$ and the
blue solid line is prepared for $K=0.1$. Here, $\kappa=g^{\prime}/|g|=0.1$ and
$g$ is set at $1$.

To excavate more on droplet-soliton transition, the chemical potential can be
regarded as an important physical parameter. Hence, we assume,
$\psi(x,t)=\phi(x)e^{-i\mu t}$ and apply in Eq.(3). This results,

$\displaystyle\mu\phi(x)=\left[-\frac{1}{2}\frac{d^{2}}{dx^{2}}+\frac{1}{2}Kx^{2}+g\phi^{2}(x)+g^{\prime}\phi^{3}(x)\right]\phi(x).$

Applying the normalization condition $\int_{-\infty}^{\infty}\phi^{2}(x)dx=1$
we obtain,

$\displaystyle\mu=\int_{-\infty}^{\infty}\left[\frac{1}{2}\left(\frac{d\psi}{dx}\right)^{2}+\phi^{2}(x)\left(\frac{1}{2}Kx^{2}+g\phi^{2}(x)+g^{\prime}\phi^{3}(x)\right)\right].$

We report the variation of the chemical potential with trap frequency for
different beyond mean-field interaction strength in Fig 4. The purple dashed-
dotted line, green dashed-double-dotted line and blue solid line depicts
$\kappa=0.1,0.5$ and $0.9$ respectively. For weaker trap frequency the system
moves more in the droplet region which is signified by the negative value of
the chemical potential. A very weak trap frequency also leads the system to a
nearly homogeneous setup where we obtain the chemical potentials for different
$\kappa$ is converging to $\sim-0.02$. Interestingly, similar result we have
already reported for $\kappa=0.33$ Debnath and Khan (2021). As noted in Fig 3,
the stronger confinement leads to the solitonic state with $\mu>0$ (right side
of the grey shaded region of Fig 4). However, we observe a non-monotonic
behaviour in the region $K=0.001$ and $0.1$ in the figure (the grey shaded
region). A comparison with Fig 3 encourages us to consider this frequency
window as the transition region as we observe a bell like density profile for
$K=0.001$.

Figure 4: (Color online) The variation of chemical potential at different
$\kappa$ is depicted here. The solid blue line corresponds to $\kappa=0.9$,
green dashed-double-dotted and purple dashed-dotted line describes variation
of chemical potential with trap frequency for $\kappa=0.5$ and $0.1$
respectively. The black dashed line indicates $\mu=0$. The shaded area
describes the crossover regime.

To investigate the transition we then calculate the first derivative of
chemical potential with respect to the trap frequency. The obtained result is
depicted in Fig 5. The main question here, is whether the transition is a
phase transition (quantum) or a crossover. We know that, for quantum phase
transition (QPT) we observe an abrupt change in the ground state of the many
body system when a controlling parameter ($\lambda$) of the Hamiltonian
crosses a critical value Khan and Pieri (2009); Dey  _et al._ (2011). Here, we
recognize $\lambda=K$, or the longitudinal trapping frequency. However, in Fig
4 and Fig 5 we do not observe any abrupt change in the chemical potential.
Thus, to our assessment the transition appears to be gradual and hence can be
termed as a crossover. Based on these evidences we can also quantify the
crossover region as $0.0001\lesssim K\lesssim 0.1$ (described in the grey
shaded region in the figures).

Figure 5: (Color online) The first derivative of chemical potential w.r.t the
trap frequency at different $\kappa$ is described. The color coding and the
line types for three different $\kappa$ follows the same definition as of Fig
4. The shaded region describes the transition.

Further, we calculate the total energy which can be defined as,

$\displaystyle E$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\left[\frac{1}{2}\left(\frac{d\psi}{dx}\right)^{2}+\phi^{2}(x)\left(\frac{1}{2}Kx^{2}+\frac{1}{2}g\phi^{2}(x)\right.\right.$
(8)
$\displaystyle\left.\left.+\frac{2}{5}g^{\prime}\phi^{3}(x)\right)\right].$

The variation of energy corresponding to the trap frequency modulation is
captured in Fig 6. The figure reveals that the energy is negative for weak
confinement suggesting the accumulation of droplet like bound pairs. The
energy crosses the zero line in the vicinity of $K=0.001$ indicating the
breakdown of droplet-like bound pairs. However, the energy flattening in the
shaded area might be suggestive of an equilibrium of dropleton-soliton
mixture. The bell shaped density profile can be noted as a signature of mixed
type of states. Further tightening of the longitudinal frequency leads to the
total destruction of the droplets.

Figure 6: (Color online) Change in energy with variation in the trap frequency
is plotted here for different BMF interaction strengths. Again, the solid blue
line, green dashed-double-dotted and purple dashed-dotted line denotes the
energy variation as a function of trap frequency for $\kappa=0.9,0.5$ and
$0.1$ respectively. The black dashed line indicates $E=0$. The grey area
corresponds to the region of transition from droplet to soliton.

Next, we calculate the size variation in the droplet-soliton transition as
depicted in Fig 7 by calculating the root-mean-square (rms) ($\sqrt{\langle
x^{2}\rangle}$) size of the condensate. The rms size is depicted in units of
the rms size ($\sqrt{\langle x^{2}_{0}\rangle}$) for $K=0$. The figure clearly
reflects on the fact that the variation of the BMF interaction does not have
much contribution on the droplet size of the condensate as we report the
variation for $\kappa=0.1$ (solid red line), $0.5$ (solid green squares) and
$0.9$ (solid blue circles) respectively. We can also conclude that the weak
confinement allows the larger pair size (with asymptotic limit extending to
$1$) suggesting the droplet formation and gradually the pair size smoothly
falls off as we increase the longitudinal frequency. We also observe that the
pair size drops much faster in the shaded region, which we have defined as the
crossover region. The rms curve flattens thereafter and asymptotically reaches
$1/10$ of the droplet size for a tightly confined one dimensional geometry.
The effect of the variation in $\kappa$ is apparently negligible for the
normalized pair size, as noted in the figure.

Figure 7: (Color online) The figure describes the variation of the rms size of
the droplets with the modulation of the trapping potential. The rms size is
defined in units of $\sqrt{\langle x_{0}^{2}\rangle}$ where $\langle
x_{0}^{2}\rangle$ is the rms size for $K=0$. The solid red line, the green
solid squares and the blue solid circles denotes $\kappa=0.1$, $0.5$ and $0.9$
respectively. The shaded region describes the crossover area.

It is also instructive now to investigate the role of particle number of the
condensate. We calculate the energy and chemical potential for different
particle numbers (from 25 to 100). We draw a phase diagram based on our
finding which is depicted in Fig. 8. The blue region corresponds to the
negative energy region, thereby suggesting of the existence of the bound pairs
and droplet formation. The red region is describes the unbound solitonic
region, and the whitish region corresponds to region of crossover where the
droplets are breaking down into solitons. The figure also indicates that, for
large particle number, even with very weak $K$ there is less possibility of
droplet formation. As already been discussed in Ref.Debnath and Khan (2021),
there exists a critical density beyond which droplet does not form.
Corroborating Fig 8 with Fig 4 and Fig 6 we can conclude that, in the region
where the chemical potential is negative for relatively higher value of $K$,
the energy of the system reaches a plateau and behaves like a latent energy
where a mixture of droplet and soliton do coexist.

Figure 8: (Color online) The figure describes the energy as a function of
longitudinal trap strength and particle number. A clear phase separation can
be noted by following the whitish region.

##  V Conclusion

The recent experimental and theoretical developments related to the newly
observed liquid like state in ultra-cold atomic BEC motivates us to study the
role of harmonic confinement in the investigation of droplets. Already we have
shown the existence of analytical solution pertaining to a droplet-like state
in a quasi one dimensional quasi homogeneous system Debnath and Khan (2021).
Here, we start from where we left in the previous analytical description and
explore the droplet state numerically in inhomogeneous quasi one dimensional
system. We find that the longitudinal trap frequency can be an appropriate
controlling parameter in experiments to access a crossover from a droplet to a
solitonic mode. We observe, that at very weak trap frequencies, the droplet
states do exists just like as they do in the homogeneous case. However, with
systematic increase in the trapping potential we move from flat top droplet
solution to bell shaped solution and later to the emergence of solitons. In
this regard, we especially emphasize on the numerical calculation of the
chemical potential and energy. For very weak trapping potential, the chemical
potential converges to the analytical value reported earlier Debnath and Khan
(2021). At moderate values of the trapping frequency, the behavior of $\mu$ is
nonlinear (the shaded regions in figure) and we observe a plateau in energy.
Here, the droplet and soliton appear to be coexisting, corresponding to a bell
type structure in the density. At relatively tighter confinement, the bound
droplets completely disappear and the solitonic mode prevails. We have
repeated our investigation for different BMF interactions strengths
($\kappa=0.1$, $0.5$ and $0.9$) while noting a similar trend in the
transition. To characterize the transition, we calculate the derivative of
chemical potential however, we do not observe any significant change in the
derivative as a function of the change of the Hamiltonian parameter ($K$
here). We also observe a smooth transition of the pair size across a weak to a
strong confinement domain. This allows us to conclude that the droplet to
soliton transition denotes a crossover phenomena.

As we report a confinement driven transition from droplet to soliton, we
expect to open up new avenues in understanding the droplet state. We believe a
study of fidelity of the numerically obtained solutions can shed more light in
understanding the transition Khan and Pieri (2009). We also expect that our
prediction can soon be verified in experiments.

## Acknowledgement

Authors acknowledge insightful discussion with P. Sadhukhan. AK also thanks
Department of Science and Technology (DST), India for the support provided
through the project number CRG/2019/000108.

## References

  * Böttcher  _et al._ (2020) F. Böttcher, J.-N. Schmidt, J. Hertkorn, K. Ng, S. Graham, M. Guo, T. Langen, and T. Pfau, Reports on Progress in Physics (2020).
  * Anderson  _et al._ (1995) M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science 269, 198 (1995).
  * Bradley  _et al._ (1995) C. C. Bradley, C. Sackett, J. Tollett, and R. G. Hulet, Physical Review Letters 75, 1687 (1995).
  * Davis  _et al._ (1995) K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. Rev. Lett. 75, 3969 (1995).
  * Dalfovo  _et al._ (1999) F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999).
  * Giorgini  _et al._ (2008) S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 80, 1215 (2008).
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * Cabrera  _et al._ (2018) C. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, Science 359, 301 (2018).
  * Ferrier-Barbut (2019) I. Ferrier-Barbut, Physics Today 72, 46 (2019).
  * Ferrier-Barbut  _et al._ (2016) I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel, and T. Pfau, Phys. Rev. Lett. 116, 215301 (2016).
  * Kadau  _et al._ (2016) H. Kadau, M. Schmitt, M. Wenzel, C. Wink, T. Maier, I. Ferrier-Barbut, and T. Pfau, Nature 530, 194 (2016).
  * Cappellaro  _et al._ (2018) A. Cappellaro, T. Macrì, and L. Salasnich, Phys. Rev. A 97, 053623 (2018).
  * Lee  _et al._ (1957) T. D. Lee, K. Huang, and C. N. Yang, Phys. Rev. 106, 1135 (1957).
  * Gross (1961) E. P. Gross, Il Nuovo Cimento (1955-1965) 20, 454 (1961).
  * Pitaevskii (1961) L. Pitaevskii, Sov. Phys. JETP 13, 451 (1961).
  * Petrov (2015) D. S. Petrov, Phys. Rev. Lett. 115, 155302 (2015).
  * Wächtler and Santos (2016) F. Wächtler and L. Santos, Phys. Rev. A 93, 061603 (2016).
  * Ferioli  _et al._ (2019) G. Ferioli, G. Semeghini, L. Masi, G. Giusti, G. Modugno, M. Inguscio, A. Gallemí, A. Recati, and M. Fattori, Phys. Rev. Lett. 122, 090401 (2019).
  * Cheiney  _et al._ (2018) P. Cheiney, C. R. Cabrera, J. Sanz, B. Naylor, L. Tanzi, and L. Tarruell, Phys. Rev. Lett. 120, 135301 (2018).
  * D’Errico  _et al._ (2019) C. D’Errico, A. Burchianti, M. Prevedelli, L. Salasnich, F. Ancilotto, M. Modugno, F. Minardi, and C. Fort, Phys. Rev. Research 1, 033155 (2019).
  * Salasnich (2018) L. Salasnich, Applied Sciences 8 (2018), 10.3390/app8101998.
  * Li  _et al._ (2018) Y. Li, Z. Chen, Z. Luo, C. Huang, H. Tan, W. Pang, and B. A. Malomed, Phys. Rev. A 98, 063602 (2018).
  * Astrakharchik and Malomed (2018) G. Astrakharchik and B. A. Malomed, Physical Review A 98, 013631 (2018).
  * Tylutki  _et al._ (2020) M. Tylutki, G. E. Astrakharchik, B. A. Malomed, and D. S. Petrov, Phys. Rev. A 101, 051601 (2020).
  * Edmonds  _et al._ (2020) M. Edmonds, T. Bland, and N. G. Parker, arXiv preprint arXiv:2002.07958 (2020).
  * Shamriz  _et al._ (2020) E. Shamriz, Z. Chen, and B. A. Malomed, Phys. Rev. A 101, 063628 (2020).
  * Mithun  _et al._ (2020) T. Mithun, A. Maluckov, K. Kasamatsu, B. A. Malomed, and A. Khare, Symmetry 12, 174 (2020).
  * Debnath and Khan (2021) A. Debnath and A. Khan, Ann. Phys. (Berlin) 533, 2000549 (2021).
  * Dey  _et al._ (2012) P. Dey, A. Khan, S. Basu, and B. Tanatar, in  _International Journal of Modern Physics: Conference Series_ , Vol. 11 (World Scientific, 2012) pp. 127–132.
  * Atre  _et al._ (2006) R. Atre, P. K. Panigrahi, and G. S. Agarwal, Physical Review E 73, 056611 (2006).
  * Abramowitz and Stegun (1972) M. Abramowitz and I. A. Stegun,  _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ , Vol. 55 (US Government printing office, 1972).
  * Debnath and Khan (2020) A. Debnath and A. Khan, The European Physical Journal D 74, 1 (2020).
  * Debnath  _et al._ (2021) A. Debnath, J. Tarun, and A. Khan, arXiv preprint arXiv:2107.11779 (2021).
  * Lieb and Liniger (1963) E. H. Lieb and W. Liniger, Physical Review 130, 1605 (1963).
  * Popov (1972) V. N. Popov, Theoretical and mathematical physics 11, 565 (1972).
  * Petrov and Astrakharchik (2016) D. S. Petrov and G. E. Astrakharchik, Phys. Rev. Lett. 117, 100401 (2016).
  * Khaykovich  _et al._ (2002) L. Khaykovich, F. Schreck, G. Ferrari, T. Bourdel, J. Cubizolles, L. D. Carr, Y. Castin, and C. Salomon, Science 296, 1290 (2002).
  * Strecker  _et al._ (2002) K. E. Strecker, G. B. Partridge, A. G. Truscott, and R. G. Hulet, Nature 417, 150 (2002).
  * Ilg  _et al._ (2018) T. Ilg, J. Kumlin, L. Santos, D. S. Petrov, and H. P. Büchler, Physical Review A 98, 051604 (2018).
  * Edler  _et al._ (2017) D. Edler, C. Mishra, F. Wächtler, R. Nath, S. Sinha, and L. Santos, Phys. Rev. Lett. 119, 050403 (2017).
  * Salasnich  _et al._ (2002) L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 65, 043614 (2002).
  * Muruganandam and Adhikari (2009) P. Muruganandam and S. K. Adhikari, Computer Physics Communications 180, 1888 (2009).
  * Khan and Pieri (2009) A. Khan and P. Pieri, Phys. Rev. A 80, 012303 (2009).
  * Dey  _et al._ (2011) P. Dey, D. Sarkar, A. Khan, and S. Basu, The European Physical Journal B 81, 95 (2011).