Turbulence-induced melting of a nonequilibrium vortex crystal in a forced thin fluid film

To develop an understanding of recent experiments on the turbulence-induced melting of a periodic array of vortices in a thin fluid film, we perform a direct numerical simulation of the two-dimensional Navier-Stokes equations forced such that, at low Reynolds numbers, the steady state of the film is a square lattice of vortices. We find that, as we increase the Reynolds number, this lattice undergoes a series of nonequilibrium phase transitions, first to a crystal with a different reciprocal lattice and then to a sequence of crystals that oscillate in time. Initially the temporal oscillations are periodic; this periodic behaviour becomes more and more complicated, with increasing Reynolds number, until the film enters a spatially disordered nonequilibrium statistical steady that is turbulent. We study this sequence of transitions by using fluid-dynamics measures, such as the Okubo-Weiss parameter that distinguishes between vortical and extensional regions in the flow, ideas from nonlinear dynamics, e.g., \Poincare maps, and theoretical methods that have been developed to study the melting of an equilibrium crystal or the freezing of a liquid and which lead to a natural set of order parameters for the crystalline phases and spatial autocorrelation functions that characterise short- and long-range order in the turbulent and crystalline phases, respectively.


Introduction
A crystal melts into a disordered liquid when the temperature is raised beyond the melting point T M ; and when the liquid is cooled below T M it freezes again.This melting or freezing transition, one of the most common equilibrium phase transitions, has been studied extensively; and the Ramakrishnan-Yussouff density-functional theory [1,2,3,4,5], which uses a variational free energy, has led to a good understanding of such freezing.It is natural to ask whether there are nonequilibrium analogues of this transition.One example is shear-induced melting of colloidal crystals [6].We investigate another example of a nonequilibrium transition in which dynamically generated turbulence plays the role of temperature and disorders a crystal consisting of an array of vortices, of alternating sign, imposed on a thin fluid film by an external force.
Recent experiments [7,8] have explored this transition from a low-Reynolds-number nonequilibrium vortex crystal, imposed on a thin fluid film by a force that is periodic in space, to a high-Reynolds-number, disordered, nonequilibrium liquid-type phase.This problem has been studied earlier by using linear-stability analysis [9,10] and direct numerical simulations (DNS) [11,12].The former is well-suited to the study of the first instability of the vortex crystal with increasing Reynolds number Re; the latter have (a) studied the route to chaos, by using techniques from nonlinear dynamics [11,12] to analyse the temporal evolution of this system, or (b) have mimicked [13] the experiments of Ref. [7] to study the creation and annihilation of vortices.
We revisit the problem of turbulence-induced melting of a vortex crystal by conducting a direct numerical simulation (DNS).Our study combines methods from turbulence, nonlinear dynamics, and statistical physics to elucidate the nature of nonequilibrium phases and transitions in this system.We use the two-dimensional (2D) Navier-Stokes (NS) equation with air-drag-induced Ekman friction to model the thin fluid films used in experiments [7,8,14]; this is a good model for flows in such films so long as the Mach number is small and the corrections arising from the finite thickness of the film and from the Marangoni effect can be neglected [15,16,17].We force this 2D NS equation in a manner that mimics the forcing used in experiments and which yields, at low Re, a stationary, periodic array of vortices of alternating signs; we refer to this array as the vortex crystal.
We then investigate the stability of this crystal as we increase the amplitude of the force and, therefore, the Reynolds number.Measures from nonlinear dynamics, including time-series, power spectra, and Poincaré-type maps, can be used fruitfully here to examine the temporal behaviour of this system as it undergoes a sequence of transitions.This part of our study complements the work of Refs.[11,12].
Furthermore, we elucidate the natures of the transitions from the spatially ordered crystal to the disordered, turbulent state by borrowing ideas from the density-functional theory of the freezing of a liquid into a crystal [1,3,4,5].Since the density field of a crystal is periodic, this theory uses the coefficients in the Fourier decomposition of the density as order parameters.Specifically, in a conventional crystal, ρ admits the Fourier decomposition where the sum is over the vectors G of the reciprocal lattice; in the density-functional theory of freezing [1,3] the Fourier coefficients ρ G are identified as the order parameters of the liquid-to-crystal transition since their mean values vanish in the liquid phase for all nonzero reciprocal lattice vectors G.Moreover, ρ G=0 is the mean density and, because most liquids are nearly incompressible, it undergoes a very small change at this transition.Spatial correlations can be characterized conveniently by the autocorrelation function g(r) = [ρ(x + r)ρ(x)] , where the angular brackets denote Gibbsian thermal averages and the overline denotes spatial averaging over x; the Fourier transform of g(r) is related [1,3] to the static structure factor S(k).In an isotropic liquid phase the arguments of g and S are, respectively, r =| r | and k =| k |.
We use the Okubo-Weiss field Λ ≡ det(A) as the analogue of the density ρ(r) in the density-functional theory of freezing; here A is the velocity-derivative matrix that has components A ij ≡ ∂ i u j , with u j the j th component of the velocity.For an inviscid, incompressible 2D fluid the local flow topology can be characterized via the Okubo-Weiss criterion [18,19]; this provides a useful measure of flow properties even if viscosity and Ekman friction are present [14,17,20]: Regions with Λ > 0 and Λ < 0 correspond, respectively, to centres and saddles [18,19].Thus, in the nonequilibrium vortex crystal, Λ(r) is a periodic function; so, like ρ(r) in a conventional crystal, it admits the Fourier decomposition where the sum is over the reciprocal-lattice vectors k; it is natural to think of Λk as the order parameters that characterise the vortex crystal.In terms of these order parameters we can define the analogue of the static structure factor S(k) for a conventional crystal; for the vortex crystal this is the two-dimensional spectrum here the angular brackets do not imply a Gibbsian thermal average, as in equilibrium melting, but denote an average over the nonequilibrium state of our system; as we show below, this state can be steady in time, or it can oscillate periodically or quasiperiodically, or it can be a chaotic state that is statistically steady.The autocorrelation function is related to E Λ (k) by a spatial Fourier transform.The turbulent phase is nearly isotropic so, to a good approximation, G depends only on r ≡| r | in this phase; and it characterises the short-range order in the system exactly as g(r) does in an isotropic liquid.
In addition to identifying a natural set of order parameters for the turbulenceinduced melting of a two-dimensional vortex crystal in a thin fluid film, our study yields several interesting results that we describe qualitatively below: In particular, we find that, as we increase the Reynolds number, a rich sequence of transitions leads from the steady ordered vortex crystal to the disordered, turbulent state.The third column in Table 1 shows the types of nonequilibrium phases we encounter: There is SX, the original, steady square crystal imposed by the force; this is followed by SXA, steady crystals that are distorted, via large-scale spatial undulations, relative to SX; these give way to distorted crystals that oscillate in time, either periodically (OPXA) or quasiperiodically (OQPXA); finally the system becomes disordered and displays spatiotemporal chaos and turbulence (SCT).OPXA and OQPXA are actually a collection of several (perhaps an infinity) of nonequilibrium crystalline phases that oscillate in time; given the resolution of our calculation we can identify only some of these, as we describe in detail below.The combination of techniques that we use helps us to elucidate the natures of these phases in far greater detail than has been attempted hitherto.We believe that, by using the ideas we develop here, it should be possible to extend recent experiments [7] on turbulence-induced melting to uncover the rich series of nonequilibrium phase transitions that we have summarised in Table 1.
The remaining part of this paper is organised as follows.In Sec. 2 we present the equations we use to model turbulence-induced melting in thin fluid films and the numerical methods we use.Section 3 contains a detailed description of our results.Section 4 is devoted to a discussion and to conclusions.

