# Melting of a nonequilibrium vortex crystal in a fluid film with polymers: elastic versus fluid turbulence

###### Abstract

We perform a direct numerical simulation (DNS) of the forced, incompressible two-dimensional Navier-Stokes equation coupled with the FENE-P equations for the polymer-conformation tensor. The forcing is such that, without polymers and at low Reynolds numbers Re, the film attains a steady state that is a square lattice of vortices and anti-vortices. We find that, as we increase the Weissenberg number Wi, a sequence of nonequilibrium phase transitions transforms this lattice, first to spatially distorted, but temporally steady, crystals and then to a sequence of crystals that oscillate in time, periodically, at low Wi, and quasiperiodically, for slightly larger Wi. Finally, the system becomes disordered and displays spatiotemporal chaos and elastic turbulence. We then obtain the nonequilibrium phase diagram for this system, in the plane, where , and show that (a) the boundary between the crystalline and turbulent phases has a complicated, fractal-type character and (b) the Okubo-Weiss parameter provides us with a natural measure for characterizing the phases and transitions in this diagram.

###### pacs:

47.50.Cd,47.27.ek,47.57.Ng,05.70.Fhalso at ]Jawaharlal Nehru Centre For Advanced Scientific Research, Jakkur, Bangalore, India.

## I Introduction

The equilibrium melting transition, from a spatially periodic crystal to a homogeneous liquid, has been studied extensively ramakrishnan1979first; haymet1987theory; chaikin2000principles; oxtoby1991liquids; singh1991density; perlekar2010turbulence. Nonequilibrium analogs of this transition have been explored in, e.g., the melting of colloidal crystals by shear das2002melting and the turbulence-induced melting of a periodic array of vortices and anti-vortices in a forced, two-dimensional (2D), fluid film perlekar2010turbulence. We elucidate, via extensive direct numerical simulations (DNSs), the nonequilibrium melting of a periodic array of vortices and anti-vortices in a forced, 2D, fluid film with polymer additives; we refer to this periodic array as a nonequilibrium, 2D crystal; this is also called a cellular flow gotoh1984instability; gotoh1987instability in some studies.

We show that such a crystal can be melted either (a) by increasing the Reynolds numbers Re^{1}^{1}1
Strictly speaking, we should use the Grashof number , the nondimensionalized force doering1995applied, because we control the forcing, which leads eventually to the Reynolds number of the statistically steady of our fluid system. However, it has been conventional to use in studies of vortex crystals or cellular flows perlekar2010turbulence; gotoh1984instability; therefore, we use too. or (b) by increasing
the Weissenberg number Wi. Case (a) can be thought of as a
nonequilibrium melting transition, which is driven by fluid
turbulence perlekar2010turbulence; ouellette2007curvature; ouellette2008dynamic
and in which the disordered phase is a turbulent fluid that shows dissipation
reduction, because of the polymer
additives perlekar2010direct; perlekar2006manifestations; gupta2015two.
Case (b) is also a nonequilibrium melting transition; in this case the vortex
crystal melts into a disordered state, which is a polymeric fluid that shows
elastic turbulence groisman2000elastic; groisman2001efficient or
rheochaos ganapathy2006intermittency; majumdar2011universality. We show
that this system of driven vortices and anti-vortices provides an excellent
laboratory for studying the crossover from the dissipation-reduction to the
elastic-turbulence regimes in a fluid film.

Case (a), namely, turbulence-induced melting, has been explored in experiments ouellette2007curvature; ouellette2008dynamic, which use a spatially periodic force to impose a nonequilibrium vortex crystal on a fluid film at low Re; this crystal melts, at high Re, into a disordered, nonequilibrium, liquid-type phase. This melting problem has been studied by using linear-stability analysis gotoh1984instability; gotoh1987instability, for the first instability of the vortex crystal, and direct numerical simulations (DNS) perlekar2010turbulence; braun1997bifurcations; molenaar2007attractor. These DNSs have explored the nonequilibrium phase diagram of this driven, dissipative system by using techniques from nonlinear dynamics perlekar2010turbulence; braun1997bifurcations; molenaar2007attractor and generalizations of order parameters that are used to characterize equilibrium crystals ramakrishnan1979first; haymet1987theory; chaikin2000principles; perlekar2010turbulence.

To the best of our knowledge, case (b), namely, elastic-turbulence-induced melting of a vortex crystal, has not been explored in experiments, even though there have been several experimental studies of elastic turbulence groisman2000elastic; groisman2001efficient and the related phenomenon of rheochaos ganapathy2006intermittency; majumdar2011universality. Direct numerical simulations of 2D fluid films with polymer additives have begun to explore the berti2008two; berti2010elastic elastic-turbulence-induced melting of spatially periodic flows; the Kolmogorov-flow pattern investigated in Ref. berti2008two; berti2010elastic is a one-dimensional crystal.

