---
abstract: 'Extending the Feynman-Kleinert variational approach, we calculate the temperature-dependent effective classical potential governing the quantum statistics of a hydrogen atom in a uniform magnetic at all temperatures. The zero-temperature limit yields the binding energy of the electron which is quite accurate for all magnetic field strengths and exhibits, in particular, the correct logarithmic growth at large fields.'
address: 'Institut für Theoretische Physik, Freie Universität Berlin, Arnimallee 14, 14195 Berlin'
author:
- 'M. Bachmann, H. Kleinert, and A. Pelster'
title: |
     Variational Approach to Hydrogen Atom in Uniform\
    Magnetic Field of Arbitrary Strength
---

Introduction
============

The quantum statistical and quantum mechanical properties of a hydrogen atom in an external magnetic field are not exactly calculable. Perturbative approaches yield good results only for weak uniform fields as discussed in detail by Le Guillou and Zinn-Justin [@zinn-justin], who interpolated with analytic mapping techniques the ground state energy between weak- and strong-field. Other approaches are based on recursive procedures in higher-order perturbation theory [@cizek1; @cizek2; @gani]. Zero-temperature properties were also investigated with the help of an operator optimization method in a second-quantized variational procedure [@feranshuk]. The behaviour at high uniform fields was inferred from treatments of the one-dimensional hydrogen atom [@landau; @loudon; @haines]. Hydrogen in strong magnetic fields is still a problem under investigation, since its solution is necessary to understand the properties of white dwarfs and neutron stars, as emphasized in Refs. [@cohen; @ivanov; @heyl].

A compact and detailed presentation of the bound states and highly accurate numerically values for the energy levels is given in Ref. [@wunner].

Equations for a first-order variational approach to the ground-state energy of hydrogen in a uniform magnetic field based on the Jensen-Peierls inequality were written down a long time ago [@dev0], but never evaluated. Apparently, they merely served as a preparation for attacking the more complicated problem of a polaron in a magnetic field [@dev0; @dev1; @dev2].

In our approach, we calculate the quantum statistical properties of the system by an extension of variational perturbation theory [@PI]. The crucial quantity is the effective classical potential. In the zero-temperature limit, it yields the ground state energy. Our calculations in a magnetic field require an extension of the formalism in Ref. [@PI] which derives the effective classical potential from the phase space representation of the partition function.

Variational perturbation theory has an important advantage over other approaches: The calculation yields a good effective classical potential for [*all*]{} temperatures and coupling strengths. The quantum statistical partition function is obtained from a simple integral over a Boltzmann-factor involving the effective classical potential. The ground state energy is then obtained from its zero-temperature limit. The asymptotic behaviour in the strong-coupling limit is emerging automatically and does not have to be derived from other sources.

