Coexisting crystal and liquid-like properties in a 2D long-range self-consistent model

A two-dimensional class of mean-field models serving as a minimal frame to study long-range interaction in two space dimensions is considered. In the case of an anisotropic mixed attractive-repulsive interaction, an initially spatially homogeneous cold fluid is dynamically unstable and evolves towards a quasi-stationary state in which the less energetic particles get trapped into clusters forming a Bravais-like lattice, mimicking a crystalline state. Superimposed to this, one observes in symplectic numerical simulations a flux of slightly more energetic particles channeling through this crystalline background. The resultant system combines the rigidity features of a solid, as particles from a displaced core are shown to snap back into place after a transient, and the dynamical diffusive features of a liquid for the fraction of channeling and free particles. The combination of solid and liquid properties is numerically observed here within the classical context. The quantum transposition of the model may be experimentally reached using the latest ultracold atoms techniques to generate long-range interactions.

at minimal energy densities, the system organizes in a Bravais-like lattice forming cores playing the role of the atoms in a crystal structure. In the present study, the energy density of the system is not minimal yet the initial temperature is vanishingly small and we use molecular dynamics simulations. As the system follows its natural evolution, it slightly warms up and organizes into a 'cold' crystal structure that is gone through by a coherent flux of the most energetic particles. After introducing the physical model, evidence is given of the dual solid and liquid nature of the quasi-stationnary state. Its quantum realization could be produced with cold atoms driven by laser light where multiple scattering of photons by atoms gives rise to infinite-range forces when the atoms couple to a single mode high-finesse cavity 19 .

The Two-Dimensional Hamiltonian Mean-Field Model
Derivation. The 2D-Hamiltonian Mean Field (HMF) model was first proposed by Antoni and Torcini 20,21 as a two dimensional generalization of the Hamiltonian Mean Field (HMF) model in the fully attractive case for the study of N-body self-gravitating systems. A generic two-body potential in a two dimensional square box of side 2π with periodic boundary conditions can be Fourier expanded as x y Retaining only the most long-range terms, with |k| = 1 and = k 2 , in the expansion and requiring that the potential be invariant under rotations of multiples of π/4 22 yields the following truncation of the potential in Eq. (1) ( cos cos ) cos cos , (2) where a is an arbitrary scaling constant, c and d are coupling constants, and due to the rotation invariance c is related to the energy scaling and d is the only free parameter 22 .
Considering N particles interacting through the potential in Eq. (2) and setting 2c + d = −a, one gets the following Hamiltonian Here the p ix and p iy are the conjugate momenta to the space positions x i and y i . The 1/N prefactor corresponds to the Kac's prescription 23 . It recovers the extensivity of the pair potential and is equivalent to a time rescaling of the type ′ = t N t. The first term in Eq. (5) is the potential of two uncoupled (one-dimensional) HMF models, while the second term couples the x and y directions.
Defining the four self-consistent mean-field variables as N N with M 1 and M 2 playing the role of the magnetization fields and P + and P − of the polarization fields and making use of the notation 〈⋅〉 N for the average over the N particles, the Hamiltonian in Eqs (3-5) simply reads The dynamics of any single particle i can be easily shown to obey the equations of motion At any time t, the equations (11 and 12) are obtained from the following one-particle Hamiltonian (defined up to a constant of integration C) The constant of integration can be fixed by using the conservation of the total energy yielding There is another constant of motion in this system: the total momentum In the numerical simulations, this is taken to be identically zero.
Equilibrium and out-of-equilibrium emergence of a Bravais lattice structure. The phase space trajectories of a Hamiltonian system such as Eq. (10) are constrained on a constant energy surface in phase space. Consequently, the time averages computed from the numerical solutions of the equations of motion are expected to converge to microcanonical ensemble averages. For some long-range interacting systems, the ensembles may not be equivalent and so it turns out that the microcanonical ensemble is the natural ensemble to derive equilibrium statistical mechanics.
In a previous work 18 , we studied the mixed case with attractive polarization mean-fields and repulsive magnetization mean-fields by choosing c = −1 and d = 1. Due to the invariance of the Hamiltonian under rotations by multiples of π/4, a simple interchange in the "charge" of the fields would provide the same results. The equilibrium statistical mechanics studies led to the following results obtained in the microcanonical ensemble: along the repulsive direction, the system behaves as the antiferromagnetic-like HMF model, whereas along the attractive directions, it behaves as the ferromagnetic-like HMF model, except that for the low energy phase, it exhibits a bicluster, instead of a single cluster structure. Due to the periodicity of the potential, one can increase the space period length by a multiple of 2π, and turn the bicluster into a periodic structure which can be regarded as a Bravais-like lattice. This space ordering into a Bravais lattice structure is not only present in the equilibrium states but also in the quasi-stationary states (QSS). Indeed, since the model (10) is long-range, its dynamics, computed here using a fourth-order symplectic scheme 24 , exhibits ergodicity breaking features. As already observed, e.g. in the HMF model, after an initial phase of violent relaxation, the system settles in an out-of-equilibrium quasi-steady state, the lifetime of which diverges with N. Figure 1 displays the typical Bravais lattice structure emerging in the QSS phase starting from cold and space homogeneous initial conditions. Contrarily to Monte Carlo simulations 18 that do not capture the real dynamics of the system, in molecular dynamical simulations one observes that some particles are hopping from one cluster to the others so that there is a flux of particles hopping from one cluster to another, the amount of the flux depending on the system energy.
We shall now report the results on the formation of the Bravais lattice structure and QSS's evolution and characteristics in the mixed attractive/repulsive 2D-HMF model.

