text:
\section{Introduction}

The standard atomistic approach to tribology, i.e. the study on friction
dissipation and wear, usually focuses on indestructible ``crushproof''
particles, typically atoms or molecules \cite{Muser03, Muser07,
  VanossiRMP13, Vigentini14, Manini15, Manini16, Manini17, Apostoli17,
  Panizon18}.
Recent works \cite{Bohlein12, Hasnain13, Hasnain14, Brazda17, Brazda18}
have brought colloids into the realm of tribology, by letting repulsive
hard-core particles interact with a periodic "corrugation potential",
generated by means of optical forces, and driven by viscous drag.
The conditions of that experiment match the single indestructible particle
paradigm.
Colloids and, more generally, soft matter systems would however allow one
to investigate situations where complex objects, carrying an internal
structure, under the stress produced by external interactions and driving,
can alter their internal structure both in shape and even in the number of
the component sub-particles.
Thus one can expect the presence of a varied set of regimes of response of
the system to the external perturbation.
In the present work, we address precisely this last situation by means of
molecular-dynamics (MD) simulations of a system forming micro phases.
Micro phase formation can take place under a number of conditions on the
interparticle interactions and on the state variables in three dimensions
\cite{Sears99, Sciortino04, Mladek06, Archer07, Archer08}
as well as in two
\cite{Seul95, Imperio04, Imperio06, Diaz-Mendez15}.
The microphases can be disordered
\cite{Sears99, Archer07, Imperio04, Imperio06}
as well as ordered in a crystalline-like state or other patterned
structures like stripes and lamellae
\cite{Leibler80, Fredrickson89, Imperio04, Imperio06, Mladek06, Diaz-Mendez15}
or even as an ordered bicontinous state
\cite{Pini15}
in which both components of a two-component system span the space.
Examples of micro-phase forming systems are block-polymers
\cite{Leibler80, Fredrickson89}
or hard colloids with competing interactions, e.g. short range
attraction-long range repulsion like for some colloid-polymer mixtures
\cite{Seul95, Sears99, Archer07, Imperio04, Sciortino04, Imperio06}.
Certain soft-matter systems like star polymers or dendrimers in a good
solvent can interpenetrate each other to a large extent so that the
effective interaction between the centers of mass of two such entities has
a soft-core character.
Also such systems can form microphases in which the colloids spontaneously
aggregate into clusters, and the clusters are organized in a crystalline
state \cite{LoVerso06, Mladek08}.
Rheological properties of cluster fluid phase \cite{Imperio08} and of
cluster crystal phase \cite{Nikoubashman11} have been investigated but no
study has yet addressed the tribology of such systems.

In the present work, we simulate a system of interacting colloidal
particles, whose mutual interaction potential generates the {\em
  spontaneous} formation of clusters.
In particular, we concentrate on a 2D geometry with repulsive
interactions, like in many experiments carried out with colloids
\cite{Mangold03, Mikhael08, Pertsinidis08, Mikhael10, Bohlein12,
  Bohlein12PRL}, but the concepts introduced here may be relevant for
vortexes in superconductors \cite{Meng17} too.
The novelty compared to these experiments is the adoption of an
interparticle interaction which supports the spontaneous formation of
cluster phases.
Cluster dynamics has been investigated in the past \cite{Reichhardt12,
  Reichhardt09b, Reichhardt09c}: in these works clustering was forced by an
external potential, while in the present work clustering is the intrinsic
result of the particle-particle interaction, and is retained even when any
external potentials are turned off.

Section~\ref{model:sec} introduces and motivates the model.
As a first step, in Sec.~\ref{gs:sec} we characterize the $T=0$ equilibrium
states of the model free from external interactions.
We study the ground-state energy of several regular arrangements of
clusters formed by $n$ particles (with $n$ from 1 to 6, that we will refer
to as $n$-clusters) as a function of the mean colloids density, obtaining a
$T=0$ phase diagram.
Then Sect.~\ref{driven:sec} reports the investigation of the dynamics of
the resulting arrays of clusters interacting with a lattice-matched
external periodic potential and an homogeneous pulling force.
Specifically, we report the mobility curves for these cluster systems, as a
function of the intensity of the pulling force, for a few values of the
amplitude of the corrugation potential.
We characterize in detail the different sliding regimes, with entire
clusters advancing, or only parts of them.
The main outcome of this research is that the clusters internal dynamics is
indeed affecting significantly the depinning force and the overall
mobility.

\section{The Model}\label{model:sec}

\subsection{The Interparticle Interaction}\label{hamilt:sec}

\begin{figure}
\centerline{
\includegraphics[width=0.45 \textwidth,clip=]{phi_cc.pdf}
}
\caption{\label{inter_pot:fig}
  The adopted colloid-colloid potential $\phi_{\rm cc}(r)$,
  Eq.~\eqref{Inter_particle}, consisting on a soft-core function plus a
  weak short-range hard-core term that we introduced to prevent the
  collapse of clusters at $T=0$.
  Inset: detail of the short-distance plateau.
}
\end{figure}

Spontaneous cluster formation can be the result of competitive interactions
\cite{Seul95, Sears99, Archer07, Imperio04, Sciortino04, Imperio06}, but
clusters can also form due to soft potentials
\cite{Diaz-Mendez15,LoVerso06, Mladek08}.
Here we adopt an interaction of the latter type.
We consider the following form for the pairwise interaction potential
energy:
\begin{equation}\label{Inter_particle}
  \phi_{\rm cc}(r) = \frac{A_0}{r^6 + R_c^6} + \frac{B_0}{r^6}
\,.
\end{equation}
This interaction, depicted in Fig.~\ref{inter_pot:fig}, consists of two
terms: a soft-core repulsive interaction, and a weak hard-core repulsion.
The $B_0=0$ (purely soft-core) model has been fully characterized in two
dimensions in Ref.~\cite{Diaz-Mendez15}: its phase diagram consists in a
high-temperature low-density regular fluid phase, and a low-temperature
high-density cluster phase.
The reason for the cluster phase is that when the inter-particle distance
decreases under $r\lesssim R_c$, the potential energy $\phi_{\rm cc}(r)$
flattens out to a nearly constant energy, with the result of producing a
quite small repulsive force.
As a result, in one dimension a particle in between two other particles
kept at fixed positions, when the distance between these fixed particles is
not too large, may be better off energetically by coming closer to either
of them, rather than remaining at the middle point.
This suggests the mechanism whereby, in any dimensionality and at large
enough density, a spontaneous symmetry breaking leading to clustering can
be energetically favorable against a uniform fluid phase.
A detailed mathematical criterion has been developed for the formation of
cluster phases in the case of soft-core potentials in terms of their
Fourier transform becoming negative in a range of finite wave vector
\cite{Likos01}.
At high $T$, entropic effects and fluctuations tend to favor the uniform
fluid phase against cluster formation.
In the opposite $T\to 0$ extreme, for $B_0=0$ no mechanism keeps the
colloids apart, and so the clustering tendency would reach the extreme
limit of collapsed point-like clusters with all colloids of a cluster
sitting exactly at the same position.
To prevent this singular behavior we add the hard-core $B_0$ term of
Eq.~\eqref{Inter_particle} for whose coefficient we adopt a relatively
small value $B_0 = 5\times10^{-5}\,A_0$ \cite{Consonni16} in order to
perturb the phase diagram as little as possible, at least in the
not-too-high-density region that we are interested in.

\begin{table}
  \begin{center}
    \begin {tabular}{l|l}
      \hline
      \hline
      Physical quantity & System units \\
      \hline
      Length & $R_c$ \\
      Energy & $E_0=A_0/R_c^6$ \\
      Mass & $m$ \\
      \hline
      Number density & $R_c^{-2}$ \\
      Force & $F_0=A_0/R_c^7$ \\
      Time & $t_0 = R_c^4\sqrt{m/A_0}$ \\
      Velocity & $v_0=R_c^{-3}\sqrt{A_0/m}$ \\
      Viscous damping & $\eta_0 = R_c^{-4}\sqrt{A_0m}$ \\
      Mobility & $\mu_0 = R_c^4/\sqrt{A_0m}$ \\
      \hline
      \hline
    \end{tabular}
  \end{center}
  \caption{\label{units:tab}
    The system of units adopted in the simulations.
    Every physical quantity is expressed in terms of these natural model units.
  }
\end{table}

We take the repulsive-potential characteristic distance $R_c$, and
characteristic energy $A_0/R_c^6$, and the mass $m$ of the colloidal
particles as fundamental units.
In the following we express all physical quantities in terms of suitable
model units, as listed in Table~\ref{units:tab}.

\subsection{The Equations of Motion}\label{eq:sec}

In simulations, we let the entire system evolve following the standard
Langevin dynamics in two dimensions \cite{Gardiner}, under the effect of an
external periodic corrugation potential $U_{\rm ext}(\mathbf r)$ and a
constant driving force $F$ applied to each particle.
The equation of motion for the $j^{\rm th}$ particle is:
\begin{eqnarray} \label{langevin:eq}
  m\ddot{\mathbf r}_j &=&
  F \, \hat{\mathbf x} - \eta \dot{\mathbf r}_j +
  \\\nonumber
 && - {\mathbf \nabla}_{{\mathbf r}_j} \left[ \sum\limits_{k \ne j}^N
    \phi_{\rm cc}(|{\mathbf r}_k-{\mathbf r}_j|) + U_{\rm ext}(\mathbf
    r_j) \right] + {\mathbf \xi}_j \,.
\end{eqnarray}
Here $\eta$ is the coefficient of viscous friction associated to the
fluid where the colloids move; we generate an overdamped dynamics by
adopting a large $\eta = 28\, \eta_0$.
The ${\mathbf\xi}_j$ terms are Gaussian-distributed random forces with
amplitude $\sigma = \sqrt{2m\eta k_{\rm B}T/\Delta t}$, where $\Delta t$
is the simulation time step: they simulate the collisions of the
molecules in the fluid with the colloids, thus generating the
appropriate Brownian dynamics.
Together, the viscous and the random-force terms in Eq.~\eqref{langevin:eq}
represent a standard Langevin thermostat \cite{Gardiner}.
In simulations, the inter-particle potential $\phi_{\rm cc}$ is cutoffed
smoothly (vanishing potential and potential derivative) at a distance
$R_{\rm cutoff} = 5\,R_c$ \cite{Consonni16}.

\begin{figure}
\centerline{
\includegraphics[width=0.45 \textwidth,clip=]{w_3_blue.pdf}
}
\caption{\label{ext_potential:fig}
  The profile $W(r_x,r_y)$ of the corrugation potential which tends to
  localize the colloids.
  It is assumed to take a periodic hexagonal-lattice form, with spacing
  $a_{\rm pot}$, see Eq.~\eqref{Wfunc:eq}.
}
\end{figure}

The periodic external potential $U_{\rm ext}(\mathbf r)$ simulates the
effects of friction against a crystalline surface.
In the laboratory, this interaction has been realized by means of a
modulated light field constructed by laser interference, so that its
spacing and intensity that can be tuned with considerable freedom
\cite{Bohlein12}.
In our simulations, we assume a corrugation potential of the form:
\begin{equation}
 U_{\rm  ext}({\mathbf r}) = V_0\,W({\mathbf r})
 \,.
