Dynamics of Local Symmetry Correlators for Interacting Many-Particle Systems

Dynamics of Local Symmetry Correlators for Interacting Many-Particle Systems

P. Schmelcher Zentrum für Optische Quantentechnologien, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany Hamburg Centre for Ultrafast Imaging, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany    S. Krönke Zentrum für Optische Quantentechnologien, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany Hamburg Centre for Ultrafast Imaging, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany    F. K. Diakonos Department of Physics, University of Athens, GR-15771 Athens, Greece
July 15, 2019
Abstract

Recently (PRL 113, 050403 (2014)) the concept of local symmetries in one-dimensional stationary wave propagation has been shown to lead to a class of invariant two-point currents that allow to generalize the parity and Bloch theorem. In the present work we establish the theoretical framework of local symmetries for higher-dimensional interacting many-body systems. Based on the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy we derive the equations of motion of local symmetry correlators which are off-diagonal elements of the reduced one-body density matrix at symmetry related positions. The natural orbital representation yields equations of motion for the convex sum of the local symmetry correlators of the natural orbitals as well as for the local symmetry correlators of the individual orbitals themselves. An alternative integral representation with a unique interpretation is provided. We discuss special cases, such as the bosonic and fermionic mean field theory and show in particular that the invariant two-point currents are recovered in the case of the non-interacting one-dimensional stationary wave propagation. Finally we derive the equations of motion for anomalous local symmetry correlators which indicate the breaking of a global into a local symmetry in the stationary non-interacting case.

I Introduction

Symmetries represent a cornerstone of modern quantum physics including many different fields, such as atomic, molecular or condensed matter physics Zare (); Hammermesh (); Harris (); Dresselhaus (). As a consequence the corresponding mathematical framework is very well developed and provides many powerful tools in order to obtain statements on the properties of specific systems even without solving the underlying equations of motion (EOM). Examples are the discrete and continuous groups constituted by inversion, reflection, translation or rotation transforms which leave the underlying Hamiltonian invariant. Theorems deduced from the existence of these global symmetries, such as the parity or Bloch theorem, are a key ingredient for all further analysis. The presence of symmetries has immediate consequences for many properties of quantum systems, such as the existence of conserved quantum numbers or resulting selection rules for e.g. electromagnetic transitions Hammermesh (); Harris (); Dresselhaus ().

The situation becomes much less transparent and accessible if a system does not exhibit a global symmetry which holds in complete space but possesses different symmetries in different domains of space. This occurs typically for materials of a higher degree of complexity ranging from self-assembled configurations of atoms on surfaces Rahsepar () to frustrated crystallization following a liquid-glass transition Shintani (). Local order and local symmetries are ubiquitous also in larger molecules or even disordered matter Wochner2009 (). One way to detect local order on the atomic scale uses X-ray cross correlation analysis Wochner2009 (); Altarelli (); Wochner2011 (); Lehmkuehler (). In view of the existing powerful methodology to tackle global symmetries the question arises whether the presence of local symmetries also allows for the development of a predictive formalism and corresponding mathematical framework. This is a question of major importance not only for principal and fundamental reasons but in particular since materials with local symmetries provide a bridge between ideal crystals and disordered matter and it is to be expected that they will lead to a rich and unique phenomenology.

Recently several steps in the direction of a theory of local symmetries have been performed indeed Kalozoumis2013a (); Kalozoumis2014a (); Kalozoumis2014b (); Kalozoumis2013b (); Kalozoumis2015a (); Morfonios2014 (); Kalozoumis2015b (); Zambetakis2016 (); Kalozoumis2016a (); Wulf2016 (). The focus was hereby on the case of non-interacting stationary one-dimensional single particle or wave mechanical (acoustics, optics) problems which are described by the corresponding Schrödinger or Helmholtz equation where the prime denotes differentiation with respect to the spatial variable . Here is a real function generated by an effective wave vector , which describes the inhomogeneity of the medium where the wave propagates. For a matter wave , is the scaled kinetic energy of a particle with mass and energy which moves in a potential . The complexity arises due to the presence of the potential term, which can exhibit a plethora of different local symmetries. Here, each local symmetry is associated with a local coordinate transform , , which leaves the potential invariant: . In what follows, we will focus on , i.e. on local discrete inversion (parity) symmetries () and local translation symmetries (). An illustration of such a composite potential landscape is provided in Figure 1. As it was recently found Kalozoumis2013b (); Kalozoumis2014a () local symmetries lead to the existence of invariant two-point correlators (ITPC) which possess the dimensionality of currents and are constant in the corresponding domains of local symmetry. Two types of these ITPC were explored Kalozoumis2014a () and read as follows

 Q=12i[σΨ(x)Ψ′(y)−Ψ(y)Ψ′(x)]=const.∀x∈D (1)
 ˜Q=12i[σΨ∗(x)Ψ′(y)−Ψ(y)Ψ′∗(x)]=const.∀x∈D (2)