The nonequilibrium vortex crystal, whose melting we study, is a 2D crystal. To describe the ordering in this crystal and its subsequent melting, either by elastic or fluid turbulence, we follow the treatment of Ref. perlekar2010turbulence, which has investigated this problem in the absence of polymers. In particular, we solve the incompressible, 2D, Navier-Stokes (NS) equation, for the fluid film, coupled to the FENE-P equation, for the polymer-conformation tensor ; this provides a good description of our system if we can neglect the film thickness and the Marangoni effect, and if the Mach number is low perlekar2010turbulence; chomaz1990soap; chomaz2001dynamics; perlekar2009statistically. The force we use yields, at low Re, a spatially periodic array of vortices and antivortices that is stationary in time; we call this a vortex crystal.

As we increase either (a) the forcing amplitude and, thereby Re or (b) the polymer-relaxation time and, hence, Wi, this crystal loses stability and gives way to a variety of nonequilibrium states. We show that time series of the energy of the fluid, their power spectra, and Poincarémaps help us to understand the spatiotemporal evolution of this system and its sequence of nonequilibrium transitions (cf., Refs. perlekar2010turbulence; braun1997bifurcations; molenaar2007attractor for turbulence-induced melting in a fluid film without polymers).

Our exploration of these transitions, which is inspired by the density-functional theory of the crystallization of a liquid ramakrishnan1979first; chaikin2000principles; oxtoby1991liquids; singh1991density; perlekar2010turbulence, is based on order parameters for our vortex crystal. The crystal density is periodic in space; the density-functional theory uses its Fourier coefficients as order parameters.

We follow Ref. perlekar2010turbulence and use the Okubo-Weiss field in the vortex-crystal instead of the crystal density , in the density-functional theory mentioned above; has elements , with the component of the velocity. The Okubo-Weiss criterion okubo1970horizontal; weiss1991dynamics, first developed to characterize the local topology of inviscid, incompressible, 2D, fluid flow, is also useful if viscosity and air-drag-induced friction are present perlekar2009statistically; rivera2001universal; elhmaidi1993elementary; this criterion states that , where the flow is vortical, and , where it is extensional or strain-dominated okubo1970horizontal; weiss1991dynamics. In a nonequilibrium vortex crystal, has the Fourier expansion

(1) |

where are the reciprocal-lattice vectors and the are the natural order parameters for the vortex crystal. In our 2D vortex crystal, the counterpart of the structure factor , for a crystal ramakrishnan1979first; chaikin2000principles; oxtoby1991liquids; singh1991density; perlekar2010turbulence, is the 2D spectrum

(2) |

where the angular brackets denote the average over the nonequilibrium state of our system (and not a Gibbsian thermal average); this state can (a) be independent of time, (b) exhibit periodic or quasiperiodic oscillations, or (c) be spatiotemporally chaotic and turbulent, but statistically steady. We define

(3) |

an autocorrelation function whose spatial Fourier transform is related to . This characterizes the spatial correlations in our system (the overbar in Eq. 3 indicates the average over the origin ). As we show below, the near isotropy of the turbulent phase implies that, to a very good approximation, depends only on here; and it characterizes the short-range vortical order in our 2D fluid film, exactly as the density correlation function does in an isotropic liquid at equilibrium ramakrishnan1979first; chaikin2000principles; oxtoby1991liquids; singh1991density; perlekar2010turbulence. We also follow the spatiotemporal evolution of the polymer-conformation tensor and relate it to the spatiotemporal evolution of .

We use the order parameters, defined above, to obtain many interesting results for the elastic-turbulence-induced melting of such a 2D vortex crystal in a fluid film with polymer additives. We first describe our principal results qualitatively. We show that, as Re or Wi increase, the steady, ordered, vortex crystal proceeds to the disordered, turbulent state via a sequence of nonequilibrium transitions between different nonequilibrium phases, which we list in Table 1 (column four): SX is the original, temporally steady, crystal, with a square unit cell that is imposed by the force; this is followed by phases of type SXA, which are temporally steady crystals that are distorted relative to SX; we then find distorted crystals that displays periodic (OPXA) or quasiperiodic (OQPXA) oscillations in time; finally, the system exhibits spatiotemporal chaos and turbulence (SCT). We conjecture that OPXA and OQPXA are actually several (perhaps an infinity) different nonequilibrium spatiotemporal crystals, which are periodic in space and time; given the spatial and temporal resolutions of our calculation, we can identify only some of these. In addition to the order parameters and correlation functions mentioned above, we also examine the spatiotemporal evolution of the polymer-conformation tensor in these phases. Our work leads to a rich, nonequilibrium phase diagram for our system. This has not been studied hitherto.

