We consider now the anharmonic oscillator problem. Effectively we are then dealing with a $1\times 1$ matrix $\phi$. The theory is given by the action \begin{eqnarray} S=\int_{0}^{\beta} dt \big[\frac{1}{2}\dot{\phi}^2+V(\phi)\big]~,~V=\frac{1}{2}\mu^2\phi^2+\lambda \phi^4. \end{eqnarray} The lattice action reads \begin{eqnarray} S[\phi_n]=\frac{1}{a}\sum_{n=1}^{\Lambda}(\phi_n^2-\phi_{n+1}\phi_n)+a\sum_{n=1}^{\Lambda}(\frac{1}{2}\mu^2\phi_n^2+\lambda \phi_n^4). \end{eqnarray} The Euclidean path integral (partition function) is given explicitly by \begin{eqnarray} Z=\int \prod_{n=1}^{\Lambda}d\phi_n\exp(-S[\phi_n]). \end{eqnarray} This lattice partition function was studied analytically and numerically in the beautiful work \cite{Creutz:1980gp}. For the harmonic oscillator ($\lambda=0$) we find, using the transfer matrix technique, the closed formula \begin{eqnarray} Z={\rm Tr}T^{\Lambda}={\rm Tr}\sqrt{2\pi a}R^{(a^{\dagger}a+\frac{1}{2})}=(2\pi a R)^{\Lambda/2}\frac{1}{1-R^{\Lambda}}. \end{eqnarray} The constant $R$ is given in terms of the lattice spacing $a$ and the mass $\mu^2$ by the expression \begin{eqnarray} R=\frac{\lambda_n}{\lambda_{n-1}}=1+\frac{a^2\mu^2}{2}-a\mu(1+\frac{a^2\mu^2}{4})^{1/2}~,~T|n\rangle=\lambda_n|n\rangle. \end{eqnarray} The ground-state energy $E_0$ is given, using the virial theorem $\langle v_i^2\rangle=\langle \phi.V^{\prime}(\phi)\rangle$, by the formula \begin{eqnarray} E_0&=&\langle \bigg(\frac{1}{2}\phi.V^{\prime}(\phi)+V(\phi)\bigg)\rangle=\mu^2\langle \phi^2\rangle+3\lambda \langle \phi^4\rangle. \end{eqnarray} The virial theorem allows us to define the mean square velocity $\langle v_i^2\rangle$ which otherwise does not really exist since the paths of the quantum particle are irregular (non-differentiable). Another definition of the mean square velocity $\langle v_i^2\rangle$ can be given using a split-point method following Feynman and Hibbs \cite{FeynmanHibbs}. In both cases the expectation value of the velocity-dependent part of the action, which otherwise diverges as $1/a$ when $a\longrightarrow 0$, can be defined. In fact, we have using the virial theorem $E_0=\langle S\rangle$.
The mean squared position $\langle \phi^2\rangle$ for the harmonic oscillator ($\lambda=0$) is found to be given by \begin{eqnarray} \langle\phi^2\rangle=\frac{1}{2\omega}\frac{1+R^{\Lambda}}{1-R^{\Lambda}}~,~\omega^2=\mu^2(1+\frac{a^2\mu^2}{4}). \end{eqnarray} In fact, the $2-$point correlation function (or propagator) for the harmonic oscillator ($\lambda=0$) is found to be given by (with $t=(n-1)a$ and $\Delta t=ma$) \begin{eqnarray} \langle\phi(t)\phi(t+\Delta t)\rangle&=&\langle \phi_{n}\phi_{n+m}\rangle\nonumber\\ &=&\frac{1}{Z}{\rm Tr}\phi T^m\phi T^{\Lambda-m}\nonumber\\ &=&\frac{1}{2\omega}\frac{R^m+R^{\Lambda-m}}{1-R^{\Lambda}}. \end{eqnarray} The case $m=0$ corresponds to the mean squared position $\langle \phi^2\rangle$. The energy $E_1$ of the first exited state (or more precisely the mass gap $E_1-E_0$) is given in terms of the $2-$point correlation function $ \langle\phi(t)\phi(t+\Delta t)\rangle$ by the following formula \begin{eqnarray} E_1-E_0=-\frac{1}{\Delta t}\ln\frac{\langle \phi(0)\phi(t+\Delta t)\rangle}{\langle \phi(0)\phi(t)\rangle}~,~t\longrightarrow\infty. \end{eqnarray} In the continuum limit $a\longrightarrow 0$ we have the behavior $R^m\longrightarrow \exp(-\mu m a(1-a^2\mu^2/24+...))$ and thus one must have \begin{eqnarray} \langle\phi^2\rangle=\frac{1}{2\mu}(1-\frac{a^2\mu^2}{8}+...)\frac{1+e^{-\mu\beta(1-\frac{a^2\mu^2}{24}+...)}}{1-e^{-\mu\beta(1-\frac{a^2\mu^2}{24}+...)}}=\frac{1}{2\mu}~,~\beta\longrightarrow\infty. \end{eqnarray} \begin{eqnarray} \langle\phi(t)\phi(t+\Delta t)\rangle &=&\frac{1}{2\mu}\frac{e^{-\mu \Delta t}+e^{-\mu(\beta-\Delta t)}}{1-e^{-\mu\beta}}=\frac{e^{-\mu\Delta t}}{2\mu}~,~\beta\longrightarrow\infty. \end{eqnarray} \begin{eqnarray} E_1-E_0=\mu. \end{eqnarray} The ground-state wave function of the the harmonic oscillator ($\lambda=0$) on the lattice ($a\neq 0$) is given explicitly by \begin{eqnarray} \psi_0(\phi)\equiv \langle\phi|0\rangle=(\frac{\omega}{\pi})^{1/4}\exp(-\frac{1}{2}\omega \phi^2). \end{eqnarray} Numerically we will also follow \cite{Creutz:1980gp} in studying the anharmonic oscillator on the lattice. Thus, we will choose the lattice spacing $a$ small enough to approximate the continuum limit and the inverse temperature $\beta$ large enough to isolate the ground-state physics of our model. Thus, we choose $a/\beta_E\sim 1/10~{\rm to}~1/20$ and $\beta/\beta_E\sim 3~{\rm to}~ 10$ where $\beta_E=2\pi\hbar/E_0$ is the timescale of the oscillator. We will still work with periodic boundary condition. The ground-state wave function is constructed numerically as the probability distribution (or histogram) of the lattice field $\phi_n$. Indeed, by using Feynman's path integral we know that the probability for finding the particle between $\phi-\Delta\phi$ and $\phi+\Delta \phi$ is the time average over transition amplitudes, viz \begin{eqnarray} P(\phi;T)=\frac{1}{T}\int_0^T dt^{\prime}\int_{\phi-\Delta \phi}^{\phi+\Delta \phi}d\phi^{\prime}\frac{Z(\phi_f,T;\phi^{\prime},t^{\prime})Z(\phi^{\prime},t^{\prime};\phi_i,0)}{Z(\phi_f,T;\phi_i,0)}. \end{eqnarray} This computes the number of times the particle passes through $\phi$ with error $\Delta \phi$. In the limit $\Delta\phi\longrightarrow 0$ and $T\longrightarrow \infty$ we obtain \cite{Creutz:1980gp} \begin{eqnarray} \frac{P(\phi;T)}{\Delta \phi}=|\psi_0(\phi)|^2+O(\frac{1}{T(E_1-E_0)}). \end{eqnarray} Thus, for $T\gg 1/(E_1-E_0)$ we can isolate the ground-state wave function. In other words, the histogram is given exactly by $|\psi_0(\phi)|^2$, viz \begin{eqnarray} \frac{P(\phi;T)}{\Delta \phi}=|\psi_0(\phi)|^2=\frac{1}{T_{\rm MC}\Delta\phi}\sum_{i=1}^{T_{\rm MC}}\theta(\Delta\phi-|\phi_n^{(i)}-\phi|). \end{eqnarray} The step function is equal $1$ only if $\phi-\Delta\phi\leq \phi_n\leq \phi+\Delta\phi$. The $2$-nd moment (the mean squared position) $\langle \phi^2\rangle$, the $4$-th moment $\langle \phi^2\rangle$ and the $2-$point function (propagator) $\langle \phi(0)\phi(t)\rangle$ are numerically defined by \begin{eqnarray} \langle \phi^2\rangle=\langle \frac{1}{\Lambda}\sum_{i=1}^{\Lambda}\phi_{i}\phi_{i}\rangle. \end{eqnarray} \begin{eqnarray} \langle \phi^4\rangle=\langle \frac{1}{\Lambda}\sum_{i=1}^{\Lambda}\phi_{i}\phi_{i}\phi_{i}\phi_{i}\rangle. \end{eqnarray} \begin{eqnarray} \langle \phi(0)\phi(t)\rangle=\langle \frac{1}{\Lambda}\sum_{i=1}^{\Lambda}\phi_{i}\phi_{n-1+i}\rangle. \end{eqnarray} The analytical results for the anharmonic oscillator with potential $V=\lambda (\phi^2-f^2)$, i.e. with $\mu^2=-4\lambda f^2$ can be found in \cite{Blankenbecler:1979pa}. The analytical solution of both the harmonic and the anharmonic oscillators will be used to calibrate the numerical simulations. Furthermore, we will study the anharmonic and the harmonic oscillators with both the Metropolis and the hybrid Monte Carlo algorithms which will allow us to compare the tunning and performance of the much more complicated hybrid Monte Carlo algorithm against those of the simpler Metropolis algorithm and compare both algorithms against the theory. Naturally, the numerical study of this model in \cite{Creutz:1980gp} was conducted with the Metropolis algorithm. Indeed, there is no need in this case to the hybrid Monte Carlo algorithm which is more suited to highly non-local theories such as matrix and supersymmetric models. Hence, our interest in the hybrid Monte Carlo algorithm here is essentially for tunning purposes. The Metropolis and hybrid Monte Carlo algorithms as applied to the anharmonic oscillator are as follows:
- Metropolis algorithm:
- The variation of the action $S[\phi_n]$ under the change $\phi_n\longrightarrow\phi_n^{\prime}=\phi_n+\epsilon$ is given by \begin{eqnarray} \Delta S(\phi_n;n,\epsilon)&=&\frac{\epsilon}{a}(\epsilon+2\phi_n-\phi_{n-1}-\phi_{n+1})+\frac{a\mu^2}{2}(\epsilon^2+2\epsilon\phi_n)\nonumber\\ &+&a\lambda(\epsilon^2+2\epsilon\phi_n)(\epsilon^2+2\epsilon\phi_n+2\phi_n^2). \end{eqnarray}
- The Metropolis step: We accept the configuration $\phi_n^{\prime}=\phi_n+\epsilon$ with the Boltzmann weight, viz \begin{eqnarray} P(\phi_n\longrightarrow \phi_n^{\prime}=\phi_n+\epsilon)={\rm min}(1,\exp(-\Delta S(\phi_n;n,\epsilon))). \end{eqnarray}
- Hybrid Monte Carlo algorithm:
- The Hamiltonian and the force in a pseudo-time $\tau$ are given respectively by \begin{eqnarray} H[p_n,\phi_n]=\frac{1}{2}\sum_{n=1}^{\Lambda}p_n^2+S[\phi_n]. \end{eqnarray} \begin{eqnarray} F_n=-\frac{\partial S}{\partial \phi_n}=\frac{1}{a}(\phi_{n+1}+\phi_{n-1}-2\phi_n)-a(\mu^2\phi_n+4\lambda\phi_n^3). \end{eqnarray}
- Hamilton's equations of motion are solved, starting from the initial configuration $\phi(0)=\phi$, by the molecular dynamics (leap-frog) algorithm given by \begin{eqnarray} &&F_n(\tau)=F_n(\phi_n(\tau))\nonumber\\ &&p_n(\tau+\frac{\delta\tau}{2})=p_n(\tau)+\frac{\delta\tau}{2}F_n(\tau)\nonumber\\ &&\phi_n(\tau+\delta\tau)=\phi_n(\tau)+\delta\tau p_n(\tau+\frac{\delta\tau}{2})\nonumber\\ &&F_n(\tau+\delta\tau)=F_n(\phi_n(\tau+\delta\tau))\nonumber\\ &&p_n(\tau+\delta\tau)=p_n(\tau+\frac{\delta\tau}{2})+\frac{\delta\tau}{2}F_n(\tau+\delta\tau). \end{eqnarray}
- The Metropolis step: We accept the configuration $\phi^{\prime}=\phi(T)$ where $T=\delta\tau.N_{\tau}$ with the Boltzmann weight, viz \begin{eqnarray} P(\phi\longrightarrow \phi^{\prime})={\rm min}(1,\exp(-\Delta H))~,~\Delta H=H(\phi^{\prime})-H(\phi). \end{eqnarray}
- The heat bath: In order to be able access the full phase space (ergodicity) we refresh the momentum $p_n$ using the Gaussian distribution $\exp(-p_n^2/2)$. Thus, if $v_1$ and $v_2$ are two uniform random numbers we write \begin{eqnarray} p_n=\sqrt{-2\ln v_1}\cos(2\pi v_2). \end{eqnarray}
No comments:
Post a Comment