# Cell Division and Motility Enable Hexatic Order in Biological Tissues

Yiwen Tang  Department of Physics, Northeastern University, Boston, MA 02115
Siyuan Chen  Department of Physics, University of California, Santa Barbara,
Santa Barbara, CA 93106  Mark J. Bowick  Department of Physics, University of
California, Santa Barbara, Santa Barbara, CA 93106  Kavli Institute of
Theoretical Physics, University of California, Santa Barbara, Santa Barbara,
CA 93106  Dapeng Bi  Department of Physics, Northeastern University, Boston,
MA 02115

###### Abstract

Biological tissues transform between solid-like and liquid-like states in many
fundamental physiological events. Recent experimental observations further
suggest that in two-dimensional epithelial tissues these solid-liquid
transformations can happen via intermediate states akin to the intermediate
hexatic phases observed in equilibrium two-dimensional melting. The hexatic
phase is characterized by quasi-long-range (power-law) orientational order but
no translational order, thus endowing some structure to an otherwise
structureless fluid. While it has been shown that hexatic order in tissue
models can be induced by motility and thermal fluctuations, the role of cell
division and apoptosis (birth and death) has remained poorly understood,
despite its fundamental biological role. Here we study the effect of cell
division and apoptosis on global hexatic order within the framework of the
self-propelled Voronoi model of tissue. Although cell division naively
destroys order and active motility facilitates deformations, we show that
their combined action drives a liquid-hexatic-liquid transformation as the
motility increases. The hexatic phase is accessed by the delicate balance of
dislocation defect generation from cell division and the active binding of
disclination-antidisclination pairs from motility. We formulate a meanfield
model to elucidate this competition between cell division and motility and the
consequent development of hexatic order.

Organ surfaces are often covered with 2D confluent monolayers of epithelial or
endothelial cells, which provide functional separation from the surrounding
environment. During development these cells grow, divide and move, dynamically
reorganizing the entire tissue. Regulated by a complex set of chemical and
mechanical signaling pathwaysLecuit  _et al._ (2011); Janmey and Miller
(2011); Mammoto and Ingber (2010); Tenney and Discher (2009), tissue
frequently undergoes a transition from a structureless fluid-like state to a
state capable of supporting a variety of stresses, most notably elastic
stressesAtia  _et al._ (2018); Park  _et al._ (2015); Mongera  _et al._
(2018); Malinverno  _et al._ (2017); Angelini  _et al._ (2011); Garcia  _et
al._ (2015); Puliafito  _et al._ (2012). Such transformations have recently
been analyzed as a crossover from a liquid to an amorphous solidBi  _et al._
(2015, 2016). In two-dimensional (2D) systems in equilibrium, however, liquids
can develop rigidity via two consecutive transitions, the first corresponding
to the development of orientational order without translational order and the
second adding translational order to the existing orientational orderHalperin
and Nelson (1978); Nelson and Halperin (1979). The intermediate phase with
(quasi-long-range) orientational order but translational disorder is known as
the hexatic phase and has been shown to occur in a very wide variety of
physical systemsChou  _et al._ (1998); Murray (1992); Zahn  _et al._ (1999);
Zahn and Maret (2000); Han  _et al._ (2008); Seshadri and Westervelt (1991);
Knobler and Desai (1992); Huang (1992); Huang and Stoebe (1993); Podgornik
_et al._ (1996); Van Winkle  _et al._ (1997); Rill  _et al._ (1991); Bishop
_et al._ (1994); Murray  _et al._ (1990). The hexatic is a particular type of
structured fluid since it flows like a fluid but has orientational rigidity.

Previous theoretical and computational models of dense tissues have studied
the emergence of hexatic order, with focus on the effects of thermal
fluctuationsLi and Ciamarra (2018); Durand and Heu (2019); xue Guo  _et al._
(2023) and motilityPasupalak  _et al._ (2020); Loewe  _et al._ (2020);
Paoluzzi  _et al._ (2021); jian Li and quan Ai (2021). Modeling typically
studies the inverse process of disordering by melting from the crystalline
state. Realistic tissues, however, are very rarely crystalline with a few
exceptionsCohen  _et al._ (2020); Hilgenfeldt  _et al._ (2008). Cell division
and apoptosis almost always destroy the crystalline stateA. Matoz-Fernandez
_et al._ (2017) and yet there has been no direct observation of the hexatic
phase in in vitro biological tissues, including those undergoing a solid-
liquid transitionPark  _et al._ (2015); Malinverno  _et al._ (2017); Atia  _et
al._ (2018); Mitchel  _et al._ (2020).

Recent in vivo experiments on Drosophila embryos have uncovered hexatic order
during development with cell divisionKanesaki  _et al._ (2011); Wu  _et al._
(2021), along with the associated increase of orientational correlationsKaiser
_et al._ (2018). The mechanism behind the emergence of this orientational
order has remained unclear.

Here we analyze whether biological systems can exhibit this rather subtle
phase by studying numerically and analytically the self-propelled Voronoi
(SPV) model of cellular tissue including cell division and deathBi  _et al._
(2016). We compare a variety of structural properties with and without
division, including translational and orientational order parameters, order
field correlation functions in space, order field susceptibility, and
topological defect densities. We find that the interplay of division/apoptosis
and cell motility does indeed give rise to a hexatic regime. In the absence of
cell division, the model undergoes a crystal-hexatic and hexatic-liquid
transition. With both cell division and motility, the model is driven through
distinct liquid-hexatic and hexatic-liquid transitions with a re-entrant
state, or phase, diagram. While cell motility is typically thought to
disorder, we show that the combined effect of cell division and cell motility
allows access to the hexatic state. A key role in this process is played by
topological defects, both disclinations and dislocations.

#### Model

We model a 2D cell layer using the Self-Propelled Voronoi (SPV)Bi  _et al._
(2016) version of the vertex modelFarhadifar  _et al._ (2007); Li  _et al._
(2019, 2018); Yan and Bi (2019); Mitchel  _et al._ (2020); Das  _et al._
(2021). The cell shapes and the cellular network are determined based on the
Voronoi tesselationDirichlet (1850); Honda (1978) of the cell centers
$\\{\vec{r}_{i}\\}$. Here mechanical interactions in the tissue are controlled
by the energy functional
$E=\sum_{i=1}^{N}[K_{A}(A_{i}-A_{0})^{2}+K_{P}(P_{i}-P_{0})^{2}]$. The first
term, quadratic in the cell areas $\\{A_{i}\\}$, originates from the
incompressibility of cell volume, giving rise to a 2D area elasticity constant
$K_{A}$ and preferred area $A_{0}$Farhadifar  _et al._ (2007); Staple  _et
al._ (2010). The second term quadratic in the cell perimeters $\\{P_{i}\\}$
arises from the contractility of the cell cortex, with an elastic constant
$K_{P}$Farhadifar  _et al._ (2007). Here $P_{0}$ is the target cell
perimeterBi  _et al._ (2015), representing the interfacial tension set by the
competition between the cortical tension and the adhesion between adjacent
cellsStaple  _et al._ (2010). The target shape index
$p_{0}=P_{0}/\sqrt{A_{0}}$ effectively characterizes the competition between
cell-cell adhesion and cortical tension, acting as a signature for the solid-
liquid phase transitionBi  _et al._ (2015); Yan and Bi (2019). Apart from the
effective mechanical interaction force $\mathbf{F}_{i}=-\nabla_{i}E$, cells
are self-propelled. A self-propulsion force is exerted along the cell polarity
direction $\mathbf{\hat{n}_{i}}=(\cos\theta_{i},\sin\theta_{i})$, where
$\theta_{i}$ is the polarity angle. The self-propulsion has a constant
magnitude $v_{0}/\mu$, with the inverse of a frictional drag $\mu$. The
equation of motion for each cell is given by

$\dot{\vec{r_{i}}}=\mu\vec{F_{i}}+v_{0}\mathbf{\hat{n}_{i}}.$ (1)

The polarity angle obeys rotational diffusion: $d\theta_{i}/dt=\eta_{i}(t)$,
where $\eta_{i}(t)$ is white-noise
($\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=2D_{r}\delta(t-t^{\prime})\delta_{ij}$),
with $D_{r}$ the rotational diffusion rate.