The remainder of this paper comprises the following Sections. Section 2, which is devoted to the equations that we use to model 2D fluid films with polymer additives and the numerical methods we employ. Section 3 contains our results. Section 4 is devoted to a discussion of our results and to conclusions.

## Ii Model and Numerical Methods

To study a 2D fluid film with polymer additives, we write the 2D, incompressible, NS and FENE-P equations in terms of the streamfunction and the vorticity , where is the fluid velocity at the point and time , as follows (we use the nondimensional form suggested in Ref. platt1991investigation for 2D fluid films without polymer additives):

(4) | |||||

(5) | |||||

(6) |

here, , the uniform solvent density , is the non-dimensionalized, air-drag-induced friction coefficient, , the Weissenberg number , is the kinematic viscosity, the zero-shear viscosity parameter for the solute (FENE-P), the polymer relaxation time, and is the non-dimensionalized, spatially periodic force, whose injection wave vector is related to , , and , where is the forcing amplitude, and is the coefficient of friction, Re is the Reynolds number. We non-dimensionalize lengths by a factor , with a wave number or inverse length platt1991investigation. The and components of the velocity are and , respectively. The superscript denotes a transpose, are the elements of the polymer-conformation tensor (angular brackets indicate the average over polymer configurations), the identity tensor has the elements , is the FENE-P potential that ensures finite extensibility, and and are, respectively, the length and the maximal possible extension of the polymers.

Our direct numerical simulation (DNS) of Eqs. (4)-(6) uses the MPI code that we have developed to study fluid turbulence with polymer additives in fluid films gupta2015two. We use periodic boundary conditions, a square simulation domain of side , and collocation points; these boundary conditions are well suited for studying vortex crystals with square unit cells. We use a fourth-order, Runge-Kutta scheme, with time step , for time marching and an explicit, fourth-order, central-finite-difference scheme in space and the Kurganov-Tadmor (KT) shock-capturing scheme kurganov2000new for the advection term in Eq. (6), because this scheme (see Eq. (7) of Ref. perlekar2010direct) successfully resolves steep gradients in and, thereby, minimizes dispersion errors. We solve the Poisson equation (5) in Fourier space by using a pseudospectral method and the FFTW library frigo1999fftw. We choose to be small enough to prevent from becoming larger than (if , there is a numerical instability). To preserve the symmetric-positive-definite (SPD) nature of at all times we adapt to 2D the Cholesky-decomposition scheme of Refs. perlekar2010direct; perlekar2006manifestations; gupta2015two; vaithianathan2003numerical. In particular, we define , so Eq. (6) becomes

(7) |

where , , , and Note that , and hence , are SPD matrices, so , where is a lower-triangular matrix with elements , such that for . We can now use Eq.(7) to obtain, for and

(8) | |||||

Equation(8) preserves the SPD property of if ; we enforce this, as in Refs. perlekar2010direct; perlekar2006manifestations; gupta2015two, by considering the evolution of instead of .

In most of our studies, we use a time step to for ; its precise value depends on the value of the polymer-relaxation time ; for , and otherwise (). Our results do not change significantly if we use or (we have checked this in representative cases). We obtain long time series for several variables (see below) to make sure that the temporal evolution of our system is captured accurately; most of our runs are at least as long as . We monitor the time-evolution of the total kinetic energy and calculate the stream function , vorticity , polymer-conformation tensor , the Okubo-Weiss parameter , and the kinetic energy spectrum , where indicates the time average over the statistically steady state. From these we obtain and , by averaging over configurations of , which are separated from each other by ; to eliminate the effects of initial transients, we remove data from the first time steps before we collect data for these averages. From we obtain the velocity and thence the component of the Fourier transform of the component of . We use the time series of to obtain its temporal Fourier transform and from that the power spectrum , which we use to differentiate between periodic, quasiperiodic, and chaotic temporal evolutions. We also use Poincaré-type sections, namely, plots of versus at successive times (for the Kolmogorov flow see Ref. platt1991investigation).

The equation of motion for a small, rigid particle (henceforth, an inertial particle) in an incompressible flow maxey1983equation; gatignol1983faxen assumes the following simple form bec2006acceleration; boffetta2004large, if the particle is much heavier than the advecting fluid:

(9) |

