Continuous Matrix Product States for Quantum Fields:
an Energy Minimization Algorithm
Abstract
The generalization of matrix product states (MPS) to continuous systems, as proposed in the breakthrough paper [F. Verstraete, J.I. Cirac, Phys. Rev. Lett. 104, 190405(2010)], provides a powerful variational ansatz for the ground state of strongly interacting quantum field theories in one spatial dimension. A continuous MPS (cMPS) approximation to the ground state can be obtained by simulating an Euclidean time evolution. In this Letter we propose a cMPS optimization algorithm based instead on energy minimization by gradient methods, and demonstrate its performance by applying it to the Lieb Liniger model (an integrable model of an interacting bosonic field) directly in the thermodynamic limit. We observe a very significant computational speedup, of more than two orders of magnitude, with respect to simulating an Euclidean time evolution. As a result, much larger cMPS bond dimension can be reached (e.g. with moderate computational resources) thus helping unlock the full potential of the cMPS representation for ground state studies.
Over the last 25 years, progress in our understanding of quantum spin chains and other strongly interacting quantum manybody systems in one spatial dimension has been dominated by a variational ansatz: the matrix product state (MPS) Fannes et al. (1992); White (1992); Schollwöck (2011); Verstraete et al. (2009). The wave function of a quantum spin chain made of spin degrees of freedom depends on complex parameters ,
(1) 
Accordingly, an exact numerical simulation has a computational cost that grows exponentially with the size of the chain. In an MPS, the coefficients are expressed in terms of the trace of a product of matrices. For instance, in a translation invariant system the MPS reads
(2) 
where and are complex matrices. Thus, the state of spins is specified by just variational parameters, allowing for the study of arbitrarily large, even infinite, systems McCulloch (2008); Vidal (2007).
A generic state of the spin chain can not be expressed as an MPS, because the bond dimension limits how entangled can be. However, ground states of local Hamiltonians happen to be weakly entangled (e.g. they obey an entanglement area law Hastings (2007); Eisert et al. (2010)) in a way that allows for an accurate approximation by an MPS. Given a Hamiltonian , White’s revolutionary density matrix renormalization group (DMRG) White (1992, 1993) algorithm provided the first systematic way of obtaining a ground state MPS approximation by minimizing the energy, see also White (1993). Subsequently, Refs. Vidal (2004, 2003) proposed an algorithm to simulate time evolution with an MPS, which in Euclidean time also produces a ground state approximation, see also Refs. Daley et al. (2004); White and Feiguin (2004); Feiguin and White (2005). An improved formulation of the time evolution simulation by MPS was obtained in terms of the timedependent variational principle (TDVP) Haegeman et al. (2011).
The continuous version of an MPS (cMPS), introduced by Verstraete and Cirac Verstraete and Cirac (2010); Haegeman et al. (2013), has the potential of duplicating, in the context of quantum field theories in the continuum, the enormous success of the MPS on the lattice. A cMPS expresses the wave function of a quantum field on a circle of radius as a path ordered exponential of the fields that define the theory. For a bosonic, translation invariant system it reads
(3) 
where is the bosonic field creation operator,
(4) 
is the empty state, i.e. , and and are complex matrices. Again, the wave function is parameterized by just parameters. A cMPS approximation to the ground state of a continuum Hamiltonian can then be obtained by simulating an Euclidean time evolution with TDVP adapted to cMPS Haegeman et al. (2011). While this algorithm and its variations work reasonably well for small up to Haegeman et al. (2010); Draxler et al. (2013); Rincón et al. (2015); Draxler et al. (2017), their performance is poor compared to lattice MPS techniques.
In this Letter we propose an energy minimization algorithm to find a cMPS approximation for ground states, based on gradient descent techniques, and demonstrate its performance with the Lieb Liniger model in the thermodynamic limit (). We also propose a useful cMPS initialization scheme, of interest on its own, based on lattice MPS algorithms. These proposals result in a very significant computational speedup with respect to Euclidean time evolution – e.g. converging a cMPS with bond dimension requires less time than a computation with TDVP. For simplicity we consider a single bosonic field. Generalization to a fermionic field and to multiple fields is straightforward.
Continuum limit and central canonical form.— In order to describe the algorithm, we must first adjust the notation in two ways. Firstly, following Verstraete and Cirac (2010), we discretize the interval in Eq.(3) into a regular lattice made of sites and with intersite spacing , and produce an MPS with matrices and mor () given, in vectorized form, by
(5) 
such that the original cMPS is recovered in the limit Verstraete and Cirac (2010). Here, and corresponds to having 0 or 1 particle at the lattice site. This lattice visualization is useful in order to manipulate the cMPS with regular MPS techniques, provided the latter have a welldefined continuum limit (). Secondly, we use the lattice visualization to reexpress the cMPS of an infinite system () in the central canonical form Stoudenmire and White (2013), Eq.(7) below. For this purpose, we consider the Schmidt decomposition of according to a left/right partition of the resulting infinite lattice Inf (),
(6) 
and denote by a diagonal matrix with the Schmidt coefficients in its diagonal. In the central canonical form, the MPS is expressed as the infinite product of (vectorized) matrices
(7) 
The matrices and are chosen such that
(8) 
are in the left and right canonical form Schollwöck (2011), namely
(9)  
(10) 
From Eqs.(7)(8) the standard MPS form Eq.(2) (in the limit) for e.g. a left normalized MPS is recovered,
(11) 
In the central canonical form, familiar to DMRG and MPS practitioners working with socalled singlesite updates, a change in the matrices and on a single site produces an equivalent change in , in the sense that the scalar product in the lattice Hilbert space and in the effective onesite Hilbert space are equivalent (they are related by an isometry). This is important when applying gradient methods, because two gradients, calculated in two different gauges of the same state, are in general not related by a gauge transformation and are not equivalent. The importance of the central gauge has been realized early on in DMRG White (1992) and also time evolution methods Vidal (2003, 2007); Stoudenmire and White (2013); Haegeman et al. (2016); van (); Halimeh and ZaunerStauber (2016).
Gradient descent.— Given a quantum field Hamiltonian , see e.g. Eq.(15), our goal is to iteratively optimize the cMPS in such a way that the energy
(13) 
is minimized. Each iteration updates a triplet and is made of two steps. (i) First, keeping fixed, we update and in the direction of steepest descent given by the gradient, namely
(14) 
where is some adjustable parameter and denotes complex conjugation. Crucially, the gradients and can be efficiently computed using standard cMPS contraction techniques. We dynamically choose the largest possible factor by requiring consistency with some simple stability conditions (alternatively, can be determined by a line search). (ii) Then, from we obtain by bringing the cMPS representation back into the central canonical form. This completes an iteration, which has a cost comparable to one time step in TDVP. We emphasize that all manipulations are implemented directly in the continuum limit i.e. is treated as an analytic parameter throughout the optimization, and the limit can be taken exactly due to exact cancellation of all divergencies.
Overall, the proposed energy minimization algorithm proceeds as follows (see sup () for technical details).
(a) Initialization: An initial triplet of matrices is obtained, either from a random initialization or, as in this Letter, through Eq.(12) from an MPS optimized on the lattice.
(b) Iteration: The above update is iteratively applied until attaining a suitably converged triplet .
(c) Final output: A standard cMPS representation as in Eq.(3) is recovered by transforming the result into . For instance, as in Eq.(8) for a final cMPS in the left canonical form (see also sup ()).
As usual in such optimization methods, convergence can be accelerated by replacing the gradient descent in Eq.(14) with e.g. a nonlinear conjugate gradient update, which reuses the gradient computed in previous steps (see Milsted et al. (2013); sup ()).
Example.— To benchmark the above algorithm, we have applied it to obtain a cMPS approximation to the ground state of the Lieb Liniger model Lieb and Liniger (1963); Lieb (1963),
(15) 
which is both of theoretical and of experimental interest and has been realized in several cold atom experiments Jaksch et al. (1998); Bloch (2005); Kinoshita et al. (2006); Cazalilla et al. (2011); Haller et al. (2010); Meinert et al. (2016). This integrable Hamiltonian has a critical, gapless ground state that can be described by Luttinger liquid theory Cazalilla et al. (2011) and can be exactly solved by Bethe ansatz Lieb and Liniger (1963); Lieb (1963); Prolhac (2016); Lang et al. (2016); Ristivojevic (2014); Zvonarev (2010).
Fig. 1 (a) (blue dots) illustrates the fast and robust convergence of the cMPS with the number of iterations of steepest descent, by showing the energy density , particle density , and reduced energy density ,
(16) 
for bond dimension and the choice of parameters . For comparison, we also show the same quantities when the cMPS is optimized instead by an Euclidean time evolution using the TDVP algorithm (green crosses), starting from the same initial state and using values for TDVP and for the steepest descent optimization sup (), respectively. These values for are typically used in common TDVP calculations for cMPS Draxler (). Fig. 1 (b) then shows the convergence of the energy to the exact value obtained from the Bethe ansatz solution eco () as a function of iteration number, again for a steepest descent (blue dots) and TDVP (green crosses) optimization. In this example, energy minimization converges towards the ground state roughly a hundred times faster than TDVP. The difference in performance is even bigger for larger bond dimension , and/or when no lattice optimization is used to initialize the cMPS, in which case TDVP may even fail to converge.
Fig. 2(a) illustrates the performance of the proposed energy minimization algorithm as a function of the bond dimension . For , we computed the reduced energy density for several values of the dimensionless interaction strength in the range and observed a uniform pattern of convergence towards the exact . For reference, a optimization employing a nonlinear conjugate gradient optimization sup () (stopped once the energy has converged to 9 digits) takes 6 minutes on a desktop computer myp (), including both the lattice initialization (2 minutes) and the nonlinear conjugate gradient optimization in the continuum (4 minutes). This value of the bond dimension is the largest reported so far using TDVP Draxler et al. (2013); Stojevic et al. (2015).
Fig. 2(b)(c) specializes to , , and considers even larger values of the bond dimensions , up to , to reproduce wellunderstood finite effects of the cMPS representation Stojevic et al. (2015); Pirvu et al. (2012); Tagliacozzo et al. (2008); Pollmann et al. (2009). Fig. 2(b) shows the relative error in the reduced energy density and the entanglement entropy across a left/right bipartition, Eq.(6). As expected, vanishes with as a powerlaw, , whereas the entanglement entropy diverges logarithmically, . Fig. 2(c) shows the superfluid correlation function , which is seen to saturate to a finite value at some distance , another wellunderstood artefact of the (c)MPS representation at finite bond dimension Stojevic et al. (2015); Pirvu et al. (2012); Tagliacozzo et al. (2008); Pollmann et al. (2009). This artificial finite correlation length is seen to diverge with growing as a powerlaw, .
Once we have established that the optimized cMPS is an accurate approximation to the ground state, we can move to exploring other properties of the model. Fig. 3 shows the superfluid correlation function and pair correlation function , respectively, for and different values of the dimensionless interaction strength . With growing we observe an increasingly rapid decay in the superfluid correlation function. The pair correlation function develops typical oscillations that are related to the fermionic nature of the ground state of the TonksGirardeau gas Girardeau (1960) at .
We can also estimate both the central charge and the Luttinger parameter , which can be used to uniquely identify the conformal field theory that characterizes the universal low energy / large distance features of the model. The central charge can be estimated from the slope of (see Ref. Stojevic et al. (2015)). For we obtain a value of , to be compared with the exact value . The Luttinger parameter Cazalilla (2004); Cazalilla et al. (2011) is obtained from fitting vs Rincón et al. (2015), where we choose to lie in the region where exhibits powerlaw decay. For and we obtain . A value of was obtained in Cazalilla (2004) from the weakcoupling approximation of the Bethe ansatz solution. The relative difference to our result is .
Discussion.— The cMPS is a powerful variational ansatz for strongly interacting quantum field theories in 1+1 dimensions Verstraete and Cirac (2010). In this Letter we have proposed a cMPS energy minimization algorithm with much better performance, in terms of convergence and the attainable bond dimension , than previous optimization algorithms based on simulating an Euclidean time evolution. For benchmarking purposes, we have applied it to the exactly solvable Lieb Liniger model, but it performs equally well for a large variety of (nonexactly solvable) field theories sup (). We envisage that this algorithm will play a decisive role in unlocking the full potential of the cMPS representation for ground state studies in the continuum.
Our algorithm works best by initializing the cMPS through an energy optimization on the lattice and by translating the resulting MPS from the lattice to the continuum through Eq.(5). A natural question is then whether the continuum algorithm is needed at all. That is, perhaps –one may wonder– an MPS algorithm working at finite lattice spacing can already provide a cMPS representation (through Eq.(5)) that can be made arbitrarily close to the one obtained with the continuum algorithm by decreasing sufficiently. The answer is that this is not possible: lattice algorithms necessarily become unstable as the lattice spacing is reduced. This can be understood from a simple scaling argument. In discretizing e.g. Hamiltonian of Eq.(15) into a lattice, the nonrelativistic kinetic term is seen to diverge with as , while the rest of terms in the Hamiltonian have a milder scaling. For small this creates a large range of energy scales that lead to numerical instability. This effect is compounded with a second fact, revealed by Eq.(5). For small the MPS matrix is made of two pieces: a constant part made of ’s and ’s and the variational parameters , which are of order . Thus the first part shadows the second one, in that the numerical precision on the variational parameters is reduced by a factor when embedded in matrix . The observant reader may then wonder if these problems could be prevented by just changing variables, to work instead with . This is indeed the case, and also the essence of working with the cMPS representation directly, as we do in the proposed energy minimization algorithm. Notice that lattice MPS techniques can be succesfully applied to ground states Dolfi et al. (2012); Stoudenmire et al. (2012); Wagner et al. (2013, 2014); Baker et al. (2015) and real time evolution Knap et al. (2014); Muth and Fleischhauer (2010) of discretized field theories. However, these simulations are conducted at sufficiently large and are often plagued with finite scaling analysis, which is not necessary when working directly with a cMPS.
We have seen that the cMPS energy minimization algorithm drastically outperforms TDVP at the task of approximating the ground state (we emphasize that the TDVP remains an extremely useful tool e.g. to simulate real time evolution, for which no other method exists). This result did not come as a surprise: on the lattice, MPS energy minimization algorithms, including DMRG, have long been observed to converge to the ground state much faster than time evolution simulation algorithms McCulloch (2008). We expect the new algorithm to also produce a significant speedup both for inhomogeneous Hamiltonians (where matrices and depend on space Verstraete and Cirac (2010); Haegeman et al. (2015)), or for a theory of multiple fields Haegeman et al. (2010); Quijandría et al. (2014); Quijandría and Zueco (2015); Chung et al. (2015); Chung and Bolech (2016); Draxler (2016); Haegeman (2011)), as we will discuss in future work.
Acknowledgements The authors thank D. Draxler, V. ZaunerStauber, A. Milsted, J. Haegeman, F. Verstraete, and E. M. Stoudenmire for useful comments and discussions. The authors also acknowledge support by the Simons Foundation (Many Electron Collaboration). Computations were made on the supercomputer Mammouth parallèle II from University of Sherbrooke, managed by Calcul Québec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), the ministère de l’Économie, de la science et de l’innovation du Québec (MESI) and the Fonds de recherche du Québec  Nature et technologies (FRQNT). Some earlier computations were conducted at the supercomputer of the Center for Nanophase Materials Sciences, which is a DOE Office of Science User Facility. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development & Innovation.
References
 Fannes et al. (1992) M. Fannes, B. Nachtergaele, and R. F. Werner, Communications in Mathematical Physics 144, 443 (1992).
 White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
 Schollwöck (2011) U. Schollwöck, Annals of Physics 326, 96 (2011).
 Verstraete et al. (2009) F. Verstraete, J. I. Cirac, and V. Murg, Adv. Phys. 57,143 (2008) (2009).
 McCulloch (2008) I. P. McCulloch, arxiv:0804.2509 [condmat] (2008).
 Vidal (2007) G. Vidal, Phys. Rev. Lett. 98, 070201 (2007).
 Hastings (2007) M. B. Hastings, Journal of Statistical Mechanics: Theory and Experiment 2007, P08024 (2007), ISSN 17425468.
 Eisert et al. (2010) J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. 82, 277 (2010).
 White (1993) S. R. White, Phys. Rev. B 48, 10345 (1993).
 Vidal (2004) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
 Vidal (2003) G. Vidal, Phys. Rev. Lett. 91, 147902 (2003).
 Daley et al. (2004) A. J. Daley, C. Kollath, U. SchollwÃ¶ck, and G. Vidal, Journal of Statistical Mechanics: Theory and Experiment 2004, P04005 (2004).
 White and Feiguin (2004) S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
 Feiguin and White (2005) A. E. Feiguin and S. R. White, Phys. Rev. B 72, 020404 (2005).
 Haegeman et al. (2011) J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Phys. Rev. Lett. 107, 070601 (2011).
 Verstraete and Cirac (2010) F. Verstraete and J. I. Cirac, Phys. Rev. Lett. 104, 190405 (2010).
 Haegeman et al. (2013) J. Haegeman, J. I. Cirac, T. J. Osborne, and F. Verstraete, Phys. Rev. B 88, 085118 (2013).
 Haegeman et al. (2010) J. Haegeman, J. I. Cirac, T. J. Osborne, H. Verschelde, and F. Verstraete, Phys. Rev. Lett. 105, 251601 (2010).
 Draxler et al. (2013) D. Draxler, J. Haegeman, T. J. Osborne, V. Stojevic, L. Vanderstraeten, and F. Verstraete, Phys. Rev. Lett. 111, 020402 (2013).
 Rincón et al. (2015) J. Rincón, M. Ganahl, and G. Vidal, Phys. Rev. B 92, 115107 (2015).
 Draxler et al. (2017) D. Draxler, J. Haegeman, F. Verstraete, and M. Rizzi, Phys. Rev. B 95, 045145 (2017).
 (22) In Eq. (5) we neglect other possible matrices for that are not relevant to the current discussion.
 Stoudenmire and White (2013) E. M. Stoudenmire and S. R. White, Phys. Rev. B 87, 155137 (2013).
 (24) In the limit of an infinite periodic lattice, , we can use the Schmidt decomposition of as in Eq.(6), which assumes a lattice with open boundary conditions. Further, in the case of a translation invariant system on an infinite chain, all left/right partitions are equivalent.
 Haegeman et al. (2016) J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Phys. Rev. B 94, 165116 (2016).
 (26) L. Vanderstraeten, Online Lecture Notes, http://quantumtensor.pks.mpg.de/index.php/school/.
 Halimeh and ZaunerStauber (2016) J. C. Halimeh and V. ZaunerStauber, arXiv:1610.02019 [condmat] (2016), arXiv: 1610.02019.
 (28) Supplemetary Material.
 Milsted et al. (2013) A. Milsted, J. Haegeman, and T. J. Osborne, Physical Review D 88, 085030 (2013).
 Lieb and Liniger (1963) E. H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963).
 Lieb (1963) E. H. Lieb, Phys. Rev. 130, 1616 (1963).
 Jaksch et al. (1998) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
 Bloch (2005) I. Bloch, Nat. Phys. 1, 23 (2005).
 Kinoshita et al. (2006) T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 900 (2006).
 Cazalilla et al. (2011) M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, and M. Rigol, Rev. Mod. Phys. 83, 1405 (2011).
 Haller et al. (2010) E. Haller, R. Hart, M. J. Mark, J. G. Danzl, L. Reichsollner, M. Gustavsson, M. Dalmonte, G. Pupillo, and H.C. Nägerls, Nature 466, 597 (2010).
 Meinert et al. (2016) F. Meinert, M. Knap, E. Kirilov, K. JagLauber, M. B. Zvonarev, E. Demler, and H.C. Nägerl, arXiv:1608.08200 [condmat] (2016).
 Prolhac (2016) S. Prolhac, J. Phys. A, 50, 144001 (2017).
 Lang et al. (2016) G. Lang, F. Hekking, and A. Minguzzi, arXiv:1609.08865 [condmat, physics:mathph] (2016).
 Ristivojevic (2014) Z. Ristivojevic, Phys. Rev. Lett. 113, 015301 (2014).
 Zvonarev (2010) M. Zvonarev, Notes on Bethe Ansatz (2010), URL http://cmt.harvard.edu/demler/TEACHING/Physics284/LectureZvonarev.pdf.
 (42) D. Draxler, private communication.
 (43) As reference value we take , where is the converged particle density at large iteration number of the TDVP or the gradient optimization, respectively.
 (44) Run times were measured on a 2014 Intel i7vPro quad core CPU with 8GB RAM.
 Stojevic et al. (2015) V. Stojevic, J. Haegeman, I. P. McCulloch, L. Tagliacozzo, and F. Verstraete, Phys. Rev. B 91, 035120 (2015).
 Pirvu et al. (2012) B. Pirvu, G. Vidal, F. Verstraete, and L. Tagliacozzo, Phys. Rev. B 86, 075117 (2012).
 Tagliacozzo et al. (2008) L. Tagliacozzo, T. . R. de Oliveira, S. Iblisdir, and J. I. Latorre, Phys. Rev. B 78, 024410 (2008).
 Pollmann et al. (2009) F. Pollmann, S. Mukerjee, A. M. Turner, and J. E. Moore, Phys. Rev. Lett. 102, 255701 (2009).
 Girardeau (1960) M. Girardeau, J. Math. Phys. 1, 516 (1960).
 Cazalilla (2004) M. A. Cazalilla, J. Phys. B: Atomic, Molecular and Optical Physics 37, S1 (2004), ISSN 09534075.
 Dolfi et al. (2012) M. Dolfi, B. Bauer, M. Troyer, and Z. Ristivojevic, Phys. Rev. Lett. 109, 020604 (2012).
 Stoudenmire et al. (2012) E. M. Stoudenmire, L. O. Wagner, S. R. White, and K. Burke, Phys. Rev. Lett. 109, 056402 (2012).
 Wagner et al. (2013) L. O. Wagner, E. M. Stoudenmire, K. Burke, and S. R. White, Physical Review Letters 111, 093003 (2013).
 Wagner et al. (2014) L. O. Wagner, T. E. Baker, E. M. Stoudenmire, K. Burke, and S. R. White, Phys. Rev. B 90, 045109 (2014).
 Baker et al. (2015) T. E. Baker, E. M. Stoudenmire, L. O. Wagner, K. Burke, and S. R. White, Phys. Rev. B 91, 235141 (2015).
 Knap et al. (2014) M. Knap, C. J. M. Mathy, M. Ganahl, M. B. Zvonarev, and E. Demler, Phys. Rev. Lett. 112, 015302 (2014).
 Muth and Fleischhauer (2010) D. Muth and M. Fleischhauer, Phys. Rev. Lett. 105, 150403 (2010).
 Haegeman et al. (2015) J. Haegeman, D. Draxler, V. Stojevic, J. I. Cirac, T. J. Osborne, and F. Verstraete, arXiv:1501.06575 (2015).
 Haegeman (2011) J. Haegeman, PhD thesis, Ghent University, Faculty of Sciences, Ghent, Belgium. (2011).
 Quijandría et al. (2014) F. Quijandría, J. J. GarcíaRipoll, and D. Zueco, Phys. Rev. B 90, 235142 (2014).
 Quijandría and Zueco (2015) F. Quijandría and D. Zueco, Phys. Rev. A 92, 043629 (2015).
 Chung et al. (2015) S. S. Chung, K. Sun, and C. J. Bolech, Phys. Rev. B 91, 121108 (2015).
 Chung and Bolech (2016) S. S. Chung and C. J. Bolech, arXiv:1612.03149 (2016).
 Draxler (2016) D. Draxler, PhD thesis, University of Vienna (2016).
 (65) M. Ganahl, J. Rincón, and G. Vidal, in progress.
 (66) Note that special care has to be taken when using the scipy.linalg.eigs function to calculate left and right eigenvectors of the transfer operator, due to erroneous behaviour of the routine.
 (67) Wikipedia entry on Nonlinear conjugate gradient method.
 (68) Ashley Milsted, private communication.