where lies in the co-domain of the local symmetry (see figure 1). is the wave function i.e. a solution of the stationary Schrödinger equation and the star denotes complex conjugation. The above quantities are defined for each domain of local symmetry separately and when considering complete real space and a complete local symmetry (see Figure 1(d)) they are piecewise constant. Looking back at the special case of a global symmetry which is respected also by the corresponding boundary conditions, one can show Kalozoumis2014a () that vanishes and represents the eigenvalue (parity or Bloch phase) under the symmetry operation applied to the wave field, where is the well-known globally conserved probability current

 J=12i[Ψ∗(x)Ψ′(x)−Ψ(x)Ψ′∗(x)] (3)

Therefore, finite values of indicate the symmetry breaking if we convert a global to a local symmetry, valid only in a limited spatial domain. This leads Kalozoumis2014a () to the generalization of the parity and Bloch theorem in the presence of local symmetries

 Ψ(y)=Ψ(F(x))=1J[˜QΨ(x)−QΨ∗(x)]∀x∈D (4)

The above relation (4) provides the mapping of the wave field between a spatial point and its symmetry transform which is solely determined by the ITPC . We remark that the above considerations can be extended beyond the linear transforms of inversion and translation symmetry Zambetakis2016 ().

The ITPC are the key quantities characterizing local discrete symmetries in one spatial dimension as has been shown in several works Kalozoumis2014b (); Kalozoumis2013b (); Kalozoumis2015a (); Morfonios2014 (); Kalozoumis2015b (); Zambetakis2016 (); Kalozoumis2016a (); Wulf2016 (). Applications to photonic multilayers have demonstrated that perfectly transmitting resonances can be completely classified according to sum rules imposed on the ITPC Kalozoumis2013b (). Even more, a construction principle has been derived which allows for a systematic design of locally symmetric multilayered structures that exhibit perfectly transmitting resonances at desired energies Kalozoumis2013b (). Very recently a local basis approach has been developed Zambetakis2016 () which provides a classification of the solutions of the wave equation in terms of ITPC and consequently a computational strategy for the determination of the wave field in complex locally symmetric potential landscapes (see Figure 1). Further very recent extensions of the theory of ITPC include time-dependent periodically driven Wulf2016 () and parity-time symmetric setups with gain and loss Kalozoumis2014b (); Kalozoumis2016a () where e.g. the ITPC provide an order parameter for the phase diagram of the scattering states Kalozoumis2014b (). The first experimental detection of the ITPC has been achieved in acoustic waveguides Kalozoumis2015b () where both the phase and magnitude of the pressure field yield a complete reconstruction of the complex ITPC. Here the theory has been extended to include losses and the experiment has also verified the occurence of perfectly transmitting resonances Kalozoumis2015b ().

While the above-mentioned theoretical framework and applications focus on one-dimensional stationary and non-interacting wave mechanical systems as ideally realized in acoustics and optics we develop here the theoretical framework for local symmetries (local inversion and local translation) in interacting many-body systems including higher dimensions and dynamics. The key idea is to define the (canonical) local symmetry correlators of an interacting many-body system as coherences of the reduced density matrices (RDM) between local symmetry related points in space. Starting from the well-known Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy of equations of motion for the RDM Born (); Bogoliubov (); Kirkwood (); Yvon (); Bonitz (); Akbari () and focusing on the dynamics of the one-body RDM local symmetry correlators, we arrive at a continuity-like equation for the latter, involving the divergence of non-local two-point current densities and a local source term stemming from the interaction with the surrounding particles. These equations become more transparent if one employs the diagonal natural orbital representation of the one-body RDM such that the local symmetry correlators and their corresponding non-local current densities are represented as convex sums over local symmetry correlators of the individual natural orbitals and their respective non-local current densities. Corresponding, continuity-like equations of motion for the local symmetry correlators of the individual natural orbitals are derived. Moreover, an integral representation for the equations of motion for the one-body local symmetry correlators is presented, which allows for a novel interpretation of the interplay between non-local currents and local symmetry correlators. Furthermore, we discuss the impact of the range of the interaction potential on the local source term of the derived continuity-like equations of motion. This provides us with the fundamental framework of local symmetries of interacting many-body systems be it fermions or bosons. Several special cases are subsequently discussed and it is, in particular, shown that the previously developed non-interacting theory in one spatial dimension is indeed contained in the current general approach. Besides the above we also establish the EOM for the anomalous symmetry correlators specializing to the ITPC for the non-interacting case.

