\section{Introduction}

Active particle models have been introduced to investigate the
dynamics of individuals or objects able to convert energy
from the environment into directed motion~\cite{ramaswamy2010mechanics, marchetti2013hydrodynamics, bechinger2016active, fodor2018statistical}.
The active systems  most studied so far are
bacteria \cite{berg2008coli}, protozoa \cite{blake1974mechanics}, spermatozoa \cite{woolley2003motility}, cells, living tissues \cite{poujade2007collective}, actin filaments \cite{kohler2011structure}, active nematics \cite{guillamat2016control},
and the so-called motor-proteins \cite{sanchez2012spontaneous}.
Recently, active microswimmers have been synthesized in
labs \cite{bechinger2016active}. Typical examples are the Janus
particles \cite{takatori2016acoustic},
spherical objects containing two faces of distinct properties, such as
hydrophobicity and hydrophilicity.
Such an asymmetry produces a self-propulsion in a given direction.
The ``activation'' of a complex microswimmer could have
important applications, for instance in drug delivery:
a propulsion mechanism could enhance and direct the transport
process with a consequent increase of its efficiency.
Alternatively, some colloids could be activated by light or magnetic fields \cite{ebbens2010pursuit}.

A large interest in the polymers community resides in the synthesis and
study of star polymers \cite{LikosPRL1998} or
dendrimers \cite{madaan2014dendrimers, Mladek2008}, which have highly
branched structures and thus endowed with several interesting
characteristics such as globular, void-containing, shapes which
make them suitable for the delivery of anticancer drugs and
imaging agents.
As a consequence of the presence of cavities and channels in
their interior, the interaction between different dendrimers or
star polymers in solution could be modeled by means of
soft-core interactions, which do not prevent the overlap
between particles. In particular, it may be described by
Generalized Exponential Model (GEM) potentials whose properties
have been reviewed by Likos~\cite{Likos2001a,
BookLikosSciortinoZaccarelli}. Following this, we
adopt a coarse-grained model which replaces a suspension of
complex polymers with overlapping spherical particles.
Remarkably, under equilibrium conditions, these particles may
bind together to form clusters and these, in turn, may organize
periodically to form cluster crystals \cite{caprini2018cluster,
Delfau2015, CoslovichIkeda2013, Likos2007, Mladek2006}.
Because of the peculiarities of two-dimensional systems
\cite{strandburg1988twodimensional}, this would not be a true
crystal phase with long-range positional order, but it displays
a clear inhomogenous and periodic distribution of particles,
much more clustered and ordered than the homogenous state found
at high temperature. Such an interesting aggregation phase
occurs also in the presence of other soft-core potentials
besides the GEM ones, such as the ultra-soft core potentials
used for low-temperature bosons \cite{Cinti2014,
diaz2015monodisperse, wang2019melting} and vortices in
superconductors \cite{varney2013hierarchical}.

The self-propulsion is often modeled by means of
an effective stochastic  force. Perhaps, the simplest model
introduced in the literature to describe the behavior of some
bacteria populations is represented by the discrete
run-and-tumble particle dynamics \cite{tailleur2008statistical,
nash2010run}. Recently, suitable descriptions in terms of
continuous stochastic processes have received large attention
for the possibility of applying well-known tools of statistical
mechanics towards the development of the thermodynamic of
active microswimmers. In this framework, we mention the active
Brownian particles (ABPs)~\cite{ten2011brownian,
romanczuk2012active} and the active Ornstein-Uhlenbeck
particles (AOUPs)~\cite{szamel2014self, marconi2015towards,
marconi2016velocity, szamel2015glassy} models. Regarding the
observed phenomenology, there is a strong resemblance between
ABP and AOUP, but they also display important differences: a)
the ABP  active force has a fixed strength whereas its
orientation fluctuates in a diffusive fashion, while in the
AOUP there is no such a constraint and the propulsive force
fluctuates both in strength and direction; b) the  correlations
are of Gaussian type  in the AOUP and non-Gaussian in ABP, in
spite of the fact that they share the same two-time
self-correlation of the active force~\cite{caprini2019activity,
das2018confined}.
 So far, the practical
consequences of these differences have not been completely elucidated. One can say that in the majority of cases the two models
display a similar phenomenology, whose characteristics are:
\begin{enumerate}

\item[i)]  In the absence of external forces some single particle properties
such as the diffusion and the mean square displacement has the same form in ABP \cite{ten2011brownian, sevilla2015smoluchowski} and AOUP \cite{Gardiner}.

\item[ii)]  Both models undergo the so-called motility induced phase separation \cite{fily2012athermal, cates2015motility, gonnella2015motility, ginot2018aggregation}. While such a topic has been well studied in the ABP case, with steric interactions \cite{digregorio2018full, bialke2015active, speck2016collective, mandal2019motility}, in the presence of an attractive component of the potential \cite{redner2013reentrant, buttinoni2013dynamical} and for more complex interactions \cite{pu2017reentrant, fischer2019aggregation}, the result with AOUP model appeared only in \cite{fodor2016far}.

\item[iii)]  In the two models particles accumulate in the proximity of the walls
\cite{maggi2015multidimensional, wensink2008aggregation}. In particular,
the wall induces a non-uniform density profile decaying with a characteristic length-scale
\cite{caprini2019active, yan2015force, yan2018curved}.

\item[iv)] The dynamics in the presence of a convex, radial
    and non-harmonic potential shows particle
    accumulation far from the minimum of the potential
    \cite{caprini2019activity, takatori2016acoustic,
    dauchot2019dynamics} when the activity is large, in
    such a way that the particle distribution is not
    Boltzmann-like \cite{marconi2017heat,
    malakar2019exact}.

\end{enumerate}
The extension of the ABP and AOUP models of activity to
particles interacting with soft potentials is considered in
the present paper. We analyze how the morphological properties
of the system of particles depend on the specific modeling of
the active force and if some features of the
dynamics are model-independent.
What happens if the non-equilibrium forcing
is replaced by a colored
noise term as in the AOUP model?
As an example, it has been found with ABP driving that for some values of the activity parameter the clusters
deform into rings with an empty interior \cite{Delfau2017}, at variance with
the equilibrium situation where clusters are compact.
We will see here that such an empty-cluster crystal phase does not appear in the AOUP description.

Hereafter, we shall study in detail some properties, such as phase
diagram and cluster size, stressing similarities and
differences between the AOUP and ABP models. Moreover,
motivated by microfluidics applications, we confine the system in
large channel to evaluate the long-range influence of the
walls on the active cluster-crystal-phase.

After introducing the active ultrasoft
model in the ABP and AOUP versions,
we present a numerical
study, for high enough density, displaying
a traveling cluster-crystal aggregation phase. The role of the
active force is elucidated, determining at first the size of
the cluster and then the occurrence of unstable regions, where
clusters shrink and reform.
Then, an effective description of the
the system is developed, with the aim of describing both the
microscopic dynamics of a particle within a cluster and the
global dynamics of the pattern.
Finally, we confine the system in a
large channel to explore the interplay between the
cluster-crystal phase and the accumulation near the walls.
In the last Section we summarize the results discussing future
perspectives in the conclusive section.

\section{Model}\label{sec:model}
Polymers are often described as large complex structures
having many internal degrees of freedom, but in some cases
it is not necessary nor possible to take into account their internal properties.
For instance, when the resolving power of instruments is low some details of their structure
can be disregarded and the polymers can be assimilated to diffusing objects.
It is coherent with such an approach to represent the effective interaction among different
polymers by a pair-wise potential which depends on the coordinates of their centers of mass.
  In some cases, such as for dendrimers,
this interaction is repulsive and of the soft-core type \cite{Likos2001a}, basically due to non-diverging potentials
which do not prevent the overlap among such complex
structures.
An appropriate model for the effective interaction
is the so-called GEM-$\alpha$ potential which reads:
\begin{equation}
\label{eq:softcorepotential}
V(\{ \mathbf{x)} \} ) \propto \sum_{i<j} \phi(\bm{x}_{ij}), \qquad \phi(\bm{x}_{ij})
=\epsilon ~ e^{-(|\bm{x}_{ij}|/R)^{\alpha}} \,,
\end{equation}
being $\alpha$ a positive real number, and $\{|\bm{x}_{ij}|\}$ are
the relative distances between the pairs of particles. We study
a system of point particles which self-propel and interact with
this potential in two dimensions. Depending on the stiffness of
the potential (in particular, $\alpha>2$ is needed), the
passive system (i.e. without the active self-propulsion) shows
the occurrence of a peculiar aggregation phase at equilibrium
\cite{Likos2007, Mladek2006, Mladek2008, Delfau2015}: particles form stable clusters,
which arrange into a periodic configuration. In two dimensions
the triangular or hexagonal lattice is the only stable pattern,
occurring for small enough temperature $T$
at a large density.
The number of
particles, $N_c$, of each cluster depends on the global density
and on the typical interaction length, $R$, while the typical
inter-cluster distance is determined just by $R$. When the
particles are no longer point-like but of small finite size
\cite{caprini2018cluster, Glaser2007, Ziherl2011, Shin2009} similar properties are found, with a
new low-temperature phase -the so-called crystal
cluster-crystal phase- in which particles inside the clusters
also show an ordered structure. The cluster-formation
phenomenon despite the repulsive interaction between all the
particles can be physically interpreted by considering the
force balance between the intra-cluster repulsion (the force
felt by particles inside clusters) and the inter-cluster
effective interaction, that is, the force exerted by
neighboring clusters \cite{Delfau2015, caprini2018cluster}.
The increasing of temperature enlarges the typical size of
clusters, destroying any structure for $T$ sufficiently large.

As mentioned in the introduction, these complex structures
could be ``activated'' by chemical reactions or biological
mechanisms taking place inside or on the surface of such a
complex microswimmer, so that each individual self-propels in a
preferential direction in the same way as a simple rod-like
active particle in the absence of any structure.
The occurrence of a driving velocity could deform the structure
of the polymer, destroying its circular symmetry and
altering the shape of the effective interaction given by
Eq.\eqref{eq:softcorepotential}. But when the speed induced by
the self-propulsion is smaller than the typical velocities of
the ``microscopic'' components of the polymer (for instance,
the arms of star polymers) then the structural changes of the
polymer shape would be negligible and our description in
terms of point particles with effective interaction and
self-propulsion remains valid. We restrict in the
following to this regime.

Self-propulsion can be modeled by means of a force
vector, applied on the center of mass of each microswimmer with
dynamics
 independent on the particle position \cite{bechinger2016active, ramaswamy2010mechanics}.
