Matrix Product State Representation without explicit local Hilbert Space Truncation with Applications to the Sub-Ohmic Spin-Boson Model

We present an alternative to the conventional matrix product state representation, which allows us to avoid the explicit local Hilbert space truncation many numerical methods employ. Utilising chain mappings corresponding to linear and logarithmic discretizations of the spin-boson model onto a semi-infinite chain, we apply the new method to the sub-ohmic SBM. We are able to reproduce many well-established features of the quantum phase transition, such as the critical exponent 1/2 predicted by mean-field theory. Via extrapolation of finite-chain results, we are able to determine the infinite-chain critical couplings at which the transition occurs and, in general, study the behaviour of the system well into the localised phase.


Introduction
The spin-boson model (SBM) describes a single two-level system (TLS), a spin, coupled to environmental degrees of freedom represented by a continuous bath of bosonic field modes. It is one of the most important models for studying the general effects arising when a quantum system is coupled to an environment [1]. In the sub-ohmic version, it possesses a mean-fieldlike quantum phase transition between a localized and a delocalized phase. This quantum phase transition has been the subject of extensive numerical and analytical investigations [2][3][4][5][6][7][8][9]. Yet, many numerical approaches face challenges near and above the transition to the localized phase. This is due to the rapidly rising number of field excitations in the localized phase, which imply that the quantum states of the field modes span an increasingly large subspace of their full Hilbert space. Most numerical methods are however based on local Hilbert space truncation and are hence discarding increasing amounts of vital information.
In this paper, we describe a variation on the matrix product state (MPS) representation that avoids the explicit local Hilbert space truncation using a soft cut-off instead. We demonstrate the usefulness of this approach by applying it to study the properties of the second-order magnetic quantum phase transition of the sub-ohmic SBM and their comparison to an analytical approach based on a variational ansatz [2].

Matrix product state method
Generally, the compound quantum state of an M-level system (whose states are labelled k = 1, . . . , M) and an environment of N bosonic modes (with states labelled i 1 , . . . , i N ∈ N 0 ) can be described by where |0 represents the vacuum state of all modes and b † m is the creation operator of the mth mode.
This representation of the state is known as a MPS [10][11][12][13]. In this form the state of the M-level system is represented by the M matrices S (1) , . . . , S (M) and the mth field mode is represented by a semi-infinite set of matrices {B (i m ) m }. In this infinite number of matrices lies the problem of the traditional MPS ansatz. Numerical calculations are limited to a finite set of matrices, hence the local Hilbert space associated with each mode has to be truncated to a finite size by limiting the local dimension. Under certain conditions, such as high mean excitation numbers, this truncation can lead to substantial errors in numerical calculations.
Here we consider an alternative MPS-type representation that retains the ability to represent correlations between subsystems and avoids the hard truncation of the local Hilbert space dimensions. This is achieved by reducing the number of matrices per mode to a single matrix X m , which is defined such that for 0 i m < ∞, i.e. the infinite set {B (i m ) m } is now formed from powers of a single matrix X m , reducing the total number of matrices required to fully describe the state to M + N . 4 The additional factor of (i m !) −1 is chosen to simplify later calculations. In addition, instead of directly restricting the bosonic Hilbert space to a finite set of low occupation sates, it introduces a 'soft cut-off' in the number of allowed bosons, giving lower occupational states a higher weight than states with large boson number, and vanishing weight in the limit i m → ∞. It should be noted that this representation shares some semblance to a MPS representation in a variable coherent state basis which becomes transparent when diagonalizing the matrix X . As we admit arbitrary forms for X that are obtained in the optimization part of our algorithm, our approach does not require a specific choice of basis but determines the optimal choice automatically.
Substituting into (2), the MPS can be written as This form not only enables us to avoid a direct truncation of the bosonic Hilbert space, but also continues to allow for straightforward determination of the normalization as well as expectation values. and we find for the norm The exponentiation of the X m matrices in combination with the deliberate choice of the (i m !) −1 factor results in the exponential functions appearing in equation (5), which are straightforward to evaluate numerically.
As an example of an expectation value, we find for the population of the kth mode In general, we find that all quantities of interest are simple traces over products of χ 2 × χ 2 matrices.
A ground state MPS in the form of equation (4) can be found for an arbitrary Hamiltonian H by starting with a state |ψ formed of randomly chosen S and X matrices and then minimizing the energy with respect to |ψ , i.e. finding the S and X matrices which minimize equation (9). The approach we are taking here is heuristic. It is motivated by the desire to reduce the number of free parameters in the description of the ansatz wavefunction, while retaining the essential features of the physical wavefunction. Indeed, if we consider a single system only, then we find that every state admits a representation as in equation (4) for sufficiently large (possibly infinite) matrices. We expect, but have not proven, that this remains true for multipartite states too. Correlations between subsystems can be described increasingly well by using growing matrix dimensions, following the philosophy of MPSs. Besides this, the demonstrated computational efficiency and the lack of a hard cut-off are additional points in favour of this approach.
To further demonstrate the viability and usefulness of this ansatz we now apply this method to analyse the ground state properties of the sub-ohmic SBM.