In addition to the polarized self-propulsion, cell division and apoptosis
serve as another source of active forcing in living tissuesA. Matoz-Fernandez
_et al._ (2017); Doostmohammadi  _et al._ (2015); Cislo  _et al._ (2023); Saw
_et al._ (2017). In the SPV model, every cell has an equal division rate
$\gamma_{0}$. For each cell division, a daughter cell is introduced by
randomly seeding a point at a distance of $d=0.1$ (in units of the average
cell diameter) near the mother cell. In order to study the density-independent
effects of cell-division, we keep the number density of the tissue constant by
implementing apoptosis at the same rate as division. Apoptosis is then
performed on randomly chosen cells, which removes the cells from the tissue.
This simulation scheme mimics the maintenance of homeostatic balance in a
tissueBasan  _et al._ (2011); Ranft  _et al._ (2010).

The model can be nondimensionalized by expressing all lengths in units of
$\sqrt{\bar{A}}$, where $\bar{A}$ is the average cell area in the tissue and
time in units of $1/(\mu K_{A}\bar{A})$. Three independent parameters remain
the cell division/apoptosis rate $\gamma_{0}$, the magnitude of motility
$v_{0}$, and the cell shape index $p_{0}$. Throughout the simulations, we
choose $D_{r}=1$ without loss of generality.

The confluent tissue with $N$ cells is simulated in a square box with size
$L=\sqrt{N}$ under periodic boundary conditions. We numerically simulate the
model using the open-source software cellGPUSussman (2017). The simulations
start with a crystalline initial state in which cell centers form a triangular
lattice. Eq. 1 is numerically integrated for $2\times 10^{6}$ steps at a step
size of $\Delta t=0.05$. For all data presented, the analysis is based on the
steady-state regime of the simulations (final $5\times 10^{5}$ steps). In the
supplementary material (Fig.6 ), we also simulate the model starting from
amorphous states and do the cooling experiment to demonstrate that the results
are independent of the initial condition and simulation approach. We set
$p_{0}=3.6$ in our simulations. During the melting process, the tissues
undergo a transition from crystalline to hexatic to liquid as motility
increasesPasupalak  _et al._ (2020).

#### Signature for the emergence of hexatic order

Translational and orientational symmetries distinguish the three phases
crystalline, hexatic and liquid. A 2D crystalline phase has quasi-long-range
translational order and long-range orientational order while the liquid phase
has no long-range order of either kind. These two symmetries are related but
not concomitant. The system in the hexatic phase has no long-range
translational order but retains quasi-long-range orientational order Halperin
and Nelson (1978); Nelson and Halperin (1979).

We begin by quantifying translational and orientational order at the tissue
level. The translational order is quantified by
$\psi_{T}(\vec{r}_{j})=\exp{(i\vec{G}_{r}\vec{r}_{j})}$, where $\vec{G}_{r}$
represents a reciprocal vector in reciprocal space. The orientational order is
quantified
$\psi_{6}(\vec{r}_{j})=\left(1/\sum_{i=1}^{z_{j}}l_{ij}\right)\sum_{i=1}^{z_{j}}l_{ij}\exp{(i6\theta_{i}^{j})}$,
where the sum runs over the $n$ Voronoi neighbors of the cell and is weighted
by their shared edge lengthMickel  _et al._ (2013); Armengol-Collado  _et al._
(2022a, b). $\theta_{i}^{j}$ is the angle of the neighboring joint vector
$(\vec{r}_{i}-\vec{r}_{j})$ to a reference axis. In Fig. 1, we plot the
tissue-level order parameters
$\Psi_{6}=\frac{1}{N}\sum_{j=1}^{N}\psi_{6}(\vec{r}_{j})$ and
$\Psi_{T}=\frac{1}{N}\sum_{j=1}^{N}\psi_{T}(\vec{r}_{j})$ as a function of
$v_{0}$. In the absence of cell division (black lines), the tissue is a
crystal at low $v_{0}$ where both $\Psi_{T},\Psi_{6}$ are close to 1. The
order parameters decrease monotonically with increasing $v_{0}$. For
$0.35\lesssim v_{0}\lesssim 0.45$, the tissue lacks translational order but
retains orientational order, suggesting the existence of a hexatic phase
before melting into a liquid phase at higher $v_{0}$. This result is
consistent with the _crystal-hexatic-liquid_ melting scenario in the previous
study using a similar modelPasupalak  _et al._ (2020).

Figure 1: (a) The translational order parameter $\Psi_{T}$ and (b) the
orientational order parameter $\Psi_{6}$ as a function of the cell motility
$v_{0}$ at various division rates $\gamma_{0}$. The errorbar represents the
standard deviation over the ensemble of random simulations. (c) The snapshot
and structure factors $S(\bf{q})$ plotted for various $v_{0}$ at fixed
division rate $\gamma_{0}=2\times 10^{-5}$. Blue represents cells with
disclination charge $q_{i}=1$, red represents cells with $q_{i}=-1$, and dark
red represents cells with $q_{i}=-2$. The $v_{0}=0.35$ snapshot has
dislocation clusters but no disclinations.

When cells divide (color lines in Fig. 1), $\Psi_{T}$ is always close to zero
at any value of $v_{0}$. This clearly illustrates that activity due to cell
cycling (division/death) always destroys the translational order and therefore
forbids the formation of permanently frozen structuresA. Matoz-Fernandez  _et
al._ (2017). Remarkably, while an actively dividing tissue lacks translational
order, it retains orientational order for a large range of $v_{0}$ values.
This suggests the emergence of a hexatic phase at intermediate $v_{0}$ values.
A transition from liquid to hexatic to liquid is visualized by the structure
factor $S(\bf{q})$ for various $v_{0}$ at a fixed division rate.

Figure 2: The (a) translational and (b) bond-orientational correlation
functions at intermediate cell motility $v_{0}=0.35$. The (c) translational
and (d) bond-orientational correlation functions at low cell division rate
$\gamma_{0}=2\times 10^{-5}$.

In order to determine the location of the transitions between different
phases, we next compute the bond-orientational and translational correlation
functions. They are given by
$g_{\alpha}(r)=\langle\psi_{\alpha}^{*}(r)\psi_{\alpha}(0)\rangle$, with
$r=|\vec{r}_{i}-\vec{r}_{j}|$ and $\alpha=6,T$ corresponding to orientational
order and translational order, respectively. The peaks of correlations are
fitted by a power law decay $g_{\alpha}(r)\sim r^{-\eta_{\alpha}}$ (long-range
order) and an exponential decay $g_{\alpha}(r)\sim e^{-r/\xi_{6}}$ (short-
range order). KTHNY theoryHalperin and Nelson (1978); Nelson and Halperin
(1979); Young (1979); Kosterlitz and Thouless (1973); Bernard and Krauth
(2011) predicts $\eta_{6}=1/4$ at the hexatic-liquid transition point and
$\eta_{T}=1/3$ at the crystal-hexatic transition pointNelson (2002);
Deutschländer  _et al._ (2014).

The correlations are drawn and compared with reference power-laws
($\eta_{T}=1/3$ or $\eta_{6}=1/4$) in Fig.2 and Fig.7. Melting (without cell
division) allows quasi-long-range translational order at low $v_{0}$, decaying
as a power law with $\eta_{T}\leq 1/3$. The translational order with cell
division decays faster. Cell division also promotes the decay of bond-
orientational correlations, but the low $\gamma_{0}$ still allows for quasi-
long-range $g_{6}(r)$ with $\eta_{6}\leq 1/4$ at intermediate $v_{0}$ values.
A broken translational symmetry without broken orientational symmetry
characterizes the emergence of a hexatic state. Exponential decay fits the
orientational order better in both low- and high-$v_{0}$ liquid phases.