In particular, we consider a two-dimensional system of $N$
interacting active particles, whose dynamics is described by
over-damped Langevin equations for the positions,
$\mathbf{x}_i$, of each particle:
\begin{equation}
\label{eq:xdynamics}
\gamma\dot{\mathbf{x}}_i= \bm{F}_i  + \sqrt{2 \gamma T}\,\bm{\eta}_i
+ \gamma\mathbf{f}_i \ ,
\end{equation}
where $\bm{F}_i=-\nabla_i V$ is the total force exerted on the
particle $i$ by the rest of the particles due to the repulsive
potential $V$ (as given by Eq.\eqref{eq:softcorepotential}). The
term $\sqrt{2\gamma T}\boldsymbol{\eta}_i$ represents the
effect of a thermal bath at temperature $T$ where active
particles are immersed, $\bm{\eta}=(\eta_x, \eta_y)$ is a
two-dimensional Gaussian white noise vector with zero average
and correlations $\langle \bm{\eta}_i
(t)\bm{\eta}_j(t')\rangle= \mathbb{I} \delta_{ij}\delta(t-t')$,
with $\mathbb{I}$ the 2d identity matrix. The constant $\gamma$
is the drag coefficient. The last term, $\mathbf{f}_i$,
models the self-propulsion mechanism of the microswimmers, and
here is where the ABP and AOUP approaches enter into the
modelling. In the ABP, $\mathbf{f}_i$ is given by a vector of
fixed norm: $\mathbf{f}_i= U_0 \mathbf{\hat{n}}_i$, where
$\mathbf{\hat{n}}_i=(\cos{\theta_i}, \sin{\theta_i})$ is a unit
vector whose angle, $\theta_i$, evolves as a Wiener process
\begin{equation}
\label{eq:ABPactiveforce}
\dot{\theta}_i = \sqrt{2D_r}\xi_i \,.
\end{equation}
The constant $D_r$ is the rotational diffusion coefficient,
while $\xi_i$ is a Gaussian white noise with zero average and
correlations $\langle \xi_i (t)\xi_j(t')\rangle=
\delta_{ij}\delta(t-t')$. Instead, in the AOUP model
$\mathbf{f}_i$ is a noise vector whose components evolve as
independent Ornstein-Uhlenbeck processes:
\begin{equation}
\label{eq:AOUPactiveforce}
\tau\dot{\mathbf{f}}_i = -\mathbf{f}_i +\sqrt{2D_a}\mathbf{w}_i \,,
\end{equation}
where $\tau$ is a correlation time, $D_a$ an effective
diffusion constant characterizing the active force, and
$\mathbf{w}_i$ is a Gaussian white noise vector with zero
averages and correlations $\langle \mathbf{w}_i
(t)\mathbf{w}_j(t')\rangle= \mathbb{I}
\delta_{ij}\delta(t-t')$.
We point out the relevance of the ratio, $D_a/\tau$, i.e. the variance of $\mathbf{f}_i$,
whose square root gives the typical average value of the active force norm.

Despite the differences between ABP and AOUP models, a
connection line between them has been already explored in
\cite{das2018confined}. We remark that AOUP has been
originally introduced as a simplification of ABP in order to
capture its phenomenology with the aim of making analytical
predictions \cite{fily2012athermal}.
Subsequently, AOUP was also used to describe the complex behavior of a passive object
immersed in a bacterial bath \cite{maggi2014generalized, maggi2017memory, steffenoni2016interacting}.
Anyway, AOUP is the simplest Gaussian model which displays
at long times the same average and two-time activity-activity correlation
function as the ABP.
In particular, this correlation in the ABP
case has an exponential form, which in two dimensions reads \cite{Farage2015}:
\begin{equation}
\label{eq:activityactivityABP}
 \langle  \mathbf{f}_i(t) \mathbf{f}_j(0) \rangle =
 U_0^2 \langle \hat{\mathbf{n}}_i(t) \hat{\mathbf{n}}_j(0) \rangle =
 \delta_{ij}\mathbb{I} \frac{U_0^2}{2} e^{-D_r t} \,.
\end{equation}
 Instead, in the AOUP model, such a correlation is \cite{Gardiner}
\begin{equation}
\label{eq:activityactivityAOUP}
 \langle  \mathbf{f}_i(t) \mathbf{f}_j(0) \rangle =
 \delta_{ij}\mathbb{I} \frac{D_a}{\tau} e^{-t/\tau} \,.
\end{equation}
Setting $\tau=1/D_r$ and $2 D_a/\tau= U_0^2$
the two-time activity-activity correlation
of the AOUP dynamics coincides with Eq.\eqref{eq:activityactivityABP} \cite{caprini2019activity},
and the average is $ \langle  \mathbf{f}_i(t) \rangle =0$ in both cases.

As commented in the introduction two main differences appear
between ABP and AOUP active forces: the fluctuating norm of
the AOUP  force \cite{das2018confined}, and the higher order
correlations \cite{sevilla2014theory} which make the ABP force
non-Gaussian. We will explore in the following if these
differences (and which) are relevant for active cluster
crystals.

\section{Numerical Results}\label{sec:numericalresults}

In this Section, we numerically explore the dynamics of
Eq.~\eqref{eq:xdynamics} for a suspension of $N$ interacting
active particles in a two-dimensional box of size $L$ with
periodic boundary conditions. We implement the Euler-Maruyama
algorithm \cite{toral2014stochastic} using both ABP and AOUP active forces, given by
Eq.\eqref{eq:ABPactiveforce} and \eqref{eq:AOUPactiveforce},
respectively. We choose a soft-core potential of the
GEM-$\alpha$ type with $\alpha=3$, which displays the
cluster-crystal aggregation phase in the passive case ($U_0=0$)
at sufficiently low temperature and large enough particle density \cite{Delfau2015}. Throughout this paper
we take $R/L=10^{-1}$, $\gamma=1$ and $\epsilon=1$. Since for
several active systems of interest  the effective diffusion
due to the active forces is much larger than the one
due to the thermal diffusion, hereafter we fix $T=(U_0^2 D_r/\gamma
)\, 10^{-4}$.

\begin{figure}[!t]
\centering
\includegraphics[clip=true,width=0.7\columnwidth,keepaspectratio]{Figure1.pdf
\caption{Snapshot of the configurations in the plane $xy$ of a system of
$N=2\times 10^3$ particles interacting with the GEM-3 potential for both
ABP (left column) and AOUP (right column) active forces.
Panels (a) and (b) are obtained with $U_0=2.75$, panels (c) and (d) with
$U_0=3.75$ and panels (e) and (f) with $U_0=4.5$.  Other parameters: $\gamma=1$,
$D_r=1$, $L=1$, $R=10^{-1}$, $\epsilon=1$, $T=10^{-4} U_0^2/D_r \gamma$. In the AOUP case the
parameters used are $D_r=1/\tau$ and $U_0^2=2 D_a/\tau$. Black lines in
panels (a) and (b) are eye-guides to show the hexagonal pattern.
Graph (c), (d), (e) and (f) display an inset around a small square region
around a cluster. In the insets, the black arrows have a length proportional to the
active force on each particle. In panel (d), isolated particles not belonging to clusters
are highlighted by red circles.
}
\label{fig:Snapshot}
\end{figure}

\begin{figure}[!t]
\centering
\includegraphics[clip=true,width=0.9\columnwidth,keepaspectratio]{Figure2.pdf
\caption{Cluster size as a function of $U_0$ for two different values of the density:
green  and blue  for $N=10^3, 2\times 10^3$, respectively, at fixed $L$.
Other parameters as in Fig. \ref{fig:Snapshot}.
Continuous lines are obtained by linear numerical fits. In the ABP case (Panel (a))
we plot the diameter, $d$, of the ring cluster, while the AOUP case (Panel (b)) we
show $2\sigma$. In the
AOUP case the label $U_0$ is a short-cut for $\sqrt{2 D_a/\tau}$. The appearance of
configurations which exchange particles among clusters is indicated by violet symbols.
Blue shades denote the occurrence of unstable clusters in the case $N=2\times10^3$,
while the same analysis for $N=10^3$ is not shown for the sake of clarity.
The left vertical dashed line indicates the
smallest value of the active force which displays the
instability of at least one cluster in the steady state.
The right vertical dashed line shows the beginning of the disordered phase.
In panel (c) we show the radial probability density,  $\rho(r)$.
Red diamonds refer to ABP, yellow triangles to  AOUP.
Black lines are the result of a numerical fit.
In panel (d) we plot the pair correlation function, $g(r)$,
for different values of the active force both for ABP and AOUP dynamics, as shown in the legend.
}
\label{fig:clustersize}
\end{figure}

In Fig.\ref{fig:Snapshot} we show long-time patterns ( $t/\tau\approx 10^2$) for different values of $U_0$
for ABP (left column) and AOUP (right column).
At such a low temperature,
both ABP and AOUP form clusters
arranged into a hexagonal pattern (see panels (a), (b), (c) and
(d)). This scenario resembles the one of passive particles
($U_0=0$) \cite{Delfau2015, caprini2018cluster}.

Notwithstanding that for low values of $U_0$ the active force
does not change the static macroscopic properties of the
pattern, i.e. the hexagonal cluster arrangement, strong
differences appear at the dynamical level. While clusters
occupy stable equilibrium positions in the passive case
($U_0=0$), without displaying any macroscopic motion, this
situation changes in the active case ($U_0 >0$). The active
forces produce a macroscopic coherent motion of the pattern,
which maintains the hexagonal cluster-crystal phase. Clusters,
stuck in the hexagonal pattern, drift persistently before
changing direction after a time which grows as $1/D_r$. In
movie 1 of the Supplementary Information, we compare the time
evolution of systems with $N=2\times 10^3$ particles for ABP
and AOUP. Despite the total active force acting on each cluster
is directed randomly in space (green arrows in movie 1),
clusters move coherently in one direction maintaining the
hexagonal arrangement. This global drift follows the average
active force of the whole system (black arrow in movie 1). The
decreasing of $D_r$ enlarges the persistence of the pattern
dynamics, for both cases, without producing any significant
change on the particles configurations.
Traveling
crystals \cite{menzel2013traveling} occur also for highly
packed suspensions of self-propelled particles interacting
by hard-core potentials, i.e. interactions diverging at $r=0$
in such a way that a finite size is attributed to the
particles. Experimental evidence of this effect has been
recently studied by means of a suspension of micro-disks
subjected to vertical vibrations
\cite{briand2018spontaneously}.
The use of density functional theory predicts such a
phenomenology and, in particular, the transition to rhombic,
quadratic, and lamellar patterns as the active force is
increased \cite{menzel2013traveling,menzel2014active}. Such
transitions do not occur in our system, where the only stable
pattern is the hexagonal one. This statement is confirmed
by movie 1 in Supplementary Information and by the study
of the pair correlation, $g(r) = \left\langle \sum_{i \neq 0}
\delta\left(x - x_i \right) \right\rangle  L^2/N$, where a target particle is at the origin, the sum runs over the other
particles, and
the brackets indicate
a circular average over positions $\mathbf{x}$ with the same
modulus $|\mathbf{x}| = r$. In Fig.\ref{fig:clustersize}d), we
compare the $g(r)$ for several values of the active force both
for ABP and AOUP dynamics, revealing the occurrence of the
typical peaks of a hexagonal pattern, indicating the
presence of first, second, third, etc. neighbors at distances
$1$, $\sqrt{3}$, $2$, $\sqrt{7}$, ... times the basic
periodicity of the pattern. This periodicity does not change
with $U_0$ and remains at the value determined by this
potential in the passive case
\cite{Delfau2015,caprini2018cluster}, namely $1.4~R\approx
0.14$. The increasing of $U_0$ leads to larger and lower peaks
towards the occurrence of the liquid shape of the $g(r)$. The
presence of the initial peak near $r\sim 0$ is due to the $N_c$
particles belonging to each cluster. We remark that the $g(r)$
remains essentially unchanged when computed at different
times during the motion (this is not shown in the figure).

The specific modeling of the active force influences the structure of a
single cluster as shown in the insets of
Fig.~\ref{fig:Snapshot}, panels (c), (d), for ABP and AOUP, respectively.
In particular, ABP's arrange in ring-like clusters
while AOUP's form Gaussian-like clusters as in the
passive case ($U_0=0$).
In Fig.\ref{fig:clustersize}c) the radial probability distribution, $\rho(r)$, of a
single cluster is numerically measured.
 $\rho(r)dr$ is the fraction
of particles of a cluster which are at a distance between $r$
and $r+dr$ from its center.
In particular, $\rho(r)=2\pi r p(\mathbf{x})$, where $p(\mathbf{x})$ is the
single-particle probability density for a specific particle to
be at $\mathbf{x}$, which turns out to depend only on the
norm $r=|\mathbf{x}|$.

In the AOUP case,
$p(\mathbf{x})$ is fitted by a Gaussian centered at
$\mathbf{x}=0$, while in ABP case $p(\mathbf{x})$ is close to a
Gaussian in $r=|\mathbf{x}|$ centered at a typical radius
larger than zero, and displays a smaller variance. For both   ABP
and AOUP  the angular location of the particles
in the clusters, i.e. the angular coordinates $\phi_i$ computed
for each particle $i$ with respect to the center of its cluster, is
approximatively equal to the orientational angle of the active
force, $\theta_i$. This feature is illustrated in the insets
 of Fig.\ref{fig:Snapshot}c) and d), where black
arrows represent the vectors of the active forces for each
particle, which point radially with respect to the center of
the cluster. Despite the apparent analogies  between the clusters obtained with the AOUP model
and the passive systems \cite{Delfau2015} (both are
non-empty), we
note the difference between these two cases: in
the large persistence regime, $D_r \ll \gamma$
an AOUP particle is stuck to its radial position
with respect to the center of its cluster only for a time $\sim
\tau$. After this correlation time, the active force is
consistently modified and, as a consequence, a large variation
in the radial coordinate of the particle occurs. In fact,
the non-existence of ring-like clusters in the AOUP
is due to the fluctuating norm of the active force at variance with the ABP where it remains constant. This fact is confirmed by simulations (not shown) where
a force with the angular time-dependence of the AOUP but
a fixed norm is used in analogy with \cite{das2018confined}.
 In this case, empty clusters are obtained.
Also, for some particular initial conditions in Eq.
(\ref{eq:AOUPactiveforce}), ring-like clusters are observed for
the standard AOUP, but only as a short-lived (for a time of the
order of $D_r^{-1}$) transient state.

In Fig. \ref{fig:clustersize} we study the size of a single
cluster, i.e. the average diameter of the ring in the ABP case
(Panel (a)) or twice its standard deviation in the AOUP case
(Panel (b)), as a function of $U_0$ ($=\sqrt{2D_a/\tau}$ in the
AOUP case). In particular, linear scaling of the standard
deviation with $U_0$ emerges for both types of active drivings
and its slope depends on the mean number of particles in each
cluster $N_c$: larger values of  $N_c$ produce larger
confinements. Such a  feature is due to the stronger repulsion
from the neighboring clusters, which is more important than the
intracluster repulsion, and needs the contribution from the
active forces to get balanced at a particular ring
radius~\cite{Delfau2017}. We remark that such size does not
depend on the persistence time of the active force, $\tau$. A
second difference between the two models emerges at large
values of $U_0$ (or of $\sqrt{2 D_a/\tau}$). As shown by the
comparison between Panels (c) and (d) of
Fig.\ref{fig:Snapshot}, at some values of $U_0$  the AOUP
reveals a different phase which does not occur in the ABP
model: In panel (d) we observe that some particles, which are
marked by red circles, leave their clusters, eventually joining
other ones, even if the cluster-crystal phase is preserved.
Quantitatively, this particle exchange is a rare event and does
not affect the cluster-size as shown in
Fig.\ref{fig:clustersize}. The number of particles which
migrate from a cluster to another, without destroying the
stable hexagonal pattern, grows with increasing $\sqrt{2
D_a/\tau}$ in the AOUP but is always a very small fraction of
$N_c$. For the ABP case, the exchange of particles could
eventually occur only for larger values of $T$ as found in
\cite{Delfau2017}, but is always independent of $U_0$.

A further increase of $U_0$ leads for both models to the
instability of some clusters. In such cases, clusters disappear
and immediately reform. The number of unstable clusters
increases as $U_0$  grows, without destroying the global order
of the pattern. While in the AOUP case clusters simply blow
up, in the ABP case we can clearly see the deformation of the
ring-like aggregates into lines after which the cluster can
reform. Then, the growth of $U_0$  produces the formation of
continuous flows of particles while other groups of particles
arrange in the hexagonal pattern (Fig.\ref{fig:Snapshot}c) and f)).
 In other words, we observe the coexistence of the