Spin-boson model
The Hamiltonian of the (unbiased) SBM is given by (h = 1) where σ i are the usual Pauli matrices describing a TLS with tunnelling amplitude . a l and a † l are the bosonic annihilation and creation operators of the environment, which consists of bath modes with frequency ω l . The key quantity in the description of the system-environment interaction is the spectral function J (ω) = π l g 2 l δ(ω − ω l ). Here we consider a spectral function of the form as given in [14], where ω c is the maximum cut-off frequency of the spectrum and (ω) = 1 − (−ω) = 1 for ω > 0. In the following we focus exclusively on the sub-ohmic case for which 0 < s < 1, in particular on the case s < 0.5, which we will compare to existing results in the literature. In [3,15] it was shown that a Hamiltonian of the form equation (10) can be mapped exactly onto a semi-infinite chain of bosonic modes that experience nearest neighbour interaction only, with the system only coupling to the first chain site. The transformed Hamiltonian can be written as where the coupling strength between the TLS and the first site is given by and the local energies and tunnelling amplitudes of the sites are and respectively. This mapping brings several advantages, which include the analytical forms for the parameters of the resulting chain model given above, an intuitive picture of how irreversibility emerges [16] and the ready applicability of the MPS method [2,3]. Utilizing the MPS equation (4), we can now find the ground state of the spin-boson system by minimizing equation (9). Specifically, we have to find the S and X which minimize with the TLS's local energy the system-chain interaction energy and the chain energy subject to the constraint ψ|ψ = 1. In terms of the MPS description we obtain, after truncating the chain length to N sites, the total energy This minimization can be carried out numerically to yield the full ground state MPS of the TLS and the N -site chain. In the following we will present some results for the ground state properties of the sub-ohmic SBM. The minimizations in this work were carried out using MATLAB's fminunc function.  (3) c as N → ∞. The dotted line for s = 0.5 corresponds to logarithmic discretization ( = 1.5) and the short dashed lines are the critical valuesα c predicted by the polaron ansatz in [2]. The discrepancy between our results and the predicted values hint at a likely inaccuracy of the polaron ansatz, particularly at low s.

