# A new direct detection electron scattering experiment to search for the X17
particle.

D. Dutta  Mississippi State University  H. Gao  Duke University  A. Gasparian
North Carolina A&T State University  T. J. Hague  North Carolina A&T State
University  Lawrence Berkeley National Laboratory  N. Liyanage  University of
Virginia  R. Paremuzyan  Thomas Jefferson National Accelerator Facility  C.
Peng  Argonne National Laboratory  W. Xiong  Shandong University  P. Achenbach
Thomas Jefferson National Accelerator Facility  A. Ahmidouch  North Carolina
A&T State University  S. Ali  University of Virginia  H. Avakian  Thomas
Jefferson National Accelerator Facility  C. Ayerbe-Gayoso  Mississippi State
University  X. Bai  University of Virginia  M. Battaglieri  Isituto Nazionale
di Fisica Nucleare, Sezione di Genova  H. Bhatt  Mississippi State University
A. Bianconi  Dipartimento di Ingegneria dell’Informazione, Università di
Brescia  Isituto Nazionale di Fisica Nucleare, sezione di Pavia  J. Boyd
University of Virginia  D. Byer  Duke University  P. L. Cole  Lamar University
G. Costantini  Dipartimento di Ingegneria dell’Informazione, Università di
Brescia  Isituto Nazionale di Fisica Nucleare, sezione di Pavia  S. Davis
North Carolina A&T State University  M. De Napoli  Isituto Nazionale di Fisica
Nucleare, Sezione di Catania  R. De Vita  Isituto Nazionale di Fisica
Nucleare, Sezione di Genova  B. Devkota  Mississippi State University  B.
Dharmasena  University of Virginia  J. Dunne  Mississippi State University  L.
El Fassi  Mississippi State University  V. Gamage  University of Virginia  L.
Gan  University of North Carolina, Wilmington  K. Gnanvo  Thomas Jefferson
National Accelerator Facility  G. Gosta  Dipartimento di Ingegneria
dell’Informazione, Università di Brescia  Isituto Nazionale di Fisica
Nucleare, sezione di Pavia  D. Higinbotham  Thomas Jefferson National
Accelerator Facility  C. Howell  Duke University  S. Jeffas  University of
Virginia  S. Jian  University of Virginia  A. Karki  Mississippi State
University  B. Karki  Duke University  V. Khachatryan  Duke University  M.
Khandaker  Energy Systems  V. Kubarovsky  Thomas Jefferson National
Accelerator Facility  I. Larin  University of Massachusetts  M. Leali
Dipartimento di Ingegneria dell’Informazione, Università di Brescia  Isituto
Nazionale di Fisica Nucleare, sezione di Pavia  V. Mascagna  DiSAT, Università
dell’Insubria  Isituto Nazionale di Fisica Nucleare, sezione di Pavia  G.
Matousek  Duke University  S. Migliorati  Dipartimento di Ingegneria
dell’Informazione, Università di Brescia  Isituto Nazionale di Fisica
Nucleare, sezione di Pavia  R. Miskimen  University of Massachusetts  P.
Mohanmurthy  Mississippi State University  H. Nguyen  University of Virginia
E. Pasyuk  Thomas Jefferson National Accelerator Facility  A. Rathnayake
University of Virginia  J. Rittenhouse West  Lawrence Berkeley National
Laboratory  A. Shahinyan  Yerevan Physics Institute  A. Smith  Duke University
S. Stepanyan  Thomas Jefferson National Accelerator Facility  E. van
Nieuwenhuizen  Duke University  L. Venturelli  Dipartimento di Ingegneria
dell’Informazione, Università di Brescia  Isituto Nazionale di Fisica
Nucleare, sezione di Pavia  B. Yu  Duke University  Z. Zhao  Duke University
J. Zhou  Duke University  Lawrence Berkeley National Laboratory

(January 2023)

###### Abstract

A new electron scattering experiment (E12-21-003) [1] to verify and understand
the nature of hidden sector particles, with particular emphasis on the so-
called X17 particle, has been approved at Jefferson Lab. The search for these
particles is motivated by new hidden sector models introduced to account for a
variety of experimental and observational puzzles: excess in $e^{+}e^{-}$
pairs observed in multiple nuclear transitions, the 4.2$\sigma$ disagreement
between experiments and the standard model prediction for the muon anomalous
magnetic moment, and the small-scale structure puzzle in cosmological
simulations. The aforementioned X17 particle has been hypothesized to account
for the excess in $e^{+}e^{-}$ pairs observed from the 8Be M1, 4He M0, and,
most recently, 12C E1 nuclear transitions to their ground states observed by
the ATOMKI group. This experiment will use a high resolution electromagnetic
calorimeter to search for or set new limits on the production rate of the X17
and other hidden sector particles in the $3-60$ MeV mass range via their
$e^{+}e^{-}$ decay (or $\gamma\gamma$ decay with limited tracking). In these
models, the $1-100$ MeV mass range is particularly well-motivated and the
lower part of this range still remains unexplored. Our proposed direct
detection experiment will use a magnetic-spectrometer-free setup (the PRad
apparatus) to detect all three final state particles in the visible decay of a
hidden sector particle allowing for an effective control of the background and
will cover the proposed mass range in a single setting. The use of the well-
demonstrated PRad setup allows for an essentially ready-to-run and uniquely
cost-effective search for hidden sector particles in the $3-60$ MeV mass range
with a sensitivity of 8.9$\times$10-8 \- 5.8$\times$10-9 to $\epsilon^{2}$,
the square of the kinetic mixing interaction constant between hidden and
visible sectors.

##  I Introduction

On large distance scales, the structure of the Universe as inferred from
cosmological data is consistent with dark matter particles that are cold,
collision-less, and interact with ordinary matter purely via gravity [2].
Decades of cosmological data have converged on the standard model of
cosmology, dubbed $\Lambda\rm CDM$, with cold dark matter (CDM) as a crucial
ingredient [3] along with dark energy (the cosmological constant of general
relativity, $\Lambda$) and ordinary matter. Weakly interacting massive
particles (WIMPs) have long been a primary candidate to explain dark matter
(DM) [4, 5]. However, recent searches have ruled out a large segment of the
phase space for WIMPs [6, 7] leading to new theoretical developments [8].
Additionally, several anomalous results in recent years have pointed to holes
in our understanding of the nature of matter, discussed further in Sec. II.

These null results and the anomalous results have led to an increased interest
in dark mediator models. In typical models these mediator particles ($X$)
interact with the standard model (SM) through a ‘kinetic mixing’ mechanism
creating a portal between the hidden sector and the SM. In most models the $X$
particle couples to the electric charges of SM matter. Some newer models
propose that the coupling to SM is tied to flavor of quarks or leptons and not
via the electric charge. A common thread between many of the observed
anomalies is that they can be explained by a new MeV-scale dark mediator
particle. Given the convergence of such models on the $1-100$ MeV mass range
for a new particle, it is critical to focus experimental searches to this mass
range in order to resolve the anomalies.

The availability of a high duty factor, high luminosity electron beam at
Jefferson Lab provides an ideal setup to search for MeV-scale dark mediators
with small coupling constants. The well tested PRad setup in Hall B [9] will
be used in this experiment to reach our physics goals. Using a magnetic-
spectrometer-free setup allows the experiment to be sensitive to the full mass
range in a single experimental setting. The detection of all three final state
particles in the PbWO4 calorimeter along with tracking with GEM chambers
allows for an effective control of the backgrounds. Moreover, it provides an
essentially ready-to-run and uniquely cost effective search for hidden-sector
particles in the 3 - 60 MeV mass range.

##  II Physics Motivation

The primary motivation for this new experiment is to resolve the X17 anomaly.
A 2015 experiment [10] at the ATOMKI 5 MV Van de Graaff accelerator sought to
study a previous 9 MeV anomaly [11, 12, 13]. While this experiment nullified
and explained the previous result, it saw an excess of $e^{+}e^{-}$ pairs
beyond the expectation of internal pair creation (IPC) for a different mass. A
subsequent analysis of these results has shown that the $6.8\sigma$ anomaly is
consistent with a new particle with a mass of 16.84 MeV, dubbed X17. A
followup experiment by the ATOMKI group was conducted on the 20.01 MeV
$0^{-}\to 0^{+}$ transition in 4He. The results for 4He also reports an
$e^{+}e^{-}$ excess consistent with the so-called X17 particle [14, 15]. In a
new 2022 experiment, the group has reported the anomaly in the 12C E1
transition [16]. Of note, this most recent result agrees with the previous
masses but disagrees with the measured branching ratios.

Figure 1: Current constraints on a fifth force explanation of the 8Be anomaly.
The vertical axis is the leptonic coupling strength relative to
$\alpha_{QED}$, with horizontal axis the mass of the mediator. Excluded
regions, in gray, are taken from measurements that depend solely on leptonic
interactions. Dark photon exclusions via hadronic measurements are not shown.
Reproduced from Ref. [17].