cluster-crystal phase with the disordered phase. Finally, the
formation of any stable structure for both models is prevented
by a further increase of $U_0$, a regime which is not under
investigation in the present manuscript.

The comparison between ABP and AOUP shows the
weaker stability of the cluster-crystal phase in the latter,
due to the particle-exchange
mechanism which does not occur in the ABP case. In addition,
clusters start to break down for AOUP at smaller values of $U_0$
as shown in Fig.\ref{fig:clustersize} (for $N=2\times 10^3$).
 The
disordered phase in  ABP  appears for larger values of
$U_0$ than for  AOUP. We remark that in both cases the
active force plays a role which resembles the one of an
effective temperature since its increasing induces a transition
towards the disordered phase. In any case, the role of the
active force reveals other interesting features, like the
the coherent motion of the pattern and the shape of a single
cluster, which cannot be understood in terms of an effective
temperature approach.

\section{Effective description}\label{sec:effectivedescription}

In this Section, we elaborate an effective description of the
dynamics of the system in the cluster-crystal phase by
separating the dynamics of the single particle from the
dynamics of the clusters. This has already been done for
equilibrium hard-core passive repulsive particles in
\cite{caprini2018cluster} and is based on the fact that the
typical distance between clusters remains approximately
constant. As commented before, such assumption holds also
in the presence of the active force (see the movie in
Supplementary Information). These observations allow us to
separate the effective dynamics of the single particle from the
effective dynamics of the clusters. In particular, the
effective equations for the $i$-th particle in the $j$-th
cluster reads
\begin{equation}
\label{eq:effectivedynamicsposition}
\gamma\dot{\mathbf{x}}_i^{(j)} = \boldsymbol{F}_{eff}\left( \mathbf{x}_i^{(j)} - \mathbf{R}^{(j)} \right)  +\gamma \mathbf{f}_i^{(j)} + \sqrt{2\gamma T} \boldsymbol{\eta}_i^{(j)} \,,
\end{equation}
where $\mathbf{F}_{eff}$ is the effective confining force due to
the neighboring clusters \cite{Delfau2017}.
The mean cluster
positions, $\mathbf{R}^{(j)}$, are located on a hexagonal structure
as in the purely Brownian case ($U_0=0$) since the
inter-cluster distance does not change significantly by the
presence of the active force.
The complex dynamics of interacting microswimmers is approximated by a set
 of independent particles in the presence of a grid of confining potential wells.
The validity of this approximation has been discussed in \cite{Mladek2006, Delfau2015,caprini2018cluster}
in the passive case, and basically follows by a Taylor expansion of the GEM-$\alpha$
 potential truncated at the second order, since the inter-cluster distance
is always larger than the typical cluster-size in the cluster
crystal phase. Within this approximation, particles belonging
to the same cluster are treated as independent and only
experience the effective force generated by the particles in
the neighboring clusters \cite{Delfau2017}, which is described
by the linear shape, $\mathbf{F}_{eff}\left(\mathbf{x}\right)
\approx -k \left( \mathbf{x}-\mathbf{R}^{(j)}  \right)$.
We stress that the same equations are obtained if using
interaction potentials different from GEM, the only difference
being the particular value of the spring constant $k$. At
variance with the equilibrium case where the pattern does not
move and each cluster fluctuates around its equilibrium
position, in the presence of the self-propulsion the positions
$\mathbf{R}^{(j)}$ change. However, since particle relative
distances remain constant within the cluster we shall neglect
this  movement and only consider the particle dynamics inside
each cluster and study a system of independent particles
confined in a harmonic well in the presence of the
self-propulsion.
 Taking the center of the cluster ($\mathbf{R}^{(j)}=\mathbf{0}$) as the origin of coordinates, the effective particle dynamics reads
\begin{equation}
\label{eq:xeffectivedynamics}
\gamma\dot{\mathbf{x}}= -k\bm{x}  + \sqrt{2  \gamma T}\,\bm{\eta} + \gamma\mathbf{f} \,,
\end{equation}
where we dropped the indices $i$ and
$j$ for the sake of simplicity.

In the AOUP case Eq.~\eqref{eq:xeffectivedynamics} is
linear, so that we can solve its associated  Fokker-Planck equation
for the steady state joint probability distribution function, $f(\mathbf{x},\mathbf{f})$, which reads \cite{caprini2019active}:
\begin{equation}
\label{eq:fxfa_AOUP}
f(\mathbf{x}, \mathbf{f})= \mathcal{C}
\exp{\left(- \dfrac{k}{2}\dfrac{\left|\mathbf{x}\right|^2}{D_a\gamma} \dfrac{D_a \Gamma}{D_a+\Gamma T/\gamma}  \right)}
\exp{\left(-\dfrac{\tau \Gamma}{2D_a}\left| \mathbf{f} - \dfrac{k}{\gamma} \dfrac{\Gamma D_a}{D_a+T/\gamma} \mathbf{x}  \right|^2  \right)} \,,
\end{equation}
where $\mathcal{C}$ is a normalization and $\Gamma=1+k \tau/\gamma$ is a numerical factor
 which depends on $k\tau/\gamma \gg1$, i.e. the ratio between the correlation
 time of the activity, $\tau$, and the relaxation time, $\gamma/k$, due to the harmonic potential.
The result is a multivariate Gaussian distribution with
non-zero correlations between each component of $\mathbf{x}$
and $\mathbf{f}$. In other words, a non-zero conditioned
first moment of the spatial distribution appears, so
 that $\langle \mathbf{x} \rangle \propto \mathbf{f}$, meaning that particles
prefer to spend their life far from the minimum of the potential in a fixed position determined by the value of $\mathbf{f}$.
Since the active force is an Ornstein-Uhlenbeck process, $\mathbf{f}$ can explore
 large values depending on its variance, $D_a/\tau$, even if the most probable values remains $\mathbf{f}=0$.
Instead, its persistence time, $\tau$, rules how long the particle remains close to
 the particular position determined by the value of $\mathbf{f}$.

Integrating out the active force in Eq.\eqref{eq:fxfa_AOUP}, we can easily find the probability
density for the position of a given particle, $p(\mathbf{x})$, which reads:
\begin{equation}
\label{eq:AOUPdensity}
p(\mathbf{x}) =  \mathcal{N} e^{-\dfrac{k (x^2+y^2)}{2 \gamma} \dfrac{\Gamma}{ D_a+\Gamma T /\gamma} } \,,
\end{equation}
where $\mathcal{N}$ is a normalization factor.
Formula \eqref{eq:AOUPdensity} implies that the clusters have a Gaussian
shape with  $\langle x\rangle=\langle y\rangle =0$ and defines an effective temperature of the system \cite{szamel2014self}:
\begin{equation}
\label{eq:varianceAOUPharmonic}
T_{e} =k \langle x^2 \rangle = k\langle y^2 \rangle=  D_a \gamma  \left(\dfrac{T}{\gamma D_a} + \dfrac{1}{1+k\tau/\gamma}   \right) \approx  \dfrac{D_a \gamma}{1+k\tau/\gamma} .
\end{equation}
 The last approximation holds
if $T \ll \gamma D_a$, the regime considered in this manuscript.
Note that the result for $T_e$  is in agreement with the scaling with $\tau$
recently observed in a dense suspension of active particles with hard-core
repulsive interaction, specifically Lennard Jones potentials \cite{mandal2016active, mandal2019extreme}.
We also note that Eq.\eqref{eq:varianceAOUPharmonic} approaches
 $T+\gamma D_a$, in the equilibrium limit, $\tau \to 0$, i.e. the effective temperature due