here, and are, respectively, the velocity and position of an inertial particle, is the Stokes or response time of the particle, whose ratio with the Kolmogorov dissipation time is the Stokes number . We assume that the radius of a particle , where is the dissipation scale of the carrier fluid, and that the particle density is so low that we can neglect interactions between particles. We also assume, as in Ref. de2012control, that the polymer additives affect the inertial-particle velocity only by virtue of their effect on the Eulerian velocity of the fluid at the position of the inertial particle. We show below that the spatial distribution of such inertial particles in the flow can also be used to distinguish the nonequilibrium crystalline and turbulent states in our system.

## Iii Results

We now present the results of our DNS. We begin with an overview, in the first subsection. In the next subsection, we present the nonequilibrium phase diagram that we obtain from our calculations. We then present the results of our numerical simulations for and the ranges of parameters given in Table 1. In the third and fourth subsections, respectively, we present our results for and . In the fifth subsection, we report our results for the order parameters, spatial correlation functions, and kinetic-energy spectra that we have defined above; and the last subsection contains our results for the spatial patterns of polymers and inertial particles.

### iii.1 Overview

We demonstrate below how our crystal of vortices melts, as we increase Wi for a fixed, low value of , through a sequence of transitions. At somewhat higher values of , at which the vortex crystal has already reached a melted state because of fluid turbulence, we find that this turbulent state moves first to a frozen state, as we increase Wi; this frozen state melts again as Wi increases some more. The air-drag-induced friction delays these transitions; we have checked this in some representative cases. However, to make contact with earlier studies of turbulence-induced melting of a vortex crystal perlekar2010turbulence; braun1997bifurcations; molenaar2007attractor without polymers, we do not include air-drag-induced friction here; our conclusions are not altered qualitatively by this.

The steady-state solution platt1991investigation; perlekar2010turbulence, of Eq. (4) is , if and ; here the subscript stands for steady state. We examine the destabilization of this solution, by increasing Wi, while we hold fixed, for one representative value of , namely, for which (see Ref. perlekar2010turbulence). The initial vorticity is taken to be ; and the system then evolves under the force . For a fixed value of , we increase Wi in steps of ; and we carry out simulations for , and from then to in steps of . To track the nonequilibrium transitions in our system, say from SXA to SCT, we have also carried out runs in which we increase Wi in steps of . We have validated our numerical scheme by comparing our results with those of Refs. platt1991investigation; perlekar2010turbulence. Reference platt1991investigation studies a Kolmogorov flow ; and Ref. perlekar2010turbulence uses an external force of the form , which is exactly the same as the one we employ.

Wi | |||
---|---|---|---|

For and , the field shows alternating vortical and extensional regions, which are also referred to as centers and saddles; these are arranged in a two-dimensional square lattice, which we illustrate via pseudocolor plots of , for , in Fig. 1 (a). For , which is well beyond , the fluid becomes turbulent; but, on the addition of polymers, Wi increases, and the fluid again forms a 2D vortex lattice, which we illustrate via pseudocolor plots of , for , in Fig. 1 (b). In Figs. 2 (a) and (b) we show the corresponding pseudocolor plots of the streamfunction .

### iii.2 Nonequilibrium Phase Diagram

We have carried out a very large number () of DNS studies of our system for a wide range of parameter values in the - plane; only a very small fraction of these parameter values are listed in Table 1. For each one of these parameter sets we calculate the quantities, such as , which we have defined above; and from these calculations we obtain the rich, nonequilibrium phase diagram of Fig. 3 (a), in which blue circles represent the steady crystal SXA, green circles represent either temporally periodic (OPXA) or quasiperiodic (OQPXA) states, and the red circles represent the state SCT that displays spatiotemporal chaos and turbulence. The most striking feature of this nonequilibrium phase diagram is that the boundaries between the different states are very complicated, with a fractal-type interleaving of the states SXA, OPXA/OQPXA, and SCT; this is especially clear in Fig. 3(b), which shows a detailed view of the phase diagram in the vicinity of the these boundaries. Earlier studies of elastic turbulence seem to have missed the complicated nature of these boundaries because they have not been able to examine these transition in as much detail as we have done.