Feng et al. [18] analyzed the X17 signal against existing constraints. The
proposed explanation is that the signal is from the decay of a protophobic
gauge boson that mediates a fifth force with a length scale of 12 fm. This
explanation can also possibly explain the muon anomalous magnetic moment and
an excess of $\pi^{0}\to e^{+}e^{-}$ decays [19]. The current constraints from
leptonic production mechanisms, where the effective coupling to a new force-
carrier is proportional to electric charge, are shown in Fig. 1. For more
generic fifth-force models with quark flavor-dependent couplings [18], a much
wider parameter space with multiple couplings must be considered.

The fifth force based explanations are being challenged by recent reanalyses
[20, 21] and the observed discrepancies could be the result of as-yet-
unidentified nuclear reactions, excited diquark states, or experimental
effects. Nonetheless, these results have garnered a lot of attention, and must
be independently validated with the highest urgency.

The Muon $g-2$ collaboration has recently reported their measurement of the
muon magnetic moment at Fermi National Accelerator Laboratory (FNAL) which is
consistent with their previous measurement at Brookhaven National Laboratory
(BNL) [22, 23]. These two results, with their uncertainties combined, show a
$4.2\sigma$ deviation from the Standard Model prediction.

There have been several proposals as to new physics that could lead to this
deviation. One such class of explanations shows that undiscovered hidden
sector particles can directly contribute to loop corrections in the
calculation of $(g-2)_{\mu}$. The contributions of a hidden sector particle to
the $(g-2)_{\mu}$ scales as $1/m_{X}^{2}$, thus a smaller mass will have a
comparably larger contribution than a higher mass particle. There have been
numerous publications proposing new particles (such as dark photons) in the
MeV-range that, at a minimum, partially account for the muon anomalous
magnetic moment [24, 25, 26, 27]. Fig. 1 shows that there is an overlap in the
phase space for a particle that would simultaneously resolve the X17 and
$(g-2)_{\mu}$ anomalies.

Several recent models have suggested that the coupling of the X17 to SM is
tied to flavor of quarks or leptons and not via the electric charge. One such
example is the protophobic X17 proposed in Ref. [18] used to explain the X17
anomaly. If an X17 were to have flavor-dependent couplings the parameter space
for it is more open than previously thought [28]. These new possibilities
point to an urgent need for independent verification of the X17 anomaly.

##  III Experimental Method and Setup

This experiment will focus on the bremsstrahlung-like production of hidden
sector particles from the initial electron or the scattered electron (both
shown in Fig.2), in the $3-60$ MeV mass range. When searching for new
particles, it is of the utmost importance to minimize backgrounds in order to
prevent false “bumps” in the mass spectrum. The primary QED background in this
experiment (see Fig. 3) are from the radiative pair production which is an
irreducible background and the Bethe-Heitler trident reactions which can be
kinematically suppressed.

Figure 2: Bremsstrahlung-like production of a hidden sector force carrier $X$
from electron scattering. The left diagram shows production from the incoming
electron and the right diagram shows production from the scattered electron.

Figure 3: The QED background from radiative (left and middle) and Bethe-
Heitler (right) process.

The experiment plans to reuse the PRad setup (shown in Fig. 4), with the
addition of a Tantalum foil target placed 7.5 m upstream of the calorimeter.
The high resolution PbWO4 crystal part of the HyCal electromagnetic
calorimeter [29, 30] will be used together with a new flash ADC (FADC) based
readout system for the calorimeter. A 2.2 and 3.3 GeV Continuous Wave (CW)
electron beam from CEBAF will be incident on a retractable ultra-thin target
consisting of a 1 $\mu$m Tantalum foil. The scattered particles will traverse
the 7.5 m long flight path in vacuum, en route to a pair of common ionization
volume GEM chambers coupled to the HyCal EM calorimeter. The vacuum chamber
will consist of the PRad target chamber (or an appropriate diameter beam pipe)
and the 5 m long PRad vacuum chamber with a thin window. The vacuum chamber
matches the geometrical acceptance of the calorimeter. An extension piece will
be added to the upstream end to couple the PRad target chamber to the super-
harp which will now hold the target foil. A reducer ring will be attached to
the downstream exit of the PRad vacuum chamber and a new 1 mm thick Al exit
window will be used such that it matches the PbWO4 portion of the calorimeter.
Two layers of GEM detectors will add a modest tracking capability to help
reduce the photon background and to reduce the background originating from the
vacuum chamber exit window.

Figure 4: Schematic of the experimental setup.

The experiment will use a “bump hunt” technique in the direct detection search
of heretofore unknown MeV mass particles. All 3 cluster events with individual
cluster energy within (0.02-0.85)$\times E_{beam}$ and with the sum of total
energy deposited $E_{sum}>$ 0.7$\times E_{beam}$ will be recorded and examined
for “bumps” in the $M_{e^{+}e^{-}}$ (or $M_{\gamma\gamma}$) invariant mass
spectrum reconstructed from these events. This will allow for the
$X$-particles production by virtual photons over a wide energy range in the
forward solid angle coverage of the PbWO4 calorimeter. The capability of
detecting events produced in a wide energy and angle range, in a single
experimental setting, is one of the important features of our experiment which
stands in contrast to all other magnetic spectrometer methods. The tracking
provided by the pair of GEM chambers will be used to suppress background
events from the experimental setup, particularly, the large area window at the
exit of the vacuum chamber. The GEM chambers will also be used to veto any
neutral particles. Only the PbWO4 part of the HyCal calorimeter will be used
in the experiment to ensure high invariant mass resolution. The experimental
method discussed here applies directly to any spin-1 boson in the hidden
sector that decay directly to a lepton pair.

##  IV Projected Results

Figure 5: Acceptance of the spectrometer setup for both 2.2 and 3.3 GeV beam
energies. Figure 6: Invariant mass resolution of the spectrometer for 2.2 GeV
beam energy.

Figure 7: Projected coverage of the square of the coupling constant
($\epsilon^{2}$) and mass ($m_{X}$) parameter space by this proposed
experiment is shown as the thick red lines for the combined statistics of the
two beam energies. The projections are superimposed on top of a constraints
plot adapted from Ref. [8], with additional data from Refs.[31, 32, 33].

In order to estimate the projected sensitivity to $\epsilon^{2}$, a
comprehensive simulation of the experiment was carried out using the Geant4
simulation package developed for the PRad experiment [9]. The simulation used
a hybrid event generator that was built by combining the Madgraph5 generator
[34] for the Bethe-Heitler process with the Geant4 generator (with Bethe-
Heitler turned off) to utilize the best features of the two generators. The
background was simulated for 200 s of beamtime and then scaled by sampling the
number of events bin-by-bin to give the expected background for 30 days of 100
nA, 3.3 GeV beam.

The production rates for the $X$ particle are calculated using the rate
equation (Eq. 14) from Ref. [35] and rearranged to give the sentitivity as:

$\epsilon^{2}=\frac{N_{X}}{{\mbox{5}}\times
N_{e}T\frac{m_{e}^{2}}{m_{X}^{2}}},$ (1)

where $N_{X}$ is the number of $X$ particles produced, $N_{e}$ is the number
of incident electrons (1.62$\times$1018 for 30 days at 100 nA), $T$ is the
target thickness (2.5$\times$10-4 r.l.), $m_{X}$ is the mass of the produced
$X$, and $\epsilon^{2}$ is the square of the dimensionless coupling constant
of the $X$ to SM matter.

The lowest $\epsilon^{2}$ achievable is calculated using Eq. 1 with $N_{X}$
obtained by scaling the signal counts with the detection efficiency and impact
of the decay length (shown in Fig. 5). The background event count is
calculated by selecting events from the background simulation within a $\pm$
3$\sigma_{m_{X}}$ window. The typical invariant mass resolution
$\sigma_{m_{X}}$ obtained from the simulation is shown in Fig. 6. The signal
counts necessary given the background counts are calculated using the criteria
for a 2.3$\sigma$ significance used by other experiments (Eq. 2).

$\frac{N_{\mbox{signal}}}{\sqrt{N_{\mbox{signal}}+N_{\mbox{bgd}}}}\geq{\mbox{2.3}}.$
(2)

The 2.2 GeV runs will serve as reference runs to better understand the
background and the signal, so it was conservatively estimated that 50% of
their statistics can be combined into the 3.3 GeV data for the final results.
Using these range of sensitivities the bounds for the $\epsilon^{2}-m_{X}$
parameter space is plotted in Fig. 7 for the combined projected statistics.

##  V Conclusion

The X17 anomaly suggests the existence of a new hidden sector particle. It is
of critical importance to resolve this anomaly quickly, as it not only could
point to new physics but also has the potential to resolve additional
anomalies. In this white paper, we have described a nearly ready-to-run
experiment at JLab using the PRad spectrometer that has the capability to
fully address this anomaly. By using a magnet-free spectrometer, two beam
energies, and detecting all three final state the experiment will provide
unprecedented control of the experimental background, leading to a results
with uniquely low systematics.