Dynamical Features of the Attractive/Repulsive 2D-HMF Model in the Low-Temperature Regime
Initial conditions. Our previous results on this system 18 showed that the homogeneous distribution of a cold ensemble of particles is linearly unstable in the attractive directions and linearly stable in the repulsive directions. The resulting linear instability triggers the so-called violent relaxation process, according to Lynden-Bell's wording [25][26][27] . The initial distribution functions considered in the present analysis are the following water-bag distribution functions where r ≡ (x, y), p ≡ (p x , p y ) and Θ is the Heaviside unitary step function. Inasmuch as the mean square momentum for distribution (20) is given by the distribution is considered to be an ensemble of cold particles if Δp is chosen to be sufficiently small. The total energy is given by In the present Hamiltonian system, the energy is conserved and equal to , whereas the polarization mean-fields are . + −   P P 0 6. These values are consistent with the results presented in the Fig. 8 of ref. 18 . The latter shows the linearly unstable and subsequent saturation stages of the time evolution of the moduli of the mean-fields. Because of the symmetry in x and y in the expression of the Hamiltonian (10) and of the spatial homogeneity, and therefore invariance, in x and y of the initial conditions, the mean-fields are approximately equal in the N ≫ 1 limit, namely ( ≈ = ). In the three-dimensional plots of Fig. 2, each of the N = 25000 particles is represented by a dot as the function of its positions in the 2D [0:2π] × [0:2π] cell and of its energy, the color of the dot depending on the particle energy. The bottom plots are projections on the horizontal plane that help to visualize the location of the cores, namely the instantaneous location of the crystalline structure. Indeed the cold mixed attractive/ repulsive 2D-HMF model is a many-body system in which the lowest energy is a state of modulated density: this low-energy fraction of particles defines a crystal in the sense of Landau 28 . Figure 4 shows the energy distribution of particles at two different stages, in the QSS regime and later in the thermalization stage. The energy distributions were obtained from averaging over the data obtained at 40 equidistant instants within a time range of 80 time units about the times t = 2000 and t = 10 4 .
Considering the particle energy distribution in the QSS regime (about t = 2000), it is possible to distinguish three types of particles: i) low energy particles (LEP), appearing in blue on Fig. 2, constituting the majority of the system particles and forming the bulk of the energy distribution function, which stay trapped in the cores of the bicluster structure; ii) a second peak around ε = −0.3 reveals a second class of particles constituted by intermediate energy particles (IEP), these particles can hop from one core to another, being still attached to the mean-fields potential, moving coherently in two well defined directions (along the diagonals) and, finally, iii) a few high energy particles (HEP), appearing in red and yellow in the snapshots of Fig. 2, that can move freely in space.
Potential topology. It is useful to figure out the behaviour of the potential energy associated to the one-particle Hamiltonian (17). This is represented in Fig. 5 in the limit case of the Vlasov limit N → ∞ inducing   Deviation from thermodynamic equilibrium. Let us now quantify the deviation from thermodynamic equilibrium of the QSS state. A simple way to check whether the system has already reached or is going toward the thermodynamic equilibrium is to compute the reduced kurtosis, k, defined as This is the fourth standardized moment minus 3. For a Gaussian distribution, which characterizes the Maxwell-Boltzmann equilibrium, one has k = 0. The graph in Fig. 6 plots the evolution of the reduced kurtosis as a function of the time divided by the number, N, of particles. It indicates that the thermalization timescale is proportional to N. From this follows, using a proof by contradiction, that the lifetimes of the QSS are, at most, diverging as N. Such a scaling would be corroborated by the results of Chavannis 29 .
These are derived from kinetic equations for bidimensional long-range interacting systems and indicate a lifetime for QSSs that scales linearly with N. This scaling with N is also obtained from a stochastic modeling of the thermalization process involving the disintegration of coherent structures sustaining out-of-equilibrium quasistationary states in the one-dimensional attractive HMF model 6 .