\end{equation}
Here $V_0$ is the amplitude of the spatial modulation and $W(\vec{r})$ is a
2D periodic function with hexagonal symmetry, period $a_{\rm pot}$, and
unit amplitude.
Its explicit form is:
\begin{eqnarray}\nonumber
  W({\mathbf r}) &=&
  -\frac{2}{9} \left[ \frac{3}{2}
    + 2\cos\left(\frac{2\pi r_x}{a_{\rm pot}}\right)
    \cos\left(\frac{2\pi  r_y}{\sqrt{3}a_{\rm pot}}\right) \right.
    \\ \label{Wfunc:eq} & &
    \left.
    + \cos\left(\frac{4\pi r_y}{\sqrt{3}a_{\rm pot}}\right)
    \right]
  .
\end{eqnarray}
This function is illustrated in Fig.~\ref{ext_potential:fig}.

Following the $r_y = 0$ path in the $x-$direction a particle goes
through alternating minima and saddle points of $W({\mathbf r})$.
This is the energetically less costly trajectory that a colloid can
follow when pushed by an $x$-directed force.
Along this path, the external potential is a simple sinusoidal
oscillation
\begin{equation}
  U_{\rm ext}(r_x,0) =
  -V_0\,\left[
    \frac{5}{9} + \frac{4}{9}\cos\left(\frac{2\pi r_x}{a_{\rm pot}}\right)
    \right]
  \,.
\end{equation}
The resulting energy barrier separating the saddle points from the
minima is therefore $\frac{8}{9} \,V_0$.
In this situation the static friction of an isolated particle at $T =
0$, i.e.\ the minimum force needed to push that particle across the
potential corrugation when inertial effects are negligible, equals the
largest value of $\frac{\partial}{\partial r_x}U_{\rm ext}(r_x, 0)$,
namely:
\begin{equation}\label{F1s}
  F_{1s} = \frac{8\pi V_0}{9a_{\rm pot}}
\,.
\end{equation}
In the simulations and in the relative graphs we express all applied
forces by comparison to the static friction $F_{1s}$ of the single
colloid, expressed in Eq.~\eqref{F1s}.
For the lattice spacing we take $a_{\rm pot} = 1.5\, R_c$, and select
values of the particles density that favor arrays of clusters whose
spacing matches precisely this separation, as discussed in
Sect.~\ref{cell:sec}.

\subsection{The Cell and the Initial Condition}\label{cell:sec}

We simulate the 2D model in a parallelogram-shaped supercell with
periodic boundary conditions (PBC), thus fixing the average density.
To respect the hexagonal symmetry that the arrays repulsive colloids (or
clusters thereof) tend to adopt, we use a supercell generated by two
primitive vectors of equal length $L$ and forming an angle of $60^\circ$
with each other.
For $V_0=0$ we are free to adopt an arbitrary cell size: in the static
simulations of Sect.~\ref{gs:sec} we will vary $L$ freely in order to
evaluate the energetics as a function of the colloid density, and therefore
determine the $T=0$ phase diagram of the free model.

For $V_0\neq 0$, this supercell needs to match the symmetry of the
corrugation potential: we must stick to $L$ given by an integer multiple of
$a_{\rm pot}$.
We settle for a compromise size $L=12\,a_{\rm pot}$, small enough to
guarantee simulations involving a not-too-large number of particles (and
therefore manageable simulation times), but large enough for a fair
averaging over fluctuations and allowing for independent and possibly
asynchronous movements of individual clusters.
Accordingly, the cell contains $144$ corrugation minima in an area
$A=\sqrt{3} L^2/2$.

As discussed below, we find that a cluster-cluster spacing $1.5 R_c$ as
optimal for favoring a well-defined $n$-cluster phase.
In the present work we focus on a lattice-matched condition between the
cluster-cluster spacing and the spacing $a_{\rm pot}$ between adjacent
minima of the corrugation.
Therefore in all dynamical simulations of Sect.~\ref{driven:sec} we
settle for fixing $a_{\rm pot}=1.5\, R_c$, and therefore a supercell
side $L=18\,R_c$, thus generating a supercell area $A\simeq 280.6\,
R_c^2$.
In this supercell, the potential has $144$ wells, each initially hosting
a single particle or a $n = 2$ to $6$-particles cluster.
We simulate a total $N=144\, n$ particles, which thus sample the
following discrete values for the average density:
\begin{equation} \label{rhoval:eq}
  \rho = \frac NA =
  n\,\frac 2{\sqrt{3} \, a_{\rm pot}^2} =
  n\, \frac 8 { 9 \sqrt{3} \, R_c^2}
\,.
\end{equation}

For the initial condition, we start off with distances $\delta = 0.32\,
R_c$ between the colloids within each cluster suggested by the inflection
of the inter-particle potential, Fig.~\ref{inter_pot:fig}.
This initial configuration relaxes very quickly to the equilibrium
interparticle distances that depends on the state of the system.
We start with a regular lattice of preformed clusters, because the
alternative possibilities of starting, e.g., with the particles
regularly or randomly distributed across the supercell leads to
irregular patterns involving different cluster sizes, which may describe
appropriately large temperature and/or phase-coexistence conditions, but
surely do not represent the optimal single-type cluster phases at $T=0$.

\section{Ground-state Energetics for the Free Model} \label{gs:sec}

\begin{figure}
\centerline{
\includegraphics[width=0.4 \textwidth,clip=]{cluster2_orient_gs.pdf}
}
\caption{\label{cluster2_orient_gs:fig}
  Small portions of fully relaxed $T=0$ configurations representing
  examples of alternative local minima for the total repulsive energy of
  the free model ($V_0 = 0$, $F=0$) in a $n=2$ cluster phase at density
  $\rho = 16 /( 9 \sqrt{3} \, R_c^2) \simeq 1.0264\,R_c^{-2}$.
  The repulsive energy per particle of these two configurations is (a)
  $E = 1.108486 \, E_0$; (b) $E =1.108907\, E_0$.
}
\end{figure}

While at finite temperature the number of colloids in each cluster
fluctuates, at $T=0$, for every given density $\rho$, one expects that
an hexagonal lattice of clusters with the same number $n$ of particles
represents a local minimum of the total repulsive energy.
We explore the density range from $\rho= 0.4\,R_{c}^{-2}$ to
$4\,R_{c}^{-2}$.
In the absence of corrugation potential ($V_0=0$) and driving force
($F=0$), for every $n$-cluster and every density (i.e. every lattice
spacing of the cluster crystal), we execute a full relaxation by means
of a damped dynamics, for at least a total time $2000\, t_0$, but
extending the relaxation for longer time until the residual kinetic
energy is safely below $10^{-8}\, E_0$.
This method leads us very close to the precise ground state.
Depending on $n$ and $\rho$, the system may remain trapped in a local
minimum characterized by some residual energy associated to the
cluster's orientational order, see Fig.~\ref{cluster2_orient_gs:fig}.
We verified that these orientational-order energies are quite small, of
the order of $10^{-4} E_0$ per particle.
For this reason and for the irrelevance (discussed below) to the driven
dynamics, we do not investigate the details of the orientational order
in the ground state of the clustered phases.

\begin{figure}
\centerline{
\includegraphics[width=0.45 \textwidth,clip=]{ground_state.pdf}
}
\caption{\label{ground_state:fig}
  The static ($T=0$, $F=0$) free ($V_0=0$) model total colloid-colloid
  repulsive energy per particle, as a function of the number density
  $\rho$ of colloids.
  Individual curves compare different $n$-cluster configurations,
  $n=1,\dots 6$.
  The orientational order is irrelevant here, as orientational-energy
  differences are well within the size of the reported points, see
  Fig.~\ref{cluster2_orient_gs:fig}.
  Arrows identify the densities adopted for the driven dynamics studied in
  Sect.~\ref{driven:sec}.
}
\end{figure}

By running several energy minimizations starting from different cluster
phases, we can follow the evolution of the ground-state energy as a
function of the density $\rho$.
For each configuration we evaluate the total colloid-colloid repulsive
energy, divide it by the number of colloids in the supercell, and obtain
the energy per particle.
Figure~\ref{ground_state:fig} displays the curves of the total energy
per particle for each $n=1$ to $6$, as a function of $\rho$.
We see that at each given density, several cluster configurations
represent competitive local energy minima.
However, at each density one of them is usually lowest in energy.
Crossings between successive curves identify the coexistence density
values where the $n$-cluster and the $(n+1)$-cluster have the same
stability, and can therefore coexist.
By carrying out a finer numerical analysis and a Maxwell construction
one could also identify density ranges characterized by the phase
coexistence of homogeneous domains consisting of clusters of different
$n$.
We need not go into this detail: based on the obtained rough phase diagram,
we identify the densities corresponding to stable $n$-cluster states,
characterized by a cluster-cluster spacing $1.5 R_c$, and we adopt them for
the driven simulations of Sect.~\ref{driven:sec}.

\section{Driven Dynamics} \label{driven:sec}

\begin{figure}
  \centerline{
    \includegraphics[width=0.45 \textwidth,clip=]{snap_static_before_depinning.pdf}
  }
\caption{\label{cluster_pulled:fig}
   Small portions of the $T=0$ configurations in a $V_0 = 1.5\,E_0$
   corrugation, under the action of a $F=0.6\, F_{1s} < F_{\rm dep}$
   pulling force, i.e.\ these are static configurations prior to
   depinning, for:
  (a): $n=2$; (b): $n=3$; (c): $n=4$; (d): $n=5$ clusters.
}
\end{figure}

Starting from the equilibrated state obtained as discussed in
Sect.~\ref{gs:sec}, we investigate the dynamics under the competing
effects of the external periodic corrugation potential, which tends to
immobilize the clusters at its minima, and the lateral driving force $F$
which, if there was no corrugation, would tend to establish a sliding
state at speed $v=F/\eta$.
As expected of a commensurate, and even lattice-matched configuration, a
static pinning threshold is always present.
When $F$ is slowly raised from $0$ up to this pinning threshold, after an
initial transient allowing for the cluster rearrangement, the steady
sliding speed remains null.
However, nontrivial rearrangements of the clusters are observed well in
advance of depinning:
the driving force clears all orientational orderings characteristic of
the ground states.
The clusters are forced to specific orientations, which are displayed in
Fig.~\ref{cluster_pulled:fig}.
These ferro-orientational \cite{Ozaki99,Potocnik12} arrangements are made
energetically advantageous by the anisotropic energy landscape present at
the regions in between the minimum and the saddle point of the corrugation
potential of Fig.~\ref{ext_potential:fig} where $F$ pushes the clusters.
Precisely this force-induced ferro-orientational cluster arrangement
observed for all $n$ prior to depinning \cite{Reichhardt09a}
makes the orientational details of the free-model ground state
(e.g.\ Fig.~\ref{cluster2_orient_gs:fig}) irrelevant for the dynamics we
address in the present work.

\begin{figure}
\centerline{
\includegraphics[width=0.45 \textwidth,clip=]{Mobility_cl2.pdf}
}
\caption{\label{mob_cl2:fig}
  The $T=0$ mobility of the $n=2$-cluster lattice as a function of the
  driving force $F$ expressed as a fraction of the depinning force
  $F_{1s} = 8\pi V_0/(9a_{\rm pot})$ for an isolated particle, for
  corrugation amplitudes $V_0/E_0 = 1$, $1.5$, and $2$.
  The numerical error affecting the mobility due to averaging over
  noninteger numbers of periods is less than 0.3\% of the free mobility
  $\eta^{-1}$.
  The labels along the curve identify different sliding regimes, with
  groups of $m/n$ particles leaving one well of $U_{\rm ext}(\mathbf r)$,
  and moving to the next well at the right.
  Arrows identify transitions from a sliding regime to the one occurring at
  larger force.
  {\it Synchro}/{\it async} regimes refer to whether all well-to-well
  particle jumps occur simultaneously across the cell or they do not.
  In this specific $n=2$ case, even without specification, all motions
  are ``1D'', meaning that colloids advance along horizontal lines
  through the minima of $W(\mathbf r)$.
}
\end{figure}

