An active colloidal system showing parallels to a time crystal

The spontaneous breaking of symmetries is a widespread phenomenon in physics. When time translational symmetry is spontaneously broken, an exotic nonequilibrium state of matter in which the same structures repeat themselves in time can arise. This state, known as ‘time crystal’, attracted a lot of interest recently. Another relatively new research area deals with active matter. Materials consisting of colloidal particles that consume energy from their environment and propel themselves forward can exhibit intriguing properties like superfluidity that were previously known only from quantum-mechanical systems. Here, we bring together these—at first glance completely different—research fields by showing that self-propelled colloidal particles show parallels to classical continuous time crystals. We present a state diagram showing where this state of matter arises. Furthermore, we investigate its properties and analyze the interactions between the particles leading to the dynamics.


Introduction
The proposal of time crystals [1,2] as a theoretical concept in 2012 opened a new area of research [3][4][5][6][7][8][9][10][11][12][13].A time crystal is characterized by spontaneously breaking time-translational symmetry and forms an analog to spatial crystals that break translational symmetry in space.Similar to spatial crystals, time crystals form a state of matter.Soon, it was proven that time crystals cannot exist in equilibrium for local interactions [9,10], although they could be realized in non-local equilibrium systems [14].Thus, when assuming local interactions, time crystals are a non-equilibrium phenomenon and never reach thermal equilibrium.Around the same time, the concept of discrete time crystals, which describe a discrete breaking of time translational symmetry in periodically driven systems, was proposed.While continuous time crystals break continuous time-translational symmetry, discrete time crystals occur in periodically driven systems but oscillate at an integer fraction of the driving frequency [15].Thus, they show a discrete time-translational symmetry breaking.Their existence in quantum mechanical systems was experimentally proven for the first time in 2017 [3,4].Different suggestions for (discrete) classical time crystals were discussed [16][17][18][19][20]. Recently, a realization of a classical continuous time crystal in a photonic metamaterial has been proposed [21].
Active colloidal particles are nano-or microparticles characterized by the capability to consume energy from their environment and to use it for self-propulsion [22].These particles include a large variety of motile microorganisms, such as archaea, bacteria, and many eukaryota.On the other hand, many types of artificial microswimmers, that utilize, e.g., acoustophoresis, chemophoresis, thermophoresis, or shape deformations for self-propulsion, have been developed in recent years [22][23][24][25][26].Such active particles have been found to exhibit a fascinating nonequilibrium behavior where various exotic properties can arise.Some of these properties, such as superfluidity, were previously known only from quantum systems and other properties, such as a negative viscosity, were previously unknown from both classical passive (i.e., not active) systems and quantum systems [27,28].Active colloidal particles have also been shown to give rise to a complex state diagram with various unusual nonequilibrium states including active crystals [29,30], which are colloidal crystals formed by active particles.Lately, there have been initial studies investigating the formation and stability of active colloidal crystals [31][32][33][34][35][36][37].However, the properties of this state of active matter remain mostly unexplored up to now.
Recently, the interaction of colloids with a choreographic time-crystalline lattice formed by traps was studied [38].Also, a first connection between an active fluid and a time crystal was proposed [39].However, this proposal was based on an abstract field-theoretical model, that can be related to active fluids but is not applicable to active crystals, and the study did not demonstrate the emergence of a time-crystalline state for a specific and realistic system of active particles.Other researchers found dynamic states in active systems [40][41][42] but did not analyze them with respect to properties of time crystals.
In the following, we establish a concrete connection between active colloidal crystals and time crystals.We show for a specific and realistic system that active colloidal crystals show parallels to classical continuous (and potentially even living) time crystals.Furthermore, we study the state diagram of these dynamic states and investigate their properties and the interactions between the particles.

Methods
To obtain new insights into the properties of active colloidal crystals, we use computer simulations combining low-Reynolds-number hydrodynamics with Brownian particle dynamics.We investigate a small active crystallite consisting of seven spherical particles with radius R that are placed in form of a hexagon with one central particle and surrounded by a liquid.The spatial layout of our system is shown in figure 1.The hexagon was chosen for the spatial configuration because it forms a part of a hexagonal grid which is common in nature.The two-dimensionality of the system avoids adding the complexity of a spatial configuration in three dimensions.In a three-particle system, we found the dynamic effects to be less clear and, therefore, more difficult to analyze.
The particles are placed at fixed positions with a distance d to each other but can freely change their orientations due to hydrodynamic interactions.We chose a minimum distance between the particles' centers of d = 2.2R to ensure enough space between the particles' surfaces for a numerically stable computation of the hydrodynamic interactions without requiring very high resolutions.Each swimmer creates a flow field caused by its swimming.This flow field influences the other swimmers and causes them to rotate.The particles themselves are modelled as squirmers, which today is a standard model for active swimmers [43].