This work is structured as follows. Section II contains the theoretical framework. It introduces our Hamiltonian of the interacting particle system and some important properties of the resulting density matrix. The BBGKY hierarchy, being the starting-point of our theoretical investigation, is subsequently discussed and we review the EOM for the natural orbitals and their populations. In section III the EOM for the symmetry correlators are derived and discussed in particular in the presence of local symmetries. Different representations are shown and special cases, such as the non-interacting situation as well as mean field theories of bosons and fermions are addressed. The EOM for the anomalous symmetry correlator reducing to the known symmetry breaking two-point current in the stationary one-dimensional and non-interacting case are established as well, and the uniqueness of this construction is critically examined. Finally, in section IV we conclude and provide some outlook.

Ii Theoretical Framework

As described above significant progress has been achieved in terms of both the formal description of local symmetries and their impact and applications. Still, the focus so far is on a one-dimensional stationary non-interacting theory, which nicely applies to acoustics and optics and non-interacting i.e. single particle quantum mechanics. Therefore the obvious question arises whether the theoretical concept of local symmetries can be extended to the time-dependent situation, to higher dimensions and interacting many-body systems. According to observation, local symmetries and local order do in many cases occur (in nature) when many interacting degrees of freedom cooperate and a high degree of complexity or nonequilibrium is met. In view of this the pressing question arises what an adequate level of description for the local symmetries of interacting systems would be. Trying to straightforwardly generalize the above-mentioned approach of the non-interacting case but now starting with the -particle Schrödinger equation including external (locally symmetric) potentials and two-body interactions and following analogous steps in the derivation leads to a conflicting situation. While in the non-interacting case a particle feels only the external local potential and it therefore definitely belongs to a certain domain of a given local symmetry, the case of interacting particles leads to the situation that two interacting particles are in general located in two different local symmetry domains. As a consequence, the reduction and derivation of expressions similar to eq.(1,2,4) are impossible. Even if one assumes sufficiently short-ranged (or contact) interactions, such that ideally two-interacting particles always belong to the same domain of local symmetry, the resulting condition represents a -dimensional divergence of a -dimensional generalization of the ITPC in eqs.(1,2) in spatial dimensions. This condition cannot be easily used to advance, and, in particular, it cannot be straightforwardly integrated as it is the case for the non-interacting one-dimensional theory.

In view of the above an approach based on integrated effective degrees of freedom is desirable. Such an approach is provided by the corresponding density matrices and their EOM, the BBGKY hierarchy Born (); Bogoliubov (); Kirkwood (); Yvon (); Bonitz (); Akbari (). In the following we will provide our setup and Hamiltonian, as well as important properties of the RDM and their equations motion. This will be the basis for our derivation of the EOM of the symmetry correlators in the presence of local symmetries.

ii.1 Hamiltonian and Reduced Density Matrices for Interacting Particle Systems

Our starting-point is a system of identical particles governed by the Hamiltonian

 H=N∑i=1(−12Δi+U(xi))+12∑i≠jV(xi,xj) (5)

where we have assumed atomic units (a.u.) and are space-spin coordinates. is the single particle potential which we will, later on, assume to possess local symmetries. represents the interaction between the -th and -th particle. Solutions to the time-dependent Schrödinger equation are then given by the time-dependent wave packets whose th order (body) RDM are given by

 ρ(n)(Xn;X′n;t)=∫Ψ∗(X′n,Xcn,t)Ψ(Xn,Xcn,t)dXcn (6)

with the compact notation and where indicates the complement w.r.t. the complete vector . Here and the integration w.r.t. includes a spatial integration and a summation over the spin degrees of freedom. Let us mention some relevant properties of the RDM. Of course, RDM inherit the bosonic or fermionic exchange symmetry within each of their two sets of coordinates from the underlying complete wave functions. They obey a partial trace relation meaning that lower order RDM can be directly calculated from a given order RDM by setting equal and integrating out the corresponding coordinates. RDM are hermitian and possess (in our probabilistic state normalization) trace one. RDM are positive semidefinite yielding that all of their eigenvalues are equal to or greater than zero. The corresponding eigenvalue relation reads

 ∫ρ(n)(Xn;X′n;t)ϕ(n)i(X′n,t)dX′n=(N−n)!N!λ(n)i(t)ϕ(n)i(Xn,t) (7)

The eigenvectors and eigenvalues of are refered to as the natural orbitals and their occupation numbers or natural populations. For bosonic systems, we have whereas for fermions . Specifically the diagonal representation of reads

 ρ(1)(x1;x′1;t)=1N∑iλi(t)ϕi(x1,t)ϕ∗i(x′1,t) (8)

where, for simplicity, we have omitted the upper index indicating the order for the natural orbital and their occupation numbers .

ii.2 The BBGKY Hierarchy and the Equations of Motion for Natural Orbitals