\begin{figure}
\centerline{
\includegraphics[width=0.45 \textwidth,clip=]{Mobility_cl3.pdf}
}
\caption{\label{mob_cl3:fig}
  Same as Fig.~\ref{mob_cl2:fig} but for the $n=3$ cluster lattice.
}
\end{figure}

\begin{figure}
\centerline{
\includegraphics[width=0.45 \textwidth,clip=]{Mobility_cl4.pdf}
}
\caption{\label{mob_cl4:fig}
  Same as Fig.~\ref{mob_cl2:fig} but for the $n=4$ cluster lattice.
}
\end{figure}

\begin{figure}
  \centerline{
    \includegraphics[width=0.45 \textwidth,clip=]{Mobility_cl5.pdf}
  }
\caption{\label{mob_cl5:fig}
  Same as Fig.~\ref{mob_cl2:fig} but for the $n=5$ cluster lattice.
  All observed dynamical modes are synchronized here.
  See text for a description of the {\it very fast hop} and {\it $4/5$
    sync altern} regimes.
}
\end{figure}

The most straightforward dynamical indicator for a driven model of this
kind is the average mobility, i.e. the ratio of the average
$\mathbf{\hat x}$-directed velocity to the driving force, which we
report for different numbers $n$ of colloids per cluster.
If the corrugation potential was removed ($V_0=0$), then for any nonzero
force, the mobility would equal the free mobility $\eta^{-1}$, which
equals $0.0357\, \eta_0$ with the present choice of the damping
parameter.
The addition of the corrugation leads to smaller values of mobility.

To evaluate the mobility, for each $2\leq n\leq 5$ we execute three
sequences of simulations, with the following amplitudes of the periodic
external potential: $V_{0} = E_0$, $1.5\,E_0$ and $2\,E_0$.
For each sequence, we increase the applied lateral force in a range
proportional to $V_0$, namely $0.55\,F_{1s} \leq  F \leq 1.3\, F_{1s}$ in
small steps $\Delta F=0.015\,F_{1s}$.
At each step, the average velocity is evaluated by averaging over a total
time $1400\, t_0$ after an initial transient time $600\, t_0$.
Figures~\ref{mob_cl2:fig}--\ref{mob_cl5:fig} report the resulting
mobilities of the cluster lattices for the investigated numbers $n$ of
colloids per cluster as a function of the applied driving force $F$.
In the initial small-$F$ simulations the driving force is insufficient to
extract the colloid clusters from the potential wells of $U_{\rm ext}$:
thus a null mobility is obtained.
Beyond a certain static friction threshold $F_{\rm dep}$, the colloids
unpin and start to advance, resulting in a finite mobility, which grows
progressively approaching the free-particle value $\eta^{-1}$.
The small ``adiabatic''increase in $F$ allows us to evaluate $F_{\rm
  dep}$ with a relatively small uncertainty, and to identify the
features of the particle motion as mobility increases after depinning.

\begin{figure}
\centerline{
\includegraphics[width=0.45 \textwidth,clip=]{dep_force.pdf}
}
\caption{\label{dep_force:fig}
  The static friction, or depinning threshold, $F_{\rm dep}$, expressed as
  a fraction of the single-colloid depinning force $F_{1s}$ (also relevant
  for the depinning of the $n=1$ no-cluster phase), for $V_0/E_0 = 1$,
  $1.5$, and $2$, as a function of the number $n$ of colloids per cluster.
}
\end{figure}

Consider the depinning force $F_{\rm dep}$.
Although its values can be read directly from
Figs.~\ref{mob_cl2:fig}--\ref{mob_cl5:fig}, for ease of comparison
Fig.~\ref{dep_force:fig} reports them as a function of $n$.
A first observation is that the clustering of the colloids causes a strong
decrease of $F_{\rm dep}$ compared to the value $F_{\rm dep}=F_{1s}$ of a
single colloid, by about 20\% for $n=2$ and $3$ up to 33\% for $n=5$.
The static friction force peaks for $n = 3$.
Above $n = 3$ the depinning force tends to decrease as $n$ increases, as
expected due to the growing size of the clusters.

The $n = 2$ clusters mark an exception, as they depin at smaller force
than $n = 3$.
Also, in the comparison of the depinning force for three values of
$V_0/E_0$, we see that depinning for $n = 2$ occurs at a smaller $F_{\rm
  dep}/F_{1s}$ for $V_0 = E_0$ than $V_0 = 2 E_0$, contrary to all other
cluster sizes.
The reason for the special behavior of the $n = 2$ cluster lays in its
elongated conformation under traction, see Fig.~\ref{cluster_pulled:fig}.
As a result, in the pulling direction $\hat {\mathbf x}$ the $n = 2$ cluster
extends over a longer distance $\approx \delta$ than the size $\approx
\delta \sqrt{3}/2$ of $n = 3$.
For this reason, when the cluster center of mass approaches the inflection
point of the corrugation potential (the point beyond which the retaining
force begins to taper off), the rightmost particle has moved forward more
for $n = 2$ than for $n = 3$.

After depinning, the mobility, on top of an overall smooth increase as a
function of $F$, exhibits sharp changes in slope,
Figs.~\ref{mob_cl2:fig} - \ref{mob_cl5:fig}.
These changes mark transitions between different sliding patterns, which
can be examined by monitoring the advancement of the individual colloids.
The following section describes precisely these sliding regimes.

\subsection{Sliding regimes}\label{patterns:sec}

\begin{figure}
  \centerline{
    \includegraphics[width=0.95 \columnwidth,clip=]{5_5sync.pdf}
  }
\caption{\label{5_5sync:fig}
  The whole-cluster ($5/5$) advancement observed in the $n=5$, $V_0 =
  E_0$ simulation driven at $F/F_{1s} = 1$, see Fig.~\ref{mob_cl5:fig}.
  Panels (a) - (d) represent subsequent snapshots, separated by $\Delta
  t = 2\,t_0$.
  A particle is labeled with a vertical bar for ease of identification
  in the successive snapshots.
}
\end{figure}

Under driving, the internal structure of the clusters makes them advance
either as a single object, or less trivially with structural
decompositions and recompositions.
Depending on the corrugation amplitude $V_0$ and on the number of
colloid per cluster we observe (or not) several of these
decomposition-recomposition phenomena as the lateral force is increased.
For convenience, we indicate the observed dynamical patterns on the
mobility curves of Figs.~\ref{mob_cl2:fig} - \ref{mob_cl5:fig}, where
the switching between patterns can induce visible mobility anomalies.
We use $k/n$ labels indicating dynamical patterns which persist for an
entire force interval.
The $k/n$ labels in the mobility curves indicates that the advancement
from one well to the next occurs with $k$ particles out of a cluster of
$n$ jumping ahead, joining the next cluster, and leaving $n-k$ particles
behind.
Arrows point at transitions from one kind of dynamical behavior to the
next.

\begin{figure}
  \centerline{
    \includegraphics[width=0.95 \columnwidth,clip=]{2_5sync.pdf}
  }
\caption{\label{2_5sync:fig}
  Successive snapshots illustrating the $2/5$ dynamics observed in the
  $n=5$, $V_0 = 2\, E_0$ simulation driven at $F/F_{1s} = 0.7$, see
  Fig.~\ref{mob_cl5:fig}.
  Five particles are labeled for ease of identification in the
  successive snapshots.
  The time interval $\Delta t = 2\,t_0$, the same as in
  Fig.~\ref{5_5sync:fig}.
}
\end{figure}

As a general rule, a small corrugation amplitude $V_0$ tends to favor the
clusters advancing as a whole, namely in ``$n/n$'' modes.
In contrast, for larger amplitude $V_0$ the corrugation competes
favorably with the energy advantage of the clusters to remain entire,
and as a result for larger $V_0$ we do observe nontrivial
decomposition-recomposition patterns.
These decompositions can occur synchronously across the entire supercell
(and in this case they are labeled {\em synchro}) or asynchronously ({\em
  async} label).
Figures~\ref{5_5sync:fig} and \ref{2_5sync:fig} report successive
snapshots illustrating the two main kinds of dynamical patterns, as
marked on the mobility plots.
Specifically, Fig.~\ref{5_5sync:fig} illustrates the most basic $5/5$
advancement mode: whole cluster advancement from one corrugation well to
the next, as observed for weak corrugation $V_0 = E_0$.
In contrast, Fig.~\ref{2_5sync:fig} illustrates an advancement mode
involving synchronous decompositions and recompositions of the $n=5$
clusters, with $2$ particles abandoning each cluster and jumping to the
next corrugation well while leaving $3$ particles behind at each step.
Examples of similar patterned movements involving different numbers of
particles are reported in the Supplementary Material
\cite{supporting_current}, as digital videos.
Figure~\ref{scalini:fig} illustrates the effect of these patterned
movements on the $x$-coordinate of the center of mass of the colloids.
As the particle hoppings are synchronized across the entire sample, the
center-of-mass position $x_{\rm cm}$ exhibits relatively fast jumps
during the synchronized hoppings, alternating with plateaus during the
clusters rearrangements.
The heights of the observed steps are a fraction of $a_{\rm pot}$, which
matches precisely the specific dynamical ratio, e.g.\ the $2/5$ pattern
(solid line in Fig.~\ref{scalini:fig}) exhibits $\frac 25\, a_{\rm
  pot}$-high steps.

\begin{figure}
\centerline{
\includegraphics[width=0.45 \textwidth,clip=]{scalini.pdf}
}
\caption{\label{scalini:fig}
  The time-dependence of the $x_{\rm cm}$ of the system with $n = 5$
  colloids per cluster for $V_0 = 2\, E_0$ and for three different
  forces representative of the fully synchronized $2/5$, $3/5$ and $4/5$
  dynamical modes indicated in Fig.~\ref{mob_cl5:fig}: given $a_{\rm
    pot} = 1.5$, the steps heights are $0.62178 \simeq \frac 25 \,
  a_{\rm pot}$, $0.88585 \simeq \frac 35 \, a_{\rm pot}$ and $1.23788
  \simeq \frac 45 \, a_{\rm pot}$.
}
\end{figure}

When the driving force reaches relatively large values $F\gtrsim F_{\rm
  1s}$, its action on the system can force the collective dynamics into
an effectively one-dimensional (1D) dynamics, with all colloids
organized in rows parallel to the direction of the driving force
($\hat{x}$).
We noted this mode on the mobility plots, using the appropriate {\it 1D}
label, except for Fig.~\ref{mob_cl2:fig}, because for $n=2$ the motion
is quasi-1D for any $F$.
Similar quasi-1D dynamics was observed also in the different model of
Ref.~\cite{Reichhardt12, Reichhardt09a, Reichhardt09c}.
In contrast for $n = 5$, this quasi-1D advancement does not develop
completely, because of a significant transverse ($\hat{\mathbf
  y}$-directed) displacement persisting at all times due to the
