Didier MAQUIN Besma Gaddouna
email: bgaddoun@ensem.u-nancy.fr

José Ragot
email: jragot@ensem.u-nancy.fr

FAULT DETECTION OBSERVERS FOR SYSTEMS WITH UNKNOWN INPUTS

Keywords: Actuator failures, observers, state estimation, failure detection, linear systems

Many approaches for the design of unknown input observers were developed to estimate the state of a linear time-invariant dynamical system driven by both known and unknown inputs. These observers are often used for component failure detection, for instrument fault detection in systems subject to plant parameter variations or uncertainties and for fault detection and isolation for systems subject to unknown inputs. In this paper, a method for the design of fault detection observers is developed for linear time-invariant systems subject to arbitrary unknown inputs. A constructive numerical algorithm using generalized inverse of matrices is included. Moreover, an example illustrates the applicability of the proposed method.

INTRODUCTION

Many approaches for the design of unknown input observers (UIOs) were developed to estimate the state of a linear time-invariant dynamical system driven by both known and unknown inputs [START_REF] Kudva | Observers for linear systems with unknown inputs[END_REF], [START_REF] Kurek | The state vector reconstruction for linear systems with unknown inputs[END_REF], [START_REF] Guan | A novel approach to the design of unknown input observers[END_REF] and [START_REF] Hou | Design of observers for linear systems with unknown inputs[END_REF]. Generally, the estimation is established by eliminating the unknown inputs from the system. This technique may be also applied for singular systems [START_REF] Miller | On designing reducedorder observers for linear time-invariant systems subject to unknown inputs[END_REF], [START_REF] Chen | Robust detection of faulty actuators via unknown input observers[END_REF] and [START_REF] Ragot | Analytical redundancy for system with unknown inputs -Application to fault detection[END_REF]. Our principal aim, in this communication, consists to compare the actual behaviour of a system to a reference one in order to detect and localize the faults when they occur.

Generally, the generation of residuals is issued from the process knowledge; for example analytical models, knowledge base, sensor and actuator signals, ... The most popular approach is based on the analytical redundancies which use the mathematical model of the process. At each moment, the residuals represent the inconsistency between the actual plant variables and the model variables; they are ideally zero but become non-zero if the actual system differs from the ideal one (this may be due to sensor or actuator faults, modelling errors, non exact parameters of the model, ...). When only few sensors are available, the analytical redundancies have to be generated on a long time interval. In this case one commonly use system observers to estimate state time histories for processes with incomplete set of state measurements [START_REF] Gertler | Analytical redundancy methods in fault detection and isolation[END_REF]. Then, the fault detection is achieved by testing the magnitude of these residuals.

The detection and the isolation of sensor and actuator failures have received during the last decade, a lot of attention in the literature in the field of automatic control. Moreover, the use of multiple observers has been extended in order to improve the localization of the faulty elements of a process. The first development in this area is probably due to [START_REF] Clark | A simplified instrument failure detection scheme[END_REF] which uses as many observers as there are outputs, the i th one is driven by the i th output of the process and all the inputs. The same strategy may be applied for inputs [START_REF] Viswanadham | Fault detection using unknown-input observers[END_REF]; the method consists in constructing as many observers as there are inputs, the i th one is driven by the i th input and all the outputs. When it is desired to isolate faults on actuators, some algorithms based on unknown input observers have been proposed; such strategies have also been reviewed by several authors and the reader should refer to the papers given in reference [START_REF] Frank | Robust fault diagnosis using unknown input observer schemes[END_REF] and [START_REF] Chen | Robust detection of faulty actuators via unknown input observers[END_REF]. Some results dealing with the computational aspects were given by O'reilly (1983) who proposed an approach which consists in reducing the problem of the observer matrices computation by expressing all of them in terms of an arbitrary one; [START_REF] Ge | Detection of faulty components via robust observation[END_REF] who treated these aspects for component failure detection using robust observers; [START_REF] Saif | A new approach to robust fault detection and identification[END_REF], for instrument fault detection and identification in linear dynamical systems subject to plant parameter variations or uncertainties. [START_REF] Chang | Fault detection observer design for linear systems with unknown inputs[END_REF] gave also, a solution to the problem of designing an observer for fault detection and isolation for linear time-invariant system subject to arbitrary unknown inputs by using the vector projection theory and propose an algorithm for the matrix calculation. The aim of this paper is also to point out the problem of solving the equations of the observer and to propose an appropriate algorithm for linear systems subject to unknown inputs. This paper is organized as follows: the section 2 deals with the general structure of an unknown input observer giving the constraint conditions that the observer matrices have to satisfy; the section 3 derives an algorithm for the matrix computations from the constraint conditions; the section 4 deals with residual analysis; section 5 presents an application of the design of the unknown input observer in detecting actuator faults and the section 6 illustrates, by an example, the matrix computations and the procedure of actuator failure detection.

UNKNOWN INPUT OBSERVER APPROACH

Consider the linear time-invariant system:

x . (t) = A x(t) + B u(t) + D d(t) (1a) y(t) = C x(t) (1b)
where x(t) ∈ R n is the state vector, u(t) ∈ R r is the input vector, y(t) ∈ R m is the output vector and d(t) ∈ R k+1 (the dimension k + 1 will be justified in the section 5) is the disturbance vector whose elements are unknown functions of t. A, B, C and D are constant matrices with appropriate dimensions; without lost of generality, C is assumed to be full row rank and D to be full column rank.

Our aim is to estimate the state vector x(t) in spite of the presence of the unknown perturbation d(t). So, let us use the full order observer described by:

z . (t) = N z(t) + L y(t) + G u(t) (2a) x ^(t) = z(t) -E y(t) (2b)
where z(t) ∈ R n and x ^(t) ∈ R n is an estimation of the system state vector x(t).

It is desired that, at steady state, x ^ approaches asymptotically x. Therefore, we define the state estimation error:

e(t) = x(t) -x ^(t) (3)
The dynamic equation governing this error can be reduced to e . (t) = N e(t) if one purposely chooses:

P = I + E C (4a) L C = P A -N P (4b) G = P B (4c) P D = 0 (4d) N stable (4e)
Our aim is to establish a systematic procedure in order to determine the matrices defining the observer (2) and verifying the above conditions (4).

Necessary and sufficient conditions for the existence of this type of observers are very classical and will not be discussed here. They are as follows :

rank(C D) = rank(D) = k+1 with m > k+1

THE UNKNOWN INPUT OBSERVER DESIGN

Since C has full row rank, suitable coordinate transformation on the states can be found such that the given system (1) is restricted system equivalent to (O'Reilly, 1983):

x . (t) = A x(t) + B u(t) + D d(t) (5a) y(t) = (I m 0 m(n-m) ) x(t) (5b)
In order to simplify the presentation, the notation used in equation ( 1) remain unchanged although the matrix A and the state x are different. According to the reduced form of C, the following equalities can be deduced:

E C = (E | 0 n(n-m) ) (6a) L C = (L | 0 n(n-m) ) (6b)
These expressions (6a) and (6b) will be used in the following to determine respectively the structure of the matrices P and L. Let us determine the matrices characterizing the observer (4).

Determination of P

We remark, taking into account the equations (6a), ( 4a) and (4d), that (E | 0 n(n-m) ) = P -I and P has to be orthogonal to D. Therefore, we propose the following structure of P which allows us to determine the other matrices:

P =       P 1 0 m.(n-m) P 2 I (n-m).(n-m) ( 7
)
where P 1 ∈ R m.m , P 2 ∈ R (n-m).m . Taking into account (4d) and ( 7), D can be partitioned as follows:

D =       D 1 D 2 (8)
where k+1) , one may write the equation (4d) under the form:

D 1 ∈ R m.(k+1) , D 2 ∈ R (n-m).(
P 1 D 1 = 0 m.(k+1) (9a) P 2 D 1 + D 2 = 0 (n-m).(k+1) (9b)
As D 1 and D 2 are known and rank(D) = k+1, we remark respectively from the equations (9a) and (9b) that m.(m-(k+1)) and (n-m).(m-(k+1)) parameters in P 1 and P 2 can be chosen and that the others can be deduced. More precisely, from the equations (9a) and (9b), the matrices P 1 and P 2 may be computed as:

P 1 = K 1 (I -D 1 D 1 + ) (10a) P 2 = -D 2 D 1 + + K 2 (I -D 1 D 1 + ) (10b)
where D 1 + is a generalized inverse matrix of D 1 , K 1 and K 2 are arbitrary matrices (we will define this choice later) with K 1 ∈ R m.m and K 2 ∈ R (n-m).m . If D 1 is a non singular matrix then P 1 is null and P 2 = -D 2 D 1 -1 .

Determination of E

From (4a) and using (7), one may write:

(E | 0) =       P 1 -I m.m 0 m.(n-m) P 2 0 (n-m).(n-m) (11)
From which we deduce:

E =       P 1 -I m.m P 2 (12)

Determination of N

The matrices A and N can be respectively partitioned as:

A =       A 11 A 12 A 21 A 22 N =       N 11 N 12 N 21 N 22 with A 11 ∈ R m.m , A 21 ∈ R (n-m).m , A 12 ∈ R m.(n-m) , A 22 ∈ R (n-m).(n-m) , N 11 ∈ R m.m , N 21 ∈ R (n-m).m , N 12 ∈ R m.(n-m)
and N 22 ∈ R (n-m).(n-m) . The matrices P A and N P are then expressed as: From ( 13) and ( 14) and taking into account the expression of C, the following equations have to be satisfied:

P A =       P 1 A
P 1 A 12 = N 12 (15a) P 2 A 12 + A 22 = N 22 (15b)
Taking into account (10a) and (10b), equations ( 15a) and (15b) can be written:

K 1 (I -D 1 D 1 + ) A 12 = N 12 (16a) A 22 + (-D 2 D 1 + + K 2 (I -D 1 D 1 + )) A 12 = N 22 (16b)
So, with the above expression ( 16) of N 12 and N 22 , the matrix N can be written:

N =       N 11 K 1 (I -D 1 D 1 + ) A 12 N 21 K 2 (I -D 1 D 1 + ) A 12 + A 22 -D 2 D 1 + A 12 (17) Otherwise: N = R + (N _ K) Q = R + F Q (18)
where:

R =       0 0 0 A 22 -D 2 D 1 + A 12 N _ =       N 11 N 21 Q =       I 0 0 (I -D 1 D 1 + ) A 12 K =       K 1 K 2
The matrix F is chosen such that N, which fix the dynamic of the observer, has a desired set of stable eigenvalues. Consequently, a simple partition of F gives the matrices N _ and K.

Determination of L

The expression of L can be deduced from (4b). Indeed, taking into account (6b), this equation is reduced to:

(L | 0) = P A -N P ( 19
)
Summarizing the different matrix computations the design of the observer is achieved by the following algorithm:

Algorithm System equations x . (t) = A x(t)+ B u(t) + D d(t) y(t) = C x(t) Observer equations z . (t) = N z(t) + L y(t) + G u(t)
x ^(t) = z(t) -E y(t)

-Partition D and A (according to the dimension of y and d)

-Compute R =       0 0 0 A 22 -D 2 D 1 + A 12 and Q =       I 0 0 (I -D 1 D 1 + ) A 12 .
-Determine F by the method of eigenvalues assignment such that N = R + F Q has a desired set of stable eigenvalues.

-Deduce K (K 1 and K 2 ) and N _ (N 11 and N 21 ) .

-

Compute P 1 = K 1 (I -D 1 D 1 + ) P 2 = -D 2 D 1 + + K 2 (I -D 1 D 1 + ).
-Deduce P, N 12 and N 22 .

-Compute G = P B and

E =       P 1 -I m.m P 2 .
-Deduce L from (L | 0) = P A -N P.

Remark : according to all the previous partitions of matrices, it is easy to verify that if K 1 is a null matrix, the first m states of z do not depend on u and y. That implies that output estimations are exactly equal to the measurements.

ANALYSIS OF THE RESIDUALS

Let us analyse the aptitude of the observer to reconstruct the state of the process and its ability to detect faults on actuators and sensors. In the presence of sensor faults p(t) and actuator faults a(t), the system (1) is equivalent to the one shown by the fig. 1.

a(t) u(t) System Observer up(t) d(t) p(t) y(t) yp(t)
y(t) Fig.

UIO in the presence of faults

The signals u p (t) and y p (t) are respectively the perturbed input and output. In this condition, the system (1) and the observer (2) become respectively:

x . (t) = A x(t) + B u(t) + D d(t) + F a a(t) (20a) y(t) = C x(t) + F p p(t) (20b) z . (t) = N z(t) + L y(t) + G u(t) (21a) x ^(t) = z(t) -E y(t) (21b)
where F p and F a are matrices with appropriate dimensions. The output estimation error is defined by:

ε(t) = y(t) -y ^(t) (22)
Taking into account the remark formulated in the previous section, it is obvious that, if K 1 is a null matrix, the output estimation error will not be very useful for detecting actuator faults. So, in order to design efficient fault detection observers it is necessary to introduce a supplementary constraint: the computation of the matrix F using eigenvalues assignment must be done taking into account that the matrix K 1 must be different from a null matrix. Therefore, using Laplace transform, equation ( 22) can be expressed as:

ε(s) = T p p(s) + T a a(s) ( 23
)
where:

T p = (I -C (s I -N) -1 L) F p and T a = C (s I -N) -1 F a
We remark that this estimation error is insensitive to the unmeasured perturbation d(t) and depends on the two types of faults which can be correctly detected if the elements of the matrices T a and T p are not null and isolated if the columns of the matrices T a and T p are independent. Thus, once the observer has been designed, it is then easy to test the isolability condition of faults. In fact, the failure signatures are closely related to the structure of the matrices T a and T p As suggested by [START_REF] Gertler | An evidential reasoning extension to quantitative model-based failure diagnosis[END_REF], equation ( 23) may be characterized by an occurrence matrix which contains 0 and 1, at any given position, whether the matrices T a and T p contain zero or non zero elements. The examination of the columns of the occurrence matrix allows us to define which fault can be detected and isolated.

APPLICATION TO THE ISOLATION OF ACTUATOR FAULTS

Our aim here is to detect and identify actuator faults using a bank of UIOs. The detection and identification scheme is constructed to be sensitive to these faults and to isolate them [START_REF] Chen | Robust detection of faulty actuators via unknown input observers[END_REF]. The scheme is not affected by uncertainties (parameter variations, process noises, ...): their influences are eliminated from the state estimation error [START_REF] Viswanadham | Fault detection using unknown-input observers[END_REF]. So, the idea is to associate a part of the input u to the uncertainties (these inputs will be considered as unknown inputs) and to construct estimators which are insensitive to failures in the actuators corresponding to the unknown inputs. We can formulate this proposition by writing the system (1) as follows:

x

. (t) = A x(t) + B i u p i (t) + B i u pi (t) + H v(t) i = 1, ..., r (26a) y
(t) = C x(t) (26b)
where v ∈ R k is an unknown input (the dimension of v justifies that of the vector d(t) given in (1a)), H is its distribution matrix, B i is the i th column of B and B i is the n(r-1) matrix obtained from B by deleting B i , u i is the i th entry of u and u i is the input obtained from u by deleting u i .

Denoting d i = (u pi v T ) T and D i = (B i H), the system (26) becomes:

x

. (t) = A x(t) + B i u p i (t) + D i d i (t) (27a) y(t) = C x(t) (27b)
where

D i ∈ R n.(k+1)
As the estimation of the system state is not affected by failures in the unknown inputs, we assume that rank(D i ) = k + 1 and that rank(C D i ) = k + 1 (see section 3) to construct r UIOs given by:

z . i (t) = N i z i (t) + L i y(t) + G i u i (t) (28a) x ^ i (t) = z i (t) -E i y(t) (28b)
where

N i ∈ R n.n , L i ∈ R n.m , G i ∈ R n.(r-1) , E i ∈ R n.m for i = 1, .
.., r. The matrices of each observer that might satisfy the following conditions

P i = I + E i C (29a) L i C = P i A -N i P i (29b) G i = P i B i (29c) P i D i = 0 (29d) N i stable (29e)
are computed according to the algorithm given in the third section.

Then, we construct a set of residuals (ε i (t) = y(t) -y ^ i (t)) using the proposed elimination strategy by considering u i as an unknown input. These residuals are used to detect and identify the actuator faults. By nature of construction, the residuals are insensitive to u i whereas variations and failures in u i will affect the outputs y(k).

By monitoring these sets of residuals, the actuator failures can be isolated. Indeed, if all sensors are well functioning and all residuals are statistically null except the i th one, then the i th actuator is faulty (fig. 2). The residual ε i (t) may be used to detect and identify the faulty actuators by comparing its magnitude with a threshold. Let us denote s i the signature of the residual ε i (t), we can write:

ε i (t)  < threshold: no failure => s i = 0. (30a) ε i (t)  > threshold: the i th actuator has failed => s i = 1.

(30b)

If the failure affects the known input u i then the residual is sensitive to this anomaly; on the other hand, when the failure affects the unknown input u i , the residual remains close to zero. If two or more actuators fail simultaneously, the above decision logic is not able to isolate faults; for such situation, [START_REF] Chen | Robust detection of faulty actuators via unknown input observers[END_REF] have proposed a scheme in which each UIO is driven by r-t inputs if t actuators have simultaneously failed.

NUMERICAL RESULTS

In this section, we give a numerical example to illustrate the technique of the matrix computation and the detection and isolation of actuator faults. The system is driven by two known inputs u 1 and u 2 (figures 3a and 3b) and one perturbation; its state has four components and the dimension of the measurement vector is equal to three (figure 4). The matrices of the system (27) are: A fault with constant magnitude has been added to the first input between the instants 100 and 150. For this example, the eigenvalue assignment problem has enough degrees of freedom. So, the matrix K may entirely be imposed. In the following, we have chosen:

A =         1.
C =         1 0 0 0 0 1 0 0 0 0 1 0 H =         1
K =         1 0 1 1 1 0 1 1 1 0 1 1 .
One of the two inputs is alternately considered as unknown. In both cases, the observer matrices and the output error estimation ε i have been determined and computed.

-Case 1: u We remark from the figure 5 that an anomaly occurs between the moments 100 and 150. The fault has affected all the components of the residual. This is due to the fact that P 1 and P 2 are not null; so, to the choice of K.

p i = u p 1 N 1 =         0
-Case 2: u p i = u p In this case the residual is null. It can be interpreted by the fact that u 2 is not affected by faults and that the anomaly is associated to the unknown input. So, we can conclude that the first actuator has failed between the moments 100 and 150.

CONCLUSION

We have developed a method which allows one to determine systematically the matrices describing a fault detection observer in the presence of unknown inputs. The detection and isolation of faulty actuators have been done using a classical extension, taking into account unknown inputs, of the generalized observer scheme. This method has been experimented on various different systems and seems to give interesting results.

  Fig. 2. The decision logic

  Fig. 3a. First input u 1

  Fig. 5. Output error estimation ε 1 with u p i = u p 1


P Clark

Di Leoni

A Alexakis

L Biferale

M Buzzicotti

Phase transitions and flux-loop metastable states in rotating turbulence

By using direct numerical simulations of up to a record resolution of 512x512x32768 grid points we discover the existence of a new metastable out-of-equilibrium state in rotating turbulence. We scan the phase space by varying both the rotation rate (proportional to the inverse of the Rossby number, Ro) and the dimensionless aspect ratio, λ = H/L, where L and H are the sizes of the domain perpendicular and parallel to the direction of rotation, respectively. We show the existence of three turbulent phases. For small Ro but finite λ, we have a split cascade where the injected energy is transferred to both large and small scales. For large λ and finite Ro there is no inverse cascade and the energy is transferred forward in Fourier space only. Surprisingly, between these two regimes, a third phase is observed as reported here for the first time. Consequently, for certain intervals of Ro and λ, energy is no longer accumulated at arbitrarily large scales, rather it stops at some characteristic intermediate length-scales from where it is then redistributed forward in Fourier space, leading to a flux-loop mechanism where the flow is out of equilibrium with vanishing net flux, and non-vanishing heterochiral and homochiral sub-fluxes. The system is further characterized by the presence of metastability and critical slowing down, explaining why previous experiments and numerical simulations were not able to detect this phenomenon, requiring extremely long observation time and huge computational resources.

Introduction. Statistical systems can develop critical behaviour, where abrupt macroscopic changes happen when varying some control parameter, like temperature or magnetic field [START_REF] Goldenfeld | Lectures on phase transitions and the renormalization group[END_REF]. Averaged quantities can show discontinuous or continuous variations across the critical lines/points in the phase space where the transition occurs. In many cases, experimental and numerical realizations are affected by long transients that are generated due to the presence of metastable states corresponding to local minima of the (free) energy. There have been many attempts to transfer such descriptions to out-of-equilibrium systems quantitatively [START_REF] Zwanzig | Nonequilibrium statistical mechanics[END_REF][START_REF] Berthier | [END_REF]. The critical behaviour of stationary systems in the presence of energy injection mechanisms, dissipation and nonvanishing fluxes remains a major topic of current research in fluid-dynamics, granular and active matter and lacks systematic theoretical understanding [4][5][6]. A paradigmatic example of (phase) transition is the sudden jump from laminar to turbulent dynamics that is observed in Poiseuille and Couette flows when changing the forcing intensity [7][8][9]. In this paper, we investigate the important case of rotating turbulence where the control parameter, given by the intensity of the Coriolis force, affects the symmetries of the macroscopic flow but not the energy injection rate [START_REF] Davidson | Turbulence in rotating, stratified and electrically conducting fluids[END_REF][START_REF] Campagne | [END_REF][12][13][START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF]. In this set-up, it is known that for sufficiently weak rotation rate, Ω, the system behaves as 3D homogeneous and isotropic turbulence transferring energy to small scales only (forward energy cascade), while for Ω above a critical value, Ω c , 3D fluctuations are sufficiently suppressed, the flow becomes quasi-two dimensional and energy is transferred with a split cascade to large scales also [5,[START_REF] Smith | [END_REF]. As as the domain size H in the direction of rotation becomes larger, Ω c increases [16][17][18]. Arguments based on wave-turbulence theory suggest that the inverse flux vanishes in the infinite volume limit, predicting that lim H→∞ Ω c → ∞ [19]. More recently, considering an asymptotic form of the governing equations for large Ω, it was shown that Ω c ∝ H [20]. However, for finite Ω we do not know the precise functional dependence of the boundary Ω c (H) nor the nature of the transition. Despite of the importance of such questions for geophysical and engineering applications we still do not have a satisfying understanding of any of them. While trying to address these problems, we discovered that the physics is even richer due to the existence of a new region in the (λ, Ro) phase-space, where turbulence develops a third, non-trivial, macroscopic out-of-equilibrium state, characterized by a flux-loop cascade, where the flow organises to spontaneously contain the tendency to condensate energy in larger and larger scales, entering a stationary regime and producing a quasi-ordered array of turbulent columnar vortices, akin to a vortex crystal [START_REF] Aref | Vortex crystals[END_REF]. The aim of this letter is to report about this new state of turbulent rotating flows, characterizing its peculiar statistical and dynamical properties, including metastability and critical slowing down as well as commenting about its potential importance for other turbulent realizations.

Setup. We begin by considering the flow in a rectangular periodic domain with aspect ratio λ ≡ H/L and dimensions 2πL × 2πL × 2πH in a rotating reference frame where the rotation is along the direction with dimension 2πH. The governing equations for the incompressible velocity field, u, can be written as

∂ t u + u • ∇u + 2Ω × u = -∇P + ν∆u + f (1)
where ν is the kinematic viscosity, f is an external forcing and 2Ω × u is the Coriolis force produced by the frame rotating with intensity Ω. We choose the forcing f to be Gaussian and delta-correlated in time, acting on a spherical Fourier shell with |k| ∈ [k f , k f + 2] and with amplitude f 0 = f 1 /λ, so that the total injection rate = f • u remains fixed, where • means an average over the forcing realization or on the whole fluid volume. In order to reduce viscous effects, most of the results shown have been obtained by adopting a modified hyperviscosity (see later and SM for details). In the present work, we keep k f L = 20 fixed, using λ as control parameter. Besides λ the other two non-dimensional quantities are given by Reynolds, Re = 1/3 k -4/3 f ν -1 and Rossby,

Ro = 1/3 k 2/3
f Ω -1 , numbers. The equations are solved using a parallel pseudo spectral code (see details in [START_REF] Biferale | [END_REF]) using grids as big as 512 × 512 × 32768 for the largest aspect ratio λ = 64. It is important to stress that accessing high aspect ratios is key to attack the infinite volume limit in the direction parallel to rotation and to assess potential singular effects induced by a finite separation of the 2D plane at k = 0 from the 3D modes with k > 0 in Fourier space [23]. We also fixed L = 1 and = 1.

RESULTS. Figure 1(a) summarizes the main results of our paper, showing the existence of three different macroscopic phases of the rotating flow in the (λ, Ro) space, consisting in (i) a pure forward cascade, (ii) a new flux-loop regime (the choice of name will become obvious later) and (iii) a split cascade regime. In Fig. 1b we show for the sub-set of simulations with λ = 1 and at various Ro the isotropic energy spectra, defined as

E k (t) = 1 2 k≤|k|<k+1 | ûk (t)| 2
where ûk (t) are the Fourier coefficients of the velocity field. The direct cascade regime (black dashed lines) does not develop any large scale fluctuations and peaks at the forcing scale. The split cascade regime (blue dashdotted curves) showcases both forward and inverse cascades and the simulations are stopped when the peak at the largest horizontal scale, k 1, is well developed.

The novelty here is given by the flux-loop phase (solid green lines) showing an intermediate spectral behavior. In this case, the energy spectra are much more irregular, break self-similarity and have their largest peak at an intermediate wavenumber, k 5. Figure 1(c) shows the evolution of the total energy, E(t) = k E k (t), for some of the most characteristic (λ, Ro) values. As expected, in the case with a forward cascade we have a constant (small) total kinetic energy as all energy input is dissipated at high wavenumbers by the viscosity. In the two data-sets showing a split cascade regime, the energy increases constantly as it is transferred to large scales without important dissipative effects. In the three flux-loop data-sets, the total energy saturates for very long times, indicating that the flux to the large scales is halted. In Figs. 2(a)-(c) we show visualisations of the vorticity projection in the direction of the rotation axis for three characteristic data-sets representing the three different phases at late times, and in the SM we show movies comparing their time evolutions. In the split cascade regime (a), the system forms many co-rotating columnar vortices which eventually merge into one. In the forward cascade regime (c), no large scale coherent vortical structures are formed, as expected. In the new flux-loop regime (b), the columnar vortices form, but do not merge and get quasi stuck in a lattice-like structure that persists in time. Similar structures, deemed "vortex crystals", have been observed in systems like 2D point vortices [24,25], 2D turbulence [26][27][28], Bose-Einstein condensates [29,30], and even Jupiter's atmosphere [31].

In particular, asymmetric states where the system resembles a crystal with defects, like what is shown in Fig. 2b, have been shown to be equilibria of co-rotating point vortex systems [24,25]. Note that in 2D randomly forced turbulence there is a symmetry between positive and negative vorticity. As a result, the kind of structures that we observe here are connected to the asymmetry between co-rotating and counter rotating vortices introduced by rotation and 3D effects.

It is worth noting that the vortex crystal state is formed in the absence of any large scale damping term to suppress the inverse cascade. The stationarity of the energy spectrum then implies that the total inverse energy flux at k < k f has to be zero. Nonetheless, the spectrum is far from the E k ∝ k 2 shape predicted by a simple equilibrium distribution [32,33]. To resolve this puzzle we show in Figs. 2(d)-(f) the total energy flux, Π k = -i |k|≤k p+q=k ( û-k • ûp )(k • ûq ), and its exact decomposition in homochiral and heterochiral sub-components, Π k = Π hom k + Π het k , built in terms of Fourier triads including modes with the same or opposite helicity signatures (see SM and [23,34] for a discussion about the importance of hetero-and homo-chiral properties for the energy cascade direction). In order to reduce fluctuations, fluxes are averaged on stationary or quasistationary time windows. As one can see, the forward cascade (panel f) does not show any exotic behaviour, both the total flux and it sub-components are zero in the k < k f range [35,36]. For the split cascade phase in the range k < k f (panel d) we have the usual negative total flux as a result of the negative contributions from the two helical sub-components (in a quasi 2D regime helicity does not play any role and homo-and hetero-chiral channels are expected to be identical [37,38]). The interesting -and non trivial-result is shown in panel (e) where the total flux for k < k f is zero, as it must be if the statistics are stationary, but it is the result of a balance between the forward, Π het k > 0, and the inverse, Π hom k < 0 sub-fluxes contributions. Hence the name of a flux-loop state [5]. This highly intricate flux-loop balance is an out-of-equilibrium effect and has already been observed in two-dimensional but three-component flows [39] and in rotating flows with only three-component mo- tions [40] where similar peaked spectra were found.

From previous figures it is clear that the inverse cascade and the flux-loop phase are the results of a competition between a tendency to become 2D-like contrasted by some residual 3D structures that push energy forward. It is therefore interesting to assess the dynamical effects of these contrasting forces. To do that, we measured the typical time it take the energy to become concentrated on the k ≈ 0 2D plane, defined as the instant of time, t * , when the ratio E (2D) (t * )/E(t * ) hits 0.5, where

E (2D) (t) = |k•Ω|<1 1 2 | ûk | 2
is the total energy in the Fourier plane perpendicular to the rotation direction. In the inset of Fig. 3 we show the evolution of E (2D) (t)/E(t) as a function of time for four different values of Ro. It is clear that in both the split and the flux-loop cascades, the flow approaches asymptotically a quasi-2D state, the main difference is the time it takes to reach it. The main panel of Fig. 3 shows the time t * /τ f vs Ro in a log-lin plot at fixed λ, where τ f = ( k 2 f ) -1/3 is the characteristic time associated with the forcing. It is interesting to notice that when entering the flux-loop region, t * increases with Ro faster than exponentially, indicating the presence of a possible divergence at the critical value of Ro c 0.65. The transition from the split to the forward cascade as a function of the Rossby number has been also analyzed in a recent important study [18], but with runs that evolved until t/τ f = 30 only, and thus miss the development of the flux-loop regime. In fact, the critical Rossby number reported in [18] is around the same value where we see the transition from split to flux-loop cascades, and the growth rates they report (measured by the rate of change of the vertical correlations) depend exponentially on Ro, the same way t * does for Ro < 0.4.

Two important remarks are now in order. First, the transition between split and flux-loop cascades exhibits hysteresis. If one takes a stable simulation under the flux-loop regime and sufficiently decreases Ro the vortexcrystal is destabilized, the columns merge, and the system switches to a split cascade regime. Conversely, the opposite does not happen. We have checked this explicitly. Second, and even more important, the vortex crystals formed in the flux-loop regime are metastable. For example, we observed that by evolving a simulation with parameters (λ, Ro) close to the transition from flux-loop to split cascade the vortex crystal structure is destabilized after a long time and the inverse cascade starts again (see the central panel in the SM movie). Metastability is consistent with the fact that the vortex crystal is observed at the boundary between the inverse and the forward cascade and with the existence of a critical slowing down.

Conclusions. By using huge high-performancecomputing resources we have studied the (λ, Ro) phase space of rotating turbulence up to record resolution of 512 × 512 × 32768 grid points. We found the existence of a new metastable flux-loop regime, where the inverse energy cascade is stopped by a delicate balance between hetero-and homo-chiral triadic non-linear interactions, leading to multiple metastable vortex-crystal like states. These states are stable for very long times but can transition to the inverse cascade regime if perturbed strongly enough.

Observations of multiple large-scale, self-organised turbulent state are becoming more and more common in the turbulent literature, having been observed in bounded flows [41], in anisotropic sheared turbulence [42], in swirling flows [43] and in magnetohydrodamic flows [44]. They play a key role for both fundamental aspects, suggesting the existence of multiple attractors in the system and applied ones, leading to huge variations in the global energetic balance for tiny changes in the control parameters. The present study relates these metastable states with the boundaries of two different states (forward and split cascading) in a phase-diagram making the connection with classical phase transitions. Some important questions remain open, connected to the asymptotic behavior of the critical lines between different phases and the robustness of the metastable states in the limit Ro → 0, λ → ∞ and λ → 0 (see dashed red lines in the phase-space summary of Fig. 1(a)). The presence of the flux-loop condensate prevented us from extracting a precise functional behaviour in the (λ, Ro) plane for the transition from forward to split cascade. Although the asymptotic scaling λ ∝ Ro -1 that was suggested in [5] and found in [20] is still plausible, the present results can not confirm nor reject this hypothesis, because of the presence of the flux-loop phase. Finally, we note that metastable properties can depend on the horizontal domain size L. If L is large enough and the size of the vortex-crystal is increased the possibility of having a destabilizing local defect increases.

with s k = ±. The homo-chiral energy fluxes is made out of triads with the same chirality: ( û±

-k • û± p )(k • û± q ) (6)
while the heterochiral is given by all resulting triads with two Fourier modes of opposite chirality:

Π het k = Π k -Π hom k , (7)
where Π k is the total energy flux defined in the main text as

Π k = -i |k|≤k p+q=k ( û-k • ûp )(k • ûq ) . ( 8
)
Video material

A visualization of an inverse and a flux-loop cascade can be found in uploaded video, where we show the volume rendering of the time evolution of the vorticity component along the direction of rotation axis, ω z , for three different simulations obtained with (Ro, λ): (0.25, 1) split cascade (left panel); (0.51, 1) split/flux-loop (center panel); and (0.58, 1) flux-loop (right panel). Notice the metastable regime shown in the center panel where the inverse cascade first stops in to a quasi vortex crystal state and then suddenly restarts, thanks to 3D vortex merging. In the bottom row of the same video we present the energy spectra and total energy evolution for the same three simulations.

FIG. 1 .

 1 FIG. 1. (a): (λ, Ro) Phase space. Different symbols represent the three macroscopic turbulent cascade phases: forward, flux-loop and split. Red continuos line are a guide to the eye to distinguish the three phases. Red dotted-lines represents possible asymptotics behaviour in the limit Ro → 0 and λ → ∞ (see discussion in the text). (b): Instantaneous energy spectra for fixed aspect ratio, λ = 1 at different Rossby numbers 0.1 ≤ Ro ≤ 1. For flux-loop and direct cascade cases the spectra are plotted in the stationary regime, for the split-cascade regime we used the final time when we stopped the simulation. (c): time evolution of the total energy, E(t), for some characteristic (λ, Ro) values (represented with empty symbols in panel (a). Line colours distinguish the three phases following the same color code of symbols in panel (a).

FIG. 2 .

 2 FIG. 2. (a)-(c): Visualizations of the parallel component of the vorticity for a split cascade case, a flux loop case, and a forward cascade case, respectively. (d)-(f): Total and chirally decomposed fluxes, for the three regimes shown in the top row.

FIG. 3 .

 3 FIG. 3. Time to bidimensionalization t * as a function of Rossby number for aspect ratio λ = 1. Inset: Ratio of the energy in the 2D modes a function of time.

  The authors acknowledge partial funding from the Eu-ropean Research Council under the European Community's Seventh Framework Program, ERC Grant Agreement No. 339032. The simulations have been done using resources provided through the PRACE initiative Pra17_4374 at CINECA. The authors acknowledge Pablo Mininni and Gregory Eyink for useful discussions. where h ± (k) are the orthogonal eigenmodes of the curl operator, hence each Fourier modes of the velocity field satisfies ik × ûs k k = s k ûs k k ,

 I. All simulations were produced with a horizontal resolution of Nx = Ny = 512, box size L = 2π, hyperviscosity of order α = 2, forcing wavenumber in the range range [k f , k f + 2] with k f = 20, forcing amplitude f0 = f1/λ with f0 = 1.66, viscosity ν = 4 × 10 -7 , energy injection rate of = 2.02 and Reynolds number defined on the forcing scale Re = 125 (see main text). The parameters shown in the table are the box aspect ratio λ = H/L, the number of collocation points in the vertical direction Nz, rotation rate Ω, Rossby number defined in terms of the energy injection properties, Ro = (ε f k 2 f ) 1/3 /Ω, and the regime the simulation is in.

SUPPLEMENTAL MATERIAL Numerical Simulations

We performed a series of direct numerical simulations (DNS) of the incompressible Navier-Stokes equations (see Eqs. [START_REF] Goldenfeld | Lectures on phase transitions and the renormalization group[END_REF] in the main text) where the normal viscosity has been replaced with hyper-viscous dissipation:

with α = 2. The motivation to use hyper-viscosity comes from request to minimize viscous effects at the forcing scale at the expense of generating spectral bottlenecks at high wavenumbers [45], allowing the inverse energy cascade to develop with minimal Reynolds effects. For a few cases, we made sure to test that increasing the Reynolds numbers do not introduce important differences with what reported here. In Table I we provide a list of all parameters used in each of the simulations presented in this paper.

Helical decomposition

In this section we provide the definition of the helical-fluxes calculated in Fig. 2(a-c) of the main text. For this purpose we exploit the decomposition of any incompressible 3D flow into helical modes proposed by [34,46]. From the incompressibility assumption it follows that u(x) is a solenoidal vector field, hence its Fourier modes û(k) depends only on two linearly independent degrees of freedom and we can decompose the velocity field as follows


Nicolás P Müller

Marc-Etienne Brachet

Alexandros Alexakis

Pablo D Mininni

Abrupt Transition between Three-Dimensional and Two-Dimensional Quantum Turbulence

whether they are published or not. The documents may come

L'archive ouverte

We present numerical evidence of a critical-like transition in an out-of-equilibrium mean-field description of a quantum system. By numerically solving the Gross-Pitaevskii equation we show that quantum turbulence displays an abrupt change between three-dimensional (3D) and two-dimensional (2D) behavior. The transition is observed both in quasi-2D flows in cubic domains (controlled by the amplitude of a 3D perturbation to the flow), as well as in flows in thin domains (controlled by the domain aspect ratio) in a configuration that mimics systems realized in laboratory experiments. In one regime the system displays a transfer of the energy towards smaller scales, while in the other the system displays a transfer of the energy towards larger scales and a coherent self-organization of the quantized vortices.

The phenomena of condensation and phase transitions in statistical mechanics has traditionally been associated with equilibria. However, observations of turbulence in experiments of gaseous Bose-Einstein condensates (BECs) [1][2][3] and of superfluid 4 He [4][5][6] have shown that these out-of-equilibrium systems can also display multiple phases. In particular, recent BEC experiments close to a two-dimensional (2D) regime, instead of a tendency towards disorder, display an intriguing out-ofequilibrium self-organization and the nucleation of quantized vortices [7][8][9] (see [10,11] for numerical studies).

In classical turbulence, a reminiscent process can take place when flows are 2D. Under certain conditions, the kinetic energy can undergo an inverse cascade (moving to larger scales), and eventually create a condensate [12]. This condensation is of a different nature than a BEC as it involves the kinetic energy of the system instead of its mass density. In classical three-dimensional (3D) turbulence, recent developments indicate that this farfrom-equilibrium system can change its behavior as its dimensionality is changed [13][14][15][16] (or, equivalently, as one of its spatial dimensions is compactified, see [13], and [17] for an example of a transition under compactification in gravitational theories). In classical fluids, when the flow is 3D energy undergoes a direct cascade (moving to smaller scales), while as the domain that contains the fluid is made thiner, the system becomes 2D and displays an inverse cascade after a critical second-order transition.

Both classical and quantum turbulence involve nonlinear and complex spatio-temporal dynamics of fields, and cascade-like solutions can develop in many different cases. In this letter we address the following questions: Is there a transition in the behavior of quantum turbulence as its dimensionality is changed as reported in recent quantum turbulence experiments [7][8][9]? And is this transition associated with the emergence of different outof-equilibrium self-similar regimes (i.e., associated with a change in the direction of the energy cascade)? To this end, we study numerically 3D condensates in periodic boundary conditions using the Gross-Pitaevskii equation (GPE), exploring two configurations. In one, we solve the equations in a cubic domain and perturb an initial 2D random array of quantized vortices with a 3D perturbation, varying the amplitude of the perturbation as a control parameter. In the other, we consider a quasi-2D array of quantized vortices and vary the aspect ratio of the domain, compactifying one of its dimensions. In both cases we find evidence of an abrupt transition towards a regime that displays two-dimensionalization, spatial aggregation of quantized vortices, and inverse energy flux.

To describe the dynamics of weakly interacting bosons of mass m at zero-temperature we solve numerically the GPE [18], i ∂ t ψ = -2 ∇ 2 ψ/(2m) + g|ψ| 2 ψ, where ψ is the condensate wave function and g is proportional to the scattering length. The fluid density, velocity, and quantized vortices can be obtained from ψ using Madelung's transformation (see [18]). The GPE is solved using a parallel pseudospectral method [19,20]. To achieve the largest possible scale separation (at a fixed spatial resolution), we resort to periodic boundary conditions in a 3D domain of size

L x × L y × L z , with spatial resolution N x × N y × N z . The size of the domain is L x = L y = L ⊥ = 2π in dimensionless
units in all cases, and L z = γL ⊥ where γ is the domain aspect ratio. In these domains, we prepared a set of randomly distributed 2D vortices with a small 3D perturbation of amplitude A z , such that the wavefunction is a solution of the GPE, and that the incompressible kinetic energy of the system peaks at an intermediate wavenumber k 0 ≈ 10 (i.e., the correlation length of the flow is 0 ≈ L ⊥ /10; see [18] for more details on the preparation of the initial conditions and for the definition of the incompressible kinetic energy). This results on quantized vortices with a random separation, and that are perfectly 2D for A z = 0 while they display stronger curvature in z for increasing A z .

10 -5 10 -3 E i k (k) (a) A z =0 A z =0.0003 A z =0.0006 A z =0.001 A z =0.01 A z =0.05 A z =0.4 A z =1
As previously mentioned, we consider two ways to observe a transition between 2D and 3D flows using these initial conditions. One of them consists on varying the amplitude of the 3D perturbation A z between 0 and 1 in a cubic domain. The other, is to vary the aspect ratio of the domain for fixed A z (the 2D limit case being that in which γ = L z /L ⊥ → 0, and the 3D case when γ = 1). In each case, when varying the control parameters between their two limits, classically we can expect an inverse cascade of energy in the 2D regime, and the absence thereof in the 3D case. To identify the direction of the cascades we consider two quantities: the incompressible kinetic energy spectrum E i k (k) (see [18,21,22] for a detailed description of energy components in the GPE) and the total energy flux Π(k) = -dE < (k)/dt, where E < (k) is the total energy of the system integrated up to wavenumber k, E < (k) = k 0 E(k )dk , and where E(k) is the total energy spectrum [18]. A direct cascade of energy corresponds to the development of a power law in E i k (k) for k > k 0 and with Π(k) > 0 constant in a range of wavenumbers, while an inverse cascade of energy corresponds to a growth of E i k (k) for k < k 0 and with Π(k) < 0. As the system has no external steering force (but no dissipation either), an inverse cascade can only develop for a transient time [23], and in the following we will focus on time averages of these quantities between t = 1 and 10 flow turnover times, as well as on their time evolution over the same time span (with the turnover time defined as τ = 0 /U , with U the r.m.s. initial flow velocity).

In cubic domains (γ = 1) we performed two sets of simulations, with spatial resolutions of N x × N y × N z = 256 3 and 512 3 grid points, varying the amplitude of the 3D perturbation A z . For large values of A z the flow quickly evolves into a 3D regime, with quantized vortices rapidly being deformed, while for small A z there is a long transient in which the flow remains quasi-2D (see the videos in [18]). Figure 1(a) shows the time average of E i k (k) for the simulations with 512 3 grid points, and for different values of A z . For large values of A z initial vertical gradients in the quantized vortices are large, and the spectrum peaks at k ≈ k 0 followed by a spectrum compatible with a direct energy cascade and with the emission of Kelvinwaves at wavenumbers smaller than the inverse mean intervortex distance [20]. The energy fluxes in Fig. 1(b), specially for A z = 1, are positive for all k and remain approximately constant for a range of wavenumbers k > k 0 . But for small values of A z initial vertical gradients are small, and the energy spectrum grows for k k 0 , developing a power law compatible with Kolmogorov scaling, and with negative total energy flux for k k 0 (albeit the negative flux does not remain constant with k, as a result of limited spatial resolution and of the inverse cascade being only transient in the absence of external forcing). In spite of this, the system develops a strong inverse transfer of energy, at least up to t = 10. For longer times the flow eventually becomes unstable and 3D. However, we verified that the 2D behavior is not simply due to an absence of 3D motions for A z 1. For t 10, when the systems display an inverse transfer of energy, the energy in the 3D modes for all A z = 0 is significant enough to nonlinearly act back to the 2D part of the flow and saturate its initial exponential growth, but not strong enough yet to suppress the inverse transfer. with a direct energy cascade is generated when γ = 1), but with different aspect ratios γ, using a spatial resolution N x = N y = 512, and with N z varied between 512 and 32 grid points to keep the vertical resolution ∆z fixed or over-resolved as L z is decreased, so that vertical gradients in the flow are always correctly resolved. Although the amplitude of the perturbation A z is fixed, by decreasing γ we also increase the wavenumber of the vertical perturbation (i.e., vertical variations of quantized vortices increase as the domain becomes thiner). As in the cubic domain, we observe an increase in E i k (k) for k k 0 and a range of wavenumbers with Π(k) < 0 but now for small values of γ, and a direct cascade of energy for large values of γ. But, unlike the case of the cubic domain, when γ is sufficiently small the flow remains quasi-2D for very long times, and quantized vortices aggregate in physical space creating larger structures (see movies in [18]).

Π < /Π tot Z = A z /A c z (N 3 = 256 3 ) Z = A z /A c z (N 3 = 512 3 ) Z = γ/γ c
To quantify the transition between the direct and inverse cascade regimes, we need to use (as an order parameter) an observable that measures the relative strength of the inverse energy cascade. To this end we first introduce the mean inverse flux at small wavenumbers (which is zero when the flux is positive), and the mean direct flux at large wavenumbers, respectively defined as

Π < = min 0, 1 k 0 k0 k=0 Π(k) , (1)
Π > = 1 k max -(k 0 + 1) kmax k=k0+1 Π(k), (2)
where k max = N x /3 is the maximum resolved wavenumber in the simulations, and k 0 is as before the wavenumber where the incompressible kinetic energy is initially concentrated. We can then compute the total energy flux (in both directions) Π tot = Π < + Π > , and define the normalized ratio of inverse energy flux to total en- ergy flux Π < /Π tot . Figure 2 shows the behavior of this quantity for all cases studied, as a function of the amplitude of the 3D perturbation normalized by its critical value A z /A c z (for spatial resolutions of 256 3 and 512 3 grid points), and as a function of the aspect ratio normalized by its critical value γ/γ c (for fixed A z ). In all cases we see an abrupt change as the control parameter is varied. For A z /A c z or γ/γ c > 1 there is no inverse energy flux, while for A z /A c z or γ/γ c < 1 it grows rapidly (albeit differently in each case). In the thin domains, from Fig. 1 it can be seen that γ c ≈ 0.1, corresponding to a domain with L z = L ⊥ /10 ≈ 11ξ (where ξ is the healing length of the condensate, proportional to the vortex core radius). This implies that the 2D behavior occurs when the height of the domain is of the same order as the correlation of the initial conditions, γ c 0 /L ⊥ , a similar condition for the layer height and the forcing lengthscale found for the compactified case in classical flows [13]. In other words, a transition towards 2D behavior does not require 2D domains or very slim films. Even moderate aspect ratios are enough to trigger an inverse energy cascade.

Energy fluxes, although they give direct indication of the presence of an inverse cascade, are not easily measurable in laboratory experiments. There are however other global quantities that are tractable experimentally and can also give an indication of a transition from 3D to 2D behavior in the flows as the control parameters are varied. Figure 3 shows the incompressible kinetic energy E inc in these flows as a function of time, both for 512 3 simulations in cubic boxes with different A z as well as for simulations in domains with different γ, and for each case, also the total length of the vortices as a function of time. In the simulations with large A z or γ, E inc decays in time after t ≈ 1, as the direct cascade of energy transfers the incompressible kinetic energy to smaller scales where it dissipates into phonons [20,21,24]. However, for small A z or γ, E inc remains constant in time or decays very slowly, indicating energy remains at large scales as in classical 2D turbulence. The same behavior is seen in the total vortex length [18, 21], which grows and reaches a maximum in the 3D regime as a result of vortex stretching (later decaying as a result of vortex reconnection), but which remains approximately constant for all times in the cases of small A z or γ, pointing to the absence of vortex stretching as expected in 2D flows. Vortex reconnection also plays an important role at early times for large A z or γ, to speed-up the three-dimensionalization of the flow, after which vortex stretching can become more efficient. Finally, it is also important to note that in the simulations in cubic domains the total length of the vortices remains approximately constant at early times in all cases, and that the time when vortex stretching starts increases as A z decreases. This is consistent with our previous observations: In the cubic box, for smaller values of A z the flow remains quasi-2D for longer times, and the observed transient inverse cascade delays the growth of 3D excitations in the flow.

Given the above we can consider various quantities that can indicate the presence of a sharp transition. Here we focus only on the thin layer case that is somehow closer to what is experimentally realizable. Figure 4 shows the time to reach the maximum vortex stretching t max , the maximum length of the vortices l max , and the flow integral scale L as a function of γ, which is obtained from the incompressible kinetic energy spectrum as L = 2π k -1 E i k (k)dk/ E i k (k)dk, and provides an estimation of the flow energy containing scale. When L ≈ L x (the domain size), the flow has self-organized at the largest available scale in the domain. These quantities display an abrupt change near the critical value γ c as γ is varied. The time t max is larger when γ < γ c , while the maximum vortex length is larger when γ > γ c . Both behaviors are to be expected when the flow is 3D and displays vortex stretching, or when the flow is 2D and as a result does not. Finally, the flow integral scale L becomes larger (and close to L x ) when γ < γ c . This indicates that the inverse transfer of energy leads to the concentration of kinetic energy at large scales, and implies the formation of large structures in the flow (e.g., resulting from spatial aggregation of vortices). In the simulations varying A z , we also verified that the overall shape of the quantities in the curves in Fig. 4 remain the same when changing the spatial resolution of the simulations, although the actual values (e.g., the time t max or the maximum vortex length l max ) depend on the resolution: at larger resolution the flow becomes more turbulent and vortex stretching is more efficient.

The numerical results show the existence of a transition between 2D and 3D behavior in quantum turbulence. This transition can be obtained by varying the dimensionality of the flow (in a 3D cubic domain), or by changing the aspect ratio of the domain and compactifying one spatial dimension. Below critical values of the controlling parameters the flows display an inverse transfer of energy which results in the growth of the incompressible kinetic energy at large-scales, and in the aggregation of quantized vortices. For the quasi-2D regimes the system suffers an interesting double condensation: the BEC, and the out-of-equilibrium inverse energy cascade which can result in a condensation of the kinetic energy at the largest available scales in the system [12]. This behavior is compatible with critical transitions reported in classical turbulence [13][14][15][16], and reminiscent of recent observations in experiments of gaseous BECs [7][8][9]. For the 3D cubic domain, the critical amplitude of the 3D perturbation is A c z ≈ 10 -2 (for the 512 3 simulations), while in the thin domains the critical aspect ratio is γ c ≈ 1/10. As our system is not forced, the inverse energy cascade can only develop as a transient (see, e.g., [23] for a discussion of the equivalent configuration in the classical case), a configuration which is comparable to experiments of gaseous BECs where the flow is let to freely decay after initially stirring it [1][2][3]. However, in experiments of gaseous BECs the condensate is trapped inside a potential, which we are not considering in our simulations to increase the scale separation between the domain size and the vortex radius as much as possible. The study of the effect of trapping potentials in these cascades is left for future work. In this work we study a system of weakly interacting bosons of mass m at zero-temperature that is described by the GPE

i ∂ψ ∂t = - 2 2m ∇ 2 ψ + g|ψ| 2 ψ, (3)
where ψ is the wave function of the condensate and g = c/( √ 2ρ 0 ξ) is proportional to the scattering length (with c the speed of sound, ρ 0 the mean mass density, and ξ the healing length); in terms of these variables m = /( √ 2cξ). In dimensionless units, all simulations have ρ 0 = 1, c = 2, and ξ such that the vortex cores are well resolved by the spatial resolution considered. This results in ξ = (40 √ 2) -1 in all simulations with N x = N y = 256 spatial grid points, and ξ = (80 √ 2) -1 in all simulations with N x = N y = 512 grid points so that ξk max = 1.5 where k max is the maximum resolved wavenumber.

The total energy E tot is a conserved magnitude in the GPE, and can be decomposed into

E tot = E k + E int + E q , (4)
where E k is the kinetic energy, E int is the internal energy, and E q is the quantum energy, which are defined respectively as

E k = ( √ ρu) 2 2 d 3 r, E int = gρ 2 2m 2 d 3 r, E q = 2 2m 2 (∇ √ ρ) 2 d 3 r, (5)
where ρ is the fluid density and u the fluid velocity, obtained from Madelung's transformation with ρ = |ψ| 2 and u = ∇φ/m. In this description, quantized vortices correspond to lines with ρ = 0, with quantum of circulation Γ = h/m. The Helmholtz decomposition √ ρu = ( √ ρu) i + ( √ ρu) c can be applied to the kinetic energy to further decompose it into incompressible E i k and compressible E c k kinetic energy components [20,21]. As these energies are quadratic, it is straightforward to construct power spectra from them as

E i,c k (k) = 1 2 ( √ ρu) i,c 2 k 2 dΩ k , E int (k) = g | ρ| 2 2m 2 k 2 dΩ k , E q (k) = 2 2m 2 (∇ √ ρ) 2 k 2 dΩ k , (6)
where the hat denotes the Fourier transform, and Ω k is the solid angle in Fourier space. The GPE was evolved in time using a fourth-order Runge-Kutta method, and a pseudospectral method to compute spatial derivatives and nonlinear terms [19]. Time steps were chosen to satisfy the CourantFriedrichsLewy condition, and resulted in ∆t = 10 -3 in dimensionles units in simulations with N x = N y = 256 grid points, and in ∆t = 5 × 10 -4 in simulations with N x = N y = 512 grid points. With these choices, total energy is conserved in all simulations up to the sixth significant digit at t = 10.

PREPARATION OF THE INITIAL CONDITIONS

An initial random two-dimensional (2D) flow with a three-dimensional (3D) perturbation is constructed using a Clebsch representation of the incompressible velocity field u = λ∇µ -∇φ, [21], where the Clebsch potentials are a superposition of modes

λ = 1 2k λ 2k λ ki=1 cos x k λ cos πk i 2k λ + y k λ sin πk i 2k λ + φ ki × 1 + A z cos 2πz L z + ϕ ki , (7)
µ = 1 2k µ 2kµ kj =1 cos x k µ cos πk j 2k µ + y k µ sin πk j 2k µ , (8)
where φ ki and ϕ ki are random phases, the brackets [ . ] indicate the integer part of the argument (to satisfy periodicity of each mode), and k z is the wavenumber of the perturbation in the z direction. The Clebsch potential φ is determined by the condition ∇ • u = 0. The parameters k λ and k µ control the initial correlation length of the field, and A z controls the amplitude of the 3D perturbation. The initial conditions are designed to generate a disordered flow with quantized vortices that have a dominant 2D component, a 3D perturbation (when A z = 0), and a correlation length at intermediate scales (or wavenumbers) such as both direct or inverse cascades can develop. As described in [20,21], with these potentials an associated initial wavefunction ψ(x, y, z) can be constructed as the product of wavefunctions ψ e (λ(x, y, z), µ(x, y, z)) where the ψ e have zeros at the zeros of the Clebsch potentials (and thus quantized vortices in the corresponding x, y, and z coordinates). To reduce the contribution of compressible modes and the initial emission of phonons, before solving the GPE these initial conditions are integrated to convergence using the advective real Ginzburg-Landau equation [20,21], which is the imaginary-time propagation of the GPE Galilean transformed to preserve the velocity field u. The final result is a wavefunction compatible with the flow u and with minimal sound emission, and which is used as the actual initial condition of the GPE. As examples of the resulting initial conditions, Fig. 5 shows slices in an xy plane of the mass density using 128 3 spatial grid points, for k λ = 4 and different values of k µ (and with A z = 0, and thus 2D). Points with zero density are defects that correspond to quantized vortices. As the number of excited modes increases (i.e., as k µ increases), more vortices are generated, and the average distance between them decreases. Figure 6 shows the spectrum of the incompressible kinetic energy for each of these initial conditions, and it can be seen that the maximum of the spectrum takes place at a wavenumber that increases with k µ (i.e., the initial correlation of the flow changes as this quantity is varied) in accordance with the relation k ∼ M k 0 k ξ with M the Mach number and k ξ ∼ 1/ξ, both of them fixed values, and k ∼ 1/ the wave number associated with the intervortex distance [21]. Similar results are obtained when k λ is varied. To consider an initial flow with scale separation such that both direct and inverse energy cascades can develop, all simulations in this study are done with k λ = 4 and k µ = 10, such that the initial energy peaks at k 0 ≈ 10.

The effect of varying A z is illustrated in Fig. 7, which shows slices in the xz plane of the mass density for initial conditions with fixed k λ and k µ , and with different values of A z . For A z = 0 vortices are parallel and straight in the z direction, and thus generate a purely 2D flow. As A z increases the vortices curve until in some cases they can even close on themselves forming rings, generating an initially 3D flow. Videos of these initial conditions in cubic boxes and in thin domains, as well as of their time evolution under the GPE, can be seen as supplemental material.

For the study of the transition in cubic boxes, we considered in simulations with N x × N y × N z = 512 3 grid points values of A z = 1, 0.4, 0.1, 5 × 10 -2 , 4 × 10 -2 , 3 × 10 -2 , 10 -2 , 8 × 10 -3 , 6 × 10 -3 , 10 -3 , 8 × 10 -4 , 7 × 10 -4 , 6 × 10 -4 , 3 × 10 -4 , 10 -5 , and 0 (even more values of A z , in the same range, were considered in the simulations with N x × N y × N z = 256 3 , for a total of 33 simulations at this resolution). In the thin domain case, the 3D perturbation was fixed at A z = 0.1, and the aspect ratio was varied to take values γ = 1, 1/2, 1/4, 1/8, 1/10, 1/12.5, 1/16, 1/20, 1/32, and 1/64. As a result, a total of 59 simulations with different parameters was considered for the analysis.

ENERGY FLUXES

Under the GPE, the dynamics of the system conserves the total energy

dE dt = 0, (9)
which as a result implies that a detailed balance equation can be written in spectral space as dE dt (k) = T (k), (10) where T (k) is the transfer function [12,15,23]. In other words, the change of energy at any given wavenumber must correspond to a transfer of this energy to or from this wavenumber to all other wavenumbers. By integrating this equation up to some wavenumber, an energy flux can be defined as

Π(k) = - k 0 T (k )dk = - d dt k 0 E(k )dk = - dE < (k) dt . (11)
Using the decomposition of the energy in Eq. ( 4) and the Helmholtz decomposition, this flux can be further decomposed as

Π(k) = Π i k (k) + Π c k (k) + Π int (k) + Π q (k), (12)
where each component of the flux corresponds to the different energy components. We verified that similar results are obtained when the total energy flux Π(k) is used to measure the direction of the cascades, and when the flux of incompressible kinetic energy Π i k (k) is considered instead.

ESTIMATION OF VORTEX LENGTH

In a similar fashion as with the energy, one can define an incompressible momentum power spectrum. The high wavenumber components of this spectrum can be approximated as the sum of the momenta of all the vortices present in the flow, counted individually. This provides an easy way to estimate the total line length of the vortices in the flow. The method is detailed in references [21,22].

FIG. 1 .

 1 FIG. 1. (a) Spectrum of the incompressible kinetic energy averaged between t = 1 and 10 for simulations in cubic domains (Nx = Ny = Nz = 512) and different values of Az. Kolmogorov power laws ∼ k -5/3 are indicated as a reference by solid black lines. The vertical dashed line indicates the inverse mean intervortex distance. Note the growth of energy and a ∼ k -5/3 scaling for k 10 when Az is small. (b) Total energy fluxes for the same simulations. For small Az the flux becomes negative for k 10, and the positive flux for k > 10 decreases. References for the values of Az in (a) and (b) are provided in the inset. (c) Same as in (a) for simulations in thin domains, for different values of γ. (d) Same as in (c) for simulations in thin domains; the inset gives the values of γ. In all panels, dashed curves highlight simulations for which movies are available in [18].

Figures 1 (

 1 Figures 1(c) and (d) show similar results for simulations with a fixed value of A z = 0.1 (such that a 3D flow

FIG. 2 .

 2 FIG. 2. Ratio of inverse energy flux to total energy flux as a function of the normalized control parameter Z (either Az or γ, normalized by their respective critical values A c z or γ c , see inset). For cubic domains two curves are shown, corresponding to spatial resolutions of 256 3 and 512 3 grid points.

1 FIG. 3 .

 13 FIG.3. Time evolution of the incompressible kinetic energy (a) in simulations in cubic domains with different perturbation amplitudes Az (512 3 runs), and (b) in simulations in domains with different aspect ratios γ. The insets show the total vortex length [18] as a function of time for each case. References are as in Fig.1; a few labels are provided as guidelines.

FIG. 4 .

 4 FIG. 4. (a) Time to reach the maximum length of the vortices as a function of γ/γ c . (b) Normalized maximum total length of the vortices, lmax/Lx, as a function of the control parameter. (c) Flow integral length scale normalized by the domain size, L/Lx, as a function of the same control parameter.

  NPM and PDM acknowledge financial support from grants UBACYT No. 20020170100508BA and PICT No. 2015-3530.

  SUPPLEMENTAL MATERIAL: ABRUPT TRANSITION BETWEEN THREE AND TWO-DIMENSIONAL QUANTUM TURBULENCE THE GROSS-PITAEVSKII EQUATION

FIG. 5 .

 5 FIG. 5. Slices of the mass density in an xy plane for initial conditions with different values of k λ and kµ and with Az = 0. Dark regions correspond to vortex cores. Note the random spatial distribution of vortices, and the change in their number and mean separation as kµ is varied (similar results are obtained when k λ is changed).

k = 4 ; k = 10 FIG. 6 .

 4106 FIG.6. Incompressible kinetic energy spectrum for the initial conditions in Fig.5. The wavenumber corresponding to the maximum of the kinetic energy spectral density changes with k λ and kµ.

FIG. 7 .

 7 FIG. 7. Slices of the mass density in an xz plane for initial conditions with different values of Az. Dark regions correspond to vortex cores. The configuration is independent of z for Az = 0, while the vortex cores become more deformed in the vertical direction as Az is increased.

MOVIES

The movies provided as supplemental material span the entire time evolution of the flow (from t = 0 to 10). The 3D renderings of quantized vortices in these movies provide examples of the behavior above and below the critical parameter A c z or γ c (or, in other words, 2D-like and 3D-like behavior), and correspond to the following cases: • Files side 0003.mp4 and side 4.mp4 are two examples of vortex evolution in the 3D cubic domain (at 512 3 resolution), respectively with A z = 0.0003 (A z < A c z ) and with A z = 0.4 (A z > A c z ). In the case with A z < A c z , note the system remains quasi-2D for a long time, until eventually 3D perturbations grow and dominate the dynamics.

• Files aniso 125.mp4 and aniso 03125.mp4 are two examples of vortex evolution in the thin domain, one with γ = 0.125 (γ > γ c ) and the other with γ = 0.03125 (γ < γ c ). In the case with γ < γ c the flow remains quasi-2D at all times, showing no vortex reconnection and spatial aggregation of quantized vortices.