Hydrodynamics
The active particles generate a flow field in the surrounding liquid that leads to hydrodynamic interactions of the particles, which affect their motion.Since the size of microswimmers typically lies in the range of micrometers, the Reynolds number r h = RU Re 0 with the liquid's mass density ρ, particle radius R, maximal flow speed U 0 , and dynamic shear viscosity η is much smaller than one.This allows to approximate the Navier-Stokes equation describing the flow field by the Stokes equation ( ) where v denotes the liquid's velocity field and p its pressure field.Any external forces applied to the system, such as gravity, can be neglected in this work.Additionally, the fluid is treated as incompressible: From the solution of the Stokes equation (1), one can calculate the stress tensor [22] T where I is the identity matrix, ⊗ is the dyadic product, and T denotes transposition.
Since the positions of the particles are fixed, they can only rotate but not translate.Therefore, we can neglect the forces and need to consider only the torques exerted on the particles by hydrodynamic interactions, when studying their motion.Based on the stress tensor (3), the hydrodynamic torque T h on a particle can be calculated as [22] where r is the position of a surface element with area ds and normal vector n on the particle surface S.

Model for microswimmers
We model the active particles as spherical squirmers as proposed by Lighthill [43] and Blake [44].This model is widely used for various microswimmers ranging from biological ones like Paramecia to artificial ones like Janus particles.Depending on the flow field generated by the squirmers, they are distinguished into pullers, which propel themselves mainly with the front similar to the algae Chlamydomonas reinhardtii, neutral squirmers, where the swimming motion happens everywhere on the body like for Paramecia, or as pushers, which propel themselves with the back part of their body like Escherichia coli [22].The squirming parameter β = B 1 /B 2 , where B 1 and B 2 are the first and second velocity mode, respectively, allows the distinction between pushers (β < 0), pullers (β > 0), and neutral squirmers (β = 0) [45].In this work, we choose the values β = − 3, 0, 3 for the squirming parameter to be consistent with [45].To facilitate the treatment of the spherical particles, spherical coordinates r, θ, f are used.The particle surface is located at r = R, where R is the radius of the particles.On the surface of the particles, the velocity of the surrounding liquid is prescribed as Here, ŝ is the tangential vector on the surface and equivalent to the unit vector ˆq e .The analytical solution of the Stokes equation around a particle [45] leads to a swimming velocity . 6 0 1