hard-core repulsion not leaving enough room for 5 colloids in a row of
length $a_{\rm pot}$.
For $n = 5$, the large-$F$ dynamics can consist of a rapid and
uninterrupted transfer of particles from one cluster to the next, as
identified by the label {\it very fast hop} in Fig.~\ref{mob_cl5:fig}.
Alternatively we have observed a {\it 4/5 sync altern} motion, with groups
of 4 colloids jumping across clusters arranged in an alternating pattern (3
above + 2 below) -- (2 above + 3 below) along the rows, as illustrated in
video {\tt cl5\_4.44V2.avi} in the Supplementary Material
\cite{supporting_current}.

The synchronization of the well-to-well jumps can be promoted or
disfavored by weak cluster-cluster interactions.
As a consequence, we occasionally observe a de-synchronization of the
well-to-well hoppings for different clusters.
In the Supplementary Material \cite{supporting_current} we report examples
of asynchronous dynamics, where the advancements at different rows occur at
subsequent times.
Contrasted to the synchronous advancement steps of Fig.~\ref{scalini:fig},
in asynchronous dynamics as a consequence of averaging no (or very small)
steps are observed in the overall center-of-mass coordinate.
The asynchronous jumps can occur in a regular sequence at different rows or
with apparently chaotic hoppings as for example in the video
{\tt cl4\_2.71V2.avi} provided as Supplementary Material
\cite{supporting_current}.
Note that any regularity in asynchronous motions (and possibly the
overall synchronization itself when present) is favored by the specific
supercell periodic boundary conditions adopted in simulation: as such,
the synchronized or regular patterns may be considered artifacts of the
simulated model.

This zero-temperature dynamics, often dominated by unique regular or
chaotic attractors \cite{Alligood00,Ott02} does not represent a
realistic condition that one could put in relation to actual
experiments.
Indeed, in the more realistic finite-temperature simulations of
Sect.~\ref{finiteT:sec} below, irregular asynchronous depinning occurs
in most cases, due to the Brownian forces randomly anticipating or
retarding the advancement of individual particles or clusters.

By cycling the driving force $F$ up and down, we investigated whether these
dynamical systems retain some kind of dynamics memory, i.e. the sliding
state is affected not only by the competition of the applied confinement
potential and driving, but also by the initial sliding configuration.
Indeed in these $T=0$ simulation we do find some small difference of the
boundaries between specific dynamical modes between the decreasing-$F$ path
and the increasing-$F$ path.
To rule out non-physical underdamping effects, we verified that such
hysteretic effects persist even under a doubled damping rate $\eta = 56
\,\eta_0$, indicating that a competition among multiple dynamical
attractors is indeed present in this model.

In the present study the adopted $\hat {\mathbf x}$-directed orientation
of the driving force tends to drive the colloids straight across the
saddle points of the corrugation potential.
By exploring other more general pushing directions, depending on the
angle formed by the force with the $\hat{\mathbf x}$-direction, the
colloids would be driven toward the corrugation maxima too.
This kind of investigation has been carried out in the past for
different models or experimental setups, typically finding locking to
energetically and/or dynamically favorable advancement directions
\cite{Reichhardt04, Reichhardt11, Bohlein12PRL}.
We expect similar directional-locking phenomena in the present model.
However, possible novel regimes realized by the current cluster model
may include concurrent multiple decomposition paths for the clusters,
with several particles leaving simultaneously a potential well to reach
different neighboring ones.
We leave the investigation of these regimes to future research.

\subsection{Thermal Effects}\label{finiteT:sec}

\begin{figure}
\centerline{
\includegraphics[width=0.36 \textwidth,clip=]{s_of_q_combined.pdf}
}
\caption{\label{structure:fig}
  The scattering intensity $S({\mathbf q})$ obtained assuming that each
  colloid acts as structureless point-like scatterer.
  The simulations refer to $n=4$, at four temperatures expressed in units
  of energy (i.e.\ assuming $k_{\rm B}=1$): (a) $T=0$ -- the ordered state
  in the static condition depicted in Fig.~\ref{cluster_pulled:fig}c, (b)
  $T=0.01\,E_0$ -- rotational melting of the cluster lattice, (c)
  $T=0.15\,E_0$ -- crystal lattice phase just below structural melting, and
  (d) $T=0.17\,E_0$ -- a liquid phase just above melting.
  The $T>0$ simulations of panels (b-d) are carried out for the free model
  ($V_0 = 0$, $F=0$).
  $S({\mathbf q})$ is averaged over 30 snapshots widely spaced in time
  along an equilibrium simulation.
}
\end{figure}

We simulate finite temperature by restoring the random forces
${\mathbf\xi}_j$ in the equations of motion \eqref{langevin:eq}, and
therefore sampling the canonical ensemble.
The starting point of the $T>0$ simulations is the appropriate $T=0$
static configuration (a perfect cluster lattice) obtained in
Sect.~\ref{gs:sec}.
In the simulated scattering intensity reported in
Fig.~\ref{structure:fig}a, the modulation of the Bragg peaks reflects
the structure factor of the orientationally ordered lozenge-shaped
clusters at $T=0$, see Fig.~\ref{cluster_pulled:fig}c.
In the finite-temperature simulations of the free model ($V_0 = 0$, $F=0$),
it is straightforward to verify that even a low temperature, such as
$T=0.01\, E_0$ induces the melting of the orientational order of the
clusters, as seen in the diffuse rings at integer multiples of $2\pi/\delta
\simeq 20 \, R_c^{-1}$ reported in Fig.~\ref{structure:fig}b.
Further raising temperature leads rapidly to more disordering of the
cluster crystal, and eventually to melting, which occurs between
$T=0.15\, E_0$ and $T=0.17\, E_0$, see Fig.~\ref{structure:fig}c,d.

According to this analysis, we focus the investigation of the driven
model on moderately-low temperature $T\leq 0.05\, E_0$.
In the simulations under driving, as we did in the $T=0$ protocol,
successive runs at larger and larger driving force $F$ are started from the
final configuration of the previous step in $F$, which is increased
``adiabatically'' by $\Delta F = 0.01\,F_{1s}$ at each step.
For each step, the simulation duration is $500\,t_0$: the first 30\% is
dropped to prevent transient effects, and the remainder is used for the
determination of the mobility, providing a fair averaging over thermal
fluctuations.

Strictly speaking, at finite temperature there is no static friction,
because thermal fluctuations would lead to a diffusive displacement of the
layer even for $F=0$, and to a slow systematic drift in the force direction
for $F>0$, if one had the patience to wait long enough for these
thermally-activated rare events to occur.
However, at the relatively small considered temperatures $T\ll V_0$, the
rate of well-to-well thermal-assisted hopping is extremely small, and as
long as the applied force is small $F \ll F_{\rm 1s}$, we observe no such
hoppings at all for the duration of our simulations.
This allows us to define $F_{\rm dep}$ even for $T>0$, as the threshold
force beyond which a significant center-of-mass displacement is observed
in the applied force direction within the simulation duration.
For forces immediately before this threshold, the mobility is not
precisely null, but it is really tiny, with overall advancements of the
entire system by far less than one lattice spacing $a_{\rm pot}$ over
the entire simulation time.
We note also that all particle well-to-well displacements occur within
horizontal rows: although {\it a-priori} possible, we never observe any
particle abandoning their initial $\hat{\mathbf x}$-oriented row and
transferring to a well in the next row.
These cross-row hopping would of course become frequent for larger
$T/V_0$.
We do however observe particle jumps within rows, leading to a limited
but nonzero random concentration of $n-1$ and $n+1$ ``defective''
clusters in the pre-depinning quasi-immobile state.

\begin{figure}
\centerline{
\includegraphics[width=0.45 \textwidth,clip=]{cl3V2.pdf}
}
\caption{\label{cl3V2:fig}
  The mobility, relative to the value $\eta^{-1}$ of free colloids, as a
  function of the driving force, relative to the $T=0$ depinning threshold
  $F_{1s}$ of isolated colloids, at three relatively small but finite
  temperatures, expressed in units of energy $E_0$.
  The $T=0$ solid curve is the same as in Fig.~\ref{mob_cl3:fig}.
  Simulations are carried out for $n=3$-particle clusters, for corrugation
  $V_0 = 2 E_0$.
  Statistical fluctuations account for error bars (not drawn) on mobility
  mostly smaller than $5\times 10^{-3} \eta^{-1}$.
}
\end{figure}

The mobility resulting in the finite-temperature simulations is exemplified
in Fig.~\ref{cl3V2:fig}.
This figure exhibits the following characteristic features, common to
analogous simulations that we carried out for clusters of sizes $n=2$ to
$5$, and reported in Appendix:
(i) The depinning threshold $F_{\rm dep}$ decreases for increasing
temperature, because thermal fluctuations help the particles anticipate the
barrier hopping, thus activating sliding even when the driving force is
nominally insufficient to overcome the barrier.
(ii) For $F\gtrsim F_{\rm dep}$, immediately above the threshold, thermal
fluctuations tend to favor the sliding state, as shown by the higher
mobility proving a sort of thermolubric effect \cite{Gnecco00, Sang01,
  Dudko02, Riedo03, Reimann04, Krylov05, Jinesh08, Franchini08, Pierno14}.
However, depending on $n$ and the corrugation amplitude $V_0$, for forces
larger than the $T=0$ depinning threshold $F_{1s}$, increasing temperature
can lead to either a slightly decreasing mobility, as in
Fig.~\ref{cl3V2:fig}, or essentially $T$-independent mobility, quite
similar to the values obtained for $T=0$, see Appendix~\ref{thermal:app}.

The main effect of finite temperature on the cluster dynamics is to
suppress the $T=0$ synchronized advancements, especially those
associated to the fractional-advancement modes.
As a result, at $T>0$ the center-of-mass advancement is generally
smoother, with no, or at least substantially smeared, steps like those
occurring for the $T=0$ fractional sliding regimes reported in
Fig.~\ref{scalini:fig}.
Nonetheless, some amount of approximate synchronization is preserved in the
special case of whole-cluster hopping, e.g.\ at $T>0$ the {\it
  $3/3$~synchro} pattern of Fig.~\ref{mob_cl3:fig} even extends down to
smaller $F$ than for $T=0$.
Relatedly, the sharp transitions between different sliding modes observed
at $T=0$, at finite $T$ are replaced by smooth crossovers.

\section{Discussion and Conclusions}

In the present work we introduce and study a driven model for friction at
the micro scale where the protagonists are spontaneously formed clusters,
rather than simple structureless objects\cite{Thompson90b, Shinjo93,
  Smith96, He99, Muser00, Consoli00, Vanossi00, Muser01, Consoli01,
  Vanossi04a, Guerra10, BraunManini11, ManiniBraun11, vandenEnde12,
  Hasnain14, Varini15, Ma15, Koren16a} or structured but unbreakable
units\cite{molshapeLin, Sivebaek12}, as usually considered in previous
investigations.
We focus on a lattice-matched situation, with the external corrugation
periodicity matching exactly the equilibrium spacing between the clusters.
In the language of the Frenkel-Kontorova model \cite{Floria96, Braun98,
  Braunbook, Popov11, VanossiRMP13}, this is the typical situation which
maximizes the pinning effect of the corrugation, leading to a finite static
friction regardless of how small the corrugation amplitude may be.
We do find indeed a finite static friction in our simulations of the
current model, with clustering affecting the static-friction threshold and
the sliding mechanisms for driving force exceeding this threshold.