Results
The sub-ohmic SBM is believed to possess a mean-field-like continuous phase transition in the magnetization for 0 < s < 0.5 at a critical coupling strength α c between system and environment. For small coupling strengths α < α c the TLS is in a delocalized phase, having no net magnetization. Above the critical point α > α c the environment induces a spontaneous magnetization on the TLS, which then exhibits a doubly degenerate localized phase.
In the delocalized phase the mean site population is expected to diverge along the chain. This has so far made it difficult to study the system close to and above the phase transition with great accuracy employing numerical methods that rely on Hilbert space truncation, such as numerical renormalization group (NRG) [4] and density matrix renormalization group (DMRG) methods [17,18]. We propose that this truncation and many of the associated problems can be avoided using the MPS of form equation (4) and the minimization equation (20).
However, in line with the other approaches mentioned above, one still has to perform a different kind of truncation to make the numerical simulation feasible, namely truncation of the chain length N , which effectively amounts to an infrared cut-off that is neglecting low frequency components of the environmental bath. Figure 1 shows the critical coupling strengths α (χ ) c , defined as the value of α at which the TLS develops a non-zero magnetization, plotted against the inverse chain length 1/N for matrices of dimension χ = 3 and several values of s, which we determined using our ansatz. A good fit to the data was found for an ansatz of the form  which we then fitted for each s, enabling us to extract an extrapolated value for α (χ ) c in the limit N → ∞. Table 1 shows the results of the extrapolation. In [2] a variational ansatz was used to predict some of the properties of the sub-ohmic SBM ground state. For the critical coupling strength they predict a valuẽ These predicted values are indicated in figure 1 by dashed lines and are also listed in table 1 along with the fractional deviation between our extrapolation results and the predicted values. We find that the predictions agree reasonably well with our result for large s but show significant deviations at smaller s. This suggests that the mean-field type ansatz for the environment holds well for larger s but fails for decreasing s, possibly as a result of the increasing correlations in the environment. In figure 2 we show the behaviour of the critical coupling strength for different matrix dimensions χ in the case of s = 0.2. Even for scalars, χ = 1, the system still exhibits a The growing error with increasing χ is due to the fact that we used the same computation time for all matrix sizes, and hence employed less stringent convergence criteria for higher χ, leading to larger uncertainties on the individual data points. Our value for α (4) c is in excellent agreement with the critical coupling α c = 0.0175 ± 0.0002 found in [8] using quantum Monte Carlo simulations for the same set of parameters, showing that, even for very moderate matrix dimensions χ , our method agrees with previous studies. In the following analysis we use χ = 2 and 3 (to speed up simulations) since in the rest of this paper we are mainly interested in the qualitative features of the system, which as mentioned above show no significant deviations for all χ 2 in the mean-field regime 0 < s 0.5.
Being able to find an MPS representation for the ground state, we were also able to analyse the general properties of the state in both the delocalized and the localized phase. Figure 3 shows the magnetization M = | σ z | of the TLS for some representative values of s. This and the following results were obtained with a chain length N = 50. For s 0.3 we used χ = 3, whereas the results for higher s were obtained using χ = 2 matrices due to slower convergence just above the phase transition. Figure 3 clearly shows the two phases, separated by a second-order transition. In the delocalized phase α < α c the order parameter M is zero. Above the critical coupling strength α c the TLS obtains a finite magnetization with a tendency to full localization M = 1 as α grows large. This localized phase is two-fold degenerate with M = ± σ z both being solutions. where we expect a deviation of the exponent from the mean-field value. In this plot the solid straight line represents the exponent 0.286 our result converges to, which appears to be in reasonable agreement with the results presented in [9]. This suggests that the method presented in this paper is also applicable outside the mean-field regime.
Mean-field theory predicts a second-order magnetic transition at α c with Our simulations do indeed reproduce the correct critical mean-field exponent 1 2 well, as can be seen in figure 4, where we have plotted the magnetization for α > α c on a log-log-plot. The mean-field result is indicated in the figure by the solid straight line. To show that the method is also valid in the non-mean-field regime of the sub-ohmic SBM, 0.5 < s < 1, we consider the specific case of s = 0.75, using matrices of dimension χ = 3. The results for this case are shown in the inset of figure 4. From figure 2 in [9] we expect to find an exponent of approximately 0.25. Our result predicts a scaling according to an exponent of 0.286, which is in reasonable agreement with [9], and certainly shows a clear deviation from the mean-field result. This suggests that our method is not limited to the mean-field regime, but can also be applied in other cases.
In [2] a variational ansatz was used to predict the amount of entanglement in the TLS, defined as the von Neumann entropy of its reduced density matrix. Our numerical results are presented in figure 5. Despite the slight deviations in the values for α c , which most likely arise as a combination of the finite chain length we considered as well as the inherent differences between the two approaches (cf table 1), our results are in excellent qualitative agreement with the analytical predictions [2]. In the delocalized phase entanglement increases. At the critical coupling α c the entanglement exhibits a cusp and then decays rather rapidly in the localized phase due to the system evolving into a product state.
In addition to the entanglement of the TLS, we also looked at the entanglement of the individual sites in the chain. As a representative example, the results for s = 0.3 are shown  . Entanglement E = −tr[ρ n log ρ n ] of the individual chain sites as a function of α for the first ten sites for s = 0.3 using χ = 3. The entanglement of site n is here defined as the von Neumann entropy of the reduced density matrix ρ n of site n. Only the first few sites show significant amounts of entanglement. This shows that the analytical ansatz in [2] does have a sound basis in this regime but fails to be accurate near the critical point and perhaps for other choices of s when entanglement is larger. in figure 6 for the first ten chain sites. We find that only the first few sites carry significant amounts of entanglement and that the general behaviour with changing α closely resembles the entanglement properties of the TLS in figure 5. A substantial spread of entanglement along the chain is only observable very close to the phase transition.
Another observable of interest is the coherence σ x which is shown in figure 7 as a function of α. The results are again in excellent qualitative agreement with the results from the variational approach [2]. The coherence is continuously decreasing, with a faster decay above the transition, at which we observe a cusp.
As mentioned above, the main reason other numerical approaches have failed to return accurate results near and above the critical coupling is that, whereas in the delocalized phase the mean occupation b † n b n of chain site n rapidly decreases along the chain, it rises considerably in the localized phase. In fact, [2] predicts that b † n b n diverges along the chain for α > α c . In figure 8 we plot b † n b n as a function of α and n. Below the transition we find that indeed the average population of the sites decreases with n. At the transition we observe a sudden increase in b † n b n and the maximum begins to shift away from the first site, further along the chain. This rapid rise in the occupancy shows why methods such as DMRG, which rely on Hilbert space truncation, are challenged in this regime, since the information about the system is spread over an increasing number of basis states, only a finite number of which are retained by these methods.
An alternative method to the linear chain mapping we have thus far considered, is provided by logarithmic discretization of the spectrum [3,4,14,19], which is extensively used in NRG. It does not linearly subdivide the bath spectral function, but instead splits it in intervals [ −(n+1) , −n ], where > 1 is the discretization parameter and n ∈ N 0 . This new Hamiltonian has again the same form equation (12), but with different site frequencies and transition amplitudes and respectively, where ζ s , A n , C n and N n are given in [3]. To be applicable for numerical methods it is again necessary to truncate the resulting chain Hamiltonian at a finite number of sites n = N . A challenge for many numerical methods with the logarithmic discretization is the fact that the mean occupation of the chain sites is on average considerably larger than on the linearly discretized chain. This quickly leads to the breakdown of these methods. However, as equation (4) avoids any direct truncation, we were again able to establish numerical results for the sub-ohmic SBM ground state. The dotted line in figure 1 shows the same extrapolation for s = 0.5 as carried out for the linear discretization, now using the logarithmically discretized Hamiltonian with discretization parameter = 1.5. Despite the missing data for large N , we see that the critical coupling converges to a similar value as found before, α LD c = 0.0725 ± 0.0029, with a fractional deviation |α c − α LD c |/α LD c = 0.033, where the superscript LD refers to logarithmic discretization. We also find that with the same chain length N , the logarithmically discretized Hamiltonian results in a value for α c that is generally lower and hence closer to the limiting value for N → ∞ found via extrapolation than in the linearly discretized case. However, for larger values of N , the simulation takes considerably longer to converge near the phase transition due to the comparatively larger Hilbert space that is populated in this scheme. Hence, using the same computational time as for the linear discretization, we did not acquire reliable data points for N > 25. The dotted magenta lines for s = 0.5 in figures 3, 5 and 7 also show results using the logarithmic discretization with N = 50, using exactly the same computational time as the results for linear discretization. The results are almost identical to those obtained via linear discretization, except just above the transition at around 0.09 α 0.11 where the convergence issues come into play.