###### Acknowledgements.

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Nuclear Physics under contract DE-AC05-06OR23177.

## References

  * Ahmidouch  _et al._ [2021] A. Ahmidouch _et al._ , A Direct Detection Search for Hidden Sector New Particles in the 3-60 MeV Mass Range, (2021), arXiv:2108.13276 [nucl-ex] .
  * S. and H.-B. [2018] T. S. and Y. H.-B., Dark Matter Self-interactions and Small Scale Structure, Phys. Rep. 730, 1 (2018), arXiv:1705.02358 [hep-ph] .
  * A.  _et al._ [1999] B. N. A., O. J. P., P. S., and S. P. J., The Cosmic Triangle: Revealing the State of the Universe, Science 284, 1481 (1999), arXiv:9906463 [astro-ph] .
  * Jungman  _et al._ [1996] G. Jungman, M. Kamionkowski, and K. Griest, Supersymmetric dark matter, Physics Reports 267, 195 (1996).
  * Feng [2010] J. L. Feng, Dark Matter Candidates from Particle Physics and Methods of Detection, Ann. Rev. Astron. Astrophys. 48, 495 (2010), arXiv:1003.0904 [astro-ph.CO] .
  * Battaglieri  _et al._ [2017] M. Battaglieri _et al._ , US Cosmic Visions: New Ideas in Dark Matter 2017: Community Report,  _U.S. Cosmic Visions: New Ideas in Dark Matter_ , arXiv:1707.04591 (2017), arXiv:1707.04591 [hep-ph] .
  * Aprile  _et al._ [2018] E. Aprile _et al._ , Dark Matter Search Results from a One Ton-Year Exposure of XENON1T, Phys. Rev. Lett. 121, 111302 (2018).
  * Alexander  _et al._ [2016] J. Alexander _et al._ , Dark Sectors 2016 Workshop: Community Report, (2016), arXiv:1608.08632 [hep-ph] .
  * Xiong  _et al._ [2019] W. Xiong _et al._ , A small proton charge radius from an electron–proton scattering experiment, Nature 575, 147 (2019).
  * Krasznahorkay  _et al._ [2016] A. J. Krasznahorkay _et al._ , Observation of Anomalous Internal Pair Creation in Be8 : A Possible Indication of a Light, Neutral Boson, Phys. Rev. Lett. 116, 042501 (2016), arXiv:1504.01527 [nucl-ex] .
  * de Boer  _et al._ [1996] F. de Boer, O. Fröhlich, K. Stiebing, K. Bethge, H. Bokemeyer, A. Balanda, A. Buda, R. van Dantzig, T. Elze, H. Folger, J. van Klinken, K. Müller, K. Stelzer, P. Thee, and M. Waldschmidt, A deviation in internal pair conversion, Phys. Lett. B 388, 235 (1996).
  * de Boer  _et al._ [1997] F. W. N. de Boer, R. van Dantzig, J. van Klinken, K. Bethge, H. Bokemeyer, A. Buda, K. A. Müller, and K. E. Stiebing, Excess in nuclear $e^{+}e^{-}$ parin near 9 MeV/c2 invariant mass, J. Phys. G 23, L85 (1997).
  * de Boer  _et al._ [2001] F. W. N. de Boer, K. Bethge, H. Bokemeyer, R. van Dantzig, J. van Klinken, V. Mironov, K. A. Muller, and K. E. Stiebing, Further search for a neutral boson with a mass around 9-MeV/c2, J. Phys. G 27, L29 (2001), arXiv:hep-ph/0101298 .
  * Krasznahorkay  _et al._ [2019] A. J. Krasznahorkay _et al._ , New evidence supporting the existence of the hypothetic X17 particle, (2019), arXiv:1910.10459 [nucl-ex] .
  * Krasznahorkay  _et al._ [2021] A. J. Krasznahorkay, M. Csatlós, L. Csige, J. Gulyás, A. Krasznahorkay, B. M. Nyakó, I. Rajta, J. Timár, I. Vajda, and N. J. Sas, New anomaly observed in He4 supports the existence of the hypothetical X17 particle, Phys. Rev. C 104, 044003 (2021), arXiv:2104.10075 [nucl-ex] .
  * Krasznahorkay  _et al._ [2022] A. J. Krasznahorkay _et al._ , New anomaly observed in 12C supports the existence and the vector character of the hypothetical X17 boson, Phys. Rev. C 106, L061601 (2022), arXiv:2209.10795 [nucl-ex] .
  * dar [2020] Search for new physics in e+e-final states near an invariant mass of 17 mev using the cebaf injector, http://www.jlab.org/exp_prog/proposals/20/PR12-20-001_Proposal.pdf (2020).
  * Feng  _et al._ [2016] J. L. Feng, B. Fornal, I. Galon, S. Gardner, T. M. P. Smolinsky, J.and Tait, and P. Tanedo, Protophobic Fifth-Force Interpretation of the Observed Anomaly in 8Be Nuclear Transitions, Phys. Rev. Lett. 117, 071803 (2016), arXiv:1604.07411 [hep-ph] .
  * Abouzaid  _et al._ [2007] E. Abouzaid _et al._ (KTeV), Measurement of the Rare Decay $\pi^{0}\to e^{+}e^{-}$, Phys. Rev. D 75, 012004 (2007), arXiv:hep-ex/0610072 .
  * Zhang and Miller [2021] X. Zhang and G. A. Miller, Can a protophobic vector boson explain the ATOMKI anomaly?, Phys. Lett. B 813, 136061 (2021), arXiv:2008.11288 [hep-ph] .
  * Kubarovsky  _et al._ [2022] V. Kubarovsky, J. R. West, and S. J. Brodsky, Quantum Chromodynamics Resolution of the ATOMKI Anomaly in ${\rm{{}^{4}He}}$ Nuclear Transitions, (2022), arXiv:2206.14441 [hep-ph] .
  * Albahri  _et al._ [2021] T. Albahri _et al._ (Muon g-2), Magnetic-field measurement and analysis for the Muon $g-2$ Experiment at Fermilab, Phys. Rev. A 103, 042208 (2021), arXiv:2104.03201 [hep-ex] .
  * Bennett  _et al._ [2006] G. W. Bennett _et al._ (Muon g-2), Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL, Phys. Rev. D 73, 072003 (2006), arXiv:hep-ex/0602035 .
  * Arcadi  _et al._ [2021] G. Arcadi, A. S. De Jesus, T. B. De Melo, F. S. Queiroz, and Y. S. Villamizar, A 2HDM for the g-2 and Dark Matter, arXiv:2104.04456 (2021), arXiv:2104.04456 [hep-ph] .
  * Athron  _et al._ [2021] P. Athron, C. Balázs, D. H. Jacob, D. Kotlarski, W.and Stöckinger, and H. Stöckinger-Kim, New physics explanations of $a_{\mu}$ in light of the FNAL muon $g-2$ measurement, arXiv:2104.03691 (2021), arXiv:2104.03691 [hep-ph] .
  * Borah  _et al._ [2021] D. Borah, M. Dutta, S. Mahapatra, and N. Sahu, Muon $(g-2)$ and XENON1T Excess with Boosted Dark Matter in $L_{\mu}-L_{\tau}$ Model, arXiv:2104.05656 (2021), arXiv:2104.05656 [hep-ph] .
  * Ge  _et al._ [2021] S.-F. Ge, X.-D. Ma, and P. Pasquini, Probing the Dark Axion Portal with Muon Anomalous Magnetic Moment, arXiv:2104.03276 (2021), arXiv:2104.03276 [hep-ph] .
  * Feng  _et al._ [2022] J. L. Feng _et al._ , The Forward Physics Facility at the High-Luminosity LHC, (2022), arXiv:2203.05090 [hep-ex] .
  * Gasparian [2002] A. Gasparian (PrimEx), Performance of PWO crystal detectors for a high resolution hybrid electromagnetic calorimeter at Jefferson Lab, in  _10th International Conference on Calorimetry in High Energy Physics (CALOR 2002)_ (2002) pp. 208–214.
  * Gasparian [2004] A. Gasparian (PrimEx), A high performance hybrid electromagnetic calorimeter at Jefferson Lab, in  _11th International Conference on Calorimetry in High-Energy Physics (Calor 2004)_ (2004) pp. 109–115.
  * Aaij  _et al._ [2018] R. Aaij _et al._ (LHCb), Search for Dark Photons Produced in 13 TeV $pp$ Collisions, Phys. Rev. Lett. 120, 061801 (2018), arXiv:1710.02867 [hep-ex] .
  * Aaij  _et al._ [2020] R. Aaij _et al._ (LHCb), Search for $A^{\prime}\to\mu^{+}\mu^{-}$ Decays, Phys. Rev. Lett. 124, 041801 (2020), arXiv:1910.06926 [hep-ex] .
  * Baltzell  _et al._ [2022] N. Baltzell _et al._ , The Heavy Photon Search Experiment, (2022), arXiv:2203.08324 [hep-ex] .
  * Conte  _et al._ [2013] E. Conte, B. Fuks, and G. Serret, MadAnalysis 5, A User-Friendly Framework for Collider Phenomenology, Comput. Phys. Commun. 184, 222 (2013), arXiv:1206.1599 [hep-ph] .
  * Bjorken  _et al._ [2009] J. D. Bjorken, R. Essig, P. Schuster, and N. Toro, New Fixed-Target Experiments to Search for Dark Gauge Forces, Phys. Rev. D 80, 075018 (2009), arXiv:0906.0580 [hep-ph] .


