EPJ Web of Conferences Lattice2017  11institutetext: Institute of Physics,
LPPC, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
22institutetext: Fakultät für Physik, Universität Bielefeld, D-33615,
Bielefeld, Germany

# Thermal Simulations, Open Boundary Conditions and Switches

Yannis Burnier  Adrien Florio Speaker, 11<EMAIL_ADDRESS>Olaf Kaczmarek
22 Lukas Mazur 22

###### Abstract

$SU(N)$ gauge theories on compact spaces have a non-trivial vacuum structure
characterized by a countable set of topological sectors and their topological
charge. In lattice simulations, every topological sector needs to be explored
a number of times which reflects its weight in the path integral. Current
lattice simulations are impeded by the so-called freezing of the topological
charge problem. As the continuum is approached, energy barriers between
topological sectors become well defined and the simulations get trapped in a
given sector. A possible way out was introduced by Lüscher and Schaefer using
open boundary condition in the time extent. However, this solution cannot be
used for thermal simulations, where the time direction is required to be
periodic. In this proceedings, we present results obtained using open boundary
conditions in space, at non-zero temperature. With these conditions, the
topological charge is not quantized and the topological barriers are lifted. A
downside of this method are the strong finite-size effects introduced by the
boundary conditions. We also present some exploratory results which show how
these conditions could be used on an algorithmic level to âreshuffleâ the
system and generate periodic configurations with non-zero topological charge.

##  1 Introduction

$SU(N)$ gauge theories on compact and orientable manifolds admit a set of
distinct vacua, labeled by an integer topological charge $Q$ DeWitt1979 . When
computing the path integral, all topological sectors need to be included
(9780521670548, , chap. 23.6). Conventional lattice QCD simulations use
periodic boundary conditions; they take place on the torus $\mathbf{T}^{4}$.
Hence, different lattice configurations are characterized by different $Q$’s
(for a review of the topology of $SU(N)$ theories on the torus, see Gonzalez-
Arroyo1997 ).

As the continuum limit is approached, the topological sectors become well
defined Luscher1982c , they decouple and small deformations of the
configurations do not allow to travel from one sector to the other. The
simulations get stuck in a given sector; this problem is often referred to as
"topological freezing" and has been investigated in the literature, see for
example Schaefer2011 ; Fritzsch2014 .

Different approaches to solve this issue have been or are studied Luscher2011b
; Mages2015a ; Czaban2014 ; Bietenholz2016 ; thisContriHasen . In particular,
the use of open boundary conditions (OBC) in the time extent has been proposed
in Luscher2011b . In this context, the base manifold of the simulations is not
compact anymore, $Q$ is not quantized and the topological freezing disappears.
An important drawback is the appearance of strong finite-size effects
Mages2015a ; thiscontribSommer .

This contribution is constructed as follow. In section 2, we discuss the
topological freezing at non-zero temperature and present quenched results
obtained with OBC in the spatial extents. This is relevant for the computation
of the topological susceptibility at non-zero temperatures, which is of utter
interest for cosmology and axion physics thiscontriMoore . It is also needed
for a precise determination of transport coefficients, whose correlators are
known to be sensitive to the topological sampling Fukaya2014b . We show that,
as expected, the topological charge is not quantized and the topological
freezing disappears. However, atop of the topological charge not being anymore
a topological invariant, OBC configurations are plagued by strong finite-size
effects.

Thus, in section 3, we try to circumvent these problems and present quenched
results where OBC are used purely as an algorithmic tool111A very related idea
was presented by M. Hasenbusch thisContriHasen  at this conference.. The idea
is to generate a certain number of configurations with PBC, then switch to OBC
in order to "reshuffle" the system, before going back to PBC. It manages to
produce higher topological charge configurations but oversamples the topology.
We conclude by discussing outlooks on how this method could be improved.

##  2 Topological Freezing and Spatial OBC

To overcome the problem of topological freezing, the idea of using OBC was
introduced in Luscher2011b . By making the manifold non-compact, one loses the
quantization of the topological charge; there are no more distinct topological
sectors DeWitt1979 ; Luscher2011b . However, to be able to simulate thermal
configurations, one needs periodicity in time Gross1981 . This enforces the
use of OBC in spatial directions. To preserve cubic symmetry, we impose them
in all three directions, thus setting222As explained in the original paper
Luscher2011b , one also needs to introduce some weights for the boundary
plaquettes in the action.

$F_{i\mu}(x)|_{x_{i}=-a}=F_{i\mu}(x)|_{x_{i}=N_{i}a}=0,\ \ \ i=1,\dots,3$ (1)

for a lattice with lattice spacing $a$ and sizes $N_{t}\times N_{x}\times
N_{y}\times N_{z}$.

We start by doing some measurements on $T=0$ configurations, $20^{4}$ lattices
at $\beta=6.2$, where the topological transitions are not damped by the quark-
gluon plasma. To compute the topological charge, we use Wilson’s flow to
evolve the configurations towards the continuum together with the clover
discretization of the field-strength tensor. When the continuum theory has
distinct topological sectors, the flow dynamically interpolates every
configuration to a given sector Luescher2010 . This is illustrated by the
curves in the left hand-side part of figure 1. The upper plot is a semi-log
plot of the topological charge’s evolution as a function of the dimensionless
flow time for given configurations. The topological charge is not quantized at
$\tau=0$. First, in the case of PBC (green and red curves), when $\tau$ is
increased, the discrete configurations are flown towards continuous bundles.
They are representative of different topologies; they carry different integer
topological charges. We see that this quantity is unambiguously determined and
very robust under flow evolution. The less stable $Q=-1$ charge decays at
about $\tau=600$, which corresponds to a smoothing range of roughly
$\rho\approx 70$ lattice sites ($\rho\approx\sqrt{8\tau}$, see Luscher2011b ),
more than three times the actual time extent of the lattice. The observed fact
that the $|Q|=2$ charge is more stable than the $|Q|=1$ is a general feature,
which had already been observed on the lattice in BilsonThompson:2003zi ,
using some smearing methods. On the torus, there exists no self-dual
configuration for $|Q|<2$, see Braam1989  and references therein. The $|Q|=1$
sectors only have infima of their action but no minima; they are less stable.

On the other hand, when dealing with OBC (blue curves), no plateau is reached;
the charge continuously flows towards zero. It reflects the non-existence of
distinct topological sectors; this is the feature which allows for a good
sampling of the topological fluctuations.

To pursue this idea, we simulated some toy lattices (small aspect ratios and
volumes) at $T=1.28T_{c}$ for different lattice spacings. In the right hand-
side part of figure 1, we show the ratio of the topological charge variance by
the lattice volumes, both for PBC and OBC configurations. For PBC, as the
coupling is reduced, the topological susceptibility drops; different sectors
are not sampled anymore; the simulation freezes. For OBC, as expected, it does
not occur. No topological barriers impede the Monte-Carlo (MC) evolution; the
topological fluctuations are well sampled.

Figure 1: Left: Evolution of the topological charge with the dimensionless
flow time, $T=0$. For PBC, it gives rise to a well-defined integer, robust to
Wilson’s flow evolution. In the case of OBC, there is no clear mean of
singling out a flow time at which the topological charge should be defined.
The lower plot shows the derivatives of the curve and displays a clear
difference between the two types of conditions. Upper right: Topological
charge variance for different couplings, $T=1.28T_{c}$. In the case of PBC, we
see that the topology freezes. This does not happen in the OBC case. Lower
right: Topological charge distributions for OBC (black) and PBC (blue). This
clearly illustrates the fact that OBC configurations do not possess an integer
topological charge. Figure 2: Finite size effects. We restrict the spatial
extent to sub-cubes of volume $L_{s}^{3}$ (drawing on the right-hand side) and
plot the expectation values of Wilson loops as a function of $L_{s}$. As
expected, the OBC bulk results match the PBC expectation value. Note however
the large extent where the finite size effects dominate.

This is of course not the whole story. The OBC come with their own
disadvantages, the main one being the strength of finite size effects (for a
careful discussion in the case of time-like OBC, the reader is referred to
thiscontribSommer ). This problem is illustrated for the spatial OBC in figure
2. There, as pictured on the right hand-side, we look at the average $4\times
4$ Wilson loops on sub-spatial volumes of $32^{4}$ lattices with OBC,
$\beta=6.872$. As expected, the bulk value matches with the PBC value, where
the finite-size effects are minimal. However, we observe that to get rid of
them, we need to restrict the OBC configurations to $12^{3}\times 32$ sub-
lattices, effectively losing a substantial part of the volume. This makes the
use of small OBC lattices unattractive, while being still practical for larger
lattices. A milder disadvantage of OBC is related to the topology. The
determination of the topological susceptibility is made harder; one does not
have anymore a clear stopping criterion for the Wilson flow (for PBC, it is
given by the quantization of the topological charge). Moreover, it needs to be
computed directly from the two point functions of the topological charge
density, see Bonati:2017woi  for example.

In the following section, we present some exploratory results on how these
problems might be overcome by using the OBC as an algorithmic tool.

##  3 Switching the Boundaries

In order to cure the problems related to OBC, we try to combine them in the
process of generating PBC configurations; we want to use them solely to
"reshuffle" the topology. As we may see in figure 3 plot, which depicts the
evolution of the topological charge with OBC as a function of MC sweeps, OBC
lead to a very efficient sampling of the topology. Our proposed algorithm is
illustrated in the right-hand side of the same figure. We start a normal MC
run with PBC. After $n_{PBC}/n_{sweeps}$ configurations, we switch PBC for OBC
in one spatial direction, perform $n_{OBC}$ sweeps in order to change
topological sector and then switch back to PBC. We repeat this $n_{trans}$
times, to have the required number of topological transitions.