Since the Schrödinger equation possesses a unique solution for a given initial wave function this holds also for the EOM of the RDM. Working with the RDM however has the advantage that typical observables are of one- or two-body character and therefore the one-body or two-body RDM are sufficient to determine their expectation values. In essence, employing the -th order RDM reduces the number of relevant degrees of freedom enormously if which is typically the case for many-body systems of identical particles. Since the procedure of deriving the BBGKY hierarchy will be applied also later in order to obtain the equation of motions for the symmetry correlators which represent the analogue of the ITPC for interacting many-particle system (see eq.(1,2)), let us briefly summarize the main steps. The BBGKY hierarchy is derived by taking the time derivative of the RDM of arbitrary degree , employing the underlying time-dependent Schrödinger equation, and dividing the total Hamiltonian into and as well as coupling terms which emerge from the two-body interaction. Here indicate the degrees of freedom of whereas are the ones traced out to obtain . Employing hermiticity of the Hamiltonian and particle exchange symmetry of either purely bosonic or fermionic wave functions one arrives at the hierarchy of EOM

 i∂tρ(n)(Xn;X′n;t)=(H1...n−H1′...n′)ρ(n)(Xn;X′n;t) (9) +(N−n)n∑i=1∫(V(xi,xn+1)−V(x′i,xn+1))ρ(n+1)(Xn,xn+1;X′n,xn+1;t)dxn+1

The so-called collision integral on the r.h.s. of this EOM stems from the two-body interaction coupling the two above-mentioned subspaces. Obviously this is a system of coupled equations and constitutes a hierarchy in the sense that every equation for couples via the integral expression to the next order . Solving this hierarchy, in particular while keeping the advantage has to introduce a cutoff at comparatively low order. At the same time a functional approximation to typically involving the lower order RDM has to be employed. For the reconstruction of as a functional of lower order RDM as well as for relevant properties of the truncated hierarchy of equations, such as particle number and energy conservation as well as compatibility, we refer the reader to Bonitz (); Akbari () and references therein. These features are not in the focus of the present work which aims at establishing the fundamental formalism in the presence of local symmetries. The resulting framework can be employed both in the sense of a propagation method (if an appropriate closure approximation for the hierarchy is utilized) and for obtaining analytical insights into the interplay between local symmetries and interactions.

It is instructive to inspect the lowest order equation separately which takes on the following appearance

 i∂tρ(1)(x1;x′1;t)=(H1−H1′)ρ(1)(x1;x′1;t) +(N−1)∫(V(x1,x2)−V(x′1,x2))ρ(2)(x1,x2;x′1,x2;t)dx2 (10)

where is the single particle part of the total Hamiltonian. Focusing on one can derive the standard continuity equation for the particle density which reflects particle number conservation

 −∇J(x1,t)=∂tn(x1,t)withJ(x,t)=(12i(∇1−∇′1)ρ(1)(x1;x′1;t))∣∣x′1=x1 (11)

Here denotes the particle current density. The EOM (10) for the single particle density matrix coupled to the integrated two-particle density matrix can be exploited to obtain an EOM for the natural orbitals and their occupations Pernal (); Appel (); Brics (); Rapp (); Jansen (); Manthe (). This goes as follows. First one inserts the representation (8) of into the EOM (10). The time-derivative of the natural orbitals is then expanded in the complete set of natural orbitals with . Extracting the global phase from the natural orbital removes the diagonal term . If one subsequently projects onto the natural orbitals one arrives at the EOM for the natural populations and the natural orbitals which take on the following appearance

 i∂tλn=N(N−1)∫dx′dx′′dzϕ∗n(x′,t)[V(x′,z)−V(x′′,z)]ρ2(x′;z;x′′;z;t)ϕn(x′′,t), (12)
 i∂tϕn(x,t)=∑p≠n( ⟨ϕp|H1|ϕn⟩+N(N−1)λn−λp∫dx′dx′′dzϕ∗p(x′,t)× ×[V(x′,z)−V(x′′,z)]ρ2(x′;z;x′′;z;t)ϕn(x′′,t))ϕp(x,t) (13)

Eq.(12) shows that the rate of change of the natural populations is triggered by interaction with the surrounding particles. The time dependence of the natural orbitals, however, is influenced by the single particle Hamiltonian (see first term in eq.(13)) as well as the interaction. The latter terms are multiplied by the inverse of the difference of the relevant occupations indicating that degeneracies111In the case of a degeneracy, any linear combination of the degenerate natural orbitals is also a natural orbital with the same natural occupation number. Thus, the singularity means that the equations of motion cannot “decide” which linear combination is appropriate for a continuous natural orbital trajectory and other criteria have to be employed. Whenever we refer to these equations of motion in the following, we tacitly assume the absence of degeneracies for simplicity. of the occupation numbers lead to corresponding singularities in the expression (13). In particular the rate of change of a given natural orbital is ’driven’ by all other natural orbitals as indicated by the sum over all .