# Role of tissue fluidization and topological defects in epithelial
tubulogenesis

Richard D.J.G. Ho  Njord Centre, Physics Department, University of Oslo,
Norway  Stig Ove Bøe  Department of Microbiology, Oslo University Hospital,
0372, Oslo, Norway  Dag Kristian Dysthe  Njord Centre, Physics Department,
University of Oslo, Norway  Luiza Angheluta  Njord Centre, Physics Department,
University of Oslo, Norway

###### Abstract

Cellular rearrangements, as primary sources of tissue fluidization, facilitate
topological transitions during tissue morphogenesis. We study the role of
intrinsic cell properties such as cell polarity and cell-cell adhesion in
shaping epithelial tissues using a minimal model of interacting polarized
cells. The presence of a vortex in the cell polarization poses the topological
constraint that induces an inwards migration with the formation of a conical
shape. Local rearrangements at the tip of the cone lead to the onset of tube
formation. Switching between collective migration and structural
rearrangements is key for balancing the contrasting tendencies, such as the
tissue rigidity needed to preserve shape and the tissue fluidity allowing for
topological transitions during tissue morphogenesis.

tissue dynamics, topological defects, morphogenesis

##  I  Introduction

Collective structural arrangements and migration of cells within epithelial
tissues are important physical processes that underlie various phenomena such
as embryonic or tissue development, wound healing, homeostasis and cancer
metastasis. Understanding the complexity of cellular interactions and the
resulting emergent collective behaviors presents formidable challenges both
experimentally and theoretically. Epithelial tissues are a ubiquitous type of
biological tissue found throughout the body of multicellular organisms.
Epithelia are organized as mono- or multi-layered cell sheets that function as
protective layers covering internal organs, line body cavities, and form the
outer skin layer.

On the structural level, epithelial monolayers can exhibit two distinct types
of cell polarity: i) the apical-basal (AB) polarity, which provides the
surface orientation, and ii) the planar cell polarity (PCP), orthogonal to the
AB polarity axis that refers to the structural arrangement of cells within the
tissue surface [1, 2]. Apical-basal polarity represents the main polarity
axis, defining the “top” (apical surface) and “bottom” (basal surface) of the
epithelial layer. Typically, the apical side faces the lumen or the outer
surface of the monolayer, while the basal side is in contact with the basal
membrane, which faces the underlying connective tissue [3]. Planar cell
polarity, on the other hand, is specified by a group of proteins specifically
localized to the “front” and “back” (or “left” and “right”) of each cell,
leading to structural alignment of cells within the monolayer [4]. The
functional role of PCP is evident in fully developed skin tissues, where PCP
helps align hair follicles [5, 6]. In addition, PCP proteins may also play a
role in dynamic processes during tissue development [7, 8, 9, 10, 11].
Epithelial monolayers can also undergo spontaneous polarization and collective
migration in a manner that seem to be independent of PCP proteins. This
process typically involves structural rearrangements of actin filaments,
forming protrusions, called lamellopodia, at the leading edge of the cell.
Such spontaneous transition is fundamental to the epithelial-to-mesenchymal
transition (EMT), which plays a crucial role in wound healing, development,
and tumor progression[12].

Another characteristics of epithelial monolayers is their ability to maintain
tight inter-cellular connections through cell-cell adhesion. This is essential
for their barrier function and for their ability to control transport across
the tissue. However, under specific conditions involving dynamic tissue
remodeling, these epithelial sheets exhibit an ability to modulate cell-cell
adhesion. By temporarily loosening cell-cell adhesion, they facilitate
increased neighbor exchange and rearrangement, enabling necessary tissue
movements. Examples of this kind of epithelial fluidization occur during wound
closure [13] and tissue morphogenesis and development [14, 15, 16].

Epithelial monolayers have the ability to form intricate three-dimensional
morphological structures such as folds, tubes, and branching networks [7, 8].
During morphogenesis, topological defects with full-integer charges emerge as
key organizational centers for morphological events guiding topological
transformation in evolving shapes. This is evidenced by experimental studies
in the model organism Hydra, where long-lived $+1$ topological defects are
formed through the nematic ordering of actin filaments and facilitate
epithelial morphogenesis [17]. In addition, a growing body of empirical
evidence from various tissue systems suggests that similar topological defects
also appear in the velocity patterns or structural alignment of cells within
epithelial tissues during development and homeostasis. One prominent example
is found in the small intestinal epithelium, characterized by the presence of
crypts and villi. During homeostasis, epithelial cells within the crypt
converge outwards and form a $+1$ defect in their structural arrangement.
Similarly, at the villi tips, where cells migrate radially before shedding,
effectively forming $+1$ defects by cells inwardly migrating towards the
defect core [18, 19, 20]. Furthermore, many epithelial tissues, such as those
in the kidney, lung, and glands, exhibit branching and tubular morphologies.
The initiation of these structures is thought to involve the formation of buds
that extend orthogonal from the parent tissue, potentially creating a $+1$
defect within the velocity field at the budding site [21]. In addition,
endothelial tissue, a specialized epithelium lining the vascular system, also
exhibits similar branching morphology, indicating a potential role of $+1$
defects in sprouting and blood vessel generation during angiogenesis [21, 22].
Finally, the role of $+1$ defects in shaping epithelial morphogenesis is
supported by studies showing that epithelial monolayers cultured on flat
surfaces spontaneously form vortices and swirls, signifying flow patterns that
generate $+1$ defects in the velocity field [23, 24, 25, 26, 27].

Despite strong evidence suggesting a role for $+1$ defects in epithelial
morphogenesis, several key questions remain unanswered, in particular what
types of collective behaviors induce the formation of these defects and their
feedback on morphological changes and tissue flows.

In this paper, we study the role of topological constraints and tissue
dynamics on the tissue morphology. We approach this problem using a
theoretical model of a tissue monolayer formed by polarized cells interacting
with each other through cell-cell adhesion forces. We predict that during
morphogenesis, the epithelial tissue evolves towards balancing competing
tendencies, namely tissue rigidity needed for its functional integrity and
sustaining shapes, and tissue fluidity enabled by structural rearrangements
and needed for shape transformations. Within our model, the tissue rigidity is
intrinsically rendered by the collective ordering of planar cell polarities,
while in biological tissues there might also be environmental factors such as
the cytoskeleton which may lead to the same effect. We find that a vortex in
the tissue polarization induces a relaxational dynamics towards a conical
shape. However, since the cone harbors a singular point, it is a source of
structural rearrangements thus local fluidization which triggers the onset of
tube formation. The straight tube is sustained by the concentric pattern in
the planar cell polarities with a fixed vortex. Fluctuations can break this
radial symmetry and induce the formation a flap-like structure instead.

The rest of the paper is structured as follows. In Section II, we introduce a
minimal model of polarized cells with adhesive interactions and discuss the
finite size effects of the disk geometry. The cell migration patterns emerging
around isolated full-integer defects in the orientation of planar cell
polarities within a flat tissue are studied in Section III. We show that
asters and vortices induce transient outward and inward migration,
respectively, whereas spirals sustain persistent spiraling migration. In
Section IV, we discuss the role of localized tissue fluidization in the
transition from a cone to tubular structures. Concluding remarks and a summary
are presented in Section V.

##  II  Minimal model of polarized epithelium

Figure 1:  A) Two neighboring cells, $i$ and $j$, since for all other third
cells, $k$, $r_{ij}<r_{k-k}$, thus $r_{ij}^{2}<r_{ik}^{2}+r_{jk}^{2}$ and the
angle $\angle IKJ<\pi/2$. B) Two cells that are not neighbors, since there
exists a third cell, $k$, where the angle $\angle IKJ>\pi/2$
($\lambda_{0}=2$).

To explore the interplay between tissue polarization and fluidization, we
introduce a minimal model of polarized cells that interact with their
neighbors through adhesion forces modulated by the cell polarities [28, 29].
The Apical-Basal (AB) polarity $\mathbf{p}$ is a unit vector which is along
the up-down axis of a cell and represents the cell elongation out of the
tissue surface. Hence, $\mathbf{p}$ tends to point in the same direction as
the outward normal of the tissue surface. Orthogonal to this is the planar
cell polarity (PCP) represented by the unit vector $\mathbf{q}$ which lies
within the tissue surface. The cells interact with their nearest neighbors
(n.n) through a potential energy that contains two contributions: i) an
isotropic interaction $V_{ij}^{(0)}$ that depends only on the pair separation
distance, $r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|$, and ii) an anisotropic
part $V_{ij}^{(1)}$ that account for soft repulsion/attraction forces
modulated by the local orientation of the cell polarities relative to each
other and to their pair separation vector. Hence,
$V_{ij}=V_{ij}^{(0)}+V_{ij}^{(1)}$, where [28]

