# Monte Carlo Hamiltonian

## 1 Motivation

The motivation for constructing a Monte Carlo Hamiltonian comes from different directions.   
(i) The renormalization group à la Kadanoff-Wilson  aims to construct a renormalized Hamiltonian, which describes physics at a critical point, based on the assumption of scale invariance. Such Hamiltonian is supposed to have much less degrees of freedom than the original Hamiltonian. A recent further development of those ideas is White’s density matrix renormalization group technique . Our goal is similar to the above in the sense that we aim at an effective Hamiltonian, which has ”less” degrees of freedom than the ”original” Hamiltonian. But it differs in describing physics in the low energy domain instead of doing so at the critical point.   
(ii) When one tries to solve field theory in the Hamiltonian formulation, the standard way to proceed is by constructing a Fock space, parametrized by some high momentum cut-off and some occupation number cut-off (Tamm-Dancoff approximation). When increasing those parameters, which means increasing the upper bound of the energy, then typically the density of states increases in an exponential manner, which renders the system beyond any control. In contrast to that, the Monte Carlo Hamiltonian is governed by a ”small” number of low-energy degrees of freedom and the spectral density decreases with increasing energy.   
(iii) The enormous success of lattice field theory over the last quarter of the century is certainly due to the fact that the Monte Carlo method with importance sampling is an excellent technique to solve high dimensional (and even ”infinite” dimensional) integrals. Conventionally, one computes a transition amplitude of an operator and evaluates it numerically via Monte Carlo (e.g. Metropolis algorithm ),

$`<O>`$ $`=`$ $`{\displaystyle \frac{[dx]O[x]\mathrm{exp}(\frac{1}{\mathrm{}}S[x])}{[dx]\mathrm{exp}(\frac{1}{\mathrm{}}S[x])}}`$ (1)   
$``$ $`{\displaystyle \frac{1}{N_c}}{\displaystyle \underset{C}{}}O[C].`$   
Here $`C`$ stands for a path configuration drawn from the distribution $`P[x]=\frac{1}{Z}\mathrm{exp}(\frac{1}{\mathrm{}}S[x])`$. The vitue of the Monte Carlo method lies in the property of yielding very good numerical results. E.g., solving a field theory model on a lattice of size $`20^4`$ and measuring the observable from a number of configurations $`N_c`$ in the order of a few hundred typically yields results with statistical errors in the order of a few percent. In this way it has been possible to determine low lying baryon and meson masses quite precisely .

On the other hand, one can express a transition amplitude in imaginary time via the Hamiltonian

$`<x_{fi},T|x_{in},0>`$ $`=`$ $`<x_{fi}|e^{HT/\mathrm{}}|x_{in}>`$ (2)   
$`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}<x_{fi}|E_n>e^{E_nT/\mathrm{}}<E_n|x_{in}>`$   
$``$ $`<x_{fi}|e^{H_{ef}T/\mathrm{}}|x_{in}>`$   
$`=`$ $`{\displaystyle \underset{\nu =1}{\overset{N}{}}}<x_{fi}|E_\nu ^{eff}>e^{E_\nu ^{eff}T/\mathrm{}}<E_\nu ^{eff}|x_{in}>.`$   
In the last two lines we have approximated the Hamiltonian $`H`$ by an effective Hamiltonian $`H_{eff}`$, which has less degrees of freedom, e.g., it has only $`N`$ eigenstates. The idea of the Monte Carlo Hamiltonian is that an effective Hamiltonian can be found via use of Monte Carlo, such that transition amplitudes become a finite sum over $`N`$ eigenstates, where $`N`$ is in the order of magnitude of $`N_c`$, i.e. the number of equilibrium configurations, sufficient to closely approximate the path integral of Eq.(1).