Using now the resolution of the identity for the complete set of natural orbitals and extracting the phase eliminates the diagonal single particle term proportional to . We therefore finally arrive at the following EOM for the natural orbitals

 i∂tϕn(x,t)=H1ϕn(x,t)+N(N−1)∑p≠n(λn−λp)−1ϕp(x,t)∫ϕ∗p(x′,t)[V(x′,z)−V(x′′,z)] ρ(2)(x′,z;x′′,z;t)ϕn(x′′,t)dx′dx′′dz (14)

Eq.(14) will be the relevant equation below to derive the EOM for symmetry correlators on the level of individual natural orbitals. From here on we will focus exclusively on a spin-independent formalism, thereby keeping the same notation as before, for reasons of simplicity. However, we note that the following derivations and results can straightforwardly be transfered and applied to the spin-dependent case.

Iii Equation of Motion for Symmetry Correlators

Having established the theoretical background we are now aiming at performing the relevant steps to derive the corresponding EOM adapted to the presence of local symmetries. First we remind the reader of the fact that for the non-interacting theory in one spatial dimension it is the constancy of the ITPC (see eqs.(1,2)) that reflect the local symmetry in a given spatial domain. The ITPC however depend on two symmetry-related spatial positions for a single particle. The key idea now for going from a local symmetry formalism for non-interacting particles to interacting many-body systems is to define the (canonical) local symmetry correlators as coherences of RDMs between local symmetry related positions. A note of caution is appropriate here: we use the term ’correlator’ in the sense that a corresponding quantity contains symmetry-related spatial points which has to be distinguished from the notion of correlations in many-body physics as a beyond mean field correlated motion of particles. Here, we concentrate on the local symmetry correlators of the one-body RDM. In the following, we will first introduce a local symmetry-adapted representation of the equations of motion for the local symmetry correlators or in other words derive equations of motion for symmetry related coherences of the one-body RDM. Here, the crucial step consists in expressing the kinetic part as a divergence of an appropriately chosen non-local two-point current density. These equations are most transparent in the natural orbital representation of the one-body RDM. The resulting equation of motion generalizes the ITPC (see eq.(2)) to the dynamics of many-body interacting systems, and for this reason, we refer to these local symmetry correlators as the canonical ones. In a second step we will introduce an anomalous local symmetry correlator whose equations of motion are derived and we will show that the non-interacting one-dimensional stationary case of these equations of motion leads to the ITPC (see eq.(1)).

iii.1 The Canonical Form of the Symmetry Correlator

We start by rewriting the single particle part of the EOM for (10) in the following way

 (H1−H1′)ρ(1)(x1;x′1;t)=12(∇′1+∇1)(∇′1−∇1)ρ(1)(x1;x′1;t)+(U(x1)−U(x′1))ρ(1)(x1;x′1;t) (15)

So far we have assumed that the two spatial arguments are independent. Aiming at the presence of local symmetries in and knowing about the structure of the ITPC for the non-interacting one-dimensional case we take . The symmetry mapping corresponds, in general in three dimensions, for to translations by and for to inversions (parity) at the center . Now the kinetic term is reformulated as a divergence of an appropriately defined non-local current density, which can be achieved in a representation-free manner w.r.t. the one-body RDM as summarized in Appendix A. Much more transparent, however, is the reformulation of the kinetic energy on the r.h.s. of eq.(15) in the natural orbital representation:

 (12N)∞∑i=1λi(t)∇x[σ(∇ϕ∗i)(y,t)ϕi(x,t)−ϕ∗i(y,t)(∇ϕi)(x,t)] (16)

where and the derivative in front of the brackets acts on both arguments and . We have therefore succeeded in writing the kinetic energy as a convex sum of divergences of two-point correlator current densities belonging to the natural orbitals

 ji(x,y,t)=[σ(∇ϕ∗i)(y,t)ϕi(x,t)−ϕ∗i(y,t)(∇ϕi)(x,t)] (17)

which is reminescent of the structure of the ITPC obtained in the noninteracting single particle case (see eqs.(1,2)). In total we therefore arrive at the local symmetry correlator EOM

 i∂t[∞∑i=1λi(t)(ϕi(x,t)ϕ∗i(y,t))]=12∇x∞∑i=1λi(t)ji(x,y,t) (18)
 +(U(x)−U(y))[∞∑i=1λi(t)(ϕi(x,t)ϕ∗i(y,t))]+N(N−1)∫(V(x,z)−V(y,z))ρ(2)(x,z;y,z;t)dz