Figure 3: Right: Topological charge history during a MC run. Note that the
$x$-axis represents the number of MC sweeps and not a configuration number
(configurations are separated by 500 sweeps). We see that the autocorrelation
time of the topological charge with OBC is very short; the sampling is very
efficient. Left: Illustration of how the "switch" algorithm works.
Measurements are performed on periodic configurations. Re-thermalization is
required after the OBC sweeps.

We keep only the PBC configurations which can be considered to have re-
thermalized (see the end of this section for a discussion on this issue) after
the OBC sweeps. Moreover, every time we proceed to a new set of OBC sweeps, we
"open-up" a different direction, in order to homogenize the potential finite-
size biases hence introduced. Our first results on $64^{3}\times 16$ lattices
at $\beta=6.872$ ($1.28T_{c}$) are shown in figure 4. The results obtained by
doing the switches are displayed in blue. We see that an efficient sampling of
the different topological sectors is achieved. It allows us to compute the
topological susceptibility, to find $\frac{\chi_{switch}}{T_{c}^{4}}\approx
0.16$. Actually, the sampling seems to be too efficient. The PBC results are
shown in black. The measured topological susceptibility
$\frac{\chi_{PBC}}{T_{c}^{4}}\approx 0.0045$ is in agreement with the value of
Borsanyi2016 . The low sampling of the topology seems to be a physical effect
due to the damping of the quark-gluon plasma and does not seem to be a sign of
topological freezing.

This potential problem of oversampling is most certainly explained by the fact
that brutally "opening-up" our lattice torus to allow the topological charge
to flow in and out disturbs durably the systems. To reduce this issue, as
mentioned before, we re-thermalized our configurations after every OBC
sequence.

This operation is delicate for the following reason: it is not an easy task to
disentangle the finite size-effects related to fixed topology from spurious
effects coming from badly thermalized lattices. To proceed, we looked at
Wilson loops, computed on subset of configurations generated continuously
after an OBC transition. We define re-thermalized lattices to be the ones
which, in a part between two OBC sequences, lead to a set of Wilson loops
expectations values which consistently comes from a single probability
distribution. It is not completely clear how appropriate this procedure is;
gaining a better understanding on the re-thermalization is a potential way of
improving these results.

As suggested to us during this conference (and following the idea of
thisContriHasen ), a seemingly more rigorous way of building such an algorithm
could involve some parallel tempering between OBC and PBC configurations. This
would remove the re-thermalization step and provide an algorithm whose
properties are better defined. Implementing such an algorithm would also allow
us to adjust or rule out our naive switch algorithm.

Figure 4:  Topological charge distribution, for the $64^{3}\times 16$ PBC
(black) and "switch" (blue) configuration sets. The PBC results agree with the
results of Borsanyi2016 ; the "switch" algorithm seems to lead to some
oversampling. Only re-thermalized configurations are shown. Left: Histogram.
Right: MC-History.

##  4 Conclusion

In lattice simulations, a good sampling of the different topological sectors
is desirable, both for thermal and non-thermal simulations. In order to try to
achieve such a sampling, we presented some results obtained with OBC imposed
in spatial directions, having in mind to use them in the former case. We
showed that, as expected, the topological charge is not quantized and has a
very small auto-correlation time. We also showed that these conditions suffer
from large finite-size effects. Building on these observations, we presented
some exploratory results where we made use of OBC as a tool to change
topological sectors between different PBC configurations. Our naive algorithm
seems to lead to some oversampling. It is worth noting that even if it also
leads to a wrong topological sampling, it has the merit of providing a new way
to generate higher charge configurations, which could for example be used in
some re-weighting procedures.

All in all, these first results are encouraging, but work remains to be done.
First, the naive algorithm should be further adjusted, by gaining a better
understanding of the re-thermalization procedure. Then, one should try to
design a more subtle algorithm, using some parallel tempering methods for
example. The next step would be to study fermionic configurations, one of the
final objectives being to be able to use it as a different way to extract the
topological susceptibility from lattice simulations.

### Acknowledgments

A.F. would like to thank A. Rothkopf, who kindly accepted to share his code
Rothkopf:2011db . It was modified and used to generate configurations in the
initial part of this project. O. K. and L. M. acknowledge support by the
Deutsche Forschungsgemeinschaft (DFG) through the grant CRC-TR 211 "Strong-
interaction matter under extreme conditions". The work of Y.B. and A.F. was
supported by the Swiss National Science Foundation.

## References

  * (1) B.S. DeWitt, C.F. Hart, C.J. Isham, Physica A96, 197 (1979)
  * (2) S. Weinberg, _The Quantum Theory of Fields, Volume 2: Modern Applications_ (Cambridge University Press, 2005), ISBN 0521670543
  * (3) A. Gonzalez-Arroyo, _Yang-Mills fields on the four-dimensional torus. Part 1.: Classical theory_ , in _Nonperturbative quantum field physics. Proceedings, Advanced School, Peniscola, Spain, June 2-6, 1997_ (1997), pp. 57–91, hep-th/9807108
  * (4) M. Lüscher, Commun. Math. Phys. 85, 39 (1982)
  * (5) S. Schaefer, R. Sommer, F. Virotta (ALPHA), Nucl. Phys. B845, 93 (2011), 1009.5228
  * (6) P. Fritzsch, A. Ramos, F. Stollenwerk, PoS Lattice2013, 461 (2014), 1311.7304
  * (7) M. Lüscher, S. Schaefer, JHEP 07, 036 (2011), 1105.4749
  * (8) S. Mages, B.C. Toth, S. Borsanyi, Z. Fodor, S. Katz, K.K. Szabo (2015), 1512.06804
  * (9) C. Czaban, A. Dromard, M. Wagner, Acta Phys. Polon. Supp. 7, 551 (2014), 1404.3597
  * (10) W. Bietenholz, K. Cichy, P. de Forcrand, A. Dromard, U. Gerber, _The Slab Method to Measure the Topological Susceptibility_ , in _Proceedings, 34th International Symposium on Lattice Field Theory (Lattice 2016): Southampton, UK, July 24-30, 2016_ (2016), 1610.00685
  * (11) M. Hasenbusch, _Fighting topological freezing in the two-dimensional $CP^{N-1}$ model _, in _Proceedings,35th International Symposium on Lattice Field Theory (Lattice2017): Granada, Spain_, to appear in EPJ Web Conf., 1709.09460
  * (12) N. Husung, M. Koren, P. Krah, R. Sommer, _$SU(3)$ Yang Mills theory at small distances and fine lattices_, in _Proceedings,35th International Symposium on Lattice Field Theory (Lattice2017): Granada, Spain_, to appear in EPJ Web Conf.
  * (13) G. Moore, _Axion dark matter and the Lattice_ , in _Proceedings,35th International Symposium on Lattice Field Theory (Lattice2017): Granada, Spain_, to appear in EPJ Web Conf., 1709.09466
  * (14) H. Fukaya, S. Aoki, G. Cossu, S. Hashimoto, T. Kaneko, J. Noaki (JLQCD), PoS LATTICE2014, 323 (2014), 1411.1473
  * (15) D.J. Gross, R.D. Pisarski, L.G. Yaffe, Rev. Mod. Phys. 53, 43 (1981)
  * (16) M. Lüscher, JHEP 08, 071 (2010), [Erratum: JHEP03,092(2014)], 1006.4518
  * (17) S.O. Bilson-Thompson, D.B. Leinweber, A.G. Williams, G.V. Dunne, Annals Phys. 311, 267 (2004), hep-lat/0306010
  * (18) P.J. Braam, P. van Baal, Commun. Math. Phys. 122, 267 (1989)
  * (19) C. Bonati, M. D’Elia (2017), 1709.10034
  * (20) S. Borsanyi, M. Dierigl, Z. Fodor, S.D. Katz, S.W. Mages, D. Nogradi, J. Redondo, A. Ringwald, K.K. Szabo, Phys. Lett. B752, 175 (2016), 1508.06917
  * (21) A. Rothkopf, T. Hatsuda, S. Sasaki, Phys. Rev. Lett. 108, 162001 (2012), 1108.1579


††thanks: Corresponding Auhtor: Dibakar Ghosh

# Directional synchrony among self-propelled particles under spatial influence

Suvam Pal  Physics and Applied Mathematics Unit, Indian Statistical Institute,
203 B. T. Road, Kolkata 700108, India  Gourab Kumar Sar  Physics and Applied
Mathematics Unit, Indian Statistical Institute, 203 B. T. Road, Kolkata
700108, India  Dibakar Ghosh  Physics and Applied Mathematics Unit, Indian
Statistical Institute, 203 B. T. Road, Kolkata 700108, India
<EMAIL_ADDRESS>(D. Ghosh) Arnab Pal  The Institute of Mathematical
Sciences, CIT Campus, Taramani, Chennai 600113, India & Homi Bhabha National
Institute, Training School Complex, Anushakti Nagar, Mumbai 400094, India

###### Abstract