One might ask: What is the virtue of such a Hamiltonian? A list of physics problems, where progress has been slow with conventional methods including standard lattice techniques, and where such a Hamiltonian might bring progress are the following topics:   
\- Non-perturbative computation of cross sections and decay amplitudes in many-body systems.   
\- Low-lying but excited states of the hadronic spectrum and the related question of quantum chaos in such a system.   
\- Hadron wave functions and the related question of hadron structure functions, in particular for small $`x_B`$ and $`Q^2`$. The Hamiltonian formulation is suited to compute wave functions, which is quite difficult in the Lagrangian lattice formulation.   
\- Finite temperature and in particular finite density in baryonic matter. This is crucial for the quark-gluon plasma phase transition, the physics of neutron stars and cosmology. The Hamiltonian formulation is suited to compute the mean value of the energy (average energy). This is difficult to compute in the Lagrangian lattice formulation where one usually computes the expectation value of the action. Finite density $`QED`$ and $`QCD`$ in the Lagrangian lattice formulation is hampered by the notorious complex action problem.   
\- Atomic physics: study of spectra and the question of quantum chaos.   
\- Condensed matter physics: study of spin systems (computation of dynamical structure factors), and high $`T_c`$ superconductivity models (search for electron pair attraction at very small energy). In the following we will outline how to construct such a Hamiltonian.

## 2 Construction of $`H_{eff}`$

In contrast to the statistical mechanics concept of the transfer matrix, which describes the time-evolution (we consider imaginary time) when advancing the system by a small discrete time step $`\mathrm{\Delta }t=a_0`$ and from which one can infer the Hamiltonian $`(a_00)`$, here we consider transition amplitudes $`<\psi |e^{HT/\mathrm{}}|\varphi >`$ corresponding to a finite, long time $`T`$ ($`T>>a_0`$), for the purpose to reconstruct the spectrum in some finite low energy domain. Let us start from a complete orthonormal basis of Hilbert states $`|e_i>,i=1,2,3\mathrm{}`$ and consider the matrix elements for a given fixed $`N`$

$$M_{ij}(T)=<e_i|e^{HT/\mathrm{}}|e_j>,i,j1,\mathrm{},N.$$
(3)   
Under the assumption that $`H`$ is Hermitian, $`M(T)`$ is a positive, Hermitian matrix. Elementary linear algebra implies that there is a unitary matrix $`U`$ and a real, diagonal matrix $`D`$ such that

$$M(T)=U^{}D(T)U.$$
(4)   
On the other hand, projecting $`H`$ onto the the subspace $`S_N`$ generated by the first $`N`$ states of the basis $`|e_i>`$, and using the eigenrepresentation of such Hamiltonian, one has

$$M_{ij}(T)=\underset{k=1}{\overset{N}{}}<e_i|E_k^{eff}>e^{E_k^{eff}T/\mathrm{}}<E_k^{eff}|e_j>,$$
(5)   
and we can identify

$$U_{ik}^{}=<e_i|E_k^{eff}>,D_k(T)=e^{E_k^{eff}T/\mathrm{}}.$$
(6)   
Let us assume for the moment that the matrix elements $`M_{ij}(T),i,j=1,\mathrm{},N`$ would be known. Then algebraic diagonalization of the matrix $`M(T)`$ yields eigenvalues $`D_k(T),k=1,\mathrm{},N`$, which by Eq.(6) gives the spectrum of energies,

$$E_k^{eff}=\frac{\mathrm{}}{T}\mathrm{ln}D_k(T),k=1,\mathrm{},N.$$
(7)   
The corresponding k-th eigenvector can be identified with the k-th column of the matrix $`U_{ik}^{}`$. From Eq.(6) we then know the wave function of the k-th eigenstate expressed in terms of the basis $`|e_i>`$. Thus starting from the matrix elements $`M_{ij}(T)`$ we have explicitly constructed an effective Hamiltonian

$$H_{eff}=\underset{k=1}{\overset{N}{}}|E_k^{eff}>E_k^{eff}<E_k^{eff}|.$$
(8)   
## 3 Computation of matrix elements by Monte Carlo