where we assumed that lies in the interior of the domain in order to avoid difficulties regarding the spatial derivatives at the boundary of and . Let us pause to think what we have achieved so far by discussing the above continuity-like equation (18) which is one of the central results of the present work. The l.h.s. of eq.(18) represents the time-derivative of the sum of two point symmetry correlators for the natural orbitals weighted by the natural occupation numbers which, as the natural orbitals do, vary in time. The first term on the r.h.s. is the above-introduced divergence of the sum of natural orbital two-point currents (see eq.(16)) also weighted by the natural populations. The last term on the r.h.s. couples to the second order RDM via the two-body interactions and and will be subject to (functional) approximations in a truncation scheme if one aims at closed equations of motion on this order. The collision integral may be interpreted as a source term in this continuity-like equation for the local symmetry correlators. Finally the second term on the r.h.s. involves the single particle potential . It is exactly this term which disappears domainwise in the presence of corresponding local translation or inversion symmetries i.e. with . This case, which is of central importance here, leads immediately to the simplified equation of motion

 i∂t[∞∑i=1λi(t)(ϕi(x,t)ϕ∗i(y,t))]=12∇x∞∑i=1λi(t)ji(x,y,t) (19)
 +N(N−1)∫(V(x,z)−V(y,z))ρ(2)(x,z;y,z;t)dz

So for any local inversion or translation symmetry, the dynamics of the local symmetry correlator, i.e. one-body coherences between local symmetry related positions, is not directly driven by the external locally symmetric potential but only, if at all, indirectly via the time-evolution of the two-body RDM entering the collision integral. Let us now discuss the role of the range of the interaction potential by decomposing the collision integral in eqs.(18,19). With , we now define the union and as its complement. The collision integral then reads

 T=TD+TE=N(N−1){∫DT+∫E}(V(x,z)−V(y,z))ρ(2)(x,z;y,z;t)dz (20)

describes the two-body interaction () within the local symmetry domain (note that the case of a gapped local symmetry (see Fig.1) needs special attention) whereas is responsible for the interaction from outside (). Let us specifically address which can, by using the symmetry mapping, be rewritten as follows

 TD=N(N−1)(∫D(V(x,z)−V(y,z))ρ(2)(x,z;y,z;t)dz) +N(N−1)(∫D(V(x,F(z′))−V(y,F(z′)))ρ(2)(x,F(z′);y,F(z′);t)dz′) (21)

If the interaction potential depends only on the distance between the particles i.e. we have and obtain

 TD=N(N−1)∫DV(|x−z|)(ρ(2)(x,z;y,z;t)−ρ(2)(x,F(z);y,F(z);t))dz +N(N−1)∫D(V(|x−F(z)|)ρ(2)(x,F(z);y,F(z);t)−V(|y−z|)ρ(2)(x,z;y,z;t))dz (22)

In the above expression the first integral contains the intradomain particle interactions and the second one describes the particle interactions between the domains and . This second term provides only a contribution if the range of the interaction exceeds the minimal distance between and . A special case is the local inversion symmetry in which case we have . Then the contribution reads

 TD=N(N−1)∫D(V(|x−z|)−V(|y−z|))(ρ(2)(x,z;y,z;t)−ρ(2)(x,F(z);y,F(z);t))dz (23)

In ultracold quantum gases Pethick () the interaction in the course of the ultracold collisions can be modeled by a so-called contact potential of the form which is the extreme limit of a short-range interaction. In this case vanishes and the collision integral in total takes the following appearance

 T=N(N−1)g(ρ(2)(x,x;y,x;t)−ρ(2)(x,y;y,y;t)) (24)

This concludes our discussion of the collision integral and we turn now back to our equations of motion for the symmetry correlators. Eq.(18) possesses an integral representation which can be obtained by integrating it over a volume and employing Gauss (divergence) theorem which yields finally

 i∂t∫V[∞∑i=1λi(t)(ϕi(x,t)ϕ∗i(y,t))]dx=12∮S(V)(∞∑i=1λi(t)ji(x,y,t))dS (25)
 +∫V(U(x)−U(y))[∞∑i=1λi(t)(ϕi(x,t)ϕ∗i(y,t))]dx
 +N(N−1)∫V∫(V(x,z)−V(y,z))ρ(2)(x,z;y,z;t)dzdx

where is the surface of the volume and denotes the normal vector of the surface element222A word of caution is in order here: Eq.(25) is valid if lies in the interior of . Otherwise, one has to extend the local coordinate transformation to a global one and carefully inspect whether the kinetic term results in possibly even distribution valued singularities on the surface of .. The time-dependent symmetry correlator change integrated over a volume involves therefore the flux out of this volume due to the convex sum of the weighted currents . Both, the external single particle potential and the interaction occur in their equally volume-integrated form. Again, in the presence of corresponding local symmetries the term due to the external potential vanishes per construction.

In the above EOM (18,19) the symmetry correlators appear in the form of convex sums, i.e. for both the time-derivative as well as the term involving the currents they are weighted by the natural populations and summed over all natural orbitals. It is however possible to derive EOM for the symmetry correlators for individual natural orbitals. Our starting-point is here eq.(14). Constructing the quantity where again are symmetry related via and reformulating the action of the Laplacian as a total divergence as well as exploiting the hermiticity of the density matrix we arrive at the desired EOM for a single natural orbital-based symmetry correlator

 i∂t(ϕn(x,t)ϕ∗n(y,t))=12∇xjn(x,y,t)+(U(x)−U(y))ϕn(x,t)ϕ∗n(y,t) +N(N−1)∞∑p=1,p≠n(λn−λp)−1(Ipnϕp(x,t)ϕ∗n(y,t)+Inpϕn(x,t)ϕ∗p(y,t)) (26)