Synchronization is one of the emerging collective phenomena in interacting
particle systems. Its ubiquitous presence in nature, science, and technology
has fascinated the scientific community over the decades. Moreover, a great
deal of research has been, and is still being, devoted to understand various
physical aspects of the subject. In particular, the study of interacting
active particles has led to exotic phase transitions in such systems which
have opened up a new research front-line. Motivated by this line of work, in
this paper, we study the directional synchrony among self-propelled particles.
These particles move inside a bounded region, and crucially their directions
are also coupled with spatial degrees of freedom. We assume that the
directional coupling between two particles is influenced by the relative
spatial distance which changes over time. Furthermore, the nature of the
influence is considered to be both short and long-ranged. We explore the phase
transition scenario in both the cases and propose an approximation technique
which enables us to analytically find the critical transition point. The
results are further supported with numerical simulations. Our results have
potential importance in the study of active systems like bird flocks, fish
schools and swarming robots where spatial influence plays a pertinent role.

> Swarming and alignment are two rich collective behaviors in statistical
> physics and nonlinear dynamics. Understanding these behaviors has been a
> challenge to the researchers over the years. These phenomena primarily
> depend on the interactions among the agents or potential landscape
> (attractive/repulsive force). In VicsekVicsek  _et al._ (1995) model, the
> vision range of each mobile agent plays a crucial role to control the
> alignment among them in the presence of noise. In this paper, we have
> replaced the vision range of the particles by distance dependent coupling
> functions. This type of influence of spatial positions on the directions of
> the particles is generic and allows us to analytically calculate the
> critical transition point from a desychronized configuration to a
> synchronized one by performing mean-field approximation.

##  I Introduction

The study of collective dynamics comes into picture when one tries to
comprehend the behavior of interacting entities that stay or move together.
Phenomena like swarming Carrillo  _et al._ (2010); Okubo (1986) and flocking
Reynolds (1987); Toner and Tu (1998) are some of the most widespread
ramifications of collective dynamics with multi-agent systems. A myriad of
applications can be found in many biological systems like bacterial colonies
Czirók  _et al._ (1996); Kearns (2010), bird flocks Cavagna  _et al._ (2010);
Hemelrijk and Hildenbrandt (2012), shoaling behavior in fishes Huth and Wissel
(1994); Partridge (1982) where the length scales span from micro (microns) to
macro (kilometers). Specifically, it is known that local interaction among
large number of individuals drives them to a state of collective motion, which
is commonly known as swarming. Generally, these systems have been demonstrated
using self-propelled particles (SPPs) where the particles consume energy from
the surrounding environment to drive themselves. They move with nonzero self
propulsion velocity and orientation. The situation where the velocity and
direction of propagation are synchronized among the agents, is regarded as
flocking in literature.

Several multi-agent systems in nature can show synchronization along with
collective spatial dynamics. Simply put, synchronization can be perceived as a
state of unison. In 1967, Winfree Winfree (1967) proposed a limit cycle
oscillator model that shows critically synchronized behavior. Later, Kuramoto
proposed a general model which was analytically tractable Kuramoto (1975).
Kuramoto’s seminal paper on synchronization opened up a new avenue of research
that lead to further studies of this model with different network topologies
Rodrigues  _et al._ (2016); Dörfler and Bullo (2011), distributed frequencies
Sakaguchi and Kuramoto (1986), time-varying interactions Ghosh  _et al._
(2022), attractive-repulsive couplings Hong and Strogatz (2011a), multi-layer
structures Gambuzza, Frasca, and Gomez-Gardenes (2015), higher-order
interactions Millán, Torres, and Bianconi (2020) etc. Synchronization
phenomena have also been rigorously studied with various amplitude oscillator
models Pikovsky  _et al._ (2001); Pecora and Carroll (1990); Boccaletti  _et
al._ (2002); Anwar and Ghosh (2023) other than the regular Kuramoto phase
oscillator. However, in most of these studies, the spatial positions of
particles played a minimal role in the sense that they did not have an effect
on the dynamics of synchronization. This gap was curtailed to some extent once
researchers started studying mobile oscillators Frasca  _et al._ (2008);
Fujiwara, Kurths, and Díaz-Guilera (2011); Chowdhury  _et al._ (2019); Majhi,
Ghosh, and Kurths (2019). In there, the particles, which are more often termed
as oscillators, move in space and interact among them depending on their
spatial locations. More recently, extensive works have been done on model
systems where the above-mentioned unidirectional influence of spatial
positions on internal dynamics has been replaced by a bidirectional mutual
feedback of spatial dynamics and internal dynamics on each other. These model
systems, termed as swarmalators O’Keeffe, Hong, and Strogatz (2017); Sar and
Ghosh (2022); O’Keeffe and Bettstetter (2019); Sar, Ghosh, and O’Keeffe
(2023); O’Keeffe, Ceron, and Petersen (2022); Sar  _et al._ (2022); Sar,
O’Keeffe, and Ghosh (2023), are still very much in their infancy and require
extensive study.

In this work, we are interested in the collective dynamics of SPPs with a
sharp distinction that the spatial positions influence the direction of
propagation. Similar study in this direction was first carried out by Vicsek
et al. Vicsek  _et al._ (1995) while studying phase transition among
interacting active particles. The agents switched from a gas-like incoherent
state to one in which the agents moved as a coherent flock beyond a critical
density. It was assumed that there is a vision range for each particle and the
direction of that particle is the mean direction of all the particles lying
inside that vision range. At later times, many studies adopted this approach
and considered short-range interaction of particles Peruani, Nicola, and
Morelli (2010); Peruani (2017); Großmann, Peruani, and Bär (2016). However,
the coupling strengths inside this interaction region were not affected by the
spatial distance between the particles. We replace this type of implicit
dependence by incorporating explicit functional dependence of spatial
positions to the direction of propagation or the orientation of the particles.

To interconnect the spatial movement and the orientation of the particles, we
consider SPPs which are moving on a plane and introduce a spatial dependent
Kuramoto-like phase equation. The orientation dynamics is influenced by the
distance between the particles in two ways. In one case, the distance between
two particles enhances their directional coupling and the effect is inverse in
the other case. Here, the directional coupling essentially indicates the
coupling among particles’ direction. Due to the noise present in the system,
particles move in randomized directions throughout the region for low coupling
strength, but for sufficiently large coupling their orientations tend to
align. We study this transition from a randomized state to an ordered state
with a mix and match of theoretical and/or numerical analysis. In the
thermodynamic limit, we analytically derive the effective coupling strength
and critical noise amplitude for the synchronized state. We also study the
rate of change of the order parameter near the critical noise strength.

The rest of the paper is arranged as follows: In section Sec. II, we introduce
and describe the model for SPPs. We proceed with the theoretical analysis of
our model in Sec. III, where we elaborately discuss the nature of spatial
influence and its dependence on the exponent $\alpha$. In Sec. IV, we study
the behavior of macroscopic order parameter near the critical transition point
by performing finite size scaling analysis on the model. Finally, we summarize
our paper with concluding remarks in Sec. V.

##  II Model set-up of SPP

We consider $N$ number of identical self-propelled particles moving inside a
square region with area $L\times L$ under periodic boundary conditions with
speed $v_{0}$. The spatial dynamics is governed by

$\dot{\vec{R}}_{i}=v_{0}\left(\cos\theta_{i}\hat{x}+\sin\theta_{i}\hat{y}\right),$
(1)

for $i=1,2,\cdots,N$, where $\vec{R}_{i}\equiv(x_{i},y_{i})$ stands for the
position of the $i$-th particle, and $\hat{x}$ and $\hat{y}$ are the unit
vectors along $x$ and $y$ axes, respectively. Here $\theta_{i}$’s represent
the internal dynamics (in our case, the direction of propagation) of the
particles. The coupled internal dynamics is given by the following equation

$\dot{\theta}_{i}=\omega_{i}+\dfrac{K}{N}\sum_{j=1,j\neq
i}^{N}\left|\vec{R}_{j}-\vec{R}_{i}\right|^{\alpha}\sin\left(\theta_{j}-\theta_{i}\right)+\zeta_{i},$
(2)

where $\omega_{i}$ signifies the angular frequency of the ith particle and
which is identical for all the particles, i.e., $\omega_{i}=\omega$, for
$i=1,2,\cdots,N$. By moving to a proper frame of reference, we can set
$\omega=0$ Hong and Strogatz (2011b). The sinusoidal coupling among the
directional component of the particles is inspired by the Kuramoto coupling
Kuramoto (1975). The coupling strength $K$ ($>0$) determines whether the
directions behave asynchronously or synchronously Acebrón  _et al._ (2005).
The magnitude of $K$ is further influenced by the spatial positions of the
particles through the term $\left|\vec{R}_{j}-\vec{R}_{i}\right|$, which is
the Euclidean distance between the $i$-th and $j$-th particle. The exponent
$\alpha$ in Eq. (2) dictates the nature of influence that is exerted by the
spatial degrees of the particles on the directional coupling. For a positive
value of $\alpha$, the effective coupling strength between two particles is
proportional to the distance, i.e., the coupling is high for far away
particles compared to the ones which are nearby. The opposite scenario occurs
for a negative $\alpha$. With $\alpha=0$, we get back the stochastic Kuramoto
dynamics without the effect of spatial positions on directions. In the
subsequent sections, we will study these two cases ($\alpha=0$ and $\alpha\neq
0$) in details. In addition, we have also assumed that the SPPs are
interacting in some thermal background. This is considered by adding a thermal
noise $\zeta(t)$ to the dynamics. The thermal noise can be understood as a
random force being generated from the random collisions of the background
particles. Here, we take this to be a uniform white noise with the following
statistical properties