We suggest to compute the matrix elements $`M_{ij}(T)`$ directly from the action via Monte Carlo with importance sampling. For the sake of simplicity, let us consider $`D=1`$. We choose basis states $`|e_i>`$ in position space by introducing a lattice with nodes $`x_i`$ and define $`e_i(x)`$ (unnormalized) by $`e_i(x)=1`$ if $`x_ixx_{i+1}`$, zero else. $`\mathrm{\Delta }x_i=x_{i+1}x_i`$. In numerical calculations we have used a regular lattice, $`\mathrm{\Delta }x_i=\text{const}`$. The matrix elements read

$`M_{ij}(T)`$ $`=`$ $`{\displaystyle _{x_i}^{x_{i+1}}}𝑑y{\displaystyle _{x_j}^{x_{j+1}}}𝑑z<y,T|z,0>`$ (9)   
$`=`$ $`{\displaystyle _{x_i}^{x_{i+1}}}𝑑y{\displaystyle _{x_j}^{x_{j+1}}}𝑑z{\displaystyle [dx]\mathrm{exp}[S[x]/\mathrm{}]}|_{z,0}^{y,T}.`$   
Here $`S`$ denotes the Euclidean action for a given path $`C`$,

$$S[C]=_0^T𝑑t\frac{1}{2}m\dot{x}^2+V(x)|_C.$$
(10)   
The Monte Carlo method with importance sampling is suited and conventionally applied to estimate a ratio of integrals, like in Eq.(1). Here we suggest to estimate the matrix elements $`M_{ij}`$ by splitting the action

$$S=S_0+S_V_0^T𝑑t\frac{1}{2}m\dot{x}^2+_0^T𝑑tV(x),$$
(11)   
and to express $`M_{ij}`$ as

$$M_{ij}(T)=M_{ij}^{(0)}(T)\frac{_{x_i}^{x_{i+1}}𝑑y_{x_j}^{x_{j+1}}𝑑z[dx]\mathrm{exp}[S_V[x]/\mathrm{}]\mathrm{exp}[S_0[x]/\mathrm{}]|_{z,0}^{y,T}}{_{x_i}^{x_{i+1}}𝑑y_{x_j}^{x_{j+1}}𝑑z[dx]\mathrm{exp}[S_0[x]/\mathrm{}]|_{z,0}^{y,T}},$$
(12)   
where $`O\mathrm{exp}[S_V/\mathrm{}]`$ is treated as an observable. The ratio can be treated by standard Monte Carlo methods with importance sampling. The matrix elements $`M_{ij}^{(0)}`$, corresponding to the free action $`S_0`$, are almost known analytically,

$$M_{ij}^{(0)}(T)=_{x_i}^{x_{i+1}}𝑑y_{x_j}^{x_{j+1}}𝑑z\sqrt{\frac{m}{2\pi \mathrm{}T}}\mathrm{exp}\left[\frac{m}{2\mathrm{}T}(yz)^2\right].$$
(13)   
## 4 Test of $`H_{eff}`$

### 4.1 Free system

In order to test the effective Hamiltonian, we have computed the energy spectrum, its wave functions and thermodynamic observables like the average energy $`U`$ and the specific heat $`C`$ as well as the partition function $`Z`$. They are defined by

$`Z(\beta )`$ $`=`$ $`Tr[e^{\beta H}],`$   
$`U(\beta )`$ $`=`$ $`{\displaystyle \frac{1}{Z}}Tr[He^{\beta H}]={\displaystyle \frac{\mathrm{log}Z}{\beta }},`$   
$`C(\beta )`$ $`=`$ $`{\displaystyle \frac{F}{𝒯}}=k_B\beta ^2{\displaystyle \frac{^2\mathrm{log}Z}{\beta ^2}},`$ (14)   
where $`\beta =(k_B𝒯)^1`$, $`𝒯`$ is the temperature, and we identify $`\beta `$ with the imaginary time $`T`$ by $`\beta =T/\mathrm{}`$. For the free system one obtains the following analytical expressions for $`Z`$, $`U`$ and $`C`$,