For small corrugation amplitude, the clustering tendency is dominating the
dynamics, so that clusters tend to advance as a individual objects, with
the partial exception of single-particle hopping induced by thermal
fluctuations near the depinning threshold.
In contrast, at larger corrugation amplitude, the stability of clusters is
challenged when the driving force moves them in regions where the curvature
of the external corrugation potential is negative: this leads to cluster
decomposition and partial hoppings.
In this regime, the internal structure of the clusters decomposes in
several ways as a function of the driving force, thus producing a quite
rich dynamics with complicate fractional advancement patterns for clusters
of various size.
We investigate the impact of these phenomena on the mobility of the system.
At $T=0$ we observe and characterize well-defined regimes of sliding of
entire clusters or of decomposition-recomposition processes, with
deformation of the clusters.
These processes can take place in a synchronous or in asynchronous way
across the system.
Thermal effects tend to destroy such distinct regimes: we observe a gradual
change of sliding under increasing driving force, already for temperature
as low as $\sim 5$\% of the temperature at which the transition from a
cluster crystal to a uniform fluid occurs.

The phenomenon of cluster decomposition/recomposition under driving is not
expected to be specific of the interparticle interaction adopted here.
We therefore predict that experiments carried out in conditions where
spontaneous clustering occurs should observe this kind of dynamical
behavior.

In the future we plan to extend the investigation to lattice-mismatched
configurations.
In incommensurate conditions, for increasing corrugation amplitude $V_0$
there is room for an Aubry-type transition from a ``superlubric'' dynamics
characterized by "solitonic" sliding modes \cite{Aubry78, Peyrard83,
  Floria96, Braun98, Santoro06, Vanossi07Hyst, Vanossi07PRL, Manini07PRE,
  Cesaratto07, Castelli08Lyon, Castelli09, Vanossi12PNAS, Vanossi13,
  MandelliPRL15, MandelliPRB15, Paronuzzi16, Mandelli17, Bonfanti17}.
A similar condition could be realized by taking for the substrate a
symmetry other than hexagonal, as e.g.\ in Ref.~\cite{Persson93b,
  ManiniBraun11, McDermott13a, McDermott13b}.
In the small-$V_0$ superlubric ground state, a finite fraction of clusters
is supposed to sit near one of the maxima of the incommensurate corrugation
potential, rather than all in a potential well as in the fully commensurate
setup of the present paper.
The novelty for a cluster-supporting condition is that, as $V_0$ is raised,
two competing effects may arise: either an entire cluster near the maximum
slides down toward a nearby potential well, thus leading to a relatively
standard Aubry-type transition \cite{MandelliPRL15,Mandelli17}, or
alternatively it decomposes under the action of divergent corrugation
forces, leading to a structurally inequivalent ground state.
Which of these mechanisms turns in earlier is open for investigation.
These competing possibilities are likely to lead to a rich phase diagram.

\begin{acknowledgments}
  We acknowledge useful discussion with C. Bechinger, R. Guerra, and
  A. Vanossi.  The COST Action MP1303 is also gratefully acknowledged.
\end{acknowledgments}


text:
\section{Introduction}
Ultracold atoms provide a testing ground for many-body physics~\cite{Bloch}. Utilizing Feshbach resonance management, the interatomic interaction can be tuned to any desired value~\cite{interactions}. This tunability allows the production of weakly or strongly correlated Bose Einstein condensates (BECs). A variety of different geometries and topologies is realizable by suitably manipulating the magneto-optical trapping potentials~\cite{arbpot}. For instance, interacting bosons in optical lattices, that is, a spatially periodic potential, have been shown to exhibit a quantum phase transition from a superfluid state (SF) to a Mott insulator (MI)~\cite{Jaksch,Greiner,Haller}.

One-dimensional optical lattices loaded with strongly interacting bosons have been shown to be experimentally more challenging than their three-dimensional counterparts: quantum fluctuations are not negligible~\cite{haldane}, trigger correlations \cite{haller1,kinoshita} and are the focus of the present paper.

The zero-temperature SF to MI transition is commonly described by the Bose-Hubbard model~\cite{Fisher}. For weak interparticle interactions, a Bose gas in an optical lattice of moderate depth is in the superfluid phase: the many-body state features long-range coherence. This SF phase persists, as long as the interatomic interaction is small compared to the tunneling coupling -- i.e. the parameter of the Bose-Hubbard Hamiltonian that determines the lattice depth and hence tunneling between neighboring sites. The SF can be characterized as a state in which all particles are delocalized between all sites of the lattice. When the repulsive interaction is large compared to the tunneling coupling, each site is filled up with an identical number of bosons and the coherence between different sites is completely lost: the many-body state enters the Mott-insulating phase.

Even though the Bose-Hubbard model can aptly describe the above SF to MI transition in lattices, its regime of validity is restricted to the case where site-localized Wannier states are an appropriate basis set~\cite{Grigory}. It has been demonstrated that a general quantum many-body description, valid at all interaction strengths, is necessary for the regime beyond the Bose-Hubbard model~\cite{Sakmann2009}. One such theory is the  the multiconfigurational time-dependent Hartree (MCTDH) method \cite{mctdh}. In Ref.~\cite{Brouzos} MCTDH was applied to small lattices and it was found that correlations depend both on the commensurability and the strength of the interparticle interactions. In the present paper, we solve the full many-body Schr\"odinger equation at a high level of accuracy by using the multiconfigurational time-dependent Hartree method for bosons (MCTDHB)~\cite{Ofir,OfirJCP} implemented in the MCTDH-X software~\cite{axel1,axel2,ultracold}.
The motivation for our present work is to investigate the pathway from superfluid through the Mott insulating to the fermionized phase and explain it in terms of production of many-body information entropy.

Entropies of quantum systems were long ago introduced and used as a measure of the degree of order of a given quantum state~\cite{Massen1,Massen2,Massen3}. Entropy measures for many-body states were shown to saturate in time to the values given by estimates of Gaussian orthogonal ensembles (GOE) of random matrices for the time evolution triggered by an interaction quench~\cite{SIE_1}. This situation, when different entropy measures take on the GOE values, is referred to as \emph{statistical relaxation} and was shown to also reflect in the correlation functions of the many-body state~\cite{SIE_1}.

Below, we analyze this connection of many-body entropies and correlation functions for the quantum phases of many-body states in an optical lattice, commensurately filled with bosons that interact with a repulsive contact interaction. We characterize the many-body state by calculating the produced entropy and the spatial correlations as a function of the interparticle interaction strength. We demonstrate that the production of different types of entropies -- defined below -- and the correlation functions can be used to identify the phase of the system.

The \emph{Shannon information entropy (SIE)} can be computed from the one-body density, i.e., the diagonal of the first-order reduced density matrix, in coordinate or momentum space~\cite{Guevara1,Guevara2,Sagar,Sen}. In this case, SIE is a measure of the spatial delocalization of the one-body probability distribution of the system. Since SIE of the one-body density is computed from a \emph{single-particle quantity}, it is insensitive to the many-body physics of the considered system and its correlations.

In this work, in order to remedy this shortcoming of the SIE of one-body quantities, we generalize the concept and introduce SIEs on the many-body level. In particular, we define and investigate: $a)$ the \emph{two-body Shannon information entropy} calculated from the second-order reduced density matrix, $b)$ the \emph{occupation SIE (O-SIE)}, i.e. the entropy of the eigenvalues of the reduced first-order density matrix~\cite{Brezinova2012,Tsatsos2015,SIE_1} and $c)$ the \emph{Shannon information entropy of the coefficients (C-SIE)} of the expansion of the state on some many-body basis (see Sec.\ref{Hamiltonian})~\cite{SIE_1}. These entropies -- as opposed to the single-particle SIE -- directly relate to the many-body physics of the system. The O-SIE and C-SIE are many-body measures that allow us to get a quantitative criterion for the applicability of theories like, for instance, the Gross-Pitaevskii (GP) mean-field theory; the C-SIE and O-SIE take on extremal values for mean-field GP states while any other value heralds beyond-mean-field physics.
In the following, we study the many-body entropy production of interacting ultracold atoms in a one-dimensional optical lattice as a function of the strength of the interparticle interactions and the number of bosons per site. We establish a relation of the phase diagram to the C-SIE and O-SIE entropy measures. We complement the emergent picture with an analysis of Glauber's normalized first- and second-order correlation functions~\cite{RJG}. Refs.~\cite{SIE_1,Brezinova2012} have already demonstrated that a fundamental connection between the correlation functions and many-body entropy measures exists, at least for the temporal evolution following a quench. In our present study, we demonstrate that this relation between correlations and many-body entropy also holds for the ground states of bosons in optical lattices.

For small barriers and weak interactions, the system is a superfluid, aptly described within a mean-field theory and characterized by a smooth increase (decrease) of the Shannon information entropy of the one-body density in coordinate (momentum) space~\cite{Massen1,Massen2,Massen3} for increasing interaction strength.

For larger barriers or stronger interactions and one atom per site, i.e. commensurate filling, the system becomes a Mott insulator with no phase coherence between distinct sites. For this MI phase, the reduced one-body density matrix acquires multiple eigenvalues: fragmentation emerges \cite{Oliver,spekkens,mueller,lode17,PV,PNAS}. We characterize the many-body state through the superfluid-to-Mott-insulator phase transition using our many-body entropy measures.  For sufficiently strong interparticle interactions we observe a complete saturation of entropies to values predicted by the Gaussian Orthogonal Ensemble of random matrix theory~\cite{Kota2001}.

If the commensurate filling of the lattice is larger than one atom per site, a second phase transition emerges for very strong interparticle interactions: the fragmented Mott insulator transmogrifies into a so-called fermionized state~\cite{Gir,Ofir2}. The onsite repulsion forces the density of the system to develop an intrasite structure. The C-SIE and O-SIE measures of many-body information entropies also saturate for large interactions in the case of the larger filling factor with, however, no limiting value from a Gaussian orthogonal ensemble existing for this case.

We quantify the coherence properties of all three above phases -- condensed superfluid, fragmented Mott insulator and fermionized phase -- by computing Glauber's normalized first- and second-order correlation functions. We find that the classification of phases using correlation functions agrees with the classification of phases using many-body entropy measures. The fundamental relation between many-body entropies and coherence properties in ultracold bosonic atoms, found in Refs.~\cite{SIE_1,Brezinova2012} in the time evolution following a quench, is extended to the present case of stationary states in lattices.

The paper is structured as follows. In Sec.~\ref{Hamiltonian} we expose the Hamiltonian and introduce the many-body theory and our quantities of interest: Glauber correlation functions and many-body entropy measures. Section~\ref{results} deals with the results of our analysis of states for various interaction parameters and filling factors. We conclude in Sec.~\ref{Conclusion}.

\section{Theory} \label{Hamiltonian}

