Order Indices and Entanglement Production in Quantum Systems

The review is devoted to two important quantities characterizing many-body systems, order indices and the measure of entanglement production. Order indices describe the type of order distinguishing statistical systems. Contrary to the order parameters characterizing systems in the thermodynamic limit and describing long-range order, the order indices are applicable to finite systems and classify all types of orders, including long-range, mid-range, and short-range orders. The measure of entanglement production quantifies the amount of entanglement produced in a many-partite system by a quantum operation. Despite that the notions of order indices and entanglement production seem to be quite different, there is an intimate relation between them, which is emphasized in the review.


Introduction
Many-body systems can be described by different characteristics representing the main system properties. One of the pivotal notions is that of the type of order associated with the system state. The order is usually quantified by order parameters (see, e.g., [1][2][3]). The order parameters, as is well known, are rigorously defined in the thermodynamic limit, while for finite systems they, strictly speaking, degenerate to zero [3]. Here we mean a mathematically rigorous definition of order parameters, although in practical experiments and in numerical modeling one can still see the order parameters behave as if the thermodynamic limit is reached, when the system size is much larger than the correlation length.
At the same time, large many-body systems, even being finite, can posses a kind of order that is not of a long-range type, but rather it is a quasi-long-range or algebraic order. In addition, some two-dimensional systems exhibit quasi-long-range order below the Berezinskii-Kosterlitz-Thouless transition [4,5].
Order indices for density matrices were introduced in Refs. [6][7][8][9] as a general tool for quantifying all types of order, whether long-range, mid-range, or short-range. These characteristics can be defined for any statistical system, whether finite or in the thermodynamic limit. The notion of order indices can be generalized for arbitrary operators or matrices [10].
Another important notion characterizing the state of a multi-partite system is entanglement widely employed in quantum information processing and quantum computing [11][12][13][14][15][16][17], as well as in the theory of quantum measurements and quantum decision theory [18][19][20][21][22]. To be more precise, one has to distinguish between the entanglement of a state, or a statistical operator, or generally of an arbitrary operator, and the entanglement production describing the action of an operator on the given Hilbert space. The state entanglement characterizes the structure of a statistical operator, while the entanglement production by an operator describes the result of an operator action. These notions will be concretized below.
The aim of the present review is two-fold. First, a survey of the notions of order indices and of entanglement production, and their applications for treating many-body systems, will be given. Second, we shall emphasize the interrelation between these characteristics. It turns out that the latter are intimately related with each other, so that the qualitative change of a state leads to the quantitative changes of order indices as well as of the entanglement production.
Throughout the paper, the system of units is used where the Planck and Boltzmann constants are set to one.

Order Indices
Below, we consider systems composed of N parts, or particles, with N being finite, although it can be rather large. Suppose a collection {Â} of trace-class operators acts on a Hilbert space H, so that The order index of an operatorÂ is defined as where the base of the logarithm can be taken according to convenience, since the above definition does not depend on the choice of the base due to the property log a (x) = log b (x) log a (b). In other words, the order index is the exponent connecting the norm and trace of an operator, ||Â || = | Tr HÂ | ω(Â) .
Dealing with Hermitian operators, we get Comparing two operatorsÂ 1 andÂ 2 , we say that the operator having a larger order index is better ordered. In physical applications, the operator or matrix order indices characterize the type of order associated with the considered operators. Examples will be given below.

Entangled Structures
One has to distinguish two different notions, entangled structures and entangling operations. The first notion characterizes the property of such structures as wave functions, statistical operators (quantum states), and which can be generalized to arbitrary operators. The second notion describes the action of quantum operations on the members of a Hilbert space. To better explain the difference of these notions, in the present section we recall the main definitions concerning entangled structures, and in the next section we shall elucidate the meaning of the entanglement production by quantum operations.