$`Z(\beta )`$ $`=`$ $`\sqrt{{\displaystyle \frac{m}{2\pi \mathrm{}^2\beta }}}I,I={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x(\text{being infinite}),`$   
$`U(\beta )`$ $`=`$ $`{\displaystyle \frac{1}{2\beta }}={\displaystyle \frac{1}{2}}k_B𝒯,`$   
$`C(\beta )`$ $`=`$ $`{\displaystyle \frac{1}{2}}k_B.`$ (15)   
Note that $`U(\beta )_\beta \mathrm{}0`$, i.e. it tends to the ground state energy of the free system (Feynman-Kac formula).

The partition function corresponding to the effective Hamiltonian is obtained via its spectrum,

$$Z_{eff}(\beta )=Tr[e^{\beta H_{eff}}]=\underset{k=1}{\overset{N}{}}e^{\beta E_k^{eff}}.$$
(16)   
Via Eq.(14) one obtains the corresponding average energy $`U_{eff}`$ and the specific heat $`C_{eff}`$. One should keep in mind that $`H_{eff}`$ has been constructed for a specific value of the time parameter, $`T=1`$ corresponding to the temperature $`𝒯=1`$ (we use $`\mathrm{}=k_B=1`$). Fig. shows a plot of the average energy, comparing the exact result with that from the effective Hamiltonian. One observes that the agreement is better where $`𝒯0`$, i.e. in the low energy regime. A similar behavior is found for the specific heat, shown in Fig..

### 4.2 Harmonic oscillator

The Euclidean action of the harmonic oscillator is given by

$$S=_o^T𝑑t\frac{1}{2}m\dot{x}^2+\frac{1}{2}m\omega ^2x^2.$$
(17)   
The energy spectrum is

$$E_n=\mathrm{}\omega (n+1/2),n=0,1,2,\mathrm{}$$
(18)   
A comparison of the spectrum of the effective Hamiltonian with the exact one is shown in Tab.. As can be seen, the error is small in the low energy domain. A more stringent test is that of the wave functions. Fig. shows a comparison for the wave functions of the three lowest states. We have also verified the low energy behavior of the effective Hamiltonian by computing the partition function, average energy and specific heat as a function of temperature. For the harmonic oscillator those are analytically known,

$`Z(\beta )`$ $`=`$ $`\left[2\mathrm{sinh}(\beta \mathrm{}\omega /2)\right]^1,`$   
$`U(\beta )`$ $`=`$ $`{\displaystyle \frac{\mathrm{}\omega }{2}}\text{ctgh}(\beta \mathrm{}\omega /2),`$   
$`C(\beta )`$ $`=`$ $`k_B\left[{\displaystyle \frac{\beta \mathrm{}\omega /2}{\mathrm{sinh}(\beta \mathrm{}\omega /2)}}\right]^2.`$ (19)   
In the limit $`\beta \mathrm{}`$ the average energy tends to the ground state energy, $`U\mathrm{}\omega /2`$ (Feynman-Kac formula). A plot of the average energy and the specific heat is shown in Figs.. The effective Hamiltonian, constructed at $`T_c=\beta _c=𝒯_c=1`$, describes well thermodynamic observables in the range $`\beta _c\beta `$ (it works also for $`\beta >10`$, not shown in the figure). However, it breaks down for $`\beta <\beta _c`$, i.e. $`𝒯>𝒯_c`$. This is due to the small dimension $`N=20`$ of the matrix. Agreement in a larger $`\beta `$-region, i.e. lowering $`\beta _c`$ can be obtained by increasing $`N`$. This can be seen, e.g. for the free system in Fig., where $`N=200`$ and $`\beta _c<0.1`$.

### 4.3 Other local potentials