The fitted exponents $\eta_{6}$ and $\xi_{6}$ at fixed division rate
$\gamma_{0}=2\times 10^{-5}$ are shown in Fig.3 and in Fig. 8 for the case of
no cell division. These results confirm the emergence of two distinct liquid-
hexatic and hexatic-liquid transitions when there is cell division. The
correlations in the hexatic indeed display quasi-long-range order, well-fitted
by power-law decays, $g_{6}(r)\sim r^{-\eta_{6}}$, while outside the hexatic
region correlations are short-range and well-fitted by exponential decays
$g_{6}(r)\sim e^{-r/\xi_{6}}$. As the hexatic phase is approached from either
side, $\xi_{6}$ grows rapidly, consistent with a diverging correlation length.

Despite excellent agreement with the KTHNY model, the correlation functions
and the associated quantities ($\xi_{6},\eta_{6}$) near the onset of hexatic
states suffer from large sample-to-sample variations as shown in Fig. 3(a,b).
We have confirmed that this is not due to finite-size effects since, even at
large system sizes, the behavior of $g_{6}(r)$ can range from exponential
decay to a power-law decay (Fig. 9). Consequently, ($\xi_{6}$, $\eta_{6}$)
cannot be used to pinpoint the precise location of the liquid-hexatic and
hexatic-liquid transitions.

Figure 3:  (a) The correlation length $\xi_{6}$ and (b) the power-law decay
exponent $\eta_{6}$ of the orientational correlation function are shown as
functions of $v_{0}$ at constant $\gamma_{0}=2\times 10^{-5}$. Circles
represent the fitting exponents for different seeds, and the solid lines
average the seeds. (c) The hexatic order parameter susceptibility $\chi_{6}$.
(d) Phase diagram as a function of cell division rate $\gamma_{0}$ and
motility $v_{0}$.

We next take advantage of the large fluctuations that arise near critical
points by using the order parameter susceptibility to pinpoint the
transitions. The susceptibility is given by
$\chi_{\alpha}=N(\langle|\Psi_{\alpha}|^{2}\rangle-\langle|\Psi_{\alpha}|\rangle^{2})$,
which characterizes the fluctuations in the translational ($\chi_{T}$) and
orientational ($\chi_{6}$) order parameters. Since $\chi_{\alpha}$ is
essentially an integral of the correlation function, it is expected to be more
robust to finite-size or finite-time effectsHan  _et al._ (2008); Durand and
Heu (2019).

In the melting process without cell division (shown in Fig.10), there is a
sharp divergence of $\chi_{T}$ at $v_{0}=0.35$, indicating a crystal-hexatic
transition. On the other hand, $\chi_{6}$ diverges at $v_{0}=0.46\pm 0.01$,
which corresponds to the hexatic-liquid transition. By analyzing system sizes
from $N=2430$ to $38880$, we show that these divergences are robust to finite-
size effects (Supplemental Material).

In contrast, $\chi_{6}$ with cell division (at $\gamma_{0}=2\times 10^{-5}/s$)
generate two peaks ( Fig.3(c)). The divergence of $\chi_{6}$ determines two
distinct transition points at $v_{0}=0.25\pm 0.01$ and at $v_{0}=0.45\pm
0.01$. Whereas the second point is a vestige of the hexatic-liquid transition
in the absence of cell division, the first transition point arises solely from
cell division. Here, a state that would be crystalline in the absence of cell
division becomes hexatic when cells divide.

Exploring various cell division rates $\gamma_{0}$ and active motilities
$v_{0}$ at fixed $p_{0}=3.6$, we plot the $v_{0}-\gamma_{0}$ phase diagram in
Fig.3(d). Color indicates the mean magnitude of the global orientational order
over tens of thousands of frames. Black dots mark the peaks of $\chi_{6}$ at
various division rates. The two transition points approach each other and
disappear as the division rate increases. The blue dashed lines indicate the
liquid-hexatic and hexatic-crystal transition points in the absence of cell
division. We also investigate the $p_{0}$ dependence of the phase diagram in
Fig.11 and compare it to previous results by Pasupalak et al.Pasupalak  _et
al._ (2020) which do not have cell division.

#### Disclinations and Dislocations

Figure 4:  (a) The volume densities of dislocations and disclinations are
plotted as functions of $v_{0}$ at constant cell division rate
$\gamma_{0}=2\times 10^{-5}$. (b) The same quantities are plotted at a
constant $v_{0}=0.35$ and varying $\gamma_{0}$.

According to KTHNY theoryHalperin and Nelson (1978); Nelson and Halperin
(1979); Young (1979); Kosterlitz and Thouless (1973); Bernard and Krauth
(2011), the distinct phases crystalline, hexatic and liquid are characterized
by the distributions of the basic topological defects known as disclinations
and dislocations. Whereas the pure crystalline phase is defect free, or
equivalently all defects are tightly bound in defect-antidefect pairs, the
hexatic phase has a non-vanishing density of free dislocations and the liquid
phase has a non-vanishing density of free disclinations.

We can define a charge $q_{i}=6-z_{i}$Bowick and Giomi (2009) associated with
disclinations, where $z_{i}$ is the coordination number (number of neighbors)
of the $i$th cell. Hexagonal cells are thus "neutral", pentagonal cells have
charge +1, heptagonal cells charge -1 and so on. Dislocations, the defects
that disrupt translational order but preserve orientational order, correspond
to tightly bound $5-7$ pairs. They are neutral as disclinations but possess a
net vectorial charge, the Burgers vector. We approximate the Burgers vector by
the displacement vector separating the 5 and the connected 7. In general there
will be clusters of connected defects and one must measure the associated
disclination and dislocation charges of the entire cluster. The density of
disclinations and dislocations are calculated by their volume fraction
averaged over time.

As shown in Fig.4(a), cell division creates dislocations at a rate dependent
on motility. Division tends to disorder, favoring a liquid. What about
motility? At low motility, division disordering wins. At high motility both
processes generate disorder, leading again to a liquid. But for a significant
range of intermediate motilities, we see that the number density of free
disclinations falls to zero whereas the free dislocation density is finite.
How is this possible? In this intermediate region we hypothesize that
disclinations are able to explore sufficient configuration space to access
local free energy minima at which all disclinations find their anti-
disclinations and bind into dislocations, thus leading to a hexatic. Fig.4(b)
shows the defect density dependence on cell division rate at a fixed motility
in the hexatic regime, showing that sufficiently high division rates lead to a
non-zero density of free disclinations, thus melting the hexatic to a liquid.
The middle snapshot in Fig.1(c) is a representative snapshot of a hexatic
state ($v_{0}=0.35$, $\gamma_{0}=2\times 10^{-5}$). Note the presence of
dislocation complexes but no isolated disclinations. Movie. S1shows a dynamic
evolution of states with various values of cell motility at a fixed division
rate. The densities of dislocations and disclinations are indicated by color
as a function of cell division rate $\gamma_{0}$ and motility $v_{0}$ at fixed
$p_{0}=3.6$ in Fig. 12. Black dots mark the same data in Fig. 3(d).

#### Meanfield Model

Figure 5: A meanfield description for defect dynamics in a tissue. (a) The
energy landscape and transition diagram between states. (b) The phase diagram
of meanfield model as the function of the dimensionless division rate
$\gamma_{0}/R$ and motility $(T/\Delta\epsilon_{X})^{1/2}$. The green(red)
dashed line indicates the liquid-hexatic and hexatic-crystal transition in the
absence of cell division.

To further understand the emergence of hexatic order through cell division we
develop a simple meanfield model (MF) (Fig.5) incorporating the competition
between cell division and motility. We simplify the state of a small cell
cluster ($\sim 4$ cells) using a meanfield approximation that allows three
states: (a) crystalline solid state ("ordered"), (b) an isolated single
dislocation, and (c) an isolated single disinclination. Transitions between
states arise from fluctuations over the energy barriers $\Delta\epsilon_{i}$,
as illustrated in Fig. 5(a). Fluctuations arise from both Brownian motility
forces and cell division, leading, in the low temperature/velocity limit, to
an equal density distribution of states rather than only the "ordered" state
(see Supplementary Text). In the steady-state limit(Fig.5(b)), the asymmetric
boundary of the hexatic region (determined by the threshold of $\rho_{H}$ and
$\rho_{D}$) is remarkably similar to our simulations and the phase diagram
shows a re-entrant liquid-hexatic-liquid transition with changing motility, as
in Fig. 3(d). The MF model predicts unique behavior for ultra-low division
rates ($\gamma_{0}/R\ll 10^{-10}$), where tissues undergo a complex
temperature-dependent transition sequence, following a Liquid-Hexatic-Crystal-
Hexatic-Liquid path at constant $\gamma_{0}/R$ (Fig.5(b) and also see
Supplementary Text for a detailed discussion). Remarkably, our phase diagram
closely mirrors that of 2D melting on a random substrate Sachdev and Nelson
(1984). In both models, temperature drives phase transitions, and in our case,
cell division plays a role analogous to substrate disorder, introducing
persistent, random spatial distortions.

