Introduction

Structures of the observable Universe, such as galaxies and globular clusters, appear to be macroscopically stationary but are not at thermodynamic equilibrium, i.e., the distribution of the velocity of their stars is not Maxwellian1. Indeed, Chandrasekhar pointed out in 1941 that the time necessary for these objects to reach thermal equilibrium is actually much larger than their age2. This has been confirmed by observations determining that these astrophysical structures are indeed far from thermal equilibrium3. In 1967 Lynden-Bell proposed a mechanism, violent relaxation that leads to the formation of these out-of-equilibrium structures, called quasi-stationary states. These structures evolve towards the quasi-stationary state on a time much faster than that required for full thermodynamic equilibrium4. Importantly, it has been subsequently understood that this mechanism is generic in Hamiltonian systems with a long range interacting potential, i.e., a potential that is not integrable as a result of its extension over large scales5. This phenomenon is similar also to what arises in plasmas subject to Landau damping, in which there is an exchange of energy between the electromagnetic wave generated by the particles of the plasma and the particles themselves6. Landau damping has been observed in plasma experiments7,8,9,10,11,12,13 and in space plasma turbulence14. Contrary to Landau damping, violent relaxation is more elusive and has not been observed to date, neither in a repeatable or controllable experiment, nor in situ. Indeed, experimental observation of the dynamics of the formation of quasi-stationary states via violent relaxation is hindered mostly for two reasons. First, there are systems in which it is potentially present but it is destroyed by the stochastic noise generally present in these systems. An example of such systems are cold atoms confined in optical traps, in which the noise is produced by the interaction of the photons of the lasers with the atoms15,16. Second, there are systems in which violent relaxation is actually present, but the associated timescales are too large to observe it. This is the case of astrophysical systems such as galaxies, independently if it is constituted by classical (non-quantum) dark matter particles17,18,19,20, or composed by quantum matter commonly called fuzzy dark matter21,22,23,24,25. In these systems violent relaxation occurs on time scales of the order of millions of years1.

Violent relaxation has however, been studied numerically for example, by simulating classical N-body systems with nonlocal (e.g., gravitational) interactions that are governed by Vlasov-Poisson equations26,27,28,29. Whilst these simulations provide confirmation of the process proposed by Lynden-Bell, they provide no guidance of how to observe violent relaxation experimentally.

Here we report the experimental observation of violent relaxation in an optical setting where we observe the evolution of an optical beam in the presence of a self-generated long-range interaction that leads to phase-mixing in the presence of a varying potential and final relaxation towards a quasi-stationary state. Furthermore, the observed optical dynamics can be viewed as an analogue of (dark matter) galaxy formation via violent relaxation. Indeed the connecting point is the underlying Vlasov-Poisson equation—this describes the dynamics and violent relaxation of an N-body system of self-gravitating particles of dark matter30,31. The Vlasov-Poisson equation is also the semi-classical limit of the quantum description of dark-matter evolution based on the Newton-Schrödinger equation (NSE)—the latter therefore also describes the evolution of classical dark matter32 and importantly, it is also the same equation that describes our optical experiments. The NSE has been experimentally realised in nonlinear optical experiments that have been used to probe gravitational lensing, tidal forces and analogue quantum processes such as Boson star evolution33,34. By choosing the appropriate system parameters, we show that it is possible to experimentally observe violent relaxation and the formation of a quasi-stationary state in the form of a table-top galaxy that bears a close resemblance to the result of an N-body numerical simulation.

Results

Self-gravitating systems

The temporal evolution of self-gravitating particles of dark matter, of mass m, defined by a wavefunction ψ(r, t), is described in 3D by the Newton–Schrödinger equations (NSE):

$${{{{{{{\rm{i}}}}}}}}\hslash {\partial }_{t}\psi +\frac{{\hslash }^{2}}{2m}{\nabla }^{2}\psi +m\phi \psi =0$$
(1a)
$${\nabla }^{2}\phi =-4\pi G| \psi {| }^{2},$$
(1b)

