Quantum Gross-Pitaevskii Equation

We introduce a non-commutative generalization of the Gross-Pitaevskii equation for one-dimensional quantum gasses and quantum liquids. This generalization is obtained by applying the time-dependent variational principle to the variational manifold of continuous matrix product states. This allows for a full quantum description of many body system ---including entanglement and correlations--- and thus extends significantly beyond the usual mean-field description of the Gross-Pitaevskii equation, which is known to fail for (quasi) one-dimensional systems. By linearizing around a stationary solution, we furthermore derive an associated generalization of the Bogoliubov -- de Gennes equations. This framework is applied to compute the steady state response amplitude to a periodic perturbation of the potential.

In 1961, Gross and Pitaevskii developed the mean-field theory description of cold Bose gasses [1][2][3], which resulted in the ubiquitous Gross-Pitaevskii equation (GPE) where φ(x, t) is the order parameter of the Bose-Einstein condensate in a trapping potential v(x). Ever since, this equation constitutes the cornerstone for the theoretical description of cold atomic gasses [4][5][6]. Its success can be explained by the fact that the GPE agrees with the full quantum solution for the three-dimensional problem in the weak-density limit, which corresponds with the typical experimental setup of trapped dilute gasses (see Ref. 7 and references therein). The GPE without trapping potential [v(x) = 0] is also known as the nonlinear Schrödinger equation and appears in several areas of theoretical physics as it offers a canonical description for slowly varying, quasi-monochromatic wave packets in dispersive, weakly nonlinear media [8,9]. As such, it has also stimulated an abundance of mathematical research towards showing the stability of its (solitary wave) solutions [10][11][12], as well as towards the development of numerical integrators [13][14][15][16][17].
The goal of the current Letter is to formulate a generalization of the one-dimensional GPE which allows to take quantum correlations into account. One approach to derive the GPE is by applying the Dirac-Frenkel time-dependent variational principle (TDVP) [35][36][37] to a variational mean field ansatz |Ψ[φ] in the canonical is the bosonic field creation operator in second quantization and |Ω is the vacuum state of the Fock space with no particles. Our generalization is based on the variational set of continuous matrix product states (cMPS) where the functions Q and R take values in C D×D and v 1,2 are D-dimensional boundary vectors for a system defined on the interval [x 1 , x 2 ] ⊂ R (either boundary could be at infinity) [38,39]. This ansatz 'generates' a physical state by acting with 'operators' Q(x) and R(x) on a virtual boundary system C D [40], whose dimension D controls the amount of quantum correlations in the variational ansatz. The path-ordered exponential Pe arises because its argument does not commute at different positions and reduces to the standard exponential for D = 1, where this ansatz reproduces the mean field ansatz with φ(x) = R(x) and e x 2 x 1 Q(x) dx some global scalar factor (norm and phase). This ansatz was conceived as a continuum limit of the matrix product state (MPS) ansatz [41][42][43], which underlies the highly successful density matrix renormalization group [44] for the description of one-dimensional quantum spin systems. By enlarging the refinement parameter D, the exact quantum state can be increasingly well approximated. Indeed, the cMPS ansatz was shown to capture both the ground state [39] and the two types of elementary excitations [45] of the Lieb-Liniger model for moderate values of D. In this Letter, we use the TDVP to generalize the GPE and to approximate the full quantum dynamics by a time-evolving cMPS.
For a complex manifold, the TDVP can be understood geometrically as a replacement of Schrödinger's equation with an additional projector P Ψ that projects the right hand side to the tangent space of the variational manifold at the point |Ψ . Whereas the Schrödinger equation for the exact evolution would immediately take an initial state away from the variational manifold, this extra projector assures that the evolution remains within the manifold at any time.