Model and Numerical Methods
We begin with the 2D Navier-Stokes (NS) equations that can be written, as follows, in the nondimensional form of Ref. [21]: Here u ≡ (−∂ y ψ, ∂ x ψ), ψ, and ω ≡ ∇×u are, respectively, the velocity, stream function, and vorticity at the position x and time t; we choose the uniform density ρ = 1; α is the non-dimensionalised Ekman-friction coefficient, ν is the kinematic viscosity, and F ω ≡ −n 3 [cos(nx) + cos(ny)]/Ω, is the non-dimensionalised force with injection wave vector n, Ω = nF amp /(ν 2 k 3 ), and α = nνα k/F amp , where F amp is the forcing amplitude, α is the Ekman friction, and lengths are non-dimensionalised via a factor k/n with k a wave vector or inverse length [21].We denote the x and y components of the velocity as u 1 ≡ u and u 2 ≡ v, respectively.The spatially periodic force F ω yields, at low Ω, a vortex crystal that is also referred to as a cellular flow.A linear-stability analysis of this flow indicates that it has a primary instability [10] at a critical Reynolds number Re c ≡ √ 2 which translates into a threshold value Ω s,n ≡ nRe c .This primary instability yields another vortex crystal, which is steady in time but whose unit cell is larger than that of the original vortex crystal [10,22].
We solve Eqs. ( 5) and (6) numerically by using a pseudo-spectral method with a 2/3 dealiasing cut-off and a second-order Runge-Kutta scheme for time marching [23,24] with a time step δt = 0.01.We use N 2 collocation points; in most of our studies we use N = 128; we have checked in representative cases that our results are unchanged if we use N = 256.Our main goal has been to obtain long time series for several variables (see below) to make sure that the temporal evolution of our system is obtained accurately; most of our runs are at least as long as 3 × 10 6 δt.We monitor the time-evolution of (a) the total kinetic energy E(t) ≡ u 2 , (b) the stream function ψ, (c) the vorticity ω, (d) the Okubo-Weiss parameter Λ, and (e) the k = (1, 0) component of the Fourier transform v of the y component v of u.Given these time series we obtain E Λ (k) at representative times and G(r), which is obtained by averaging over 20 configurations of Λ(r) separated from each other by 10 5 δt, after transients in the first 10 6 time steps have been removed.From the time series of E(t) we obtain its temporal Fourier transform E(f ) and thence the spectrum | E(f ) | that helps us to distinguish between periodic, quasiperiodic, and chaotic temporal behaviours.We also augment this charaterisation by using Poincaré-type sections in which we plot v(1,0) versus v(1,0) at successive times (see, e.g., Ref. [21] for the Kolmogorov flow).
We show below that the vortex crystal melts, as we increase Ω, via a rich sequence of transitions.The principal effect of the Ekman friction is to delay the onsets of these transitions; we have checked this explicitly in some cases.However, to make contact with earlier linear-stability and DNS studies of this problem, the results we present below have been obtained with no Ekman friction.Our qualitative conclusions are not affected by this.
So long as Ω < Ω s,n , the steady-state solution [21], indicated by the subscript s, of Eq. ( 5) is ω s,n = −n[cos(nx) + cos(ny)].We examine the destabilisation of this state, with increasing Ω, for two representative values of n, namely, n = 4 and n = 10 for which Ω s,4 5.657 and Ω s,10 14.142; in our DNS we choose the initial vorticity field to have the form 2 ); and then we let the system evolve under the application of the force F ω .We increase Ω from Ω s,n to 6.85Ω s,n in steps of 0.15Ω s,n , for n = 4 (runs R4-1 to R4-7 in Table 1), and from Ω s,n to 15.9Ω s,n in steps of 0.1Ω s,n , for n = 10 (runs R10-1 to R10-4 in Table 1).To trace the transition to chaos accurately from SXA to SCT we have also conducted runs where Ω is increased in steps of 0.08Ω s,n for n = 4 (runs R4-4 to R4-6 in Table 1) and 0.07Ω s,n for n = 10 (runs R10-2 and R10-3 in Table 1).We have benchmarked our numerical scheme by comparing results from our code with those of Ref. [21], which deals with a Kolmogorov flow imposed by an external force of the form F ω = n cos(ny).
For Ω < Ω s,n , the Λ field shows alternating centres and saddles, arranged in a two-dimensional square lattice, which we illustrate via pseudocolour plots for n = 4 and n = 10 in Figs. 1 (a plots of ψ in Figs. 2 (a) and (b), respectively.

Results
In this section we present the results of our numerical simulations for n = 4 and n = 10 and the ranges of parameters given in Table 1.In the first subsection we present our results for n = 4; the next subsection contains our results for n = 10; in the last subsection we describe our results for the order parameters and spatial correlation functions.

The case n = 4
When we increase Ω beyond Ω s,n=4 , the steady-state solution ω s,n=4 becomes unstable.
From the range of values of Ω in our runs R4 − 2 (Table 1) we observe that a new steady state is attained, which we illustrate, for Ω = 6.5, via pseudocolour plots of ψ and Λ in Figs.3(a) and (b), respectively.The new steady state is also a vortex crystal; however, it is different from the original vortex crystal as can be seen especially clearly by comparing the pseudocolour plots of ψ in Figs. 2 (a) and 3(a).This difference also shows up as a very slight distortion of the crystalline structure in the pseudocolour plot of Λ shown in Fig. 3(b).To use the language of solid state physics, this is an example of a very weak structural phase transition.Normally such a phase transition is mirrored in new superlattice peaks that appear in the reciprocal-space spectrum E Λ in addition  to the dominant peaks associated with the original crystal structure; however, given the weakness of the distortion, such superlattice peaks are not visible, with our resolution in Fig. 3(c).Clear examples of such superlattice peaks appear as we increase Ω (see below).
At Ω = 8.202, a new regime appears (run R4 − 3).The time-series of the energy E(t) now shows a periodic array of spikes.This regime has no analogue in a conventional crystal; it is a crystal that oscillates periodically in time and, to that extent, it can be thought of as a spatiotemporal crystal.The time between successive spikes is very large ( 10 4 δt) as shown by the plot of E(t) in Fig. 4(a); this is why our DNS runs must be very long to distinguish such states from one that is steady; we have also checked that the time between successive spikes is the same (to three-figure accuracy) for N = 64 and N = 128.In Figs.4(b) and (c) we show pseudocolour plots of ψ and Λ, respectively; the former shows a large-scale spatial undulation and the latter some deformation relative to the original crystal.This deformation is also mirrored in the distortion, relative to Fig. 3(c), of the dominant peaks in the reciprocal-space spectrum E Λ (k) shown in Fig. 4(d).
For runs R4 − 4, i.e., 9.05 ≤ Ω < 15.3, we find a new crystalline state that is steady in time.It has a large-scale spatial undulation relative to the original vortex crystal as illustrated, for Ω = 11.3, by the pseudocolour plots of ψ and Λ in Figs.5(a) and (b), respectively.This undulation leads to a distortion of the dominant peaks in the reciprocal-space spectrum E Λ (k) of Fig. 5(c), which also shows new superlattice peaks that occur at smaller values of (k x , k y ) relative to the dominant peaks.
On further increasing Ω we enter a new regime (runs R4 − 5, i.e., 15.3 ≤ Ω < 17.3) in which we have a spatiotemporal crystal, i.e., a spatially periodic Λ that oscillates in time.The time-series E(t) displays a periodic array of spikes as shown in Fig. 6(a); this leads to the frequency-space spectrum | E(f ) | of Fig. 6(b).The peaks in this spectrum can be labelled as f 0 , where is a positive integer and f 0 is the fundamental frequency that can be obtained from the inverse of the temporal separation between successive spikes in Fig. 6(a); this is a clear signature of periodic temporal evolution.The Poincarétype map in the ( [v(1, 0)], [v(1, 0)]) plane, Fig. 6(c), shows that the projection of the attractor on this plane is a closed loop in this case.Pseudocolour plots of ψ and Λ [ Figs. 7(a)-(b)] are similar, respectively, to those in Figs.5(a)-(b) if we look at their spatial patterns; however, they oscillate in time as can be seen most clearly from their animated versions [mpeg files psi movie R5.mpeg-lam movie R5.mpeg].The associated reciprocal-space spectrum E Λ also oscillates between the spectra shown in Figs.7(c) and (d) as can be seen clearly from its animated version [avi file lamf movie R5.avi].
The time between successive spikes in E(t) decreases as Ω increases.To quantify this, we define the inter-spike interval T i as follows: T i ≡ t i+1 − t i , where t i is the time at which E(t) crosses, for the i th time, its mean value, E(t) , from below; we can think of i as the spike index.In Fig. 8(a) we plot T i versus i; this shows that the mean value of T i decreases as Ω increases [Fig.8(b)]; furthermore, T i oscillates slightly about its mean value for any given value of Ω.The magnitude of these oscillations, which we have  used for the error bars in Fig. 8(b), decreases as Ω increases.We have checked explicitly that our results here do not change when we increase the resolution of our DNS from N = 128 to N = 256.
At Ω = 17.8, i.e., run R4 − 6, another transition occurs: The time series of the energy E(t) and its frequency spectrum | E(f ) | are shown, respectively, in Figs.9(a) and (b).The latter displays peaks superimposed on a noisy background signal; these peaks can be indexed as f 0 , f 1 , f 1 − 2f 0 , f 0 − 2f 1 , and 3f 0 − 2f 1 , within our numerical accuracy and with f 0 0.001653 and f 1 0.001707.Since f 0 /f 1 is not a simple rational number, we conclude that Fig. 9(b) indicates principally quasiperiodic temporal evolution with a small chaotic admixture, the former associated with the peaks indexed above and the  latter with the noisy background signal.We believe the chaotic part of the signal comes from transitions between the elliptical islands in the Poincaré-type section of Fig. 9(c).The plot of the inter-spike interval T i versus the spike index i in Fig. 9(d) confirms the complicated temporal evolution of this state.
As we increase Ω further (R4 − 7) the temporal evolution of the system becomes ever more chaotic; this is associated with a disordered pattern of vortices in space too.Thus we obtain a state with spatiotemporal chaos and turbulence, which is our analogue of the liquid state.We illustrate this for Ω = 50.We begin with the time series of E(t) and the spectrum | E(f ) | in Figs.10(a) and (b), respectively; the latter shows clearly a broad background that is indicative of temporal chaos.This is further confirmed by the nearly uniform spread of points in the Poincaré-type section [Fig.10(c)] in the ( [v(1, 0)], [v(1, 0)]) plane.The disordered spatial organisation of this state is illustrated by the pseudocolour plots of ψ and Λ [Figs.11(a)-(b), respectively] and the reciprocal-space spectrum E Λ (k) of Fig. 11(c) that shows several new modes in addition to the original peaks, whose vestiges are still visible since F ω continues to act on the system.
Thus we see that the turbulence-induced melting of our nonequilibrium vortex crystal is far richer than its equilibrium counterpart.For the case n = 4 investigated above it proceeds as described in column 4 in Table 1.Before we present an analysis of the disordered state in terms of the spatial autocorrelation function G(r) and the evolution of the order parameters Λk with Ω, we give below a short summary of our results for n = 10; the route to turbulence is slightly different for this case than for n = 4.

The case n=10
Our results for n = 10 are based on the runs R10 − 1 to R10 − 5 in Table 1.
For Ω < Ω s,n the steady vortex crystal is shown by the pseudocolour plot of Λ in Fig. 1(b).As we increase Ω beyond Ω s,n we find in runs R10 − 2, i.e., the range Ω s,n < Ω < 22.6, a new steady state in which pseudocolour plots of ψ and Λ show largescale spatial undulations caused by small deformations of the original vortex crystal [Figs.12 (a) and (b), respectively]; consequently the dominant peaks in the reciprocalspace spectrum E Λ are slightly distorted.
Around Ω = 24 (run R10 − 3) another transition occurs: the time series of E(t) is periodic [Fig.13(a)] and its spectrum | E(f ) | [Fig.13(b)] shows one dominant peak, i.e., higher harmonics are nearly absent.Thus the Poincaré-type plot in the ( v [1,0], v [1,0]) plane [Fig.13(c)] displays a simple attractor.The spatial structure of this state is illustrated by the pseudocolour plots of ψ and Λ shown, respectively, in Figs.14(a) and (b); the associated reciprocal-space spectrum E Λ is shown in Fig. 14(c).Given the temporal behaviour of this state, these structures, in real or reciprocal space, oscillate in time at the frequency given by the temporal evolution of E(t).For Ω = 24 the large-scale spatial structures in ψ oscillate around their mean positions; but, for 25 ≤ Ω ≤ 28, we find a travelling-wave type pattern, which reenters our simulation domain by virtue of the periodic boundary conditions that we use in our pseudo-spectral method.If we compare the frequency spectra | E(f ) | for the cases Ω = 24 and Ω = 28, we find higher harmonics in the latter but they are all multiples of one fundamental frequency.
In Fig. 15(a) we show plots of T i versus i (cf.Fig. 8(a) for n = 4) for various values of Ω. From these we obtain the plot of the mean value T i versus Ω shown in Fig. 15.This first decreases, as we increase Ω, and then increases mildly at Ω = 28.
For Ω ≥ 29 the time-series of E(t) appears chaotic and the associated frequency spectrum | E(f ) | displays a broad background, as we show in the illustrative Figs.16 (a) and (b) for Ω = 225.The associated Poincaré-type section in Fig. 16(c) confirms that the temporal behaviour is chaotic.The spatial patterns are also disordered as we show via the pseudocolour plots of ψ and Λ in Figs.17(a) and (b), respectively.The corresponding reciprocal-space spectrum E Λ shows that a large number of modes are excited.Thus we have both spatial disorder and temporal chaos; as in the case n = 4, the analogue of the liquid state is a turbulent one with spatiotemporal chaos.However, given the resolution in Ω that we have been able to obtain in our calculations, the route to this state of spatiotemporal chaos is different for n = 10 and n = 4 as can be seen from column 4 in Table 1.

Order parameters and spatial autocorrelation functions
We now return to ideas borrowed from the density-functional theory [1,3,4,5] of freezing by examining the behaviour of the order parameters Λk as functions of Ω. Recall that, in equilibrium melting, ρ G jumps discontinuously from a nonzero value in the crystal to zero in the liquid at the first-order melting transition.As we have noted above, the turbulence-induced melting of our vortex crystal is far more complicated; it proceeds via a sequence of transitions.In turbulence-induced melting of a vortex crystal Λk , obtained by summing Λk over the four forcing wave vectors, is the equivalent of For small values of Ω the steady state is SX so, in Fourier space, only modes with the forcing wave vector are significant.On increasing Ω, the crystal undergoes a series of transitions that take it from the crystal SX to the disordered, turbulent state SCT; as this happens the spectral weight slowly shifts away the forcing wave numbers.This explains the trend observed in Figs.18.The spatial correlations in the crystalline and disorderd states can be characterised by the spatial autocorrelation function G(r) defined in Eq. ( 4).Representative plots are shown, respectively, for n = 4 and n = 10 in Figs.19 (a)-(d).For the crystalline case we evaluate G(r) along the line connecting r = (π/2, π/2) and r = (π/2, π); this shows a periodic array of peaks [Figs.19 (a) and (b) for n = 4 and n = 10, respectively]; the widths of these peaks are related to the widths of vortical or strain-dominated regions.
In the turbulent phase we present data obtained by a circular average of G [Figs.19  (c) and (d) for n = 4 and n = 10, respectively]; here the peaks decay over a length scale that indicates the degree of short-range order.This decay is similar to the decay of spatial correlation functions in a disordered liquid.

Conclusions
We have carried out a detailed numerical study of turbulence-induced melting of a nonequilibrium vortex crystal in a forced, thin fluid film.We use ideas from the density-functional theory of freezing [1,3,4,5], nonlinear dynamics, and turbulence to characterise this.Correlation functions, similar to those in liquid-state theory, have been used by some recent experiments [7,8] to analyse the short-range order in the turbulent phase; nonlinear-dynamics methods, such as Poincaré-type maps, have been used in the numerical studies of Ref. [11]; experimental studies have used the curvature of Lagrangian trajectories to identify extrema in vortical and strain-dominated regimes.To the best of our knowledge, there is no study that brings together a variety of methods, as we do, to analyse turbulence-induced melting.
The advantages of our approach are as follows: (a) it helps us to identify the order parameters for turbulence-induced melting and thus contrast it with conventional melting; (b) the sequence of transitions can be characterised completely in terms of the Eulerian fields ψ and Λ, the total energy E(t), and suitable Fourier transforms of these; (c) spatial correlations in crystalline and turbulent phases can be studied conveniently in terms of G.
Equilibrium phase transitions occur strictly only in the thermodynamic limit that is, roughly speaking, the limit of infinite size [25].It is interesting to ask how we might take the thermodynamic limit for the vortex crystals we have studied here.There seem to be at least two ways to do this: (a) in the first the system size can be taken to infinity in such a way that the areal density of the vortical and strain-dominated regimes remains the same in the ordered, crystalline phase; (b) in the second way we can increase the parameter n in the forcing F ω so that more and more unit cells occupy the simulation domain (see, e.g., Figs.1(a) and (b) for n = 4 and n = 10, respectively).Such issues have not been addressed in detail by any study, partly because, for large system sizes, it is not possible to obtain the long time series that are required to characterise the temporal evolution of the system (especially in the states we have referred to as spatiotemporal crystals).In particular, it is quite challenging to investigate the system-size dependence of the transitions summarised for n = 4 and n = 10 in Table 1.From Figs. 19 (c) and (d) we can extract a correlation length; this length is much smaller than the linear size of our simulation domain so we expect that our results in the turbulent phase (SCT) will not change if we increase the size of our system.However, subtle size dependence might occur in the ordered phase as follows: as we increase Ω, the original crystal is distorted by large-scale spatial undulations that are associated with the inverse cascade of energy in two-dimensional turbulence; if these undulations lead to crystalline phases with larger and larger unit cells, then our finite-size calculations will become unreliable when the size of the unit cell becomes comparable to the size of our simulation domain.A systematic study of such subtle finite-size effects is a very challenging task that requires much more detailed simulations than have been attempted so far.
As we have shown above, the sequence of transitions that comprise turbulenceinduced melting of a vortex crystal is far richer than conventional equilibrium melting.There is another important way in which the former differs from the latter: To maintain the steady states, statistical or strictly steady, we always have a force F ω ; thus, in the language of phase transitions, we always have a symmetry-breaking field, both in the ordered and disordered phases.Strictly speaking, therefore, there is no symmetry difference between the disordered, turbulent state and the ordered vortex crystal, as can be seen directly from the remnants of the dominant peaks in the reciprocal-space spectra E Λ (k) in Figs.11 and 17(c) for n = 4 and n = 10, respectively.One consequence of this is that the order parameters Λk , with k equal to the forcing wave vectors, do not vanish identically in the disordered, turbulent phase; however, they do assume very small values.Moreover, in the case of turbulence-induced melting the crystal undergoes a transition from an ordered state to an undulating crystal to a fully turbulent state.Thus there is no noise and hence no fluctuations in the crystalline state; i.e., it is equivalent to a crystal at zero temperature.This has no analogue in the equilibrium melting of a crystal.
In equilibrium, different ensembles are equivalent; we can, e.g., use either the canonical or the grand-canonical ensemble to study the statistical mechanics of a system and, in particular, the phase transitions in it.However, this equivalence cannot be taken for granted when we consider nonequilibrium statistical steady states (see, e.g., Ref. [26]).We have seen an example of this in Ref. [17] where certain PDFs show slightly different behaviours depending on whether we keep the Grashof number (i.e., the nondimensionalised amplitude of the force) fixed or whether we keep the Reynolds number fixed.Turbulence-induced melting offers another example of the inequivalence of dynamical ensembles: the precise sequence of transitions that we encounter in going from the vortex crystal to the turbulent state depends on whether we do so by changing the Grashof number as in Ref. [11] or whether we do so by changing Ω as we have done here.We have checked explicitly that we can reproduce the sequence of transitions in Ref. [11] if we tune the Grashof number rather than Ω to obtain the turbulent state.
Investigations of similar transitions, such as in the Kolmogorov flow [21,27], can benefit by using the combination of methods we have used above.Detailed studies of the effects of confinement, air-drag induced Ekman friction on turbulence-induced melting, initiated, e.g., in Refs.[12,22], can also make use of our methods, but that lies beyond the scope of this paper.We hope, therefore, that our study will encourage experimental groups to analyse turbulence-induced melting by using the set of techniques and ideas that we have described above.

Figure 3 .
Figure 3. (Colour online) Pseudocolour plots for Ω = 6.5 of (a) the streamfunction ψ and (b) the Okubo-Weiss parameter Λ with superimposed contour lines.(c) A filled contour plot of the reciprocal-space energy spectrum E Λ showing clear, dominant peaks at the forcing wave vectors.

Figure 4 .
Figure 4. (Colour online) Plot for Ω = 8.202 of (a) the time evolution of the energy E(t) for Ω = 8.202.Pseudocolour plots of (b) the streamfunction ψ and (c) the Okubo-Weiss parameter Λ with superimposed contour lines.(d) A filled contour plot of the reciprocal-space energy spectrum E Λ showing the distortions of peaks at the forcing wave vectors.

Figure 5 .
Figure 5. (Colour online) Pseudocolour plots for Ω = 11.3 of (a) the streamfunction ψ and (b) the Okubo-Weiss parameter Λ with superimposed contour lines.(c) A filled contour plot of the reciprocal-space energy spectrum E Λ showing clear, dominant peaks, at the forcing wave vectors, and subdominant superlattice peaks at smaller wave vectors.

Figure 7 .
Figure 7. (Colour online) Pseudocolour plots for Ω = 15.3 of (a) the streamfunction ψ and (b) the Okubo-Weiss parameter Λ with superimposed contour lines.(c) and (d) Filled contour plots of E Λ showing the two spectra between which it oscillates in time.

Figure 11 .
Figure 11.(Colour online) Pseudocolour plots for Ω = 50 of (a) the streamfunction ψ and (b) the Okubo-Weiss parameter Λ with superimposed contour lines.(c) A filled contour plot of the reciprocal-space energy spectrum E Λ which shows that a large number of modes are excited.

Figure 12 .
Figure 12. (Colour online) Pseudocolour plots for Ω = 22.62 of (a) the streamfunction ψ and (b) the Okubo-Weiss parameter Λ with superimposed contour lines.(c) A filled contour plot of the reciprocal-space energy spectrum E Λ showing clear, but slightly distorted, dominant peaks at the forcing wave vectors.

Figure 14 .
Figure 14.(Colour online) Pseudocolour plots for Ω = 24 of (a) the streamfunction ψ and (b) the Okubo-Weiss parameter Λ with superimposed contour lines.(c) A filled contour plot of the reciprocal-space energy spectrum E Λ showing clear, dominant peaks at the forcing wave vectors.

Figure 17 .
Figure 17.(Colour online) Pseudocolour plots for Ω = 225 of (a) the streamfunction ψ and (b) the Okubo-Weiss parameter Λ with superimposed contour lines.(c) A filled contour plot of the reciprocal-space energy spectrum E Λ showing that a large number of modes are excited.

Table 1 .
Table indicating the number of the Run (e. ) and (b), respectively.[These patterns are reminiscent of a two-dimensional version of a perfectly ordered binary alloy, with two kinds of atoms, whose analogues here are centres and saddles.]We show corresponding pseudocolour g., R4-1), the values of n and Ω, and the type of order.Here SX denotes the original, steady square crystal imposed by the force (the precise pattern depends on n); SXA denotes steady crystals that are distorted, via a large-scale undulation, relative to SX; OPXA indicates a crystal that is distorted slightly with respect to SX and which oscillates periodically in time; OQPXA is like OPXA but with quasiperiodic oscillations; SCT denotes the disordered phase with spatiotemporal chaos and turbulence.