\subsection{Hamiltonian}
Consider a system of $N$ bosons interacting by a contact interparticle interaction potential in one spatial dimension. In dimensionless units, such a system is governed by the following Hamiltonian:
\begin{equation}
H = \sum_{i=1}^N \left( \frac{1}{2} \frac{\partial^2}{\partial x_i^2} + V(x_i) \right) + \lambda \sum_{i<j}^N \delta (x_i - x_j). \label{eq:HAM}
\end{equation}
Here, $V(x)= V_0 \sin^2(kx)$ is the lattice potential, where $k=\frac{\pi}{d}$; $V_0$ is the depth and $d$ is the periodicity of the lattice. The strength of the two-body interactions, $\lambda$, can be experimentally tuned almost at will in quasi-one-dimensional systems by manipulating the strength of the transversal confinement~\cite{1d_scattering}. In the remainder of the paper, we set a depth of $V_0 = 12.0$, $x_{\text{min}}=-4.7124$, $x_{\text{max}}=4.7124$ and $d=3$ wells are considered. The depth $V_0$ is chosen such that it allows superfluidity, for appropriately chosen boson number $N$ and interaction strength $\lambda$. We find the stationary solutions of the many-body Schr\"odinger equation with periodic boundary conditions and obtain the observables defined in Sec.~\ref{QOI} as a function of the strength of the interparticle interactions $\lambda$ and the number of atoms per lattice site.

\subsection{MCTDHB}
In the multiconfigurational time-dependent Hartree for bosons (MCTDHB) approach, the wave function of the interacting $N$-boson problem is expanded over a set of permanents. Permanents are symmetrized bosonic states of $N$ particles in $M$ single-particle states. Each permanent can be constructed by acting products of $N$ boson creation operators $b^\dagger_k$ ($k=1,...,M$) onto the vacuum $\vert vac \rangle$:
\begin{eqnarray}
& \vert \Psi(t) \rangle & ~=~ \sum_{\vec{n}} C_{\vec{n}}(t) \vert  \vec{n}; t \rangle; \label{cansatz} \\
& \vert \vec{n} \rangle & ~= \vert n_1,n_2, ...  n_M \rangle = \prod_{i=1}^M \left[ \frac{\left(b_i^\dagger (t)\right)^{n_{i}}}{\sqrt{n_i !}}\right] \vert vac \rangle.
\label{ansatz}
\end{eqnarray}
Here, each operator $b^\dagger_k(t)$ creates a boson occupying the time-dependent single-particle state (orbital) $\phi_{k}(x,t)$. The number of possible configurations of $N$ bosons in $M$ orbitals is equal to $\binom{N+M-1}{N}$ and defines the number of complex-valued coefficients $C_{\vec{n}}(t)$ in Eq.~\eqref{cansatz}. A formal variational treatment with the above ansatz leads to the MCTDHB equations of motion~\cite{Ofir,Streltsov2007,OfirJCP}. The solution of the latter yields the time evolution of the coefficients $C_{\vec{n}}(t)$ and orbitals $\phi_k(x,t)$ that built-up our solution: an approximation to the solution of the time-dependent Schr\"odinger equation at a desired (arbitrary) degree of accuracy. In this work, we find the eigenstates of the MCTDHB equations by propagating the MCTDHB equations in imaginary time using the MCTDH-X package~\cite{axel1,axel2,ultracold}.
In the limit of infinite orbitals, $M\rightarrow\infty$, the set of permanents $\vert n_1,n_2, \dots  n_M; t \rangle$ in Eq.~\eqref{ansatz} spans the complete $N$-particle Fock space and MCTDHB becomes exact~\cite{exact,exact2,exact3}.

\subsection{Quantities of interest}\label{QOI}
We now introduce the quantities that we use to characterize the solutions of the MCTDHB equations.
\paragraph{\bf Entropies:}
The Shannon information entropies of the one-body density in coordinate space $\rho^{(1)}(x)=\langle \Psi \vert \hat{\Psi}(x)^\dagger \hat{\Psi}(x) \vert \Psi \rangle$ is given by
\begin{equation}
S_x(t) = - \int dx \rho^{(1)}(x,t) \ln [\rho^{(1)}(x,t)].
\label{SIE1x}
\end{equation}
Analogously, for the momentum density, $\rho^{(1)}(k)=\langle \Psi \vert \hat \Psi (k)^\dagger \hat \Psi(k) \vert \Psi \rangle$ we have:
\begin{equation}
S_k(t) = - \int dk \rho^{(1)}(k,t) \ln [\rho^{(1)}(k,t)].\label{SIE1k}
\end{equation}
The SIEs in Eqs.~\eqref{SIE1x},\eqref{SIE1k} are a measure of the delocalization of the corresponding distributions $\rho^{(1)}(x)$ and $\rho^{(1)}(k)$. Refs.~\cite{Massen1,Massen2,Massen3} establish a universal relation between entropy and the number of interacting particles for diverse systems like atoms, nuclei, and atomic clusters.
Since the distributions in Eqs.~\eqref{SIE1x},\eqref{SIE1k} are relate to the one-body density, they are insensitive to correlations that may be present in the state $\vert \Psi \rangle$~\cite{PNAS,PV,Tsatsos2015,Brezinova2012}.
We can, however, formulate an SIE using the two-body density distributions $\rho^{(2)}(x_1,x_2)=\langle \Psi \vert \hat\Psi(x_1)^\dagger \hat\Psi(x_2)^\dagger \hat \Psi(x_1) \hat{\Psi(x_2)} \vert \Psi \rangle$ and $\rho^{(2)}(k_1,k_2)=\langle \Psi \vert \hat\Psi(k_1)^\dagger \hat\Psi(k_2)^\dagger \hat\Psi(k_1)\hat\Psi(k_2) \vert \Psi \rangle$:
\begin{equation}
S_{\rho-x}(t) = - \int dx_1 dx_2 \rho^{(2)}(x_1,x_2;t) \ln [\rho^{(2)}(x_1,x_2;t)],
\label{SIE2x}
\end{equation}
and
\begin{equation}
S_{\rho-k}(t) = - \int dk_1 dk_2 \rho^{(2)}(k_1,k_2;t) \ln [\rho^{(2)}(k_1,k_2;t)]
\label{SIE2k}
\end{equation}
Here $\rho^{(2)}(x_1,x_2)$ [$\rho^{(2)}(k_1,k_2;t)$] is the diagonal part of the two-body reduced density matrix in position (momentum) space. Since the SIEs in Eqs.~\eqref{SIE2x},\eqref{SIE2k} are computed from two-body quantities, they are measures sensitive to correlations in the many-body state $\vert \Psi \rangle$. By comparing the SIE based on the one-body density [Eqs.~\eqref{SIE1x} and \eqref{SIE1k}] with the SIE based on the two-body density [Eqs.~\eqref{SIE2x},\eqref{SIE2k}] the presence of correlations in the state $\vert \Psi \rangle$ can be inferred.

Since the state we consider is expanded in the MCTDHB theory as $\vert \Psi \rangle = \sum_{\vec{n}} C_{\vec{n}} (t) \vert \vec{n};t \rangle$, we can define an alternative SIE using the coefficients $C_{\vec{n}}$ that characterize the distribution of the state $\vert \Psi \rangle$ in the underlying Fock space [cf. Eq.~\eqref{ansatz}]:
\begin{equation}
S_c(t) = -\sum_{\vec{n}} |C_{\vec{n}}(t)|^2  \ln \left[|C_{\vec{n}}(t)|^2\right]. \label{SIEC}
\end{equation}
We term this entropy \emph{coefficient Shannon information entropy} (C-SIE), or simply coefficient entropy.
A mean-field state is a single-configuration state, i.e. only a single coefficient contributes in Eq.~\eqref{ansatz}. For such a state $S_c(t)=0$ holds at all $t$. Coefficient entropy $S_c$ thus cannot be produced in a mean-field theory.
When the state $\vert \Psi \rangle$ spreads across several configurations $\vert \vec{n} ; t \rangle$, several expansion coefficients contribute [cf. Eq.~\eqref{ansatz}] and the coefficient entropy $S_{c}$ gradually increases.
In the limiting case, when the complete $N$-body Fock space is populated by the state $\vert \Psi \rangle$, all coefficients are equally large and $S_c$ saturates to its maximal value.

Last, we define an SIE measure related to the emergence of fragmentation \cite{Brezinova2012,Tsatsos2015,SIE_1}, i.e., the emergence of multiple significant eigenvalues $n_i;\;i=1,\dots,M$ of the reduced one-body density matrix of the state $\vert \Psi \rangle$~\cite{Oliver}. These eigenvalues are also referred to as natural occupations and thus we term the following measure of entropy the occupation Shannon information entropy (O-SIE):
\begin{equation}
S_n = - \sum_i^M n_i \ln n_i .
\label{SIEO}
\end{equation}
For a state described by a single-orbital mean-field theory, the reduced density matrix is characterized by only a single eigenvalue and hence $S_n=0$ holds. For single-configuration states with multiple contributing orbitals as well as multiconfigurational states, there maybe several occupation numbers and hence $S_n\neq 0$.
For the condensed gas only a single occupation dominates and the occupation entropy is zero. For an increase in interaction strength, the O-SIE gradually increases and saturates only for a maximally fragmented state.

\paragraph{\bf Correlation functions:}
The normalized $p$-th order correlation function is defined by
\begin{equation}
 g^{(p)}(x_1^{\prime},...,x_p^{\prime},x_1, ...,x_p)=
\frac{\rho^{(p)}(x_1,...,x_p \vert x_1^{\prime},...,x_p^{\prime})}
{\sqrt{\prod_{i=1}^p\rho^{(1)}(x_i \vert x_i )\rho^{(1)}(x_i^{\prime} \vert x_i^{\prime} )}}
\end{equation}
and is the key quantity to define spatial $p$-th order coherence. Here, $\rho^{(p)}(x_1,...,x_p\vert x_1^{\prime},...,x_p^{\prime}; t)$ is the $p$-th order reduced density matrix of the state $\vert \Psi \rangle$~\cite{RJG}.
In the case of $\vert g^{(p)}(x_1...,x_p,x_1 ...,x_p;t) \vert > 1$ ($<1$), the detection probabilities of $p$ particles at positions $x_1,...,x_p$ are referred to as (anti-)correlated.
Recent progress in experiments in quantum gases has been remarkable and the measurement of higher-order correlation functions is now possible \cite{Hodgman2011,Dall2013,Langen2015,Schweigler2017,Tsatsos2017}. In particular, in the work of Ref.~\cite{Schweigler2017} it has been explicitly shown how a many-body system is characterized via its higher-order correlations.

\section{Results}\label{results}

In this section, we display our numerical results for the SIE measures and spatial correlations in the eigenstates of the Hamiltonian of Eq.~\eqref{eq:HAM} as a function of the strength of repulsive interactions, $\lambda$ . We choose three different sets of system parameters representative for the three distinct phases that the system may be in. Results for the superfluid phase, the Mott insulating phase, and the fermionized phase are discussed in the following Sections~\ref{SF}, \ref{MI}, and \ref{FERM}, respectively.

We define the filling factor as the ratio of the number of atoms $N$ and the number of lattice sites $W$: $\nu= \frac{N}{W}$. For the entire manuscript we focus on commensurate filling factors $\nu=1, 2$ and $7$. We keep the lattice depth $V_0 =12.0$ and lattice site number $d=3$ fixed. We thus change the number of particles to change the filling factor and the interaction strength is gradually tuned to cover all the emergent phases for each filling factor. The transitions between the emergent phases are discriminated using information entropies and correlation functions.

We choose to present the results in the order of increasing complexity of the found many-body state: $\nu=7,1,2$. Indeed, we find that the $\nu=7$ case is captured within a mean-field theory for the range of the interparticle interaction strengths which we investigate. The $\nu=7$ and $\nu=1$ cases can qualitatively be described within the Bose-Hubbard model as there is no structure that forms within sites. Notably, the fermionization that emerges for the $\nu=2$ case at stronger interparticle interactions is a many-body phase that neither of the aforementioned approaches can describe.