#### Discussion

The subtle balance required to establish hexatic order in equilibrium means
that it is often confined to a rather narrow region of the relevant parameter
space. Our findings suggest that cell division provides a new way of exploring
the configuration space of physical systems, as noted above. In particular,
the dynamics of dislocation defects generated by cell-division, both self-
propelled and relaxational, promote fluctuations over barriers separating the
hexatic phase from crystalline or liquid phases. This phenomenon, which we may
call defect-driven structure development, may well have implications beyond
biological systems. In terms of the configuration space explored by the vertex
model, cell division and apoptosis correspond to adding T2 moves (or
interstitial insertion/deletion) to the allowed lattice updates – this yields
a more efficient exploration of the space of all Voronoi tesselations and thus
better routes to local hexatic minima Bowick  _et al._ (2007); Irvine  _et
al._ (2012); J. Bowick  _et al._ (2007). It is remarkable that the early work
of Swope and AndersenSwope and Andersen (1995) found the hexatic phase by
employing a grand canonical ensemble in which particles are added and removed.
The mechanism we find here is very different from that found in colloidsHan
_et al._ (2008) and models of active particlesDigregorio  _et al._ (2018),
where packing density plays a crucial role.

We have taken cell division to be isotropic. The inclusion of oriented cell
divisions, however, would only enhance hexatic order. Recent workCislo  _et
al._ (2023) has shown that oriented cell divisions can give rise to novel
four-fold orientational order in vivo through active defect climb, where
defects introduced into the nascent lattice by cell divisions are healed by
subsequent divisions along a well-defined global polarity axis. The effect of
oriented divisions on hexatic order is a subject for the future.

###### Acknowledgements.

This work was supported in part by NSF DMR-2046683 (Y.T. and D. B.),
PHY-1748958 (D. B. and M. J. B.), the Center for Theoretical Biological
Physics NSF PHY-2019745 (Y. T. and D. B.), Alfred P. Sloan Foundation (Y. T.
and D. B.) and The Human Frontier Science Program (Y. T. and D. B.))