Quantum Gross-Pitaevskii Equation -
We now specialize to the set of cMPS with fixed virtual dimension D and use the Lieb-Liniger Hamiltonian [26] The energy expectation value was computed in Refs. 38 and 39 and is given by with the left and right density matrices ρ L (x) and ρ R (x) defined by and where can be interpreted as a covariant derivative, as explained below.
To facilitate the rest of the derivation, we also introduce the notationM (y, z) = Pe z y dx Q(x)⊗1+R(x)⊗ψ † (x) , using which a general cMPS tangent vector can be written as In particular, we find for the left-hand side of Eq.
where we use the dot notation to denote the time derivative and omit the explicit t-dependence of the variational parameters for the sake of brevity. The orthogonal projection in the right hand side of Eq. (1) corresponds to finding the tangent vector which minimises |Φ[V, W, w 1 , w 2 ] −Ĥ |Ψ . The solution to this problem can be found explicitly and is constructed in the Supplementary Material. This then gives rise to the following set of differential equations to which we henceforth refer as the Quantum Gross-Pitaevskii Equation (QGPE).
In these equations, P (x) is an arbitrary x-dependent (and t-dependent) D × D matrix, which originates from the gauge freedom that is inherent to the cMPS parameterization. Indeed, the physical cMPS which we henceforth refer to as a gauge transformation. Invariance of the energy expectation value in Eq. (3) follows trivially from the induced transformation laws , where the last transformation validates the identification of D x with a covariant derivative (in the adjoint representation).
thus explaining the appearance of the free parameters P (x) in the QGPE [Eq. (11)]. The QGPE is then invariant under arbitrary xand t-dependent gauge transformations of the form of Eq. (8), provided that we also transform P (x) as where, as before, we have omitted the time-dependence throughout these equations. Indeed, identifying P and Q with respectively the A 0 and A 1 components of a gauge potential, the left hand side of Eq. (7a) can be written as covariant time derivative in the adjoint representation, whereas the left hand side of Eq. (7b) can be recognized as the only non-zero of the antisymmetric field tensor. Note that in all of this the gauge gauge group is the general linear group GL(D) and that this gauge invariance is purely virtual.
Translation invariance -For translation invariant systems with v(x) = v, we can immediately take the thermodynamic limit and approximate the quantum state of the system with a cMPS with space-independent matrices Q and R. The left and right density matrices correspond to the left and right fixed point (eigenvector corresponding to the eigenvalue with largest real part) of the transfer matrix Q ⊗ 1 1 + 1 1 ⊗ Q + R ⊗ R, and the time evolution is given by the set of first order non-linear differential equations for some suitable choice of P . These equations can rather straightforwardly be integrated and the gauge freedom can help to further simplify the problem to some extent.
Boundary Conditions -For finite systems, the QGPE needs to be supplemented with appropriate boundary conditions to fully specify the problem. These will also affect the evolution equation for the boundary vectors v 1,2 , which have so far not been written down yet. As the TDVP can be obtained from extremizing an action, there are only two types of self-consistent boundary conditions (similar to e.g. the classical wave equation for a vibrating string), unless explicit boundary terms are included in the Hamiltonian from Eq. (2). These can be derived by considering the quantized field operatorψ(x) and expressing stability with respect to variations ofψ † (x). Either the value of ψ(x) is fixed at the boundaries, resulting in the Dirichlet conditions: Alternatively, ifψ(x) is left free at the end points, then variations with respect toψ(x) † will -after partial integration in the kinetic energy-enforce the homogenous Neumann conditions: We could of course also have mixed boundary conditions at x 1 and x 2 . While the resulting boundary conditions for the variational parameters R are inherently gauge invariant, they only correspond to D instead of D 2 equations each. In e.g. the case of Dirichlet conditions, they only fix one eigenvalue and eigenvector of the matrix R. In fact, the other directions of R at the boundary do not appear in physical expectation values and can thus not be fixed from physical considerations. In order to eliminate any interplay with the gauge transformation, we will 'promote' the boundary conditions in a gauge invariant manner by imposing them as identity matrix, e.g. R(x 1 ) = a1 D in the case of Eq. (12a). Note that there are no separate boundary condition on Q(x), as these degrees of freedom can be interpreted as pure gauge degrees of freedom and can be set to zero by inserting in Eq. (8) the transformation G(x) = Pe x 2 x Q(z) dz .
The boundary conditions then also affect the TDVP equation for the boundary vectors. For the case of Dirichlet where the left hand side contains the covariant time derivative in the conjugate and fundamental representation, respectively. These equations are also valid for the Neumann conditions, where the right hand side becomes zero.
Discussion -Let us now discuss in more detail the mathematical structure of the QGPE. Since it contains the quantities ρ L,R (x), which are defined by integrating Eq. (4), it forms a set of coupled non-linear partial integro-different equations containing first order time derivatives and second order space derivatives. It is a noncommutative generalization of the normal GPE in that it is defined in terms of matrix variables. Indeed, the normal GPE is recovered from Eq. (7a) in the limit D = 1 by setting R(x) = φ(x) and observing that multiplications then commute and commutators thus vanish. In that limit, x2 x1 Q(x) dx acts as on overall scalar factor (norm and phase) that can be absorbed in the boundaries.
Just like the normal GPE and essentially any TDVP equation, the QGPE for real-time evolution takes the form of a classical Hamiltonian system. The expectation value Ψ[Q, R, v 1 , v 2 ]|Ĥ|Ψ[Q, R, v 1 , v 2 ] plays the role of the classical Hamiltonian function and will be a constant of motion for solutions of the QGPE when the Hamiltonian H is time-independent. The symplectic structure comes for free, since any complex submanifold of Hilbert space is automatically a Kähler manifold using the induced metric. This Kähler metric corresponds to the physical overlap of two tangent vectors, as computed in the Supplementary Material. When using the QGPE with imaginary time evolution t → −iτ to find a cMPS approximation for the quantum ground state, this symplectic structure is of course lost and the energy expectation value decreases monotonically until convergence.
We have already discussed the gauge invariant structure of the QGPE as it was of key importance in correctly deriving and interpreting the equations. So far, we have not discussed how to fix the gauge degrees of freedom, which cor-responds to specifying P (x) in the QGPE. While P (x) = 0 represents a valid choice, there might be reasons to make a different choice. The gauge freedom in the cMPS allows to transform the parameters into e.g. a left-canonical form where Q(x)+ Q(x) † + R(x) † R(x) = 0 and consequently ρ L (x) = 1 1 [46]. We can then choose P such that this form is preserved at any point in time. This implicitly determines P fromQ(x) +Q(x) † + R(x) †Ṙ (x) +Ṙ(x) † R(x) = 0. One can for example check that that by choosing P (x) = −iR(x) † D x R(x) + iK(x) we obtain (now also omitting the argument x) which is purely antihermitian for any hermitian matrix K(x) and thus satisfies the above gauge fixing condition. The freedom of K originates from the fact that the left canonical form does not fix the gauge completely but only reduces the gauge group from GL(D) to U(D). We can further choose K such that the above expression becomes zero and the left canonical form is even maintained under imaginary time evolution. Developing a stable numerical integration scheme for the QGPE is equally challenging. Firstly, there are the the inherent complexities associated with solving a set of nonlinear partial integro-differential equations. Because of the second order spatial derivative and first order time derivative, the Courant-Friedrichs-Levy condition limits the time step of explicit schemes. A typical workaround for the GPE is to use a splitting scheme [13,15,17], where the evolution is decomposed into the local terms (external potential and interaction) and the kinetic term. The linearity of the latter allows to solve it using a Crank-Nicolson method [47] or in Fourier space [48]. While the QGPE is still linear in the second order spatial derivative, it has non-linear terms containing first-order spatial derivatives, which cannot easily be integrated in Fourier space. Another complication of the QGPE, as formulated in Eq. (11), is that it depends on the inverses of the density matrices ρ L (x) and ρ R (x). These become rank deficient near respectively the left and right boundary, as is clear from the definition in Eq. (4). It is well known in the tensor network community that the gauge degrees of freedom in the underlying matrix product state have to be exploited to transform these density matrices into identity matrices [49]. This was implemented only recently for the TDVP equation for matrix product states [50], using a non-trivial decomposition of the tangent space projector that allows to split the non-linear differential equation for all variables into a set of linear differential equations for the individual MPS tensors [51], which is made possible by the fact that the MPS parameterization is multilinear. Here too, we have to face complications introduced by the intrinsic nonlinearity in the cMPS parameterization. A final challenge is to develop suitable continuum analogous of the factorization routines such as the QR-or singular value decomposition, which are exploited in the MPS simulations to robustly implement the required gauge transformations. Note, in addition, that in order to exploit gauge freedom, the discretization of the QGPE should not break gauge invariance. Hereto, ideas from lattice gauge theory can serve as inspiration. We develop a possible numerical integrator for the QGPE that deals with these complexities in a forthcoming publication.
Outlook and conclusion -We have developed a natural generalization of the GPE based on the formalism of cMPS, which should be better suited to study onedimensional quantum systems where entanglement plays an important role and the mean-field ansatz underlying the GPE is not justified. While it would be interesting to have a complete mathematical study of the existence, uniqueness (up to gauge transformations) and stability of the solutions of the QGPE for given boundary conditions, this is beyond the goal and scope of this paper. It would also be enlightening to investigate whether there exist quantum generalization of certain exact solutions of the GPE, such as the dark soliton solution, and how they relate to the integrability of the full Lieb-Liniger Hamiltonian or to the topological excitations constructed using the related ansatz of Ref. 45. In addition, it would be instructive to compare the predictions of the QGPE to existing beyond-mean field studies such as e.g. Ref. 52.
As the cMPS ansatz is a versatile variational ansatz that can readily be extended to bosonic and fermionic systems with e.g. multiple particle species [38,[53][54][55], the QGPE equations generalize straightforwardly to such systems. Moreover, the approach described in this paper is in no way restricted to the Lieb-Liniger Hamiltonian, but is applicable to arbitrary Hamiltonians in one spatial dimension. cMPS methods have already been used to study (1+1) dimensional relativistic theories for fermions [53] and bosons [56], in thermodynamic limit and in the translationally invariant setting. Using the regularisation method described in Ref. 56, the derivation of the QGPE presented in this paper extends straightforwardly to such systems, enabling in principle the study of general bosonic non-linear σ-models with boundaries. Non-linear σ-models are of significant interest for the high energy physics community, for example, providing the underlying description of bosonic strings [57,58] propagating on curved backgrounds. Given that cMPS methods are intrinsically non-perturbative, the approach of this paper has the potential to be particularly useful in the study of such models in the limit of large curvature, when perturbative quantization schemes fail. It is furthermore very encouraging that the implementation of Dirichlet and Neumann conditions in the quantum theory is straightforward [Eqs. (12) - (13)], and we thus expect that the boundary condition implementation described in this paper can be used "as is" for the description of strings with either freely propagating endpoints, or with endpoints restricted to lie on D-branes [59,60].
Finally, a challenging open problem would be to construct a higher-dimensional generalization of this QGPE. This would arise as the TDVP equation for the variational set of continuous projected-entangled pair states [61], which are less well studied and understood.
We acknowledge discussions with C. Lubich. Research supported by the Research Foundation Flanders (JH), the EPSRC under grant number EP/L001578/1 (VS), the Austrian FWF SFB grants FoQuS and ViCoM, the cluster of excellence EXC 201 "Quantum Engineering and Space-Time Research" and the European grants SISQ, QUTE and QFTCMPS.

SUPPLEMENTARY MATERIAL: DERIVATION OF THE QGPE
To properly deal with the boundary conditions, it is useful to derive the QGPE following the general recipe of the TDVP, i.e. as the Euler-Lagrange equations corresponding to extremizing the classical action We can easily derive the Euler-Lagrange equations by considering variations with respect to the complex conjugates of the variational parameters, which are treated as independent and appear only in the bra. Using our definition of tangent vectors, this immediately leads to the condition for any possible variation. This is indeed equivalent to the geometric formulation of the TDVP in Eq. (1). More generally, it tells us that we can find the tangent space projection |Φ[V, W, w 1 , for all possible V ′ (x) = δQ(x), W ′ (x) = δR(x) and w ′ 1,2 = δv 1,2 . The left hand side contains the overlap of two different tangent vectors and is given by where the terms on the first 5 lines corresponds to all contributions of non-zero V and W , and the terms on lines 6 and 7 correspond to non-zero w 1 and non-zero w 2 , respectively. Here, we have introduced a new notation E(x, y) = P exp y x Q(z) ⊗ 1 1 + 1 1 ⊗ Q(z) + R(z) ⊗ R(z) dz .
It is the continuum equivalent of the (product of) MPS transfer matrices and allows to e.g. write Let us now first computeĤ |Ψ[Q, R, v 1 , v 2 ] itself. Using the rules from Ref. 38 and by applying partial integration to the kinetic energy term, we obtain 20) Given the linearity of the tangent space projector, we can compute the projection of the 4 different terms separately and add the result: Hence, for the boundary conditions here considered, only the second term of Eq. (20) needed to be projected onto the tangent space and, when considering the full Schrödinger equation, would be responsible for taking the exact evolution out of the manifold.
Grouping everything together gives rise to |Φ[V, W, w 1 , which we have to equate to i |Φ[Q,Ṙ,v 1 ,v 2 ] . Together with the gauge invariance of Eq. (9), this gives rise to the QGPE in Eq. (11) and the corresponding evolution of the boundary vectors Eq. (14).