$V_{ij}^{(0)}=\lambda_{0}\left(e^{-r_{ij}}-e^{-r_{ij}/a_{0}}\right),$ (1)

with $a_{0}>1$, and minimized by an equilibrium pair distance
$r_{0}=\frac{a_{0}}{a_{0}-1}\ln a_{0}$. Adhesion forces mediated by polarities
are anisotropic and described by the second contribution

$V^{(1)}_{ij}=e^{-r_{ij}}-(\lambda_{1}S_{1}+\lambda_{2}S_{2}+\lambda_{3}S_{3})e^{-r_{ij}/a_{0}},$
(2)

which is minimized by the same equilibrium distance when
$\lambda_{1}+\lambda_{2}+\lambda_{3}=1$. The strength of the attraction is
modulated by the orientation of polarities relative to each other and the pair
separation vector,

$\displaystyle S_{1}$ $\displaystyle=$
$\displaystyle(\mathbf{p}_{i}\times\mathbf{r}_{ij})\cdot(\mathbf{p}_{j}\times\mathbf{r}_{ij})$
(3) $\displaystyle S_{2}$ $\displaystyle=$
$\displaystyle(\mathbf{p}_{i}\times\mathbf{q}_{i})\cdot(\mathbf{p}_{j}\times\mathbf{q}_{j})$
(4) $\displaystyle S_{3}$ $\displaystyle=$
$\displaystyle(\mathbf{q}_{i}\times\mathbf{r}_{ij})\cdot(\mathbf{q}_{j}\times\mathbf{r}_{ij}),$
(5)

such that orthogonal vectors are favoured for equilibrium configurations. We
may notice that $S_{1}$ couples the surface normal with the cell positions,
and maintains cells within the tissue surface whereby the AB polarities align
along the normal to the tissue surface. On the other hand, the $S_{3}$
coupling the cell positions with the planar cell polarities may contribute to
the normal migration of the cells. The $S_{2}$ term is important for
maintaining orthogonality between the AB and PCP polarities, i.e. having two
distinct polarities. The coefficient corresponding to $S_{1}$ needs to be
greater than that of $S_{3}$ in order to identity the AB polarity with the
surface orientation as discussed later.

Figure 2:  A) Metastable orientation of $\mathbf{q}$. B) Envelope of the
deviation of $\mathbf{q}$ orientation from the uniform bulk for different
system sizes. Semi-log inset shows the exponential relaxation with the
characteristic lengthscale fitted to $\xi\approx 0.2R$.

Cells migrate and reorient their polarities due to these interactions,
following the overdamped dynamics

$\displaystyle\dot{\mathbf{r}}_{i}$ $\displaystyle=$
$\displaystyle-\sum_{j=\textrm{n.n. of }i}\frac{\partial
V_{ij}}{\partial\mathbf{r}_{i}}+\bm{\eta}^{(1)}_{i}$ (6)
$\displaystyle\dot{\mathbf{p}}_{i}$ $\displaystyle=$
$\displaystyle-\sum_{j=\textrm{n.n. of }i}\frac{\partial
V_{ij}}{\partial\mathbf{p}_{i}}+\bm{\eta}^{(2)}_{i}$ (7)
$\displaystyle\dot{\mathbf{q}}_{i}$ $\displaystyle=$
$\displaystyle-\sum_{j=\textrm{n.n. of }i}\frac{\partial
V_{ij}}{\partial\mathbf{q}_{i}}+\bm{\eta}^{(3)}_{i},$ (8)

with small random fluctuations $\bm{\eta}^{(\alpha)}$ from the environment,
modelled as uncorrelated Gaussian noise of zero mean and small standard
deviation

$\langle{\eta}_{i}^{(\alpha)}(t){\eta}_{j}^{(\beta)}(t^{\prime})\rangle=2\sigma\delta(t-t^{\prime})\delta_{ij}\delta_{\alpha\beta}.$
(9)

The deterministic evolution preserves the unit norm of the polarities.
However, due to noise, the norms can fluctuate. Therefore, we normalize the
polarity vectors after each timestep.

Figure 3: A) Radial dependence of the average cell density at a late time
$t=1000$ for different types of defects. B) Instantaneous cell migration
profile (blue) for an aster (B) and a vortex (C). Averages are taken over the
transient time, $10$ initial configurations, and noise seeds. The
corresponding polarity profiles are plotted in red. ($\lambda_{0}=2$)

#### Cell rearrangements:

We use an Euclidean metric to construct the network of cell neighbors. Cells
$i$ and $j$ are considered neighbors if and only if they are closest in
distance to their midpoint than any other third cell. This rule is achieved
when the distances between cells $i$ and $j$ and any other cell $k$ satisfy
the cosine rule $r_{ij}^{2}<r_{ik}^{2}+r_{jk}^{2}$, as illustrated in Fig.
(1). This condition is similar to the Voronoi construction, but it is
computationally more efficient and can be generalized to three-dimensions, and
beyond single surfaces.

We distinguish two different dynamical regimes: i) cells remain in contact
with their initial neighbors thus rendering tissue rigidity, and ii) cells can
undergo neighbor exchanges to maintain the minimum pair distance, resulting in
tissue fluidization. Numerically, it suffices to update the cell neighbor list
at regular time intervals which may be longer than the discretization
timestep.

The evolution equations from (6)-(8) with fixed cell-cell connectivity
represent a relaxational dynamics towards minimising the total potential
energy stored by surface bending (curvature) or orientational distortions in
the alignment of planar cell polarities. In this regime, the collective cell
migration tends to get arrested in metastable configurations due to
frustrations in the cell connectivity network. This is, however, circumvented
by relative migration corresponding to local rearrangements to preserve the
minimum neighbor distance.

In this computational study, we set the parameters associated with the
anisotropic interactions to fixed values given by $\lambda_{1}=5/11$,
$\lambda_{2}=4/11$ and $\lambda_{3}=2/11$ such that they sum up to one. The
relative values of $\lambda_{1}$ and $\lambda_{3}$ determine which of the
polarities point out of the surface and which remain in the surface. For
$\lambda_{1}>\lambda_{2}$, we ensure that the $\mathbf{p}$’s represent the AB
polarities pointing out of the surface. We have included the specific value of
$\lambda_{0}$ for the isotropic interactions in each figure caption. The value
of $\lambda_{0}$ determines the strength of the isotropic interactions in
comparison to the polar forces. Since we keep
$\lambda_{1}+\lambda_{2}+\lambda_{3}=1$, a value $\lambda_{0}\ll 1$ means that
isotropic adhesion is negligible, whilst $\lambda_{0}\gg 1$ leads to clumping.
The isotropic interactions are important for maintaining the tissue as a
simply connected surface under outwards migration, like the one induced by an
aster (no holes).

We simulate the dynamical equations using as initial condition a disk geometry
with hexagonal packing of the cells. The cell positions and polarities have
free boundary conditions. For a flat tissue, we fixed the AB polarities along
the z-axis, i.e. $\mathbf{p}=(0,0,1)$. This suffices to maintain the evolution
of the cells within the $(x,y)$ plane. The different configurations of cell
polarities are discussed separately in the next sections.

#### Finite-size effects:

We use free boundary conditions for the dynamical variables. To quantify the
boundary effect on the PCP polarities, we consider a metastable configuration
where the $\mathbf{q}$ polarities uniformly point in the same direction for
the bulk cells and rotate in different directions for the cells within a
boundary layer to ensure that there is a net $+2\pi$ rotation imposed by disk
geometry through its Euler characteristic $\chi=1$. This is also illustrated
in Fig. (2 a)). We measure the local deviation in the orientation of the
$\mathbf{q}$ polarities from their uniform (bulk) orientation
$\mathbf{q}_{b}$, i.e. $1-\mathbf{q}_{b}\cdot\mathbf{q}_{n}$ as a function of
the radial distance from the center of the disk domain. In Fig. (2 B), the
envelope of this deviation is plotted as function of the radial distance $r/R$
(relative to the radius $R$). The deviation decreases exponentially
$e^{-(R-r)/\xi}$ with the distance from the boundary and its corresponding
characteristic length $\xi\approx 0.2R$ is independent of model parameters.
The hexagonal packing of cells introduces a ragged edge of the disc which also
influences the perpendicular orientation of the polarity at the disc edge.
However, this effect is much more localized and negligible compared to the one
introduced by the free boundary conditions.

Figure 4:  A) Transient cell migration induced by $-1$ defect deforms the disc
into an "X" shape. Averages are taken over the transient time, initial
conditions and noise seeds. The structure of the polarities is shown in red.
B) Snapshot of the tissue configuration at $t=5000$ induced by the $-1$
defect. ($\lambda_{0}=2$)