where is the single natural orbital symmetry correlator current (see eq.(17)) and denote the matrix elements of the collision integral in the natural orbital basis

 (27)

Local symmetry implies again and one should note that in this case the EOM for the individual symmetry correlators depend only indirectly on the external potential, namely via as we have already observed above for the total local symmetry correlator. Eq.(26) shows that the dynamics of the symmetry correlator of a single natural orbital involves apart from the corresponding correlator current in particular all off-diagonal correlators involving the other natural orbitals. We remark that a corresponding EOM can be derived for the dynamics of these off-diagonal correlators as done above for the diagonal term. In passing, we note that by considering the trivial map , eq.(26) turns into a continuity equation with a local source term for the individual natural orbital densities. This source term is induced by the collision integral, which couples the natural orbital densities to the other natural orbitals, and vanishes when being integrated over the whole space, which reflects the conservation of probability.

In the following subsections we will derive and discuss several special cases for the above equation. This includes first the non-interacting case where we derive a continuity equation for the two-point symmetry correlators and recover the case of the ITPC for a time-independent theory. The two subsequent subsections discuss the mean field theory for bosons and fermions, respectively.

iii.1.1 Non-interacting Particle Systems

Let us first focus on the case of non-interacting particles in the presence of local symmetries which leads via eq.(19) to

 i∂t[∞∑i=1λi(t)(ϕi(x,t)ϕ∗i(y,t))]=12∇x∞∑i=1λi(t)ji(x,y,t) (28)

and via eq.(26) to

 i∂t(ϕi(x,t)ϕ∗i(y,t))=12∇xji(x,y,t) (29)

This is a generalized continuity equation for the two-point symmetry correlators. Reducing to a single particle in one dimension which is the case previously treated extensively Kalozoumis2013a (); Kalozoumis2014a (); Kalozoumis2014b (); Kalozoumis2013b (); Kalozoumis2015a (); Morfonios2014 (); Kalozoumis2015b (); Zambetakis2016 (); Kalozoumis2016a (); Wulf2016 () we obtain

 i∂t[ϕ(x,t)ϕ∗(y,t)]=12∂xj(x,y,t) (30)

with where one should always keep in mind that and the derivative acts therefore on both and . Focusing on the stationary i.e. time-independent Schrödinger (or Helmholtz) equation with the single particle energy eigenvalue and the stationary single particle bound or scattering state the time derivative vanishes due to the time independence of the product . Integrating the r.h.s. of eq.(30) yields then

 12ij(x,y,t)=12i[σ(∂xϕ∗)(y,t)ϕ(x,t)−ϕ∗(y,t)(∂xϕ)(x,t)]=C (31)

which is identical to the complex conjugate of the previously discovered Kalozoumis2013b (); Kalozoumis2014a () invariant two-point current correlator given in eq.(2). The latter is constant in each corresponding domain of local symmetry. Our general theoretical framework therefore contains the already known special case of stationary states in one dimension. This remark refers to the ITPC . The case of the ITPC will be treated separately below since it goes beyond the canonial ket bra structure of the density matrix formalism.

iii.1.2 Mean-field theory: Bosons

Let us focus next on the case of a bosonic many-body ensemble where we have in particular in mind ultracold quantum gases which can in many situations successfully be described by a mean field Gross-Pitaevskii equation that assumes the particles to be condensed in a single orbital Pitaevskii (); Pethick () i.e. that there is a natural orbital with . The many-body wave function then reads which trivially obeys the particle exchange symmetry. The resulting single particle and two particle density matrices read and respectively. The interaction among the bosons can in the ultracold regime be described by a so-called contact potential which depends on a single parameter, the -wave scattering length contained in . Assuming again a locally symmetric potential the above EOM (19) for the symmetry correlator simplifies to

 i∂t(ϕ(x,t)ϕ∗(y,t))=12∇xj(x,y,t)+(N−1)g[|ϕ(x,t)|2−|ϕ(y,t)|2]ϕ(x,t)ϕ∗(y,t) (32)

with given in eq.(17) for a single natural orbital. Now, due to the interactions on the mean field level a ’source term’ of a nonlinear character emerges. This term is proportional to the difference which shows the remarkable fact that if the magnitude of the wave function follows the local symmetry then the symmetry correlator behaves as if there would be no interactions present on the mean field level. This is of particular importance for the case of a stationary state for which eq.(32) leads to

 12∇xj(x,y)=(N−1)g[|ϕ(y)|2−|ϕ(x)|2]ϕ(x)ϕ∗(y) (33)