Entangled Functions
The notion of wave-function entanglement was introduced by Schrödinger [23,24] with respect to quantum systems that can be separated into several subsystems. Then the system Hilbert space can be represented as the tensor product of the subsystem Hilbert spaces. Wave functions of the whole system, pertaining to the space H can be separated into two classes, separable and entangled functions. Separable functions have the form of the product Entangled functions can be represented as the linear combinations where at least two coefficients c α are not zero. From the physical point of view, separable wave functions are rather exceptional, being attributed to quantum subsystems that have never interacted and are distinguishable [25]. Generally, subsystems are characterized by entangled wave functions [26,27]. The collection of all separable functions forms a disentangled set All entangled functions compose an entangled set Considering wave functions, one keeps in mind the functions normalized to one, but in principle, the same definitions are valid for non-normalized functions that are the members of the Hilbert space H.

Entangled States
Generally, quantum states are represented by statistical operators. For pure states, statistical operators can be written asρ where the standard bra-ket notation is used. Separable wave functions compose separable stateŝ while entangled wave functions form entangled stateŝ More generally, a state is separable if and only if it has the structurê in which andρ iα is a statistical operator acting on a partial Hilbert space H i . A state that cannot be represented in the above form is entangled.

Entangled Operators
The notion of entangled states can be generalized to trace-class entangled operators [19][20][21]. Let us consider an algebra A = {Â} of trace-class operators acting on a Hilbert space H A . For any two operatorsÂ 1 andÂ 2 from the algebra A one can introduce the Hilbert-Schmidt scalar product The triple of the operator algebra A, acting on the Hilbert space H A , and the Hilbert-Schmidt scalar product form the Hilbert-Schmidt space Similarly, one can define the Hilbert-Schmidt space for another trace-class operator algebra B as The composite Hilbert-Schmidt space is given by the tensor product An operatorĈ sep from the composite Hilbert-Schmidt space in Equation (15) is separable if and only if it can be represented in the form Otherwise, the operator is entangled. In general, an operatorĈ sep , defined on the Hilbert space in Equation (6) is separable if and only if it can be represented asĈ where the operatorsÂ iα act on H i . These definitions lift the notion of entanglement to the operator level. In the Hilbert-Schmidt space, operators are isomorphic to the space members, so that the distinction between entangled and separable states in the Hilbert-Schmidt space becomes similar to that in the Hilbert space [19][20][21][28][29][30][31].

Entangling Operations
The other notion is the entanglement production by quantum operations. Considering the operators acting on a Hilbert space H, it is possible to distinguish two types of the operator actions, entangling and nonentangling. If an operatorÂ, acting on any function from the disentangled set D, leaves the function in this set, it is called a nonentangling operator [32,33], However, if there exists at least one function of the disentangled set D that becomes entangled under the action of the operator, this operator is termed entangling [34,35]. The strongest type of an entangling operator is a universal entangling operator that makes all disentangled functions entangled [36], A principal problem is how to measure the entangling power of operators. When one is interested in just a few wave functions, it is admissible to analyze the action of a given operator on all the functions of interest. In the case of a bipartite system, one can check the amount of the produced entanglement by studying the entanglement entropy for the considered few wave functions [28,30,31,37]. This, however, does not allow for quantifying the entangling power of the examined operator on the whole Hilbert space.
A general measure of entanglement production for arbitrary trace-class operators acting on a Hilbert space was advanced in Refs. [38,39]. The idea behind the definition of this measure is to compare the action of the given operatorÂ with the action of its nonentangling counterpart in whichÂ i are partially traced operatorsÂ The coefficient in Equation (20) is defined so that to preserve the trace normalization The measure of entanglement production [38,39] by an operatorÂ on a Hilbert space H is This measure is based on the comparison of the action of an operator on the whole Hilbert space with the action of its nonentangling counterpart that leaves invariant the disentangled set, It is useful to emphasize that, as has been proved [40][41][42][43][44][45], the only operators preserving separability are the operators having the form of tensor products of local operators and a swap operator permuting Hilbert subspaces in the tensor product of the total Hilbert space of a composite system. However, the action of the swap operator is trivial, merely permuting the indices labeling the subspaces. Up to the enumeration of subspaces, the product operators are the sole operators preserving the separability of functions. The tensor-product operators, as is evident, do not produce entanglement, Thus, the entanglement-production measure (23) is zero for nonentangling operators, and also it is continuous, additive, and invariant under local unitary operations [39,46]. As it should be for being a measure, it is semipositive. The sketch of the proof of this important property is as follows [39,46].
The set of trace-class operatorsÂ acting on the Hilbert space H, with a given operator norm, forms the Banach space which is a complete normed linear space. Similarly, the set of the product operatorsÂ ⊗ , leaving invariant the disentangled set D, composes the Banach space The latter space, by definition, is a subspace of the Banach space seen in Equation (24), Then it is admissible to define a projector transforming the members of space in Equation (24) into the members of space in Equation (25),P with the standard projector propertieŝ Therefore we have ||Â ⊗ || = ||P ⊗Â || ≤ ||P ⊗ || · ||Â || = ||Â || , from where the semi-positivity of the measure follows: It is important to stress that the entangled structure of an operator and entanglement production by this operator are quite different notions. An operator can be separable but entangling. Thus, the action of the separable operator in Equation (17) on the separable function in Equation (7) results in an entangled function,Ĉ It is possible to notice that the measure of entanglement production in Equation (23) and the order indices in Equation (2), although having rather different meanings, but are connected with each other through the relations The relations between these quantities will be considered in more detail below.