to the joint effect of self-propulsion and thermal noise for a
free microswimmer \cite{ten2011brownian, winkler2015virial, kurzthaler2016intermediate}.
The Gaussianity of the density
agrees with the shape of the clusters obtained in
Fig.~\ref{fig:clustersize}c). Moreover, the standard
deviation of the above distribution, obtained as  the square
root of Eq.~\eqref{eq:varianceAOUPharmonic}, confirms the
linear cluster-size scaling with $\sqrt{2 D_a/\tau} = U_0$,
numerically measured in Fig.~\ref{fig:clustersize} (b). We note
that the authors of ~\cite{caprini2019activity}
showed that significant
anharmonicity of the trap (for instance
$\mathbf{F}_{eff}=-\nabla U(\mathbf{x})$, with $U(\mathbf{x})\propto
|\mathbf{x}|^{2n}$ and $n\geq 2$) would lead, under AOUP dynamics, to a
``delocalization'' phenomenon, in which particles accumulate
far from the minimum of the potential,
 displaying a non-Boltzmann distribution
\cite{marconi2017heat, fodor2016far, bonilla2019active, caprini2018linear}. The absence of this
phenomenon here confirms that the effective trap potential
induced by the neighboring clusters is harmonic to a good
approximation (see also \cite{Delfau2017}).

\begin{figure}[!h]
\centering
\includegraphics[clip=true,width=0.9\columnwidth,keepaspectratio]{Figure3.pdf
\caption{A single ABP  particle in a harmonic trap. Panel (a): radial
probability density $\rho(r)$  for different values of
$D_r$. The inset shows a comparison between the Gaussian
approximation with some effective temperature, Eq.\eqref{eq:varianceAOUPharmonic},
and numerical data for $D_r=500$.
Panel (b): $\rho(r)$  for $D_r=0.3$ for three values of $T$: colored symbols from
numerical simulations, black lines from
Eq.\eqref{eq:radialprobT}. The red line is  for $T=0$.
Panel (c): angular probability distribution functions, $g(\theta-\phi)$, for two
different values of $T$ (colored points) compared with Eq.\eqref{eq:angularprobT} (black lines).
Panel (d): Variances of $g(\theta-\phi)$ vs $D_r$ for different values of $T$.
 Points are obtained by simulations while
continuum lines correspond to the variances predicted by the distribution \eqref{eq:radialprobT}.
The dashed lines indicate the limiting variances obtained for $D_r \to 0$.
Other parameters: $k=10$ and $U_0=3$.
}
\label{fig:harmonicprob}
\end{figure}

In the case of the ABP active force, in this effective dynamics description
in which interparticle forces are replaced by an external confining
potential, the solution of
the single-particle Fokker-Planck equation is not an easy task even
under the simple harmonic potential \cite{pototsky2012active, hennes2014self, dauchot2019dynamics, malakar2019exact}.
Recently, such a system has been studied in the presence of
hydrodynamic and steric interactions \cite{rana2019tuning} by means of the density functional theory \cite{hoell2019multi}.
In Fig.~\ref{fig:harmonicprob}(a), we numerically study the radial
density, $\rho(r)=2\pi r p(\mathbf{x})$, associated to
Eq.~\eqref{eq:xeffectivedynamics} with ABP self-propulsion. We
identify two regimes: i) when $D_r \gg k/\gamma$, the
activity plays the role of effective temperature and the
particle position density $p(\mathbf{x})$, has a Gaussian form
as in the case of AOUP and passive suspensions (see the inset of
Fig.\ref{fig:harmonicprob} (a)). In this regime, the increasing
of $D_r$ changes only the effective temperature, $\gamma
U_0^2/2D_r +  T$
( which is equivalent to $\gamma D_a + T$ discussed below Eq~\eqref{eq:varianceAOUPharmonic}), which determines the variance of the
distribution. It is straightforward to check that for large $D_r$ the variances of the ABP
and AOUP positions coincide. ii) When  $D_r \lesssim \gamma/k$,
particles arrange on a circular crown and the region near the
minimum of the harmonic potential becomes empty as in the
clusters in Fig.~\ref{fig:Snapshot}. The decrease of $D_r$
enhances the accumulation of particles in the proximity of $r^*\approx \gamma U_0/k$ (shown in Fig.~\ref{fig:harmonicprob}(a)), where
$r^*$ corresponds to the radius at which active and
the confining harmonic force balance.
 For
$T=0$ (or $T$ small enough with respect to $\gamma U_0^2/D_r$), the
distribution is strongly non-Gaussian and can be  approximated
by a Dirac $\delta$-function centered at $r-r^*$ in the limit
$D_r \to 0$. At $T=0$ the asymmetry of $\rho(r)$ is
quite evident (see Fig.~\ref{fig:harmonicprob}b)),
since as the norm of the
active force is fixed at $U_0$, particles cannot explore
regions with $r>r^*$. When $D_r \lesssim k/\gamma$,
particles arrange on the circular crown in such a way that
their orientational angles, $\theta$, can be approximated by
$\theta\approx \phi=\text{arctan}(y/x)$, i.e.
their angular coordinate with respect to the
minimum of the potential. In this case, the
probability distribution, $g(\theta-\phi)$, is a Gaussian which
becomes narrower as $D_r$ decreases (Fig.~\ref{fig:harmonicprob}c)). In
Fig.~\ref{fig:harmonicprob}(d) we study the variance of
$g(\theta-\phi)$ vs $D_r$ for several values of $T$, showing a
linear scaling (black triangles).

The effect of $T$ is shown in Fig.~\ref{fig:harmonicprob} b),
c) and d). On one hand, the larger thermal fluctuations
symmetrize the shape of $\rho(r)$ leading the system to explore the region $r> r^*$,
otherwise inaccessible. In this regime, $p(\mathbf{x})$ is
well-approximated by a Gaussian concentrated on a ring of radius
$|\mathbf{x}| \approx r^*$, so that the radial density is:
\begin{equation}
\label{eq:radialprobT}
\rho (r) \approx\mathcal{N} r \exp{\left(-\frac{k}{2T} \left( r-r^*\right)^2   \right)} \,,
\end{equation}
as confirmed in Fig.~\ref{fig:harmonicprob}(b). On the other
hand, even if $T$ does not change the Gaussian shape of
$g(\theta-\phi)$, its value can determine the variance of the
distribution. For large $D_r$ the variance is independent of
the value of $T$ (Fig.\ref{fig:harmonicprob}d)). A decrease of $D_r$ determines
a deviation from the linear behavior until a $T$-dependent
plateau is reached.
The shape of $g(\theta-\phi)$ is approximatively described by a Gaussian distribution:
\begin{equation}
\label{eq:angularprobT}
g(\theta-\phi)\approx \mathcal{N}\exp{\left(\frac{r^*}{2}\frac{\gamma U_0 }{T+\gamma D_r (r^*)^2 }\left(\phi - \theta  \right)^2  \right)} 	\,,
\end{equation}
whose variance excellently agrees with numerical simulations as shown in Fig.3 (d).
In the Supplementary Information
we show an analytical argument,
which comes from the analysis of the Fokker-Planck
equation, to derive Eqs.~\eqref{eq:radialprobT} and
\eqref{eq:angularprobT}.

The study of the distribution of particles inside each cluster,
within the present approximation of independent particles in a
confining potential allows to further understand the origin of
some differences between the AOUP and ABP active forces,
occurring at large $U_0$.
We already showed that
the constancy of the norm $U_0$ makes the difference. Another
way to see this is by noting that in the AOUP case the
fluctuations in the positions are  ruled by
$T_e/k$ given by Eq.\eqref{eq:varianceAOUPharmonic},
 which scales as
$\sim \gamma^2 D_a/\tau k^2$ for
large $\tau$. In this case, we can find a particle at $r>r^*$
with a finite probability, which is controlled by the strength
of the active force in the regime of very small $T$. A large
value of $D_a/\tau$ increases the diameter of the
cluster and the probability of finding a particle far away
from its most probable value. On the other hand, in the ABP
case fluctuations are mostly ruled by $T$, while $\gamma U_0/k$ only
determines the maximal cluster radius, without increasing the
particle positions fluctuations. Only  fluctuations induced
by $T$ could lead a particle to explore radial distances larger
than $\gamma U_0/k$ from the center of its cluster. This
difference explains why in the regime of negligible $T$, ABP
does not display the particle-exchange phase, at variance with
the AOUP model. This is also the main reason for which the AOUP
active force has a region of cluster instability for
smaller values of $U_0$ ($\sqrt{2 D_a/\tau}$). This description
directly agrees with the radial distribution measured in panel
Fig.\ref{fig:clustersize}c), where, despite the same value
of $U_0$, the probability of finding an AOUP-particle for $r
\gtrsim r^*$ is consistently larger than the ABP counterpart.
The particle-exchange phase could occur only if a ``small''
fraction of particles has an active force large enough to overcome
the effective barrier induced by the neighboring clusters,
which for $T\to 0$ can be provided only by the norm
fluctuations of $\mathbf{f}$, a mechanism which is absent in
the ABP model with $T=0$.
We remark that ``small'' has to be
considered in relation to the average number of particles
inside each cluster \cite{caprini2019transport}.

\section{Active cluster crystal in a channel}\label{sec:confining}

In this Section, we return to the original system of $N$ active
particles interacting through a repulsive GEM-$\alpha$
potential, and consider its behavior in a long
channel. The aim is studying the interplay between
cluster-crystal active aggregation and geometrical confinement.

The cluster crystal phase  obtained with soft-core interactions
and in the presence of a confining mechanism has been already
studied in \cite{wang2019melting} in the passive case. The hexagonal pattern is
stable if the cross-section of the channel is larger than the
typical interaction length of the soft-core potential, $L/R \gg
1$, a regime which will be considered hereafter.

In addition, it is well-known that active particles, also in
the non-interacting case, manifest the tendency to accumulate
near the walls  \cite{MinoClement2018EColi, li2009accumulation},
in the regime of large persistence. Both for AOUP and ABP active
forces, a microswimmer maintains its direction roughly during
the active force correlation time ($1/D_r$ or $\tau$).
 This persistence induces a profile in the space-density of the system,
 producing an anomalous maximum in front of each wall higher
 than the typical bulk density
\cite{maggi2015multidimensional, wittmann2016active, Elgeti2013wall, wensink2008aggregation, wagner2017steady}.
We, now, evaluate the relation
between the cluster crystal phase and the wall-accumulation.

