Matrix product state representation without explicit local Hilbert space truncation with applications to the sub-ohmic spin-boson model

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-field-like 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–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].

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,\ldots ,i_N \in \mathbb {N}^0$) can be described by

$$\begin{equation} |{\psi}\rangle = \sum_{k, i_1,\ldots, i_N} c_{k, i_1,\ldots ,i_N}|{k, i_1,\ldots ,i_N}\rangle. \end{equation} \tag{ 1 }$$

Alternatively, it is possible to express the coefficients c_k,i1,...,iN as products of χ × χ matrices, where the M-level system as well as each mode have a unique set of matrices associated with them

$$\begin{eqnarray} |{\psi}\rangle &=& \sum_{k, i_1,\ldots , i_N} \mathrm{tr}\left[S^{(k)}\prod_{m=1}^N B_m^{(i_m)}\right] \prod_{m=1}^N(b_m^{\dagger})^{i_m}|{k, 0}\rangle \nonumber \\ &=& \sum_{k, i_1,\ldots , i_N} \mathrm{tr}\left[S^{(k)}\prod_{m=1}^N B_m^{(i_m)}\sqrt{i_m!}\right] |{k, i_1,\ldots ,i_N}\rangle, \end{eqnarray} \tag{ 2 }$$

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–13]. In this form the state of the M-level system is represented by the M matrices S⁽¹⁾,...,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

$$\begin{equation} B_k^{(i_m)} = \frac{X_m^{i_m}}{i_m!} \end{equation} \tag{ 3 }$$

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.⁴ The additional factor of (i_m!)⁻¹ 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

$$\begin{equation} |{\psi}\rangle = \sum_{k=1}^M\sum_{i_1,\ldots , i_N = 0}^{\infty} \mathrm{tr}\left[S^{(k)}\prod_{m=1}^N \frac{X_m^{i_m}}{\sqrt{i_m!}}\right] |{k, i_1,\ldots ,i_N}\rangle. \end{equation} \tag{ 4 }$$

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. Introducing the new notation

$$\begin{equation} \Gamma_a^b \equiv \prod_{m=a}^b \mathrm{e}^{\bar{X}_m \otimes X_m} \end{equation} \tag{ 5 }$$

and

$$\begin{equation} \Xi \equiv \sum_{k=1}^M \bar{S}^{(k)} \otimes S^{(k)}, \end{equation} \tag{ 6 }$$

we find for the norm

$$\begin{equation} \langle{\psi | \psi}\rangle = \mathrm{tr}[\Xi \Gamma_1^N]. \end{equation} \tag{ 7 }$$

The exponentiation of the X_m matrices in combination with the deliberate choice of the (i_m!)⁻¹ 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

$$\begin{equation} \langle{\psi | b_k^{\dagger} b_k | \psi}\rangle = \mathrm{tr}[\Xi \Gamma_1^k (\bar{X}_k \otimes X_k) \Gamma_{k+1}^N]. \end{equation} \tag{ 8 }$$

In general, we find that all quantities of interest are simple traces over products of χ² × χ² 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

$$\begin{equation} E = \frac{\langle{\psi | H | \psi}\rangle}{\langle{\psi | \psi}\rangle} \end{equation} \tag{ 9 }$$

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.

The Hamiltonian of the (unbiased) SBM is given by (ℏ = 1)

$$\begin{equation} H_{\mathrm{SB}} = -\frac{1}{2}\Delta\sigma_x + \frac{1}{2}\sigma_z\sum_l g_l(a_l + a_l^{\dagger}) + \sum_l\omega_l a_l^{\dagger} a_l, \end{equation} \tag{ 10 }$$

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(\omega )=\pi \sum _lg_l^2\delta (\omega -\omega _l)$. Here we consider a spectral function of the form

$$\begin{equation} J(\omega) = 2\pi\alpha\omega_{\mathrm{c}}^{1-s}\omega^s\Theta(\omega_{\mathrm{c}}-\omega) \end{equation} \tag{ 11 }$$

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

$$\begin{eqnarray} &&H = -\frac{1}{2}\Delta\sigma_x + c_0\sigma_z(b_0 + b_0^{\dagger})+\sum_{n=1}^{\infty}\omega_n b_n^{\dagger}b_n + t_n(b_{n+1}^{\dagger}b_n + b_n^{\dagger}b_{n+1}), \end{eqnarray} \tag{ 12 }$$

where the coupling strength between the TLS and the first site is given by

$$\begin{equation} c_0 = \sqrt{\frac{\alpha}{2(s+1)}}\omega_{\mathrm{c}} \end{equation} \tag{ 13 }$$

and the local energies and tunnelling amplitudes of the sites are

$$\begin{equation} \omega_n = \frac{\omega_{\mathrm{c}}}{2}\Bigl(1 + \frac{s^2}{(s+2n)(2+s+2n)}\Bigr) \end{equation} \tag{ 14 }$$

and

$$\begin{equation} t_n = \frac{\omega_{\mathrm{c}}(1+n)(1+s+n)}{(s+2+2n)(3+s+2n)}\sqrt{\frac{3+s+2n}{1+s+2n}}, \end{equation} \tag{ 15 }$$

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

$$\begin{equation} E = \langle{\psi | H | \psi}\rangle \equiv E_{\mathrm{loc}} + E_{\mathrm{int}} + E_{\mathrm{chain}} \end{equation} \tag{ 16 }$$

with the TLS's local energy