Fractal-type boundaries between different nonequilibrium states have been obtained in nonequilibrium phase diagrams for some other extended dynamical systems. Examples of such boundaries can be found in the transition to turbulence in pipe flow schneider2007turbulence, the transitions between different spiral-wave patterns in mathematical models for cardiac tissue shajahan2007spiral; shajahan2009spiral, and the transition from no-dynamo to dynamo regimes in a shell model for magnetohydrodynamic (MHD) turbulence sahoo2010dynamo. In many dynamical systems, a fractal basin boundary can separate the basin of attraction of a fixed point or a limit cycle from the domain of attraction of a strange attractor; thus, a small change in the initial condition can lead either to the simple dynamical evolutions, associated with fixed points or limit cycles, or to chaos, because of a strange attractor. The 2D Navier-Stokes (NS) and FENE-P equations that we study here comprise a -dimensional dynamical systems, where is the number of grid points; the factor of appears here because we have the 2D vorticity and the three independent components of the matrix at each grid point. It is not feasible to find the complete basin boundaries for such a high-dimensional dynamical system (we use so ); however, we can reasonably assume that complicated, fractal-type boundaries separate the basins of attraction of the SXA, OPXA/OQPXA, and SCT states. Note that we do not change the initial condition; instead, we change our dynamical system by changing the parameters and and, hence, and Wi. This change affects the long-time evolution of our 2D NS and FENE-P system as sensitively as does a change in the initial conditions because the fractal basin boundary itself changes with these parameters. Reference schneider2007turbulence has found such a sensitive dependence on model parameters in the transition to turbulence and the edge of chaos in pipe flow (see their Fig. 3, which shows that the border, in parameter space, between the laminar and turbulent regions is very complicated, in much the same way as the phase boundary is in our study). In Refs. shajahan2007spiral; shajahan2009spiral, it has been shown that the spatiotemporal evolution of spiral waves of electrical activation depends sensitively on inhomogeneities in partial-differential-equation models for ventricular tissue; changes in the positions of these inhomogeneities in models for cardiac tissue correspond in our system to changes in parameters such as and Wi. A shell-model study of the dynamo transition in MHD has also found such a fractal-type boundary, separating dynamo and no-dynamo regimes, in the plane of the magnetic Reynolds number and the magnetic Prandtl number sahoo2010dynamo.

### iii.3 The case

If , which lies in region , the steady-state solution is , in the absence of polymers. When we add FENE-P polymers to the 2D Navier-Stokes fluid and increase the Weissenberg number Wi, the steady-state solution becomes unstable. For the range of values of Wi in our runs (Table 1), we obtain a steady vortex crystal, which is shown by the pseudocolor plots of and in Figs. 4(a) and (b), respectively, for the illustrative value ; for this set of parameters, Fig. 4(c) shows the reciprocal-space spectrum , which has clear, dominant peaks at the forcing wave vectors.

As we increase Wi beyond , new dynamical regimes appear in the range of values for our runs . In particular, we observe, for and , that the time series of is periodic (Fig. 5 (a)) and its power spectrum (Fig. 5(b)) shows one dominant peak, with hardly any sign of higher harmonics. Thus, the Poincaré-type map in the plane (Fig. 5(c)) displays a simple attractor. In Figs. 5(d) and (e), we show, for these parameter values, pseudocolor plots of and , respectively; there is spatial undulation in the former and the latter is deformed relative to SX (Fig. 1(a)). In the reciprocal-space spectrum (Figure 5(f)), we do not see any major peaks, other than the ones at the forcing wave vectors, because the deformation of the original crystal is weak; the amplitudes of these dominants peaks are smaller than those of their counterparts, at the forcing wave vectors, in Fig.4(c).

An additional increase of Wi (runs with ) leads to chaotic temporal evolution and a disordered set of vortices and antivortices in space. Thus, we have states with spatiotemporal chaos and turbulence; for our nonequilibrium system, turbulent states with spatiotemporal chaos are the analogs of the disordered liquid state, which appears on the melting of an equilibrium crystal ramakrishnan1979first; haymet1987theory; chaikin2000principles; oxtoby1991liquids; singh1991density; perlekar2010turbulence. We illustrate a state with spatiotemporal chaos for and by the plots in Fig. 6. In particular, Fig. 6(a) shows the chaotic energy time series , whose power spectrum (Fig 6(b)) shows a broad background, which indicates temporal chaos. The scattered points in the Poincaré-type section in the plane confirm the presence of chaos. The disorder of this state is shown in the pseudocolor plots of and in Figs. 9(d) and (e), respectively; and the associated reciprocal-space spectrum (Fig. 9(f)) shows many new Fourier modes in addition to the residual peaks at the forcing wave vectors.

Thus, we see that the onset of elastic turbulence can lead to vortex-crystal melting, through a very rich sequence of transitions, even when is itself too small to melt this nonequilibrium crystal. The spatial autocorrelation function and the evolution of the order parameters with Wi are given in the fifth subsection.

### iii.4 The case