Effective Classical Representations for the Quantum Statistical Partition Function {#effrep}
==================================================================================

A point particle in $D$ dimensions with a potential $V({\bf x})$ and a vector potential ${\bf A}({\bf x})$ is described by a Hamiltonian $$\label{ham00}
H({\bf p},{\bf x})=\frac{1}{2M}\left[{\bf p}-\frac{e}{c}{\bf A}({\bf x})\right]^2+V({\bf x}).$$ The quantum statistical partition function is given by the euclidean phase space path integral $$\label{ham01}
Z=\oint {\cal D'}^Dx {\cal D}^Dp\,e^{-{\cal A}[{\bf p},{\bf x}]/\hbar}$$ with an action $$\label{ham02}
{\cal A}[{\bf p},{\bf x}]=\int_0^{\hbar\beta}d\tau\left[-i{\bf p}(\tau)\cdot\dot{\bf x}(\tau) 
+H({\bf p}(\tau),{\bf x}(\tau))\right],$$ and the path measure $$\label{ham03}
\oint{\cal D'}^Dx {\cal D}^Dp=\lim_{N\to\infty}\prod\limits_{n=1}^{N+1}
\left[\int\frac{d^Dx_nd^Dp_n}{(2\pi\hbar)^D} \right].$$ The parameter $\beta=1/k_BT$ denotes the usual inverse thermal energy at temperature $T$, where $k_B$ is the Boltzmann constant. From $Z$ we obtain the free energy of the system: $$\label{ham04}
F=-\frac{1}{\beta}{\rm ln}\,Z.$$ In perturbation theory, one treats the external potential $V({\bf x})$ as a small quantity, and expands the partition function into powers of $V({\bf x})$. Such a naive expansion is applicable only for extremely weak couplings, and has a vanishing radius of convergence. Convergence is achieved by variational perturbation theory [@PI], which yields good approximations for all potential strengths, as we shall see in the sequel.

Effective Classical Potential
-----------------------------

All quantum-mechanical systems studied so far in variational perturbation theory were governed by a Hamiltonian of the standard form $$\label{ham05}
H({\bf p},{\bf x})=\frac{{\bf p}^2}{2M}+V({\bf x}).$$ The simple quadratic dependence on the momenta makes the momentum integrals in the path integral (\[ham01\]) trivial. The remaining configuration space representation of the partition function is used to define an effective classical potential $V_{\rm eff}({\bf x}_0)$, from which quantum mechanical partition function is found by a classically looking integral $$\label{ham06}
Z=\int \frac{d^Dx_0}{\lambda_{\rm th}^{D}}\,\exp\left[-\beta V_{\rm eff}({\bf x}_0)\right],$$ where $\lambda_{\rm th}=\sqrt{2\pi\hbar^2\beta/M}$ is the thermal wavelength. The Boltzmann factor plays the role of a [*local partition function*]{} $Z^{{\bf x}_0}$, which is calculated from the restricted path integral $$\label{ham07}
e^{-\beta V_{\rm eff}({\bf x}_0)}\equiv
Z^{{\bf x}_0}=
\lambda_{\rm th}^{D}\oint {\cal D}^Dx\,
\delta({\bf x}_0-\overline{{\bf x}(\tau)})\,e^{-{\cal A}[{\bf x}]/\hbar},$$ with the action $$\label{ham08}
{\cal A}[{\bf x}]=\int_0^{\hbar\beta}d\tau\,\left[\frac{M}{2}\dot{\bf x}^2(\tau)+
V({\bf x}(\tau)) 
\right],$$ and the path measure $$\label{ham09}
\oint{\cal D}^Dx=\lim_{N\to\infty}\prod\limits_{n=1}^{N+1}
\left\{\int\frac{d^Dx_n}{[2\pi\hbar^2\beta/M(N+1)]^{D/2}} \right\}.$$ The special treatment of the temporal average of the Fourier path $$\label{ham10}
{\bf x}_0=\overline{{\bf x}(\tau)}=\frac{1}{\hbar\beta}\int_0^{\hbar\beta}d\tau\,{\bf x}(\tau)$$ is essential for the quality of the results. It subtracts from the harmonic fluctuation width $\langle{\bf x}^2\rangle^{\rm cl}$ the classical divergence proportional to $T=1/k_B\beta$ of the Dulong-Petit law [@PI; @density]. Such diverging fluctuations cannot be treated perturbatively, and require the final integration in expression (\[ham06\]) to be done numerically.

For the Coulomb potential $V({\bf x})=-e^2/4\pi\varepsilon_0\,|{\bf x}|$ in three dimensions, the effective classical potential in Eq. (\[ham07\]) can be approximated well by variational perturbation theory [@PI; @density; @hatom].

Effective Classical Hamiltonian
-------------------------------

In order to deal with Hamiltonians like (\[ham00\]) which contain a ${\bf p}\cdot{\bf A}({\bf x})$-term, we must generalize the variational procedure. Extending (\[ham07\]), we define an [*effective classical Hamiltonian*]{} by the phase space path integral $$\label{ham12}
e^{-\beta H_{\rm eff}({\bf p}_0,{\bf x}_0)}\equiv
Z^{{\bf p}_0,{\bf x}_0}=
(2\pi\hbar)^D\oint {\cal D'}^Dx{\cal D}^Dp\,\delta({\bf x}_0-\overline{{\bf x}(\tau)})
\delta({\bf p}_0-\overline{{\bf p}(\tau)})\,e^{-{\cal A}[{\bf p},{\bf x}]/\hbar},$$ with the action (\[ham02\]) and the measure (\[ham03\]). This allows us to express the partition function as the classically looking phase space integral $$\label{ham11}
Z=\int\frac{d^Dx_0d^Dp_0}{(2\pi\hbar)^D}\,\exp\left[-\beta H_{\rm eff}({\bf p}_0,{\bf x}_0)
\right],$$ where ${\bf p}_0$ is the temporal average of the momentum: $$\label{ham13}
{\bf p}_0=\overline{{\bf p}(\tau)}=\frac{1}{\hbar\beta}\int_0^{\hbar\beta}d\tau\,{\bf p}(\tau).$$ The fixing of ${\bf p}_0$ is done for the same reason as that for ${\bf x}_0$, since the classical expectation value $\langle {\bf p}^2\rangle^{\rm cl}$ is diverging linearly with $T$, just as $\langle{\bf x}^2 \rangle^{\rm cl}$.

In the special case of a standard Hamiltonian (\[ham05\]), the effective Hamiltonian in Eq. (\[ham11\]) reduces to the effective classical potential, since the momentum integral in Eq. (\[ham12\]) can then be easily performed, and the resulting restricted partition function becomes $$\label{ham14}
Z^{{\bf p}_0,{\bf x}_0}=\exp\left(-\beta\frac{{\bf p}_0^2}{2M} \right)\,
Z^{{\bf x}_0}$$ with the local partition function $Z^{{\bf x}_0}=\exp[-\beta V_{\rm eff}({\bf x}_0)]$ of Eq. (\[ham07\]). Thus the complete quantum statistical partition function is given by (\[ham11\]), with an effective classical Hamilton function $$\label{ham15}
H_{\rm eff}({\bf p}_0,{\bf x}_0)=\frac{{\bf p}_0^2}{2M}+V_{\rm eff}({\bf x}_0).$$ As a consequence of the purely quadratic momentum dependence of $H({\bf p},{\bf x})$ in (\[ham05\]), the ${\bf p}_0$-integral in (\[ham11\]) can be done, thus expressing the quantum statistical partition function as a pure configuration space integral over the Boltzmann factor involving the effective classical potential $V_{\rm eff}({\bf x}_0)$, as in Eq. (\[ham06\]).

Exact Effective Classical Hamiltonian for an Electron in a Constant Magnetic Field {#exham}
----------------------------------------------------------------------------------

The effective classical Hamiltonian for the electron moving in a constant magnetic field can be calculated exactly. We consider a magnetic field ${\bf B}=B {\bf e}_z$ pointing along the positive $z$-axis. The only nontrivial motion of the electron is in the $x\!\!-\!\!y$-plane. In symmetric gauge the vector potential is given by $$\label{ham16}
{\bf A}({\bf x})=\frac{B}{2}(-y,x,0).$$ The choice of the gauge does not affect the partition function since the periodic path integral (\[ham01\]) is gauge invariant. Ignoring the trivial free particle motion along the $z$-direction, we may restrict our attention to the two-dimensional Hamiltonian $$\label{ham17}
H({\bf p},{\bf x})=\frac{{\bf p}^2}{2M}-\frac{1}{2}\omega_c l_z({\bf p},{\bf x})+
\frac{1}{8}M\omega_c^2{\bf x}^2$$ with ${\bf x}=(x,y)$ and ${\bf p}=(p_x,p_y)$. Here, $\omega_c=eB/Mc$ is the Landau frequency, and $$\label{ham18}
l_z({\bf p},{\bf x})=({\bf x}\times{\bf p})_z=xp_y-yp_x$$ the third component of the orbital angular momentum. The partition function of the problem is given by Eq. (\[ham11\]), with $D=2$. Being interested in an effective classical formulation, we have to calculate the path integral (\[ham12\]). First we express the $\delta$-function for the averaged momentum as a Fourier integral $$\begin{aligned}
\label{ham19}
\delta({\bf p}_0-\overline{{\bf p}(\tau)})&=&
\int \frac{d^2\xi}{(2\pi\hbar)^2}\,\exp\left(-\frac{i}{\hbar}{\mbox{\boldmath$\xi$}}\cdot{\bf p}_0
\right)\exp\left[-\frac{1}{\hbar}\int_0^{\hbar\beta}d\tau\,
{\bf v}_0({\mbox{\boldmath$\xi$}})\cdot {\bf p}(\tau)\right]\end{aligned}$$ involving an auxiliary source $$\label{ham20}
{\bf v}_0({\mbox{\boldmath$\xi$}})=-\frac{i}{\hbar\beta}\,{\mbox{\boldmath$\xi$}}$$ which is constant in time. Substituting the $\delta$-function in Eq. (\[ham12\]) by this source representation, the partition function reads $$\begin{aligned}
\label{ham21}
Z^{{\bf p}_0,{\bf x}_0}&=&\int d^2\xi\,\exp\left(-\frac{i}{\hbar}{\mbox{\boldmath$\xi$}}\cdot{\bf p}_0\right)
\oint {\cal D'}^2x{\cal D}^2p\,\delta({\bf x}_0-\overline{{\bf x}(\tau)})\nonumber\\
&&\times\exp\left\{-\frac{1}{\hbar}\int_0^{\hbar\beta}d\tau\,
\left[-i{\bf p}(\tau)\cdot\dot{\bf x}(\tau)+H({\bf p}(\tau),{\bf x}(\tau))
+{\bf v}_0({\mbox{\boldmath$\xi$}})\cdot{\bf p}(\tau)\right]\right\}.\end{aligned}$$ Evaluating the momentum integrals and utilizing the periodicity property ${\bf x}(0)={\bf x}(\hbar\beta)$, we obtain the configuration space path integral $$\begin{aligned}
\label{ham22}
Z^{{\bf p}_0,{\bf x}_0}=\lim_{\Omega\to 0}
&&\int d^2\xi\,\exp\left(-\frac{i}{\hbar}{\mbox{\boldmath$\xi$}}\cdot{\bf p}_0-\frac{M}{2\hbar^2\beta}
{\mbox{\boldmath$\xi$}}^2\right)
\oint {\cal D}^2x\,\delta({\bf x}_0-\overline{{\bf x}(\tau)})\nonumber\\
&&\times\exp\left\{-\frac{1}{\hbar}\int_0^{\hbar\beta}d\tau\,\left[\frac{M}{2}\dot{\bf x}^2(\tau) 
+\frac{1}{2}M\Omega^2{\bf x}^2(\tau)
-\frac{i}{2}M\omega_c({\bf x}(\tau)\times \dot{\bf x}(\tau))_z+
{\bf x}(\tau)\cdot {\bf j}_1
({\mbox{\boldmath$\xi$}})\right]\right\},\end{aligned}$$ where the source ${\bf v}_0$ coupled to the momentum in (\[ham21\]) has turned to a source ${\bf j}_1$ coupled to the path in configuration space [@correlation], with components $$\label{ham23}
{\bf j}_1({\mbox{\boldmath$\xi$}})=\frac{M}{2}\omega_c\left(\,v_{0y}({\mbox{\boldmath$\xi$}}),-v_{0x}({\mbox{\boldmath$\xi$}}) 
\,\right)=\frac{i\omega_cM}{2\hbar\beta}\left(\,-\xi_y,\xi_x \,\right).$$ We have introduced an additional harmonic oscillator in Eq. (\[ham22\]) which will turn out to be useful at intermediate stages of the development. At the end of the calculation, only the limit $\Omega\to 0$ will be relevant.

Expressing the $\delta$-function in the path integral of Eq. (\[ham22\]) by the Fourier integral $$\label{ham23b}
\delta({\bf x}_0-\overline{{\bf x}(\tau)})=\int \frac{d^2\kappa}{(2\pi)^2}\,
\exp\left(\frac{i}{\hbar}{\mbox{\boldmath$\kappa$}}\cdot {\bf x}_0 \right)
\exp\left[-\frac{1}{\hbar}\int_0^{\hbar\beta}d\tau\,{\bf j}_2({\mbox{\boldmath$\kappa$}})\cdot
{\bf x}(\tau) \right]$$ with the new source $$\label{ham23c}
{\bf j}_2({\mbox{\boldmath$\kappa$}})=\frac{i{\mbox{\boldmath$\kappa$}}}{\beta},$$ the partition function (\[ham22\]) can be written as $$\label{ham23d}
Z^{{\bf p}_0,{\bf x}_0}=\lim_{\Omega\to 0}
\int d^2\xi\,\exp\left(-\frac{i}{\hbar}{\mbox{\boldmath$\xi$}}\cdot{\bf p}_0-\frac{M}{2\hbar^2\beta}
{\mbox{\boldmath$\xi$}}^2\right)\int\frac{d^2\kappa}{(2\pi)^2}\,\exp\left(\frac{i}{\hbar}
{\mbox{\boldmath$\kappa$}}\cdot {\bf x}_0 \right)\,Z_\Omega[{\bf J}({\mbox{\boldmath$\xi$}},{\mbox{\boldmath$\kappa$}})].$$ The functional $Z_\Omega[{\bf J}({\mbox{\boldmath$\xi$}},{\mbox{\boldmath$\kappa$}})]$ is defined as the configuration space path integral $$\label{ham23e}
Z_\Omega[{\bf J}({\mbox{\boldmath$\xi$}},{\mbox{\boldmath$\kappa$}})]=\oint{\cal D}^2x\,
\exp\left[-\frac{1}{2}\int_0^{\hbar\beta}d\tau\int_0^{\hbar\beta}d\tau'\,
{\bf x}(\tau)\,{\bf G}^{-1}(\tau,\tau')\,{\bf x}(\tau')-\frac{1}{\hbar}\int_0^{\hbar\beta}
d\tau\,{\bf J}({\mbox{\boldmath$\xi$}},{\mbox{\boldmath$\kappa$}})\cdot{\bf x}(\tau) \right],$$ where we have introduced the combined source ${\bf J}({\mbox{\boldmath$\xi$}},{\mbox{\boldmath$\kappa$}})={\bf j}_1({\mbox{\boldmath$\xi$}})+{\bf j}_2({\mbox{\boldmath$\kappa$}})$. Formally, the solution reads $$\label{ham23ee}
Z_\Omega[{\bf J}({\mbox{\boldmath$\xi$}},{\mbox{\boldmath$\kappa$}})]=Z_\Omega[0]
\exp\left[-\frac{1}{2\hbar^2}\int_0^{\hbar\beta}d\tau\int_0^{\hbar\beta}d\tau'\,
{\bf J}({\mbox{\boldmath$\xi$}},{\mbox{\boldmath$\kappa$}})\,{\bf G}(\tau,\tau')\,{\bf J}({\mbox{\boldmath$\xi$}},{\mbox{\boldmath$\kappa$}}) \right],$$ where ${\bf G}(\tau,\tau')$ is the matrix of Green functions obtained by inverting $$\label{ham23f}
{\bf G}^{-1}(\tau,\tau')=\frac{M}{\hbar}\,\left(\begin{array}{cc}
-\frac{d^2}{d \tau^2}+\Omega^2 & -i\omega_c\frac{d}{d\tau}\\
i\omega_c\frac{d}{d\tau} & -\frac{d^2}{d \tau^2}+\Omega^2
\end{array}\right)\,\delta(\tau-\tau').$$ The inversion is easily done in frequency space after spectrally decomposing the $\delta$-function into the Matsubara frequencies $\omega_m=2\pi m/\hbar\beta$, $$\label{ham23g}
\delta(\tau-\tau')=\frac{1}{\hbar\beta}\sum\limits_{m=-\infty}^\infty\,e^{i\omega_m(\tau-\tau')}.$$ The result is $$\label{ham23h}
\tilde{\bf G}(\omega_m)=\frac{\hbar}{M}\frac{1}{{\rm det}\,\tilde{\bf G}}
\left(\begin{array}{cc}\omega_m^2+\Omega^2 & -\omega_c\omega_m\\
 \omega_c\omega_m & \omega_m^2+\Omega^2\end{array}\right).$$ At this point, the additional oscillator in Eq. (\[ham22\]) proves useful: It ensures that the determinant $$\label{ham23i}
{\rm det}\,\tilde{\bf G}(\omega_m)=(\omega_m^2+\Omega^2)^2+\omega_c^2\omega_m^2$$ is nonzero for $m=0$, thus playing the role of an infrared regulator. The Fourier expansion $$\label{ham23j}
{\bf G}(\tau,\tau')=\frac{1}{\hbar\beta}\sum\limits_{m=-\infty}^\infty\,\tilde{\bf G}(\omega_m)
e^{-i\omega_m(\tau-\tau')}$$ yields the matrix of Green functions $$\label{ham23k}
{\bf G}(\tau,\tau')=\left(\begin{array}{cc}G_{xx}(\tau,\tau') & G_{xy}(\tau,\tau')\\
G_{yx}(\tau,\tau') & G_{yy}(\tau,\tau')  \end{array}\right)$$ which inherits the symmetry properties from the kernel (\[ham23f\]): $$\label{ham23l}
G_{xx}(\tau,\tau')=G_{yy}(\tau,\tau'),\qquad G_{xy}(\tau,\tau')=
-G_{yx}(\tau,\tau').$$ A more detailed description of these Green functions is given in Apps. \[appgen\] and \[appprop\].

Since the current ${\bf J}$ does not depend on the euclidean time, the expression (\[ham23ee\]) simplifies therefore to $$\label{ham23m}
Z_\Omega[{\bf J}({\mbox{\boldmath$\xi$}},{\mbox{\boldmath$\kappa$}})]=Z_\Omega[0]
\exp\left[-\frac{1}{\hbar^2}{\bf J}^2({\mbox{\boldmath$\xi$}},{\mbox{\boldmath$\kappa$}})
\int_0^{\hbar\beta}d\tau\int_0^{\hbar\beta}d\tau'\,G_{xx}(\tau,\tau')\right].$$ The Green function has the Fourier decomposition $$\label{ham23n}
G_{xx}(\tau,\tau')=\frac{1}{M\beta}\sum\limits_{m=-\infty}^\infty\,
\frac{\omega_m^2+\Omega^2}{(\omega_m^2+\Omega_+^2)(\omega_m^2+\Omega_-^2)}
\,e^{-i\omega_m(\tau-\tau')},$$ where $\Omega_\pm$ are the frequencies $$\label{ham23o}
\Omega_\pm=\sqrt{\Omega^2+\frac{1}{2}\omega_c^2\pm\omega_c\sqrt{\Omega^2+\frac{1}{4}\omega_c^2}}.$$ The ratios in the sum of (\[ham23n\]) can be decomposed into two partial fractions, each of them representing a single harmonic oscillator with frequency $\Omega_+$ and $\Omega_-$, respectively. The analytic form of the periodic Green function of a single harmonic oscillator is well known (see Chap. 3 in [@PI]), and we obtain for the present Green function (\[ham23o\]): $$\label{ham23p}
G_{xx}(\tau,\tau')=\frac{1}{M\beta}\left(\frac{\hbar\beta}{2\Omega_+}\frac{\Omega_+^2-\Omega^2}
{\Omega_+^2-\Omega_-^2}\frac{\cosh{\Omega_+(|\tau-\tau'|-\hbar\beta/2)}}
{\sinh{\hbar\beta\Omega_+/2}}- \frac{\hbar\beta}{2\Omega_-}\frac{\Omega_-^2-\Omega^2}
{\Omega_+^2-\Omega_-^2}\frac{\cosh{\Omega_-(|\tau-\tau'|-\hbar\beta/2)}}
{\sinh{\hbar\beta\Omega_-/2}}\right).$$ By writing the determinant (\[ham23i\]) as $$\label{ham23ii}
{\rm det}\,\tilde{\bf G}(\omega_m)=(\omega_m^2+\Omega_+^2)(\omega_m^2+\Omega_-^2)$$ and summing over the logarithms of this, we calculate the partition function as a product of two single harmonic oscillators: $$\label{ham23q}
Z_\Omega=Z_\Omega[0]=
\frac{1}{2\sinh{\hbar\beta\Omega_+/2}}\frac{1}{2\sinh{\hbar\beta\Omega_-/2}}.$$ The results (\[ham23p\]) and (\[ham23q\]) determine the generating functional (\[ham23m\]). The euclidean time integrations are then easily done, and subsequently the ${\mbox{\boldmath$\kappa$}}$- and ${\mbox{\boldmath$\xi$}}$-integrations in (\[ham23d\]). As a result, we obtain the restricted partition function $$\label{ham25}
Z^{{\bf p}_0,{\bf x}_0}=\lim_{\Omega\to 0}\exp\left\{-\beta\left(-\frac{1}{\beta}{\rm ln}\,
\frac{\sinh \hbar\beta\Omega_+/2}{\hbar\beta\Omega_+/2}
\frac{\sinh \hbar\beta\Omega_-/2}{\hbar\beta\Omega_-/2}
 +\frac{{\bf p}_0^2}{2M}-\frac{1}{2}\omega_c l_z({\bf p}_0,{\bf x}_0)
+\frac{1}{8}M\omega_c^2{\bf x}_0^2+\frac{1}{2}M\Omega^2{\bf x}_0^2
\right)\right\}.$$ If we now remove the additional oscillator by taking the limit $\Omega\to 0$, we find from (\[ham23o\]): $\Omega_+\rightarrow \omega_c$, $\Omega_-\rightarrow 0$, and therefore $$\label{ham26}
\lim_{\Omega\to 0}\frac{\sinh{\hbar\beta\Omega_+/2}}{\hbar\beta\Omega_+/2}=
\frac{\sinh{\hbar\beta\omega_c/2}}{\hbar\beta\omega_c/2},\quad 
\lim_{\Omega\to 0}\frac{\sinh{\hbar\beta\Omega_-/2}}{\hbar\beta\Omega_-/2}=1.$$ Recalling the definition (\[ham12\]), we identify the exact effective classical Hamiltonian for an electron in a magnetic field as $$\label{ham27}
H_{\rm eff}({\bf p}_0,{\bf x}_0)=\frac{1}{\beta}\,{\rm ln}\,\frac{\sinh{\hbar\beta\omega_c/2}}
{\hbar\beta\omega_c/2}+\frac{{\bf p}_0^2}{2M}-\frac{1}{2}\omega_c\,l_z({\bf p}_0,{\bf x}_0)
+\frac{1}{8}M\omega_c^2{\bf x}_0^2.$$ Integrating out the momenta in Eq. (\[ham11\]), the configuration space representation (\[ham06\]) for the partition function contains the effective classical potential for a charged particle in the plane perpendicular to the direction of a uniform magnetic field $$\label{ham28}
V_{\rm eff}({\bf x}_0)=\frac{1}{\beta}\,{\rm ln}\,\frac{\sinh{\hbar\beta\omega_c/2}}
{\hbar\beta\omega_c/2}.$$ Note that this is a constant potential.

Denoting the area $\int d^2x_0$ by $A$, we find the exact quantum statistical partition function $$\label{ham29}
Z=\frac{A}{\lambda_{\rm th}^2}\frac{\hbar\beta\omega_c/2}{\sinh{\hbar\beta\omega_c/2}}.$$ After these preparations, we can turn our attention to the system we want to study in this paper: the hydrogen atom in a uniform magnetic field, where the additional Coulomb interaction prevents us from finding an exact solution for the effective classical Hamilton function.

Hydrogen Atom in Constant Magnetic Field
========================================

The zero-temperature properties of the hydrogen atom without external fiels are exactly known. For the quantum statistics at finite temperatures, an analytic expression exists, but it is hard to evaluate. It is easier to find an accurate approximative result with the help of variational perturbation theory [@hatom]. Similar calculations have been performed for the electron-proton pair distribution function which can be interpreted as the unnormalized density matrix [@density].

Here we extend this method of calculation to the hydrogen atom in a constant magnetic field. This extension is quite nontrivial since the weak- and strong-field limits will turn out to exhibit completely different asymptotic behaviours. Let us first generalize variational perturbation theory to an electron in a constant magnetic field and arbitrary potential.

Generalized Variational Perturbation Theory
-------------------------------------------

We consider once more the effective classical form (\[ham11\]) of the quantum statistical partition function, which requires the path integration (\[ham12\]) in phase space. Fluctuations parallel and vertical to the magnetic field lines are now both nontrivial, and we must deal with the full three-dimensional system and the components of the electron position and momentum are now denoted by ${\bf x}=(x,y,z)$ and ${\bf p}=(p_x,p_y,p_z)$. For the uniform magnetic field pointing along the $z$-axis, the vector potential ${\bf A}({\bf x})$ is used in the gauge (\[ham16\]). Thus the Hamilton function of an electron in a magnetic field and an arbitrary potential $V({\bf x})$ is $$\label{vpt00}
H({\bf p},{\bf x})=\frac{{\bf p}^2}{2M}-\frac{1}{2}\omega_cl_z({\bf p},{\bf x})+\frac{1}{8}
M\omega_c^2{\bf x}^2+V({\bf x}).$$ The orbital angular momentum $l_z({\bf p},{\bf x})$ was introduced in Eq. (\[ham18\]), and the Landau frequency $\omega_c$ below Eq. (\[ham17\]). The importance of the separation of the zero frequency components ${\bf x}_0$ and ${\bf p}_0$ was discussed in Sect. \[effrep\]. Their divergence with the temperature $T$ prevents a perturbative treatment. Thus it is essential to set up the perturbation theory only for the fluctuations around ${\bf x}_0$ and ${\bf p}_0$. For this we rewrite the action functional (\[ham02\]) associated with the Hamiltonian (\[vpt00\]) as $$\label{vpt01}
{\cal A}[{\bf p},{\bf x}]=
{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf p},{\bf x}]
+{\cal A}_{\rm int}[{\bf p},{\bf x}],$$ where we have introduced the fluctuation action $$\begin{aligned}
\label{vpt02}
{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf p},{\bf x}]=
\int_0^{\hbar\beta}d\tau\, \Big\{&&-i[{\bf p}(\tau)-{\bf p}_0]\cdot\dot{\bf x}(\tau)+
\frac{1}{2M}[{\bf p}(\tau)-{\bf p}_0]^2+\frac{1}{2}{\Omega_{\perp 1}}l_z({\bf p}(\tau)-{\bf p}_0,{\bf x}(\tau)-{\bf x}_0)\nonumber\\
&&+\frac{1}{8}M{\Omega_{\perp 2}}^2
\left[{\bf x}^\perp(\tau)-{\bf x}_0^\perp\right]^2+\frac{1}{2}M{\Omega_{\parallel}}^2[z(\tau)-z_0]^2 \Big\},\end{aligned}$$ in which ${\bf x}^\perp=(x,y)$ denotes the transverse part of ${\bf x}$. The interaction is now $$\label{vpt03}
{\cal A}_{\rm int}[{\bf p},{\bf x}]=\int_0^{\hbar\beta}d\tau\,
V_{\rm int}({\bf p}(\tau),{\bf x}(\tau))
={\cal A}[{\bf p},{\bf x}]-
{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf p},{\bf x}]$$ with the interaction potential $$\begin{aligned}
\label{vpt03b}
V_{\rm int}({\bf p}(\tau),{\bf x}(\tau))&=&\frac{1}{2M}\left\{{\bf p}^2(\tau)-
\left[{\bf p}(\tau)-{\bf p}_0 \right]^2 \right\}+\frac{1}{2}\omega_c\,
{\bf p}^\perp(\tau)\times
{\bf x}^\perp(\tau)
\nonumber\\
&&-\frac{1}{2}{\Omega_{\perp 1}}({\bf p}^\perp(\tau)-{\bf p}^\perp_0)\times
({\bf x}^\perp(\tau)-{\bf x}^\perp_0)
+\frac{1}{8}M\omega_c^2{{\bf x}^\perp}^2(\tau)\nonumber\\
&&-\frac{1}{8}M{\Omega_{\perp 2}}^2
\left[{\bf x}^\perp(\tau)-{\bf x}_0^\perp\right]^2
-\frac{1}{2}M{\Omega_{\parallel}}^2[z(\tau)-z_0]^2+V({\bf x}(\tau)),\end{aligned}$$ where ${\bf p}^\perp=(p_x,p_y)$. The frequencies ${\bf \Omega}=({\Omega_{\perp 1}},{\Omega_{\perp 2}},{\Omega_{\parallel}})$ are for the moment arbitrary. The decomposition (\[vpt01\]) forms the basis for the variational approach, where the first term in the action (\[vpt01\]) allows an exact treatment. The transverse part was given in Sec. \[exham\] and the longitudinal part is trivial, since it is harmonic with frequency ${\Omega_{\parallel}}$. The associated partition function is given by the path integral $$\label{vpt04}
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}=
\oint{\cal D'}^3x{\cal D}^3p\,\delta({\bf x}_0-\overline{{\bf x}(\tau)})
\delta({\bf p}_0-\overline{{\bf p}(\tau)})
e^{-{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf p},{\bf x}]/\hbar},$$ which can be performed. Details are given in Appendix \[appgenrp\]. The result is $$\label{vpt04b}
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}=
\frac{\hbar\beta\Omega_+/2}{\sinh{\hbar\beta\Omega_+/2}}\,
\frac{\hbar\beta\Omega_-/2}{\sinh{\hbar\beta\Omega_-/2}}\,
\frac{\hbar\beta{\Omega_{\parallel}}/2}{\sinh{\hbar\beta{\Omega_{\parallel}}/2}},$$ where auxiliary frequencies are composed of the frequencies ${\Omega_{\perp 1}},{\Omega_{\perp 2}}$ in the action (\[vpt02\]) as $$\label{vpt04c}
\Omega_{\pm}({\Omega_{\perp 1}},{\Omega_{\perp 2}})=\frac{1}{2}\,|{\Omega_{\perp 1}}\pm{\Omega_{\perp 2}}|.$$ This partition function serves in the subsequent pertubation expansion as trial system which depends explicitly on the frequencies ${\bf \Omega}$. The correlation functions are a straightforward generalization of (\[ham23k\]) to three dimensions: $$\label{hyd05c}
{\bf G}^{{\bf x}_0}(\tau,\tau')=\left(\begin{array}{ccc}
G_{xx}^{{\bf x}_0}(\tau,\tau') & G_{xy}^{{\bf x}_0}(\tau,\tau') & 0\\
G_{yx}^{{\bf x}_0}(\tau,\tau') & G_{yy}^{{\bf x}_0}(\tau,\tau') & 0\\
0 & 0 & G_{zz}^{{\bf x}_0}(\tau,\tau')\end{array}\right),$$ whose explicit form is derived in App. \[appgenrp\].

The ${\bf \Omega}$-dependent action in Eq. (\[vpt01\]) is treated perturbatively. Writing the partition function (\[ham12\]) as $$\label{vpt05}
Z^{{\bf p}_0,{\bf x}_0}=
(2\pi\hbar)^3\oint {\cal D'}^3x{\cal D}^3p\,\delta({\bf x}_0-\overline{{\bf x}(\tau)})
\delta({\bf p}_0-\overline{{\bf p}(\tau)})\,
\exp\left\{-\frac{1}{\hbar}{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
[{\bf p},{\bf x}]\right\}
\exp\left\{-\frac{1}{\hbar}\int_0^{\hbar\beta}d\tau\,
V_{\rm int}({\bf p}(\tau),{\bf x}(\tau))\right\},$$ the second exponential is expanded into a Taylor series, yielding $$\begin{aligned}
\label{vpt06}
Z^{{\bf p}_0,{\bf x}_0}&=&
(2\pi\hbar)^3\oint {\cal D'}^3x{\cal D}^3p\,\delta({\bf x}_0-\overline{{\bf x}(\tau)})
\delta({\bf p}_0-\overline{{\bf p}(\tau)})\,
\exp\left\{-\frac{1}{\hbar}{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
[{\bf p},{\bf x}]\right\}\nonumber\\
&&\times\left[1-\frac{1}{\hbar}\int_0^{\hbar\beta}d\tau\,V_{\rm int}({\bf p}(\tau),{\bf x}(\tau))
+\frac{1}{2!\hbar^2}\int_0^{\hbar\beta}d\tau_1\int_0^{\hbar\beta}d\tau_2\,
V_{\rm int}({\bf p}(\tau_1),{\bf x}(\tau_1))V_{\rm int}({\bf p}(\tau_2),{\bf x}(\tau_2))-
\ldots\right].\end{aligned}$$ Defining harmonic expectation values with respect to the restricted path integral as $$\label{vpt07}
{\left\langle\,\ldots\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=
\frac{(2\pi\hbar)^3}{Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}\oint{\cal D'}^3x {\cal D}^3p
\,\ldots\,\delta({\bf x}_0-\overline{{\bf x}(\tau)})\delta({\bf p}_0-\overline{{\bf p}(\tau)})
\exp\left\{-\frac{1}{\hbar}{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf p},{\bf x}]
\right\},$$ the perturbation expansion for the partition function (\[vpt06\]) reads $$\label{vpt08}
Z^{{\bf p}_0,{\bf x}_0}=Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}\,\sum\limits_{n=0}^\infty
\frac{(-1)^n}{\hbar^nn!}\,
{\left\langle\,\left(\int_0^{\hbar\beta}d\tau\,V_{\rm int}({\bf p}(\tau),{\bf x}(\tau))\right)^n\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}.$$ This power series expansion can be rewritten in the exponential form $$\label{vpt09}
Z^{{\bf p}_0,{\bf x}_0}=Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}\,\exp\left\{
\sum\limits_{n=0}^\infty
\frac{(-1)^n}{\hbar^nn!}\,
{\left\langle\,\left(\int_0^{\hbar\beta}d\tau\,V_{\rm int}({\bf p}(\tau),{\bf x}(\tau))\right)^n\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{{\bf \Omega},c}}\right\},$$ where the subscript $c$ on the expectation values indicates cumulants. The lowest cumulants are related to the full expectation values as follows: $$\begin{aligned}
\label{vpt10}
{\left\langle\,O_1({\bf p}(\tau_1),{\bf x}(\tau_1))\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{{\bf \Omega},c}}&=&{\left\langle\,O_1({\bf p}(\tau_1),{\bf x}(\tau_1))\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}},
\nonumber\\
{\left\langle\,O_1({\bf p}(\tau_1),{\bf x}(\tau_1))O_2({\bf p}(\tau_2),{\bf x}(\tau_2))\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{{\bf \Omega},c}}
&=&{\left\langle\,O_1({\bf p}(\tau_1),{\bf x}(\tau_1))O_2({\bf p}(\tau_2),{\bf x}(\tau_2))\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}\nonumber\\
&&-{\left\langle\,O_1({\bf p}(\tau_1),{\bf x}(\tau_1))\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}
{\left\langle\,O_2({\bf p}(\tau_2),{\bf x}(\tau_2))\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}},\nonumber\\
&\vdots&\quad,\end{aligned}$$ where $O_i({\bf p}(\tau_j),{\bf x}(\tau_j))$ denotes any observable depending on position and momentum. Recalling the relation (\[ham12\]) between partition function (\[vpt09\]) and effective classical Hamiltonian $H_{\rm eff}({\bf p}_0,{\bf x}_0)$, we obtain from (\[vpt09\]) the effective classical Hamiltonian as a cumulant expansion: $$\label{vpt11}
H_{\rm eff}({\bf p}_0,{\bf x}_0)=-\frac{1}{\beta}{\rm ln}\,Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
+\frac{1}{\beta}\sum\limits_{n=1}^\infty
\frac{(-1)^n}{\hbar^nn!}\,
{\left\langle\,\left(\int_0^{\hbar\beta}d\tau\,V_{\rm int}({\bf p}(\tau),{\bf x}(\tau))\right)^n\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{{\bf \Omega},c}}.$$ Up to now, we did not make any approximation. The expansion on the right-hand side is an exact expression for the effective classical Hamiltonian for any ${\bf \Omega}$.

For systems with a nontrivial interaction, we are capable of calculating only some initial truncated part of the series (\[vpt11\]), say up to the $N$th order, leading to the approximate effective classical Hamiltonian $$\label{vpt12}
{\cal H}_{\bf \Omega}^{(N)}({\bf p}_0,{\bf x}_0)=-
\frac{1}{\beta}{\rm ln}\,Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
+\frac{1}{\beta}\sum\limits_{n=1}^N
\frac{(-1)^n}{\hbar^nn!}\,
{\left\langle\,\left(\int_0^{\hbar\beta}d\tau\,V_{\rm int}({\bf p}(\tau),{\bf x}(\tau))\right)^n\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{{\bf \Omega},c}}.$$ This depends explicitly on the three parameters ${\bf \Omega}$. Since the exact expression (\[vpt11\]) is independent of ${\bf \Omega}$, the best approximation for ${\cal H}_{\bf \Omega}^{(N)}({\bf p}_0,{\bf x}_0)$ should depend on ${\bf \Omega}$ [*minimally*]{}. Thus the optimal solution will be found by determining the parameters from the conditions $$\label{vpt13}
\nabla_{\bf \Omega} {\cal H}_{\bf \Omega}^{(N)}({\bf p}_0,{\bf x}_0)\stackrel{!}{=}0.$$ Let us denote the optimal variational parameters to $N$th order by $$\label{vpt14}
{\bf \Omega}^{(N)}=\left({\Omega_{\perp 1}}^{(N)}({\bf p}_0,{\bf x}_0), 
{\Omega_{\perp 2}}^{(N)}({\bf p}_0,{\bf x}_0), {\Omega_{\parallel}}^{(N)}({\bf p}_0,{\bf x}_0)\right).$$ Inserting these into Eq. (\[vpt12\]) yields the optimal effective classical Hamiltonian ${\cal H}^{(N)}({\bf p}_0,{\bf x}_0)$.

First-Order Effective Classical Potential
-----------------------------------------

The first-order approximation of the effective classical Hamiltonian (\[vpt12\]) reads $$\label{fir00}
{\cal H}_{\bf \Omega}^{(1)}({\bf p}_0,{\bf x}_0)=-
\frac{1}{\beta}{\rm ln}\,Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
-{\left\langle\,V_{\rm int}({\bf p},{\bf x})\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}.$$ In writing the last term we have used the fact that, as a consequence of the time translation invariance of the system, the first-order expectation value of $V_{\rm int}({\bf x})$ is independent of the euclidean time $\tau$.

In order to calculate ${\cal H}_{\bf \Omega}^{(1)}({\bf p}_0,{\bf x}_0)$, we use the two-point correlation functions derived in App. \[appgenrp\], and the vanishing of the linear expectations, e.g. $$\label{fir01}
{\left\langle\,p_x(\tau)-{p_0}_x\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=0$$ to find $$\label{fir02}
{\cal H}_{\bf \Omega}^{(1)}({\bf p}_0,{\bf x}_0)=
\frac{{\bf p}_0^2}{2M}-\frac{1}{2}\omega_c
l_z({\bf p}_0,{\bf x}_0)+\frac{1}{8}M\omega_c^2(x_0^2+y_0^2)+
W_{\bf \Omega}^{(1)}({\bf x}_0),$$ where we have collected all terms depending on the variational parameters ${\bf \Omega}$ in the potential $$\label{fir03}
W_{\bf \Omega}^{(1)}({\bf x}_0)=-\frac{1}{\beta}{\rm ln}\,
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
+(\omega_c-{\Omega_{\perp 1}})\,b^2_\perp({\bf x}_0)-\frac{1}{4}\left({\Omega_{\perp 2}}^2-\omega_c^2 \right)
\,a^2_\perp({\bf x}_0)-\frac{1}{2}M{\Omega_{\parallel}}^2a^2_\parallel({\bf x}_0)+
{\left\langle\,V({\bf x})\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}.$$ The quantities $a^2_\perp({\bf x}_0)$ and $a^2_\parallel({\bf x}_0)$ are the transverse and longitudinal fluctuation widths $$\label{fir03b}
a^2_\perp({\bf x}_0)=G^{{\bf p}_0,{\bf x}_0}_{xx}(0),
\quad a^2_\parallel({\bf x}_0)=G^{{\bf p}_0,{\bf x}_0}_{zz}(0),
\quad b^2_\perp({\bf x}_0)=G^{{\bf p}_0,{\bf x}_0}_{xp_y}(0).$$ Note that the potential (\[fir03\]) is independent of ${\bf p}_0$. This means that the approximation (\[fir02\]) to the effective classical Hamiltonian contains no coupling of the momentum ${\bf p}_0$ to a variational parameter ${\bf \Omega}$, such that the optimal ${\bf \Omega}^{(1)}$ determined by minimizing ${\cal H}_{\bf \Omega}^{(1)}({\bf p}_0,{\bf x}_0)$ is independent of ${\bf p}_0$. We may therefore integrate out ${\bf p}_0$ in the phase space representation of the first-order approximation for the partition function $$\label{fir04}
Z^{(1)}=\int\frac{d^3x_0d^3p_0}{(2\pi\hbar^3)}\,
e^{-\beta{\cal H}_{\bf \Omega}^{(1)}({\bf p}_0,{\bf x}_0)}$$ to find the pure configuration space integral $$\label{fir05}
Z^{(1)}=\int\frac{d^3x_0}{\lambda_{\rm th}^3}\,
e^{-\beta W_{\bf \Omega}^{(1)}({\bf x}_0)},$$ in which $W_{\bf \Omega}^{(1)}({\bf x}_0)$ is the first-order approximation to the effective classical potential of an electron in a potential $V({\bf x})$ and a uniform magnetic field.

Application to the Hydrogen Atom in a Magnetic Field
----------------------------------------------------

We now apply the formulas of the preceding section to the Hamiltonian (\[vpt00\]) with an attracting Coulomb potential $$\label{hyd00}
V({\bf x})=-\frac{e^2}{4\pi\varepsilon_0\,|{\bf x}|}\,,$$ where $|{\bf x}|$ is the distance between the electron and the proton. The only nontrivial problem is the calculation of the expectation value ${\left\langle\,V({\bf x}(\tau))\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}$ in Eq. (\[fir03\]). This is done using the so-called [*smearing formula*]{}, which is a Gaussian convolution of $V({\bf x})$. This formula was first derived by Feynman and Kleinert [@fk], and exists now also in an extension to arbitrary order [@density; @hatom]. The generalization to position and momentum dependent observables was given in the phase space formulation [@correlation]. We briefly rederive the first-order smearing formula. The expectation value is defined by $$\label{hyd01}
{\left\langle\,V({\bf x}(\tau'))\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=\frac{(2\pi\hbar)^3}{Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}
\oint {\cal D'}^3x{\cal D}^3p\,V({\bf x}(\tau'))\,\delta({\bf x}_0-\overline{{\bf x}(\tau)})
\delta({\bf p}_0-\overline{{\bf p}(\tau)})
e^{-{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf p},{\bf x}]/\hbar},$$ Now we substitute the potential by the expression $$\label{hyd02}
V({\bf x}(\tau'))=\int d^3x\,V({\bf x})\delta({\bf x}-{\bf x}(\tau'))
=\int d^3x\,V({\bf x})\int\,\frac{d^3\kappa}{(2\pi)^3}\,\exp\left[i{\mbox{\boldmath$\kappa$}}\cdot
({\bf x}-{\bf x}_0) \right]
\exp\left\{-\frac{1}{\hbar}
\int_0^{\hbar\beta}d\tau\,{\bf j}(\tau)\cdot[{\bf x}(\tau)-{\bf x}_0]
\right\},$$ where we have introduced the source $$\label{hyd03}
{\bf j}(\tau)=i\hbar{\mbox{\boldmath$\kappa$}}\delta(\tau-\tau').$$ Inserting the expression (\[hyd02\]) into Eq. (\[hyd01\]) we obtain $$\label{hyd04}
{\left\langle\,V({\bf x}(\tau'))\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=\frac{1}{Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}
\int d^3x\,V({\bf x})\int\,\frac{d^3\kappa}{(2\pi)^3}\,\exp\left[i{\mbox{\boldmath$\kappa$}}\cdot
({\bf x}-{\bf x}_0) \right]\,Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf j}],$$ with the harmonic generating functional $$\label{hyd05}
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf j}]=
(2\pi\hbar)^3\oint {\cal D'}^3x{\cal D}^3p\,\delta({\bf x}_0-\overline{{\bf x}(\tau)})
\delta({\bf p}_0-\overline{{\bf p}(\tau)})
\exp\left\{-\frac{1}{\hbar}{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf p},{\bf x}]
-\frac{1}{\hbar}\int_0^{\hbar\beta}d\tau\,{\bf j}(\tau)\cdot[{\bf x}(\tau)-{\bf x}_0]\right\}.$$ The solution is $$\label{hyd05b}
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf j}]=
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
\exp\left[\frac{1}{2\hbar^2}\int_0^{\hbar\beta}d\tau\int_0^{\hbar\beta}d\tau'\,
{\bf j}(\tau)\,{\bf G}^{{\bf x}_0}(\tau,\tau')\,{\bf j}(\tau')\right]$$ with the $3\times 3$-matrix of Green functions of Eq. (\[hyd05c\]). The properties of the Green functions are discussed in the Appendices \[appgen\] and  \[appprop\]. Expressing the source ${\bf j}(\tau)$ in terms of ${\mbox{\boldmath$\kappa$}}$ via Eq. (\[hyd03\]) and performing the $\tau$-integrations, we arrive at $$\label{hyd06}
{\left\langle\,V({\bf r}(\tau'))\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=\int d^3x\,V({\bf x})\int\,\frac{d^3\kappa}{(2\pi)^3}\,
\exp\left\{i{\mbox{\boldmath$\kappa$}}\cdot[{\bf x}-{\bf x}_0] \right\}\,
\exp\left[-\frac{1}{2}{\mbox{\boldmath$\kappa$}}\,{\bf G}^{{\bf x}_0}(0)\,{\mbox{\boldmath$\kappa$}} \right].$$ Recognizing that $G^{{\bf x}_0}_{yx}(0)=G^{{\bf x}_0}_{xy}(0)$ vanish, the ${\mbox{\boldmath$\kappa$}}$-integral is easily calculated and leads to the first-order smearing formula for an arbitrary position dependent potential $$\label{hyd07}
{\left\langle\,V({\bf x}(\tau'))\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=\frac{1}{(2\pi)^{3/2}\,a^2_\perp({\bf x}_0)
\sqrt{a^2_\parallel({\bf x}_0)}}\int d^3x
\,V({\bf x})\,
\exp\left[-\frac{(x-x_0)^2+(y-y_0)^2}{2a^2_\perp({\bf x}_0)}
-\frac{(z-z_0)^2}{2a^2_\parallel({\bf x}_0)} \right],$$ the right-hand side containing the Gaussian fluctuation widths (\[fir03b\]).

For the Coulomb potential (\[hyd00\]) that we are interested in, the integral in the smearing formula (\[hyd07\]) can not be done exactly. An integral representation for a simple numerical treatment is $$\begin{aligned}
\label{hyd08}
{\left\langle\,-\frac{e^2}{4\pi\varepsilon_0\,|{\bf x}|}\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}&=&-\frac{e^2}{4\pi\varepsilon_0}
\sqrt{\frac{2}{\pi}\,a^2_\parallel({\bf x}_0)}\int\limits_0^1\frac{d\xi}
{a^2_\parallel({\bf x}_0)+\xi^2[a^2_\perp({\bf x}_0)-a^2_\parallel({\bf x}_0)]}\nonumber\\
&&\times\exp\left\{-\frac{\xi^2}{2}\left(\frac{x_0^2+y_0^2}{a^2_\parallel({\bf x}_0)+
\xi^2[a^2_\perp({\bf x}_0)-a^2_\parallel({\bf x}_0)]}+\frac{z_0^2}{a^2_\parallel({\bf x}_0)} 
\right) \right\}.\end{aligned}$$ With this expression we know the entire first-order effective classical potential (\[fir03\]) for an electron in a Coulomb potential and a uniform magnetic field which has to be optimized in the variational parameters ${\bf \Omega}$.

Results {#res}
=======

We are now going to optimize the effective classical potential by extremizing it in ${\bf \Omega}$ at different temperatures and magnetic field strengths. In the zero-temperature limit this will produce the ground state energy.

Effective Classical Potential for Different Temperatures and Magnetic Field Strengths
-------------------------------------------------------------------------------------

The optimization of $W_{{\bf \Omega}}^{(1)}({\bf x}_0)$ proceeds by minimization in ${\bf \Omega}$ and must be done for each value of ${\bf x}_0$. Reinserting the optimal parameters ${\bf \Omega}^{(1)}({\bf x}_0)$ into the expressions (\[fir03\]) and (\[hyd08\]), we obtain the optimal first-order effective classical potential $W^{(1)}({\bf x}_0)$. The calculations are done numerically, where we used natural units $\hbar=e^2/4\pi\varepsilon_0=k_B=c=M=1$. This means that energies are measured in units of $\epsilon_0=Me^4/(4\pi\varepsilon_0)^2\hbar^2\equiv 2\,{\rm Ryd}\approx 27.21\,{\rm eV}$, temperatures in $\epsilon_0/k_B\approx 3.16\times 10^5\,{\rm K}$, distances in Bohr radii $a_B=(4\pi\varepsilon_0)^2\hbar^2/Me^2\approx
0.53\times 10^{-10}\,{\rm m}$, and magnetic field strengths in $B_0=e^3M^2/\hbar^3(4\pi\varepsilon_0)^2\approx 2.35\times 10^5\,{\rm T}=2.35\times10^9\,
{\rm G}$. Figure \[wofq1\] shows the resulting curves for various magnetic field strengths $B$ and an inverse tempature $\beta=1/T=1$. Examples of the lower temperature behaviour are shown in Fig. \[wofq100\] for $\beta=100$. To see the expected anisotropy of the curves in the magnetic field direction and in the plane perpendicular to it, we plot simultanously the curves for $W^{(1)}({\bf x}_0)$ transversal to the magnetic field as a function of $\rho_0=\sqrt{x_0^2+y_0^2}$ at $z=0$ (solid curves) and parallel as a function of $z_0$ at $\rho_0=0$ (dashed curves). The curves become strongly anisotropic for low temperatures and increasing field strengths (Fig. \[wofq100\]). At a given field strength $B$, the two curves converge for large distances from the origin, where the proton resides, to the same constant depending on $B$. This is due to the decreasing influence of the Coulomb interaction which shows the classical $1/r$-behaviour in each direction. When approaching the classical high-temperature limit, the effect of anisotropy becomes less important since the violent thermal fluctuations do not have a preferred direction (see Fig. \[wofq1\]). For $\rho_0\to\infty$ or $z_0\to\infty$, the expectation value of the Coulomb potential (\[hyd08\]) tends to zero. The remaining effective classical potential $$\label{res01}
W_{\bf \Omega}^{(1)}({\bf x}_0)\longrightarrow -\frac{1}{\beta}{\rm ln}\,
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
+(\omega_c-{\Omega_{\perp 1}})\,b^2_\perp-\frac{1}{4}\left({\Omega_{\perp 2}}^2-\omega_c^2 \right)
\,a^2_\perp-\frac{1}{2}M{\Omega_{\parallel}}^2a^2_\parallel$$ is a constant with regard to the position ${\bf x}_0$, and the optimization yields ${\Omega_{\perp 1}}^{(1)}={\Omega_{\perp 2}}^{(1)}=\omega_c$ and ${\Omega_{\parallel}}^{(1)}=0$, leading to the asymptotic constant value $$\label{res02}
W^{(1)}({\bf x}_0)\longrightarrow -\frac{1}{\beta}{\rm ln}\,
\frac{\hbar\beta\omega_c/2}{\sinh{\hbar\beta\omega_c/2}}.$$ The $B=0$ -curves are of course identical with those obtained from variational perturbation theory for the hydrogen atom [@hatom].

Ground State Energy of the Hydrogen Atom in Uniform Magnetic Field {#vptlow}
------------------------------------------------------------------

In what follows we investigate the zero-temperature behaviour of the theory. Figures \[wofq1\] and \[wofq100\] show that the minimum of each potential curve lies at the origin. This means that the first-order approximation to the ground state energy for a fixed magnitude of the magnetic field $B$ is found by considering the zero-temperature limit of the first-order effective classical potential in the origin $$\label{gr00} 
E^{(1)}=\lim_{\beta\to\infty}W^{(1)}(0).$$ Thus we obtain from Eq. (\[fir03\]) the variational expression for the ground state energy: $$\label{gr01}
E^{(1)}_{\bf \Omega}(B)=\frac{\hbar}{4{\Omega_{\perp 2}}}\left({\Omega_{\perp 2}}^2+\omega_c^2\right)+\frac{\hbar{\Omega_{\parallel}}}{4}-
\frac{e^4}{4\pi\varepsilon_0}{\left\langle\,\frac{1}{|{\bf x}|}\,\right\rangle^{\bf 0}_{\bf \Omega}},$$ where the expectation value for the Coulomb potential (\[hyd08\]) can now be calculated exactly since the exponential in the integral simplifies to unity: $$\label{gr02}
{\left\langle\,\frac{1}{|{\bf x}|}\,\right\rangle^{\bf 0}_{\bf \Omega}}=2\sqrt{\frac{M}{\pi\hbar}}\times\left\{
\begin{array}{cl}  
{\displaystyle}\sqrt{\frac{{\Omega_{\parallel}}{\Omega_{\perp 2}}}{{\Omega_{\parallel}}-{\Omega_{\perp 2}}}}\,{\rm arctan}\,\sqrt{\frac{2{\Omega_{\parallel}}}{{\Omega_{\perp 2}}}-1}, & \quad 2{\Omega_{\parallel}}> {\Omega_{\perp 2}},\\
{\displaystyle}\sqrt{{\Omega_{\parallel}}}, & \quad 2{\Omega_{\parallel}}= {\Omega_{\perp 2}},\\
{\displaystyle}\frac{1}{2i}\sqrt{\frac{{\Omega_{\parallel}}{\Omega_{\perp 2}}}{{\Omega_{\parallel}}-{\Omega_{\perp 2}}}}\,{\rm ln}\,
\frac{1+i\sqrt{2{\Omega_{\parallel}}/{\Omega_{\perp 2}}-1}}{1-i\sqrt{2{\Omega_{\parallel}}/{\Omega_{\perp 2}}-1}}, & \quad 2{\Omega_{\parallel}}< {\Omega_{\perp 2}}.
\end{array}\right.$$ The equations (\[gr01\]) and (\[gr02\]) are independent of the frequency parameter ${\Omega_{\perp 1}}$ such that the optimization of the first-order expression for the ground state energy (\[gr01\]) requires the satisfying of the equations $$\label{gr03}
\frac{\partial E^{(1)}_{\bf \Omega}(B)}{\partial {\Omega_{\perp 2}}}\stackrel{!}{=}0,\qquad
\frac{\partial E^{(1)}_{\bf \Omega}(B)}{\partial {\Omega_{\parallel}}}\stackrel{!}{=}0.$$ Reinserting the resulting values ${\Omega_{\perp 2}}^{(1)}$ and ${\Omega_{\parallel}}^{(1)}$ into Eq. (\[gr01\]) yields the first-order approximation for the ground state energy $E^{(1)}(B)$. In the absence of the Coulomb interaction the optimization with respect to ${\Omega_{\perp 2}}$ yields ${\Omega_{\perp 2}}^{(1)}=\omega_c$, rendering the ground state energy $E^{(1)}(B)=\omega_c/2$, which is the zeroth Landau level. An optimal value for ${\Omega_{\parallel}}$ does not exist since the dependece of the ground state energy of this parameter is linear in Eq. (\[gr01\]) in this special case. To obtain the lowest energy, this parameter can be set to zero (all optimal frequency parameters used in the optimization procedure turn out to be nonnegative). For a vanishing magnetic field, $B=0$, Eq. (\[gr01\]) exactly reproduces the first-order variational result for the ground state energy of the hydrogen atom, $E^{(1)}(B=0)\approx -0.42\, {\rm [2\,Ryd]}$, obtained in Ref. [@hatom].

To investigate the asymptotics in the strong-field limit $B\to\infty$, it is useful to extract the leading term $\omega_c/2$. Thus we define the binding energy $$\label{gr04}
\varepsilon(B)\equiv\frac{\omega_c}{2}-E(B)$$ which possesses an characteristic strong-field behaviour to be discussed in detail subsequently. The result is shown in Fig. \[bofb\] as a function of the magnitude of the magnetic field $B$, where it is compared with the high-accuracy results of Ref. [@zinn-justin]. As a first-order approximation, this result is satisfactory. It is of the same quality like other first-order results, for example those from the operator optimization method in first order of Ref. [@feranshuk]. The advantage of variational perturbation theory is that it yields good results over the complete range of the coupling strength, here the magnetic field. Moreover, as a consequence of the exponential convergence [@PI Chap. 5], higher orders of variational perturbation theory push the approximative result of any quantity very rapidly towards the exact value.

### The Weak-Field Case

We investigate now the weak-field behaviour of our theory starting from the expression (\[gr04\]) and the expectation value of the Coulomb potential (\[gr02\]) in natural units: $$\label{we00}
\varepsilon^{(1)}_{\eta,\Omega}(B)=
\frac{B}{2}-\frac{\Omega}{4}\left(1+\frac{\eta}{2}\right)-\frac{B^2}{4\Omega}
-\sqrt{\frac{\eta\Omega}{2\pi}}h(\eta)$$ with $$\label{we01}
h(\eta)=\frac{1}{\sqrt{1-\eta}}\,{\rm ln}\,\frac{1-\sqrt{1-\eta}}{1+\sqrt{1-\eta}}.$$ In comparison with Eq. (\[gr01\]) we introduced new variational parameters $$\label{we02}
\eta\equiv\frac{2{\Omega_{\parallel}}}{{\Omega_{\perp 2}}},\qquad \Omega\equiv{\Omega_{\perp 2}}$$ and utilized, as the calculations for the binding energy showed, that always $\eta\le 1$. Performing the derivatives with respect to these variational parameters and setting them zero yields conditional equations which can be written after some manipulations as $$\begin{aligned}
\label{we03}
\frac{\Omega}{8}+\sqrt{\frac{\Omega}{2\pi\eta}}\frac{1}{1-\eta}\left(1+\frac{1}{2}
\frac{1}{\sqrt{1-\eta}}{\rm ln}\,\frac{1-\sqrt{1-\eta}}{1+\sqrt{1-\eta}}\right)&\stackrel{!}{=}&0,
\nonumber\\
\label{we04}
\frac{1}{4}+\frac{\eta}{8}-\frac{B^2}{4\Omega^2}+\frac{1}{2}\sqrt{\frac{\eta}{2\pi\Omega}}
\frac{1}{\sqrt{1-\eta}}{\rm ln}\,\frac{1-\sqrt{1-\eta}}{1+\sqrt{1-\eta}}&\stackrel{!}{=}&0.\end{aligned}$$ Expanding the variational parameters into perturbation series of the square magnetic field $B^2$, $$\label{we05}
\eta(B)=\sum\limits_{n=0}^{\infty}\,\eta_n B^{2n},\qquad 
\Omega(B)=\sum\limits_{n=0}^{\infty}\,\Omega_n B^{2n}$$ and inserting these expansions into the self-consistency conditions (\[we03\]) and (\[we04\]) we obtain order by order the coefficients given in Table \[tab1\]. Inserting these values into the expression for the binding energy (\[we00\]) and expand with respect to $B^2$, we obtain the perturbation series $$\label{we06}
\varepsilon^{(1)}(B)=\frac{B}{2}-\sum\limits_{n=0}^\infty\,\varepsilon_n B^{2n}.$$ The first coefficients are also given in Table \[tab1\]. We find thus the important result that the first-order variational perturbation solution possesses a perturbative behaviour with respect to the square magnetic field strength $B^2$ in the weak-field limit thus yielding the correct asymptotics. The coefficients differ in higher order from the exact ones but are improved in higher orders of the variational perturbation theory [@PI Chap. 5].

  ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
                  $n$                                         0                                                    1                                                             2                                                                       3
  ----------------------------------- ------------------------------------------------- -------------------------------------------------------- ----------------------------------------------------------------- -----------------------------------------------------------------------------
               $\eta_n$                                     $1.0$                        ${\displaystyle}-\frac{405\pi^2}{7168}\approx -0.5576$   ${\displaystyle}\frac{16828965\pi^4}{1258815488}\approx 1.3023$   ${\displaystyle}-\frac{3886999332075\pi^6}{884272562962432}\approx -4.2260$

          \[1mm\] $\Omega_n$           ${\displaystyle}\frac{32}{9\pi}\approx 1.1318$       ${\displaystyle}\frac{99\pi}{224}\approx 1.3885$      ${\displaystyle}-\frac{1293975\pi^3}{19668992}\approx -2.03982$     ${\displaystyle}\frac{524431667187\pi^5}{27633517592576}\approx 5.8077$

        \[1mm\] $\varepsilon_n$        ${\displaystyle}-\frac{4}{3\pi}\approx -0.4244$      ${\displaystyle}\frac{9\pi}{128}\approx 0.2209$         ${\displaystyle}-\frac{8019\pi^3}{1835008}\approx -0.1355$                  ${\displaystyle}\frac{256449807\pi^5}{322256764928}
                                                                                                                                                                                                                                                  \approx 0.2435$

   \[1mm\] $\varepsilon_n$ [@cizek1]                       $-0.5$                                                $0.25$                                   ${\displaystyle}-\frac{53}{192}\approx -0.2760$                        ${\displaystyle}\frac{5581}{4608}\approx 1.2112$
  ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

  : \[tab1\] Perturbation coefficients up to order $B^6$ for the weak-field expansions of the variational parameters and the binding energy in comparison to the exact ones of Ref. [@cizek1].

### Asymptotical Behaviour in the Strong-Field Regime

In the discussion of the pure magnetic field below Eq. (\[gr03\]) we have mentioned that the variational calculation for the ground state energy which is thus associated with the zeroth Landau level yields a frequency ${\Omega_{\perp 2}}\propto B$ while ${\Omega_{\parallel}}=0$. Therefore we use the assumption (with ${\Omega_{\perp}}\equiv {\Omega_{\perp 2}}$) $$\label{st00}
{\Omega_{\perp}}\gg 2{\Omega_{\parallel}},\qquad {\Omega_{\parallel}}\ll B$$ for the consideration of the ground state energy (\[gr01\]) of the hydrogen atom in a strong magnetic field. In a first step we expand the last expression of the expectation value (\[gr02\]) which corresponds to the condition (\[st00\]) in terms of $2{\Omega_{\parallel}}/{\Omega_{\perp}}$ and reinsert this expansion in the equation of the ground state energy (\[gr01\]). Then we omit all terms proportional to $C/{\Omega_{\perp}}$ where $C$ stands for any expression with a value much smaller than the field strength $B$. In natural units, we thus obtain the strong-field approximation for the first-order binding energy (\[gr04\]) $$\label{st01}
\varepsilon_{{\Omega_{\perp}},{\Omega_{\parallel}}}^{(1)}=\frac{B}{2}-\left(\frac{{\Omega_{\perp}}}{4}+\frac{B^2}{4{\Omega_{\perp}}}+\frac{{\Omega_{\parallel}}}{4}
+\sqrt{\frac{{\Omega_{\parallel}}}{\pi}}\,{\rm ln}\,\frac{{\Omega_{\parallel}}}{2{\Omega_{\perp}}}\right).$$ As usual, we consider the zeros of the derivatives with respect to the variational parameters $$\label{st02}
\frac{\partial \varepsilon_{{\Omega_{\perp}},{\Omega_{\parallel}}}^{(1)}}{\partial {\Omega_{\parallel}}}\stackrel{!}{=}0,\qquad
\frac{\partial \varepsilon_{{\Omega_{\perp}},{\Omega_{\parallel}}}^{(1)}}{\partial {\Omega_{\perp}}}\stackrel{!}{=}0,$$ which lead to the self-consistence equations $$\begin{aligned}
\label{st03}
\sqrt{{\Omega_{\parallel}}}&=&-\frac{2}{\sqrt{\pi}}\left({\rm ln}\,{\Omega_{\parallel}}-{\rm ln}\,{\Omega_{\perp}}+2-{\rm ln}\,2 \right),\\
\label{st04}
{\Omega_{\perp}}&=&2\sqrt{\frac{{\Omega_{\parallel}}}{\pi}}+B\sqrt{1+4\frac{{\Omega_{\parallel}}}{\pi B^2}}\end{aligned}$$ Let us first consider the last equation. Utilizing the second of the conditions (\[st00\]) we expand the second root around unity yielding the expression $$\label{st05}
{\Omega_{\perp}}=B+2\sqrt{\frac{{\Omega_{\parallel}}}{\pi}}+2 \frac{{\Omega_{\parallel}}}{\pi B}-4\frac{{\Omega_{\parallel}}^2}{\pi^2 B^3}+\ldots,$$ where the terms are sorted with regard to their contribution starting with the biggest. Since we are interested in the strong $B$ limit, we can obviously neglect terms suppressed by powers of $1/B$. Thus we only consider the following terms for the moment: $$\label{st06}
{\Omega_{\perp}}\approx B+2\sqrt{\frac{{\Omega_{\parallel}}}{\pi}}.$$ Inserting this into the other condition (\[st03\]), expanding the corresponding logarithm, and, once more, neglecting terms of order $1/B$, we find $$\label{st07}
\sqrt{{\Omega_{\parallel}}}\approx\frac{2}{\sqrt{\pi}}\left({\rm ln}\,B-{\rm ln}\,{\Omega_{\parallel}}^{(1)}+{\rm ln}\,2-2 \right).$$ To obtain a tractable approximation for ${\Omega_{\parallel}}$, we perform some iterations starting from $$\label{st08}
\sqrt{{\Omega_{\parallel}}^{(1)}}=\frac{2}{\sqrt{\pi}}{\rm ln}\,2Be^{-2}$$ Reinserting this on the right-hand side of Eq. (\[st07\]), one obtains the second iteration $\sqrt{{\Omega_{\parallel}}^{(2)}}$. We stop this procedure after an additional reinsertion which yields $$\label{st09}
\sqrt{{\Omega_{\parallel}}^{(3)}}=\frac{2}{\sqrt{\pi}}\left({\rm ln}\,2Be^{-2}-2{\rm ln}\left[\frac{2}{\sqrt{\pi}}
\left\{{\rm ln}\,2Be^{-2}-2{\rm ln}\,\left(\frac{2}{\sqrt{\pi}}{\rm ln}\,2Be^{-2}\right) 
\right\}\right] \right).$$ The reader may convince himself that this iteration procedure indeed converges. For a subsequent systematical extraction of terms essentially contributing to the binding energy, the expression (\[st09\]) is not satisfactory. Therefore it is better to separate the leading term in the curly brackets and expand the logarithm of the remainder. Then this proceeding is applied to the expression in the angular brackets and so on. Neglecting terms of order ${\rm ln}^{-3} B$, we obtain $$\label{st10}
\sqrt{{\Omega_{\parallel}}^{(3)}}\approx \frac{2}{\sqrt{\pi}}\left({\rm ln}\,2Be^{-2}+{\rm ln}\frac{\pi}{4} 
-2{\rm ln}{\rm ln}\,2Be^{-2}\right).$$ The double-logarithmic term can be expanded in a similar way as described above: $$\label{st11}
{\rm ln}{\rm ln}\,2Be^{-2}={\rm ln}\left[{\rm ln}\,B\left(1+\frac{{\rm ln}\,2-2}{{\rm ln}\,B} 
\right) \right]= {\rm ln}{\rm ln}\,B+\frac{{\rm ln}\,2-2}{{\rm ln}\,B}-\frac{1}{2}
\frac{({\rm ln}\,2-2)^2}{{\rm ln}^2B}+{\cal O}({\rm ln}^{-3} B).$$ Thus the expression (\[st10\]) may be rewritten as $$\label{st12}
\sqrt{{\Omega_{\parallel}}^{(3)}}= \frac{2}{\sqrt{\pi}}\left( {\rm ln}\,B -2{\rm ln}{\rm ln}\,B+
\frac{2a}{{\rm ln}\,B}+\frac{a^2}{{\rm ln}^2B}+b \right) +{\cal O}({\rm ln}^{-3} B)$$ with abbreviations $$\label{st13}
a=2-{\rm ln}\,2 \approx 1.307,\qquad b= {\rm ln}\frac{\pi}{2}-2\approx -1.548.$$ The first observation is that the variational parameter ${\Omega_{\parallel}}$ is always much smaller than ${\Omega_{\perp}}$ in the high $B$-field limit. Thus we can further simplify the approximation (\[st06\]) by replacing $$\label{st14}
{\Omega_{\perp}}\approx B\left(1+\frac{2}{B}\sqrt{\frac{{\Omega_{\parallel}}}{\pi}}\right)\longrightarrow B$$ without affecting the following expression for the binding energy. Inserting the solutions (\[st12\]) and (\[st14\]) into the equation for the binding energy (\[st01\]) and expanding the logarithmic term once more as described, we find up to the order ${\rm ln}^{-2}B$: $$\begin{aligned}
\label{st15}
\varepsilon^{(1)}(B)&=&\frac{1}{\pi}\left({\rm ln}^2B -4\,{\rm ln}\,B\; {\rm ln}{\rm ln}\,B
+4\,{\rm ln}^2{\rm ln}\,B -4b\,{\rm ln}{\rm ln}\,B
+2(b+2)\,{\rm ln}\,B
+b^2
-\frac{1}{{\rm ln}\,B}\left[8\,{\rm ln}^2{\rm ln}\,B-8b\,{\rm ln}{\rm ln}\,B+2b^2 
\right]\right)\nonumber\\
&&+{\cal O}({\rm ln}^{-2}B)\end{aligned}$$ Note that the prefactor $1/\pi$ of the leading ${\rm ln}^2B$-term differs from a value $1/2$ obtained by Landau and Lifschitz [@landau]. Our different value is a consequence of using a harmonic trial system. The calculation of higher orders in variational perturbation theory would improve the value of the prefactor.

At a magnetic field strength $B=10^5 B_0$, which corresponds to $2.35\times 10^{10}\,{\rm T}=
2.35\times10^{14}\,{\rm G}$, the contribution from the first six terms is $22.87\,[2\,{\rm Ryd}]$. The next three terms suppressed by a factor ${\rm ln}^{-1} B$ contribute $-2.29\,[2\,{\rm Ryd}]$, while an estimate for the ${\rm ln}^{-2} B$-terms yields nearly $-0.3\,[2\,{\rm Ryd}]$. Thus we find $$\label{st16}
\varepsilon^{(1)}(10^5)=20.58\pm 0.3\,[2\,{\rm Ryd}].$$ This is in very good agreement with the value $20.60\,[2\,{\rm Ryd}]$ obtained from the full treatment described in Sec. \[vptlow\].

Table \[asymp\] lists the values of the first six terms of Eq. (\[st15\]). This shows in particular the significance of the second-leading term $-(4/\pi){\rm ln}\,B\;{\rm ln}{\rm ln}\,B$, which is of the same order of the leading term $(1/\pi){\rm ln}^2B$ but with an opposite sign. In Fig. \[bofb\], we have plotted the expression $$\label{st17}
\varepsilon_L(B)=\frac{1}{2}\,{\rm ln}^2B$$ from Landau and Lifschitz [@landau] to illustrate that it gives far too large binding energies even at very large magnetic fields, e.g. at $2000 B_0\propto 10^{12}\,{\rm G}$.

This strength of magnetic field appears on surfaces of neutron stars ($10^{10}-10^{12}\,{\rm G}$). A recently discovered new type of neuton star is the so-called magnetar. In these, charged particles such as protons and electrons produced by decaying neutrons give rise to the giant magnetic field of $10^{15}\,{\rm G}$. Magnetic fields of white dwarfs reach only up to $10^6-10^8\,{\rm G}$. All these magnetic field strengths are far from realization in experiments. The strongest magnetic fields ever produced in a laboratory were only of the order $10^5\,{\rm G}$, an order of magnitude larger than the fields in sun spots which reach about $0.4\times10^4\,{\rm G}$. Recall, for comparison, that the earth’s magnetic field has the small value of $0.6\,{\rm G}$.

As we see in Fig. \[bofb\], the nonleading terms in Eq. (\[st15\]) give important contributions to the asymptotic behaviour even at such large magnetic fields. It is an unusual property of the asymptotic behaviour that the absolute value of the difference between the Landau-expression (\[st17\]) and our approximation (\[st15\]) diverges with increasing magnetic field strengths $B$, only the relative difference decreases.

   $(1/\pi){\rm ln}^2B$   $-(4/\pi){\rm ln}\,B\;{\rm ln}{\rm ln}\,B$   $(4/\pi)\,{\rm ln}^2{\rm ln}\,B$   $-(4b/\pi)\,{\rm ln}{\rm ln}\,B$   $[2(b+2)/\pi]\,{\rm ln}\,B$   $b^2/\pi$
  ---------------------- -------------------------------------------- ---------------------------------- ---------------------------------- ----------------------------- -----------
        $42.1912$                         $-35.8181$                               $7.6019$                           $4.8173$                        $3.3098$             $0.7632$

  : \[asymp\] Example for the competing leading six terms in Eq. (\[st15\]) at $B=10^5B_0\approx 2.35\times 10^{14}\,{\rm G}$.

Summary
=======

We have calculated the effective classical potential for the hydrogen atom in a magnetic field. For this we have generalized variational perturbation theory to make it applicable to physical systems with uniform external magnetic field.

The effective classical potential containing the complete quantum statistical information of the system was determined in first-order variational perturbation theory. For zero-temperature, it gave the energy of the system. Our result consists of a single analytic expression which is quite accurate at all temperatures and magnetic field strengths.

Acknowledgments {#acknowledgments .unnumbered}
===============

We thank Prof. J. Čížek and Dr. J. Weniger for useful hints and references of a perturbative treatment of the ground state properties of the hydrogen atom in a magnetic field. The authors also thank Dr. J. Ortner and M. Steinberg for discussions of the finite temperature behaviour of this system. For interesting discussions we would also like to thank Prof. J.T. Devreese and Prof. G. Wunner.

One of us (M.B.) is grateful for support by the Studienstiftung des deutschen Volkes.

Generating Functional for Particle in Magnetic Field and Harmonic Oscillator Potential {#appgen}
======================================================================================

For the determination of the correlation functions of a system, we need to know the solution of the two-dimensional generating functional in the presence of an external source ${\bf j}=(j_x,j_y)$: $$\label{app00}
Z^{{\bf x}_0}[{\bf j}]=\lambda_{\rm th}^2\oint {\cal D}^2x\,
\delta({\bf x}_0-\overline{{\bf x}(\tau)})\,
e^{-{\cal A}^{{\bf x}_0}[{\bf x};{\bf j}]/\hbar}.$$ The action of a particle in a magnetic field in $z$-direction and a harmonic oscillator reads $$\label{app01}
{\cal A}^{{\bf x}_0}[{\bf x};{\bf j}]=\int_0^{\hbar\beta}d\tau\,\left[
\frac{M}{2}\dot{\bf x}^2(\tau)-\frac{i}{2}M\omega([{\bf x}(\tau)-{\bf x}_0]
\times\dot{\bf x}(\tau))_z 
+\frac{1}{2}M\Omega^2[{\bf x}(\tau)-{\bf x}_0]^2+{\bf j}(\tau)\cdot({\bf x}(\tau)-{\bf x}_0)  
\right].$$ The position dependent terms are centered around ${\bf x}_0=(x_0,y_0)$, which is the temporal average of the path ${\bf x}(\tau)$, and thus equal to the zero frequency component of the Fourier path, is $$\label{app02}
{\bf x}(\tau)={\bf x}_0+\sum\limits_{m=1}^\infty\,\left({\bf x}_me^{i\omega_m\tau}+
{\bf x}^\star_me^{-i\omega_m\tau} \right)$$ with the Matsubara frequencies $\omega_m=2\pi m/\hbar\beta$ and complex Fourier coefficients ${\bf x}_m={\bf x}_m^{\rm re}+i{\bf x}_m^{\rm im}$. Introducing a similar Fourier decomposition for the current ${\bf j}(\tau)$ with Fourier components ${\bf j}_m$ and using the orthonormality relation $$\label{app03}
\frac{1}{\hbar\beta}\int_0^{\hbar\beta}d\tau\,e^{i(\omega_m-\omega_n)\tau}=\delta_{m\,n},$$ the generating functional can be written as $$\label{app04}
Z^{{\bf x}_0}[{\bf j}]=\prod\limits_{m=1}^\infty\left[\int\frac{dx_m^{\rm re}dx_m^{\rm im}
dy_m^{\rm re}dy_m^{\rm im}}{(\pi/M\beta\omega_m^2)^2} 
e^{-{\cal A}_m({\bf x}_m,{\bf x}^\star_m;{\bf j}_m,{\bf j}^\star_m)/\hbar} \right]$$ with $$\begin{aligned}
\label{app05}
{\cal A}_m({\bf x}_m,{\bf x}^\star_m;{\bf j}_m,{\bf j}^\star_m)&=&
\hbar\beta M(\omega_m^2+\Omega^2)
([x_m^{\rm re}]^2+[x_m^{\rm im}]^2+[y_m^{\rm re}]^2+[y_m^{\rm im}]^2)
+2i\hbar\beta M \omega\omega_m(x_m^{\rm re}y_m^{\rm im}-x_m^{\rm im}y_m^{\rm re})\nonumber\\
&&+2\hbar\beta(x_m^{\rm re}{j_x}_m^{\rm re}+x_m^{\rm im}{j_x}_m^{\rm im}+
y_m^{\rm re}{j_y}_m^{\rm re}+y_m^{\rm im}{j_y}_m^{\rm im}).\end{aligned}$$ Expression (\[app04\]) is equivalent to the path integral (\[app00\]) and we obtain after performing the integrations and retransforming the currents $$\label{app06}
{\bf j}_m=\frac{1}{\hbar\beta}\int_0^{\hbar\beta}d\tau\,{\bf j}(\tau)e^{-i\omega_m\tau}$$ the resulting generating functional $$\label{app07}
Z^{{\bf x}_0}[{\bf j}]=Z^{{\bf x}_0}\exp\left\{\frac{1}{2\hbar^2}
\int_0^{\hbar\beta}d\tau \int_0^{\hbar\beta}d\tau'\,
{\bf j}(\tau){\bf G}^{{\bf x}_0}(\tau,\tau'){\bf j}(\tau') \right\}$$ with the partition function $$\label{app08}
Z^{{\bf x}_0}\equiv Z^{{\bf x}_0}[0]=
\prod\limits_{m=1}^\infty\,\frac{\omega_m^4}{\omega^2\omega_m^2+
(\omega_m^2+\Omega^2)^2}$$ and the $2\times 2$-matrix of Green functions $$\label{app09}
{\bf G}^{{\bf x}_0}(\tau,\tau')=\left(\begin{array}{cc}
G_{xx}^{{\bf x}_0}(\tau,\tau') & G_{xy}^{{\bf x}_0}(\tau,\tau')\\
G_{yx}^{{\bf x}_0}(\tau,\tau') & G_{yy}^{{\bf x}_0}(\tau,\tau') \end{array}\right).$$ The elements of this matrix are position-position correlation functions what can be easily proved by applying two functional derivatives with respect to the desired component of the current to the functional (\[app00\]), for example $$\label{app10}
G_{xx}^{{\bf x}_0}(\tau,\tau')={\left\langle\,(x(\tau)-x_0)\,(x(\tau')-x_0)\,\right\rangle^{{\bf x}_0}}
=\left[\hbar^2\frac{1}{Z^{{\bf x}_0}[{\bf j}]}
\frac{\delta^2}{\delta j_x(\tau)\delta j_x(\tau')}Z^{{\bf x}_0}[{\bf j}]\right]_{{\bf j}=0},$$ where we have defined expectation values by $$\label{app11}
{\left\langle\,\ldots\,\right\rangle^{{\bf x}_0}}=\frac{\lambda_{\rm th}^2}
{Z^{{\bf x}_0}}\oint {\cal D}^2x\,\ldots\,\delta({\bf x}_0-\overline{{\bf x}(\tau)})
e^{-{\cal A}^{{\bf x}_0}[{\bf x};0]/\hbar}.$$ From the above calculation we find the following expressions for the Green functions in Fourier space ($0\le\tau,\tau'\le\hbar\beta$): $$\begin{aligned}
\label{app12a}
G_{xx}^{{\bf x}_0}(\tau,\tau')&=&
{\left\langle\,{\tilde{x}}(\tau)\,{\tilde{x}}(\tau')\,\right\rangle^{{\bf x}_0}} = G_{yy}^{{\bf x}_0}(\tau,\tau')={\left\langle\,{\tilde{y}}(\tau)\,{\tilde{y}}(\tau')\,\right\rangle^{{\bf x}_0}}
\nonumber\\
&=&\frac{2}{M\beta}\sum\limits_{m=1}^\infty\,
\frac{\omega_m^2+\Omega^2}{\omega^2\omega_m^2+(\omega_m^2+\Omega^2)^2}\,
e^{-i\omega_m(\tau-\tau')},\\
\label{app12b}
G_{xy}^{{\bf x}_0}(\tau,\tau')&=&{\left\langle\,{\tilde{x}}(\tau)\,{\tilde{y}}(\tau')\,\right\rangle^{{\bf x}_0}}=
-G_{yx}^{{\bf x}_0}(\tau,\tau')=-{\left\langle\,{\tilde{y}}(\tau)\,{\tilde{x}}(\tau')\,\right\rangle^{{\bf x}_0}}\nonumber\\
&=&\frac{2\omega}{M\beta}\sum\limits_{m=1}^\infty\,
\frac{\omega_m}{\omega^2\omega_m^2+(\omega_m^2+\Omega^2)^2}\,
e^{-i\omega_m(\tau-\tau')},\end{aligned}$$ where, for simplicity, ${\tilde{\bf x}}(\tau)={\bf x}(\tau)-{\bf x}_0$. It is desirable to find analytical expressions for the Green functions and the partition function (\[app08\]). All these quantities possess the same dominator which can be decomposed as $$\label{app13}
\omega^2\omega_m^2+(\omega_m^2+\Omega^2)^2=(\omega_m^2+\Omega_+^2)(\omega_m^2+\Omega_-^2)$$ with frequencies $$\label{app14}
\Omega_\pm(\omega,\Omega)=
\sqrt{\Omega^2+\frac{1}{2}\omega^2\pm\omega\sqrt{\Omega^2+\frac{1}{4}\omega^2}}.$$ Therefore the partition function (\[app08\]) can be split into two products, each of which known from the harmonic oscillator [@PI Chap. 5]: $$\label{app15}
Z^{{\bf x}_0}=\prod\limits_{m=1}^\infty\left[\frac{\omega_m^2}{\omega_m^2+\Omega_+^2}\right]
\prod\limits_{m=1}^\infty\left[\frac{\omega_m^2}{\omega_m^2+\Omega_-^2}\right]
=\frac{\hbar\beta\Omega_+/2}{\sinh{\hbar\beta\Omega_+/2}}\,
\frac{\hbar\beta\Omega_-/2}{\sinh{\hbar\beta\Omega_-/2}}.$$ Now we apply the property (\[app13\]) to decompose the Green functions (\[app12a\]) into partial fractions, yielding $$\label{app16}
G_{xx}^{{\bf x}_0}(\tau,\tau')=G_{yy}^{{\bf x}_0}(\tau,\tau')=\frac{1}{M\beta}
\left(\alpha_1\,\sum\limits_{m=-\infty}^\infty\,\frac{1}{\omega_m^2+\Omega_+^2}\,
e^{-i\omega_m(\tau-\tau')}+\alpha_2\,\sum\limits_{m=-\infty}^\infty\,
\frac{1}{\omega_m^2+\Omega_-^2}\,
e^{-i\omega_m(\tau-\tau')}-\frac{1}{\Omega^2} \right)$$ with coefficients $$\label{app17}
\alpha_1=\frac{\Omega_+^2-\Omega^2}{\Omega_+^2-\Omega_-^2},\qquad 
\alpha_2=-\frac{\Omega_-^2-\Omega^2}{\Omega_+^2-\Omega_-^2}.$$ Following Ref. [@PI Chap. 3], sums of the kind occuring in expression (\[app16\]) are spectral decompositions of the correlation function for the harmonic oscillator and can be summed up: $$\label{app18}
\sum\limits_{m=-\infty}^\infty\,\frac{1}{\omega_m^2+\Omega_\pm^2}\,e^{-i\omega_m(\tau-\tau')}=
\frac{\hbar\beta}{2\Omega_\pm}
\frac{\cosh{\Omega_\pm(|\tau-\tau'|-\hbar\beta/2)}}{\sinh{\hbar\beta\Omega_\pm/2}}.$$ Thus, the $xx$- and $yy$-correlation functions can be expressed by $$\begin{aligned}
\label{app19}
&&G_{xx}^{{\bf x}_0}(\tau,\tau')=G_{yy}^{{\bf x}_0}(\tau,\tau')\nonumber\\
&&\hspace*{5mm}=\frac{1}{M\beta}\left(\frac{\hbar\beta}{2\Omega_+}\,
\frac{\Omega_+^2-\Omega^2}{\Omega_+^2-\Omega_-^2}\,
\frac{\cosh{\Omega_+(|\tau-\tau'|-\hbar\beta/2)}}{\sinh{\hbar\beta\Omega_+/2}} 
-\frac{\hbar\beta}{2\Omega_-}\,
\frac{\Omega_-^2-\Omega^2}{\Omega_+^2-\Omega_-^2}\,
\frac{\cosh{\Omega_-(|\tau-\tau'|-\hbar\beta/2)}}{\sinh{\hbar\beta\Omega_-/2}}
-\frac{1}{\Omega^2}
\right),\end{aligned}$$ where, from Eq. (\[app14\]), $\Omega_\pm=\Omega_\pm(\omega,\Omega)$ are functions of the original frequencies $\omega$ from the magnetic field and $\Omega$ from the additional harmonic oscillator (\[app01\]). It is obvious that expression (\[app19\]) reduces to the Green function of the harmonic oscillator for $\omega\to 0$: $$\label{app19a}
\lim_{\omega\to 0} G_{ii}^{{\bf x}_0}(\tau,\tau')=\frac{1}{M\beta\Omega^2}\left(
\frac{\hbar\beta\Omega}{2}\,\frac{\cosh{\Omega(|\tau-\tau'|-\hbar\beta/2)}}
{\sinh{\hbar\beta\Omega/2}}-1\right)$$ with $i\in\{x,y\}$. In this limit, the partition function (\[app15\]) turns out to be the usual one [@PI Chap. 5] for such a harmonic oscillator $$\label{app19aa}
\lim_{\omega\to 0} Z^{{\bf x}_0}=\frac{\hbar\beta\Omega/2}{\sinh{\hbar\beta\Omega/2}}.$$ It is worth mentioning that with the last term in Green function (\[app19\]) the classical harmonic fluctuation width $$\label{app19b}
G^{\rm cl}_{xx}=\left\langle x^2\right\rangle^{\rm cl}=\frac{1}{M\beta\Omega^2}$$ is subtracted. This is the consequence of the exclusion of the zero frequency mode of the Fourier path (\[app02\]) in the generating functional (\[app00\]). The necessity to do this has already been discussed in Sect. \[effrep\]. The other terms in Eq. (\[app19\]) are those which we would have obtained [*without*]{} separation of the $x_0$-component. Thus these terms represent the quantum mechanical Green function containing all quantum as well as thermal fluctuations. It is a nice property of all Green functions discussed in this paper that $$\label{app19c}
G^{{\bf x}_0}_{xx}(\tau,\tau')=G^{\rm qm}_{xx}(\tau,\tau')-G^{\rm cl}_{xx}.$$ Such a relation exists for all other Green functions appropriately, including momentum-position correlations which we consider subsequently.

The knowledge of relation (\[app18\]) makes it quite easy to determine the algebraic expression for the mixed $xy$-correlation functions. Rewriting Eq. (\[app12b\]) as $$\label{app20}
G_{xy}^{{\bf x}_0}(\tau,\tau')=-G_{yx}^{{\bf x}_0}(\tau,\tau')
= \frac{i\omega}{M\beta}\,\frac{1}{\Omega_+^2-\Omega_-^2}\frac{\partial}{\partial \tau}
\left(\sum\limits_{m=-\infty}^\infty\,\frac{1}{\omega_m^2+\Omega_+^2}\,e^{-i\omega_m(\tau-\tau')}
+\sum\limits_{m=-\infty}^\infty\,\frac{1}{\omega_m^2+\Omega_-^2}\,e^{-i\omega_m(\tau-\tau')}
\right)$$ and applying the derivative with respect to $\tau$ to relation (\[app18\]), we obtain the following expression for the mixed Green function: $$\begin{aligned}
\label{app21}
G_{xy}^{{\bf x}_0}(\tau,\tau')&=&-G_{yx}^{{\bf x}_0}(\tau,\tau')\nonumber\\
&=&\frac{\hbar\omega}{2iM}\,\frac{1}{\Omega_+^2-\Omega_-^2}\left\{  
\Theta(\tau-\tau')[g_{\Omega_+}(\tau,\tau')-g_{\Omega_-}(\tau,\tau')]-
\Theta(\tau'-\tau)[g_{\Omega_+}(\tau',\tau)-g_{\Omega_-}(\tau',\tau)]
\right\},\end{aligned}$$ where we have used the abbreviation $$\label{app22}
g_{\Omega_\pm}(\tau,\tau')=\frac{\sinh{\Omega_\pm(\tau-\tau'-\hbar\beta/2)}}
{\sinh{\hbar\beta\Omega_\pm/2}},\qquad \tau,\tau'\in(0,\hbar\beta).$$ Note that classically $\left\langle xy\right\rangle^{\rm cl}=0$ such that Eq. (\[app19c\]) reduces to $$\label{app22a}
G_{xy}^{{\bf x}_0}(\tau,\tau')=G_{xy}^{\rm qm}(\tau,\tau').$$ The Heaviside function in Eq. (\[app21\]) is defined symmetrically: $$\label{app22b}
\Theta(\tau-\tau')=\left\{\begin{array}{cc} 1 & \quad\tau>\tau',\\ 1/2 & \quad\tau=\tau', \\
0 & \quad\tau<\tau'. \end{array}\right.$$ In the quantum mechanical limit of zero-temperature ($\beta\to \infty$), the Green function (\[app19\]) simplifies to $$\label{app23a}
\lim_{\beta\to\infty} G_{xx}^{{\bf x}_0}(\tau,\tau')=
\lim_{\beta\to\infty} G_{yy}^{{\bf x}_0}(\tau,\tau')=\frac{\hbar}{2M}
\left(\frac{1}{\Omega_+}\,\frac{\Omega_+^2-\Omega^2}{\Omega_+^2-\Omega_-^2}\,
e^{-\Omega_+|\tau-\tau'|}-\frac{1}{\Omega_-}\,\frac{\Omega_-^2-\Omega^2}{\Omega_+^2-\Omega_-^2}\,
e^{-\Omega_-|\tau-\tau'|}\right),$$ while in Eq. (\[app21\]) only $g_{\Omega_\pm}(\tau,\tau')$ changes: $$\label{app23b}
\lim_{\beta\to\infty} g_{\Omega_\pm}(\tau,\tau')=-e^{-\Omega_\pm(\tau-\tau')}.$$

Properties of Green Functions {#appprop}
=============================

In this section we list properties of the Green functions (\[app19\]) and (\[app21\]) which are important for the forthcoming consideration of the generating functional with sources coupling linearily to position or momentum in Appendix \[appgenrp\]. For all relations we suppose that $0\le\tau,\tau'\le\hbar\beta$.

General Properties
------------------

A first observation is the temporal translational invariance of the Green functions: $$\label{prop00}
G_{ij}^{{\bf x}_0}(\tau,\tau')=G_{ij}^{{\bf x}_0}(\tau-\tau'),$$ where each of the indices $i,j$ stands for $x$ or $y$, respectively. For equal times we find $$\label{prop01}
G_{ij}^{{\bf x}_0}(\tau,\tau)=\frac{1}{M\beta}
\left(\frac{\hbar\beta}{2\Omega_+}\,\frac{\Omega_+^2-\Omega^2}{\Omega_+^2-\Omega_-^2} 
\coth{\hbar\beta\Omega_+/2}-\frac{\hbar\beta}{2\Omega_-}\,
\frac{\Omega_-^2-\Omega^2}{\Omega_+^2-\Omega_-^2} \coth{\hbar\beta\Omega_-/2}-\frac{1}{\Omega^2}
\right)\times
\left\{\begin{array}{cc}
1 & \quad i=j,\\
0 & \quad i\neq j.
\end{array}\right.$$ Moreover we read off the following symmetries from the expressions (\[app19\]) and (\[app21\]): $$\label{prop02}
G_{ij}^{{\bf x}_0}(\tau,\tau')=G_{ij}^{{\bf x}_0}(\tau',\tau)\times\left\{\begin{array}{cc}
1 & \quad i=j,\\
-1 & \quad i\neq j.
\end{array}\right.$$ Otherwise, $$\label{prop03}
G_{ij}^{{\bf x}_0}(\tau,\tau')=G_{ji}^{{\bf x}_0}(\tau',\tau).$$ Throughout the paper we always use periodic paths. Hence it is obvious that all Green functions are periodic, too: $$\label{prop04}
G_{ij}^{{\bf x}_0}(0,\tau')=G_{ij}^{{\bf x}_0}(\hbar\beta,\tau'),\qquad 
G_{ij}^{{\bf x}_0}(\tau,0)=G_{ij}^{{\bf x}_0}(\tau,\hbar\beta).$$

Derivatives of Green Functions
------------------------------

We now proceed with derivatives of the Green functions (\[app19\]) and (\[app21\]), since these are essential for the derivation of the generating functional of position and momentum dependent correlations in the forthcoming Appendix \[appgenrp\].

Before considering the concrete expressions we introduce a new symbol indicating uniquely to which argument the derivative is applied. A dot on the left-hand side means to perform the derivative with respect to the first argument and the dot on the right-hand side indicates that to differentiate with respect to the other argument. Having a dot on both sides the Green function is derived with respect to both arguments: $$\label{der00}
{{^\bullet G}}_{ij}^{{\bf x}_0}(\tau,\tau')=
\frac{\partial G_{ij}^{{\bf x}_0}(\tau,\tau')}{\partial \tau},\quad
{{G^\bullet}}_{ij}^{{\bf x}_0}(\tau,\tau')=
\frac{\partial G_{ij}^{{\bf x}_0}(\tau,\tau')}{\partial \tau'},\quad
{{^\bullet G^\bullet}}_{ij}^{{\bf x}_0}(\tau,\tau')=
\frac{\partial^2 G_{ij}^{{\bf x}_0}(\tau,\tau')}{\partial \tau \partial \tau'}.$$ Applying such derivatives to the Green functions (\[app19\]), we obtain ($i\in \{x,y\}$): $$\label{der01}
{{^\bullet G}}_{ii}^{{\bf x}_0}(\tau,\tau')=\frac{\hbar}{2M}\frac{1}{\Omega_+^2-\Omega_-^2}
\left[\Theta(\tau-\tau')g(\tau,\tau') 
-\Theta(\tau'-\tau)g(\tau',\tau)\right]=-{{G^\bullet}}_{ii}^{{\bf x}_0}(\tau,\tau')$$ with $$\label{der02}
g(\tau,\tau')=(\Omega_+^2-\Omega^2)g_{\Omega_+}(\tau,\tau')
-(\Omega_-^2-\Omega^2)g_{\Omega_-}(\tau,\tau'),$$ where $g_\pm(\tau,\tau)$ was defined in Eq. (\[app22\]). Performing the derivatives to both arguments leads to the expression $$\label{der03}
{{^\bullet G^\bullet}}_{ii}^{{\bf x}_0}(\tau,\tau')={{^\bullet \tilde{G}^\bullet}}_{ii}^{{\bf x}_0}(\tau,\tau')+
\frac{\hbar}{M}\delta(\tau-\tau'),$$ where we have introduced the partial function $$\label{der04}
{{^\bullet \tilde{G}^\bullet}}_{ii}^{{\bf x}_0}(\tau,\tau')=-\frac{\hbar}{2M}\left[\Omega_+
\frac{\Omega_+^2-\Omega^2}{\Omega_+^2-\Omega_-^2}
\frac{\sinh{\Omega_+(|\tau-\tau'|-\hbar\beta/2)}}{\sinh{\hbar\beta\Omega_+/2}}-
\Omega_-
\frac{\Omega_-^2-\Omega^2}{\Omega_+^2-\Omega_-^2}
\frac{\sinh{\Omega_-(|\tau-\tau'|-\hbar\beta/2)}}{\sinh{\hbar\beta\Omega_-/2}} \right]$$ which is finite for equal times.

Applying derivatives with respect to the first respective second argument to the mixed correlation function (\[app21\]), we find: $$\label{der05}
{{^\bullet G}}_{xy}^{{\bf x}_0}(\tau,\tau')=\frac{\hbar\omega}{2iM}\frac{1}{\Omega_+^2-\Omega_-^2}
\left[\Omega_+\frac{\cosh{\Omega_+(|\tau-\tau'|-\hbar\beta/2)}}{\sinh{\hbar\beta\Omega_+/2}}
-\Omega_-\frac{\cosh{\Omega_-(|\tau-\tau'|-\hbar\beta/2)}}{\sinh{\hbar\beta\Omega_-/2}}\right]
=-{{G^\bullet}}_{xy}^{{\bf x}_0}(\tau,\tau')$$ and $$\label{der05b}
{{^\bullet G}}_{yx}^{{\bf x}_0}(\tau,\tau')=-{{^\bullet G}}_{xy}^{{\bf x}_0}(\tau,\tau').$$ Differentiating each argument of the mixed Green function results in $$\label{der06}
{{^\bullet G^\bullet}}_{xy}^{{\bf x}_0}(\tau,\tau')=\frac{i\hbar\omega}{2M}\frac{1}{\Omega_+^2-\Omega_-^2}
\left[\Theta(\tau-\tau')h(\tau,\tau')-\Theta(\tau'-\tau)h(\tau',\tau) \right]=
-{{^\bullet G^\bullet}}_{yx}^{{\bf x}_0}(\tau,\tau')$$ with $$\label{der07}
h(\tau,\tau')=\Omega_+^2g_{\Omega_+}(\tau,\tau')-\Omega_-^2g_{\Omega_-}(\tau,\tau').$$ An additional property we read off from Eqs. (\[der01\]) and (\[der05\]) is ($i,j\in \{x,y\}$): $$\begin{aligned}
\label{der08a}
{{^\bullet G}}_{ij}^{{\bf x}_0}(\tau,\tau')={{^\bullet G}}_{ij}^{{\bf x}_0}(\tau',\tau)\times
\left\{\begin{array}{cc}
-1 &\quad i=j,\\
1 & \quad i\neq j,
\end{array}\right.\\
\label{der08b}
{{G^\bullet}}_{ij}^{{\bf x}_0}(\tau,\tau')={{G^\bullet}}_{ij}^{{\bf x}_0}(\tau',\tau)\times
\left\{\begin{array}{cc}
-1 &\quad i=j,\\
1 & \quad i\neq j.
\end{array}\right.\end{aligned}$$ The double-sided derivatives (\[der03\]), (\[der04\]), and (\[der06\]) imply $$\label{der09}
{{^\bullet G^\bullet}}_{ij}^{{\bf x}_0}(\tau,\tau')={{^\bullet G^\bullet}}_{ij}^{{\bf x}_0}(\tau',\tau)\times
\left\{\begin{array}{cc}
1 &\quad i=j,\\
-1 & \quad i\neq j.
\end{array}\right.$$ The derivatives (\[der01\]), (\[der04\]), (\[der05\]), and (\[der06\]) are periodic: $$\begin{aligned}
\label{der10a}
&&{{^\bullet G}}_{ij}^{{\bf x}_0}(\tau,0)={{^\bullet G}}_{ij}^{{\bf x}_0}(\tau,\hbar\beta),\quad 
{{^\bullet G}}_{ij}^{{\bf x}_0}(0,\tau')={{^\bullet G}}_{ij}^{{\bf x}_0}(\hbar\beta,\tau'),\\
\label{der10b}
&&{{G^\bullet}}_{ij}^{{\bf x}_0}(\tau,0)={{G^\bullet}}_{ij}^{{\bf x}_0}(\tau,\hbar\beta),\quad 
{{G^\bullet}}_{ij}^{{\bf x}_0}(0,\tau')={{G^\bullet}}_{ij}^{{\bf x}_0}(\hbar\beta,\tau'),\\
\label{der10c}
&&{{^\bullet \tilde{G}^\bullet}}_{ii}^{{\bf x}_0}(\tau,0)={{^\bullet \tilde{G}^\bullet}}_{ii}^{{\bf x}_0}(\tau,\hbar\beta),\quad 
{{^\bullet \tilde{G}^\bullet}}_{ii}^{{\bf x}_0}(0,\tau')={{^\bullet \tilde{G}^\bullet}}_{ii}^{{\bf x}_0}(\hbar\beta,\tau'),\\
\label{der10d}
&&{{^\bullet G^\bullet}}_{ij}^{{\bf x}_0}(\tau,0)={{^\bullet G^\bullet}}_{ij}^{{\bf x}_0}(\tau,\hbar\beta),\quad 
{{^\bullet G^\bullet}}_{ij}^{{\bf x}_0}(0,\tau')={{^\bullet G^\bullet}}_{ij}^{{\bf x}_0}(\hbar\beta,\tau'),\qquad (i\neq j).\end{aligned}$$

Generating Functional for Position- and Momentum-Dependent Correlation Functions {#appgenrp}
================================================================================

With the discussion of the generating functional for position-dependent correlation functions and, in particular, the Green functions in Appendix \[appgen\] and their properties in Appendix \[appprop\], we have layed the foundation to derive the generating functional for correlation functions depending on both, position and momentum. Following the framework presented in an earlier work [@correlation], such a functional involving sources coupled to the momentum can always be reduced to one containing position-coupled sources only.

We start from the three-dimensional effective classical representation for the generating functional $$\label{xp00}
Z_{\bf \Omega}[{\bf j},{\bf v}]=
\int\frac{d^3x_0d^3p_0}{(2\pi\hbar)^3}Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf j},{\bf v}]$$ with zero frequency components ${\bf x}_0=(x_0,y_0,z_0)={\rm const.}$ and ${\bf p}_0=({p_x}_0,{p_y}_0,{p_y}_0)={\rm const.}$ of the Fourier path separated. The reduced functional is $$\label{xp01}
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf j},{\bf v}]=
(2\pi\hbar)^3\oint {\cal D'}^3x{\cal D}^3p\,\delta({\bf x}_0-\overline{{\bf x}(\tau)})
\delta({\bf p}_0-\overline{{\bf p}(\tau)})\,
\exp\left\{-\frac{1}{\hbar}{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
[{\bf p},{\bf x};{\bf j},{\bf v}] \right\},$$ where the path integral measure is that defined in Eq. (\[ham03\]). Extending the action (\[ham02\]) by source terms, considering a more general Hamilton function than (\[ham16\]), and introducing an additional harmonic oscillator in $z$-direction, the action functional in Eq. (\[xp01\]) shall read $$\begin{aligned}
\label{xp02}
{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf p},{\bf x};{\bf j},{\bf v}]=
\int_0^{\hbar\beta}d\tau\,\Big\{&&-i{\tilde{\bf p}}(\tau)\cdot \dot{\bf x}(\tau)+
\frac{1}{2M}{\tilde{\bf p}}^2(\tau)-\frac{1}{2}{\Omega_{\perp 1}}l_z({\tilde{\bf p}},{\tilde{\bf x}}) 
+\frac{1}{8}M{\Omega_{\perp 2}}^2\left[ 
{\tilde{x}}^2(\tau)+{\tilde{y}}^2(\tau)\right]+\frac{1}{2}M{\Omega_{\parallel}}^2{\tilde{z}}^2(\tau)\nonumber\\
&&+{\bf j}(\tau)\cdot{\tilde{\bf x}}(\tau)+{\bf v}(\tau)\cdot{\tilde{\bf p}}(\tau) 
\Big\}\end{aligned}$$ with shifted positions and momenta $$\label{xp03}
{\tilde{\bf x}}={\bf x}(\tau)-{\bf x}_0,\qquad {\tilde{\bf p}}={\bf p}(\tau)-{\bf p}_0.$$ The orbital angular momentum $l_z({\bf p},{\bf x})$ is defined in Eq. (\[ham18\]) and is used in Eq. (\[xp02\]) with the shifted phase space coordinates (\[xp03\]). We have introduced three different frequencies in (\[xp02\]), ${\bf \Omega}=({\Omega_{\perp 1}},{\Omega_{\perp 2}},{\Omega_{\parallel}})$, where the first both components are used in regard to the oscillations in the plane perpendicular to the direction of the magnetic field which shall be considered here to point into $z$-direction. The last component, ${\Omega_{\parallel}}$, is the frequency of a trial oscillator parallel to the field lines.

Due to the periodicity of the paths, we suppose that the sources might also be periodic: $$\label{xp03b}
{\bf j}(0)={\bf j}(\hbar\beta),\qquad {\bf v}(0)={\bf v}(\hbar\beta).$$ Since we want to simplify expression (\[xp01\]) such that we can use the results obtained in Appendix \[appgen\], the momentum path integral is solved in the following. In a first step we reexpress the momentum $\delta$-function in (\[xp01\]) by $$\label{xp04}
\delta({\bf p}_0-\overline{{\bf p}(\tau)})=\int\frac{d^3\xi}{(2\pi\hbar)^3}\,
\exp\left\{-\frac{1}{\hbar}\int_0^{\hbar\beta} d\tau \,{\bf v}_0\cdot[{\bf p}(\tau)-{\bf p}_0]
\right\},$$ where $$\label{xp05}
{\bf v}_0({\mbox{\boldmath$\xi$}})=\frac{i}{\hbar\beta}{\mbox{\boldmath$\xi$}}$$ is an additional current which is coupled to the momentum and is constant in time. Defining the sum of all sources coupled to the momentum by $$\label{xp06}
{\bf V}({\mbox{\boldmath$\xi$}},\tau)={\bf v}(\tau)+{\bf v}_0({\mbox{\boldmath$\xi$}}),$$ the functional (\[xp01\]) can be written as $$\begin{aligned}
\label{xp07}
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf j},{\bf v}]&=&
\int d^3\xi\oint {\cal D'}^3x{\cal D}^3p\,\delta({\bf x}_0-\overline{{\bf x}(\tau)})
\exp\Bigg\{-\frac{1}{\hbar}\int_0^{\hbar\beta}d\tau\,\Big[-i{\bf p}(\tau)
\cdot\dot{\bf x}(\tau) 
+\frac{{\bf p}^2(\tau)}{2M}-\frac{1}{2}{\Omega_{\perp 1}}l_z({\bf p}(\tau),{\tilde{\bf x}}(\tau))\nonumber\\
&&+\frac{1}{8}M{\Omega_{\perp 2}}^2\left\{{\tilde{x}}^2(\tau)+{\tilde{y}}^2(\tau)\right\}+\frac{1}{2}M{\Omega_{\parallel}}^2{\tilde{z}}^2(\tau)
+{\bf j}(\tau)\cdot{\tilde{\bf x}}(\tau)+{\bf V}({\mbox{\boldmath$\xi$}},\tau)\cdot{\bf p}(\tau)
\Big] \Bigg\},\end{aligned}$$ where we have used the translation invariance ${\tilde{\bf p}}\to {\bf p}$ of the path integral. To solve the momentum path integral, it is useful to express it in its discretized form. Performing quadratic completions such that the momentum path integral separates into an infinite product of simple Gaussian integrals which are easily calculated, the remaining functional is reduced to the configuration space path integral $$\begin{aligned}
\label{xp08}
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf j},{\bf v}]=
\int d^3\xi\,\exp\left[\frac{M}{2\hbar}\int_0^{\hbar\beta}d\tau\,{\bf V}^2({\mbox{\boldmath$\xi$}},\tau) 
\right]
\oint{\cal D}^3x\,\delta({\bf x}_0-\overline{{\bf x}(\tau)})
\exp\left\{-\frac{1}{\hbar}{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
[{\bf x};{\bf j},{\bf V}]\right\}\end{aligned}$$ with the measure (\[ham09\]) for $D=3$. The action functional is $$\begin{aligned}
\label{xp09}
{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf x};{\bf j},{\bf V}]&=&
\int_0^{\hbar\beta}d\tau\,\Bigg[\frac{M}{2}\dot{\bf x}^2(\tau)  
+\frac{1}{2}iM{\Omega_{\perp 1}}\left\{\dot{x}(\tau){\tilde{y}}(\tau)-\dot{y}(\tau){\tilde{x}}(\tau)\right\}
+\frac{1}{8}M\left({\Omega_{\perp 2}}^2-{\Omega_{\perp 1}}^2\right)\left\{{\tilde{x}}^2(\tau)+{\tilde{y}}^2(\tau)\right\}\nonumber\\
&&\hspace*{-3pt}
+\frac{1}{2}M{\Omega_{\parallel}}^2{\tilde{z}}^2(\tau)+{\tilde{x}}(\tau)\left[j_x(\tau)+\frac{1}{2}M{\Omega_{\perp 1}}V_y({\mbox{\boldmath$\xi$}},\tau)
\right]\nonumber\\
&&\hspace*{-3pt}+{\tilde{y}}(\tau)\left[j_y(\tau)-\frac{1}{2}M{\Omega_{\perp 1}}V_x({\mbox{\boldmath$\xi$}},\tau)\right]+{\tilde{z}}(\tau)j_z(\tau)
\Bigg]
-\frac{iM}{\hbar}\int_0^{\hbar\beta}d\tau\,\dot{\bf x}(\tau)\cdot
{\bf V}({\mbox{\boldmath$\xi$}},\tau),\end{aligned}$$ where the last term simplifies by the following consideration. A partial integration of this term yields $$\label{xp10}
\int_0^{\hbar\beta}d\tau\,\dot{\bf x}(\tau)\cdot{\bf V}({\mbox{\boldmath$\xi$}},\tau)
=-\int_0^{\hbar\beta}d\tau\,({\bf x}(\tau)-{\bf x}_0)\cdot\dot{\bf V}({\mbox{\boldmath$\xi$}},\tau).$$ The surface term vanishes as a consequence of the periodicity of the path and the source. This periodicity is also the reason why we could shift ${\bf x}(\tau)$ by the constant ${\bf x}_0$ on the right-hand side of Eq. (\[xp10\]). Obviously, the importance of this expression lies in the coupling of the time derivative of ${\bf V}({\mbox{\boldmath$\xi$}},\tau)$ to the path ${\bf x}(\tau)$. Thus, $\dot{\bf V}({\mbox{\boldmath$\xi$}},\tau)$ can be handled like a ${\bf j}(\tau)$-current [@correlation] and the action (\[xp09\]) can be written as $$\label{xp11}
{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf x};{\bf j},{\bf V}]=
{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf x};{\bf J},0]=
{\cal A}^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf x};0,0]-\frac{1}{\hbar}
\int_0^{\hbar\beta}d\tau\,{\tilde{\bf x}}(\tau)\cdot{\bf J}({\mbox{\boldmath$\xi$}},\tau)$$ with the new current vector ${\bf J}({\mbox{\boldmath$\xi$}},\tau)$ which has the components $$\begin{aligned}
\label{xp12}
J_x({\mbox{\boldmath$\xi$}},\tau)&=&j_x(\tau)+\frac{1}{2}M{\Omega_{\perp 1}}V_y({\mbox{\boldmath$\xi$}},\tau)-
iM\dot{V}_x({\mbox{\boldmath$\xi$}},\tau),\nonumber\\ 
J_y({\mbox{\boldmath$\xi$}},\tau)&=&j_y(\tau)-\frac{1}{2}M{\Omega_{\perp 1}}V_x({\mbox{\boldmath$\xi$}},\tau)-
iM\dot{V}_y({\mbox{\boldmath$\xi$}},\tau),\\
J_z({\mbox{\boldmath$\xi$}},\tau)&=&j_z(\tau)-\frac{1}{2}M{\Omega_{\parallel}}V_z({\mbox{\boldmath$\xi$}},\tau)\nonumber.\end{aligned}$$ and couples to the path ${\bf x}(\tau)$ only. With the expression (\[xp08\]) for the generating functional and the action (\[xp11\]), we have derived a representation similar to Eq. (\[app00\]) with the action (\[app01\]), extended by an additional oscillator in $z$-direction. We identify $$\label{xp13}
\omega\equiv {\Omega_{\perp 1}},\quad \Omega\equiv \frac{1}{4}\left({\Omega_{\perp 2}}^2-{\Omega_{\perp 1}}^2\right),
\quad j_x\equiv J_x,\quad
j_y\equiv J_y.$$ Thus the auxiliary frequencies $\Omega_{\pm}$ (\[app14\]) become $$\label{xp13b}
\Omega_{\pm}({\Omega_{\perp 1}},{\Omega_{\perp 2}})=\frac{1}{2}\,|{\Omega_{\perp 1}}\pm{\Omega_{\perp 2}}|.$$ Inserting the substitutions (\[xp13\]) into the solution (\[app07\]) for the generating functional in two dimensions and performing the usual calculation for a harmonic oscillator with external source [@PI Chaps. 3,5] in $z$-direction, we obtain an intermediate result for the generating functional in three dimensions (\[xp01\]): $$\label{xp14}
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf j},{\bf v}]=
\lambda_{\rm th}^{-3}\,Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
\int d^3\xi\,\exp\left\{\frac{M}{2\hbar}\int_0^{\hbar\beta}d\tau\,{\bf V}^2({\mbox{\boldmath$\xi$}},\tau) 
\right\}\,\exp\left\{\frac{1}{2\hbar^2}\int_0^{\hbar\beta}d\tau\int_0^{\hbar\beta}d\tau'\,
{\bf J}({\mbox{\boldmath$\xi$}},\tau)\,{\bf G}^{{\bf x}_0}(\tau,\tau')\,{\bf J}({\mbox{\boldmath$\xi$}},\tau')\right\}.$$ The partition function follows from Eqs. (\[app15\]) and (\[app19aa\]) $$\label{xp15}
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}=Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[0,0]=
\frac{\hbar\beta\Omega_+/2}{\sinh{\hbar\beta\Omega_+/2}}\,
\frac{\hbar\beta\Omega_-/2}{\sinh{\hbar\beta\Omega_-/2}}\,
\frac{\hbar\beta{\Omega_{\parallel}}/2}{\sinh{\hbar\beta{\Omega_{\parallel}}/2}}$$ and ${\bf G}^{{\bf x}_0}(\tau,\tau')$ is the $3\times 3$-matrix of Green functions $$\label{xp16}
{\bf G}^{{\bf x}_0}(\tau,\tau')=\left(\begin{array}{ccc}
G_{xx}^{{\bf x}_0}(\tau,\tau') & G_{xy}^{{\bf x}_0}(\tau,\tau') & 0\\
G_{yx}^{{\bf x}_0}(\tau,\tau') & G_{yy}^{{\bf x}_0}(\tau,\tau') & 0\\
0 & 0 & G_{zz}^{{\bf x}_0}(\tau,\tau')\end{array}\right).$$ Except $G_{zz}^{{\bf x}_0}(\tau,\tau')$, the Green functions are given by the expressions in Eqs. (\[app19\]) and (\[app21\]) with frequencies (\[xp13b\]). The Green function of the pure harmonic oscillator in $z$-direction $$\label{xp17}
G_{zz}^{{\bf x}_0}(\tau,\tau')=\frac{1}{M\beta{\Omega_{\parallel}}^2}\left(
\frac{\hbar\beta{\Omega_{\parallel}}}{2}\,\frac{\cosh{{\Omega_{\parallel}}(|\tau-\tau'|-\hbar\beta/2)}}
{\sinh{\hbar\beta{\Omega_{\parallel}}/2}}-1\right)$$ follows directly from the limit (\[app19a\]). Since the current ${\bf J}$ (\[xp12\]) still depends on time derivatives of ${\bf V}$, we have to perform some partial integrations in the functional (\[xp14\]). This is a very extensive but straightforward work and thus we only present an instructive example. For that we apply the properties and the time derivatives of the Green functions which we presented in Appendix \[appprop\]. Consider the integral $$\label{xp18}
I=-\frac{M^2}{2\hbar^2}\int_0^{\hbar\beta}d\tau\int_0^{\hbar\beta}d\tau'\,\dot{V}_i
({\mbox{\boldmath$\xi$}},\tau)\,G_{ii}^{{\bf x}_0}(\tau,\tau')\,\dot{V}_i({\mbox{\boldmath$\xi$}},\tau')$$ occuring in the second exponential of Eq. (\[xp14\]) with $i\in\{x,y,z\}$. A partial integration in the $\tau'$-integral leads to $$\begin{aligned}
\label{xp18a}
I&=&-\frac{M^2}{2\hbar^2}\int_0^{\hbar\beta}d\tau\,\dot{V}_i({\mbox{\boldmath$\xi$}},\tau)
\left(G_{ii}^{{\bf x}_0}(\tau,\tau')\,V_i({\mbox{\boldmath$\xi$}},\tau')\Big|_{\tau'=0}^{\tau'=\hbar\beta}
-\int_0^{\hbar\beta}d\tau'\,\frac{\partial G_{ii}^{{\bf x}_0}(\tau,\tau')}{\partial \tau'}
\,V_i({\mbox{\boldmath$\xi$}},\tau')\right)\nonumber\\
&=&\frac{M^2}{2\hbar^2}\int_0^{\hbar\beta}d\tau\int_0^{\hbar\beta}d\tau'\,
\dot{V}_i({\mbox{\boldmath$\xi$}},\tau)\,{{G^\bullet}}_{ii}^{{\bf x}_0}(\tau,\tau')\,V_i({\mbox{\boldmath$\xi$}},\tau').\end{aligned}$$ The surface term in the first line vanishes as a consequence of the periodicity of the current (\[xp03b\]) and the Green function (\[prop04\]). A second partial integration, now in the $\tau$-integral, results in $$\begin{aligned}
\label{xp18b}
I&=&-\frac{M^2}{2\hbar^2}\int_0^{\hbar\beta}d\tau\int_0^{\hbar\beta}d\tau'\,
V_i({\mbox{\boldmath$\xi$}},\tau)\,{{^\bullet G^\bullet}}_{ii}^{{\bf x}_0}(\tau,\tau')\,V_i({\mbox{\boldmath$\xi$}},\tau')\nonumber\\
&=&-\frac{M^2}{2\hbar^2}\int_0^{\hbar\beta}d\tau\int_0^{\hbar\beta}d\tau'\,
V_i({\mbox{\boldmath$\xi$}},\tau)\,{{^\bullet \tilde{G}^\bullet}}_{ii}^{{\bf x}_0}(\tau,\tau')\,V_i({\mbox{\boldmath$\xi$}},\tau')-
\frac{M}{2\hbar}\int_0^{\hbar}d\tau\,V_i^2({\mbox{\boldmath$\xi$}},\tau).\end{aligned}$$ Here we have applied the periodicity property of the right-hand derivative of the Green function (\[der10b\]), leading to a vanishing surface term in this case, too. In the second line, we have used the decomposition (\[der03\]) of the double-sided differentiated Green function. Note that the last term just cancels the appropriate term in the first exponential of the right-hand side of Eq. (\[xp14\]). Eventually, after performing all such partial integrations, we reexpress Eq. (\[xp14\]) by $$\label{xp19}
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf j},{\bf v}]=
\lambda_{\rm th}^{-3}\,Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
\int d^3\xi\,\exp\left\{\frac{1}{2\hbar^2}\int_0^{\hbar\beta}d\tau\int_0^{\hbar\beta}d\tau'\,
\tilde{\bf s}({\mbox{\boldmath$\xi$}},\tau)\,{\bf H}^{{\bf x}_0}(\tau,\tau')\,\tilde{\bf s}({\mbox{\boldmath$\xi$}},\tau')  
\right\}$$ with six-dimensional sources $$\label{xp20}
\tilde{\bf s}({\mbox{\boldmath$\xi$}},\tau)=\left(\,{\bf j}(\tau),{\bf V}({\mbox{\boldmath$\xi$}},\tau)\,\right).$$ and the $6\times 6$-matrix ${\bf H}^{{\bf x}_0}(\tau,\tau')$ which has no significance as long as we have not done the ${\mbox{\boldmath$\xi$}}$-integration. We explicitly insert the decomposition (\[xp06\]) into expression (\[xp20\]) of the source vector $\tilde{\bf s}$. Since ${\bf v}_0({\mbox{\boldmath$\xi$}})$ from Eq. (\[xp05\]) is constant in time, some temporal integrals in the exponential of Eq. (\[xp19\]) can be calculated and we obtain $$\begin{aligned}
\label{xp21}
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf j},{\bf v}]&=&
\lambda_{\rm th}^{-3}\,Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
\exp\left\{\frac{1}{2\hbar^2}\int_0^{\hbar\beta}d\tau\int_0^{\hbar\beta}d\tau'\,
{\bf s}(\tau)\,{\bf H}^{{\bf x}_0}(\tau,\tau')\,{\bf s}(\tau')\right\}\nonumber\\
&&\times\int d^3\xi\,
\exp\left\{-\frac{M}{2\hbar^2\beta}{\mbox{\boldmath$\xi$}}^2+i\frac{M}{\hbar^2\beta}{\mbox{\boldmath$\xi$}}\cdot
\int_0^{\hbar\beta}d\tau\,{\bf v}(\tau) \right\}\end{aligned}$$ with the new $6$-vector $$\label{xp22}
{\bf s}(\tau)=(\,{\bf j}(\tau),{\bf v}(\tau)\,)$$ consisting of the original sources ${\bf j}$ and ${\bf v}$ only. The Gaussian $\xi$-integral in Eq. (\[xp21\]) can be easily solved and the terms appearing from quadratic completion modify the above matrix ${\bf H}^{{\bf x}_0}(\tau,\tau')$. The final result for the generating functional of all position and momentum dependent correlations is given by $$\label{xp23}
Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}[{\bf j},{\bf v}]=Z^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}
\exp\left\{\frac{1}{2\hbar^2}\int_0^{\hbar\beta}d\tau\int_0^{\hbar\beta}d\tau'\,
{\bf s}(\tau)\,{\bf G}^{{\bf p}_0,{{\bf x}_0}}(\tau,\tau')\,{\bf s}(\tau')\right\}.$$ The complete $6\times 6$-matrix ${\bf G}^{{\bf p}_0,{{\bf x}_0}}(\tau,\tau')$ contains all possible Green functions describing position-position, position-momentum, and momentum-momentum correlations. As a consequence of separating the fluctuations into those perpendicular and parallel to the direction of the magnetic field, all correlations between $x,y$ on the one and $z$ on the other hand vanish as well as those for the appropriate momenta. The symmetries for the Green functions and their derivatives were investigated in detail in Appendix \[appprop\] and lead to a further reduction of the number of significant matrix elements. It turns out that only 9 elements are independent of each other. Therefore we can write the matrix $$\label{xp24}
{\bf G}^{{\bf x}_0,{{\bf p}_0}}(\tau,\tau')=\left(\begin{array}{cccccc}
G_{xx}^{{\bf p}_0,{\bf x}_0}(\tau,\tau') & G_{xy}^{{\bf p}_0,{\bf x}_0}(\tau,\tau') & 0 & 
G_{xp_x}^{{\bf p}_0,{\bf x}_0}(\tau,\tau') & G_{xp_y}^{{\bf p}_0,{\bf x}_0}(\tau,\tau') & 0 \\
G_{xy}^{{\bf p}_0,{\bf x}_0}(\tau',\tau) & G_{xx}^{{\bf p}_0,{\bf x}_0}(\tau,\tau') & 0 & 
-G_{xp_y}^{{\bf p}_0,{\bf x}_0}(\tau,\tau') & G_{xp_x}^{{\bf p}_0,{\bf x}_0}(\tau,\tau') & 0 \\
0 & 0 & G_{zz}^{{\bf p}_0,{\bf x}_0}(\tau,\tau') & 0 & 0 & 
G_{zp_z}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')\\
G_{xp_x}^{{\bf p}_0,{\bf x}_0}(\tau',\tau) & -G_{xp_y}^{{\bf p}_0,{\bf x}_0}(\tau',\tau) & 0 & 
G_{p_xp_x}^{{\bf p}_0,{\bf x}_0}(\tau,\tau') & G_{p_xp_y}^{{\bf p}_0,{\bf x}_0}(\tau,\tau') & 0 \\
G_{xp_y}^{{\bf p}_0,{\bf x}_0}(\tau',\tau) & G_{xp_x}^{{\bf p}_0,{\bf x}_0}(\tau',\tau) & 0 & 
G_{p_xp_y}^{{\bf p}_0,{\bf x}_0}(\tau',\tau) & G_{p_xp_x}^{{\bf p}_0,{\bf x}_0}(\tau,\tau') & 0 \\
0 & 0 & G_{zp_z}^{{\bf p}_0,{\bf x}_0}(\tau',\tau) & 0 & 0 & 
G_{p_zp_z}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')
\end{array}\right).$$ The matrix decays into four $3\times 3$-blocks, each of the which describing another type of correlation: the upper left position-position, the upper right position-momentum (as well as the lower left one), and the lower right momentum-momentum correlations. The different elements of the matrix are $$\begin{aligned}
\label{xp25a}
G_{xx}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')&=&{\left\langle\,{\tilde{x}}(\tau){\tilde{x}}(\tau')\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=
G_{xx}^{{\bf x}_0}(\tau,\tau'),\\
\label{xp25b}
G_{xy}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')&=&{\left\langle\,{\tilde{x}}(\tau){\tilde{y}}(\tau')\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=
G_{xy}^{{\bf x}_0}(\tau,\tau'),\\
\label{xp25c}
G_{zz}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')&=&{\left\langle\,{\tilde{z}}(\tau){\tilde{z}}(\tau')\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=
G_{zz}^{{\bf x}_0}(\tau,\tau'),\\ 
\label{xp25d}
G_{xp_x}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')&=&{\left\langle\,{\tilde{x}}(\tau){\tilde{p}_x}(\tau')\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=
iM{{G^\bullet}}_{xx}^{{\bf x}_0}(\tau,\tau')-\frac{1}{2}M{\Omega_{\perp 1}}G_{xy}^{{\bf x}_0}(\tau,\tau'),\\
\label{xp25e}
G_{xp_y}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')&=&{\left\langle\,{\tilde{x}}(\tau){\tilde{p}_y}(\tau')\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=
iM{{G^\bullet}}_{xy}^{{\bf x}_0}(\tau,\tau')+\frac{1}{2}M{\Omega_{\perp 1}}G_{xx}^{{\bf x}_0}(\tau,\tau'),\\
\label{xp25f}
G_{zp_z}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')&=&{\left\langle\,{\tilde{z}}(\tau){\tilde{p}_z}(\tau')\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=
iM{{G^\bullet}}_{zz}^{{\bf x}_0}(\tau,\tau'),\\
\label{xp25g}
G_{p_xp_x}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')&=&{\left\langle\,{\tilde{p}_x}(\tau){\tilde{p}_x}(\tau')\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=
-M^2{{^\bullet \tilde{G}^\bullet}}_{xx}^{{\bf x}_0}(\tau,\tau')-iM^2{\Omega_{\perp 1}}{{^\bullet G}}_{xy}^{{\bf x}_0}(\tau,\tau')+\frac{1}{4}
M^2{\Omega_{\perp 1}}^2G_{xx}^{{\bf x}_0}(\tau,\tau')-\frac{M}{\beta},\\
\label{xp25h}
G_{p_xp_y}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')&=&{\left\langle\,{\tilde{p}_x}(\tau){\tilde{p}_y}(\tau')\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=
iM^2{{^\bullet G}}_{xx}^{{\bf x}_0}(\tau,\tau')-M^2{{^\bullet G^\bullet}}_{xy}^{{\bf x}_0}(\tau,\tau')+\frac{1}{4}
M^2{\Omega_{\perp 1}}^2G_{xy}^{{\bf x}_0}(\tau,\tau'),\\
\label{xp25i}
G_{p_zp_z}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')&=&{\left\langle\,{\tilde{p}_z}(\tau){\tilde{p}_z}(\tau')\,\right\rangle^{{\bf p}_0,{\bf x}_0}_{\bf \Omega}}=
-M^2{{^\bullet \tilde{G}^\bullet}}_{zz}^{{\bf x}_0}(\tau,\tau')-\frac{M}{\beta},\end{aligned}$$ where the expectation values are defined by Eq. (\[vpt07\]). Note that all these Green functions are invariant under time translations such that $$\label{xp26}
G_{\mu\nu}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')=G_{\mu\nu}^{{\bf p}_0,{\bf x}_0}(\tau-\tau')$$ with $\mu,\nu\in \{x,y,z,p_x,p_y,p_z\}$.

It is quite instructive to prove that all these Green functions can be decomposed into a quantum statistical and a classical part as we did it in Eq. (\[app19\]). Since we know that the classical correlation functions do not depend on the euclidean time, all derivative terms in Eqs. (\[xp25a\])–(\[xp25i\]) do not contribute a classical term. We can write each Green function $$\label{xp27}
G_{\mu\nu}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')=G_{\mu\nu}^{\rm qm}(\tau,\tau')-
G_{\mu\nu}^{\rm cl},$$ This relation has been already checked for Eqs. (\[xp25a\])-(\[xp25c\]) in Appendix \[appgen\]. The classical contribution is zero in Eqs. (\[xp25d\]), (\[xp25f\]), and (\[xp25h\]) following from the absence of classical terms in derivatives of the Green functions and mixed correlations like (\[app22a\]). It seems surprising that the correlation (\[xp25e\]) contains a classical term while (\[xp25d\]) possesses none. This is, however, a consequence of the cross product of the orbital angular momentum appearing in the action (\[xp02\]) and the explicit classical calculation entails $$\label{xp28}
G_{xp_x}^{\rm cl}=\langle xp_x\rangle^{\rm cl}=0,
\qquad G_{xp_y}^{\rm cl}=\langle xp_y\rangle^{\rm cl}=\frac{2}{\beta}\,\frac{{\Omega_{\perp 1}}}{{\Omega_{\perp 2}}^2-{\Omega_{\perp 1}}^2},$$ where the latter is the subtracted classical term in Eq. (\[app19\]) when considering the first two substitutions in (\[xp13\]). In Eq. (\[xp25i\]), the second term is obviously the classical one since $$\label{xp29}
G_{p_zp_z}^{\rm cl}=\langle p_zp_z\rangle^{\rm cl}=\frac{M}{\beta}.$$ The extraction of the classical terms $$\label{xp30}
G_{p_xp_x}^{\rm cl}=\langle p_xp_x\rangle^{\rm cl}=\frac{M}{\beta}\left(1+
\frac{{\Omega_{\perp 1}}^2}{{\Omega_{\perp 2}}^2-{\Omega_{\perp 1}}^2}\right)$$ in the case of the Green function $G_{p_xp_x}^{{\bf p}_0,{\bf x}_0}(\tau,\tau')$ requires the consideration of the last two terms in Eq. (\[xp25g\]). Thus we have shown that the decomposition (\[xp27\]) holds for each of the Green functions (\[xp25a\])–(\[xp25i\]). Note the necessity of subtracting the classical terms since they all diverge in the classical limit of high temperatures ($\beta\to 0$).

[199]{} J.C. Le Guillou and J. Zinn-Justin, Ann. Phys. (N.Y.) [**147**]{}, 57 (1983). J.E. Avron, B.G. Adams, J. Čížek, M. Clay, M.L. Glasser, P. Otto, J. Paldus, and E. Vrscay, Phys. Rev. Lett. [**43**]{}, 691 (1979). B.G. Adams, J.E. Avron, J. Čížek, P. Otto, J. Paldus, R.K. Moats, and H.J. Silverstone, Phys. Rev. A [**21**]{}, 1914 (1980). V.A. Gani, A.E. Kudryavtsev, and V.M. Weinberg, eprint: physics/9708005 (1997). I.D. Feranshuk and L.I. Komarov, J. Phys. A: Math. Gen. [**17**]{}, 3111 (1984). L.D. Landau and E.M. Lifschitz, [*Quantenmechanik*]{}, Sechste Auflage (Akademie-Verlag Berlin, 1979). R. Loudon, Am. J. Phys. [**27**]{}, 649 (1959). L.K. Haines and D.H. Roberts, Am. J. Phys. [**37**]{}, 1145 (1969). R. Cohen, J. Lodenquai, and M. Ruderman, Phys. Rev. Lett. [**25**]{} (1970). M.V. Ivanov and P. Schmelcher, Phys. Rev. A [**61**]{} (in press 2000). J.S. Heyl and L. Hernquist, Phys. Rev. A [**58**]{}, 3567 (1998). H. Ruder, G. Wunner, H. Herold, and F. Geyer, [*Atoms in Strong Magnetic Fields*]{} (Springer-Verlag, Berlin, 1994). J.T. Devreese and F. Brosens in [*Elementary Excitations in Solids*]{}, p. 283, Eds. J.L. Birman, C. Sébenne, and R.F. Wallis, Elsevier, 1992. J.T. Devreese and F. Brosens, Solid State Comm. [**79**]{}, 819 (1991). J.T. Devreese and F. Brosens, Phys. Rev. B [**45**]{}, 6459 (1992). H. Kleinert, [*Path Integrals in Quantum Mechanics, Statistics, and Polymer Physics*]{}, 2nd ed. (World Scientific, Singapore, 1995). R.P. Feynman and H. Kleinert, Phys. Rev. A [**34**]{}, 5080 (1986); see also Chapter 5 of the textbook [@PI]. M. Bachmann, H. Kleinert, and A. Pelster, Phys. Rev. A [**60**]{}, 3429 (1999), eprint: quant-ph/9812063. H. Kleinert, W. Kürzinger, and A. Pelster, J. Phys. A: Math. Gen. [**31**]{}, 8307 (1998), eprint: quant-ph/9806016. H. Kleinert, A. Pelster, and M. Bachmann, Phys. Rev. E [**60**]{}, 2510 (1999), eprint: quant-ph/9902051.