\subsection{Commensurate filling factor $\nu=7$: Superfluid state}\label{SF}
Here, we characterize the superfluid phase using natural occupations, normalized correlation functions, and the Shannon information entropy measures that we have introduced. In Fig.~\ref{fig.fig1}(a), we plot the occupations as a function of the interaction strength $\lambda$. For relatively weak interactions, only the first natural orbital is occupied and the population of the second and third orbital remain below $10\%$ for all values of $\lambda<0.01$. In this parameter regime the system is a condensate, according to Penrose and Onsager~\cite{Oliver}, and one can well approximate the many-body wave function with a single-orbital mean-field state $\vert N,0,... \rangle$. To achieve convergence in the occupations, i.e. the eigenvalues of the reduced one-body density matrix the interaction strengths that we consider, three orbitals [$M=3$, cf. Eq.~\eqref{ansatz}] are enough; adding more orbitals does not quantitatively change the many-body state.

We plot the one-body SIE for the density and momentum distributions as a function of the interaction strength $\lambda$ in Fig.~\ref{fig.fig1}(b). As $\lambda$ gradually increases the (momentum) density within each well is broadened (narrowed) and consequently the SIE computed from the density (momentum) distribution, $S_{x}$ ($S_{k}$) [see Eqs.~\eqref{SIE1x} and \eqref{SIE1k}] increases (decreases). Note that the curves for $S_x$ and $S_k$ cross each other for $\lambda =0.007$. In Fig.~\ref{fig.fig1}(b), we also present the C-SIE $S_{c}$ and O-SIE $S_{n}$. As anticipated, the values of $S_{c}$ and $S_{n}$ are very close to zero, because the condensate is phase coherent and can be well described by a macroscopically occupied single-particle state and the Gross-Pitaevskii mean-field equation yields a good description.

In Fig.\ref{fig.fig1}(c), we plot the two-body SIE [Eqs.~\eqref{SIE2x},\eqref{SIE2k}]. The two-body SIEs behave similarly to the one-body entropy measures; the crossing of the momentum-based and the real-space-based SIE happens at the same point $\lambda =0.007$. However the rate of increase in $S_{\rho -x}$ is greater than that in $S_x$ and, likewise, the rate of decrease of $S_{\rho -k}$ is greater than that in $S_k$. A similar trend was observed in the calculation of one-electron and two-electron entropies in Ref.~\cite{Guevara1}. The fact that the two-body density matrix can -- to a good approximation -- be written as a product of one-body densities for the coherent superfluid gas is a possible explanation for this behavior.

\begin{figure}
\includegraphics[angle=270,width=\linewidth]{Fig1.png}
\caption{{\bf Entropies and occupations in the superfluid phase.} All quantities are shown for $\nu=7$ and $N=21$ bosons as a function of the interaction strength $\lambda$. (a) Eigenvalues of the reduced density matrix -- the ground state is condensed and only a single eigenvalue is macroscopic for the superfluid phase. (b) One-body SIE $S_x$ and $S_k$, C-SIE $S_c$ and O-SIE $S_n$. (c) Two-body SIE $S_{\rho -x}$ and $S_{\rho -k}$. All quantities are dimensionless. See text for discussion.}
\label{fig.fig1}
\end{figure}