## References

  * Lecuit  _et al._ (2011) Thomas Lecuit, Pierre-François Lenne, and Edwin Munro, “Force generation, transmission, and integration during cell and tissue morphogenesis,” Annu. Rev. Cell Dev. Biol. 27, 157–184 (2011).
  * Janmey and Miller (2011) Paul A. Janmey and R. Tyler Miller, “Mechanisms of mechanical signaling in development and disease,” J Cell Sci 124, 9–18 (2011).
  * Mammoto and Ingber (2010) Tadanori Mammoto and Donald E. Ingber, “Mechanical control of tissue and organ development,” Development 137, 1407–1420 (2010).
  * Tenney and Discher (2009) Rebeca M Tenney and Dennis E Discher, “Stem cells, microenvironment mechanics, and growth factor activation,” Current Opinion in Cell Biology Cell-to-Cell Contact and Extracellular Matrix, 21, 630–635 (2009).
  * Atia  _et al._ (2018) Lior Atia, Dapeng Bi, Yasha Sharma, Jennifer A. Mitchel, Bomi Gweon, Stephan A. Koehler, Stephen J. DeCamp, Bo Lan, Jae Hun Kim, Rebecca Hirsch, Adrian F. Pegoraro, Kyu Ha Lee, Jacqueline R. Starr, David A. Weitz, Adam C. Martin, Jin-Ah Park, James P. Butler, and Jeffrey J. Fredberg, “Geometric constraints during epithelial jamming,” Nature Phys 14, 613–620 (2018).
  * Park  _et al._ (2015) Jin-Ah Park, Jae Hun Kim, Dapeng Bi, Jennifer A. Mitchel, Nader Taheri Qazvini, Kelan Tantisira, Chan Young Park, Maureen McGill, Sae-Hoon Kim, Bomi Gweon, Jacob Notbohm, Robert Steward Jr, Stephanie Burger, Scott H. Randell, Alvin T. Kho, Dhananjay T. Tambe, Corey Hardin, Stephanie A. Shore, Elliot Israel, David A. Weitz, Daniel J. Tschumperlin, Elizabeth P. Henske, Scott T. Weiss, M. Lisa Manning, James P. Butler, Jeffrey M. Drazen, and Jeffrey J. Fredberg, “Unjamming and cell shape in the asthmatic airway epithelium,” Nature Mater 14, 1040–1048 (2015).
  * Mongera  _et al._ (2018) Alessandro Mongera, Payam Rowghanian, Hannah J. Gustafson, Elijah Shelton, David A. Kealhofer, Emmet K. Carn, Friedhelm Serwane, Adam A. Lucio, James Giammona, and Otger Campàs, “A fluid-to-solid jamming transition underlies vertebrate body axis elongation,” Nature 561, 401–405 (2018).
  * Malinverno  _et al._ (2017) Chiara Malinverno, Salvatore Corallino, Fabio Giavazzi, Martin Bergert, Qingsen Li, Marco Leoni, Andrea Disanza, Emanuela Frittoli, Amanda Oldani, Emanuele Martini, Tobias Lendenmann, Gianluca Deflorian, Galina V. Beznoussenko, Dimos Poulikakos, Kok Haur Ong, Marina Uroz, Xavier Trepat, Dario Parazzoli, Paolo Maiuri, Weimiao Yu, Aldo Ferrari, Roberto Cerbino, and Giorgio Scita, “Endocytic reawakening of motility in jammed epithelia,” Nature Mater 16, 587–596 (2017).
  * Angelini  _et al._ (2011) Thomas E. Angelini, Edouard Hannezo, Xavier Trepat, Manuel Marquez, Jeffrey J. Fredberg, and David A. Weitz, “Glass-like dynamics of collective cell migration,” Proceedings of the National Academy of Sciences 108, 4714–4719 (2011).
  * Garcia  _et al._ (2015) Simon Garcia, Edouard Hannezo, Jens Elgeti, Jean-François Joanny, Pascal Silberzan, and Nir S. Gov, “Physics of active jamming during collective cellular motion in a monolayer,” Proceedings of the National Academy of Sciences 112, 15314–15319 (2015).
  * Puliafito  _et al._ (2012) Alberto Puliafito, Lars Hufnagel, Pierre Neveu, Sebastian Streichan, Alex Sigal, D. Kuchnir Fygenson, and Boris I. Shraiman, “Collective and single cell behavior in epithelial contact inhibition,” Proceedings of the National Academy of Sciences 109, 739–744 (2012).
  * Bi  _et al._ (2015) Dapeng Bi, J. H. Lopez, J. M. Schwarz, and M. Lisa Manning, “A density-independent rigidity transition in biological tissues,” Nature Phys 11, 1074–1079 (2015).
  * Bi  _et al._ (2016) Dapeng Bi, Xingbo Yang, M. Cristina Marchetti, and M.L. Manning, “Motility-driven glass and jamming transitions in biological tissues,” Physical Review X 6, 021011 (2016).
  * Halperin and Nelson (1978) B. I. Halperin and David R. Nelson, “Theory of two-dimensional melting,” Phys. Rev. Lett. 41, 121–124 (1978).
  * Nelson and Halperin (1979) David R. Nelson and B. I. Halperin, “Dislocation-mediated melting in two dimensions,” Phys. Rev. B 19, 2457–2484 (1979).
  * Chou  _et al._ (1998) Chia-Fu Chou, Anjun J. Jin, S. W. Hui, C. C. Huang, and John T. Ho, “Multiple-step melting in two-dimensional hexatic liquid-crystal films,” Science 280, 1424–1426 (1998).
  * Murray (1992) Cherry A. Murray, “Experimental studies of melting and hexatic order in two-dimensional colloidal suspensions,” Bond-Orientational Order in Condensed Matter Systems , 137–215 (1992).
  * Zahn  _et al._ (1999) K. Zahn, R. Lenke, and G. Maret, “Two-stage melting of paramagnetic colloidal crystals in two dimensions,” Phys. Rev. Lett. 82, 2721–2724 (1999).
  * Zahn and Maret (2000) K. Zahn and G. Maret, “Dynamic criteria for melting in two dimensions,” Phys. Rev. Lett. 85, 3656–3659 (2000).
  * Han  _et al._ (2008) Y. Han, N. Y. Ha, A. M. Alsayed, and A. G. Yodh, “Melting of two-dimensional tunable-diameter colloidal crystals,” Phys. Rev. E 77, 041406 (2008).
  * Seshadri and Westervelt (1991) R. Seshadri and R. M. Westervelt, “Hexatic-to-liquid melting transition in two-dimensional magnetic-bubble lattices,” Phys. Rev. Lett. 66, 2774–2777 (1991).
  * Knobler and Desai (1992) Charles M. Knobler and Rashmi C. Desai, “Phase transitions in monolayers,” Annual Review of Physical Chemistry 43, 207–236 (1992).
  * Huang (1992) C. C. Huang, “Nature of phase transitions related to stacked hexatic phases in liquid crystals,” in  _Bond-Orientational Order in Condensed Matter Systems_, Partially Ordered Systems, edited by Katherine J. Strandburg (Springer, New York, NY, 1992) pp. 78–136.
  * Huang and Stoebe (1993) C.c. Huang and T. Stoebe, “Thermal properties of ‘stacked hexatic phases’ in liquid crystals,” Advances in Physics 42, 343–391 (1993).
  * Podgornik  _et al._ (1996) R Podgornik, H H Strey, K Gawrisch, D C Rau, A Rupprecht, and V A Parsegian, “Bond orientational order, molecular motion, and free energy of high-density dna mesophases.” Proceedings of the National Academy of Sciences 93, 4261–4266 (1996).
  * Van Winkle  _et al._ (1997) D. H. Van Winkle, Arunava Chatterjee, Robert Link, and R. L. Rill, “Magnetic-field alignment of cholesteric liquid-crystalline dna,” Phys. Rev. E 55, 4354–4359 (1997).
  * Rill  _et al._ (1991) Randolph L. Rill, Teresa E. Strzelecka, Michael W. Davidson, and David H. Van Winkle, “Ordered phases in concentrated dna solutions,” Physica A: Statistical Mechanics and its Applications 176, 87–116 (1991).
  * Bishop  _et al._ (1994) David J. Bishop, Peter L. Gammel, and Cherry A. Murray, “Magnetic decoration studies of flux line lattices in the cuprate superconductors,” in  _The Vortex State_, NATO ASI Series, edited by Nicole Bontemps, Yvan Bruynseraede, Guy Deutscher, and Aharon Kapitulnik (Springer Netherlands, Dordrecht, 1994) pp. 99–123.
  * Murray  _et al._ (1990) C. A. Murray, P. L. Gammel, D. J. Bishop, D. B. Mitzi, and A. Kapitulnik, “Observation of a hexatic vortex glass in flux lattices of the high-$\mathit{T}_{\mathit{c}}$ superconductor $\mathrm{Bi}_{\mathrm{2.1}}\mathrm{Sr}_{\mathrm{1.9}}\mathrm{Ca}_{\mathrm{0.9}}\mathrm{Cu}_{\mathrm{2}}\mathrm{O}_{\mathrm{8+\delta}}$,” Phys. Rev. Lett. 64, 2312–2315 (1990).
  * Li and Ciamarra (2018) Yan-Wei Li and Massimo Pica Ciamarra, “Role of cell deformability in the two-dimensional melting of biological tissues,” Physical Review Materials 2, 045602 (2018).
  * Durand and Heu (2019) Marc Durand and Julien Heu, “Thermally driven order-disorder transition in two-dimensional soft cellular systems,” Phys. Rev. Lett. 123, 188001 (2019).
  * xue Guo  _et al._ (2023) Rui xue Guo, Jia jian Li, and Bao quan Ai, “Melting of two-dimensional deformable particle systems,” Physica A: Statistical Mechanics and its Applications 623, 128833 (2023).
  * Pasupalak  _et al._ (2020) Anshuman Pasupalak, Li Yan-Wei, Ran Ni, and Massimo Pica Ciamarra, “Hexatic phase in a model of active biological tissues,” Soft Matter 16, 3914–3920 (2020).
  * Loewe  _et al._ (2020) Benjamin Loewe, Michael Chiang, Davide Marenduzzo, and M. Cristina Marchetti, “Solid-liquid transition of deformable and overlapping active particles,” Phys. Rev. Lett. 125, 038003 (2020).
  * Paoluzzi  _et al._ (2021) Matteo Paoluzzi, Luca Angelani, Giorgio Gosti, M. Cristina Marchetti, Ignacio Pagonabarraga, and Giancarlo Ruocco, “Alignment interactions drive structural transitions in biological tissues,” Phys. Rev. E 104, 044606 (2021).
  * jian Li and quan Ai (2021) Jia jian Li and Bao quan Ai, “Melting of two-dimensional biological tissues containing active ornstein–uhlenbeck cells,” New Journal of Physics 23, 083044 (2021).
  * Cohen  _et al._ (2020) Roie Cohen, Liat Amir-Zilberstein, Micha Hersch, Shiran Woland, Olga Loza, Shahar Taiber, Fumio Matsuzaki, Sven Bergmann, Karen B. Avraham, and David Sprinzak, “Mechanical forces drive ordered patterning of hair cells in the mammalian inner ear,” Nat Commun 11, 5137 (2020).
  * Hilgenfeldt  _et al._ (2008) Sascha Hilgenfeldt, Sinem Erisken, and Richard W. Carthew, “Physical modeling of cell geometric order in an epithelial tissue,” Proceedings of the National Academy of Sciences 105, 907–911 (2008), https://www.pnas.org/doi/pdf/10.1073/pnas.0711077105 .
  * A. Matoz-Fernandez  _et al._ (2017) D. A. Matoz-Fernandez, Kirsten Martens, Rastko Sknepnek, J. L. Barrat, and Silke Henkes, “Cell division and death inhibit glassy behaviour of confluent tissues,” Soft Matter 13, 3205–3212 (2017).
  * Mitchel  _et al._ (2020) Jennifer A. Mitchel, Amit Das, Michael J. O’Sullivan, Ian T. Stancil, Stephen J. DeCamp, Stephan Koehler, Oscar H. Ocaña, James P. Butler, Jeffrey J. Fredberg, M. Angela Nieto, Dapeng Bi, and Jin-Ah Park, “In primary airway epithelial cells, the unjamming transition is distinct from the epithelial-to-mesenchymal transition,” Nat Commun 11, 5053 (2020).
  * Kanesaki  _et al._ (2011) Takuma Kanesaki, Carina M. Edwards, Ulrich S. Schwarz, and Jörg Grosshans, “Dynamic ordering of nuclei in syncytial embryos: A quantitative analysis of the role of cytoskeletal networks,” Integrative Biology 3, 1112–1119 (2011).
  * Wu  _et al._ (2021) Xiaoxuan Wu, Kakit Kong, Wenlei Xiao, and Feng Liu, “Attractive internuclear force drives the collective behavior of nuclear arrays in drosophila embryos,” PLOS Computational Biology 17, e1009605 (2021).
  * Kaiser  _et al._ (2018) Franz Kaiser, Zhiyi Lv, Daniel Marques Rodrigues, Jan Rosenbaum, Timo Aspelmeier, Jörg Großhans, and Karen Alim, “Mechanical model of nuclei ordering in drosophila embryos reveals dilution of stochastic forces,” Biophysical Journal 114, 1730–1740 (2018).
  * Farhadifar  _et al._ (2007) Reza Farhadifar, Jens-Christian Röper, Benoit Aigouy, Suzanne Eaton, and Frank Jülicher, “The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing,” Current Biology 17, 2095–2104 (2007).
  * Li  _et al._ (2019) Xinzhi Li, Amit Das, and Dapeng Bi, “Mechanical heterogeneity in tissues promotes rigidity and controls cellular invasion,” Phys. Rev. Lett. 123, 058101 (2019).
  * Li  _et al._ (2018) Xinzhi Li, Amit Das, and Dapeng Bi, “Biological tissue-inspired tunable photonic fluid,” Proceedings of the National Academy of Sciences 115, 6650–6655 (2018).
  * Yan and Bi (2019) Le Yan and Dapeng Bi, “Multicellular rosettes drive fluid-solid transition in epithelial tissues,” Phys. Rev. X 9, 011029 (2019).
  * Das  _et al._ (2021) Amit Das, Srikanth Sastry, and Dapeng Bi, “Controlled neighbor exchanges drive glassy behavior, intermittency, and cell streaming in epithelial tissues,” Phys. Rev. X 11, 041037 (2021).
  * Dirichlet (1850) G. Lejeune Dirichlet, “über die reduction der positiven quadratischen formen mit drei unbestimmten ganzen zahlen.” Journal für die reine und angewandte Mathematik (Crelles Journal) 1850, 209–227 (1850).
  * Honda (1978) Hisao Honda, “Description of cellular patterns by dirichlet domains: The two-dimensional case,” Journal of Theoretical Biology 72, 523–543 (1978).
  * Staple  _et al._ (2010) D. B. Staple, R. Farhadifar, J. C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, “Mechanics and remodelling of cell packings in epithelia,” Eur. Phys. J. E 33, 117–127 (2010).
  * Doostmohammadi  _et al._ (2015) Amin Doostmohammadi, Sumesh P. Thampi, Thuan B. Saw, Chwee T. Lim, Benoit Ladoux, and Julia M. Yeomans, “Celebrating soft matter’s 10th anniversary: Cell division: A source of active stress in cellular monolayers,” Soft Matter 11, 7328–7336 (2015).
  * Cislo  _et al._ (2023) Dillon J. Cislo, Fengshuo Yang, Haodong Qin, Anastasios Pavlopoulos, Mark J. Bowick, and Sebastian J. Streichan, “Active cell divisions generate fourfold orientationally ordered phase in living tissue,” Nat. Phys. 19, 1201–1210 (2023).
  * Saw  _et al._ (2017) Thuan Beng Saw, Amin Doostmohammadi, Vincent Nier, Leyla Kocgozlu, Sumesh Thampi, Yusuke Toyama, Philippe Marcq, Chwee Teck Lim, Julia M. Yeomans, and Benoit Ladoux, “Topological defects in epithelia govern cell death and extrusion,” Nature 544, 212–216 (2017).
  * Basan  _et al._ (2011) Markus Basan, Jacques Prost, Jean-François Joanny, and Jens Elgeti, “Dissipative particle dynamics simulations for biological tissues: Rheology and competition,” Phys. Biol. 8, 026014 (2011).
  * Ranft  _et al._ (2010) Jonas Ranft, Markus Basan, Jens Elgeti, Jean-François Joanny, Jacques Prost, and Frank Jülicher, “Fluidization of tissues by cell division and apoptosis,” Proceedings of the National Academy of Sciences 107, 20863–20868 (2010).
  * Sussman (2017) Daniel M. Sussman, “cellgpu: Massively parallel simulations of dynamic vertex models,” Computer Physics Communications 219, 400–406 (2017).
  * Mickel  _et al._ (2013) Walter Mickel, Sebastian C. Kapfer, Gerd E. Schröder-Turk, and Klaus Mecke, “Shortcomings of the bond orientational order parameters for the analysis of disordered particulate matter,” J. Chem. Phys. 138, 044501 (2013).
  * Armengol-Collado  _et al._ (2022a) Josep-Maria Armengol-Collado, Livio Nicola Carenza, and Luca Giomi, “Hydrodynamics and multiscale order in confluent epithelia,”  (2022a), arxiv:arXiv:2202.00651 .
  * Armengol-Collado  _et al._ (2022b) Josep-Maria Armengol-Collado, Livio Nicola Carenza, Julia Eckert, Dimitrios Krommydas, and Luca Giomi, “Epithelia are multiscale active liquid crystals,”  (2022b), arxiv:arXiv:2202.00668 .
  * Young (1979) A. P. Young, “Melting and the vector coulomb gas in two dimensions,” Phys. Rev. B 19, 1855–1866 (1979).
  * Kosterlitz and Thouless (1973) John Michael Kosterlitz and David James Thouless, “Ordering, metastability and phase transitions in two-dimensional systems,” Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Bernard and Krauth (2011) Etienne P. Bernard and Werner Krauth, “Two-step melting in two dimensions: First-order liquid-hexatic transition,” Phys. Rev. Lett. 107, 155704 (2011).
  * Nelson (2002) David R. Nelson,  _Defects and Geometry in Condensed Matter Physics_ (Cambridge University Press, 2002).
  * Deutschländer  _et al._ (2014) Sven Deutschländer, Antonio M. Puertas, Georg Maret, and Peter Keim, “Specific heat in two-dimensional melting,” Phys. Rev. Lett. 113, 127801 (2014).
  * Bowick and Giomi (2009) Mark J. Bowick and Luca Giomi, “Two-dimensional matter: Order, curvature and defects,” Advances in Physics 58, 449–563 (2009).
  * Sachdev and Nelson (1984) S Sachdev and D R Nelson, “Crystalline and fluid order on a random topography,” J. Phys. C: Solid State Phys. 17, 5473–5489 (1984).
  * Bowick  _et al._ (2007) Mark Bowick, Homin Shin, and Alex Travesset, “Dynamics and instabilities of defects in two-dimensional crystals on curved backgrounds,” Phys. Rev. E 75, 021404 (2007).
  * Irvine  _et al._ (2012) William T. M. Irvine, Mark J. Bowick, and Paul M. Chaikin, “Fractionalization of interstitials in curved colloidal crystals,” Nature Mater 11, 948–951 (2012).
  * J. Bowick  _et al._ (2007) Mark J. Bowick, David R. Nelson, and Homin Shin, “Interstitial fractionalization and spherical crystallography,” Physical Chemistry Chemical Physics 9, 6304–6312 (2007).
  * Swope and Andersen (1995) William C. Swope and Hans C. Andersen, “A computer simulation method for the calculation of chemical potentials of liquids and solids using the bicanonical ensemble,” The Journal of Chemical Physics 102, 2851–2863 (1995).
  * Digregorio  _et al._ (2018) Pasquale Digregorio, Demian Levis, Antonio Suma, Leticia F. Cugliandolo, Giuseppe Gonnella, and Ignacio Pagonabarraga, “Full phase diagram of active brownian disks: From melting to motility-induced phase separation,” Phys. Rev. Lett. 121, 098003 (2018).
  * Bi  _et al._ (2014) Dapeng Bi, Jorge H. Lopez, J. M. Schwarz, and M. Lisa Manning, “Energy barriers and cell migration in densely packed tissues,” Soft Matter 10, 1885–1890 (2014).