where ∣ψ∣2 is the mass density, G the gravitational constant and ∇2 the three-dimensional (3D) Laplacian. The gravitational potential, ϕ, generated by the mass distribution itself, depends on the constant G and the mass density. When the system is in the semi-classical regime, which corresponds to ℏ/m ≪ 1, the process of violent relaxation can be observed. This process leads the system towards its quasi-stationary state1.

The features of violent relaxation

Violent relaxation consists in the evolution of the energy distribution due to the variation in time of the potential ϕ(r, t). The evolution towards the quasi-stationary state is accompanied by phase mixing due to the evolution of the wavefunction in the non-harmonic potential ϕ(r, t). We underline that mixing alone is a relaxation process by itself i.e., when violent relaxation is present, mixing is also generically present, whilst the opposite is not generally true as mixing may occur in the absence of a time-varying potential. Indeed, the concomitant presence of mixing and a time-varying potential leads to significantly accelerated and abrupt relaxation dynamics, hence the naming ‘violent’ (see also Supplementary Note 3). Therefore, any demonstration of violent relaxation requires as a minimum, the presence of both of these ingredients.

In the semi-classical regime, the quasi-stationary solution for violent relaxation process corresponds to the formation of an oscillating solitonic core in the center of the system (defined as the ground state of Eq. (1),35,36,37) surrounded by the stationary solution of the classical Vlasov-Poisson equation, which is also the semi-classical limit of the NSE (i.e., the limit ℏ/m → 0 of Eq. (1))22,32. In order to be in the appropriate regime, the soliton should to be small compared to the size of the whole system. When this happens, the system can be considered to be sufficiently semi-classical (ℏ/m → 0) to observe violent relaxation. Therefore, we monitor the degree of classicality with the parameter χ = ξ/s, where ξ is the characteristic size of the soliton. ξ can be estimated by calculating the scale for which the kinetic and potential energies are of the same order of magnitude, giving ξ = ℏ2/(8πGMm2) in the case of a 3D gravitational system and where M is the total mass, and s the size of the whole system. The goal then is (differently from previous studies looking at soliton evolution) to study the behaviour of the broad semi-classical background solution, corresponding to the galaxy, in the potential generated by the total field. In order to do this, we take advantage of the formal identity of Eq. (1) to that describing our optical system.

Optical system