$\displaystyle\langle\zeta_{i}(t)\rangle=0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\langle\zeta_{i}(t)\zeta_{j}(t^{\prime}\rangle=2D\delta_{ij}\delta(t-t^{\prime}),$
(3)

where $D$ is the noise strength. The white noise signal is generated from a
uniform distribution $\left[-\eta/2,\eta/2\right]$ so that $2D=\eta^{2}/12$.
Initially, the spatial positions of the particles are uniformly distributed
inside the square region $L\times L$, and $\theta_{i}$’s are chosen randomly
from $[0,2\pi)$.

##  III Results

In this section, we will centrally focus on the analysis of directional
synchrony among the particles. Here, the coupling topology is all-to-all, and
the density of the particles does not play a factor in achieving synchrony,
unlike the Vicsek model Vicsek  _et al._ (1995) where particles move with a
unit interaction radius. The main control parameters are the coupling strength
$K$ and the noise amplitude $\eta$. The values of these two parameters have an
opposing effect on the system: $K$ brings synchrony, whereas $\eta$ brings
randomness in the directions. Moreover, the system dynamics is dependent on
the value of $\alpha$. As mentioned earlier, the nature of influence of
spatial dynamics on directional coupling varies depending on the value of
$\alpha$ – below, we list three such cases

  1. 1.

$\alpha=0$: No influence of spatial dynamics on directional coupling. This is
the conventional Kuramoto dynamics with noise Kuramoto (1975).

  2. 2.

$\alpha>0$: Coupling strength is enhanced between faraway particles and
reduced when they are within a unit distance. Attraction is dominant in long-
range than in short-range.

  3. 3.

$\alpha<0$: On the directional degree, nearby particles are coupled with a
greater strength than the distant ones. This means attraction is dominant in
short-range than in long-range.

We will discuss these individual cases in the following subsections. We start
with the case $\alpha=0$ which can be considered as simple Kuramoto network
model without the impact of spatial positions and in presence of noise.

###  III.1 Case I: $\alpha=0$ (no impact of spatial positions)

In this case, the angular degrees of freedom follow the Kuramoto dynamics and
the positions of the particles do not affect the coupling. Eq. (2) boils down
to

$\dot{\theta}_{i}=\dfrac{K}{N}\sum_{j=1}^{N}\sin\left(\theta_{j}-\theta_{i}\right)+\zeta_{i}.$
(4)

To investigate the synchronization phenomena in the directionaility, we
introduce the following complex quantity

$re^{\tilde{i}\psi}=\dfrac{1}{N}\sum_{j=1}^{N}e^{\tilde{i}\theta_{j}},\hskip
10.0pt\tilde{i}=\sqrt{-1}$ (5)

where $r$ is the well-known order parameter (a measure of synchrony) of the
Kuramoto model Acebrón  _et al._ (2005). By definition, the value of $r$
always lie within $0$ and $1$. Clearly, $r\approx 1$ means that the directions
are oriented towards the same, i.e., they are completely synchronized. On the
other hand, $r<1$ indicates asynchronous behavior among them. Here, $\psi$ is
the mean or average of $\theta_{i}$-s. In terms of the order parameter, Eq.
(4) reads

$\dot{\theta_{i}}=Kr\sin\left(\psi-\theta_{i}\right)+\zeta_{i}.$ (6)

Equation (6) shows that a mean field $Kr$ acting on each oscillator
$\theta_{i}$. Here, we assume that $r$ is time-independent, which allows us to
find a self-consistent equation for $r$ by substituting the solution of Eq.
(6) into Eq. (5). It is known that, near the phase transition point
$r\propto(K-K_{c})^{1/2}$ and $\psi=\omega_{0}t+\psi_{0}$ Sakaguchi (1988),
which means $\psi$ rotates at a uniform frequency $\omega_{0}$. By going into
the rotating frame of reference with frequency $\omega_{0}$, the mean or
average angle of propagation can be set to zero (i.e., $\psi$ = 0), without
loss of generality. Then the modified form of the Eq. (6) becomes

$\dot{\theta_{i}}=-Kr\sin\theta_{i}+\zeta_{i}.$ (7)

Equation (7) is the well-known Langevin equation in angular space with
sinusoidal force Evans (2012). Setting $\rho(\theta,t)$ as the probability
density of a particle having direction $\theta$ at time $t$, we can write the
corresponding Fokker-Planck equation as

$\dfrac{\partial}{\partial
t}\rho(\theta,t)=-\dfrac{\partial}{\partial\theta}\left[\left(-Kr\sin\theta\right)\rho(\theta,t)\right]+D\dfrac{\partial^{2}}{\partial\theta^{2}}\rho(\theta,t),$
(8)

where note that $\rho$ is a $2\pi$ periodic function and it is normalized,
i.e., $\int_{0}^{2\pi}\rho(\theta,t)d\theta=1$. In the long time, we expect
the system to reach a stationary state defined as $\rho_{st}(\theta)$ which
satisfies the following relation

$\dfrac{d}{d\theta}\left[\left(Kr\sin\theta\right)\rho_{st}(\theta)\right]+D\dfrac{d^{2}}{d\theta^{2}}\rho_{st}(\theta)=0.$
(9)

Since there is no loss of probability, Eq. (8) can also be understood as a
continuity equation $\frac{\partial\rho}{\partial t}+\frac{\partial
J}{\partial\theta}=0$, where $J$ is the probability current or the flux.
Within this description, in the stationary state, the current takes the
following form

$J(\theta)=Kr\rho_{st}\sin\theta+D\dfrac{d}{d\theta}\rho_{st}.$ (10)

Here the probability density function in stationary state $\rho_{st}$ and the
probability density current $J$ both are periodic in nature. Imposing periodic
boundary condition on $\rho_{st}$ and $J$, we can write
$\rho_{st}(\theta)=\rho_{st}(\theta+2\pi)$ and $J(\theta)=J(\theta+2\pi)$. To
solve Eq. (9) with these boundary conditions, we consider the following
function,

$f(\theta)=\exp\left[-\int_{0}^{\theta}\dfrac{Kr\sin\phi}{D}d\phi\right].$
(11)

Skipping details (see Appendix A), we can write the probability density
function for stationary state in the following form

$\rho_{st}(\theta)=\dfrac{\exp\left(\dfrac{Kr}{D}\cos\theta\right)}{2\pi
I_{0}\left(\dfrac{Kr}{D}\right)},$ (12)

where $I_{0}$ is Bessel function of first kind of zeroth order. In the
$N\to\infty$ limit, we can compute the order parameter from the stationary
probability density function

$r=\lim_{t\to\infty}r(t)=\int_{0}^{2\pi}\exp(\tilde{i}\theta)\rho_{st}(\theta)d\theta,$
(13)

where one just needs to care about the real part of the integral Sakaguchi
(1988). Performing the integral one finds

$\displaystyle
r=I_{1}\left(\dfrac{Kr}{D}\right)\bigg{/}I_{0}\left(\dfrac{Kr}{D}\right),$
(14)

where $I_{1}$ is Bessel function of first kind of first order. Near the phase
transition point, i.e., $r\rightarrow 0$ limit, $I_{0}$ and $I_{1}$ can be
expanded in series with respect to $K_{c}r/D$, where $K_{c}$ is the critical
coupling strength that marks the continuous transition from sync state to
async state. Neglecting higher order terms, Eq. (14) can be reduced to the
following form

$K_{c}=2D.$ (15)

The above analysis provides a distinct critical coupling strength for
$\alpha=0$. Drawing insights from this, we would like to explore whether a
similar analysis can be done for the $\alpha\neq 0$ case. However, the spatial
coupling on the orientations renders a new layer of complexity. Thus, moving
forward, we apply rational approximations to the overall dynamics of the
system which allows us to make some coherent inference. Next, we discuss these
in details.

###  III.2 Case II: $\alpha\neq 0$ (under the influence of spatial positions)

The self-propelled particles roam inside the bounded region following the
dynamics (1). The spatial coordinates do not interact among them. Thus, at a
given time, the spatial coordinates inside the region remain uniformly
distributed. This might not be the case with nonzero and nonidentical
frequency distributions of the particles as one can expect chiral pattern for
large frequency values Liebchen and Levis (2017, 2022). Since the frequency of
the particles is set to zero in our model, we do not encounter this issue and
numerics also confirm that the spatial distribution remains to be uniform. To
avoid the occurrence of collision among them, we fix a lower bound on the
inter-particle distance, say $\Delta(\alpha)$. This ensures that at any moment
they can not come closer than $\Delta(\alpha)$ and thus collision among them
is safely disregarded. We employ the following condition along with Eq. (1):
Suppose at time $t$ the spatial distance between the $i$-th and $j$-th
particles is $d_{ij}(t)$. Then, we take
$d_{ij}(t)=\max\\{d_{ij}(t),\Delta(\alpha)\\}$. It should be noted that the
spatial dynamics of Eq. (1) does not contain any attraction or repulsion
forces among the particles’ positions unlike other swarming models or cell
dynamical models Fetecau, Huang, and Kolokolnikov (2011). Since the particles
move inside a bounded domain, the maximum distance between any two particles
inside a square box is the diagonal of the square region. But under the
periodic boundary condition and minimum image convention, all four vertices of
the region are identical. Thus, the maximum distance between two particles
essentially reduces to half of the diagonal length, which is $L/\sqrt{2}$.

To seek the condition for which the particles are synchronized, we rewrite Eq.
(2) as

$\dot{\theta}_{i}=\dfrac{1}{N}\sum_{j=1,j\neq
i}^{N}K_{ij}\sin\left(\theta_{j}-\theta_{i}\right)+\zeta_{i},$ (16)

where $K_{ij}=K\left|\vec{R}_{j}-\vec{R}_{i}\right|^{\alpha}$ can be
considered as time-varying weighted coupling strength between the $i$-th and
the $j$-th particle. Since the distances are uniformly distributed in the
bounded region $L\times L$, so are the $K_{ij}$’s. We can next take the
annealed approximation such that $\sum_{j\neq i}K_{ij}=\langle
K_{ij}\rangle\sum_{j\neq i}$, where $\langle K_{ij}\rangle$ is the mean of the
distributed coupling strengths. We consider $\langle K_{ij}\rangle$ as the
effective coupling strength and denote it as $K^{eff}(\alpha)$ and it depends
on the choice of $\alpha$. Equation (16) in terms of effective coupling
strength becomes,

$\dot{\theta}_{i}\approx\dfrac{K^{eff}(\alpha)}{N}\sum_{j=1,j\neq
i}^{N}\sin\left(\theta_{j}-\theta_{i}\right)+\zeta_{i}.$ (17)

In the continuum limit $N\rightarrow\infty$, $K^{eff}(\alpha)$ can be
calculated in terms of a finite integral

$K^{eff}(\alpha)=\dfrac{K}{L/\sqrt{2}-\Delta(\alpha)}\int_{\Delta(\alpha)}^{L/\sqrt{2}}x^{\alpha}dx.$
(18)

Now, depending on the value of $\alpha$, the effective coupling among
particles varies. It also has a dependence on the box size $L$. Next, we
investigate the phase transition behavior for the cases $\alpha>0$ and
$\alpha<0$, one by one.

Before delving deeper into the analysis, we can qualitatively gain some
insights of these two cases. Since we are considering a bounded squared region
of length $L$, the area under consideration is $L^{2}$. For a given particle,
the strength of coupling with others depends crucially on their distances from
that particle, in particular, to the fact that whether they lie within a unit
distance or outside. The ratio of number of particles lying inside to the same
lying outside is $\dfrac{\pi}{L^{2}-\pi}$, which is less than 1 for
$L>\sqrt{2\pi}$. For $\alpha>0$, particles which are lying outside a unit
distance are coupled with higher strength than those lying inside. This
essentially means that more number of particles are coupled with higher
coupling strengths. A converse scenario appears for $\alpha<0$. In here, the
particles lying inside a unit distance are coupled with higher strengths. So,
a comparatively lesser number of particles are coupled with larger strengths.
We infer that the particles can attain synchrony faster with $\alpha>0$ than
with $\alpha<0$ as one changes the coupling strength $K$. Keeping these
insights in mind, we now proceed to study these cases both analytically and
numerically.

####  III.2.1 $\alpha>0$ : Long-range dominance of coupling

In this case, the phase coupling strength between two particles is enhanced
when their distance is greater than unity and is reduced otherwise. In other
words, the attraction between the particles is more effective in the long-
distance than in the short-distance. Here, we can safely allow the particles
to come asymptotically close to each other without colliding, i.e.,
$\Delta(\alpha)\rightarrow 0$ since the interaction term appears only in the
numerator. Within this assumption, the integral on the right hand side of Eq.
(18) can be calculated easily and the effective coupling strength becomes

$K^{eff}(\alpha)=\dfrac{K}{\alpha+1}\Big{(}\dfrac{L}{\sqrt{2}}\Big{)}^{\alpha}.$
(19)

This phase coupling strength $K^{eff}(\alpha)$ among the particles is
independent of the inter-particle distances. Then the model reduces to one
where the particles interact among themselves with a constant coupling like in
the case of $\alpha=0$ in Eq. (4). It enables us to write down the critical
coupling strength using Eq. (15) as

$K^{eff}_{c}(\alpha)=2D,$ (20)

where $2D={\eta}^{2}/12$. Substituting Eq. (19) into Eq. (20), we further get

$\dfrac{K_{c}}{\alpha+1}\Big{(}\dfrac{L}{\sqrt{2}}\Big{)}^{\alpha}=\dfrac{\eta_{c}^{2}}{12},$
(21)

where $K_{c}$ and $\eta_{c}$ are the critical coupling strength and critical
noise amplitude, respectively. Next essential step is to check the validation
of this mean field analysis comparing with simulations. To see this, we set
$\alpha=1$ in Eq. (21) to get

$\eta_{c}^{2}=3\sqrt{2}K_{c}L.$ (22)

The presence of the component $L$ in the Eq. (22) clearly indicates the effect
of box size on the synchronization condition. We validate this fact by
plotting the order parameter $r$ as a function of $\eta$ for two different
values of $L$ with a fixed $K=0.05$. In Fig. 1, the blue and red colors
demonstrate the results for $L=10$ and $L=20$, respectively. The solid lines
are the $r$ vs. $\eta$ curves calculated numerically from the exact equations
without approximation. The blue stars and red dots are those calculated from
the model after performing approximation. When we increase $\eta$, randomness
in the directionality grows, which effectively leads to a disordered state. It
is thus evident that there is a transition from the ordered state (i.e.,
synchronization) to the disordered (asynchronous) state. We refer to
Supplementary Movies 1 & 2 which illustrate the evolution of particles in the
desynchronized and synchronized states, respectively. Subsequently, the value
of order parameter $r$ is monotonically decreasing with increasing $\eta$. We
find the critical $\eta$ values from Eq. (22) for $L=10$ and $L=20$ as
$\eta_{c}\approx 1.47$ and $\eta_{c}\approx 2.06$, respectively. These match
pretty well with our numerical curves. The small fluctuations are due to the
finite particle number since the approximation is performed in the
$N\to\infty$ limit. Fig. 1 highlights three facts: (i) the obvious dependence
of the transition curve on box size $L$; (ii) the accuracy of our
approximation; and (iii) the agreement of numerical and analytical results.

Figure 1: Order parameter $r$ as a function of $\eta$ for different box sizes
$L$ for $\alpha=1$. Variation of $r$ has been shown as a function of $\eta$
for a fixed coupling strength $K=0.05$. We perform simulations with velocity
$v_{0}=0.2$ and $N=500$ over $500$ time units with step size $0.1$ using
Heun’s method. The initial positions of the particles are uniformly
distributed inside the domain and their directions are also chosen randomly
from $[0,2\pi]$. We calculate $r$ over last $5\%$ data after discarding the
transients and it is further averaged with five different initial
configurations. Solid lines (in blue and red) represent numerically calculated
values of $r$ from Eqs. (1)-(2), for $L=10$ and $20$, respectively, which we
term as exact. Blue stars and red circles are the values of $r$ calculated
numerically from the mean-field analysis and we term them as approx. We
observe that the critical noise amplitude $\eta_{c}$ is enhanced as the box
size increases. $\eta_{c}\approx 1.47$ and $\eta_{c}\approx 2.06$ for $L=10$
and $L=20$, respectively.

In Fig. 1, we check the transition phenomena with a fixed value of $K$. Having
observed the effectiveness of our approximation and the agreement between
numerics and analytical results, we intend to check if it is also valid for
arbitrary coupling strength. For that, we move to a two-dimensional phase
space (mesh) where we vary both $\eta$ and $K$ within some given intervals,
and at each of the mesh points, we numerically calculate the value of the
order parameter $r$. Then the points are colored corresponding to their
respective $r$ values. We perform this simulation on the original model (i.e.,
without approximation) by taking $K$ from the interval $[0,0.2]$ and $\eta$
from $[0,3]$. We have done this by considering a $100\times 100$ mesh. The
result is presented in Fig. 2(a). In the color bar, yellow indicates synchrony
in directions, and blue indicates asynchrony. The synchronization and
desynchronization regions separate themselves as they are colored accordingly.
We draw the solid curve (in black) from Eq. (22) which stands for the critical
point of transition. It is quite evident that the numerical result from the
original model, Eq. (2), is well supported by the analytical findings on the
approximated model, Eq. (17).

Figure 2: Contour plot of order parameter $r$ in the $K$-$\eta$ parameter
plane for $\alpha>0$. Here, we choose $v_{0}=0.2$, $L=10$, and $N=500$. We
divide the $K$-$\eta$ plane in $100\times 100$ mesh points and calculate the
value of $r$ at each of these points. Data is generated by integrating the
Eqs. (1)-(2) using Heun’s method with a step-size 0.1 for $t=2000$ time units.
And order parameter $r$ is calculated by taking the average over the last
$5\%$ data and further averaging it over $5$ realizations. Left panel (a)
$\alpha=1$: The solid curve (in black) is the analytical prediction drawn
using Eq. (22). The region in yellow is the synchronized state whereas the
blue region is the desynchronized one. Right panel (b) $\alpha=2$: The black
curve is drawn using Eq. (21) with $\alpha=2$.

To see the robustness of our analysis, we repeat our study for a different
value of $\alpha>0$. For $\alpha=2$, the critical curve Eq. (21) becomes
$\eta_{c}^{2}=2K_{c}L^{2}$. This is represented in Fig. 2(b) by the solid
curve (in black). It acts as a separatrix between the synchronized and the
desynchronized regions. Here too, we find an excellent argument between the
numerical simulation and the mean-field analysis.

####  III.2.2 $\alpha<0$ : Short-range dominance of coupling

Contrary to the previous case, for $\alpha<0$, the directional coupling
strength between two particles is enhanced when they lie within a unit
distance. On the other hand, if the distance between them is greater than
unity, the coupling strength is reduced. The coupling strength keeps
increasing as long as the distance between particles decreases. However, we
can not allow the particles to come asymptotically close to each other as the
coupling becomes unfeasible. To avoid this scenario, we set a cut-off distance
$\Delta(\alpha)$ for the particles. We specifically choose
$\Delta(\alpha)=(0.1)^{1/|\alpha|}$ so that
$\left|\vec{R}_{j}-\vec{R}_{i}\right|^{\alpha}$ always remains greater than
$0.1$ making the system well-bounded without any ambiguity.

We again begin with the approximated expression for effective coupling
strength as proposed in Eq. (18) to study the directional synchronization. In
this case, the effective coupling strength becomes

$\displaystyle
K^{eff}=\begin{cases}\dfrac{K}{L/\sqrt{2}-\Delta(\alpha)}\ln\Big{[}\dfrac{L}{\sqrt{2}\Delta(\alpha)}\Big{]},&\alpha=-1\\\
\dfrac{K}{L/\sqrt{2}-\Delta(\alpha)}\dfrac{(L/\sqrt{2})^{\alpha+1}-\Delta^{\alpha+1}(\alpha)}{\alpha+1},&\alpha\neq-1.\end{cases}$
(23)

For $\alpha=-1$, using Eq. (23), the critical value of noise amplitude becomes

$\eta_{c}^{2}=\dfrac{12K_{c}}{L/\sqrt{2}-\Delta(-1)}\ln\Big{[}{\dfrac{L}{\sqrt{2}\Delta(-1)}}\Big{]}.$
(24)

We verify this condition by performing numerical simulations. The set-up is
similar as before. At each point of the $2$d mesh in the $K$-$\eta$ plane, the
value of $r$ is calculated. The point is color coded according to it. The
result is shown in Fig. 3(a). The points in yellow denote synchrony, whereas
blue color coded points denote asynchrony. The most visible difference from
the case $\alpha=1$ (as in Fig. 2(a)) appears as a significant reduction of
the synchronized region. This is because the coupling strength is only
dominant inside a unit distance. The solid curve (black) in Fig. 3(a)
represents the analytical critical curve, Eq. (24). This is again fairly close
to the transition region and strengthens the validity of our approximated mean
field like method. We refer to Supplementary Movies 3 & 4 which respectively
show the time evolution of particles in the disordered and ordered states when
$\alpha=-1$.

We also verify our approximation with $\alpha=-2$. After putting $\alpha=-2$
in the Eq. (23) and taking $\Delta(\alpha=-2)=0.1^{1/2}$, Eq. (23) reduces to

$K^{eff}=\dfrac{K}{L/\sqrt{2}-\sqrt{0.1}}\left(\dfrac{1}{\sqrt{0.1}}-\dfrac{\sqrt{2}}{L}\right).$
(25)

From Eqs. (20) and (25), the expression for separatrix becomes,

$\eta_{c}^{2}=\dfrac{12K_{c}}{L/\sqrt{2}-\sqrt{0.1}}\left(\dfrac{1}{\sqrt{0.1}}-\dfrac{\sqrt{2}}{L}\right),$
(26)

which is again plotted by a solid curve (black) in Fig. 3(b).

Figure 3: Contour plot of order parameter $r$ in the $K$-$\eta$ parameter
plane for $\alpha<0$. We divide the $K$-$\eta$ plane in $100\times 100$ mesh
points and calculate the value of $r$ at each of these points. Data is
generated by integrating the Eqs. (1)-(2) using Heun’s method with a step-size
0.1 for $t=2000$ time units. And order parameter $r$ is calculated by taking
the average over the last $5\%$ data and further averaging it over $5$
realizations. Bifurcation diagram for (a) $\alpha=-1$, and (b) $\alpha=-2$.
Here, compared to $\alpha=1$ case in Fig. 2, the synchronized region (in
yellow) is significantly reduced. The asynchronized region in blue prevails
over most of the $K$-$\eta$ plane. The separatrices are represented with solid
lines (in black) using Eqs. (24) and (26) respectively, in panel (a) and panel
(b). Other parameters are $v_{0}=0.2$, $L=10$, and $N=500$.

So far we have established the numerical results and validated our
approximations with them for both $\alpha>0$ and $\alpha<0$. In both cases, we
find the critical transition points from the analytical expressions. The
transition points depend on the system size $L$, we verify this numerically by
varying the system size. In Fig. 4, we draw critical coupling strength $K_{c}$
as function of $L$. The theoretical curves are indicated by the blue lines
derived from Eq. (24) (for $\alpha=-1$) and Eq. (22) (for $\alpha=1$). Red
circles are the numerically calculated critical coupling strengths using
Binder’s cumulant (see Sec. IV for details). In the next section, we study the
behavior of order parameter $r$ near the transition point by performing finite
size scaling analysis.

Figure 4: Critical coupling strength $K_{c}$ as a function of the system size
$L$. The noise strength $\eta$ is kept fixed here. Solid lines (in blue)
represent the theoretical values, whereas red circles represent the numerical
values of $K_{c}$ evaluated numerically using Binder’s cumulant. (a) The
behavior of $K_{c}$ with respect to $L$ for $\alpha=-1$ with noise strength
$\eta=1.5$. (b) Represents the behavior of $K_{c}$ with respect to $L$ for
$\alpha=1$ with noise strength $\eta=2.06$. Initially, particles are placed
uniformly randomly inside a box with size $L\times L$ and directions are
chosen randomly from $[0,2\pi]$. Here, $(t,dt,N)=(500,0.1,1000)$. The
macroscopic order parameter $r$ is calculated using last $5\%$ data and is
averaged over $5$ realizations of initial conditions.

##  IV Finite size scaling analysis

In thermodynamic limit, our system exhibits a continuous phase transition.
Near the transition point, variation of the macroscopic order parameter
becomes sensitive to the noise amplitude. To examine this, we make the
following ansatz for the rate of variation

$r\sim(\eta_{c}-\eta)^{\beta},$ (27)

where $\beta$ signifies how fast $r$ varies near the critical point. We have
already calculated the critical noise amplitude analytically in the previous
section. Eqs. (22) and (24) are the expressions of $\eta_{c}$ for $\alpha=1$
and -1, respectively. We can confirm this critical value by means of the well-
known estimator Binder’s cumulant Vollmayr  _et al._ (1993) which detects the
behavior of the phase transition. This is essentially the kurtosis of the
macroscopic order parameter given by

$U_{4}=1-\dfrac{\langle r^{4}\rangle}{3\langle r^{2}\rangle^{2}},$ (28)

where $\langle\cdots\rangle$ signifies the ensemble average of initial
configurations. For small fluctuations in $r$, we can approximate $\langle
r^{4}\rangle\approx\langle r^{2}\rangle^{2}$ and thus, $U_{4}$ lies close to
$\dfrac{2}{3}$. But near the transition point, where the fluctuation in $r$ is
large, $U_{4}$ takes a jump. The value of the tuning parameter at which
$U_{4}$ takes the minimum value, is usually considered to be the critical
point Vollmayr  _et al._ (1993). In first-order transition, the variation of
macroscopic order parameter shows a discontinuous transition with respect to
any external tuning parameter and $U_{4}$ takes a sharp jump to a very large
negative value at the transition point. But in the second order phase
transition, the variation of order parameter shows a continuous transition and
$U_{4}$ takes a jump but always remains positive. We have calculated $U_{4}$
for $\alpha=1$ with $K=0.05$ while varying $\eta$. The results are accumulated
in Fig. 5. Star indicators (in red) indicate the values of $U_{4}$ and the
solid curve (in blue) stands for the order parameter $r$ numerically
calculated from the model. We find that $U_{4}$ takes a jump around
$\eta\approx 1.47$ where $r$ also goes through a visible change from a higher
value to a lower value. The fact that $U_{4}$ remains positive near the
transition indicates a second order phase transition. Using Eq. (22), we find
$\eta_{c}=1.5$ for $\alpha=1$, $K=0.05$, and $L=10$. This is satisfactorily
close to the one we achieved by performing the $U_{4}$ analysis. By doing a
similar analysis for $\alpha=-1$, we find critical noise amplitude
$\eta_{c}=0.54$ for coupling strength $K=0.05$.

Figure 5: Finding the critical noise amplitude using Binder’s cumulant. We
simulate Eqs. (1) & (2) numerically to calculate the value of $r$ at each
$\eta$ using Heun’s method. Data is generated with a step-size 0.1 for
$t=2000$ time units. $r$ is calculated by taking the average over the last
$5\%$ data and further averaging it over $5$ realizations. Star indicators
(red) are the values of $U_{4}$, given in Eq. (28). They are joined by dashed
line to show the variation of $r$ with respect to $\eta$ for $\alpha=1$. The
scale of variation is indicated in the right $y$ axis. The lower bound value
of the jump remains positive which essentially confirms that the phase
transition is of second order in nature. Lowest value occurs at
$\eta_{c}\approx 1.47$ which closely matches with the analytical result, Eq.
(22). Here, we have taken $K=0.05$ and $L=10$. Blue curve signifies the
variation of $r$ inside the same range of $\eta$. $r$ varies between 0 to 1
which is shown in the left $y$ axis. Here, $v_{0}=0.2$ and $N=500$.

Finally, we try to find out the dependency of the noise amplitude on the
macroscopic order parameter for fixed coupling strength near the phase
transition. In Fig. 6, the slopes of the curves indicate the values of
$\beta$. As mentioned earlier, for a fixed coupling strength, the system with
exponent $\alpha=1$ can sustain synchrony for larger value of noise amplitude
than $\alpha=-1$. The curves in Fig. 6 show that the slope for the system with
$\alpha=-1$ is larger than $\alpha=1$ thus conferring to the same argument.

Figure 6: Behavior of order parameter near the transition point. The decay of
the order parameter $r$ has been examined for different $\alpha$-s but for a
fixed coupling strength. The particles are moving inside a bounded region with
box length $10$ units and the coupling strength among them is fixed at
$K=0.05$. Using the description of Binder’s cumulant, we estimate the values
of critical noise amplitude $\eta_{c}$. Panels (a) and (b) respectively
describe the phenomena for $\alpha=1$ and $-1$ and the critical noise
amplitudes are found to be $\eta_{c}\approx 1.47$ and $0.54$. By finite size
scaling, we get $\beta\approx 0.29$ for $\alpha=1$ and $\beta\approx 0.44$ for
$\alpha=-1$. For simulations, we use $v_{0}=0.2$ and $N=500$.

##  V Conclusions

In summary, we have studied the directional synchrony of SPPs inside a
confined region under periodic boundary conditions. The key input here is the
introduction of the spatial dependence among the particles in the
orientational degrees of freedom and to demonstrate how they played a crucial
role in effective coupling and in the phase transition. First we benchmarked
our model in the absence of spatial influence and found the critical condition
for synchronization to be $K_{c}=2D$. Under the impact of spatial influence,
we took annealed approximation of the distributed coupling strengths which
enabled us in calculating the effective coupling strength. We also showcased
our numerical results which are in excellent agreement with the analytical
findings. According to our simulations, in a wide range of the velocities of
order from $10^{-3}$ to $10^{-1}$, the value of $v_{0}$ does not alter the
results.

In this detailed study, we considered both the cases where coupling between
two particles is strong when they are nearby and the other way around. In one
case, the spatial influence enhances the synchronization (for $\alpha>0$) and
an opposite effect is observed otherwise (for $\alpha<0$). Earlier studies
have considered the effect of mobility of the particles on achieving
synchronization when they move randomly on one-dimensional and two-dimensional
lattices Uriu, Morishita, and Iwasa (2010). Over the years, a persistent
effort has been made to understand the effect of spatial positions and
mobility on the synchronization behavior and our work is one such attempt
towards this goal. Intuitively, our model can mimic the properties of Vicsek
model Vicsek  _et al._ (1995) for a reasonable negative value of $\alpha$.
Only the particles inside a unit vision radius are shown to be coupled with
significant coupling strength, and outside it the coupling strength becomes
negligibly small. In a nutshell, our model captures the essential
characteristics of active matter systems that follow Vicsek type interactions.

Our model, albeit general, has certain limitations that open up the stage for
future studies. For example, we assume that the particles can collide among
themselves unless some external condition is imposed on them. In our work,
this drawback was taken care of by invoking a minimal distance
$\Delta(\alpha)$ among the particles. However, collision exclusion is a
necessary condition for the study of flocking and many other natural set-ups.
Thus, attempts should be made by considering intrinsic spatial attraction and
repulsion among the particles with this inherent condition. It is convenient
to mention here that, although our mean-field approximation gives fairly
satisfactory result, there can be other methods for calculating the phase
transition curves. One such method can be to consider the local order
parameter depending on the spatial position and incorporate it in the Fokker-
Planck equation Peruani, Nicola, and Morelli (2010). A field theory approach
using tools from statistical mechanics can also be useful Großmann, Peruani,
and Bär (2016). In here, we have assumed zero frequency, but many interesting
questions arise when one also considers nonzero angular frequencies in these
models namely the chiral active matter Liebchen and Levis (2017, 2022).
Another interesting direction would be to consider other physically amenable
forms for the spatial interaction. The analytical tractability of our work
provides a testing ground for attempting these questions in near future.

## Acknowledgment

Arnab Pal gratefully acknowledges research support from the Department of
Science and Technology, India, SERB Start-up Research Grant Number
SRG/2022/000080.

## Data Availability

Data sharing is not applicable to this article as no new data were created or
analyzed in this study.

##  Appendix A Derivation of the steady state in Eq. (12)

In this section we sketch out the steps to derive Eq. (12). To this end, let
us recall the set-up where we have assumed that the orientation of a particle
diffuses in the presence of a nonlinear force. The Fokker-Planck equation for
such process reads

$\dfrac{\partial\rho(\theta,t)}{\partial
t}=-\dfrac{\partial}{\partial\theta}\left[A(\theta)\rho(\theta,t)\right]+\dfrac{1}{2}\dfrac{\partial^{2}}{\partial\theta^{2}}\left[B(\theta)\rho(\theta,t)\right].$
(29)

In the steady state, setting $\dfrac{\partial\rho(\theta,t)}{\partial t}=0$,
we find

$-\dfrac{d}{d\theta}\left[A(\theta)\rho_{st}(\theta)\right]+\dfrac{1}{2}\dfrac{d^{2}}{d\theta^{2}}\left[B(\theta)\rho_{st}(\theta)\right]=0,$
(30)

where $\rho_{st}(\theta)$ is the time independent stationary state. From the
continuity equation, one can easily identify the probability current in the
steady state which reads

$J(\theta)=A(\theta)\rho_{st}(\theta)-\dfrac{1}{2}\dfrac{d}{d\theta}\left[B(\theta)\rho_{st}(\theta)\right].$
(31)

Both the stationary state and the currents are periodic in nature so that

$\displaystyle\rho_{st}(\theta)=\rho_{st}(\theta+2\pi),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
J(\theta)=J(\theta+2\pi).$ (32)

In the stationary state, $\partial J/\partial\theta=0$ which implies $J$ is
constant throughout the space. To deal with the boundary conditions, we
consider the following function

$f(\theta)=\exp\left[2\int_{0}^{\theta}\dfrac{A(\phi)}{B(\phi)}d\phi\right].$
(33)

A first differentiation with respect to $\theta$ yields

$f^{\prime}(\theta)=2\dfrac{A(\theta)}{B(\theta)}f(\theta).$ (34)

Let us now also consider the following relation

$\dfrac{d}{d\theta}\left[\dfrac{B(\theta)\rho_{st}(\theta)}{f(\theta)}\right]=-\dfrac{2J(\theta)}{f(\theta)},$
(35)

where we have used the explicit form of the current from Eq. (31) and also
used the relation above. Integrating both sides of the above equation and
using the fact that $J(\theta)=J$ in the stationary state, we find

$\rho_{st}(\theta)=\dfrac{B(0)}{B(\theta)}\dfrac{f(\theta)}{f(0)}\rho_{st}(0)-2J\dfrac{f(\theta)}{B(\theta)}\int_{0}^{\theta}\dfrac{d\phi}{f(\phi)}.$
(36)

It is now easy to evaluate $J$ from the above equation. To this end, we use
the periodic boundary condition $\rho_{st}(0)=\rho_{st}(2\pi)$ in Eq. (36) to
get

$\rho_{st}(2\pi)=\dfrac{B(0)}{B(2\pi)}\dfrac{f(2\pi)}{f(0)}\rho_{st}(0)-2J\dfrac{f(2\pi)}{B(2\pi)}\int_{0}^{2\pi}\dfrac{d\phi}{f(\phi)},$
(37)

from which one trivially arrives at the following expression of the steady
state current

$J=\dfrac{\left[\dfrac{B(0)}{f(0)}-\dfrac{B(2\pi)}{f(2\pi)}\right]\rho_{st}(0)}{2\int_{0}^{2\pi}\dfrac{d\phi}{f(\phi)}}.$
(38)

Substituting the current from Eq. (38) into Eq. (36), one finally obtains

$\rho_{st}(\theta)=\dfrac{\dfrac{B(0)}{f(0)}\int_{\theta}^{2\pi}\dfrac{d\phi}{f(\phi)}+\dfrac{B(2\pi)}{f(2\pi)}\int_{0}^{\theta}\dfrac{d\phi}{f(\phi)}}{\dfrac{B(\theta)}{f(\theta)}\int_{0}^{2\pi}\dfrac{d\phi}{f(\phi)}}\rho_{st}(0),$
(39)

which is a general expression for the stationary state. Inspecting Eq. (9),
one now easily identifies $A(\theta)=Kr\sin\theta$ & $B(\theta)=-D/2$.
Finally, substituting all these in Eq. (33) and Eq. (39) we find the desired
form of the stationary state

$\displaystyle\rho_{st}(\theta)=\dfrac{\exp\left(\dfrac{Kr}{D}\cos\theta\right)}{2\pi
I_{0}\left(\dfrac{Kr}{D}\right)},$ (40)

which was announced in the main text.

## References

  * Vicsek  _et al._ (1995) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, “Novel type of phase transition in a system of self-driven particles,” Physical Review Letters 75, 1226–1229 (1995).
  * Carrillo  _et al._ (2010) J. A. Carrillo, M. Fornasier, G. Toscani, and F. Vecil, “Particle, kinetic, and hydrodynamic models of swarming,” Mathematical Modeling of Collective Behavior in Socio-Economic and Life Sciences , 297–336 (2010).
  * Okubo (1986) A. Okubo, “Dynamical aspects of animal grouping: swarms, schools, flocks, and herds,” Advances in Biophysics 22, 1–94 (1986).
  * Reynolds (1987) C. W. Reynolds, “Flocks, herds and schools: A distributed behavioral model,” in  _Proceedings of the 14th annual conference on Computer graphics and interactive techniques_ (1987) pp. 25–34.
  * Toner and Tu (1998) J. Toner and Y. Tu, “Flocks, herds, and schools: A quantitative theory of flocking,” Physical Review E 58, 4828 (1998).
  * Czirók  _et al._ (1996) A. Czirók, E. Ben-Jacob, I. Cohen, and T. Vicsek, “Formation of complex bacterial colonies via self-generated vortices,” Physical Review E 54, 1791 (1996).
  * Kearns (2010) D. B. Kearns, “A field guide to bacterial swarming motility,” Nature Reviews Microbiology 8, 634–644 (2010).
  * Cavagna  _et al._ (2010) A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. Santagati, F. Stefanini, and M. Viale, “Scale-free correlations in starling flocks,” Proceedings of the National Academy of Sciences 107, 11865–11870 (2010).
  * Hemelrijk and Hildenbrandt (2012) C. K. Hemelrijk and H. Hildenbrandt, “Schools of fish and flocks of birds: their shape and internal structure by self-organization,” Interface Focus 2, 726–737 (2012).
  * Huth and Wissel (1994) A. Huth and C. Wissel, “The simulation of fish schools in comparison with experimental data,” Ecological Modelling 75, 135–146 (1994).
  * Partridge (1982) B. L. Partridge, “The structure and function of fish schools,” Scientific American 246, 114–123 (1982).
  * Winfree (1967) A. T. Winfree, “Biological rhythms and the behavior of populations of coupled oscillators,” Journal of Theoretical Biology 16, 15–42 (1967).
  * Kuramoto (1975) Y. Kuramoto, “Self-entrainment of a population of coupled non-linear oscillators,” in  _International symposium on mathematical problems in theoretical physics_ (Springer, 1975) pp. 420–422.
  * Rodrigues  _et al._ (2016) F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths, “The kuramoto model in complex networks,” Physics Reports 610, 1–98 (2016).
  * Dörfler and Bullo (2011) F. Dörfler and F. Bullo, “On the critical coupling for kuramoto oscillators,” SIAM Journal on Applied Dynamical Systems 10, 1070–1099 (2011).
  * Sakaguchi and Kuramoto (1986) H. Sakaguchi and Y. Kuramoto, “A soluble active rotater model showing phase transitions via mutual entertainment,” Progress of Theoretical Physics 76, 576–581 (1986).
  * Ghosh  _et al._ (2022) D. Ghosh, M. Frasca, A. Rizzo, S. Majhi, S. Rakshit, K. Alfaro-Bittner, and S. Boccaletti, “The synchronized dynamics of time-varying networks,” Physics Reports 949, 1–63 (2022).
  * Hong and Strogatz (2011a) H. Hong and S. H. Strogatz, “Kuramoto model of coupled oscillators with positive and negative coupling parameters: An example of conformist and contrarian oscillators,” Physical Review Letters 106, 054102 (2011a).
  * Gambuzza, Frasca, and Gomez-Gardenes (2015) L. V. Gambuzza, M. Frasca, and J. Gomez-Gardenes, “Intra-layer synchronization in multiplex networks,” Europhysics Letters 110, 20010 (2015).
  * Millán, Torres, and Bianconi (2020) A. P. Millán, J. J. Torres, and G. Bianconi, “Explosive higher-order kuramoto dynamics on simplicial complexes,” Physical Review Letters 124, 218301 (2020).
  * Pikovsky  _et al._ (2001) A. Pikovsky, M. Rosenblum, J. Kurths, and A. Synchronization, “Synchronization: a universal concept in nonlinear sciences,” Self 2, 3 (2001).
  * Pecora and Carroll (1990) L. M. Pecora and T. L. Carroll, “Synchronization in chaotic systems,” Physical Review Letters 64, 821 (1990).
  * Boccaletti  _et al._ (2002) S. Boccaletti, J. Kurths, G. Osipov, D. Valladares, and C. Zhou, “The synchronization of chaotic systems,” Physics Reports 366, 1–101 (2002).
  * Anwar and Ghosh (2023) M. S. Anwar and D. Ghosh, “Synchronization in temporal simplicial complexes,” SIAM Journal on Applied Dynamical Systems 22, 2054–2081 (2023).
  * Frasca  _et al._ (2008) M. Frasca, A. Buscarino, A. Rizzo, L. Fortuna, and S. Boccaletti, “Synchronization of moving chaotic agents,” Physical Review Letters 100, 044102 (2008).
  * Fujiwara, Kurths, and Díaz-Guilera (2011) N. Fujiwara, J. Kurths, and A. Díaz-Guilera, “Synchronization in networks of mobile oscillators,” Physical Review E 83, 025101 (2011).
  * Chowdhury  _et al._ (2019) S. N. Chowdhury, S. Majhi, M. Ozer, D. Ghosh, and M. Perc, “Synchronization to extreme events in moving agents,” New Journal of Physics 21, 073048 (2019).
  * Majhi, Ghosh, and Kurths (2019) S. Majhi, D. Ghosh, and J. Kurths, “Emergence of synchronization in multiplex networks of mobile rössler oscillators,” Physical Review E 99, 012308 (2019).
  * O’Keeffe, Hong, and Strogatz (2017) K. P. O’Keeffe, H. Hong, and S. H. Strogatz, “Oscillators that sync and swarm,” Nature Communications 8, 1504 (2017).
  * Sar and Ghosh (2022) G. K. Sar and D. Ghosh, “Dynamics of swarmalators: A pedagogical review,” Europhysics Letters (2022).
  * O’Keeffe and Bettstetter (2019) K. O’Keeffe and C. Bettstetter, “A review of swarmalators and their potential in bio-inspired computing,” Micro-and Nanotechnology Sensors, Systems, and Applications XI 10982, 383–394 (2019).
  * Sar, Ghosh, and O’Keeffe (2023) G. K. Sar, D. Ghosh, and K. O’Keeffe, “Pinning in a system of swarmalators,” Physical Review E 107, 024215 (2023).
  * O’Keeffe, Ceron, and Petersen (2022) K. O’Keeffe, S. Ceron, and K. Petersen, “Collective behavior of swarmalators on a ring,” Physical Review E 105, 014211 (2022).
  * Sar  _et al._ (2022) G. K. Sar, S. N. Chowdhury, M. Perc, and D. Ghosh, “Swarmalators under competitive time-varying phase interactions,” New Journal of Physics 24, 043004 (2022).
  * Sar, O’Keeffe, and Ghosh (2023) G. K. Sar, K. O’Keeffe, and D. Ghosh, “Swarmalators on a ring with uncorrelated pinning,” Chaos: An Interdisciplinary Journal of Nonlinear Science 33, 111103 (2023).
  * Peruani, Nicola, and Morelli (2010) F. Peruani, E. M. Nicola, and L. G. Morelli, “Mobility induces global synchronization of oscillators in periodic extended systems,” New Journal of Physics 12, 093029 (2010).
  * Peruani (2017) F. Peruani, “Hydrodynamic equations for flocking models without velocity alignment,” Journal of the Physical Society of Japan 86, 101010 (2017).
  * Großmann, Peruani, and Bär (2016) R. Großmann, F. Peruani, and M. Bär, “Superdiffusion, large-scale synchronization, and topological defects,” Physical Review E 93, 040102 (2016).
  * Hong and Strogatz (2011b) H. Hong and S. H. Strogatz, “Conformists and contrarians in a kuramoto model with identical natural frequencies,” Physical Review E 84, 046202 (2011b).
  * Acebrón  _et al._ (2005) J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, “The kuramoto model: A simple paradigm for synchronization phenomena,” Reviews of Modern Physics 77, 137 (2005).
  * Sakaguchi (1988) H. Sakaguchi, “Cooperative phenomena in coupled oscillator systems under external fields,” Progress of Theoretical Physics 79, 39–46 (1988).
  * Evans (2012) L. C. Evans,  _An introduction to stochastic differential equations_ , Vol. 82 (American Mathematical Soc., 2012).
  * Liebchen and Levis (2017) B. Liebchen and D. Levis, “Collective behavior of chiral active matter: Pattern formation and enhanced flocking,” Physical Review Letters 119, 058002 (2017).
  * Liebchen and Levis (2022) B. Liebchen and D. Levis, “Chiral active matter,” Europhysics Letters 139, 67001 (2022).
  * Fetecau, Huang, and Kolokolnikov (2011) R. C. Fetecau, Y. Huang, and T. Kolokolnikov, “Swarm dynamics and equilibria for a nonlocal aggregation model,” Nonlinearity 24, 2681 (2011).
  * Vollmayr  _et al._ (1993) K. Vollmayr, J. D. Reger, M. Scheucher, and K. Binder, “Finite size effects at thermally-driven first order phase transitions: A phenomenological theory of the order parameter distribution,” Zeitschrift für Physik B Condensed Matter 91, 113–125 (1993).
  * Uriu, Morishita, and Iwasa (2010) K. Uriu, Y. Morishita, and Y. Iwasa, “Random cell movement promotes synchronization of the segmentation clock,” Proceedings of the National Academy of Sciences 107, 4979–4984 (2010).