We have tested the effective Hamiltonian for other local potentials. For example,

$$V(x)=V_0\text{sech}^2(x/d)$$
(20)   
is a potential having a minimum $`V_0`$ at $`x=0`$ and rising asymptotically to zero at $`x=\pm \mathrm{}`$. It generates a bound state spectrum, being analytically known . It is given by

$`E_n`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2md^2}}\lambda _n,`$   
$`\lambda _n`$ $`=`$ $`\left[(n+1/2)\sqrt{Q+1/4}\right]^2,`$   
$`Q`$ $`=`$ $`{\displaystyle \frac{2md^2}{\mathrm{}^2}}V_0,`$   
$`n`$ $`=`$ $`0,1,2,\mathrm{},n_{max}<\sqrt{Q+1/4}1/2.`$ (21)   
The results are shown in Tabs..

## 5 Conclusion

We have proposed to construct an effective low-energy Hamiltonian from the action via use of the Monte Carlo method. We have shown that the method works for a number of systems in $`1D`$ quantum mechanics, by computing the spectrum, wave functions and thermodynamical observables.   
\- We have not given an error estimate of the statistical errors. The reason is that the statistical error of the matrix elements can be estimated easily, however, to get from that an error estimate of the energy spectrum is difficult. We defer that to a later study.   
\- We have not discussed an application to a field theory or a many-body system, although this is the area where the method should prove to be most useful. The reason is that this requires a new step, namely a stochastic (Monte Carlo) selection from the set of basis functions. This is presently under investigation.   
\- In our opinion an effective low-energy Hamiltonian will be very useful in condensed matter physics, atomic physics, nuclear physics, and high energy particle physics.

Acknowledgments   
H.K. would like to acknowledge helpful discussions with M. Creutz, B. Berg, W. Janke. and P. Amiot. H.K. and K.J.M.M. are grateful for support by NSERC Canada. X.Q.L. is supported by the National Natural Science Fund for Distinguished Young Scholars, supplemented by the National Natural Science Foundation of China, fund for international cooperation and exchange, the Ministry of Education, and the Hong Kong Foundation of the Zhongshan University Advanced Research Center.

Figure Caption

Average energy of the free system. Solid line and diamonds, respectively, represent the exact analytical result, and that from the exact matrix elements for $`\mathrm{\Delta }x=0.5`$ and $`N=100`$. The cross at $`\beta =0.1`$ corresponds to $`\mathrm{\Delta }x=0.2`$ and $`N=200`$.

Specific heat over $`k_B`$ of the free system. Symbols as in Fig..

Wave function of the harmonic oscillator, (a) ground state, (b) first excited state, (c) second excited state. Solid line, diamonds and crosses, respectively, represent the exact analytical result, that from the exact matrix elements, and that from Monte Carlo simulation.

Average energy of the harmonic oscillator. Symbols as in Fig..

Specific heat over $`k_B`$ of the harmonic oscillator. Symbols as in Fig..

Eigenvalues of the harmonic oscillator. $`m=1`$, $`\mathrm{}=1`$, $`\omega =0.6`$, $`\mathrm{\Delta }x=1`$, $`N=20`$. $`E_n^{exact}`$, $`E_n^{e.m.}`$ and $`E_n^{m.c.}`$, respectively, represent the exact analytical result, that from the exact matrix elements, and that from Monte Carlo simulation.

Bound state spectrum for potential given by Eq.(20). (a) For $`m=1.0,\mathrm{}=1.0,T=1.0,V_0=1.0,d=1.0,Q=2,\mathrm{\Delta }x=1.0,N=10`$, there is only one bound state $`n_{max}<1`$. This is confirmed by the Monte Carlo data. (b) For $`m=1.0,\mathrm{}=1.0,T=1.0,V_0=1.0,d=2.0,Q=8,\mathrm{\Delta }x=1.0,N=20`$, there are three bound states $`n_{max}<3`$. This is confirmed by the Monte Carlo data.