Density Matrices
To illustrate physical applications of the introduced notions, it is reasonable to consider such important physical quantities as reduced density matrices, which can be treated as matrix elements of reduced density operators [47]. For instance, the first-order density operator is expressed through the matrix elements Here x denotes a set of variables, such as spatial coordinates and spin, ψ(x) are field operators, the trace is over the Fock space generated by the field operators [48], andρ is a statistical operator. The second-order density operatorρ has the matrix elements Generally, the n-th order density operator is defined through the matrix elements The reduced density matrices of different orders are connected with each other by means of the relations in which N is the number of particles in the system. For example, The trace operation here implies the summation over the variable x, defined as In particular The non-entangling product operator for the second-order density matrix is expressed through the matrix elements Keeping in mind the operator norm in Equation (5), associated with the vector norm, we have where {ϕ k } is a natural ortho-normalized basis [47] and k is a multi-index labeling quantum states.
In first order, we have with the notation In the second order, we find with In this way, we can calculate the order indices of density matrices as well as the measure of entanglement production Using the properties of reduced density matrices, one can show [9,10,47,49] that the order indices, depending on statistics, satisfy the inequalities for Bose-Einstein statistics and for Fermi-Dirac statistics. When the upper boundary is reached, this signifies the occurrence of long-range order. Otherwise, there can only be mid-range or short-range order.

Correlation Matrices
Reduced density matrices are a particular case of correlation functions constructed by means of the field operators. In general, it is possible to consider correlation functions based on some other operators. Thus, taking an arbitrary operatorÂ(x), acting on the Fock space, one can set the correlation functions Then, one can introduce the correlation operator whose matrix elements are the above correlation functions. For the correlation operator, it is straightforward to define the order indices and the entanglement-production measure By choosing appropriate correlation functions one can quantify the properties of arbitrary physical systems.

Examples of Order Indices
Order indices are defined by the structure of the system state and can essentially vary under phase transformations [50]. Some examples of order indices that have been considered in literature are mentioned below.

Superconducting State
The structure of reduced density matrices have been analyzed in several works [6,8,9,47,49]. In the thermodynamic limit, the order indices of density matrices in a three-dimensional system are zero for the normal state and take the values for superconducting state. This corresponds to even long-range order.

Bose-Condensed System
Under the usual Bose-Einstein condensation into the state with zero momentum, in the thermodynamic limit in three dimensions, all order indices grow from zero to one, becoming This corresponds to the total long-range order [47,49]. Developing long-range order strongly influences the structure of correlation functions and density matrices [51][52][53][54].

Even Bose Condensate
There are theoretical speculations [55][56][57][58][59][60][61][62][63] that in Bose systems there can develop the so-called even condensate, with the formation of pairs similar to those for Fermi systems in the superconducting state, when thus exhibiting even long-range order.