\begin{figure}[!h]
\centering
\includegraphics[clip=true,width=0.9\columnwidth,keepaspectratio]{Figure4.pdf
\caption{Panels (a)-(f): Long-time snapshot for $N=2\times 10^3$
particles confined by  walls at positions $x=\pm L$. Panels (a), (c), (e) are
obtained with the AOUP model while panels (b), (d), (f) with the ABP one. From the
left to the right we increase $U_0=4,6,15$ ($\sqrt{2 D_a/\tau}$).
Grey regions are
forbidden due to the presence of the walls.
Panels (g)-(i):  spatial density for the AOUP case (similar for ABP -not shown-)
 along the transversal direction of the channel,
in correspondence of the snapshots.
Black dashed lines in correspondence of the peaks are guides to the eye.
Other  parameters are: $\epsilon=1$, $R=10^{-1}$, $L=1$, $k=10$ and $\gamma=D_r=1$.
}
\label{fig:Snapconf}
\end{figure}

We consider a channel of infinite length in the $y$-direction, and
width $2L$ in the transverse $x$-direction and model two parallel
walls by means of external repulsive potentials at positions $x=\pm
L$. For the sake of simplicity, we choose a truncated harmonic
wall-shape in such a way that the force exerted on the
particle, $\mathbf{F}_w$, is linear and directed along the
direction of the unit vector $\hat{\mathbf{x}}$. In detail:
\begin{equation}
\begin{aligned}
\mathbf{F}_w^{l,r}=-k(x \pm L) \theta(x \pm L)\hat{\mathbf{x}} \,,
\end{aligned}
\end{equation}
where the indices $l$ and $r$ refer to the left and the right
walls, respectively. $\theta(x)$ is the Heaviside function and the constant
$k$ represents the strength of the harmonic
potential. Along the $y$-direction we consider periodic
boundary conditions to mimic the infinite size of the channel.
In this way, the equations of motion read:
\begin{equation}
\label{eq:xdynamics_wall}
\gamma\dot{\mathbf{x}}_i= \bm{F}_i  + \sqrt{2 \gamma T}\,\bm{\eta}_i + \gamma\mathbf{f}_i +\mathbf{F}^l_w + \mathbf{F}^r_w \,.
\end{equation}

It is well-known that active particles show a peculiar behavior
in confined geometries \cite{wittmann2019pressure}. When the active force is strong and
persistent the microswimmers accumulate at the walls, so
 that their stationary probability distribution displays
anomalous peaks in the proximity of the walls, at $x\approx \pm
L$. The main reason is related to the time-persistence of their
motion, which keeps the direction of the active forces
roughly for a time $\sim 1/D_r$, depending on the
model employed \cite{lee2013active}. The number of particles accumulating at the
walls with respect to the particles in the bulk is controlled by
the persistence length of the active motion, $\lambda_a=U_0/D_r
\propto \sqrt{D_a \tau}$. When $\lambda_a \ll L$ the majority
of the particles moves freely in the bulk, while the number of
particles at the wall become comparable with the
bulk-particles in the opposite regime \cite{caprini2019transport}.
 In such a situation, the role of $T$ has important consequences: the accumulation
 is reduced and a bulk profile of the density occurs even far from the walls, both for ABP \cite{yan2015force, yan2018curved}
and AOUP \cite{caprini2019active} active forces.
Such a long-range effect of the wall has not a Brownian counterpart and increases with the ratio $T/\gamma  D_a$.

 In fig.~\ref{fig:Snapconf}, we display
configurations obtained for three different values of $U_0$
keeping fixed $D_r$. At small $U_0$, the cluster-crystal
 phase occurs and, for both ABP and AOUP, aligns
along the wall direction, as shown in panels a), b), c) and d).
Interestingly, we find two symmetric narrow stripes of clusters
attached to the walls whose shapes are strongly
deformed with respect to the bulk-clusters (see the insets in
panels a) and b)): in the ABP case, we observe
ellipsoidal-like clusters instead of the ring-like one, while
in the AOUP the Gaussian clusters are compressed in the
transversal direction with respect to the channel as if the
wall affects the $x$-variance of the particle distribution
inside the wall-clusters. This narrow stripe of clusters slides
along the walls changing direction with an average rate $\sim
1/D_r =\tau$, without breaking away from the walls. Clusters,
hexagonally arranged, form in the bulk at very small $U_0$,
behaving in the same way as unconfined clusters, with one of
the hexagonal directions aligned with the walls, as confirmed
by the study of the $x$-density, $\rho(x)$, in panel g) which
reveals separated peaks. This situation changes by increasing
$U_0$. In the bulk, far from the walls, clusters become
unstable: they disappear and reform continuously in time as in
the unconfined case, but maintain a clear alignment. As shown
in Fig.~\ref{fig:Snapconf} c) and d), the inner
clusters elongate along the $y$-direction until they collapse in
vertical stripes. In this case, the study of $\rho(x)$ cannot
capture this dynamical effect: the central peaks are only less
pronounced with respect to the lateral ones and overlap.
 A further increase of $U_0$ destroys the stripes order in
the inner region, creating a fluid-like homogenous bulk-phase,
whose structure is confirmed by the study of the density in
panel i) which is roughly homogenous until the occurrence of
the first stripes. In the proximity of the walls two stripes
clearly separated are stable: the one in front of the wall and
the one at the interface with the fluid-region, clearly
separated by an empty region. In practice, the wall stripe
creates an effective extra wall at distance $\sim \pm(L -R)$,
where active particles accumulate.

Finally, when $U_0$ is very large we recover the usual behavior
of active particles in the presence of walls (not shown), i.e. accumulation
near walls without any structures in the bulk region. We remark
that the suppression of any ordered structure occurs at larger
values of $U_0$ than in the case of the unconfined system. The presence
of the walls, which breaks the rotational symmetry, stabilizes
the stripe-phase which cannot form in the absence of walls.
As in the unconfined case, the drift of the pattern is
controlled by the average active force of the whole system.
Nevertheless, we observe that the confinement prevents any
motion along the transversal direction of the channel.

\section{Conclusion}

In this work, we have studied a system of active particles in
the presence of soft-core repulsive interactions,
a coarse-grained model for suspensions of
 activated complex polymers, such as dendrimers or star polymers.
We explored how the active force affects the cluster-crystal
phase of the system. We discover the existence of
traveling cluster crystals, with a speed induced by the active
force. The crystal moves coherently in space, for a
typical time which depends on the persistence of the active
force, maintaining its hexagonal structure. In addition, the
cluster-shape is deeply affected by the strength of the active
force which determines its size, until producing an unstable
region where at first clusters can exchange particles and then
destroy and reform continuously in time. Finally, for
large enough self-propulsion, the crystal melts. We have
checked in specific cases that such a phenomenology is not
restricted to our choice of interparticle potential and we
expect it to be present in a large class of soft-core
interactions. We explore two different modelizations of the
active force, both well-known in the literature, exploiting
analogies and differences between them. Besides some
differences in the particular parameter values at which
transitions occur, the only feature which distinguishes the
two descriptions is the cluster shape: for some values of the
control parameters they display a central hole in the ABP
case, but not in the AOUP. We have explained in detail the
reasons for this difference.

Finally, we confine the system into an infinitely long channel to
explore the dynamics of active soft repulsive particles. The
effective long-range effect of the wall clearly appears, deeply
influencing the structure of the pattern until to produce a
collapse into a stripe-phase aligned to the walls. Such a
phenomenon has not a passive counterpart and is entirely due to
the active force.

The consideration of more complicated confining geometries, finite-size for the particles,
and the role of repulsive and attractive potentials acting at different scales may give rise to interesting behaviors and applications of the interplay
of active forces and aggregated phases of particles.
This is material to be explored
in the close future.


\section{Introduction}

Given a smooth complex algebraic variety $X$ and a nonzero regular function $f\in\mathcal{O}_X(X)$, the \emph{Bernstein-Sato polynomial} $b_f(s)\in{\mathbf C}[s]$
is the monic polynomial of minimal degree such that
$$
b_f(s)f^s\in {\mathcal D}_X[s]\bullet f^{s+1}.
$$
Here ${\mathcal D}_X$ is the sheaf of differential operators on $X$ and we use $\bullet$ to denote the action of differential operators.
Note that $f^s$ can be treated as a symbol on which differential operators
act in the expected way. By making $s=-1$, we see that if $f$ is not invertible, then $b_f(s)$ is divisible by $(s+1)$, and the quotient
$\widetilde{b}_f(s)=b_f(s)/(s+1)$ is the \emph{reduced Bernstein-Sato} polynomial of $f$. The existence of $b_f(s)$ was proved by Bernstein
for the case when $X={\mathbf A}^n$ in
\cite{Bernstein} and a proof in the general case (in the analytic setting) is given in \cite{Bjork}. The Bernstein-Sato polynomial of $f$ is a subtle invariant
of the singularities of the hypersurface defined by $f$ and it is connected to several other invariants of singularities (for example, by
\cite{Malgrange}, its roots determine the
eigenvalues of the monodromy action on the cohomology of the Milnor fiber).

The above invariant has been extended to arbitrary (nonzero) coherent ideals $\mathfrak{a}$ in $\mathcal{O}_X$ in \cite{BMS}. Working locally, we may and will assume that
we have nonzero regular functions $f_1,\ldots,f_r\in\mathcal{O}_X(X)$ that generate the ideal $\mathfrak{a}$. In this case, the Bernstein-Sato polynomial $b_{\mathfrak{a}}(s)\in{\mathbf C}[s]$
is the monic polynomial of minimal degree such that
$$
b_{\mathfrak{a}}(s)f_1^{s_1}\cdots f_r^{s_r}\in \sum_{|u|=1}{\mathcal D}_X[s_1,\ldots,s_r]\bullet \prod_{u_i<0}{{s}\choose {-u_i}}f_1^{s_1+u_1}\cdots f_r^{s_r+u_r},
$$
where the sum is over all $u=(u_1,\ldots,u_r)\in{\mathbf Z}^r$ such that $|u|:=\sum_iu_i=1$. Here $s=s_1+\ldots+s_r$, where $s_1,\ldots,s_r$ are independent variables,
$f_1^{s_1}\cdots f_r^{s_r}$ is a symbol on which differential operators act in the expected way, and for every positive integer $m$, we put
${{s_i}\choose {m}}=\frac{1}{m!}\prod_{j=0}^{m-1}(s_i-j)$. The existence, independence of the choice of the generators $f_1,\ldots,f_r$, and some basic properties
of $b_{\mathfrak{a}}(s)$ were proved in \cite{BMS}. The main observation
 of this note is the following result. Given $f_1,\ldots,f_r$ as above, we consider the regular function
$g=\sum_{i=1}^rf_iy_i$ on $X\times{\mathbf A}^r$, where $y_1,\ldots,y_r$ are the coordinates on ${\mathbf A}^r$.

\begin{thm}\label{thm_description}
If $f_1,\ldots,f_r$ are nonzero regular functions on the smooth, complex algebraic variety $X$, generating the coherent ideal $\mathfrak{a}$, and if
$g=\sum_{i=1}^rf_iy_i$, then $b_{\mathfrak{a}}(s)=\widetilde{b}_g(s)$.
\end{thm}

In fact, this observation can be used to give a new proof of the existence of $b_{\mathfrak{a}}(s)$ and of its independence of the generators $f_1,\ldots,f_r$.
We hope that it will be useful for extending properties of Bernstein-Sato polynomials from the case of principal ideals to arbitrary ones.

By combining the above description of $b_{\mathfrak{a}}(s)$ with results in \cite{BMS}, we can relate invariants and properties of $g$ with those
of the ideal $\mathfrak{a}$. Recall that by a result of Kashiwara \cite{Kashiwara}, for every nonzero $f\in\mathcal{O}_X(X)$,
all roots of the Bernstein-Sato polynomial $b_f(s)$ are negative rational numbers. If $f$ is not invertible, then the negative of the largest root
of $\widetilde{b}_f(s)$ is the \emph{minimal exponent} $\widetilde{\alpha}_f$ of $f$ (with the convention that $\widetilde{\alpha}_f=\infty$ if
$\widetilde{b}_f(s)=1$, which is the case if and only if the hypersurface defined by $f$ is smooth). Therefore $\min\{1,\widetilde{\alpha}_f\}$
is the negative of the largest root of $b_f(s)$; by a result of Lichtin and Koll\'{a}r (see \cite[Theorem~10.6]{Kollar}), this is equal to the log canonical threshold
${\rm lct}(f)$ of $f$.

\begin{cor}\label{cor1}
With the notation in the theorem, we have $\widetilde{\alpha}_g={\rm lct}(\mathfrak{a})$.
\end{cor}

\begin{cor}\label{cor2}
With the notation in the theorem, if $\mathfrak{a}$ defines a reduced, complete intersection subscheme $W$, of pure codimension $r$, then
$W$ has rational singularities if and only if $\widetilde{\alpha}_g=r$ and $-r$ is a root of multiplicity $1$ of $\widetilde{b}_g(s)$.
\end{cor}

Finally, we apply the description of $b_{\mathfrak{a}}(s)$ in the theorem to show that the Strong Monodromy Conjecture for Igusa zeta functions
associated to hypersurfaces
implies the similar statement for arbitrary ideals. For the sake of simplicity, we work in the $p$-adic setting, though a similar result holds for the
motivic zeta function (see Remark~\ref{rmk_motivic} below).