##  Appendix A Supplementary Text

##  Appendix B Susceptibility of Order Parameters

$\chi_{6}$ is a measure of the fluctuation in tens of thousands of frames,
which come from 50 simulations. Each simulation has a different seed,
controlling the initial velocity direction and dividing cells. $\chi_{T}$ is a
measure of the fluctuations of the time-averaged $\Psi_{T}$ from the 20
simulations. The ensemble average plays a key role in the susceptibility
calculation.

##  Appendix C  Tracking Topological Defects

The video (Movie. S1) shows a dynamic evolution of states with various cell
motilities ($v_{0}=0.15,0.35,0.55$) with a fixed division rate
($\gamma_{0}=2\times 10^{-5}$). We color the dislocations in cyan($z_{i}<6$)
and magenta($z_{i}>6$) and color the disclinations in blue($z_{i}<6$) and
red($z_{i}>6$). At $v_{0}=0.15$(left), the fluctuations from motility are
insufficient to anneal the dislocations created by cell division. Dislocations
eventually unbind into disclinations. The tissue is in a liquid state.
$v_{0}=0.35$(middle) allows dislocations to reorganize and overcome barriers,
leading to a hexatic. In contrast $v_{0}=0.55$(right) generates the
dislocations and disclinations by itself.

Figure 6: The orientational order parameter is independent of simulation
protocols. (a) The orientational order parameter remains consistent across
different initial conditions. The crystal (solid lines) and random (dash
lines) initial conditions both lead to the same form for $\Psi_{6}(v_{0})$.
$\Psi_{6}$ is non-zero at intermediate $v_{0}$ values and declines for both
higher and lower $v_{0}$. (b) The orientational order parameter is unaffected
by the approach taken during the simulation process. The cooling run (green
solid line) has the same $\Psi_{6}(v_{0})$ function as the direct run (blue
dash lines) from crystal initial conditions at constant division rate
$\gamma_{0}=4\times 10^{-5}$. In the cooling run, we start from the simulation
for $v_{0}=0.60$. Subsequently, we use the final state of $v_{0}=0.6$ as the
initial state for $v_{0}=0.55$ and employ the resulting state of $v_{0}=0.55$
as the starting state for $v_{0}=0.50$, continuing this sequence until we
reach $v_{0}=0.10$. Figure 7: The (a) translational and (b) bond-orientational
correlation functions in the absence of cell division. The curves evolve from
a power-law decay to an exponential decay as $v_{0}$ increases. (a) The
reference line $g_{T}\sim r^{-1/3}$ indicates the theoretically expected for
the crystal-hexatic transition, while (b) the $g_{6}\sim r^{-1/4}$ fall-off
indicates the expected decay for a hexatic-liquid transition. The decay of
$g_{6}$ lags $g_{T}$, leaving an intermediate range of $v_{0}$ without quasi-
long-range translational order but with quasi-long-range orientational order.
Figure 8: The correlation length $\xi_{\alpha}$ and exponent $\eta_{\alpha}$
($\alpha=6,T$) as a function of $v_{0}$ in the absence of cell division. The
translational correlations in the crystal phase are quasi-long-range with a
power-law decay $g_{T}(r)\sim r^{-\eta_{T}}$, and those in the hexatic and
liquid range are short-range with an exponential decay $g_{T}(r)\sim
e^{-r/\xi_{T}}$. The orientational correlations in the hexatic are quasi-long-
range with a power-law decay $g_{6}(r)\sim r^{-\eta_{6}}$, and those in the
liquid are short-range with an exponential decay $g_{6}(r)\sim
e^{-r/\xi_{6}}$. $\xi_{6}$ and $\eta_{6}$ grow significantly at the transition
point.  Figure 9: The finite-size effects for the translational and
orientational correlations. Different colors represent various $N$ and
different lines indicate simulations with distinct seeds. (a) Near the
transition point ($v_{0}=0.25$ and $\gamma_{0}=2\times 10^{-5}$), the behavior
of $g_{6}(r)$ suffers large sample-to-sample fluctuations, ranging from
exponential decay to a power-law decay. (b) Deep in the hexatic phase
($v_{0}=0.35$ and $\gamma_{0}=2\times 10^{-5}$), the finite-size analysis over
a large range of $N$ confirms quasi-long-range order. The power-law decay is
independent of system size and seeds.  Figure 10: The susceptibility during
melting as a function of $v_{0}$ for various system sizes. (a) All the
$\chi_{T}$ in melting have a sharp peak at $v_{0}=0.35$, clearly indicating
the crystal-hexatic transition. (b) All the $\chi_{6}$ have a sharp peak at
$v_{0}=0.46\pm 0.01$, clearly indicating the hexatic-liquid transition. Figure
11: Depedence of the orientational order parameter $\Psi_{6}$ on different
choices of $p_{0}$. (a) The phase diagram as a function of cell shape index
$p_{0}$ and motility $v_{0}$ with $\gamma_{0}=2\times 10^{-5}$. The color
indicates the magnitude of $\Psi_{6}$. (b) The orientational order parameter
$\Psi_{6}$ as the function of $v_{0}$ with various $p_{0}$ and fixed
$\gamma_{0}=2\times 10^{-5}$. The high $\Psi_{6}$ values at intermediate
$v_{0}$ levels indicate the appearance of the hexatic phase. Cell division
leads to a liquid-hexatic-liquid transition sequence across different $p_{0}$
values. The peak of $\Psi_{6}$ shrinks and shifts to smaller $v_{0}$ as
$p_{0}$ increasingly approaches the $p_{0}\sim 3.75$, which is in harmony with
the no-division casePasupalak  _et al._ (2020).  Figure 12: The density of (a)
disclinations and (b) dislocations as a function of cell division rate
$\gamma_{0}$ and motility $v_{0}$. Color indicates the magnitude of the
density of (a) disclinations and (b) dislocations. White dots are the hexatic-
liquid transition points, obtained from the divergence of $\chi_{6}$. In the
hexatic phase, the density of disclinations is very low. In the liquid phase,
the density of disclinations increases rapidly at low $v_{0}$.