If we do not have any polymers and if , which is greater than , the steady-state solution is unstable, and the temporal evolution of our system is chaotic as we can see from the time series of the energy in Fig. 7(a). Its power spectrum , which is plotted in Fig. 7(b), has several peaks in a broad background, and the Poincaré-type section in the (Fig. 7 (c)) displays points that cover a two-dimensional area. The disordered spatial structure of this state, shown by the pseudocolor plots of and in Figs. 7 (d) and (e), respectively, indicates that this state is spatiotemporally chaotic and turbulent. Figure 7 (f) shows the reciprocal-space energy spectrum , which exhibits a large number of Fourier modes in addition to the ones at the forcing wave vectors.

When we add polymers to our 2D Navier-Stokes system at , these polymers lead first to a reappearance of states that are periodic in time; this occurs, e.g., for the range of Wi covered in runs . For , Fig. 8(a) shows that time series of the energy is periodic in time and its power spectrum (Fig. 8(b)) has a fundamental peak at and harmonics at , , , and . The Poincaré-type section in the plane in Fig. 8 (c) exhibits a simple attractor, which is consistent with the periodic behaviors that emerge from Figs. 8(a) and (b). The spatial organization of this state, shown by pseudocolor plots of and in Figs. 8 (d) and (e), respectively, is disordered, although some vestiges of the unit cells of the original crystal are visible in the pseudocolor plot of ; and Fig. 8 (f) shows the reciprocal-space spectrum .

If we increase the Weissenberg number to in run , we obtain a new steady state that is a vortex crystal; however, it is not exactly the same as the original vortex crystal as can be seen clearly by comparing the pseudocolor plots of in Figs. 1(a) and 9(b). This difference is even more striking if we compare the pseudocolor plots of shown in Figs. 9(a) and 2(a). Not surprisingly, then, the reciprocal-space energy spectrum shows new peaks in addition to the ones at the forcing wave vectors, which now appear with a reduced amplitude.

If we increase Wi some more (runs with ), we obtain a spatiotemporal crystal, i.e., one that is spatially periodic and which also oscillates in time perlekar2010turbulence. For instance, if , the time series of the energy is periodic in time, as shown in Fig. 10(a); and Fig. 10(b) shows that its power spectrum has only one dominant peak, which is a clear indication that the temporal evolution is periodic. Figure 10(c) shows the Poincaré-type section in the plane, which displays that the projection of the attractor on this plane is a closed loop. Figures 10(d) and (e) show pseudocolor plots of and , respectively; and the associated reciprocal-space spectrum (Fig. 10(f)) shows dominant peaks principally at the forcing wave vectors.

We now increase Wi again (run ); the temporal evolution of the system becomes chaotic and the spatial organization of the vortex crystal is distorted. Thus, we finally enter a state with spatiotemporal chaos and turbulence (SCT); this is our analog of the disordered, liquid state, which appears on the melting of an equilibrium crystal ramakrishnan1979first; chaikin2000principles; oxtoby1991liquids; singh1991density. However, some vestiges of the incipient crystalline ordering can still be found in reciprocal-space spectra, as we illustrate for in Fig. 11. From the time series of the energy (Fig. 11(a)), it is clear that system is chaotic; this is confirmed by the power spectrum (Fig. 11(b)), which displays several peaks, and the Poincaré-type section in the plane (Fig. 11(c)), in which the points fill a two-dimensional area. The disordered spatial organization of vortical structures is shown in Figs. 11(d) and (e) via pseudocolor plots of and , respectively. The reciprocal-space spectrum (Fig. 11(f)) shows residual peaks at the forcing wave vectors, which indicate the incipient vortex crystal, but the other Fourier modes show that long-range spatial periodicity has been lost in this SCT state.

Thus, our careful study of the case , for various values of Wi, shows how the turbulent SCT state, at , gives way to the states OPXA, SXA, and OPXA states, as we increase Wi (Table 1); finally these states show reentrant melting into the SCT state.

### iii.5 Order Parameters, Spatial Correlation Functions, and Kinetic-Energy Spectra

We now use statistical measures that are employed in the density-functional theory ramakrishnan1979first; chaikin2000principles; oxtoby1991liquids; singh1991density of freezing. In particular, we examine the dependence of the order parameters and the spatial autocorrelation function , defined in Eq. (3), on Wi and . Recall that equilibrium melting is a first-order transition at which jumps discontinuously from a nonzero value in the crystal to zero in the liquid. The melting of our nonequilibrium vortex crystal is far more complicated, in all parts of the parameter space of the phase diagram of Fig. 3, because there are many transitions.