In Fig.~\ref{fig.fig2}(a), we plot the absolute value of the normalized first-order correlation function, $\vert g^{(1)}(x^{\prime},x)\vert^2$ for $\lambda= 0.01$. We see that in the $(x,x')$ region where the density is localized $\vert g^{(1)}\vert^2 \approx 1$ holds. We infer that coherence within and between sites is maintained. Fig.~\ref{fig.fig2}(b) shows the two-body correlation function $g^{(2)}(x^{\prime},x,x^{\prime},x)\equiv g^{(2)}(x,x^{\prime})$ for the same parameters as in Fig.~\ref{fig.fig2}(a). We find $g^{(2)}(x,x')$ to be close to unity for all $x',x$. Remarkably, second-order coherence between different wells is perfectly maintained: $g^{(2)}(x,x')\approx 1$ for the off-diagonal $x'\neq{}x$ with $\vert x'-x \vert$ larger than the size of a single well.
The diagonal part of the two-body correlation function is slightly depleted, i.e. $g^{(2)}(x,x')\lesssim 1$ for the diagonal $x' \sim x$ if $\vert x'-x \vert$ is smaller than the size of a single well. This signifies that anti-bunching starts to develop for particles within the same site due to the repulsive interactions: second order coherence is locally decreased, even though the interaction strength $\lambda$ is relatively small.

\begin{figure}
\includegraphics[angle=270,width=\linewidth]{Fig2.png}
\caption{{\bf (a) first- and (b) second-order normalized correlation function for the superfluid phase, i.e. $\nu=7, N=21$.} The interaction strength is $\lambda=0.01$ while the rest of the parameters are as in Fig.~\ref{fig.fig1}. (a) Complete first-order coherence is observed within each well;  $\vert g^{(1)} \vert \approx 1$ for all $x\approx x'$. Significant inter-well coherence can be inferred from $\vert g^{(1)} \vert \in [0.5,1]$ for all $x\neq x'$ on the off-diagonal. (b) Second-order coherence is retained among distinct wells $g^{(2)}(x,x')\approx 1$ for all $x\approx x'$. The fact that $g^{(2)}$ is maximal among distinct wells but drops within each well indicates that the probability for detecting two particles in different wells is higher than that of finding them in the same well. This behavior is a consequence of the interparticle interactions and it is termed antibunching. All quantities are dimensionless. See text for discussion.}
\label{fig.fig2}
\end{figure}

\subsection{Commensurate filling factor $\nu=1$: transition from condensation to fragmentation}\label{MI}
We now turn our attention to the Mott insulating phase and investigate its many-body physics from the viewpoint of correlation functions, coherence and Shannon information entropies. To obtain an archetypical Mott-insulator for strong interactions, we consider $N=3$ bosons in three wells, i.e. a single atom per site. To quantify the departure of the many-body state $\vert \Psi \rangle$ from a mean-field state in the superfluid-to-Mott-insulator transition, we plot the expansion coefficients $\vert C_{n} \vert^{2}$ [Eq.~\eqref{ansatz}] as a function of the index $n$ of the basis states for two choices for the interaction strength: $\lambda=0.01$ which puts the system onto a superfluid phase and $\lambda=10.0$ that renders it a Mott-insulator [see Fig.~\ref{fig.fig3}(a) and (b)].
In the superfluid phase ($\lambda=0.01$) the number of significant coefficients is only a small portion of the total number of configurations $N_{conf}=\binom{N+M-1}{N}$. We refer to such a state as localized. We emphasize that our use of the term ``localized'' refers to many-boson Fock space and not real-space; localized states are thus close to a mean-field description where only a single coefficient would contribute. However, as $\lambda$ increases a larger amount of coefficients become significant; we refer to such a state as delocalized [cf. $\lambda=10.0$ in Fig.~\ref{fig.fig3}(b)].

\begin{figure}
\includegraphics[angle=270,width=\linewidth]{Fig3.png}\\
\caption{{\bf Localization and delocalization in Fock space in the superfluid to Mott-insulator transition at $\nu=1,N=3$.} Panels (a) and (b) show the distribution of the magnitude of the coefficients $\vert C_{n}\vert^{2}$ as a function of the index $n$ for the superfluid at interaction strength $\lambda=0.01$ in (a) and the Mott insulator at interaction strength $\lambda=10.0$ in (b). The index $n$ is computed from the vector $\vec{n}$ using the mapping described in Ref.~\cite{MAP}. In (a) the $\vec N,0,...\rangle$ coefficient and its neighbors dominate while all others are small; the state is \emph{localized} in the many-body Fock space. At (b) the larger interactions force the many-body state to spread at a larger part of the available space and many coefficients are significant; the state is \emph{delocalized} in the many-body Fock space. Panel (c): Populations of the first three natural orbitals. As $\lambda$ increases, the occupation of the first orbital gradually decreases while another two orbitals begin to contribute (green and blue curves). Eventually the state becomes three-fold fragmented ($n_1=n_2=n_3=1/3$) at $\lambda\gtrsim6.0$. All quantities are dimensionless.}
\label{fig.fig3}
\end{figure}

In Fig. \ref{fig.fig3}(c), we plot the corresponding occupation of the first, second, and third natural orbital as a function of the interaction strength. With an increase of $\lambda$, the occupation of the first natural orbital gradually decreases while the occupations of the second and third natural orbitals gradually increases. For large interactions, the first three natural occupations saturate at $33.33\%$. In our computations with $M=6$ orbitals, all but the first three orbitals had occupations smaller than $10^{-7}$.
Three-fold fragmentation in a triple well indicates loss of partial coherence of distinct sites and is a signature for the transition to the insulating phase \cite{Streltsov2011}. Thus, across the superfluid-to-Mott-insulator transition the increase of the interaction strength drives the state from a condensed to a fragmented one. Generally, in the Mott-insulating phase there are as many significant eigenvalues of the reduced density matrix as there are lattice sites.

\begin{figure}
\includegraphics[angle=270,width=\linewidth]{Fig4.png}\\
\caption{{\bf One-body density and entropies in the superfluid-to-Mott-insulator transition at unit filling $\nu=1$ and $N=3$ as a function of the interaction strength.} (a) The one-body density as a function of interaction strength $\lambda$ shows almost no variation in the transition from a condensed superfluid to a fragmented Mott insulator. $(b)$: Shannon information entropy of the one-body (momentum) density $S_x$ ($S_k$), the occupation entropy $S_n$ and the coefficient entropy $S_c$ all as a function of $\lambda$.
(c) SIE of the reduced two-body (momentum) density matrix $S_{\rho-x}$ ($S_{\rho -k}$). For the SIE measures in (b) and (c), the curved part for weak interactions signifies the transition from a condensed superfluid to a fragmented Mott insulator. Within the Mott-insulating phase, for larger interaction strengths $\lambda$, all SIE measures saturate. All quantities are dimensionless, see text for further discussion.}
\label{fig.fig4}
\end{figure}

From Fig.~\ref{fig.fig4}(a) we see that the above transition cannot be inferred from the density alone. Indeed, the density before and after the transition (in the $\lambda$ space) are identical and many-body measures, such as the entropies, are required to mark the transition. Examining the C-SIE $S_{c}$ and O-SIE $S_{n}$ we observe that for small  $\lambda$ both start at zero, as only a single coefficient contributes to the the state, which is a condensed superfluid. The behavior of the SIE measures computed from the one- and two-body density are in mutual agreement with the C-SIE and O-SIE measures, compare Figs.~\ref{fig.fig3} and Fig.~\ref{fig.fig4}(b),(c). As $\lambda$ rises, the C-SIE and O-SIE measures increase and the many-body state converges to a pure Mott-insulator while fragmentation, as well as all SIE measures saturate.

We shall now explore this saturation of entropies using an analogous Gaussian Orthogonal Ensemble (GOE). The GOE of random matrices can be formulated for isolated quantum systems of interacting particles which exhibit time reversal symmetry and rotational invariance~\cite{rigo,sant3,rigol}. The entropies evaluated for the GOE provide estimates for the entropies of many-body systems for the limiting case of maximal disorder. The GOE of random matrices with infinite interaction has entropy $S_{c}^{GOE}=\ln (0.48 D_c)$, where $D_c \times D_c$ is the dimension of the considered random matrices.

In the present $\nu=1$ case, three bosons are distributed in (essentially) three orbitals, $N_{orb}=3$ [Fig.~\ref{fig.fig3}(c)]. Hence, the size of the GOE is $D_c= \binom{N + N_{orb} -1}{N} = 10$ thus obtaining $S_{c}^{GOE} = 1.568$. This is in excellent agreement with the value at which our numerical result for the C-SIE $S_c$ saturates: $1.552$. The small discrepancy is because the interaction strength is still finite. The coefficient entropy thus saturates at the value predicted from the GOE for sufficiently large $\lambda$.
For the O-SIE, $S_n$, the dimensionality of the GOE $D_n=3$ is equal to the number of occupied orbitals in our MCTDHB treatment and $S_{n}^{GOE} = -\sum_{i=1}^{D_n}\frac{1}{D_n} \ln\left(\frac{1}{D_n}\right) = \ln (D_n) = 1.098$. Our numerical result for the saturated value of the O-SIE is identical, $S_{n} = 1.098$ \cite{M9comparison}.

\begin{figure}
\includegraphics[angle=270,width=\linewidth]{Fig5.png}\\
\caption{{\bf First-order and second-order coherence of a fragmented Mott-insulator at $\nu=1, N=3$.} (a) The normalized first-order correlation function $\vert g^{(1)}(x,x')\vert^{2}$ for $\lambda=10$ exhibits three separated regions with $\vert g^{(1)} \vert^2\approx1$ along the diagonal and $\vert g^{(1)} \vert^2\approx0$ on the off-diagonal. Thus, the first-order coherence is maintained within wells and lost between them. (b) The normalized second-order correlation function $g^{(2)}(x,x')$ also for $\lambda=10$. The ``correlation hole'' is clearly seen and $g^{(2)}\approx 0$ holds for the diagonal $x\approx x'$; two particles are thus almost never simultaneously found in the same well. Second-order coherence $g^{(2)} \approx 1$ on the off-diagonal signifies that, in a measurement, two uncorrelated particles will almost surely be detected in distinct wells. The Mott-insulating phase is characterized by independent lattice sites with strongly localized particles. All quantities are dimensionless.}
\label{fig.fig5}
\end{figure}

We analyze the first-order and second-order coherence in the fragmented Mott insulator for $\lambda = 10.0$ in Fig.~\ref{fig.fig5}. The diagonal of the first-order correlation function shows three completely separated coherent regions where $\vert g^{(1)} \vert^2\approx1$; the coherence of bosons within the same well is maintained while it is lost between distinct wells, since $ \vert g^{(1)} \vert^2\approx 0$ at the off-diagonal, as expected for the fully localized particles in a Mott insulator.
Looking at $g^{(2)}$ we see that coherence is maintained ($g^{(2)}\approx 1$) at the off-diagonal but not at the diagonal.
The vanishing diagonal part of the normalized two-body correlation function -- referred to as the correlation hole~\cite{Brouzos} -- translates to the fact that the probability of finding a double occupation of a single well is practically zero. Second order coherence between wells is maintained, because the outcome of a two-particle detection is almost always a certain and uncorrelated detection of two particles in distinct wells.

\begin{figure}
\includegraphics[angle=270,width=\linewidth]{Fig6.png}
\vspace{-9mm}
\caption{{\bf One-body density and entropy measures for the transition from the Mott-insulator to fermionization.} All panels are for $\nu=2,N=6$ and are plotted against the interaction strength $\lambda$. (a) Plot of the one-body density; the transition to fermionization is visible in the formation of a two-hump intrawell structure -- for large $\lambda$, the number of peaks in the density is equal to the number of particles $N$. The dashed contour line $\rho(x)=0.22$ is drawn to visualize the two maxima of each well. (b) SIE computed from the one-body (momentum) density $S_x$ ($S_k$), the coefficients entropy $S_c$, and the occupation entropy $S_n$. (c) SIE computed from the two-body (momentum) density, $S_{\rho-x}$ ($S_{\rho-k}$).
Compared to the system with a single boson per well, the SIE measures saturate only for comparatively large interaction strengths [compare panels (b) and (c) to Fig.~\ref{fig.fig4} panels (b) and (c), respectively]. All quantities are dimensionless, see text for further discussion.}
\label{fig.fig6}
\end{figure}

\subsection{Commensurate filling factor $\nu=2$: from a condensed superfluid via a fragmented Mott insulator to fermionization}\label{FERM}
We now consider a system with $\nu=2$ bosons per site and characterize its phases using the density, normalized correlation functions and Shannon information entropy measures. In lattices with commensurate filling factor larger than one, i.e. with more than one boson per site, two phase transitions can be observed as a function of the strength of the interparticle interactions $\lambda$: The superfluid-to-Mott-insulator transition, that is also present if the filling factor is equal to one, emerges for intermediate interaction strengths. For large interaction strengths the Mott-insulator-to-fermionization transition -- absent for $\nu=1$ -- now occurs. Since the physics of the superfluid-to-Mott-insulator transition is unchanged, we focus here on the Mott-insulator-to-fermionization transition. The pathway from the condensed superfluid via the fragmented Mott insulator to the fermionization regime depends on the shape of the trapping potential~\cite{Ofir2}. The optical lattice with filling factor two is the generic system to for this transition to be studied.

Similarly to Sec.~\ref{MI}, we see that the system transitions to MI at $\lambda\approx10$. The transition cannot be seen in the one-body density, since the latter remains unchanged across the SF-to-MI transition as $\lambda$ increases [see Figs.~\ref{fig.fig6}(a) and \ref{fig.fig4}(a)].
When the interaction strength $\lambda$ is sufficiently large ($\lambda\gtrsim200$), the Mott-insulator-to-fermionization transition takes place and the number of maxima in the one-body density becomes equal to the number of bosons of the system, see Fig.~\ref{fig.fig6}(a). In this case the bosonic distribution resembles a fermionic one where Pauli exclusion forbids the spatial overlap of the particles.
The transition to fermionization is thus manifest in the formation of intrawell structure in the density. Such intrawell structures cannot be appropriately described by the single-band Hubbard model. When a complete set of Wannier states is considered in each site, i.e. a multi-band Hubbard model, an intractably large number of bands quickly becomes necessary~\cite{optimal_lattices}. Thus, the transition to fermionization is a hallmark for the inapplicability of Hubbard models.

\begin{figure}
\includegraphics[angle=270,width=\linewidth]{Fig7.png}\\
\caption{{\bf First- and second-order coherence for $N=6$ fermionized bosons and $\nu=2$.} Both panels show normalized correlation functions for an interaction strength $\lambda=200$. (a) The normalized first-order correlation function $\vert g^{(1)}(x,x')\vert^{2}$. The coherence between wells is completely lost, i.e. $\vert g^{(1)}\vert^{2}\approx 0$ for the off-diagonal. The fermionization of the system is manifest in the loss of first-order coherence also within each of the wells [cf. Fig.~\ref{fig.fig5}(a)]. (b) The normalized second-order correlation function $g^{(2)}(x,x')$ has a fully developed correlation hole, $\vert g^{(2)}\approx0$ on the diagonal. No bunching or anti-bunching is seen between particles in distinct wells, i.e.$g^{(2)}\approx1$ on the off-diagonal. Fermionization is manifested in the fact that $g^{(2)}<1$ as long as $x,x'$ are within each of the wells [cf. Fig.~\ref{fig.fig5} (b)]. All quantities shown are dimensionless, see text for further discussion.}
\label{fig.fig7}
\end{figure}

The SIE measures plotted in Fig.~\ref{fig.fig6}(b) and \ref{fig.fig6}(c) as a function of $\lambda$ also saturate in the transition to fermionization. The saturation of SIEs in the Mott-insulator-to-fermionization transition with two bosons in each well, however, takes place at a much larger interaction strength ($\lambda\sim200$) as compared to the saturation of SIEs for the superfluid-to-Mott-insulator transition ($\lambda\sim10$) in the case of $\nu=1$. Interestingly, no GOE analogue for the fermionized state has been established so far and, consequently, we cannot provide a comparison the values that the entropies shall saturate to.

We plot the normalized one-body and two-body correlation functions $g^{(1)}$ and $g^{(2)}$ to determine how the Mott-insulator-to-fermionization transition reflects in the coherence of the many-body state, see Fig.~\ref{fig.fig7}.
The first and second order coherence between different wells is similar for the fragmented and the fermionized Mott insulator, i.e. the off-diagonal parts of the normalized correlation functions $g^{(1)}$ and $g^{(2)}$ are practically identical (cf. Figs.~\ref{fig.fig5} and \ref{fig.fig7}).
However, in contrast to the correlation functions of the fragmented Mott insulator with a single boson per site, the correlation functions for the fermionized Mott insulator feature a distinct intrawell structure (cf. Figs.~\ref{fig.fig5} and ~\ref{fig.fig7}).

\section{Conclusions}\label{Conclusion}
We have shown that various information entropies can be defined to reflect and quantify many-body properties of the ground state of interacting bosons in one-dimensional lattices. These entropies serve the dual purpose to $i)$ identify phase transitions not seen in the single-particle properties of the system and $ii)$ allow a connection to known results from random matrix theory.

We demonstrated that our SIE measures are in agreement with predictions for the Gaussian orthogonal ensemble (GOE) of random matrix theory in the case of strong interactions. For sufficiently strong interactions and a single boson per site the system is fully in the Mott insulating regime and all our SIE measures saturate to their respective GOE values. For sufficiently strong interactions and two bosons per site the system is fermionized and the SIE measures are saturated. For this case, however no GOE estimation is available.

All the emergent phases -- the superfluid, the Mott-insulator, and the fermionized state -- are identified with a distinct value of their SIE measures. Our analysis of the SIE measures is complemented by an analysis of the density and the spatial first- and second-order coherence. We showed that all emergent phases have characteristic density, correlation, and coherence patterns.
We hence demonstrated that the superfluid, the Mott-insulating and the fermionized phases of bosons in one-dimensional lattices are not only characterized uniquely by their densities and Glauber normalized correlation functions~\cite{RJG,Ofir1}, but also by their SIE values. This corroborates the fundamental connection between correlation functions and entropy measures conjectured already in Refs.~\cite{SIE_1,Brezinova2012}.

In the transition to fermionization, intrawell structure in the densities and correlations of the many-body state is formed; this marks the breakdown of the single-band Hubbard model. Since the SIE measures and Glauber normalized correlations allow to identify the fermionized state, we also demonstrated two independent ways to assess the applicability of the Hubbard description.

Studying higher-order correlations and their `holes' in the $(x,x')$ would consist a natural extension of our works. Moreover, phases beyond fermionization and the relation of our findings for the stationary states to quench dynamics are open questions.

\begin{acknowledgments}
RR acknowledges UFC fellowship and B. Chakrabarti acknowledges FAPESP (grant No. 2016/19622-0). AG and MCT acknowledge FAPESP and AG thanks CNPq for financial support. B. Chatterjee acknowledges financial support from the Department of Science and Technology, Government of India under the DST Inspire Faculty fellowship. AUJL acknowledges financial support by the Austrian Science Foundation (FWF) under grant No. F65 (SFB ``Complexity in PDEs''), and the Wiener Wissenschafts- und TechnologieFonds (WWTF) project No MA16-066 (``SEQUEX'').
\end{acknowledgments}

redpajama_set_name: RedPajamaArXiv