##  III  Migration patterns in a flat tissue

Through cell-cell adhesion forces described by the Eq. (6)-(8), the PCP
polarities $\mathbf{q}_{i}$ tend to align along a preferred orientation to
form polar order. Since the disk geometry has the Euler characteristic
$\chi=1$, there will be a net $+2\pi$ rotation of the $\mathbf{q}_{i}$. This
can be achieved as an edge effect (see Fig. 2) or through a structural
arrangement induced by the presence of a topological defect of $+1$ charge at
the center of the disk. Topological defects locally melt the orientational
order and induce long-range deformations which feed into both collective and
relative migration.

We first consider a flat tissue with fixed AB polarities along the normal to
the tissue plane, i.e. $\mathbf{p}_{n}=(0,0,1)$, and allow cells to migrate
following Eqs. (6) with or without neighbor exchanges.

As proof-of-concept and to bridge with hydrodynamic models, we consider the
setup of a single defect embedded in a uniform polarity orientation. The
profile of the PCP polarities $\mathbf{q}_{i}$ induced by the defect is fixed
instantaneously to the imposed orientational structure of an isolated defect.
Alternatively, one may include an energy contribution due to deviations from
the $\mathbf{p}_{i}$ polarities from the imposed profile, resulting in a
linear restoring force to the desired polarity configuration.

Figure 5:  A-B) Average azimuthal flow velocity as a function of the radial
distance from a spiral ($m=1$, $\theta_{0}=\pi/4$) for different disk radius,
$R$. The dashed line shows the $1/r$ scaling. The finite size effect is
determined by the thickness of the boundary layer $\xi$. C) Persistent flow
for a spiral defect. Motion of the particles in blue, embedded defect in
orange ($\lambda_{0}=2$).

The isolated defect inserted at the center of the disk imposes a polar
orientational field $\theta$ given by [30]

$\theta=\theta_{0}+m\arctan{(y/x)}\ ,\ $ (10)

where $m=\pm 1$ is the topological charge picked up by an arbitrary contour
integral $\oint\limits_{C}d\theta=2\pi m$ enclosing the defect. The constant
phase $\theta_{0}$ gives the baseline of the defect. For a defect with $m=-1$,
the baseline phase can be set to $\theta_{0}=0$ by a reorientation of the
defect. However, this baseline phase $\theta_{0}$ is an intrinsic phase for
the $m=+1$ defect, and its value distinguishes three types of defects:
$\theta_{0}=0$ corresponds to an aster, $\theta_{0}=\pi/2$ for a vortex and
$\theta_{0}\in(0,\pi/2)$ gives a spiral defect [30]. Due to their migrations,
cells acquire different orientations of their PCP polarity depending on their
distance from the defect. Thus, while the polar orientation $\theta$ is
quenched, the actual $\mathbf{q}_{i}$ polarities are obtained by evaluating
$\theta_{i}=\theta(\mathbf{r}_{i})$ at their current cell position
$\mathbf{r}_{i}$, i.e. $\mathbf{q}_{i}=(\cos(\theta_{i}),\sin(\theta_{i}),0)$.

###  III.1 Transient migration

With or without neighbor exchanges, both an aster ($\theta_{0}=0$) and a
vortex ($\theta_{0}=\pi/2$) give rise to a transient radial migration
corresponding to an outward motion for an aster (regardless on whether the
aster arrangement of cell polarities is inward or outward) and an inward
motion for a vortex. This divergent/convergent migration leads to
depletion/accumulation of cell density around an aster or a vortex,
respectively, as shown the Fig. (3, A). We also compute the typical patterns
of cell migration by averaging over different initial conditions, time and
noise seeds. These average profiles are illustrated for an aster (B panel) and
a vortex (C panel). The divergent migration from an aster is balanced out by
the isotropic attraction forces (tuned by $\lambda_{0}$ parameter) preventing
ruptures, such that the cells tend to stagnate into simply connected
configurations.

Notice that for an aster, the outward migration pattern aligns with the
divergent polarity profile, while for the vortex there is a reversal of motion
from inward in the bulk to outward on the rim of the disk. This reversal
effect is characteristic to $+1$ defects as predicted from a hydrodynamic
model in Ref. [30]. Interestingly, we find that it is also present for a
spiral defect as discussed later.

Figure 6: Time snapshots of the evolving $\mathbf{q}_{i}$ polarities for two
topological defects: A) Two initial vortices turning into spirals while
repelling each other, B) One vortex turning into a spiral as it approaches the
$m=-1$ defect, C) Two $m=-1$ defects repelling each other. Background field
shows the defect density charge $\rho$ field ($\lambda_{0}=1.5$).

For the $-1$ defect, the cell polarity alignment has a 4-fold saddle
structure. This induces an 8-fold saddle structure in the migration pattern as
shown in Fig. (4, B) and consistent with the theoretical prediction from a
hydrodynamic model of polar active matter [30]. The saddle point introduces
large frustrations in the cell-cell connectivity which are removed by local
neighbor exchanges. This spontaneous rearrangements occur along the principal
axes oriented at $\pi/4$ degrees with respect to the saddle point axes and
leads to the transient fluidization with the formation of an "X" shaped
tissue.

###  III.2 Persistent migration

The only configuration with persistent cell migration corresponds to a spiral
$\theta_{0}\in(0,\pi/2)$. Under neighbor exchanges, cells sustain a chiral
flow as shown in Fig. (4 A). The flow chirality is determined by that of the
spiral in the center of the disk and changes sign near the boundary of the
disk. This is similar to the reversal of transient divergent migration induced
by a vortex and the flow reversal predicted in Ref. [30]. To further bridge to
the hydrodynamic models, we compute the azimuthal speed for this vortical flow
as function of the radial distance to the spiral defect. In Fig. (5), we see
that it exhibits the $1/r$ scaling at intermediate scales away from the defect
core (discrete nature) and the finite boundary. Notice that the azimuthal
speed reverses its sign near the boundary layer, i.e. at $1-\xi$, which
implies that $\xi$ is the relevant lengthscale for this flow reversal.

In Fig. (6), we illustrate the typical orientational flow fields (normalized
vector field) induced by one pair of defects (two vortices, one vortex and one
$-1$ defect and two $-1$ defects). The cell migration pattern changes non-
trivially in the presence of multiple defects due to the long-range
interactions between defects.

Figure 7: Morphological response to an imprinted vortex in the absence of
neighbor exchanges. A) Initial Gaussian bump perturbation (with height
exaggerated for visualization purposes). B) Formation of a conical shape. C)
Profile of the instantaneous cell velocity along the $(x,z)$ cross-section for
different times. D) Relaxation to the asymptotic conical shape
($\lambda_{0}=2$).

###  III.3 Defect density field

To further characterize the polar order within the tissue and its topological
defects, we coarse-grain the $\mathbf{q}_{n}$ polarities to obtain a smooth
polarization field

$\mathbf{\Psi}(\mathbf{r}^{\prime})=\frac{1}{N}\sum\limits_{i=1}^{N}\mathcal{P}_{i}\cdot\mathbf{q}_{i}\delta(\mathbf{r}^{\prime}-\mathbf{r}^{\prime}_{i})$
(11)

corresponding to an $O(2)$ vector order parameter of the tissue polarization.
$\mathbf{r}^{\prime}=\mathcal{P}_{i}\cdot\mathbf{r}$ is the in-plane
coordinate and $\mathcal{P}_{i}=(\mathcal{I}-\mathbf{p}_{i}\mathbf{p}_{i})$ is
the in-plane projection operator for the $i$-th cell. Whilst the
$\mathbf{q}_{i}$ are unit vectors, the corresponding polarization field has
varying magnitude (i.e. the order parameter space is the unit disk). The
topological defects are associated with regions where the polarization
vanishes in magnitude and becomes multi-valued in phase. Thus, defects can be
tracked as zeros of the order parameter using the Mazenko-Halperin method
[31]. The defect density $\rho$ follows as

$\rho=\frac{1}{\pi\Psi_{0}^{2}}[(\partial_{x}\Psi_{x})(\partial_{y}\Psi_{y})-(\partial_{x}\Psi_{y})(\partial_{y}\Psi_{x})],\
$ (12)

where $\Psi_{0}$ is the constant magnitude of the uniform polarization.
Topological defects are initially embedded in the phase of the polarization
field as singularities

$\theta=\theta_{0}+\sum_{i}m_{i}\arctan\bigg{(}\frac{y-y_{0,i}}{x-x_{0,i}}\bigg{)}\
,$ (13)

with topological charges $m_{i}=\pm 1$. For two defects, we use the initial
position coordinates $(x_{0,i},y_{0,i})=(0,\pm 15)$ and evolve the cell
positions and the $\mathbf{q}_{i}$ polarities according to Eqs. (6) and (8),
where cells are allowed to exchange neighbors in order preserve the minimum
neighbor distance.