For our nonequilibrium vortex crystal , which is obtained by summing over the four forcing wave vectors, is the analog of the order parameters in a conventional crystal. changes with Wi as shown, respectively, for (a) and and (b) and in Figs. 12 (a) and (b). For and small values of Wi, the steady state of our system is SX, so the Fourier modes at the forcing wave vector are the most significant ones. As we increase Wi, our nonequilibrium system undergoes a series of transitions from SX to the turbulent state SCT (see above); therefore, decreases with Wi as we show in Fig. 12(a). For and , we begin with a system that is disordered and turbulent; if we then increase Wi, our system first goes to the state SX, so first increases as we show in Fig. 12(b); as we increase Wi, the system becomes turbulent and disordered again, so decreases at large Wi (Fig. 12(b)).

We give representative plots of the correlation function (see Eq.(3)) in Figs. 13 (a)-(d). For the crystal, we calculate along the line from and ; this shows periodic peaks [Figs. 13 (a) for ] with widths that are related to those of vortical or strain-dominated regions. In the turbulent phase, we use a circular average of [Figs. 13 (b), (c), and (d) for and , and , and and , respectively]; these peaks decay over a length scale that yields the degree of short-range order. This decay is similar to the decay of spatial correlation functions in a disordered liquid phase in an equilibrium system. For similar, but not the same, correlation functions in a viscoelastic Taylor-Couette flow see Refs. latrache2012transition; latrache2016defect.

We have also obtained the energy spectra in two parts of the turbulent state of our system; the first part is dominated by elastic turbulence and the other one lies in the region of fluid turbulence, where polymers lead to dissipation reduction. For a given resolution, there is one qualitative difference between the energy spectra in these two regimes. In the forward-cascade part of , the slope of the energy spectrum is higher in the elastic-turbulence regime than in the dissipation-reduction regime as we show in Fig. 14 (a), for and (red line with circles), and (green line with squares), and and (black line with triangles). The exponent of the energy spectra is , which is in reasonable agreement with the exponent in laboratory experiments of Groisman and Steinberg groisman2004elastic. Note that, in both elastic-turbulence and dissipation-reduction regimes, the energy spectra show that the energy of the polymeric fluid is spread over many decades of (or, equivalently, over many length scales); this is a well-known signature of turbulent flows, whether they are engendered by large Re or large Wi. We have also obtained the spectrum of polymer stretching, which is defined as , where denotes the time average over the statistically steady state. In Fig. 14 (b) we show the polymer-stretching spectra, for and (red line with circles), and (green line with squares), and and (black line with triangles).

### iii.6 Spatial Patterns of Polymers and Inertial Particles in the Nonequilibrium States of our System

It is instructive to study the spatial patterns of polymers and inertial particles in the nonequilibrium states of our system; we show below that such patterns also provide clear signatures of these states. In Figs. 15 (a), (b), (c), and (d) we show contours of the polymer stretching superimposed on pseudocolor plots of the Okubo-Weiss field for and (crystalline state), and (melted crystal in the elastic-turbulence regime), and (frozen crystal), and and (melted crystal in the dissipation-reduction regime), respectively. As we have seen in Ref gupta2015two, the polymers stretch preferentially in extensional regions, where ; thus, the patterns of their extension mirror the spatial periodicity or lack thereof in the crystalline and turbulent phases, respectively. The spatiotemporal evolution of the plots in Fig. 15 is given in Video crystal-rp-a, b, c, and d at http://www.physics.iisc.ernet.in/rahul/movies.html .

The analogs of Figs. 15 (a), (b), (c), and (d), with polymers replaced by inertial particles (black dots), are given in Figs. 16 (a), (b), (c), and (d), respectively. We see from these plots that inertial particles organize themselves in regions where ; thus, their spatial organization can also be used to surmise the underlying periodicity of a nonequilibrium vortex crystal. This method of visualization has been used in experiments on fluid films without polymers ouellette2007curvature; ouellette2008dynamic. All these plots are for a Stokes number . Note the wavy patterns of particles across the flow field, which are strongly influenced by the forcing field; these patterns are not permanent; and, as the elastic turbulence changes the particle positions, these wavy patterns of particles disappear; similar patterns appear later in time, but at different places in the simulation domain.

The analogs of Figs. 16 (a), (b), (c), and (d), with the field replaced by polymers, are given in Figs. 17 (a), (b), (c), and (d), respectively. We see from these plots that the correlation between the spatial organizations of polymers and inertial particles is not dramatic. This is not surprising insofar as the polymers stretch preferentially where , whereas inertial particles tend to cluster in regions where .

## Iv Conclusions