Conclusion
By modifying the traditional MPS representation, we were able to avoid the explicit local Hilbert space truncation that leads to the failure of many numerical methods in a regime of high field mode excitation. Instead we introduce a soft cut-off which gives higher weight to lower population numbers, but does not directly truncate the Hilbert space at any dimension. Using this modified state representation combined with a method of energy minimization, we were able to give a detailed study of the ground state properties of the sub-ohmic SBM. Our findings are in good agreement with previous numerical and analytical results, but extend these to new regimes of the SBM, particularly the region close to and above the phase transition. This regime poses considerable challenges to numerical investigations at present, since methods such as DMRG fail to produce reliable results due to the rapid increase of the Hilbert space dimension of the environmental modes. In addition, our method allowed us to give an analysis of the chain properties such as mode excitation and entanglement for specific sites along the chain. It also has the advantage of being comparatively easy to implement numerically. Hence the method provides a promising new tool to investigate the localized phase of the SBM near and above the transition and to test the current analytical results such as the mean-field type approaches in this regime. The results are particularly remarkable considering the brute-force approach (an inbuilt MATLAB function) used for the minimization. We believe that further study of the method will provide more specialized techniques, thus speeding up convergence and allowing for efficient simulations with larger matrix dimension χ. Some first applications to other models such as coupled harmonic systems [20] also show promising results but require further study to confirm a general applicability of the method.