Recall that if $f\in{\mathbf Z}_p[x_1,\ldots,x_n]$ is a nonzero polynomial over the ring of $p$-adic integers, the Igusa zeta function associated to $f$ is the
formal power series in $p^{-s}$ given by
$$Z_p(f;s):=\int_{{\mathbf Z}_p^n}|f(x)|_p^sd\mu_p(x),$$
where $|\cdot|_p$ is the $p$-adic absolute value on ${\mathbf Q}_p$ and $\mu_p$ is the Haar measure on ${\mathbf Q}_p^n$. This power series encodes
the numbers $a_m$ of roots of $f$ in $({\mathbf Z}/p^m{\mathbf Z})^n$ for $m\geq 1$. It was shown by Igusa \cite{Igusa1}, \cite{Igusa2} that $Z_p(f; s)$ is a rational function of $p^{-s}$,
with the candidate poles determined in terms of a log resolution of the pair $({\mathbf A}_{{\mathbf C}}^n,f)$.
The following is the outstanding open problem in this area:

\begin{conj}[Strong Monodromy Conjecture, Igusa]
Given $f\in{\mathbf Z}[x_1,\ldots,x_n]$, for every prime $p$ large enough, if
 $s_0$ is a pole of $Z_p(f;s)$, then ${\rm Re}(s_0)$ is a root of $b_f(s)$. Moreover, if the order of $s_0$ as a pole is $m$,
 then ${\rm Re}(s_0)$ is a root of $b_f(s)$ of multiplicity $\geq m$.
\end{conj}

One can study an analogue of Igusa's zeta function for arbitrary ideals $\mathfrak{a}\subseteq{\mathbf Z}_p[x_1,\ldots,x_n]$ (see \cite{Veys}).
More precisely, if $f_1,\ldots,f_r$ generate $\mathfrak{a}$, then we have a function
$\varphi_{\mathfrak{a}}\colon {\mathbf Z}_p^n\to{\mathbf Q}$ given by $\varphi_{\mathfrak{a}}(x)=\max_{i=1}^n|f_i(x)|_p$ and the corresponding Igusa zeta function
$$Z_p(\mathfrak{a};s):=\int_{{\mathbf Z}_p^n}\varphi_{\mathfrak{a}}(x)^sd\mu_p(x).$$
Again, this is a rational function of $p^{-s}$ and candidate poles can be given in terms of a log resolution of $({\mathbf A}_{{\mathbf C}}^n,\mathfrak{a})$.

\begin{thm}\label{thm_Igusa}
If $\mathfrak{a}$ is the ideal of ${\mathbf Z}_p[x_1,\ldots,x_n]$ generated by the nonzero polynomials $f_1,\ldots,f_r$
and if $g=\sum_{i=1}^rf_iy_i\in{\mathbf Z}_p[x_1,\ldots,x_n,y_1,\ldots,y_r]$, then
$$Z_p(g;s)=\frac{1-p^{-1}}{1-p^{-s-1}}Z_p(\mathfrak{a};s).$$
In particular, if $f_1,\ldots,f_r\in{\mathbf Z}[x_1,\ldots,x_n]$ and
$g$ satisfies the Strong Monodromy Conjecture, then for every prime $p$ large enough,
if $s_0$ is a pole of $Z_p(\mathfrak{a},s)$ of order $m$, then ${\rm Re}(s_0)$ is a root of $b_{\mathfrak{a}}(s)$ of multiplicity $\geq m$.
\end{thm}

In the next section we give the proof of Theorem~\ref{thm_description} and of its corollaries. The last section contains the proof of Theorem~\ref{thm_Igusa}.

\subsection*{Acknowledgement} I am indebted to Nero Budur for bringing to my attention the reference \cite{Javier} and to Wim Veys for pointing out
an inaccuracy in a previous version of this note.

\section{The description of the Bernstein-Sato polynomial of an ideal}

We begin with the formula relating the Bernstein-Sato polynomials of $\mathfrak{a}$ and $g$.

\begin{proof}[Proof of Theorem~\ref{thm_description}]
By taking an affine open cover of $X$, we see that we may assume that $X={\rm Spec}(R)$ is affine.
By definition, the Bernstein-Sato polynomial $b_g(s)$ is the monic polynomial of minimal degree such that there is $P\in \Gamma(X\times{\mathbf A}^r,
{\mathcal D}_{X\times{\mathbf A}^r})[s]$ such that
\begin{equation}\label{eq_proof1}
b_g(s)g^s=P\bullet g^{s+1}.
\end{equation}
Such $P$ can be uniquely written as $P=\sum_{\alpha,\beta\in{\mathbf Z}_{\geq 0}^r}P_{\alpha,\beta}\frac{1}{\beta !}y^{\alpha}\partial_y^{\beta}$,
with $P_{\alpha,\beta}\in\Gamma(X,{\mathcal D}_X)[s]$, only finitely many being nonzero.
Here we use the multi-index notation $y^{\alpha}=y_1^{\alpha_1}\cdots y_r^{\alpha_r}$
and $\partial_y^{\beta}=\partial_{y_1}^{\beta_1}\cdots\partial_{y_r}^{\beta_r}$ and $\beta!=\prod_{i=1}^r(\beta_i)!$
for $\alpha=(\alpha_1,\ldots,\alpha_r)$ and $\beta=(\beta_1,\ldots,\beta_r)$
in ${\mathbf Z}_{\geq 0}^r$. Furthermore, the equality in (\ref{eq_proof1}) is equivalent to
\begin{equation}\label{eq_proof2}
b_g(m)g^m=\sum_{\alpha,\beta}P_{\alpha,\beta}(m)\bullet g^{m+1}\quad\text{for all}\quad m\geq 0.
\end{equation}
Since $g=\sum_{i=1}^rf_iy_i$, we have
\begin{equation}\label{eq_proof3}
b_g(m)g^m=b_g(m)\cdot \sum_{|a|=m}{{m}\choose{a_1,\ldots,a_r}}f_1^{a_1}\cdots f_r^{a_r}y_1^{a_1}\cdots y_r^{a_r},
\end{equation}
where the sum is over all $a=(a_1,\ldots,a_r)\in{\mathbf Z}_{\geq 0}^r$ with $|a|:=a_1+\ldots+a_r=m$.
On the other hand, the right-hand side of (\ref{eq_proof2}) is equal to
$$
\sum_{\alpha,\beta}P_{\alpha,\beta}(m)\frac{1}{\beta !} y^{\alpha}\partial_y^{\beta}\bullet
\sum_{|b|=m+1}{{m+1}\choose {b_1,\ldots,b_r}}f_1^{b_1}\cdots f_r^{b_r}y_1^{b_1}\cdots
y_r^{b_r},
$$
where the second sum is over all $b=(b_1,\ldots,b_r)\in{\mathbf Z}_{\geq 0}^r$, with $|b|=m+1$.
This is further equal to
\begin{equation}\label{eq_proof4}
\sum_{\alpha,\beta}\sum_{|b|=m+1}\big(P_{\alpha,\beta}(m)\bullet f_1^{b_1}\cdots f_r^{b_r}\big)\cdot {{m+1}\choose {b_1,\ldots,b_r}}
\cdot \prod_{i=1}^r{{b_i}\choose{\beta_i}}\cdot \prod_{i=1}^ry_i^{b_i-\beta_i+\alpha_i},
\end{equation}
where we make the convention that ${{b_i}\choose{\beta_i}}=0$ if $\beta_i>b_i$.
Via the formulas in (\ref{eq_proof3}) and (\ref{eq_proof4}), the equality in (\ref{eq_proof2}) is equivalent to the fact that for every
$a=(a_1,\ldots,a_r)\in{\mathbf Z}_{\geq 0}^r$, we have
$$
b_g\big(|a|\big){{|a|}\choose {a_1,\ldots,a_r}}f_1^{a_1}\cdots f_r^{a_r}=
$$
$$
\sum_{|\beta|-|\alpha|=1}\big(P_{\alpha,\beta}\big(|a|\big)\bullet f_1^{a_1+\beta_1-\alpha_1}\cdots
f_r^{a_r+\beta_r-\alpha_r}\big)\cdot{{|a|+1}\choose {a_1+\beta_1-\alpha_1,\ldots,a_r+\beta_r-\alpha_r}}\cdot\prod_{i=1}^r{{a_i+\beta_i-\alpha_i}
\choose{\beta_i}}.
$$
An easy computation shows that this is further equivalent to
$$b_g\big(|a|\big)f_1^{a_1}\cdots f_r^{a_r}=$$
$$\big(|a|+1\big)\cdot \sum_{|\beta|-|\alpha|=1}\prod_{i=1}^r\frac{(a_i)!}{(\beta_i)!(a_i-\alpha_i)!}
P_{\alpha,\beta}\big(|a|\big)\bullet f_1^{a_1+\beta_1-\alpha_1}\cdots
f_r^{a_r+\beta_r-\alpha_r},$$
where the sum is over all $\alpha,\beta\in{\mathbf Z}_{\geq 0}^r$ with $\beta|-|\alpha|=1$ and such that $\alpha_i\leq a_i$ for all $i$.
Since it is clear that $g$ is not invertible, we know that $(s+1)$ divides $b_g(s)$, with $\widetilde{b}_g(s)=b_g(s)/(s+1)$.
It follows that $\widetilde{b}_g(s)$ is the monic polynomial of smallest degree such that we have $P_{\alpha,\beta}$ as above such that
for all $a=(a_1,\ldots,a_r)\in{\mathbf Z}_{\geq 0}^r$, we have
$$\widetilde{b}_g\big(|a|\big)f_1^{a_1}\cdots f_r^{a_r}=\sum_{|\beta|-|\alpha|=1}\prod_{i=1}^r\frac{(a_i)!}{(\beta_i)!(a_i-\alpha_i)!}
P_{\alpha,\beta}\big(|a|\big)\bullet f_1^{a_1+\beta_1-\alpha_1}\cdots
f_r^{a_r+\beta_r-\alpha_r}.$$
Equivalently, there are $P_{\alpha,\beta}\in \Gamma(X,{\mathcal D}_X)[s]$, for $\alpha,\beta\in{\mathbf Z}_{\geq 0}^r$ satisfying $|\beta|-|\alpha|=1$, with only finitely many nonzero, such that we have the equality
\begin{equation}\label{eq_proof6}
\widetilde{b}_g(s_1+\ldots+s_r)f_1^{s_1}\cdots f_r^{s_r}=\sum_{|\beta|-|\alpha|=1}\frac{\alpha !}{\beta !}\cdot \prod_{i=1}^r {{s_i}\choose {\alpha_i}}\cdot
P_{\alpha,\beta}(s_1+\ldots+s_r)\bullet f_1^{s_1+\beta_1-\alpha_1}\cdots f_r^{s_r+\beta_r-\alpha_r}.
\end{equation}
Equivalently, $\widetilde{b}_g(s)$ is the monic polynomial of minimal degree such that
$\widetilde{b}_g(s_1+\ldots+s_r)f_1^{s_1}\cdots f_r^{s_r}$ lies in
$$\sum_{|\beta|-|\alpha|=1}\prod_{i=1}^r{{s_i}\choose {\alpha_i}}{\mathcal D}_X[s_1+\ldots+s_r]\bullet f_1^{s_1+\beta_1-\alpha_1}\cdots f_r^{s_r+\beta_r-\alpha_r}.$$
This sum can be rewritten as
$$\sum_{|\gamma|=1}\sum_{\alpha}D_X[s_1+\ldots+s_r]\bullet\prod_{i=1}^r{{s_i}\choose{\alpha_i}}\cdot f_1^{s_1+\gamma_1}\cdots f_r^{s_r+\gamma_r},$$
where the first summation index runs over those $\gamma\in{\mathbf Z}_{\geq 0}^r$ such that $|\gamma|=1$ and the second summation index runs over
those $\alpha\in{\mathbf Z}_{\geq 0}^r$ such that $\alpha_i+\gamma_i\geq 0$ for all $i$. The polynomials ${{s_i}\choose{\alpha_i}}$ such that $\alpha_i+\gamma_i\geq 0$
give a basis of ${\mathbf C}[s_i]$ if $\gamma_i\geq 0$ and give a basis of ${{s_i}\choose {-\gamma_i}}\cdot{\mathbf C}[s_i]$ if $\gamma_i<0$. We thus conclude that
$\widetilde{b}_g(s)$ is the monic polynomial of smallest degree such that
$$\widetilde{b}_g(s_1+\ldots+s_r)f_1^{s_1}\cdots f_r^{s_r}\in \sum_{|\gamma|=1}{\mathcal D}_X[s_1,\ldots,s_r]\bullet
\prod_{\gamma_i<0}{{s_i}\choose {-\gamma_i}}f_1^{s_1+\gamma_1}\cdots f_r^{s_r+\gamma_r},$$
hence it is equal to the Bernstein-Sato polynomial\footnote{This is not the definition of the Bernstein-Sato polynomial
$b_{\mathfrak{a}}(s)$ in \cite{BMS}, but the definition is equivalent to this one, as explained in \cite[Section 2.10]{BMS}.} $b_{\mathfrak{a}}(s)$.
This completes the proof of the theorem.
\end{proof}