This is reminescent of the occurence of the so-called symmetric perfectly transmitting resonances Kalozoumis2013b () which indeed possess the property that the magnitude of the wave function follows the local symmetries of the underlying potential .

iii.1.3 Mean-field theory: Fermions

Within the mean field theory of fermions, namely Hartree-Fock theory, electrons occupy orbitals forming a Slater determinant that is completely antisymmetric with respect to particle exchange (as indicated previously, the spin degrees of freedom are ignored tacitly here). We have therefore natural orbitals with occupation numbers equal to one. Exploiting once again local symmetry of the external potential the corresponding EOM reads

 i∂t(N∑i=1ϕi(x,t)ϕ∗i(y,t))=12∇xN∑i=1ji(x,y,t) (34)
 +∫(V(x,z)−V(y,z))N∑i,j=1[ϕi(x,t)ϕ∗i(y,t)|ϕj(z,t)|2−ϕj(x,t)ϕ∗i(y,t)ϕi(z,t)ϕ∗j(z,t)]dz

where now due to the Pauli principle natural orbital symmetry correlator currents contribute. The interaction term contains the corresponding Coulomb and exchange terms, respectively. The stationary case follows immediately as

 12∇xN∑i=1ji(x,y)=∫(V(y,z)−V(x,z))N∑i,j=1[ϕi(x)ϕ∗i(y)|ϕj(z)|2−ϕj(x)ϕ∗i(y)ϕi(z)ϕ∗j(z)]dz

We do not address the problems which arise due to the mean field approximation in the equations of motion but refer for this to the literature (see refs. Bonitz (); Akbari () and references therein).

iii.2 The Anomalous Symmetry Correlator

The theoretical framework presented so far is derived from the BBGKY-hierarchy which has the canonical ket bra structure such that the scalar product of Hilbert space can be exploited. The stationary state case leads then to the canceling of the dynamical phase factor and single (natural) orbital orthogonality within e.g. mean field theory can be exploited. As a matter of fact, however, the ITPC in eq.(1) is not of this canonical structure but is still a relevant quantity as has been shown via the generalization of the parity and Bloch theorem Kalozoumis2014a () or in applications to parity-time symmetry breaking of scattering states Kalozoumis2014b () where plays a decisive role to construct the underlying phase diagram. So, in what follows, we ’leave safe grounds’ based on the canonical structure and establish the anomalous EOM of the generalization of to higher dimensional many-body interacting systems. Though this might seem disturbing at first glance, it is indeed motivated by the above-indicated importance of the ITPC .

Let us start by defining the quantity

 γ(N)(x1,...,xN;x′1,...,x′N;t)=Ψ(x1,...,xN,t)Ψ(x′1,...,x′N,t) (35)

which is a complex symmetric quantity as compared to the hermitian density matrix . We can then construct the corresponding reduced quantities of order

 γ(n)(Xn;X′n;t)=∫Ψ(X′n,Xcn,t)Ψ(Xn,Xcn,t)dXcn (36)

In order to derive an EOM for we define the time-derivative 333This time-derivative is reminiscent of a multiple-time calculus .

 i^∂tγ(N)=i(∂tΨ(XN,t))Ψ(X′N,t)−i(∂tΨ(X′N,t))Ψ(XN,t) (37)

and correspondingly for where the time derivative acts under the spatial integral. Following the same scheme as in the derivation of the BBGKY hierarchy, i.e. the steps which have been outlined above (see section II.2), we can again partition the Hamiltonian in the equation of motion for into degrees of freedom of and those which are traced out, i.e. integrated over. Partial integration is now sufficient to get rid of the terms where the corresponding Hamiltonian acts on the integrated spatial variables and, together with permutation symmetry, we obtain the EOM for

 i^∂tγ(n)(Xn;X′n;t)=(H1...n−H1′...n′)γ(n)(Xn;X′n;t) (38) +(N−n)n∑i=1∫(V(xi,xn+1)−V(x′i,xn+1))γ(n+1)(Xn,xn+1;X′n,xn+1;t)dxn+1

and for the lowest order equation for

 i^∂tγ(1)(x1;x′1;t)=(H1−H1′)γ(1)(x1;x′1;t) (39) +(N−1)∫(V(x1,x2)−V(x′1,x2))γ(2)(x1,x2;x′1,x2;t)dx2

Due to the fact that is complex symmetric (and not hermitian) it is not necessarily diagonalizable. However, we do assume here that it can be diagonalized in which case it has an eigendecomposition which reflects the complex symmetry Horn (). More precisely, diagonalization is possible if and only if the corresponding eigenvector matrix can be chosen such that and , where is complex orthogonal and is the diagonal eigenvalue matrix. The complex orthogonality of reflects the complex symmetry of . This leads us to the representation

 γ(1)(x1;x′1;t)=∑iμi(t)