The optical system is based on a slab of glass or crystal that exhibits strong thermo-optical nonlinearity. Briefly, when an intense laser beam propagates inside the crystal, it induces a nonlocal interaction (heating) of the medium. This heat profile in turn acts on the laser beam with a focusing action—this therefore can emulate for example the nonlocal self-attractive gravitational force. In the paraxial approximation, the propagation of a monochromatic laser beam with amplitude, \({{{{{{{\mathcal{E}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{\perp },z)\), in a thermally focusing nonlinear medium is described by33,34,38,39:

$$\begin{array}{l}{{{{{{{\rm{i}}}}}}}}{\partial }_{z}{{{{{{{\mathcal{E}}}}}}}}+\frac{1}{2{k}_{0}{n}_{{{{{{{{\rm{b}}}}}}}}}}{\nabla }_{\perp }^{\,2}{{{{{{{\mathcal{E}}}}}}}}+{k}_{0}\Delta n{{{{{{{\mathcal{E}}}}}}}}+{{{{{{{\rm{i}}}}}}}}\,\frac{\alpha }{2}\,{{{{{{{\mathcal{E}}}}}}}}=0,\\ {\nabla }_{\perp }^{\,2}\Delta n=-\frac{\alpha \beta }{\kappa }| {{{{{{{\mathcal{E}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{\perp },z){| }^{2},\end{array}$$
(2)

where r⊥ = (x, y) is the two-dimensional (2D) position in the plane transverse to the propagation direction z. The operator \({\nabla }_{\perp }^{\,2}\) is the transverse 2D Laplacian, \({k}_{0}=\frac{2\pi }{\lambda }\) the wave-number of the incident laser with nb the background refractive index of the medium. The non-local nonlinear refractive index change, Δn, is induced by the beam itself heating the medium. β is the medium thermo-optic coefficient, κ its thermal conductivity and α its absorption coefficient. The last term of Eq. (2) accounts for the absorption in the crystal and in our parameter space, has little effect on the violent relaxation dynamics40 (see discussion in Supplementary Note 6). We hence neglect it in the following discussion.

Provided that z plays the role of time t, the similarity between Eqs. (1) and (2) underpins the opportunity to directly observe 2D violent relaxation in a laboratory experiment. Especially, the presence of violent relaxation is independent of the dimension of space, and occurs also in 1D41 and, as in our case, in 2D42. The main difference between the 2D optical and 3D gravitational systems is the shape of the potential being logarithmic in the 2D system. However, the mechanism that governs the physics in both systems is the same, i.e., modes or particles that live in a confining potential that is evolving in time can undergo mode mixing and violent relaxation, leading to the formation of a quasi-stationary state. The optical equivalent of the above-mentioned semi-classical regime is obtained when χ = ξ/s ≪ 1. In the optical case, \(\xi =\sqrt{{z}_{{{{{{{{\rm{nl}}}}}}}}}/(2{k}_{0}{n}_{{{{{{{{\rm{b}}}}}}}}})}\) is the soliton size, defined as the transverse length scale for which both the linear and nonlinear effects are of the same order. znl = κ/(αβk0P) is the longitudinal length over which the effect of the nonlinear term becomes substantial and \(P=\int\,d{{{{{{{{\bf{r}}}}}}}}}_{\perp }{\left\vert {{{{{{{\mathcal{E}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{\perp },z)\right\vert }^{2}\) is the power of the laser beam. The initial beam with transverse width s dictates the propagation regime of the system. We note that we do not work in the soliton-dominated regimes studied in the past. There, the waist s was typically chosen to be close to the soliton waist ξ, so as to generate as little background (i.e., surrounding radiation that is not converted into the soliton) as possible38,43. Instead, the present work relies on the formation of a significant, surrounding ‘background’, i.e., the analogue galaxy, hence with s ≪ ξ.

We define the following spatially dependent quantity for the optical system

$${{{{{{{\mathcal{U}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{\perp },z)=\frac{{\left\vert {\nabla }_{\perp }{{{{{{{\mathcal{E}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{\perp },z)\right\vert }^{2}}{2k{\left\vert {{{{{{{\mathcal{E}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{\perp },z)\right\vert }^{2}}-{k}_{0}\Delta n({{{{{{{{\bf{r}}}}}}}}}_{\perp },z),$$
(3)

where k = k0nb. The first term is a kinetic (linear) energy density \({{{{{{{\mathcal{K}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{\perp },z)\), the second term is a potential (nonlinear) energy density \({{{{{{{\mathcal{V}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{\perp },z)\). \({{{{{{{\mathcal{U}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{\perp },z)\) is the energy needed to add or remove a particle from the system and is therefore effectively a spatially dependent chemical potential. Particles will tend to reduce the free energy of the system by moving from high to lower chemical potential regions44 and therefore, in keeping with the N-body analogy, \({{{{{{{\mathcal{U}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{\perp },z)\) quantifies the redistribution of particles and energy due to violent relaxation.

In order to then characterise and quantify violent relaxation in optical experiments, we consider two quantities: firstly, we define the evolution of the distribution of the chemical potential density \(\nu ({{{{{{{\mathcal{U}}}}}}}})\) of the optical field \({{{{{{{\mathcal{E}}}}}}}}\) (see Supplementary Note 1 for details), that captures the main signature of the violent relaxation process, i.e., the change in the distribution of the energy due to the variation in the potential \({{{{{{{\mathcal{V}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{\perp },z)=-{k}_{0}\Delta n({{{{{{{{\bf{r}}}}}}}}}_{\perp },z)\) along z. Secondly, we consider the Wigner transform45 F(r⊥, k⊥, z) of the optical field \({{{{{{{\mathcal{E}}}}}}}}\), i.e., the density of probability to find a portion of the optical beam at the position r⊥ with wavevector k⊥. We use the evolution of F with respect of z to study the mixing of the phase-space.

Experimental setup

Figure 1a shows a schematic representation of the experiment. A continuous-wave laser beam with a Gaussian profile and wavelength λ = 532 nm propagates in a thermo-optical nonlinear medium made of three aligned identical slabs of lead-doped glass for a total length L = 30 cm, represented here as a single slab.

Fig. 1: Overview of the experimental setup and results.
figure 1

a Sketch of the experiment. A Gaussian laser beam propagates in a lead-doped glass slab. Fully detailed experimental layout is shown in SM. The diffusion of heat inside the nonlinear medium is represented by the glowing red profile. Insets show measured input and output experimental light intensity profiles (analogues of particle density) at power P = 5 W. b y = 0 slice of the beam intensity profile as a function of one transverse coordinate x and power, obtained from the numerical simulation (the fields have sufficient azimuthal symmetry for these lineouts to be representative of the evolution). c Simulation of the full transverse plane distribution, \({\left\vert {{{{{{{\mathcal{E}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{\perp },z)\right\vert }^{2}\) at z = L for input power P = 5 W. We observe the soliton (red dot indicated by the black arrow) surrounded by the classical part corresponding to χ → 0. d y = 0 slice of the beam intensity profile as a function of one transverse coordinate x and power, obtained from experimental data. e N-body simulation result under the same conditions, performed evolving self-gravitating particles of dark matter. We have used 217 N-body particles and the parameters map directly on to those used in the experiments (see Supplementary Note 2 for details), i.e., this galaxy is the particle version of the optical galaxies shown in a (experiment) and c (numerical simulation). All simulations represent the optical field evolution (based on Eq. (2)) with the exception of e that shows the result of the equivalent N-body particle simulation (see Supplementary Note 2)). All simulations are implemented in Matlab or C (see Supplementary Note 2).

Full size image

The beam width s = 350 μm at the sample input facet is chosen experimentally by a system of lenses such that \(\chi =2.3\times 1{0}^{-2}/\sqrt{P}\), with P measured in Watts and is of order 10−2 ≪ 1 over the full evolution. This therefore ensures that we satisfy the semi-classical regime requirement.

The beam at the output facet of the medium is imaged onto a camera, where we collect its interference with a reference beam. By using the Off-Axis Digital Holography technique46, we measure the spatial distribution of both the intensity and the phase of the output field. To explore the full dynamics of the laser beam, we tune the initial power from 0.2 W to 5.5 W. The insets in Fig. 1a show the experimental beam intensity profile at the input and output crystal facets for an input power P = 5 W.

Experimentally, it is only possible to access the field at the output facet of the sample and not the full nonlinear propagation inside the material.

However, there is a direct mapping between the power P and the propagation length z, when χ ≪ 1 and the beam initial phase is negligible (see Supplementary Notes 1 and 5). This mapping allows us to follow the z-evolution of the amplitude \({{{{{{{\mathcal{E}}}}}}}}\) by varying the input power of the beam and measuring the intensity at fixed z = L (L is the sample length). We then re-scale the propagation coordinate z in terms of a relevant dynamical characteristic scale \({z}_{{{{{{{{\rm{dyn}}}}}}}}}=s\sqrt{{n}_{{{{{{{{\rm{b}}}}}}}}}\kappa /(\alpha \beta P)}\). Therefore, varying the initial power P and measuring the intensity at fixed z = L is equivalent to measuring the intensity at different steps z/zdyn inside the material at fixed P. Hereafter, we will use P to parameterize the system evolution along z.

Observation of violent relaxation and formation of the quasi-stationary state

Figure 1b, d depict the numerical and experimental intensity profiles (along y = 0) measured at the output of the glass sample as a function of power P, respectively. We observe good qualitative agreement: the initial beam collapse is then followed by a stabilization in the sense that a central high intensity (or N-body particle density) peak is formed that, despite some oscillations, persists for the rest of the propagation. In nonlinear optics terms, the optical beam is undergoing self-focusing and is trying to stabilize on the central solitonic peak that acts as an attractor for the dynamics, by expelling energy in the form of a broader, lower amplitude surrounding field. The presence of nonlocality prevents the light from undergoing a catastrophic collapse in this system47,48,49. The semi-classical regime chosen is not ideal for the formation of a soliton, but instead maximizes the generation of the surrounding background that indeed is the result of phase-mixing and violent relaxation (in a gravitational context, this corresponds to the galaxy). A plot of the simulated intensity distribution for an incident power P = 5 W is shown in Fig. 1c, and is in good agreement with the experimental inset in Fig. 1a).

Violent relaxation can be identified by looking at the distribution of the chemical potential density \(\nu ({{{{{{{\mathcal{U}}}}}}}}/{{{{{{{{\mathcal{U}}}}}}}}}_{0})\) and at the phase-space behaviour along the evolution. We expect a variation in the chemical potential due to a variation in the overall potential \({{{{{{{\mathcal{V}}}}}}}}({{{{{{{{\bf{r}}}}}}}}}_{\perp },z)\)1. Figure 2 shows the experimental (a) and numerical (b) distribution of the normalized chemical potential density, \({{{{{{{\mathcal{U}}}}}}}}({{{{{{{\mathcal{E}}}}}}}})/{{{{{{{{\mathcal{U}}}}}}}}}_{0}\), obtained for various input powers, P, as well as the potential \({{{{{{{\mathcal{V}}}}}}}}/{{{{{{{{\mathcal{V}}}}}}}}}_{0}\) evolution (c), computed at z = L. The initial region (0–2 W) is where the variation in the potential is strongest and corresponds to the same region in which the chemical potential density variation and mixing is also strongest, as expected for violent relaxation. In contrast, after the collapse (after P = 3 W), the distribution of the chemical potential density exhibits two characteristic ‘structures’, which persist for the whole subsequent evolution: one at smaller values, which corresponds to the centre of the beam near the solitonic core; a second ‘structure’ at higher values related to the more external rings. It is worth noticing that, after the collapse, the distribution of the chemical potential density does not vary significantly, tending asymptotically to a constant profile associated to the quasi-stationary state, formed by the narrow solitonic core plus the broad background analogue galaxy. Therefore, the region in which violent relaxation takes place can be identified between P ~ 0.5 W and P ~ 3 W. After this, the evolution tends to a quasi-stationary state where the chemical potential distribution of the system is constant, despite the fact that the intensity profile keeps evolving (see Fig. 1a). A similar behaviour is observed in the astrophysical context: galaxies can present a slowly evolving chemical potential, despite showing a continuing evolution in time of the overall shape.

Fig. 2: Distribution of the chemical potential density \(\nu ({{{{{{{\mathcal{U}}}}}}}}/{{{{{{{{\mathcal{U}}}}}}}}}_{0})\).
figure 2

Experiment a and simulation b. The chemical potential is normalised to \({{{{{{{{\mathcal{U}}}}}}}}}_{0}=\left(\alpha \beta {k}_{0}P\right)/\left(2\pi \kappa \right)\). c Numerical y = 0 slice of the normalized potential \({{{{{{{\mathcal{V}}}}}}}}/{{{{{{{{\mathcal{V}}}}}}}}}_{0}\), computed at z = L, as a function of transverse coordinate x and power P (\({{{{{{{{\mathcal{V}}}}}}}}}_{0}={k}_{0}P\)). \(t{\prime}\) is the light propagation time through the sample.

Full size image

We study the existence of mixing in the system by analysing the evolution of the phase-space. Figure 3 shows the experimental Wigner distribution F(r⊥, k⊥, z) of the full complex-valued optical field45. At the input, the system has a Gaussian spatial distribution with a very narrow dispersion along the kx-axis. As P increases, the phase mixing starts by first twisting the phase-space (indicated by the white arrows) and then forming filaments characteristic of violent relaxation1. We have also verified that in a system where only mixing is present (without violent relaxation), such as in the Snyder-Mitchell model50, the evolution of the system is significantly different (see Supplementary Note 4).

Fig. 3: Wigner Distribution.
figure 3

Results of experiment for the y = 0, ky = 0 profiles of the Wigner distribution at different input powers. a P = 0.2W, b P = 1W, c P = 2W, d P = 3W, e P = 4W, f P = 5W. The white arrows in c refer to the sense of twist of the Wigner distribution as the power is increased.

Full size image

Conclusions

Violent relaxation was first formulated by Lynden-Bell in 19674 in the context of galaxy formation but subsequent studies highlighted its very general nature that extends also to other systems. However, violent relaxation dynamics have never been observed experimentally before. We have provided experimental evidence of violent relaxation in an optical system that exhibits the required long-range interactions. The observed dynamics can also be likened to those of a gravitational system through the optical-gravitational analogy of the underlying equations. Indeed, we can directly connect our optical experimental parameters to those of a particle-based dark matter galaxy, as shown in Fig. 1e, where we plot the galaxy distribution for a particle system with parameters equivalent to those of the experiment, to be compared with Fig. 1c. This additional connection between the experiments and actual N-body particle systems arises as a result of the fact that the Schrödinger-Poisson equations (that underly our optical system) converge to the Vlasov-Poisson equations (that describe N-body system in the limit of large N).

A possible merit in this experimental evidence of violent relaxation is that these results extend beyond the original astrophysical setting (where in any case we would not be able to carry out repeatable experiments). Violent relaxation and related phenomena therefore can occur more readily than expected in other long-range systems and where these might actually dominate in the temporal or length scales that are relevant to the experiment itself. These experiments, aside from demonstrating violent relaxation for the first time in any experimental setting, to the best of our knowledge, provide a generic and easily accessible platform in which experiments can be carried out and tunes so as to recreate for example analogues of galaxy evolution or analogues of plasma evolution, where violent relaxation experiments and the underlying dynamics are less accessible.

Methods

Experimental setup

A CW laser with wavelength λ = 532 nm is split into 2 beams: a reference and a target beam. The reference beam is expanded by a system of lenses and collected by a CMOS camera. See Supplementary Note 1 for a detailed image of the setup. The target beam is shaped to have waist s = 350 μm (waist calculated where the intensity falls of 1/e2—the value has been obtained by a Gaussian fit of the beam intensity at the sample input face—see Supplementary Note 2) and shines onto three aligned identical slabs of lead-doped glass (height D = 5 mm, width W = 40 mm and length L0 = 100 mm each, hence a total length L = 300 mm).

The glass is a self-focusing nonlinear optical medium with background refractive index nb = 1.8, thermal conductivity κ = 0.7 Wm−1 K−1, absorption coefficient α = 1m−1, thermo-optic coefficient \(\beta =\frac{\partial n}{\partial T}=2.2\cdot 1{0}^{-5}\) K−1 and transmission coefficient at the sample interface T = 0.92. The value of the coefficient β is found by a fit of the experimental beam evolution and results to be 1.6 times larger than the value provided by the manufacturer. The target beam input powers range from 0.2 W to 5.5 W, with a 0.25 W step. By means of the off-axis digital holography technique46, we reconstruct the amplitudes and phases of the target beam.