As shown in Fig. (6), the defect density field $\rho$ tracks very well the
topological charges of the moving defects. The defect core size is about two
cell units, thus small enough compared to the system size to observe
hydrodynamic effects and long-range interactions between defects. We see that
defects of same sign repels while those of opposite signs attract and
eventually annihilate, as predicted in hydrodynamic models.

Figure 8: Morphological response induced by a fixed vortex in the presence of
cell rearrangements. A) Tubular growth at t=500, B) at t = 1500. C) Profile of
the instantaneous cell velocity along the $(x,z)$ cross-section for different
times. The colormap shows the density of neighbor exchanges per cell
indicating localized fluidization around the evolving tube ($\lambda_{0}=2$).

##  IV  Cone versus Tube

We now relax the zero curvature constraint by allowing the AB polarities to
evolve according to Eq. (7), which will induce cell migration normal to the
surface. We study the morphological response of the tissue due to presence of
a vortex in the $\mathbf{q}_{i}$ polarities and see how cell neighbor
exchanges induce the formation of a tubular structure.

To induce a morphological change, we perturb the disk with a Gaussian bump of
small height $h=3$ cell units as illustrated in Fig. (7, A). Using the radial
symmetry of the vortex, we can extend the method of embedding a vortex in the
configuration of the $\mathbf{q}_{i}$ polarities from Sec. (III) by fixing the
orientation of the $\mathbf{q}_{i}$ along concentric rings centered at the
vortex position. The cell migration and the evolution of tissue curvature are
governed by the dynamics of $\mathbf{r}_{i}$’s and $\mathbf{p}_{i}$’s. At each
timestep, we evaluate the orientation of the $\mathbf{q}_{i}$ polarities along
concentric rings at a distance $r$ from the vortex using the static
orientational profile of $\theta$. Embedding other types of full-integer
defects for curved tissues is more challenging and remains to be studied
separately.

Figure 9: Time snapshots during the formation of a flap under vortex dynamics
and cell rearrangements: A) $t=575$, B) $t=625$, C) $t=675$, D) $t=1500$. The
colormap shows the density of neighbor exchanges per cell and indicates that
structural reconfiguration is localized around the evolving flap
($\lambda_{0}=2$).

The vortex is a source of inward migration and accumulation of cells towards
the Gaussian bump and this leads to a normal migration. From energetic
considerations, the vortex localises the energy associated with large
orientational distortions. In the absence of neighbor exchanges, the
collective migration is towards minimizing this total energy and thus
attaining an equilibrium shape. Subsequently, the normal migration is seeded
at the vortex center by the "leading" cells that pull the rest of the tissue
into a conical shape. The relaxation towards the cone is shown in Fig. (7
A-D). The overdamped dynamics slows down as the conical shape is approached
asymptotically. This equilibrium shape is attained also when we allow the PCP
polarities to evolve unconstrained.

Figure 10: Onset of out-of-plane growth. The fraction of cells that are more
than two cell widths above the flat surface as a function of time in semilog-
scale: flap, tube, cone I (rigid, moving defect) and cone II (rigid, fixed
defect) ($\lambda_{0}=2$).  Figure 11: The number of neighbor exchanges,
$N_{c}$, as function of time for tube and the flap shapes. Inset the
dependence in the log-log scale after the initial elastic response
($\lambda_{0}=2$).

By allowing the local neighbor rule to change the cell neighbor connectivity,
we enable localized fluidization, and this fundamentally alters the normal
migration path. For the vortex configuration, where the inwards migrations
within the tissue leads to normal migration, the local fluidization induces a
topological transition from a conical shape to a tubular structure as shown in
Fig. (8).

We introduce a scalar field defined by the cumulative number of neighbor
exchanges per cell to better quantify the fluidization at the onset of tube
formation and show that it is highly localized in space. This is shown in the
colormap of Fig. (8). The density of rearrangements is zero in the beginning
as the tissue responds by collective relaxation towards the conical shape.
Eventually, the curvature at the tip of the cone is sufficiently high to
trigger a local cell neighbor exchange which blunts the tip by opening a small
hole and marks the onset of tube formation as shown in Fig. (8). The
continuous cell rearrangements at the top and on the side of the tube
maintains the normal growth. Most of the cells away from the tube remain
fixated by their neighbors and migrate collectively. However, the fluidization
initiated by the normal growth permeates through the in-tissue migration on
localized streaks or shear zones.

Interestingly, the unconstrained evolution of the $\mathbf{q}_{i}$’s
polarities renders a morphological transition from a tube to flap structure as
shown in Fig. (9). At the base of the bump and in the far-field, the
concentric ordering of the polarities induced by a vortex is still maintained,
but the vortex center may move its position due to noise. These fluctuations
are enough to spontaneously break the rotationally-symmetric shape by inducing
a bend. This leads to a difference in curvature between the inside edge of the
bending tube (higher curvature) and the outside edge of the tube (lower
curvature). As a consequence, the side of tube with higher curvature stores
more energy and stagnates in the growth, whilst the side of the tube with
lower curvature continues to grow rendering a flap-like shape. Interestingly,
the ordering of the polarities on the flap retains the vortex-like structure.
This indicates that the tissue would rather change its shape than alter its
topological ordering, pointing to the robustness of the topological
constraint. The density of structural rearrangements remains localized on the
growing flap and along narrow steaks in the far-field akin to shear zones in
granular matter (see Fig. 9).

This topological robustness is also reflected in the global properties of the
growing shape and the total number of neighbor exchanges as shown in Figs.
(10) and (11). At the onset of normal growth, we see that the density of cells
inside the growing shape tends to increase logarithmically $\sim\ln t$ for
different classes of shapes with or without neighbor exchanges. We attribute
this logarithmic growth to the presence of the vortex guiding the growth
through the inward migration pattern, but this requires further theoretical
study. On the other hand, the total number of the neighbor exchanges tends to
increase algebraically with time as $t^{3/2}$ in the asymptotic limit. This
exponent is the same for both the tube and flap and indicates that it is
determined by the topological constraint imposed by the vortex. It also
implies that the rate of rearrangements, representing the flux of cells into
the fluidized (shear) zones, scales as $\sim t^{1/2}$, from the time
derivative of the density of neighbor exchanges.

To further emphasise the topological robustness, we consider the morphological
change induced by two vortices. For an intact tissue rigidity, we find that
the two vortices lead to collective normal migration with the formation of two
conical shapes as shown in Fig. (12 A). By contrast, two flaps are formed
under tissue fluidization and when vortices are allowed to move away from each
other due to their repelling interaction as shown in Fig. (12 B). We notice
that as the vortices turn into spirals as in Fig. (6), the flaps also develop
the same chirality as that of the corresponding spiral at the base.

##  V  Discussion and conclusions

Figure 12: Morphological changes of two initial Gaussian bumps due to two
interacting vortices. A) Conical shapes in the solid-limit, B) flap-like
shapes during tissue fluidization ($\lambda_{0}=2$).

In summary, we have studied the interplay of topological defects and tissue
fluidization on the morphological transformations of a tissue monolayer using
a minimal discrete model of interacting polarized cells. Topological defects
in the ordering of planar cell polarities pose robust topological constraints
on tissue shape. However, we find that the collective migration, which is
preferred by epithelial cells, gets fragmented by the localized tissue
fluidization during shape dynamics. In particular, we have shown that under
the topological constraint imposed by a vortex, the inwards migration within
the tissue leads to a normal migration in the region of high curvature. We
find that the onset of tube formation from a conical shape is attributed with
the onset of tissue fluidization occurring at the tip of the cone. Local
rearrangements may trigger each other similar to shear transformation zones in
amorphous materials leading to the formation of shear zones where the tissue
yields and flows.

The observed flap structures, extending from the tubular protrusion of the
epithelial surface, suggest a mechanism by which cells at the tip can be
displaced from the contiguous monolayer plane. This results in a cohort of
cells at the protrusion with both sides exposed to the extracellular
environment, potentially facilitating enhanced interaction with their
surroundings. A similar phenomenon occurs during angiogenesis, the formation
of new blood vessels. Specialized tip cells, known as endothelial sprouts,
breach the contiguous endothelial monolayer through a process called
sprouting. Subsequent cell division and pulling by neighboring cells lead to
the formation of a collective of specialized endothelial cells at the sprout
tip. These tip cells exhibit invasive properties, enabling the growth of the
new blood vessel [21, 22, 32].

While several studies have shown the ability of developing tissues to regulate
their cellular fluidity during morphogenesis, the role of fluidization in
epithelial tubulogenesis remains largely unexplored. Interestingly, a recent
work highlights the potential role of regulating tissue fluidity in the
spiraling alignment of endothelial cells within the tubular structure of blood
vessels [33]. Thus, our computational model, which demonstrates an increased
propensity for local fluidity during epithelial tube formation, suggests that
the regulation of tissue fluidity might be a general and critical process
during epithelial tube morphogenesis.