Our detailed DNS elucidates the nature of melting of a vortex crystal in a nonequilibrium, forced, thin fluid film with polymer additives. This melting can be induced by (a) elastic turbulence at low Reynolds numbers Re (or ) but large Weissenberg numbers Wi (or ), (b) fluid turbulence at low Wi but large , or (c) a combination of these two types of turbulence (see Figs. 3 (a) and (b)). Our work leads to the rich, nonequilibrium phase diagram of Fig. 3. The topology of this phase diagram and the intricate natures of the boundaries between the different phase in it have neither been anticipated nor been studied hitherto.

Our work builds on the DNS study of Ref. perlekar2010turbulence, which uses ideas from the density-functional theory of freezing ramakrishnan1979first; chaikin2000principles; oxtoby1991liquids; singh1991density, nonlinear dynamics, and turbulence to characterize the phases and transitions in the melting, by turbulence, of a nonequilibrium vortex crystal in a 2D fluid film. Most studies of the turbulence-induced melting of a vortex crytal in a fluid film do not include polymer additives; for an overview of such studies we refer the reader to Ref. perlekar2010turbulence. The DNS study of 2D fluid films with polymer additives gupta2015two; berti2008two; berti2010elastic, which uses the Oldroyd-B model, is an exception; it has studied the elastic-turbulence-induced melting of spatially periodic Kolmogorov-flow pattern, which is, in our terminology, a one-dimensional crystal. To the best of our knowledge, there is no study of vortex crystals in fluid films with polymers, which brings together the variety of methods we use to analyze the melting of such crystals.

We follow Ref. perlekar2010turbulence, which does not include polymers, in (a) identifying the order parameters for the vortex crystal, (b) characterizing the series of transitions in terms of and , the energy time series , and their Fourier transforms, and (c) using the spatial correlation functions in crystalline and turbulent phases. However, we go considerably beyond the study of Ref. perlekar2010turbulence by including polymers and, therefore, (a) a new dimensionless control parameter, the Weissenberg number Wi (or ), (b) the polymer-conformation tensor field , and (c) inertial particles. Therefore, we can examine the nonequilibrium phase diagram of this system in a two-dimensional parameter space (see Figs. 3 (a) and (b)) and elucidate the organization of polymers and inertial particles in the flow, both in crystalline and turbulent phases. We show also that our system provides a natural laboratory in which we can study the crossover from elastic-turbulence to dissipation-reduction regimes in the turbulent phase.

We refer the reader to Ref. perlekar2010turbulence for a discussion of the thermodynamic limit and finite-size effects in such systems of nonequilibrium, vortex crystals. No study (including ours) has addressed these issues because, for large system sizes, we need very large DNSs to characterize the temporal evolution of the system. If we extract a correlation length from the correlation function , it is much smaller than the linear size of our simulation domain in the turbulent phase (SCT), so our results should not change if we increase the system size. However, subtle size dependence might occur in the ordered phase because of the inverse cascade, well known in 2D fluid turbulence, which might lead to undulations that give rise to crystals with ever larger unit cells, whose size can be controlled only by the inclusion of friction. A systematic study of such subtle finite-size effects lies beyond the scope of our study.

The sequence of transitions that take us from a vortex crystal to the turbulent state is far richer than conventional equilibrium melting. The former differ from the latter in another important way: To maintain the steady states, statistically or strictly steady, we must impose the force ; in the language of phase transitions, this force is a symmetry-breaking field, both in the ordered and disordered phases. Therefore, there is no symmetry difference between the disordered, turbulent state and the ordered vortex crystal, as we see directly from the vestiges of the dominant peaks in the reciprocal-space spectra ; and the order parameters , with the forcing wave vectors, do not become exactly zero in the turbulent phase; however, they do become very small. Moreover, our nonequilibrium vortex crystal undergoes a sequence of transitions from an ordered state to an undulating crystal and then to a fully turbulent state; there is no external or thermal noise here and hence no fluctuations in the perfect vortex crystal; this has no equilibrium analog.

We hope that our study will encourage experimental groups to try to obtain the rich phase diagram of the melting of a nonequilibrium vortex crystal in a fluid film with polymer additives. In particular, experiments on such a system could look for spatiotemporal crystals, which are periodic both in space and time and which cannot be obtained easily in nonequilibrium, hard-matter settings. Furthermore, the reentrant crystallization, i.e., the recrystallization, with an increase in Wi, of the melted crystal, which emerges from our detailed DNS, is certainly worth exploring in experiments. We have recently learnt of a shell-model study of elastic turbulence that is just published ray2016elastic.

## V Acknowledgments

We thank the Council for Scientific and Industrial Research, University Grants Commission, and Department of Science and Technology (India) for support, SERC (IISc) for computational resources, and P.Perlekar, D. Mitra, S.S. Ray, and D. Vincenzi for discussions. AG acknowledges gratefully funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013)/ERC Grant Agreement N. 279004.