\begin{rmk}\label{rmk_existence}
Note that in the proof of Theorem~\ref{thm_description} we did not assume the existence of $b_{\mathfrak{a}}$, hence by the theorem, we can deduce
the existence of the Bernstein-Sato polynomial associated to $f_1,\ldots,f_r$ from the existence of $b_g(s)$. Furthermore, we see that
$b_{\mathfrak{a}}(s)$ only depends on the ideal generated by $f_1,\ldots,f_r$ and not on these generators. Indeed, it is enough to show that if we consider
$f_{r+1}=\sum_{i=1}^ra_if_i$ for some $a_1,\ldots,a_r\in\mathcal{O}_X(X)$ and $h=\sum_{i=1}^{r+1}f_iy_i$, then $b_g(s)=b_h(s)$.
Note that $h=\sum_{i=1}^rf_i(y_i+a_iy_{r+1})$. We have an automorphism of $X\times{\mathbf A}^{r+1}$ over $X$ which maps $y_{r+1}$ to $y_{r+1}$
and $y_i$ to $y_i+a_iy_{r+1}$ for $1\leq i\leq r$. Since this maps $g$ to $h$, it follows that $b_g(s)=b_h(s)$.
\end{rmk}

\begin{rmk}
The hypersurface $g=\sum_{i=1}^rf_iy_i$ also appeared in \cite{Javier}, where it was shown that its Milnor fibration
(at the origin) has trivial geometric monodromy and fiber homotopic to the complement of the germ defined by the ideal $(f_1,\ldots,f_r)$.
\end{rmk}

We can now deduce the first consequences of the theorem.

\begin{proof}[Proof of Corollary~\ref{cor1}]
It is shown in \cite[Theorem~2]{BMS} that the negative of the largest root of $b_{\mathfrak{a}}(s)$ is the log canonical threshold $\lct(\mathfrak{a})$ of $\mathfrak{a}$.
Since $\widetilde{\alpha}_g$ is, by definition,  the negative of the largest root of $\widetilde{b}_g(s)$, the assertion follows from Theorem~\ref{thm_description}.
\end{proof}

\begin{proof}[Proof of Corollary~\ref{cor2}]
Since $W$ is reduced and a complete intersection of pure codimension $r$, it follows from \cite[Theorem~4]{BMS} that $W$ has rational singularities if and only if $\lct(\mathfrak{a})=r$
and $-r$ is a root of multiplicity 1 of ${b}_{\mathfrak{a}}(s)$. The assertion in the corollary thus follows from Theorem~\ref{thm_description}.
\end{proof}

\section{An application to the Strong Monodromy conjecture}

For a nice introduction to Igusa's zeta function we refer to \cite{Nicaise}.
We only recall here the definition of the $p$-adic absolute value and of the Haar measure on ${\mathbf Z}_p^n$.
Let us denote by ${\rm ord}_p$ the $p$-adic valuation on ${\mathbf Q}_p$ (so that any element $u\in{\mathbf Q}_p$ can be
 written as $u=p^{{\rm ord}_p(u)}v$, with $v$ invertible in ${\mathbf Z}_p$). With this notation, if ${\rm ord}_p(u)=m$, then the $p$-adic absolute value of $u$
 is given by $|u|_p=\frac{1}{p^m}$.

The Haar measure $\mu_p$ on ${\mathbf Z}_p^n$ is the unique translation-invariant measure such that $\mu_p({\mathbf Z}_p^n)=1$. In particular, for
every $u\in{\mathbf Z}_p^n$ and every positive integer $m$, we have
$$\mu_p(u+p^m{\mathbf Z}_p^n)=\frac{1}{p^{mn}}.$$
Note also that the Haar measure is multiplicative with respect to the Cartesian product of cylinders in ${\mathbf Z}_p^n\times{\mathbf Z}_p^r\simeq{\mathbf Z}_p^{n+r}$
(recall that a cylinder in ${\mathbf Z}_p^n$ is the inverse image of some set via a projection map ${\mathbf Z}_p^n\to ({\mathbf Z}/p^m{\mathbf Z})^n$).

Given a nonzero $f\in{\mathbf Z}_p[x_1,\ldots,x_n]$, we denote by ${\rm ord}_f$ the function ${\rm ord}_p\circ f\colon {\mathbf Z}_p^n\to{\mathbf Z}_{\geq 0}$.
It then follows by definition that
\begin{equation}\label{eq1_Igusa}
Z_p(f;s)=\sum_{m\in{\mathbf Z}_{\geq 0}}\mu_p\big({\rm ord}_f^{-1}(m)\big)p^{-ms}.
\end{equation}
Similarly, if $\mathfrak{a}=(f_1,\ldots,f_r)$ is an ideal in ${\mathbf Z}_p[x_1,\ldots,x_n]$ and if we put ${\rm ord}_{\mathfrak{a}}=\min_{i=1}^r{\rm ord}_{f_i}$, then
\begin{equation}\label{eq2_Igusa}
Z_p(\mathfrak{a};s)=\sum_{m\in{\mathbf Z}_{\geq 0}}\mu_p\big({\rm ord}_{\mathfrak{a}}^{-1}(m)\big)p^{-ms}.
\end{equation}
We can now prove the main result of this section.

\begin{proof}[Proof of Theorem~\ref{thm_Igusa}]
The key point is the computation of the $p$-adic measure of ${\rm ord}_g^{-1}(m)\subseteq{\mathbf Z}_p^{n+r}$ for each $m\geq 0$.
Since $g=\sum_{i=1}^rf_iy_i$, it follows that if $(u,v_1,\ldots,v_r)\in{\mathbf Z}_p^{n+r}$ lies in ${\rm ord}_g^{-1}(m)$, then
$${\rm ord}_{\mathfrak{a}}(u)=\min_{i=1}^r{\rm ord}_{f_i}(u)\leq m.$$

Suppose now that $u\in{\mathbf Z}_p^n$ is such that $\min_{i=1}^r{\rm ord}_p(u_i)=d\leq m$. We want to describe the set $W_u(m)$ consisting of those
$v=(v_1,\ldots,v_r)\in{\mathbf Z}_p^r$ such that ${\rm ord}_p(u_1v_1+\ldots+u_rv_r)=m$.
Suppose that $j$ is such that ${\rm ord}_p(u_j)=d$.
By assumption, we can write $u_i=t^du'_i$ for $1\leq i\leq r$ and $u'_i\in{\mathbf Z}_p$, with $u'_j$ invertible.
In this case, we have ${\rm ord}_p(u_1v_1+\ldots+u_rv_r)=m$ if and only if ${\rm ord}_p(u'_1v_1+\ldots+u'_rv_r)=m-d$.
Since $u'_j$ is invertible, this means that $v_1,\ldots,\widehat{v_j},\ldots, v_r$ can be chosen arbitrarily and then the class of $v_j$
in ${\mathbf Z}/p^{m-d+1}{\mathbf Z}$ can take precisely $(p-1)$ values
(and then every lift of this class satisfies the desired condition). We thus conclude that $W_u(m)\subseteq {\mathbf Z}_p^r$ is a cylinder whose $p$-adic measure
is $\frac{p-1}{p^{m-d+1}}$.

The projection ${\mathbf Z}_p^n\times{\mathbf Z}_p^r\to{\mathbf Z}_p^n$ onto the first component induces a map
$$\tau\colon {\rm ord}_g^{-1}(m)\to\bigsqcup_{d=0}^m {\rm ord}_{\mathfrak{a}}^{-1}(d).$$
If we decompose each ${\rm ord}_{\mathfrak{a}}^{-1}(d)$ as a disjoint union of cylinders such that on each of these cyclinders
$\min_i{\rm ord}_{f_i}$ is achieved by some fixed $i$, then for every such cylinder $C\subseteq {\rm ord}_{\mathfrak{a}}^{-1}(d)$, the subset
$\tau^{-1}(C)\subseteq {\mathbf Z}_p^n\times{\mathbf Z}_p^r$ is a cylinder with
$$\mu_p\big(\tau^{-1}(C)\big)=\mu_p(C)\cdot \frac{p-1}{p^{m-d+1}}.$$
Therefore we have
$$\mu_p\big({\rm ord}_g^{-1}(m)\big)=\sum_{d=0}^m\mu_p\big({\rm ord}_{\mathfrak{a}}^{-1}(d)\big)\cdot\frac{p-1}{p^{m-d+1}}.$$
Using the formulas (\ref{eq1_Igusa}) and (\ref{eq2_Igusa}), we obtain
$$Z_p(g;s)=\sum_{m\geq 0}\frac{1}{p^{ms}}\cdot \sum_{d=0}^m\mu_p\big({\rm ord}_{\mathfrak{a}}^{-1}(d)\big)\cdot\frac{p-1}{p^{m-d+1}}$$
$$=\frac{p-1}{p}\cdot \sum_{d\geq 0}\frac{\mu_p\big({\rm ord}_{\mathfrak{a}}^{-1}(d)\big)}{p^{ds}}\cdot\sum_{m\geq d}\frac{1}{p^{(m-d)(s+1)}}
=\frac{1-p^{-1}}{1-p^{-(s+1)}}Z_p(\mathfrak{a};s).$$
This gives the first assertion in the theorem.

The formula relating $Z_p(g;s)$ and $Z_p(\mathfrak{a};s)$ shows that if we denote by
$n_p(g; \lambda)$ and $n_p(\mathfrak{a}; \lambda)$ the order of $\lambda$ as a pole of $Z_p(g; s)$ and $Z_p(\mathfrak{a}; s)$, respectively, then
$n_p(g; \lambda)=n_p(\mathfrak{a}; \lambda)$ for $\lambda\neq -1$; moreover,
if $n_p(\mathfrak{a};-1)\geq 1$, then $n_p(g;-1)=n_p(\mathfrak{a};-1)+1$.
The second assertion in the theorem follows from this
and Theorem~\ref{thm_description}.
\end{proof}