Finite-Momentum Condensate
If the condensed state is characterized not by zero momentum but by a momentum with finite absolute value and random direction [51,[64][65][66], then the system possesses mid-range order with the order indices that, depending on dimensionality d, are For d > 1 this implies a mid-range order [10].

Two-Dimensional Systems
Two-dimensional systems below the Berezinskii-Kosterlitz-Thouless transition [4,5] temperature T K possess the order indices [6] where η is the exponent describing the behavior of the pair correlation function at large distance.
In that case so that there exists a mid-range order with the order indices Mid-range order also develops in two-dimensional antiferromagnets in a strong magnetic field [67,68].

Critical Point
At the point of a second-order phase transition, we have [6] with η being the exponent characterizing the pair correlation function. For different dimensions, this gives which implies mid-range order. In particular, for the two-dimensional Ising model, we have [69][70][71] For the three-dimensional Heisenberg model, η = 0.036, therefore

Tonks-Girardeau Gas
One-dimensional systems usually do not exhibit long-range order, although they can possess mid-range order. As an example, let us consider the Tonks-Girardeau gas [72]. This is a one-dimensional system of impenetrable bosons described by the Hamiltonian where U(x) is an external potential, and complimented by the condition on the wave function where a 0 is an effective diameter of particles. Studying the properties of reduced density matrices [73][74][75][76] one concludes that the order indices are for both homogeneous as well as for systems trapped in power-law potentials. This shows that the Tonks-Girardeau gas possesses a mid-range order.

Lieb-Liniger Model
This is a one-dimensional model described by the Hamiltonian [77] with point-like interactions of finite strength Φ 0 . The characteristic dimensionless interaction parameter is defined by the ratio of the typical potential energy to the typical kinetic energy The model is exactly solvable by Bethe ansatz [77,78], which provides the exact expression of the many-body eigenfunctions [79]. At zero temperature, the ground-state energy, the sound velocity c, and other equilibrium quantities can be expressed in terms of the solution of the Lieb-Liniger integral equations [77]. The other useful quantity is the Luttinger parameter [80] in which v F is the Fermi velocity and c is the sound velocity. Note that this parameter reminds the effective Mach number used for characterizing statistical systems [81]. Solving the Lieb-Liniger integral equations one has access to the sound velocity c and the Luttinger parameter M for any values of the coupling constant [82][83][84][85].
The order index of the first-order density matrix can be estimated [86] as For small γ, one has c ρ m While for strong interactions, when the order index tends to that of the Tonks-Girardeau gas, As we see, the Lieb-Liniger model demonstrates a mid-range order. Formally, the order index tends to one in the limit of the ideal gas, when γ = 0. However, the ideal one-dimensional uniform gas is unstable [87,88].

Finite Systems
Strictly speaking, usually long-range order can arise only in thermodynamic limit. In finite systems, order parameters not always can be rigorously defined. However, the order indices can be well defined for any finite system. It is reasonable to expect that in finite systems, instead of long-range order, there could exist only mid-range or short-range order. In that case, the order indices vary in the interval depending on the number of particles in the system, on the interaction strength, and on other system parameters [89,90]. We shall discuss the details of calculating the order indices for systems with a finite, although large, number of particles in Sections 9 and 10.

Examples of Entanglement Production
Now we shall give several examples of the entanglement-production measure calculated [39,91] according to Section 4. For bipartite states, it is also possible to find the entanglement entropy is a partially traced statistical operator. The latter can be compared with the entanglement-production measure in Equation (23). In general, measures in Equations (23) and (69) do not need to be equal.
In addition, the entanglement entropy in Equation (69) can be defined only for bipartite states, while the entanglement-production measure in Equation (23) can be defined for arbitrary states and operators.

Einstein-Podolsky-Rosen-States
The statistical operator of this entangled pure state iŝ This is a two-particle two-mode state. For the entanglement-production measure, we have This coincides with the entanglement entropy in Equation (69).

Bell States
This is also an entangled two-particle two-mode pure state with the statistical operator The measure in Equation (23) ε(ρ B ) = log 2 equals the entanglement entropy in Equation (69).

Greenberger-Horne-Zeilinger States
This is an N-particle two-mode statê For the measure in Equation (23), we get For the number of particles more than two, the entanglement entropy is not defined.

Multicat States
The statistical operator for the two-mode N-particle statê is expressed through the multicat function where |c 1 | 2 + |c 2 | 2 = 1 .
The entanglement-production measure becomes Multicat states can be realized in multiparticle systems where each particle possesses two internal energy states, for instance, for trapped ions subject to the action of resonant laser beams [92,93] or for Bose-condensed neutral atoms under coherent Raman scattering [94]. Instead of internal single-particle states, one can use collective nonlinear states created by means of the resonant excitation of topological coherent modes in trapped Bose-Einstein condensates [95][96][97][98].

Multimode States
The state of N particles possessing N M modes each, is formed by the multimode function This is a generalization of the multicat state to the case of N particles. The measure of entanglement production is The multimode functions can be represented by coherent states [99].

Hatree-Fock States
The state with the statistical operator ρ HF = | HF HF | is formed by N-particle symmetrized or antisymmetrized functions Such functions represent N identical particles, bosons or fermions, in N different energy states. The entanglement-production measure reads as In the case of two particles, the measure reduces to the that for the Einstein-Podolsky-Rosen state, ε(ρ HF ) = log 2 (N = 2) .
For macroscopic systems, it tends to the expression ε(ρ HF ) N log e (N 1) .

Reduced Statistical Operators
In the above examples, pure states are considered. The measure of entanglement production, being general, can be found for mixed states as well. To this end, we can study reduced statistical operators defined through partial traces of the Hartree-Fock statê ρ n ≡ Tr H n+1 Tr H n+2 . . . Tr H Nρ HF , with n = 1, 2, . . . , N − 1. The entanglement-production measure for this operator is For a large number of particles, this yields

Gibbs States
An equilibrium system of N particles, characterized by a Hamiltonian H acting on a Hilbert space is described by the Gibbs statistical operator in which β ≡ 1/T is inverse temperature. The entanglement-production measure can be calculated as is explained in Section 4. For two-particle registers, represented by Ising and Heisenberg Hamiltonians, this was done in Refs. [39,91].

Werner Operator
The case of the Werner operator [100] is interesting since, depending on parametrization, it can represent a separable or entangled state. The operator has the form with the unity operator1 The name of the latter comes from the property of the flip operator to flip functions: The parametrization of the Werner operator is done by the parameter γ defined by the equation γ = Tr HρWσ (Tr HρW = 1) .
The partially traced operators arê These reduced operators are normalized to one, The nonentangling counterpart of the Werner operator iŝ which gives its norm In view of the matrix elements Thence the norm of the Werner operator is Therefore for the entanglement-production measure we find The measure is positive for all values of γ, which means that the Werner operator is always entangling. At the same time, the positive partial transpose criterion [101,102] tells us that the Werner state is separable if and only if γ ≥ 0, while it is entangled for γ < 0. However, even being separable, the Werner operator is entangling, that is, producing entanglement. The Werner operator is an explicit example of an operator that can be separable, although being entangling, as is discussed in Section 4. The entanglement entropy for the Werner operator is

Bose-Einstein Condensation
The aim of the present section is twofold. First, it shows explicitly how the order indices and the entanglement-production measure can be calculated for a system with a phase transformation for a large, however finite, number of particles. Second, it is demonstrated that both these characteristics, order indices as well as entanglement production, simultaneously change under phase transitions.
The correct description of a Bose-condensed system requires the validity of the following important stipulations. (i) The theory has to respect conservation laws and thermodynamic relations.
(ii) The excitation spectrum must be gapless. (iii) The Bose-Einstein condensation has to be a phase transition of second order. (iv) The system must be stable, satisfying the necessary stability conditions. (v) Reasonable quantitative agreement with experiments is required. The validity of these stipulations can be achieved in the self-consistent approach based on representative statistical ensembles [103][104][105][106] which we use here.
In the presence of Bose-Einstein condensate, the boson field operators acquire the Bogolubov shift [107,108], being represented as the sum in which η is the condensate wave function and ψ 1 is the field operator of uncondensed particles. The field operators are also functions of time, which are not shown explicitly for the sake of notation compactness. The functional variables η and ψ are orthogonal with each other, and satisfy the conditions The condensate function is normalized to the number of condensed particles, while the number of uncondensed particles is The first-order density matrix reads as ρ(r, r ) = η * (r )η(r) + ψ † 1 (r )ψ 1 (r) .
Its eigenvalues can be represented through the integral in which ϕ k are the natural orbitals [47]. Taking account of Equation (87) leads to the sum where Thus we obtain the norm In the absence of the condensate, η = 0 and N 0k = 0.
If the system is uniform, then and it follows Considering the second-order density matrix, we shall use the simplified notation writing just j instead of r j and j instead of r j . For instance, Then the second-order density matrix reads as In what follows, let us consider a uniform system, when Equation (92) holds true. In that case, the natural orbitals are plane waves. By employing the Hartree-Fock-Bogolubov decoupling reduces the second-order density matrix to the expression The matrix element in Equation (94) becomes where Similarly, the anomalous average in Equation (91) takes the form with For the norm of the second-order density operator, we get In that way, we find the order indices for the density operators of first order, in which N k = N 0 δ k0 + n k , and of second order where For the entanglement-production measure, we obtain These quantities can be calculated [89,90] as functions of the number of particles and other system parameters.
For illustration, let us consider the limiting case of N → ∞. Then, for temperatures below the condensation temperature T c , we find that while sup k n k ∝ N 1/3 . At temperatures above the condensation temperature T c , we have Therefore the studied order indices are (103) and the entanglement-production measure is At the phase transition point, both these quantities experience noticeable change. The appearance of order is accompanied by the reduction of the entanglement production.

Magnetic Transitions
Magnetic systems are characterized by spin operators S j located at lattice sites enumerated by j = 1, 2, . . . , N. Instead of density matrices defined through field operators, we need now to introduce correlation matrices, as in Section 6, composed of spin operators [10]. Thus the first-order spin correlation matrixĈ is composed of the matrix elements The second-order spin correlation matrixĈ is formed by the correlation functions Keeping in mind that S 2 j = S(S + 1) , we have the traces and In calculating the matrix norms, we use the natural lattice orbital Then the norm of the first-order correlation matrix is Employing the properties and separating the terms with coinciding and non-coinciding lattice sites, we get Similarly, for the second-order spin correlation matrix, we have the norm which leads to Defining the magnetization for different lattice sites, one can resort to the mean-field approximation Then we have Depending on the number of coinciding and non-coinciding lattice sites, we find This leads to the matrix elements Therefore we obtain The nonentangling counterpart of the second-order correlation matrix, for which results in the norm The norms are essentially different for the magnetically ordered state, when a nonzero magnetization is present, and for the paramagnetic state without average magnetization: For the order indices we have For the entanglement-production measure we find Here a large system, with N 1 is assumed. Again we see that the entanglement production diminishes upon arising order and increases for a disordered state.

Diagonal Order
In the previous Sections 10 and 11, the cases of phase transformations with the arising off-diagonal order are analyzed. There exists an opinion that the transitions with the arising diagonal order have to be treated differently, since the related reduced density matrices behave in a different way. However, in the approach based on the order indices, there is no difference in the method of treating any type of phase transition. What needed is to define the appropriate correlation matrix [10]. In the present section, we illustrate this for the solid-liquid transition that is the most known transition exhibiting diagonal order.
Under the solidification-melting phase transition, what changes is the particle density that is the statistical average of the density operator In the liquid state, the particle density is constant in space, being equal to the average density In the solid state, the density is nonuniform, having minima and maxima. We can label the points of the density maxima by a j , enumerating them by the index j = 1, 2, . . . , N L , thus defining them by the condition max r ρ(r) = ρ(a j ) .
Generally, the points of maxima do not need to form a periodic structure, but can be randomly located, as in amorphous solids. For crystalline structures, the points of density maxima form a periodic crystalline lattice, where the set of {a j } fixes lattice sites. Below, we shall consider a crystalline lattice, where the points of density maxima are identical to each other, so that with a being any site from the set of lattice sites. This is not principal for the approach, but just simplifies some notations. Let us introduce the operatorÂ (a j ) ≡ρ(a j ) − ρ .
With these operators, it is straightforward to define the correlation functions These functions play the role of matrix elements of the correlation matrix for which The trace of this matrix is with N L being the number of lattice sites. This trace can be rewritten as where varρ(a) ≡ ρ(a) 2 − ρ(a) 2 is the density variance at a lattice site. For different sites, one can use the Hartree decoupling that is known to provide a good description for crystals [109][110][111][112][113][114]. Keeping in mind the lattice function in Equation (111), for the norm of the matrix, we find In this way, we come to the order index For a solid, this reduces to while for a liquid, it gives These expressions can be used for both, finite systems as well as for macroscopic systems with a large number of particles. Note that finite crystalline structures can be stable, or metastable, even in low dimensions, forming crystalline chains and planes [115] for which the Lindemann stability criterion [116] is valid. In the thermodynamic limit, we obtain Similarly, the characteristics of higher-order correlation matricesD n can be calculated. The overall situation is analogous to the cases of phase transitions treated in Sections 10 and 11.

Dynamical Effects
Order indices and entanglement production can vary with time. Time dependence can come through equations of motion. If the operator obeys a unitary evolution, the order indices do not change with time, although the entanglement production can vary since the non-entangling operator, generally, does not follow the unitary evolution. Reduced density matrices and correlation matrices, in general, also are not governed by unitary time dependence. Therefore the order indices and entanglement-production measure, say for reduced density matrices of nonequilibrium systems, do depend on time,

Multitrap Multimode States
As an example, let us consider the case of multiple traps filled by Bose-Einstein condensate. This can be realized with an optical lattice having deep wells, where condensate clouds are located. Thus, cold rubidium 87 Rb atoms were loaded [117,118] into an optical lattice, with adjacent sites spaced so that these sites were practically independent, with the tunneling time between sites above 10 18 s. The number of lattice sites was typically between 5 to 35. The number of condensed atoms in each site could be varied between about 200 to 10 4 . Shaking the lattice as a whole, it is possible to create in each lattice site multimode states of excited Bose-Einstein condensates [95][96][97][98][119][120][121][122][123]. Such states are described by the statistical operator where n k is the fraction of atoms in the k-th mode, N M is the number of modes in a lattice site, whose number is N L . The normalization of the statistical operator to one, requires the validity of the summation The entanglement-production is quantified by the measure varying in the range The temporal evolution of the measure was studied for two and three modes [124][125][126][127].
For the realization of the multitrap multimode states it is important that the traps be identical, so that the transition frequencies between the ground-state condensate and the generated coherent mode be the same in all traps. Only then it is feasible to shake the lattice with a frequency that would be in resonance with the transition frequency of all traps, which is necessary for the generation of the same mode in these traps. A collection of different traps, such as effective potential wells in a random matter, where the so-called Bose glass can arise [128][129][130] is not appropriate, since these effective wells are not identical, hence are characterised by different transition frequencies.

Entangling by Evolution Operators
Entanglement production by the evolution operator was studied in Refs. [46,131]. The operator norm was defined as the Hilbert-Schmidt norm. Calculations were performed for finite-site Heisenberg and Ising Hamiltonians. Depending on the system parameters, the evolution of the entanglement-production measure is periodic or quasiperiodic.

Coherence Phenomena
Coherence phenomena occurring in nonequilibrium systems can be treated as a kind of dynamic phase transitions. Therefore the temporal behavior of these coherent phenomena can also be accompanied by drastic changes in the order indices and entanglement production. Examples of these phenomena are given by the dynamics of strongly nonequilibrium spin systems [132,133] and radiating systems [134][135][136]. Being many-particle collections, these systems allow for considering order indices and entanglement production characterizing their temporal behavior. To describe nonequilibrium processes in these systems, it is convenient to resort to quasispin representation [134][135][136] based on the quasispin operators where the index j enumerates particles, or spins, and σ α j are Pauli matrices. The main observable quantities are the fractional population imbalance (or longitudinal spin polarization) and the dimensionless coherence intensity or radiation intensity that are functions of time coming from the equations of motion. Summation is over all atoms or spins. It is useful to mention that in nonequilibrium processes, such as radiation, the population imbalance s is never identically equals minus one, which would mean an equilibrium state, when the coherence intensity w would be zero. Since there are two physically different operators, one related to coherence intensity and the other describing the atomic population imbalance (spin polarization), it is convenient to consider two different types of correlation functions [137]. Dealing with the population-imbalance operators requires to define correlation functions using the z-spin operators. Thus, it is possible to introduce the first-order correlation operatorQ with the matrix elements The second-order correlation operator isQ having the matrix elements Similarly, it is straightforward to define the correlation operators composed of the ladder spin operators, the first orderR with the elements and the second orderR with the matrix elements To calculate the traces and norms of the correlation operators, we shall need the following equalities for coinciding locations, Let us start with the operatorsQ n . For the traces, we have Calculating the order indices and entanglement-production measure, we simplify the final expressions by keeping in mind a large number of atoms N 1. Then we get the order indices ω(Q 1 ) = ω(Q 2 ) = 1, s = 0 0, s = 0 (164) and the entanglement-production measure ε(Q 2 ) = 0, s = 0 log 2, The behavior of the order indices here is similar to the case of magnetic systems. However the physical nature of the arising order is rather different from the latter. In the process of coherent dynamics or radiation, the ordering occurs for a short time, since s = s(t) is a function of time. Here it is an example of a dynamical order, while in a magnet, it is a stationary order. Similarly to the stationary magnetic order, the dynamic order also reduces the measure of entanglement production. For the correlation operators related to the coherence intensity, we have the traces The partially traced operatorsR and the entanglement-production measure ε(R 2 ) = 0, w = 0 log 2, These quantities essentially depend on the coherence intensity w = w(t) that plays here the role of a dynamic order characteristic. A nonzero w signifies the appearance of a coherent dynamic order. Sometimes, the occurrence of the coherent dynamic order is paralleled with the nonequilibrium magnon condensation [138][139][140]. Note that in equilibrium systems magnons cannot condense [141,142].

Conclusions
We have surveyed the meaning and applications of two concepts, order indices and entanglement production. The concept of order indices generalizes that of order parameters. The latter usually describe long-range order and are well defined only for macroscopic systems in the thermodynamic limit. However the order indices can be defined for any system, whether in the thermodynamic limit or finite. They characterize all types of order, long-range, mid-range, or short-range, making it possible not merely distinguishing the qualitative types of orders, but prescribing a measure uniquely quantifying the level of ordering and being applicable to equilibrium as well as to nonequilibrium systems.
The entanglement production describes how much entanglement is produced by an operator. This should be distinguished from the notion of entanglement describing the operator structure. It is possible to say that entanglement, characterizing the operator structure, is a static notion, while the entanglement production, describing the operator action, is a kind of a dynamic notion. The entanglement production can be quantified by a measure that is valid for any system, bipartite or multipartite, equilibrium or not.
It turns out that order indices and the entanglement-production measure are closely related with each other. As a rule, the larger the level of order in a system, the smaller the entanglement-production measure. The use of these concepts helps to better understand the properties of the studied systems and to more efficiently employ them in applications.
In the review, we considered the cases that could be treated analytically. Dealing with finite systems, it is usually necessary to resort to numerical calculations. Thus powerful numerical methods have been developed for studying the dynamics of finite Bose systems [143][144][145]. Hopefully, the described notions could be successfully employed for the analysis of finite systems by applying numerical methods.
Funding: This research received no external funding.