Currently, the model does not integrate the influence of the micro-
environment, which comprises essential extrinsic factors like the
extracellular matrix (ECM), signaling molecules, and mechanical cues.
Additionally, the model does not consider cellular processes such as cell
division, extrusion, and apoptosis, which could be vital for remodeling and
homeostasis of the epithelial _in vivo_. Yet, we have shown that this minimal
model offers a valuable advantage by capturing the essential, intrinsic
mechanisms governing epithelial cell sheets, such as polarity and inter-
cellular adhesions. This suggests that intrinsic cell properties, combined
with topological constraints (which may also be posed by the micro-
environment) which trigger the tissue response to converge towards a $+1$
defect core, might be sufficient for building epithelial tubes. Thus, the
present study highlights the importance of fundamental cell-intrinsic
principles in shaping complex tissues. Further investigations could explore
how these intrinsic properties of epithelial monolayers cooperate with
extrinsic factors, such as mechanical cues and biochemical signaling, to
facilitate the formation of complex organ structures. These provide many
future directions of research for this model, with the present work serving as
a basis.

###### Acknowledgements.

This project was funded by the University of Oslo, UiO:Life Science through
the convergence environment ITOM.

## References

  * Drubin and Nelson [1996] D. G. Drubin and W. J. Nelson, Origins of cell polarity, Cell 84, 335 (1996).
  * Zallen [2007] J. A. Zallen, Planar polarity and tissue morphogenesis, Cell 129, 1051 (2007).
  * Apodaca  _et al._ [2012] G. Apodaca, L. I. Gallo, and D. M. Bryant, Role of membrane traffic in the generation of epithelial cell asymmetry, Nature cell biology 14, 1235 (2012).
  * Butler and Wallingford [2017] M. T. Butler and J. B. Wallingford, Planar cell polarity in development and disease, Nature reviews Molecular cell biology 18, 375 (2017).
  * Montcouquiol  _et al._ [2003] M. Montcouquiol, R. A. Rachel, P. J. Lanford, N. G. Copeland, N. A. Jenkins, and M. W. Kelley, Identification of vangl2 and scrb1 as planar polarity genes in mammals, Nature 423, 173 (2003).
  * Wang  _et al._ [2006] Y. Wang, T. Badea, and J. Nathans, Order from disorder: Self-organization in mammalian hair patterning, Proceedings of the National Academy of Sciences 103, 19800 (2006).
  * Cetera and Horne-Badovinac [2015] M. Cetera and S. Horne-Badovinac, Round and round gets you somewhere: collective cell migration and planar polarity in elongating drosophila egg chambers, Current opinion in genetics & development 32, 10 (2015).
  * Stedden  _et al._ [2019] C. G. Stedden, W. Menegas, A. L. Zajac, A. M. Williams, S. Cheng, E. Özkan, and S. Horne-Badovinac, Planar-polarized semaphorin-5c and plexin a promote the collective migration of epithelial cells in drosophila, Current Biology 29, 908 (2019).
  * Classen  _et al._ [2005] A.-K. Classen, K. I. Anderson, E. Marois, and S. Eaton, Hexagonal packing of drosophila wing epithelial cells by the planar cell polarity pathway, Developmental cell 9, 805 (2005).
  * Gong  _et al._ [2004] Y. Gong, C. Mo, and S. E. Fraser, Planar cell polarity signalling controls cell division orientation during zebrafish gastrulation, Nature 430, 689 (2004).
  * Nishimura  _et al._ [2012] T. Nishimura, H. Honda, and M. Takeichi, Planar cell polarity links axes of spatial dynamics in neural-tube closure, Cell 149, 1084 (2012).
  * Francou and Anderson [2020] A. Francou and K. V. Anderson, The epithelial-to-mesenchymal transition in development and cancer, Annual review of cancer biology 4, 197 (2020).
  * Tetley  _et al._ [2019] R. J. Tetley, M. F. Staddon, D. Heller, A. Hoppe, S. Banerjee, and Y. Mao, Tissue fluidity promotes epithelial wound healing, Nature physics 15, 1195 (2019).
  * Hannezo and Heisenberg [2022] E. Hannezo and C.-P. Heisenberg, Rigidity transitions in development and disease, Trends in Cell Biology (2022).
  * Petridou  _et al._ [2021] N. I. Petridou, B. Corominas-Murtra, C.-P. Heisenberg, and E. Hannezo, Rigidity percolation uncovers a structural basis for embryonic tissue phase transitions, Cell 184, 1914 (2021).
  * Pinheiro and Mitchel [2024] D. Pinheiro and J. Mitchel, Pulling the strings on solid-to-liquid phase transitions in cell collectives, Current Opinion in Cell Biology 86, 102310 (2024).
  * Maroudas-Sacks  _et al._ [2021] Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun, and K. Keren, Topological defects in the nematic order of actin fibres as organization centres of hydra morphogenesis, Nature Physics 17, 251 (2021).
  * Gehart and Clevers [2019] H. Gehart and H. Clevers, Tales from the crypt: new insights into intestinal stem cells, Nature reviews Gastroenterology & hepatology 16, 19 (2019).
  * Pérez-González  _et al._ [2021] C. Pérez-González, G. Ceada, F. Greco, M. Matejčić, M. Gómez-González, N. Castro, A. Menendez, S. Kale, D. Krndija, A. G. Clark, _et al._ , Mechanical compartmentalization of the intestinal organoid enables crypt folding and collective cell migration, Nature cell biology 23, 745 (2021).
  * Ritsma  _et al._ [2014] L. Ritsma, S. I. Ellenbroek, A. Zomer, H. J. Snippert, F. J. De Sauvage, B. D. Simons, H. Clevers, and J. Van Rheenen, Intestinal crypt homeostasis revealed at single-stem-cell level by in vivo live imaging, Nature 507, 362 (2014).
  * Eilken and Adams [2010] H. M. Eilken and R. H. Adams, Dynamics of endothelial cell behavior in sprouting angiogenesis, Current opinion in cell biology 22, 617 (2010).
  * Hendriks and Ramasamy [2020] M. Hendriks and S. K. Ramasamy, Blood vessels and vascular niches in bone development and physiological remodeling, Frontiers in Cell and Developmental Biology 8, 602278 (2020).
  * Doxzen  _et al._ [2013] K. Doxzen, S. R. K. Vedula, M. C. Leong, H. Hirata, N. S. Gov, A. J. Kabla, B. Ladoux, and C. T. Lim, Guidance of collective cell migration by substrate geometry, Integrative biology 5, 1026 (2013).
  * Heinrich  _et al._ [2020] M. A. Heinrich, R. Alert, J. M. LaChance, T. J. Zajdel, A. Košmrlj, and D. J. Cohen, Size-dependent patterns of cell proliferation and migration in freely-expanding epithelia, Elife 9, e58945 (2020).
  * Li and Sun [2014] B. Li and S. X. Sun, Coherent motions in confluent cell monolayer sheets, Biophysical journal 107, 1532 (2014).
  * Lång  _et al._ [2022] E. Lång, C. Pedersen, A. Lång, P. Blicher, A. Klungland, A. Carlson, and S. O. Bøe, Mechanical coupling of supracellular stress amplification and tissue fluidization during exit from quiescence, Proceedings of the National Academy of Sciences 119, e2201328119 (2022).
  * Segerer  _et al._ [2015] F. J. Segerer, F. Thüroff, A. P. Alberola, E. Frey, and J. O. Rädler, Emergence and persistence of collective cell migration on small circular micropatterns, Physical Review Letters 114, 228102 (2015).
  * Nissen  _et al._ [2017] S. B. Nissen, M. Perera, J. M. Gonzalez, S. M. Morgani, M. H. Jensen, K. Sneppen, J. M. Brickman, and A. Trusina, Four simple rules that are sufficient to generate the mammalian blastocyst, PLoS biology 15, e2000737 (2017).
  * Nissen  _et al._ [2018] S. B. Nissen, S. Rønhild, A. Trusina, and K. Sneppen, Theoretical tool bridging cell polarities with development of robust morphologies, Elife 7, e38407 (2018).
  * Rønning  _et al._ [2023] J. Rønning, J. Renaud, A. Doostmohammadi, and L. Angheluta, Spontaneous flows and dynamics of full-integer topological defects in polar active matter, Soft Matter 19, 7513 (2023).
  * Skogvoll  _et al._ [2023] V. Skogvoll, J. Rønning, M. Salvalaglio, and L. Angheluta, A unified field theory of topological defects and non-linear local excitations, npj Comput Mater 122, https://doi.org/10.1038/s41524-023-01077-6 (2023).
  * Geudens and Gerhardt [2011] I. Geudens and H. Gerhardt, Coordinating cell behaviour during blood vessel formation, Development 138, 4569 (2011).
  * Zhang  _et al._ [2024] H. Zhang, T. Rahman, S. Lu, A. P. Adam, and L. Q. Wan, Helical vasculogenesis driven by cell chirality, Science Advances 10, eadj3582 (2024).