Supplementary Material
Appendix A DMRG preconditioning
As mentioned in the Main Text, our cMPS gradient optimization works best when using porperly prepared initial states. The purpose of the following chapter is to introduce a novel method that employs DMRG energy minimization to obtain cMPS on a discretized system. These discrete states are then taken as initial states to a cMPS optimization in the continuum. The main result here is that DMRG can be used to optimize an MPS with tensors of the form
(S1) 
in such a way that the cM PS structure, i.e. the matrices , is always explicitly available. We use instead of to highlight that this is a lattice MPS. Note that even though we are treating soft core bosons here, we are restricting the local Hilbert space dimension to (see below). We start by discretizing the Hamiltonian from Eq.(15) on a grid with spacing . With the replacement , and using a first order discretization of the kinetic energy, this results in the following local Hamiltonian terms:
(S2)  
(S3)  
(S4) 
Note that we have replaced the local interaction term with a nearest neighbour interaction term. This replacement becomes exact in the limit . By comparing the cMPS expressions for the local energy density
where and are the left and right reduced density matrices (see below), with the ones obtained using the above discretization Eqs.(S2)(S4) and MPS tensors of the form Eq.(S1), it is easy to see that they can be made to coincide if in the numerical derivative Eq.(S2) one makes the replacement
(S5) 
where is a projector onto the state with no particles at position . Thus, using Eq.(S3) and Eq.(S4) and a modified kinetic energy Eq.(A) with MPS matrices of the form Eq.(S1) is almost identical to a discretized cMPS description as obtained in Haegeman et al. (2015). The difference lies in the normalization: whereas MPS tensors Eq.(S1) are normalized according to a finite , a true cMPS is normalized with respect to . The use of the projector is crucial to the restriction of the local Hilbert space dimension to . The final Hamiltonian then assumes the form
For a given , we then use standard iDMRG on a two site unit cell McCulloch (2008) to obtain the ground state. We note that our proposed gradient optimization for cMPS could, with minor modifications, also be applied to the discrete setting, and to homogeneous lattice systems in general Ganahl et al. (). Once the iDMRG has been sufficiently converged, one can use the matrices to initialize a new DMRG run for a finer discretization by simply changing to in Eq.(S1), and renormalizing the MPS tensor. Using the terminology of Ref.Dolfi et al. (2012), we call this operation a prolongation of the state. Even for as large as and for the parameter values used in this Letter, this gives already good initial states for a cMPS optimization at . Note that since the matrices depend on , so does the quality of the prolonged initial state. This proposed fine graining procedure generalizes the approach presented in Dolfi et al. (2012) to arbitrary prolongations. Prior to prolonging the state to a finer discretization, it is vital to bring it into the canonical form Vidal (2003) to remove any gauge jumps at the boundaries of the twosite unit cell. Note that in the diagonal gauge there is still a gauge freedom with diagonal, complex gauge matrices that has to be removed as well before prolongation. These gauge jumps originate from a unitary freedom in the SVD and QR decompositions, which are used heavily in lattice MPS optimizations.
Appendix B Determination of the Gradient for the Lieb Liniger model
In this section we elaborate on the determination of the gradient of the energy expectation value using cMPS methods (see also van () for details on lattice MPS gradients). The goal is to optimize the matrices to lower the energy
see Eq(13) in the main text. Similar to DMRG Schollwöck (2011), we calculate the gradient of this expression with respect to . Taking the derivative gives
We normalize our state in such a way that (see Eqs.(S8) and (S9) below), and thus the second term vanishes identically. Since the energy depends nonlinearly on , we have to apply the chain rule when calculating . However, due to translational invariance, this full gradient decomposes into a sum of identical terms, where each term corresponds to a local gradient, and is the number of discretization points:
is a local derivative, similar to DMRG, with respect to . Note that for inhomogeneous systems, this local gradient would contain derivatives of the matrix Haegeman et al. (2013). Knowledge of this local gradient is thus sufficient to calculate the gradient of . To obtain the local gradient , it is convenient to consider the state in the center site form
where . is a diagonal matrix containing the Schmidt values. and are left/right normalized cMPS matrices, i.e.
Here we have defined the left/right cMPS transfer operator Haegeman et al. (2013)
(S6) 
The energy expectation with respect to now local (at the position where we want to take the local derivative) can be evaluated to
(S7) 
Symbols represent tensor contractions. In DMRG, the matrices and in Eq.(B) are known as the left and right block Hamiltonians. Below we describe how to calculate them in our context. We use a particular renormalization of these left and right block Hamiltonians, such that
(S8)  
(S9) 
Due to this renormalization, the expectation value . The local derivative of with respect to is then given by
(S10)  
(S11)  
signs denote free matrix indices due to the removal of a matrix at this point in the tensor network. and are then used to update :
At this point, the new state is of the form
The matrices are now absorbed back into the cMPS tensors from e.g. the right side, i.e.
are then brought back into the central canonical gauge (see below), which produces a new triplet of matrices and completes one update step. Convergence is measured by the norm of the gradient
This measure is similar to, but more stringent than, the one used in the TDVP, in general (see below for a definition of , the latter does not take into account).
Next we address the calculation of the matrices and in Eqs.(S10) and (S11). In the context of DMRG, they are known as the left and right block Hamiltonians, respectively. For an infinite, translational invariant system they are given by an infinite sum of the form
(S12)  
(S13) 
where is an infinitesimal MPS transfer operation, i.e.
(see Eq.(S6)) and
is related to the energy content of an infinitesimal interval . Since has a left and right eigenvector and to eigenvalue 1 (and similarly has eigenvectors and ) the two sums in Eqs.(S12) and (S13) are divergent. These divergences can be regularized by removing the subspace from and , and from and , i.e. by replacing
The geometric series can then be summed to
(S14) 
(note the cancellation of in Eqs.(B)) which is equivalent to
(S15)  
(S16) 
with the left/right cMPS transfer operator Eq.(S20). We solve Eqs.(S15) and (S16) using the lgmres routine provided by the scipy.sparse.linalg module pyt ().
Here is a summary of our proposed optimization scheme:

Initialize a cMPS with matrices and bring them into central canonical gauge. Set a desired convergence .

Update:

Regauge into the central canonical form and measure . If stop, otherwise go back to 2.
During simulations we dynamically adapt the parameter , i.e. if we detect a large increase in we reject the update and redo the step with a smaller value of .
Appendix C Time Dependent Variational Principle for cMPS
In the following we summarize the Time Dependent Variational Principle (TDVP) Haegeman et al. (2011, 2016); Halimeh and ZaunerStauber (2016) for homogeneous cMPS (see also the Appendix in Rincón et al. (2015) for more details). To obtain the ground state of e.g. Eq.(15) in the Main Text, one employs Euclidean time evolution using the TDVP. The general strategy is the following: given a cMPS at time , one finds an approximation that can be used to evolve :
The most general ansatz ansatz for is to superimpose local perturbations on the state :
(S17) 
where we is
and are irrelevant boundary vectors at . The optimal and are found from minimizing the norm
(S18) 
with defined in Eq.(C). This amounts to the calculation of tensor network expressions similar to Eqs.(S10) and (S11). The nonlocal nature of the perturbations in leads to a complication in the minimization, that can be overcome by resorting to a new parametrization (or ) such that the norm has nonvanishing contributions only when and act at the same position in space. Two possible parametrizations are given by
(S19) 
These are called the left or right tangentspace gauges, respectively. We note that this parametrization covers the full tangentspace to the state Haegeman et al. (2013). Here, the only free parameter left is . and are the left and right reduced density matrices obtained as the left and right eigenvectors of the cMPS transfer operator
(S20) 
to eigenvalue , i.e. in braket notation. In the TDVP for cMPS, one further simplifies all expressions in the minimization by reparametrizing as
where is now the free parameter. The convergence measure used in TDVP is given by the norm
For left normalized matrices , the update step is given by
with a small time step in imaginary time. and correspond to and in the gradient optimization approach.
Our gradient calculation includes one additional computational step as compared to the TDVP. There, depending on the choice of tangentgauge, the calculation of one of the two expressions containing or in Eqs.(S10) and (S11) can be omitted. In our case we need to calculate both of them. This adds however only a small overhead in computational time.
Appendix D Obtaining the central canonical form
An important ingredient to our gradient optimization is the central canonical form, introduced in the Main Text. Here we describe how to obtain the central canonical form for cMPS. The procedure is a straight forward adaption of the corresponding lattice algorithm. Starting from unnormalized matrices , one first calculates the left and right eigenvectors ( in matrix notation) of the transfer operator Eq.(S20) to the largest real eigenvalue ,
where . Next, the cMPS is normalized by
This shifts the eigenvalue of and to . Using , compute and use a singular value decomposition to obtain
where contains the Schmidt values. Normalize by
The left or right normalized cMPS matrices or are then obtained from
with
(S21) 
Note that .
Appendix E Nonlinear conjugate gradient
The nonlinear conjugate gradient approach is an improvement of the steepest descent method that reuses gradient information from previous iterations, see e.g. nlc (), and has recently been proposed as an improvement of the TDVP Milsted et al. (2013). At iteration , the new search direction is given by a linear combination of the current gradient and the previous search direction from step :
(S22) 
The parameter is calculated from a semiempirical formula. In our case, we use the FletcherReeves form
(S23) 
In a regular implementation, one uses a line search method to find the minimum in the direction of , which is then taken to be the starting point for the next iteration. We omit this line search in our implementation. The procedure is initialized with . In our case, there is a small caveat due to a gauge change of the state between iterations and when regauging the state into the central canonical form. Due to this gauge change, the current gradient and the previous search direction can only be added after the latter has been transformed into the gauge of the current state ash (). This is achieved using the matrices and from Eq.(D). Consider an update for the state