\begin{rmk}\label{rmk_motivic}
For the sake of simplicity, we assumed in Theorem~\ref{thm_Igusa} that $\mathfrak{a}$ is an ideal in ${\mathbf Z}_p[x_1,\ldots,x_n]$. A similar formula holds, with the same proof,
if we assume that $f\in O_K[x_1,\ldots,x_n]$, where $O_K$ is the ring of integers of a $p$-adic field $K$. Moreover, the proof generalizes immediately to the case of the motivic zeta functions of Denef and Loeser \cite{DL}. In this case, we see that if $X$ is a smooth complex algebraic variety, $\mathfrak{a}$ is the coherent ideal
generated by $f_1,\ldots,f_r\in\mathcal{O}_X(X)$, and $g=\sum_{i=1}^rf_iy_i$, then the motivic zeta functions $Z_{\rm mot}(g;s)$ and $Z_{\rm mot}(\mathfrak{a}; s)$
of $g$ and $\mathfrak{a}$, respectively, are related by the following formula
$$Z_{\rm mot}(g; s)=\frac{1-{\mathbf L}^{-1}}{1-{\mathbf L}^{-(s+1)}}Z_{\rm mot}(\mathfrak{a}; s).$$
\end{rmk}

\section*{References}
\begin{biblist}

\bib{Bernstein}{article}{
   author={Bern\v{s}te\u{i}n, I. N.},
   title={Modules over a ring of differential operators. An investigation of
   the fundamental solutions of equations with constant coefficients},
   journal={Funkcional. Anal. i Prilo\v{z}en.},
   volume={5},
   date={1971},
   number={2},
   pages={1--16},
}

\bib{Bjork}{book}{
       author={Bj{\"o}rk, J.},
       title={Analytic ${\mathcal D}$-modules and applications},
       series={Mathematics and its Applications},
       publisher={Kluwer Academic Publishers},
       date={1993},
}

\bib{BMS}{article}{
   author={Budur, N.},
   author={Musta\c{t}\u{a}, M.},
   author={Saito, M.},
   title={Bernstein-Sato polynomials of arbitrary varieties},
   journal={Compos. Math.},
   volume={142},
   date={2006},
   number={3},
   pages={779--797},
}

\bib{DL}{article}{
   author={Denef, J.},
   author={Loeser, F.},
   title={Motivic Igusa zeta functions},
   journal={J. Algebraic Geom.},
   volume={7},
   date={1998},
   number={3},
   pages={505--537},
}

\bib{Javier}{article}{
   author={Fern\'{a}ndez de Bobadilla, J.},
   title={On homotopy types of complements of analytic sets and Milnor
   fibres},
   conference={
      title={Topology of algebraic varieties and singularities},
   },
   book={
      series={Contemp. Math.},
      volume={538},
      publisher={Amer. Math. Soc., Providence, RI},
   },
   date={2011},
   pages={363--367},
}

\bib{Igusa1}{article}{
   author={Igusa, J.},
   title={Complex powers and asymptotic expansions. I. Functions of certain
   types},
   note={Collection of articles dedicated to Helmut Hasse on his
   seventy-fifth birthday, II},
   journal={J. Reine Angew. Math.},
   volume={268/269},
   date={1974},
   pages={110--130},
}

\bib{Igusa2}{article}{
   author={Igusa, J.},
   title={Complex powers and asymptotic expansions. II. Asymptotic
   expansions},
   journal={J. Reine Angew. Math.},
   volume={278/279},
   date={1975},
   pages={307--321},
}

\bib{Kashiwara}{article}{
   author={Kashiwara, M.},
   title={$B$-functions and holonomic systems. Rationality of roots of
   $B$-functions},
   journal={Invent. Math.},
   volume={38},
   date={1976/77},
   number={1},
   pages={33--53},
}

\bib{Kollar}{article}{
   author={Koll\'ar, J.},
   title={Singularities of pairs},
   conference={
      title={Algebraic geometry---Santa Cruz 1995},
   },
   book={
      series={Proc. Sympos. Pure Math.},
      volume={62},
      publisher={Amer. Math. Soc., Providence, RI},
   },
   date={1997},
   pages={221--287},
}

\bib{Malgrange}{article}{
  author= {Malgrange, B.},
     title= {Polynomes de {B}ernstein-{S}ato et cohomologie \'evanescente},
 booktitle= {Analysis and topology on singular spaces, {II}, {III}
              ({L}uminy, 1981)},
    series = {Ast\'erisque},
    volume = {101},
    pages = {243--267},
 publisher = {Soc. Math. France, Paris},
      date = {1983},
      }

\bib{Nicaise}{article}{
   author={Nicaise, J.},
   title={An introduction to $p$-adic and motivic zeta functions and the
   monodromy conjecture},
   conference={
      title={Algebraic and analytic aspects of zeta functions and
      $L$-functions},
   },
   book={
      series={MSJ Mem.},
      volume={21},
      publisher={Math. Soc. Japan, Tokyo},
   },
   date={2010},
   pages={141--166},
}

\bib{Veys}{article}{
   author={Veys, W.},
   author={Z\'{u}\~{n}iga-Galindo, W. A.},
   title={Zeta functions for analytic mappings, log-principalization of
   ideals, and Newton polyhedra},
   journal={Trans. Amer. Math. Soc.},
   volume={360},
   date={2008},
   number={4},
   pages={2205--2227},
}

\end{biblist}

\end{document}


\section{SUPPLEMENTAL MATERIAL}

\subsection{NOTE 1. Resonance field of the free Py(8) layer as a function of temperature}

\begin{figure}[h]
\includegraphics[width=9 cm]{FigSM1.eps}\\
\textbf{Figure A1.} (a) Resonance field of the free Py(8) layer as a function of temperature. (a) In-plane angular profile of the resonance field ($H_\mathrm{r}$) shows the same weak uniaxial magnetic anisotropy at all temperatures. The slight increase in the magnitude of $H_\mathrm{r}$ is due to the thermal demagnetization (decrease in the effective magnetization, $M_\mathrm{eff}$) of the free Py layer, shown in panel (b). The values of $M_\mathrm{eff}$ at the given measurement temperatures were obtained from fitting the $H_\mathrm{r}$-vs-$\varphi_\mathbf{H}$ profiles [solid lines in panel (a)] using the standard FMR phenomenology [J. Smit and H. G. Beljers, Phillips Res. Rep. 10, 113 (1955)].
\label{fig_A1}
\end{figure}

\newpage
\subsection{NOTE 2. Resonance field and FMR linewidth versus temperature}

\begin{figure}[h]
\includegraphics[width=10 cm]{FigSM2.eps}\\
\textbf{Figure A2.} Resonance field and FMR linewidth versus temperature for the structure with a nonmagnetic seed layer (with FeMn/Ta bilayer). (a) In-plane angular profiles of the resonance field of the free Py layer show a similar behavior to that for the structure with a magnetic seed layer [with FeMn/Py bilayer; see Fig.~A1(a)]. (b) FMR linewidth ($\Delta H$) shows a much less pronounced angular dependence at all temperatures, which is in contrast to the significant temperature induced changes in $\Delta H$ for the structure with FeMn/Py [Figs.~\ref{fig_2}(c),~\ref{fig_3}(a)].
\label{fig_A2}
\end{figure}

\newpage
\subsection{NOTE 3. Rotation of the N\'eel vector of the AF}

\begin{figure}[h]
\includegraphics[width=10 cm]{FigSM3.eps}\\
\textbf{Figure A3.} Rotation of the N\'eel vector of the AF extracted from the angular dependences of the resonance field using the model fit detailed in Note~1. (a) In-plane angular dependences of the resonance field for the free Py layer ($H_\mathrm{r}^\mathrm{free}$) and the pinned Py layer ($H_\mathrm{r}^\mathrm{pin}$) of the Py(8)/Cu(5)/FeMn(5)/Py(5) structure, measured at 280~K. The lines are the fits performed using the model of Note~4. (b) Corresponding angular dependences of the FMR linewidth of the free Py layer ($\Delta H_\mathrm{free}$) and (c) simulated misalignment of the N\'eel vector ($\varphi_\mathbf{N}^\mathrm{sim}$) relative to the exchange-pinning direction ($\varphi_\mathbf{H} = 0\degree$) set at fabrication in the FeMn/Py bilayer. The dependence in (c) was calculated using the parameters obtained from the model fit shown in (a). The misalignment $\varphi_\mathbf{N}^\mathrm{sim}$ is linear in $\varphi_\mathbf{H}$ in the angle interval between the two maxima in $\Delta H_\mathrm{free}$-vs-$\varphi_\mathbf{H}$, shown in (b).
\label{fig_A3}
\end{figure}

\newpage
\subsection{NOTE 4. Gibbs free energy for F/AF bilayers}

To model the deviation of the N\'eel vector of the F/AF bilayer from its easy (exchange-pinning) direction, we split the energy term commonly used for describing exchange bias into two terms, corresponding to the interlayer exchange coupling at the F/AF interface and the magnetic anisotropy in AF. The Gibbs free energy per unit volume becomes

\noindent where $J_\mathrm{F}$ and $J_\mathrm{AF}$ are the constants of the uniaxial magnetic anisotropy in the pinned ferromagnetic Py and the antiferromagnetic FeMn layers, respectively; $J_\mathrm{ex}$ is the constant of exchange coupling at the Py/FeMn interface; $Z$ and $G_\mathrm{dem}$ are the Zeeman and demagnetization energies; $\mathbf{m}_\mathrm{p}$ is the unit vector along magnetization $\mathbf{M}_\mathrm{p}$ of the pinned Py.

\begin{equation}\tag{A1}
\begin{aligned}
G_\mathrm{p} = Z~+ &~G_\mathrm{F} + G_\mathrm{AF} + G_\mathrm{ex} + G_\mathrm{dem} = \\
& = -\mathbf{M}_\mathrm{p} \mathbf{H} - J_\mathrm{F}(\mathbf{x}\cdot\mathbf{m}_\mathrm{p})^2 - J_\mathrm{AF}(\mathbf{x}\cdot\mathbf{N})^2 - J_\mathrm{ex}(\mathbf{m}_\mathrm{p}\cdot\mathbf{N}) + 2\pi M_\mathrm{p}^2 (\mathbf{z}\cdot\mathbf{m}_\mathrm{p})^2.
\end{aligned}
\end{equation}

When external magnetic field $\mathbf{H}$ is oriented in the film plane, terms $G_\mathrm{F}$, $G_\mathrm{AF}$, and $G_\mathrm{ex}$ can be transformed as follows: $G_\mathrm{F} = -0.5 M_\mathrm{p}H_\mathrm{ua}\cos 2 \varphi_\mathrm{p}$, $G_\mathrm{AF} = -J_\mathrm{AF} \cos \varphi_\mathbf{N}$ and $G_\mathrm{ex} = -M_\mathrm{p}H_\mathrm{ex}\cos (\varphi_\mathrm{p}-\varphi_\mathbf{N})$, where angles $\varphi_\mathrm{p}$ and $\varphi_\mathbf{N}$ are defined in the inset to Fig.~A3(a). Minimization of (A1) gives

\begin{equation}\tag{A2}
\begin{cases}
\dfrac{J_\mathrm{AF}}{J_\mathrm{ex}} \sin( 2\varphi_\mathbf{N}) = \sin (\varphi_\mathrm{p}-\varphi_\mathbf{N}), \\[3mm]
\dfrac{Z}{J_\mathrm{ex}} \sin (\varphi_\mathrm{p}-\varphi_\mathbf{H}) = -\sin (\varphi_\mathrm{p}-\varphi_\mathbf{N}),
\end{cases}
\end{equation}

\noindent where $Z = H M_\mathrm{p}$.

Equations (A1) and (A2) can be used for fitting the angular dependence of the FMR spectra of the pinned Py layer in the general case of non-rigid $\mathbf{N}$. Curve $H_\mathrm{r}^\mathrm{pin}(\varphi_\mathbf{H})$ in Fig.~A3(a) is an example of such fitting, showing a good agreement between the fit and the experimental data.

\end{document}