##  Appendix D Meanfield model

The mean-field approximation dynamics of a single small ($\sim 4$ cells) cell
cluster can be classified into three states: (1) a fully ordered crystalline
state (X), (2) an isolated dislocation corresponding to a hexatic state (H),
(3) an isolated single disinclination corresponding to a liquid state (D). The
corresponding volume densities (or probabilities) are given by $\rho_{X}$,
$\rho_{H}$, and $\rho_{D}$, respectively.

Here, the crystalline (X) state is a global energy minimum, while the hexatic
(H) and liquid (D) states are local energy minima. The model assumes that one
state can transition to another due to thermal fluctuations. Taken together,
we model the dynamics as activated hops between energy minima (Fig. 5(a)),
where the associated energy barriers are given by $\Delta\epsilon_{X}$,
$\Delta\epsilon_{H}$, and $\Delta\epsilon_{D}$, respectively. Since
dislocations and disclinations are excitations of the crystalline state, we
write $\Delta\epsilon_{H}=c_{H}\Delta\epsilon_{X}$ and
$\Delta\epsilon_{D}=c_{D}\Delta\epsilon_{X}$ with $1>c_{H}>c_{D}$. This is
also well supported by previous work Bi  _et al._ (2014) which computed the
energy barriers between four-cell clusters with different topologies. Here,
rather than computing the precise values of the energy barriers, we make an
arbitrary choice of $c_{H}=0.8$ and $c_{D}=0.2$ without the loss of
generality.

###  D.1 In the absence of cell divisions, the meanfield model recapitulates
the two-step melting process

In the absence of cell divisions, the cell motility provides a source of
random fluctuations. In the limit of $D_{r}\gg 1$, it essentially provides a
source of uncorrelated fluctuations that can be described by an effective
temperature given by Bi  _et al._ (2016)

$T\propto v_{0}^{2},$ (2)

which provides the thermal-like activation to overcome barriers in the mean-
field model. The transition rates between states are determined by the
thermally activated process, i.e. $Re^{-\Delta\epsilon/T}$. Here the attempt
frequency $R$ between two states (assumed to be the same for all states), and
the related energy barriers as illustrated in Fig. 5(a). The states therefore
evolve according to

$\displaystyle\dot{\rho}_{X}$
$\displaystyle=Re^{-\Delta\epsilon_{H}/T}\rho_{H}-Re^{-\Delta\epsilon_{X}/T}\rho_{X},$
(3) $\displaystyle\dot{\rho}_{H}$
$\displaystyle=Re^{-\Delta\epsilon_{X}/T}\rho_{X}+Re^{-\Delta\epsilon_{D}/T}\rho_{D}-2Re^{-\Delta\epsilon_{H}/T}\rho_{H},$
$\displaystyle\dot{\rho}_{D}$
$\displaystyle=Re^{-\Delta\epsilon_{H}/T}\rho_{H}-Re^{-\Delta\epsilon_{D}/T}\rho_{D}.$