Equation of motion
Starting from Newton's equation of motion, one obtains for the considered system the Langevin equation that describes the motion of an active particle.Here, J is the angular momentum of the particle, γ R its rotational friction coefficient, ω its angular velocity, and ξ R (t) is a Gaussian white noise random torque with variance 2D R that takes Brownian rotation into account.γ R is related to the rotational diffusion coefficient of the particle by where k B is the Boltzmann constant and T the temperature.
In the regime of small Reynolds numbers, the rotational friction is large compared to  J .Thus, inertia can be neglected leading to This equation represents an energy balance between the energy lost by friction, the energy obtained by the hydrodynamic torque that acts on the particle, and the torque contribution from the Brownian motion.Using additionally the relation with the particle director û (a unit vector that describes the particle's orientation), one obtains If the swimming velocity is high and, therefore, the dimensionless Péclet number Pe = B 1 /(RD R ) is large, the Brownian motion is also negligible because T h increases with the velocity mode B 1 .Thus, we get In the simulations where an external torque T ext is taken into account, we simply add to T h in the equations above, where T ext is the scalar torque and êz the unit vector corresponding to the z axis.

Periodic motion in active crystallite
For an active crystallite formed by pushers with squirming parameter β = − 3 and a distance of d = 2.2R between the particles, we observed a periodic rotational motion after an initial relaxation (see figure 1).It is mainly characterized by the rotation of the central particle which rotates clockwise.The other six particles only slightly change their orientation, but this motion follows the same periodicity as the one of the central particle.
The rotation of the central particle can be explained by the orientation of the outer particles.They self-organize in a circular pattern where each particle has a similar angle with respect to the central particle when one takes the symmetry into account.Thus, the part of the outer particle that causes the propulsion of the pusher (colored in yellow in figure 1(a)) is oriented towards the central particle causing its acceleration.In this state of periodic rotational motion, all particles lie in the x-y-plane and the particle directors, which denote the orientations of the particles, do not leave this plane for the chosen parameter setting.Note that figure 1(b) only shows the xcomponent of the particle direction vector.Therefore, the different amplitudes are caused by the different orientations around which the periodic motion occurs.
The plane forms a fixed point of the system, as can be seen by symmetry considerations.Thus, in absence of noise the particles should not rotate out of the plane once all particles reached it.We do not impose planar swimming directions as initial conditions, but the system self-organizes into the planar motion.This suggests that this fixed point is stable.We investigate its stability in more detail in section 3.4.The orientation of the outer particles is strongly dependent on the propulsion mechanism and, therefore, on the parameter β leading to the self-organization of the dynamic state for pushers.We investigate the dependency on the different parameters in section 3.2.
For a more general analysis including all degrees of freedom of the particles, we combine their seven orientation vectors (i.e., directors) ˆ( ) u t i with i = 1,K,7 and time t to a 21-dimensional vector T that describes the state of the whole system at a certain time t.We project the 21dimensional trajectory to a lower-dimensional space using a principal component analysis (PCA) [53] (see supplementary material for more details).For this periodic state, we found that the first two principal components already contain χ = 99.3% of the variation in the data set.By visualizing these first two components (see figure 2(a)), we obtained a spherical motion.Thus, the system can be described by a simple two-dimensional spherical motion in the 21-dimensional phase space.
This periodic motion spontaneously breaks the continuous time-translational symmetry in a way that is similar to the breaking of translational symmetry occurring in an ordinary crystal in space.Analogous to the periodicity in space that is a characteristic of a normal crystal, our system is periodic in time (see figure 1), which is typical for a time crystal.This is especially interesting as the active particles are no chiral particles but the selfpropulsion of a free particle has an axis of rotational symmetry that points through the center of the spherical particle.
To call a system a time crystal and distinguish it from other well-known nonequilibrium phenomena that show periodicity, like oscillating chemical reactions or convection cells, it has to fulfill four criteria in addition to the time-translational symmetry breaking [16], namely independence on fine-tuning of parameters, stability against fluctuations, periodicity arises from interactions, and periodicity holds in larger systems.In the following, we discuss how our system fulfills these criteria.

State diagram
First, time crystallinity forms a state of matter and, thus, should not require fine-tuning the system's initial conditions or parameters.We confirmed that this is the case here by performing the simulations for different, randomly set initial conditions and various parameter combinations (see Methods for details).
We start with the initial conditions.Repeating the simulations that showed the periodic state described in section 3.1, we observed that the length of the initial relaxation varied, but in all simulations the swimmers showed the same qualitative steady-state behavior as depicted in figure 1.To verify this finding more quantitatively, we introduce a distance measure where m 1 is the closest local maximum of the r-dependent function ∥u(t) − u(r)∥ to t with m 1 < t, m 2 is the closest local maximum with m 2 > t, and t max is the last simulated point in time.Thus, δ(t) describes how similar the system's state at time t is to its state at another time.For periodic motion, the trajectory should come back to a given point for every time and, therefore, δ(t) should become 0, potentially after a transition phase.Note that the vector u(t) describes the state of the system relevant for the motion and, thus, fully characterizes the motion in the absence of perturbations.Therefore, δ(t) = 0 identifies a periodic motion.For non-periodic cases, however, δ(t) is larger than zero, as the trajectory never reaches any original point again.With this definition, we can define the end of the initial relaxation for periodic motions as the point where δ(t) becomes zero.In figure 2 of the Supplementary Information, we plot the distance measure δ(t) for 10 different initial conditions.We see that after some initial relaxation, the distance decreases to zero for the different initial conditions.This confirms that the occurrence of a periodic state does not rely on particular initial conditions.
To analyze how the state of the system depends on the values of its parameters, we now investigate the system's state diagram (see figure 3), which shows the behavior of the active crystallite for particle distances d = 2.2R, 2.4R,K,3R and squirming-parameter values β = − 3, 0, 3 that correspond to different types of swimmers (pushers for β < 0, neutral squirmers for β = 0, pullers for β > 0).
Interestingly, self-organization and periodic behavior are observed only for pushers, whereas for pullers and neutral squirmers only chaotic-looking aperiodic motion is found (see figure 3(a)).This shows that the type of interaction between the particles is crucial for the dynamic motion (see section 3.4 for details).This finding agrees with other studies of squirmers which also found that pushers show a qualitatively different behavior [27].
In the case of pushers, different states occurring for different particle distances can be distinguished.Periodic motion of the particles is observed for the smallest considered particle distance d = 2.2R (see figure 3(b)).For a larger particle distance of d = 2.4R, the motion is almost periodic but with clearly visible fluctuations (see figure 3(c)).When the particles have an even larger distance d = 2.6R or d = 2.8R, the motion looks more irregular and suggests a coexistence.At some times, it seems as if the system is converging towards a periodic state, but never reaches it completely; at other times, the motion looks chaotic (see figure 3(d) and, for a larger time frame, figure 3 of the Supplementary Information).When comparing the latter state to classical systems, we can interpret it as a dynamical coexistence between ordered and unordered (temporal) regions.This reveals parallels to a coexistence of a liquid phase and a crystalline phase.We therefore consider this state as a coexistence state.For the largest considered particle distance d = 3R, aperiodic motion similar to the behavior of pullers and neutral squirmers is observed (see figure 3(e)).The transition from a periodic to an aperiodic motion for increasing particle distance can be explained by the decrease in the hydrodynamic interaction between the particles for increasing particle distance.
To better understand how the system's full trajectory u(t) varies between the different states, we observe the PCA results for the individual states as shown in figure 2. The circular structure in phase space of the timecrystalline state (see figure 2(a)) disappears when d or β is increased.We can see that the dimensionality of the curve in phase space increases and the ordered motion in the first principal component vanishes.While we can represent χ = 97.7% of the motion in four principal components for d = 2.4R (see figure 2(b)), we already need ten components for χ = 94.1% of the variation in case of d = 2.6R (see figure 2(c)).In the latter case, even the first component shows no order.For β = 0, the first ten components only cover χ = 87.6% of the variation without any visible order (see figure 2(d)).
For studying the transitions between the observed states, we define an order parameter d as the time average of the minimum distance of the 21-dimensional system trajectory u(t) to itself: Here, δ(t) is the distance measure introduced in equation (12).Thus, d = 0 applies to periodic states and d > 0 to non-periodic states.Compared to frequency-based order parameters, this one allows to capture gradually changing behavior, if the system's trajectory is not periodic but moves in a confined region of the state space.For the calculation of the transitions, we run simulations starting with the same initial conditions but with varied values for d and β.We refine the sampling of the parameter range in regions of a transition.All simulations are run over a time span of 1000R/B 1 .For the calculation of the order parameter, we use the results of the last 500R/ B 1 to exclude the initial relaxation.To verify that we still take a sufficient number of time steps into account and avoid finite size effects, we investigate the dependence of d on t max as shown in figure 4. We can observe little variation for the four states shown figure 3 if t max is larger than 350R/B 1 .There is a small decrease for d = 2.6R and d = 3R, but this is to be expected as these states show aperiodic motion.Therefore, increasing the observed time will lead to a denser filling of the 21-dimensional space which induces a reduction of d.As the smaller distances show no further variation, we can conclude that our chosen time interval of 500R/B 1 is sufficiently large.
The simulation results for a variation of the particle distance d are shown in figure 5(a).Until a distance of 2.32R, we see a periodic state which corresponds to an order parameter of d = 0 and a very narrow distribution in δ.It is followed by a continuous transition with bare analogies to a second-order phase transition in equilibrium systems.For an increasing distance, up to d = 2.51R, we can see fluctuations in the particles' trajectories which result in a growing but still relatively small order parameter and a broader, unimodal distribution.At d = 2.51R, a jump in the order parameter occurs which seems analogous to a discontinuous, first-order phase transition in equilibrium systems.We can also see that the distribution becomes multimodal.For larger distances, there is no visible order in the particles' trajectories which results in fluctuating, larger values of d.However, when observing the individual trajectories of the system, we can distinguish between the coexistence state and the fully aperiodic state as shown in figures 3(d) and 3(e), respectively.The boundary between both states can be located at about d = 2.95R.However, this transition between both states is not visible in the curve for the order parameter, which is consistent with our previous finding that the former of the two states resembles rather a coexistence state than an individual state.
The evolution of the order parameter d, when β is varied, is depicted in figure 5(b).Here, we can identify a periodic behavior for strongly negative values of β and a discontinuity at β = − 2. This discontinuity has bare analogies to a first-order phase transition.For β = − 0.5, we reach a special case.There, the system converges towards a static state which also manifests itself in an order parameter of d = 0.However, by observing the orientations of the individual particles, we can confirm that no periodic motion but a static state occurs there.In this state, all particles but the central one are oriented in the same direction.The central particle aligns with the zaxis of the system.
The previous results show that fine-tuning either the initial conditions or the systems parameters is not required.

Influence of fluctuations
Second, the system should be stable against fluctuations.All simulations are influenced by numerical errors that are caused by the chosen granularity of the underlying mesh and by rounding errors.As is it usual, we ensured that the limited resolution has no significant influence on our simulation results by performing comparative simulations with a finer discretization.Even though the numerical errors are small, they can influence the system's behavior at large time scales.Although the system's symmetry should prevent particles from moving out of the particle plane, they do this in some cases for parameter combinations that correspond to non-periodic  states.For the periodic state, however, we did not observe such a behavior, suggesting that this state is stable against numerical fluctuations.
Heating of the system is prevented by coupling to a bath introducing friction and Brownian noise.We therefore studied also the influence of Brownian motion on the system's behavior.For this purpose, we included Brownian motion in the particles' equations of motion.Brownian motion is commonly introduced in active systems consisting of colloidal particles, since their motion is influenced by collisions with solvent molecules.We study the influence of Brownian motion as a function of the dimensionless Péclet number Pe = B 1 /(RD R ), where B 1 denotes the surface velocity mode of the active particles that also determines the swimming velocity and D R is their rotational diffusion coefficient.The results for different values of Pe are shown in figure 6.We can see that the overall motion of the system is robust against Brownian motion if it is small compared to the propulsion of the particles.The Euclidean distance d T between the system's trajectory with Brownian motion and the corresponding trajectory without Brownian motion is presented in figure 6(a).Here, we calculate the distance d T between two 21-dimensional vectors of the trajectories at time t as where u B,i (t) denotes the ith component of the 21-dimensional vector for a trajectory including Brownian motion and u i (t) denotes the ith component of the vector at time t for the unperturbed trajectory.To take shifts in time into account, we apply dynamic time warping [54] before calculating the distances.Dynamic time warping is a common technique to measure distances between time series that are not synchronized.For this purpose, the algorithm computes the optimal match between two time series.This is necessary for our comparison between the trajectories including and excluding Brownian motion as the Brownian motion induces shifts in time even though the overall structure of the system's behavior is maintained (see figures 6(b)-(d)).For our analysis, we use the implementation of an efficient approximation algorithm for dynamic time warping called FastDTW [54] that is provided by the fastdtw package [55].Overall, the trajectories get closer to the undisturbed trajectory if Pe increases.This follows the expectation as an increasing value of Pe means a lower influence of the Brownian motion.Even the motion of particles with Pe = 100 shows similarities to the motion of the undisturbed system (see figure 6(b)), albeit the influence of the Brownian motion is clearly visible.Furthermore, we can see that the trajectories for Pe 1000 only show small deviations.Comparing the motion of the individual particles confirms that the characteristic behavior is preserved with Brownian motion (cf figures 6(c), (d)).Especially the central particle still rotates even though the frequency varies, whereas the influence of the fluctuations on other particles is larger.

Particle-particle interactions
Third, the time periodicity should only arise from the interaction between the particles, i.e., the degrees of freedom should be coupled locally.This also applies to our system, where the hydrodynamic interactions are crucial.Without interactions between the particles, they would not move at all.The trajectories of single particles for Pe = 100 (not included in (a)) are strongly disturbed but the central particle (particle 1) shows a similar motion as in the undisturbed case.c, d For increasing Pe, the influence of the Brownian motion becomes smaller and the motion of the particles approximates the undisturbed motion.
We study the interactions between two particles in the system and their dependence on the orientation angle α of a particle, which is measured counterclockwise within the particle plane (i.e., the x-y-plane) relative to the position of the other considered particle, as well as on the parameters d and β to better understand how the dynamic states emerge and why transitions between the states occur.This also helps to understand how the propulsion type influences the system's behavior.We observe the influence of a particle of the system on itself and on the other particles.As the periodic behavior is a planar motion and the motion in the plane forms a fixed point in the dynamics of the system, we consider both the torque T z exerted on a particle that leads to a rotation within the particle plane and the torque T rot that rotates the particle out of or back into the plane.
In figure 7(a), we can see that the swimmers' hydrodynamic interactions mainly influence their nearest neighbors.The torque applied on a second nearest neighbor is, on average, only 1.5% of the one applied on a first nearest neighbor.Moreover, the self-interaction strength of the outer particles is comparable the interactions with their nearest neighbors.In the case of the central particle, in contrast, the self-interaction is very small.
Studying the stability of the fixed point that is formed by the motion in the plane allows understanding why the particles can rotate out of the plane for particular parameter configurations.The strength of the torque T rot that rotates a particle back into the plane (or further out of the plane) for small deviations from the plane is shown in figures 7(b) and (c).Here, positive values mean a rotation towards the plane whereas negative torques lead to a rotation further out of the plane.The dependence of this torque on the distance d is shown in figure 7(b).Here, we observe large values for small distances.For increasing d, the torque decreases and it converges towards 0. The dependence on β, which is depicted in figure 7(c), shows relatively large stabilizing torques for the time-crystalline state that decrease approximately linearly and become negative for positive values of β.The values of T rot are small compared to the torques T z that cause a rotation within the plane and not much larger than the numerical uncertainties that lead to the visible fluctuations in the curves.These observations show that the periodic motion within the plane is stable for small distances d and small values of β but the stability decreases when the parameters increase.
Figures 7(d) and (e) show our analysis of the pairwise interaction of the inner particle and one of the outer particles in the plane.We see that the distance d between the particles determines the strength of the interactions and that the choice of the propulsion mechanism via the value of β mainly determines the shape of the function T z (α).For an additional reference, we show the positions of the state transitions as red lines in figures 7(d) and (e).From the form of the curves it can be seen that the torque T z (α), which the inner particle exerts on one of the outer particles, can be described by a function with two fit parameters a and b.We use nonlinear least squares to obtain the values of the fit parameters for each combination of the values of d and β (see Tab. II in the Supplementary Information).For both parameters a and b, we find an inverse dependence on d (see figures 7(f) and (g)), which can be described by a function of the form 1/(c 0 + c 1 d).The parameter a does not depend on β, whereas the parameter b turned out to be proportional to β (see figure 7(h)).For the dependence of a and b on d and β we thus found the relations By a further nonlinear least squares fit, we obtained the following values for the fit parameters c a,0 , c a,1 , c b,0 , and c b,1 : Here, η is the viscosity of the liquid.Thus, the total torque that the inner particle exerts on one of the outer particles and its dependence on α, d, and β are given by the function To investigate whether the time-crystalline behavior also forms for other interactions, we visualize the dependence of the state of the system on the values of the parameters a and b from equation (15) in figure 8.The parameter a, in general, can take on values between 0 and 5.848PeηR 2 B 1 due to the limitation of d to the interval [2R, ∞ ].As the parameter b depends linearly on β and β is not restricted, b also is unrestricted.For figure 8, we added some further simulations that sample the parameter region where we observed a motion showing parallels to time crystals.Its occurrence depends significantly on both parameters.We found time crystalline-like behavior without fluctuations only for a > 2PeηR 2 B 1 and −4.3PeηR 2 B 1 < b < − 1.9PeηR 2 B 1 .Even though it would be interesting to investigate whether the periodic motion also occurs for a > 3.1PeηR 2 B 1 , this yields numerical problems because the space between the particles becomes very small then.This could be solved by increasing the mesh resolution to still be capable of solving the hydrodynamic equations, but this would further increase the computational expense of the simulations.Therefore, we leave this for future work.For all four states, we found extended regions, showing that the states are stable against small changes in the parameters determining the interactions between the particles.
For a further investigation of the influence of the interactions on the state of the system, we add a constant offset to the torque of the central particle.We set the parameters as d = 2.2R and β = − 3, which leads to timecrystalline behavior in the unperturbed case.Then, we apply external torques T ext to the central particle and determine in each case how the period length τ f of the oscillation changes relative to the period length τ 0 for no external torque which corresponds to an unperturbed motion.Like the torque exerted on the central particle by its interaction with the other particles, the external torque is perpendicular to the particle plane (both torques are parallel or antiparallel).When the central particle rotates in the opposite direction compared to the case without external torque, we add a negative sign to τ f .The behavior of the system considered in this work is invariant regarding the rotation of the whole system by an angle of 180°about either of the coordinate axes.Introducing an external torque, however, breaks this symmetry.Now, in case of a rotation about the x-or y-axis by 180°, the rotation direction of the particles is reversed.Thus, there are two equivalent possible setups when T ext is nonzero.These two setups lead to two curves describing the relation of T ext and (τ f − τ 0 )/τ 0 .The two curves are equivalent and can be mapped onto each other by a point reflection at (T ext , τ f ) = (0, 0).In the following, we study the properties of the state with parallels to a time crystal observed in figure 3 based on the variation of the torque applied on the central particle (see figure 9).Thereby, we focus on the curve with a clockwise rotation of the central particle for T ext = 0 (the blue curve in figure 9(a)), since the observations for the other curve are completely analogous.
For negative and decreasing T ext , τ f declines and asymptotically approaches τ f = 0. Also, for large and increasing T ext , the curve asymptotically approaches τ f = 0, but now it approaches from negative values of τ f .The asymptotic behavior for very large and very small values of T ext results from the fact that the external torque dominates the torques exerted by the other particles on the central particle in these limiting cases.Thus, it determines the rotation of the central particle for large |T ext |.When T ext is small and positive, τ f increases with T ext .For approaching intermediate positive values of T ext from above, the curve runs towards strongly negative values of τ f , whereas when approaching the intermediate values from below, the curve runs towards large positive values of τ f .This behavior is in line with the fact that at some intermediate external torque, the external torque cancels the torque exerted on the central particle by interactions with the other particles.Here, one would expect a divergence of τ f .The cancellation can be estimated to occur at T ext ≈ 35ηR 2 B 1 , since the torque exerted by the other particles on the central one for However, in the range between T ext = 20ηR 2 B 1 and T ext = 50ηR 2 B 1 , the behavior of the system is in fact much more complicated.In this range, the time-crystalline behavior changes qualitatively or vanishes.For values between T ext = 20ηR 2 B 1 and T ext = 31ηR 2 B 1 , we saw that the orientation of the central rotating particle is no longer within the particle plane.Its deviation from the plane increases with T ext (see figure 9(b)).We also observed periodic high-frequency fluctuations disturbing the angle between the orientation of the central particle and the plane.The standard deviation of this disturbance is shown in figure 9(b) as error bars.One can consider the behavior in this interval of T ext as a different state but still as similar to time crystalline behavior.Applying a larger external torque T ext = 33ηR 2 B 1 , however, leads to a completely different behavior.After a relatively long aperiodic initial relaxation, the particles rotate again periodically within the particle plane (see figure 9(c)).This state occurs at a point in the diagram that corresponds to a point lying on the second (orange) curve.Thus, applying this torque destabilizes the motion enough to flip the system to its other equivalent setup.Indeed, the central particle now rotates in the opposite direction than for lower values of T ext .The total torque acting on this particle is then 67.5ηR 2 B 1 , which equals approximately the externally applied torque T ext = 33ηR 2 B 1 plus the torque 35ηR 2 B 1 exerted on the central particle by the other ones when the system is not influenced externally.Thus, the point lies on the curve created by a point reflection of the blue curve at (T ext , τ f ) = (0, 0).From this reasoning, we can deduce that the behavior described by the blue curve looses its stability in this point.Near T ext = 35ηR 2 B 1 , the total torque acting on the central particle is approximately zero.Then, the motion is still periodic, but it has a rather large period length τ f and the motion within a period is rather complex (see figure 9(d) and supplementary video V5).Finally, in the region between T ext = 37ηR 2 B 1 and T ext = 50ηR 2 B 1 , the system shows an aperiodic motion (see figure 9(e)).This aperiodic motion is similar to the one observed in the coexistence state of the system.When T ext increases, the aperiodicity decreases until the system shows a periodic motion again for T ext = 50ηR 2 B 1 and larger external torques.However, in this case, the central particle rotates in the direction opposite to that for large negative external torques.In the region with T ext > 50ηR 2 B 1 , the external torque starts dominating the system's behavior again and a periodic state emerges.Though, it should be noted that this state is not a truly time-crystalline state as the periodicity is induced by the external torque.
To summarize our findings on the influence of adding a constant offset to the torque of the central particle, we found that for large absolute values this external torque is dominating.For smaller absolute values, we saw complex variations in the state of the system including a destruction of the periodic motion.In this analysis, we also found analogies to the mechanical properties of classical materials which can be visualized in stress-strain diagrams.See Appendix A for a more detailed description.
Thus, we can conclude that the interactions have a strong influence on the state of the system but periodic motion can be found for various different interactions.By investigating the interactions between the particles, we found a simple description of the induced torques that facilitates the transfer of the findings in this particular system to other systems.

Larger systems
Fourth, the time-crystallinity must also hold for larger systems.This criterion distinguishes a time crystal from a simple periodic motion, e.g., of a harmonic oscillator.We investigated this criterion by simulating interactions with additional crystallites.
Therefore, we extended the particle assembly described in section 3.1 to an infinite system.For this purpose, we reduced the size of the cubic domain in the y-direction and introduced periodic boundary conditions.The modified setup is shown in figure 10(a).To speed up the calculations and considering that the rotation of the particles in the system shown in figure 1 takes place completely in the plane, we restricted the swimmers' rotation to the plane.The results are presented in figure 10(b).We observed no qualitative changes compared to the simulation results described in section 3.1.However, we observed some changes in the way that the outer particles move.These changes are induced by the interactions with the newly added particles that create the larger system.However, the motion remains time-periodic even though the particles now interact with additional particles.This shows strong evidence that the observed behavior also holds in many-body systems.However, as periodic boundary conditions do not increase the degrees of freedom in the system, we do not simulate a true many-body system.For definitely stating that this criterion is fulfilled, simulations with many more particles would be needed.Due to the high computational cost for calculating the hydrodynamic interactions between the particles, we were not able to perform these simulations.Nevertheless, our simulations provided evidence that the periodic behavior also holds in many-body systems.
These considerations together with the ones in the previous sections show in which aspects the time-periodic state of the considered active colloidal crystallite resembles a time crystal.

Conclusions
We observed a periodic motion of an active colloidal crystallite, compared it to a continuous time crystal, and analyzed the occurrence and properties of this exceptional state of matter.Such colloidal systems form a fascinating new class of materials that deserves extensive exploration in future research.
Our investigation revealed requirements for the occurrence of a colloidal time crystal.It showed that an assembly of squirmers as studied in this work self-assembles such that it shows parallels to a classical continuous time crystal when the values of the particle distance and squirming parameter are sufficiently small.We found that the occurrence of this state is neither prevalent nor limited to a selection of very particular parameter values, that it is robust against fluctuations originating from numerical errors or Brownian motion, that it strongly depends on the interactions of the particles, and that it occurs if interactions as in larger systems are included.Considering different particle assemblies, we found that the spatial assembly of the particles has a strong influence on their hydrodynamic interactions and thus on their potential formation of a time-crystalline state.
In the future, this work should be continued by additional simulation studies that consider the dependence of the behavior of the system on the parameter values in greater detail, calculate the state diagram with higher resolution and for larger ranges of the parameter values, and take also qualitatively different particle assemblies and shapes into account.In the case of different assemblies, it would be interesting to study various regular and irregular spatial assemblies of active particles in two and three spatial dimensions.For investigating the influence of the particle shape, the presented model can be applied directly.However, the computations are expected to take longer as the geometry for the flow field generation changes for each individual time step.Our analytic representation of the hydrodynamic interactions of two particles provides an ansatz for replacing the computationally expensive hydrodynamic simulations performed in this work by much more efficient Brownian dynamics simulations that could strongly speed up future research.With equation (22) that we provided for the interaction torque, it should be possible to perform Brownian dynamics simulations for an infinitely large hexagonal crystal of active particles.For other particle assemblies, including those studied in the present work, a similar approach should be possible, but in this case it would be necessary to first determine also fit functions for the interaction torques exerted by a particle at the boundary of the crystallite on other particles.
While we identified different state transitions of the system, a deeper study of the nature of these transitions can be performed using large deviation theory [56,57].This theory, that has been used to study other nonequilibrium transitions [58,59], can reveal further insights and allow for finding parallels to other systems.Additionally, an analytical investigation of the system investigated in this work would be beneficial to obtain a better understanding for the interactions between the particles.This might not only provide a deeper understanding of the influence of the spatial assembly of the particles but also facilitate a better understanding of the propulsion mechanism's influence.
Another important future step in the investigation of colloidal time crystals is their experimental realization.For this, our numerical work lays a good basis by providing guidelines how a system of active particles must be designed to be able to exhibit time-crystalline behavior.The realization of the system studied in the present article is possible by arranging and holding biological or artificial microswimmers on a hexagonal lattice with optical or acoustical tweezers [60].Particles that generate suitable flow fields are those that are well described by the pusher model.This is the case for microorganisms like Escherichia coli or artificial microswimmers like Janus particles [61].Recently, it was also shown that particles driven by ultrasound create a flow field similar to the one described by the pusher model [62].These particles allow an easy tracking of the particle directions, e.g., by observing the system in a microscope [63][64][65].Regarding the assembly of the particles, we propose to use acoustic tweezers, since they are not affected by the inhomogeneous optical properties that are common for active particles.Light fields should preferably be used to supply light-propelled active particles by energy [66].The spatial configuration of the microswimmers that we studied in our work could also be achieved by selfassembly when the particles are equipped with suitable interactions.An option to equip artificial active particles with repulsive interactions that can cause a hexagonal particle assembly is to embed a magnetic core within the particles [67].By varying the strength of an external magnetic field, which is oriented perpendicular to the particle plane and induces a magnetic dipole in the particles, it is then possible to tune the repulsive interactions of the particles and thus the lattice constant.

Figure 1 .
Figure 1.Self-organized time-periodic motion of active colloidal particles.a If the particles are pushers with radius R, distance d = 2.2R, and squirming parameter β = −3, they self-organize an orientational configuration that behaves periodically in time and shows parallels to a colloidal time crystal.The snapshot shows the particle orientations ûi with i = 1,K,7 and the flow field v around the particles at time t = 475R/B 1 , where B 1 denotes the first velocity mode.b As the time evolution of the particle orientations shows, the central particle (i = 1) rotates clockwise with a constant speed, whereas the other particles oscillate periodically.

Figure 2 .
Figure 2. PCA projections of the 21-dimensional trajectory u(t) describing the whole system.The axes of each subplot range from −2 to 2 and the numbers on the diagonal indicate the principal components.χ denotes the percentage of variation covered by the chosen number of eigenvalues.a The time-periodic state is described by a two-dimensional circular trajectory.b When d increases, the trajectory leaves the two-dimensional plane.c For even larger values of d, no order is visible and more components occur.d The trajectory describing a crystallite of neutral squirmers (β = 0) also shows no order.

Figure 3 .
Figure 3. State diagram of an active-particle system showing a region with parallels to a time crystal.a The state diagram of the system from figure 1 for variable particle distance d and squirming parameter β involves four different states including one exhibiting parallels to a time crystal.b-e For each state, a representative part of the time evolution of the particle orientations is shown.The supplementary videos V1-V4 show the time evolution of the flow field for each of the four states.The disturbed periodic state shows fluctuations but an approximately regular pattern (c).The coexistence state is characterized by time spans that are similar to the disturbed periodic state but also entirely aperiodic time intervals (d).

Figure 4 .
Figure 4. Influence of the last simulated point in time t max on the order parameter d.We show the four states with different distances d/R from figures 3(b)-(d) and observe a decrease of d for increasing t max with little variation for > t R B 350 max

Figure 5 .
Figure 5. Transitions between the states from figure 1(a).The red disks represent the order parameter d while the heatmap shows the distribution of the values δ(t).a Variation of the particle distance d for β = − 3. The transition at d = 2.32R has bare analogies to a second-order phase transition and the transition at d = 2.51R has bare analogies to a first-order phase transition.At d = 2.95R, we found by visual inspection a change from the coexistence state to the fully aperiodic state that is not visible in the order parameter d. b Variation of the squirming parameter β for d = 2.2R.The system's behavior has bare analogies to a first-order phase transition at β = − 2. Note that at β = − 0.5 the order parameter vanishes due to a static state of the system.

Figure 6 .
Figure 6.Influence of Brownian motion.All plots correspond to the observed time-crystalline state with parameter values d = 2.2R and β = − 3. a The Euclidean distance d T between the undisturbed trajectories and the trajectories with Brownian motion decreases for an increasing Péclet number.To eliminate shifts in time, dynamic time warping [54] was applied before calculating the distance.bThe trajectories of single particles for Pe = 100 (not included in (a)) are strongly disturbed but the central particle (particle 1) shows a similar motion as in the undisturbed case.c, d For increasing Pe, the influence of the Brownian motion becomes smaller and the motion of the particles approximates the undisturbed motion.

Figure 7 .
Figure 7. Analysis of interactions.Unless stated otherwise, the parameters are d = 2.2R and β = − 3. a Torque T z exerted by the propulsion of an outer particle (particle 7) on the other particles and itself.The nearest neighbors (1, 2, 6) are strongly influenced whereas the second nearest neighbors (3, 4, 5) only experience very small torques.b, c Torque T rot that causes the inner particle (particle 1) or an outer particle, when its orientation deviates by 2°from the plane, to rotate itself back to the plane (positive) or further away from the plane (negative) as a function of the parameters d and β.The curves represent the torque averaged over all in-plane orientation angles α and the shaded areas show the corresponding range of values.Note that these curves only include the values of α that occurred in the simulations for figure 5.The vertical lines indicate boundaries of states and correspond to those from figure 5. d, e Torque T z exerted by the propulsion of the inner particle on one of the other particles as a function of the angle α and for different values of the parameters d and β. f-h Dependence of the fit parameters a and b from equation (15) on the parameters d and β for the functions shown in (d) and (e).

Figure 8 .
Figure 8. Occurrence of motion with parallels to a time crystals depending on particle interactions.The data show the values of the parameters a and b from equation (15) and correspond to the simulations for figure 5 and additional simulations.

Figure 9 .
Figure9.Properties of a time-periodic active crystal when applying an external torque.a The relationship between the period length τ f of the oscillation of the dynamic state and various external torques T ext applied on the central particle shows for which external torques the periodic motion is destroyed.b In the interval marked in a, the orientation û1 of the central particle deviates from the particle plane.c-e A representative part of the time evolution of the particle orientations ûi with i = 1,K,7 is shown for T ext = 33ηR 2 B 1 (c), T ext = 35ηR 2 B 1 (d), and T ext = 37ηR 2 B 1 (e).

Figure 10 .
Figure 10.Time periodicity in larger systems.a Setup for the simulations of an infinite system.b Corresponding simulation results for d = 2.2R and β = − 3.