Combination of Solid and Liquid-Like Features
Heterogeneous diffusive properties. The fact that the QSS combines solid and liquid-like features is first evidenced by the heterogeneity of its diffusion properties. Figure 7 presents the computation of the mean square displacement (MSD) for each kind of particles during the QSS regime. The MSD plotted in Fig. 7 is defined as   where the brackets 〈 ⋅ 〉 N stands for the average over N particles. The results shown were obtained for t 0 = 2500, n = 10 and for a total number of 10 4 particles. The choice for t 0 = 2500 was in order to ensure that the system had already gone into the QSS regime and formed the periodic structure. The time evolution of the MSD for the core particles shows that they have a vitreous profile, characterized by time intervals with almost no diffusion. This can be interpreted as the Bravais-like lattice having a high viscosity. Meanwhile, the diffusion profile in Fig. 7 and the snapshots of the QSS shown in Fig. 2 show the movement of the cores without changes in the periodic array of the Bravais-like lattice. This global movement of the crystal structure results from the requirement of total momentum conservation so that the heavy cores move slowly in order to counterbalance the rapid movement of the few flux particles.
Conversely, the IEP, which are the particles that support the flux between the cores, and the free particles (HEP) present almost the same profile for diffusion, markedly different from that of the particles (LEP) forming the cores. The diffusion regime for IEP and HEP appears first as diffusive then as sub-diffusive at large time. The latter behaviour is however an artefact due to the fact that the system is confined into a finite square box of sides 2nπ by 2nπ, in such a way that there exists a maximum value that the mean square displacement can reach. Figure 7 also shows that the flux particles, even trapped to the mean-field potential of the cores, diffuse in a way independent from the vitreous diffusion of the cores. Figure 7 reveals an heterogeneous diffusive behaviour: particles forming the cores of the crystal structure have a glassy behaviour whereas the rest of the particles diffuses normally as in a normal fluid.  Moreover, for the initial conditions under consideration, namely space-homogeneous with vanishing temperature and total momentum, the collective momentum of these diffusive particles flowing through the crystal-like structure is the opposite of the drifting momentum of the clustered particles, namely −P Core (t). The time behaviour of the core total momentum vector has been plotted in Fig. 8. This means that, in the reference frame of the crystal, the particles forming the liquid-like phase have a nonzero mass flow.

Response to External Perturbations
Finally, in order to test the solid character of the QSS, we studied the response of the QSS to external perturbations. A core was displaced at some given time as represented in the left plot of Fig. 9. As visible from the central plot displaying the time behaviors of the magnetization mean-fields, particles from the displaced core snap back into place after a transient oscillatory stage. Indeed, in the case of Fig. 9, all the LEP particles from the left bottom core are displaced at some given time to the right along x. This produces an increase of M 1 . The core then starts to oscillate along x, yet with a decreasing amplitude, about its original position that it recovers at time  t 500 which is captured by the damped oscillatory behaviour of M 1 . This phenomenology is a solid feature in contrast to a fluid that would be permanently displaced.

Final Remarks and Paths to Quantum Realization
The results presented here indicate that, in the N → ∞ (Vlasov) limit, the mixed attractive/repulsive 2D-HMF in the low-temperature regime remains frozen in the QSS having a periodic Bravais-like structure with a flux of particles with non-zero mass flow between the cores. According to previous Monte Carlo results in the microcanonical ensemble 18 , there exists a threshold in the energy density above which the mixed 2D-HMF is homogeneous at equilibrium. Moreover former linear theory calculations for the stability of waterbag-like initial conditions 20 indicate that the homogeneous state will no longer be unstable above some threshold in the initial temperature, forbidding the violent relaxation stage. This leads us to predict the existence of some threshold in the initial temperature above which the QSS reported here does not survive and the system remains quiescent in the homogeneous state.
In the limit of vanishing temperature, that is valid at least in the initial stage of the dynamics considered in this work, the thermal de Broglie wavelength  λ π = mk T (2 /( )) th B 1/2 (with mass m = 1 in our model) becomes larger than the interparticle distance and quantum effects are significant. The quantum transposition of this model would therefore be interesting to investigate. In order to prepare for a quantum transposition of the model, mean-field quantum descriptions have been introduced in the Supplementary material. In the case of bosons, a linear theory is presented in the zero temperature limit. This enables to fix the constant parameters in front of the mean-fields in the expression of the potential for which linear instability (violent relaxation) exists in the quantum regime. This secures the existence of QSS, due to the combination of nonlinear effects leading to the Figure 8. Instantaneous total momentum P Core of the particles trapped in the self-consistent potential wells as a function of time for N = 25000 particles (same parameters as in Fig. 2). The fraction of trapped particles forming the crystal-like structure remains equal to about 70% during the whole QSS stage for the initial conditions (20) considered in the present study. saturation of mean-fields and of ergodicity failures for quantum long-range systems 14 . It remains to be studied what is the impact of the nonlinearities entering through quantum effects (in the quantum pressure term) on the nature of the QSS and possibly modify the potential to counteract this effect, if needed. The case of fermions may display the richer phenomenology since, along the attractive direction, one may obtain a transition between a Bose-Einstein condensate with a superfluid phase for strong interaction and a supraconductor state for weak interaction. These last conjectures are left for further studies. Finally, the interaction range involved in the 2D-HMF model is infinite. As for the experimental quantum realization, it needs then to involve interactions that have a range longer than that of the dipole interaction having a space dependence as 1/r 3 . There is presently a reckless quest in the quantum experimental community to access infinite range interactions. This may be attainable using optical cavities 19 .
Finally, compared with the traditional egg-crate 2D potential represented in the right part of Fig. 5, the potential resulting from the mixed attractive-repulsive interaction, has two barriers in the potential well, one for the confinement of particles forming the crystal-like structure and another one for the IEP fluid-like particles. This offers more freedom to control and access delocalization.