The steady-state solution gives $\rho_{X}=Q^{-1}e^{\Delta\epsilon_{X}/T}$,
$\rho_{H}=Q^{-1}e^{\Delta\epsilon_{H}/T}$, and
$\rho_{D}=Q^{-1}e^{\Delta\epsilon_{D}/T}$, where the normalization factor is
given by
$Q=e^{\Delta\epsilon_{X}/T}+e^{\Delta\epsilon_{H}/T}+e^{\Delta\epsilon_{D}/T}$.
Choosing $1/R$ as the unit of time and $\Delta\epsilon_{X}$ as the unit of
energy, we plot the fraction of the three states as a function of the
effective temperature in Fig. 13(a). As the temperature increases, the model
exhibits a two-step melting process. The system first transitions from a
crystalline phase to a hexatic phase and then from hexatic to liquid. We can
define each phase by selecting a threshold value as follows: (1) Crystalline:
For values where $\rho_{X}>0.01$ and both $\rho_{H}$ and $\rho_{D}$ are less
than $0.01$. (2) Hexatic: When $\rho_{H}>0.01$ and $\rho_{D}<0.01$. (3)
Liquid: In cases where $\rho_{D}>0.01$. It is important to note that the
results are independent of the precise threshold value chosen.

Figure 13:  (a) Plot of $\rho_{X}$, $\rho_{H}$ and $\rho_{D}$ as a function of
temperature in the absence of cell division as obtained by the meanfield
model. Here, the tissue undergoes a two-step melting process, transitioning
from a crystalline structure to a liquid state through an intermediate hexatic
phase. (b) In the presence of cell divisions ($\gamma_{0}/R=10^{-5}$), the
crystalline state ceases to exist across the entire temperature range, giving
rise to a re-entrant Liquid - Hexatic - Liquid transition.

###  D.2 Cell division/apoptosis alters the melting process

Cell division and apoptosis introduce a new type of active force. The
resulting active cell-shape deformations enhance the fluctuations over energy
barriers by ‘tilting’ the energy landscape. The transition rates are modified
by adding a $\gamma_{0}$ term as shown in Fig. 5(a).The master equations then
become

$\displaystyle\dot{\rho}_{X}=$
$\displaystyle(Re^{-\Delta\epsilon_{H}/T}+\gamma_{0})\rho_{H}-(Re^{-\Delta\epsilon_{X}/T}+\gamma_{0})\rho_{X}$
(4) $\displaystyle\dot{\rho}_{H}=$
$\displaystyle(Re^{-\Delta\epsilon_{X}/T}+\gamma_{0})\rho_{X}+(Re^{-\Delta\epsilon_{D}/T}+\gamma_{0})\rho_{D}-2(Re^{-\Delta\epsilon_{H}/T}+\gamma_{0})\rho_{H}$
$\displaystyle\dot{\rho}_{D}=$
$\displaystyle(Re^{-\Delta\epsilon_{H}/T}+\gamma_{0})\rho_{H}-(Re^{-\Delta\epsilon_{D}/T}+\gamma_{0})\rho_{D}$

The steady-state solution gives
$\rho_{X}=Q^{-1}(e^{-\Delta\epsilon_{X}/T}+\gamma_{0})^{-1}$,
$\rho_{H}=Q^{-1}(e^{-\Delta\epsilon_{H}/T}+\gamma_{0})^{-1}$ and
$\rho_{D}=Q^{-1}(e^{-\Delta\epsilon_{D}/T}+\gamma_{0})^{-1}$ with
$Q=(e^{-\Delta\epsilon_{X}/T}+\gamma_{0})^{-1}+(e^{-\Delta\epsilon_{H}/T}+\gamma_{0})^{-1}+(e^{-\Delta\epsilon_{D}/T}+\gamma_{0})^{-1}$.
Fig. 13(b) shows $\rho_{X}$, $\rho_{H}$ and $\rho_{D}$ as a function of
temperature at a constant division rate ($\gamma_{0}/R=10^{-5}$). In the limit
of high temperature, the model behaves similarly to the case of no division.
However, the dominance of the division term in the transition rates
dramatically changes the steady-state solution in the low-temperature limit,
driving the tissue from a crystal phase to a liquid phase.

In Figure 5(b) of the main text, we present the phase diagram, plotted in
terms of the dimensionless variables $\gamma_{0}/R$ and
$\sqrt{T/\Delta\epsilon_{X}}$. This diagram is constructed using a uniform
threshold criterion across the three distinct phases. This illustrates that
mean-field theory effectively replicates the behavior evident in our numerical
simulations, specifically within the experimentally accessible range of cell
division rates($\gamma_{0}/R>10^{-5}$).

Figure 14:  The behavior of the three state densities as functions of
temperature for a constant, ultra-low choice of division rate
$\gamma_{0}/R=10^{-15}$.  Figure 15:  The state densities (a) $\rho_{X}$ (b)
$\rho_{H}$ (c) $\rho_{D}$ as the function of dimensionless temperature
$T/\Delta\epsilon_{X}$ under the $\gamma_{0}\rightarrow 0$ limit. As the
division decreases to ultra-low values, the state densities share the same
steady solution with the no-division case in a bigger range with a lower limit
of $T/\Delta\epsilon_{X}$. However, the tissue is inevitably in a liquid phase
at the zero temperature limit, having equal $\rho_{X}$, $\rho_{H}$, and
$\rho_{D}$. The possible crystal or hexatic phase can only happen in the
middle-temperature range. The dimensionless division rate
$10^{-4},10^{-6},...10^{-30}$ are represented by the rainbow color from red to
purple. (d) A re-plot of the phase diagram of the MF model where the cell
division rates are displayed on a “log of log" scale to emphasize the behavior
at extremely slow values.

###  D.3 Behavior at ultra-low division rates

Intriguingly, the theory also makes a prediction for ultra-low division rates
($\gamma_{0}/R\ll 10^{-10}$). Here, tissues would experience a notably complex
transition sequence as a function of temperature, following a Liquid-Hexatic-
Crystal-Hexatic-Liquid path at constant $\gamma_{0}/R$ as shown in Fig. 5(b).
We also show one such transition path for $\gamma_{0}/R=10^{-15}$ in Fig. 14.

In order to delve deeper into the behavior as $\gamma_{0}/R\to 0$, we examine
the state densities as functions of temperature for $\gamma_{0}/R$ values
smaller than $10^{-10}$, as illustrated in Fig. 15(a-c). Firstly, in the
ultra-low regime of $\gamma_{0}/R$, there exists a condition for a crystalline
phase at sufficiently low temperatures. Under this condition, the transition
rate from the Crystal phase to the Hexatic phase exceeds that from the Hexatic
phase to the Crystal phase. Simultaneously, the transition rate from the
Dislocation phase to the Hexatic phase is also higher than the rate from the
Hexatic phase to the Dislocation phase. This condition is given by

$\frac{e^{-\Delta\epsilon_{H}/T}+\gamma_{0}/R}{e^{-\Delta\epsilon_{X}/T}+\gamma_{0}/R}<1\quad\text{and}\quad\frac{e^{-\Delta\epsilon_{H}/T}+\gamma_{0}/R}{e^{-\Delta\epsilon_{X}/T}+\gamma_{0}/R}<1.$

As the ratio $\gamma_{0}/R$ decreases further, the stable crystalline phase
not only emerges at specific temperatures (as shown in Fig. 14), but also
expands its temperature range until it ultimately converges with the
crystalline phase present at $\gamma_{0}/R=0$.

Finally, we summarize these behaviors in the phase diagram (Fig. 15(d)), which
is the same as Fig. 5(b), except the $\gamma_{0}/R$ here is plotted on a "log-
of-log" scale to emphasizes the slow convergence of the quantity. The phase
diagram of the MF model exhibits a striking resemblance to that of the 2D
melting on a substrate with quenched random disorder in its topography, as
discussed in work by Sachdev and Nelson Sachdev and Nelson (1984). In both
cases, temperature acts as the driving force for the phase transition.
Interestingly, in our model, cell division serves a role analogous to the
degree of quenched disorder in Sachdev and Nelson (1984). This correlation is
logical since cell division introduces spatially random distortions in the
tissue, akin to the disorder in the substrate, and these effects persist
without annealing over time, forming a quenched random variable.