$$\begin{equation} E_{\mathrm{loc}} = -\textstyle\frac{1}{2}\Delta\langle{\sigma_x}\rangle \end{equation} \tag{ 17 }$$

the system–chain interaction energy

$$\begin{equation} E_{\mathrm{int}} = c_0\langle{\sigma_z(b_0 + b_0^{\dagger})}\rangle \end{equation} \tag{ 18 }$$

and the chain energy

$$\begin{equation} E_{\mathrm{chain}} = \sum_{n=1}^{\infty}\omega_n \langle{b_n^{\dagger}b_n}\rangle + t_n(\langle{b_{n+1}^{\dagger}b_n}\rangle + \langle{b_n^{\dagger}b_{n+1}}\rangle) \end{equation} \tag{ 19 }$$

subject to the constraint 〈ψ|ψ〉 = 1. In terms of the MPS description we obtain, after truncating the chain length to N sites, the total energy

$$\begin{eqnarray} &&\fl E = -\frac{\Delta}{2} \mathrm{tr}\Bigl[\!\Bigl(\!\bar{S}^{(1)}\otimes S^{(2)} + \bar{S}^{(2)}\otimes S^{(1)} \Bigr)\Gamma_1^N\Bigr] + c_0 \mathrm{tr}\Bigl[\!\Bigl(\!\bar{S}^{(1)}\otimes S^{(1)} - \bar{S}^{(2)}\otimes S^{(2)}\!\Bigr)\Bigl(\!\mathbbm{1}\otimes X_1 + \bar{X}_1 \otimes \mathbbm{1}\! \Bigr)\Gamma_1^N\!\Bigr] \nonumber \\ &&+ \sum_{n=1}^{N-1}\Bigl\{\mathrm{tr}\Bigl[\Xi\Gamma_1^n\bigl(\omega_n(\bar{X}_n\otimes X_n) + t_n(\bar{X}_{n+1}\otimes X_n + \bar{X}_n\otimes X_{n+1})\bigr)\Gamma_{n+1}^N\Bigr]\Bigr\}\nonumber\\ &&+\, \omega_N \mathrm{tr}\Bigl[\Xi\Gamma_1^N\bigl(\bar{X}_N\otimes X_N)\Bigr]. \end{eqnarray} \tag{ 20 }$$

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.

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

$$\begin{equation} \alpha_{\mathrm{c}}^{(\chi)}(N) = a\mathrm{e}^{b/N} \end{equation} \tag{ 21 }$$

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 value

$$\begin{equation} \tilde{\alpha}_{\mathrm{c}} = \frac{\sin(\pi s)\mathrm{e}^{-s/2}}{2\pi(1-s)}\Bigl(\frac{\Delta}{\omega_{\mathrm{c}}}\Bigr)^{1-s}. \end{equation} \tag{ 22 }$$

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.

Zoom In Zoom Out Reset image size

Table 1. Critical coupling strengths α⁽³⁾_c for various values of s, together with the respective critical values $\tilde {\alpha }_{\mathrm {c}}$ derived from the polaron ansatz in [2] and the fractional deviation.

s	α(3)c	$\tilde {\alpha }_{\mathrm {c}}$	$(\alpha ^{(3)}_{\mathrm {c}} - \tilde {\alpha }_{\mathrm {c}})/\tilde {\alpha }_{\mathrm {c}}$
0.1	0.0117±0.0002	0.0065	0.800
0.2	0.0189±0.0003	0.0168	0.125
0.3	0.0315±0.0004	0.0316	−0.003
0.4	0.0485±0.0007	0.0519	−0.066
0.5	0.0749±0.0009	0.0784	−0.045

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 phase transition, but the qualitative behaviour is quite different from χ ⩾ 2. For all χ ⩾ 2 the qualitative features appear to be the same, the only difference being the location of the phase transition as obtained from finite scaling. Via our extrapolation ansatz we find for s = 0.2

$$\begin{eqnarray} &&\begin{array}{l} \alpha^{(2)}_{\rm c} = 0.0218\pm 0.0001, \\ \alpha^{(3)}_{\rm c} = 0.0189\pm 0.0003,\\ \alpha^{(4)}_{\rm c} = 0.0179\pm 0.0005. \end{array} \end{eqnarray} \tag{ 23 }$$

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 α⁽⁴⁾_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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

Mean-field theory predicts a second-order magnetic transition at α_c with

$$\begin{equation} M \propto |\alpha_{\rm c}-\alpha|^{1/2}. \end{equation} \tag{ 24 }$$

Our simulations do indeed reproduce the correct critical mean-field exponent $\frac {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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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^†_nb_n〉 of chain site n rapidly decreases along the chain, it rises considerably in the localized phase. In fact, [2] predicts that 〈b^†_nb_n〉 diverges along the chain for α > α_c. In figure 8 we plot 〈b^†_nb_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^†_nb_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.

Zoom In Zoom Out Reset image size

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 [Λ⁻⁽ⁿ⁺¹⁾,Λ⁻ⁿ], where Λ > 1 is the discretization parameter and $n\in \mathbb {N}^0$. This new Hamiltonian has again the same form equation (12), but with different site frequencies and transition amplitudes

$$\begin{equation} \omega_n = \zeta_s(A_n + C_n) \end{equation} \tag{ 25 }$$

and

$$\begin{equation} t_n = -\zeta_s\frac{N_{n+1}}{N_n}A_n, \end{equation} \tag{ 26 }$$

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.

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.

We acknowledge discussion with Alex W Chin. This work was supported by the Alexander von Humboldt Foundation.