0$ and $x,y$ in $\mbox{\boldmath $R$}^3$. This follows at once from a similar bound on the heat kernel $\{\exp[-t(\mbox{\boldmath $p$} +\mbox{\boldmath $A$})^2]\}(x,y)$ which, in turn, follows from its representation as a path integral. This was pointed out in \cite{AHS,CSS}. Only the resolvent powers $|\mbox{\boldmath $p$} +\mbox{\boldmath $A$}|^{-s}$ enter the proof of Theorem 1 in \cite{LY}.) Using (\ref{eq:a1}), $H$ is bounded below by $\bar H = \sum_{i=1}^N h_i$ where $h$ is the one-body operator $ h = {\cal T} - \kappa^{-1} \vert \mbox{\boldmath $p$} +\mbox{\boldmath $A$} \vert$. Thus, $E$ is bounded below by $\varepsilon\int B^2+\bar E_N$, where $\bar E_N = q\sum_{j=1}^{[N/q]} \varepsilon_j$ and $\varepsilon_1 \leq \varepsilon_2 \leq...$ are the eigenvalues of $h$. For $e>0$, let $N_{-e}(h)$ be the number of eigenvalues of $h$ less than or equal to $-e$. Choose $\mu >0$ and note that \begin{equation}\label{eq:a3} \bar E_N \geq -N\mu -q\int_{\mu}^\infty N_{-e}(h) {\rm d} e\ . \end{equation} The crucial step in our proof is noting that the positivity of the operator ${\cal T}$ implies that ${\cal T}\geq \mu{\cal T}/e$ when $e\geq \mu$. Thus, ${\cal T}\geq \mu e^{-1}{\cal T} \geq \mu e^{-1}(\mbox{\boldmath $p$} +\mbox{\boldmath $A$})^2-\mu e^{-1} B(x)$ when $e\geq \mu$. By Schwarz's inequality, $\kappa^{-1}|\mbox{\boldmath $p$} +\mbox{\boldmath $A$}|\leq (1/3)e^{-1}\kappa^{-2} (\mbox{\boldmath $p$} +\mbox{\boldmath $A$})^2+(3e/4)$ and hence if we set $\mu=(4/3)\kappa^{-2}$ we obtain $$ h\geq e^{-1}\kappa^{-2}(\mbox{\boldmath $p$} + \mbox{\boldmath $A$})^2-(4/3)e^{-1}\kappa^{-2}B(x) -(3e/4)\equiv h_{e}\ . $$ Thus, $N_{-e}(h)\leq N_{-e}(h_e)$ and this can be estimated by the Cwikel-Lieb-Rozenblum (CLR) bound \cite{L2}, i.e., $N_{-e}((\mbox{\boldmath $p$} +\mbox{\boldmath $A$})^2-U(x))\leq L_3\int[U(x)-e]^{3/2}_+ {\rm d}^3x$ where $[a]_+\equiv\max(a,0)$ and $L_3=0.1156$. In our case: \begin{equation}\label{eq:a4} N_{-e}(h_e)\leq L_3\int\left[\frac{4B(x)}{3} -\frac{e^2\kappa^2}{4} \right]_+^{3/2} {\rm d}^3x\ . \end{equation} Inserting this bound in (\ref{eq:a3}), a simple calculation yields $$ \bar E_N\geq -N\mu-(2\pi/3) q\kappa^{-1} L_3 \int B(x)^2{\rm d}^3 x\ . $$ We choose $\kappa$ so that the field energy terms are non-negative, i.e., $\kappa\geq(16\pi^2/3)L_3\alpha^2 q= 6.1\alpha^2q$. We conclude, by (\ref{eq:a1}), that magnetic stability holds if \begin{equation}\label{eq:a5} q\alpha\leq 0.071\quad\hbox{and}\quad qZ\alpha^2\leq0.052\ . \end{equation} For $q=2$, the first condition is $\alpha\leq 1/28$. For $q=2$ and $\alpha=1/137$ stability occurs if $Z\leq490$. Assuming (\ref{eq:a5}) holds, we then use (\ref{eq:a1}) and choose $\kappa = \min\{0.0315q^{-1}, (\pi Z)^{-1}\}$. Our lower bound on the ground state energy per electron, by this method, is then $-\mu=-(4/3)\kappa^{-2}=-\max\{1345q^2, 13.2Z^2 \}$. {\it Remark:\/} We used the CLR bound in (\ref{eq:a4}). Since the derivation of this bound is not elementary the reader might wish to use an easier to derive bound---at the cost of worsening the final constants. A useful substitute is $$ N_{-e}\leq 0.1054 e^{-1/4}\int\left[U(x)-e/2\right]_+^{7/4}{\rm d}^3 x\ , $$ which is in (2.8) of \cite{LT} and which can be derived by means originally employed for the Lieb-Thirring inequality. This same remark also applies to our other calculations below. \paragraph*{B. The Lieb-Thirring Inequality:} As before we note that $\sum \varepsilon_i= -\int_0^\infty N_{-e}({\cal T}-U){\rm d} e$. We write $\int_0^\infty= \int_0^\mu+\int_\mu^\infty$. The parameter $\mu$ will be optimized below. Noting that ${\cal T}\geq (\mbox{\boldmath $p$} +\mbox{\boldmath $A$})^2-B(x)$ and applying the CLR bound in the same fashion as before to $\int_0^\mu$ yields \begin{equation}\label{eq:b3} L_3\int_0^\mu\int [B(x)+U(x)-e]_+^{3/2}{\rm d}^3x{\rm d}e\ . \end{equation} In $\int_\mu^\infty$ we replace ${\cal T}$ by the lower bound $\mu e^{-1}[(\mbox{\boldmath $p$} + \mbox{\boldmath $A$})^2-B(x)]$ and obtain $N_{-e}({\cal T}-U)\leq N_{-e}(\mu e^{-1}[(\mbox{\boldmath $p$} + \mbox{\boldmath $A$})^2-B]-U)$. A further application of the CLR inequality yields the bound on $\int_{\mu}^\infty$ \begin{equation}\label{eq:b4} L_3\int_{\mu}^\infty\int[B(x)+(e/\mu)U(x)- e^2/\mu]_+^{3/2}{\rm d} x{\rm d} e\ . \end{equation} It is easy to see that for any $0<\gamma<1$ the integrand in (\ref{eq:b3}) is bounded above by $$ \sqrt{2}\left([B(x)-\gamma e^2/\mu]_+^{3/2}+ [U(x)-(1-\gamma)e]_+^{3/2}\right)\ . $$ Treating the integrand in (\ref{eq:b4}) in a similar fashion and combining the inequalities we find \begin{eqnarray*} \sum|\varepsilon_i|&&\leq\sqrt{2}L_3 \int\biggl\{\int_0^\infty [B(x)- \gamma e^2/\mu]_+^{3/2}{\rm d} e \nonumber\\&& +\int_0^\mu [U(x)-(1-\gamma) e]_+^{3/2}{\rm d} e \\ && +\int_\mu^\infty\left[(e/\mu) U(x) - (1-\gamma) e^2/\mu\right]_+^{3/2}{\rm d} e \biggr\}{\rm d}^3x\ .\nonumber \end{eqnarray*} After extending the last two integrals to $\int_0^\infty$, a straightforward computation yields \begin{eqnarray*} \sum|\varepsilon_i|\leq&&\sqrt{2}L_3\int \biggl\{\frac{2}{5(1-\gamma)} U(x)^{5/2} +\frac{3\pi\mu^{1/2}}{16\gamma^{1/2}} B(x)^2\\ && \mbox{} + \frac{3\pi}{128} \mu^{-3/2}(1-\gamma)^{-5/2} U(x)^4\biggr\}{\rm d}^3x\ . \end{eqnarray*} Optimizing over $\mu$ yields (\ref{eq:lt2}) \paragraph*{C. Proof of Theorem 1:} We turn now to our third illustration of the concept of running energy scale and prove the stability directly, not relating it to the relativistic problem. By this method we get the correct dependence of the ground state energy on Z and also somewhat better critical constants than in (\ref{eq:a5}). Following \cite{LY} we first replace the Coulomb potential by a single particle potential in (\ref{eq:c1}) below. We break up $\mbox{\boldmath $R$}^3$ into Voronoi cells defined by the nuclear locations, i.e., $\Gamma_j = \{ x: |x-R_j| \le |x-R_k|\ {\rm for\ all}\ k \}$ is the j--th Voronoi cell. Each $\Gamma_j $ contains a ball centered at $R_j$ with radius $D_j=$ min$\{|R_j-R_k|:j \ne k\}/2$. The following bound on $V_{\rm c}$ is proved in \cite{LY}: Choose some $0<\lambda<1$. Then \begin{equation}\label{eq:c1} V_{\rm c}\ge - \sum_{i=1}^N W(x_i)+ \sum_{j=1}^K \frac{Z^2}{8D_j} \ , \end{equation} where $W(x)=Z|x-R_j|^{-1}+F_j(x)$ for $x \in \Gamma_j$ with $F_j(x)$ defined by \begin{eqnarray*} (2D_j)^{-1}(1-D_j^{-2}|x-R_j|^2)^{-1}&\ \ \mbox{for }\ & |x-R_j| \le \lambda D_j\\ (\sqrt {2Z}+1/2)|x-R_j|^{-1} &\ \ \mbox{for }\ & |x-R_j| > \lambda D_j\ . \end{eqnarray*} The point about this inequality is that the potential $W$ has the same singularity near each nucleus as $V_c$ has, and that the rightmost term in (\ref{eq:c1}) is repulsive. This term will be responsible for stabilizing the system. The problem is thus reduced to obtaining a lower bound on $q\sum^\prime \varepsilon_j$, where $\sum^\prime \varepsilon_j$ is the sum of the first $[N/q]$ negative eigenvalues of of ${\cal T} - W$. Note that Theorem~2 cannot be applied directly to this problem, since $W$ is neither integrable to the power $5/2$ nor to the power 4. Instead we have to do the calculations directly. For $\nu >0$ (a number that is chosen later) set $W_{\nu}(x)\equiv (W(x)-\nu)_+$ and note that $W(x)-\nu\leq W_\nu(x)$. Then, as in (\ref{eq:a3}), $ q\sum^\prime \varepsilon_j \ge-N\nu -q \int_0^\infty N_{-e}({\cal T}- W_{\nu}) {\rm d} e $. Again, \begin{eqnarray} \int_0^\infty N_{-e}({\cal T}- W_{\nu}) {\rm d} e &\le& \int_0^{\mu} N_{-e}({\cal T}- W_{\nu}) {\rm d} e \nonumber\\ &+&\int_{\mu}^\infty N_{-e}(\mu e^{-1}{\cal T}- W) {\rm d} e\ , \label{eq:c2} \end{eqnarray} where we have replaced $W_{\nu}(x)$ by $W(x)$ in the second term. Applying the CLR bound to the first expression on the right side we obtain $ L_3 \int \int_0^{\mu} [B(x)+ W_{\nu}(x)-e]^{3/2}_+{\rm d} e {\rm d} ^3x $ which can be bounded, as in part B, by \begin{eqnarray} \sqrt 2 L_3 \int \biggl\{ && \int_0^\mu[B(x)- \frac{\gamma e^2}{\mu}]^{3/2}_+{\rm d} e \nonumber\\ && +\mbox{} \frac{2}{5}(1-\gamma)^{-1}W_{\nu}(x)^{5/2} \biggl\}{\rm d} ^3x \ , \label{eq:c3} \end{eqnarray} for any $0< \gamma <1$. The difficulty in dominating the second term in (\ref{eq:c2}) comes from the Coulomb singularity of $W(x)$ which is not fourth power integrable. The singularity can be controlled using the following operator inequality which follows from the diamagnetic inequality $\int |(\mbox{\boldmath $p$}+\mbox{\boldmath $A$})\psi|^2 {\rm d} ^3x \ge \int \bigl|\mbox{\boldmath $p$}|\psi|\bigr|^2 {\rm d} ^3x$ and Lemma~2a on p.~708 of \cite{LD}. $$ (\mbox{\boldmath $p$} +\mbox{\boldmath $A$})^2-Z/|x|\geq -\cases{Z^2/4 +(3/2)ZR^{-1},& if $|x|\leq R$\cr Z|x|^{-1},& if $|x|\geq R$}\ . $$ Choose $R=\lambda D_j$ and write $(\mbox{\boldmath $p$} + \mbox{\boldmath $A$})^2 =\beta(\mbox{\boldmath $p$} +\mbox{\boldmath $A$})^2+ (1-\beta)(\mbox{\boldmath $p$} +\mbox{\boldmath $A$})^2$ for some $0<\beta<1$. Then, by scaling, $$ (\mu/e){\cal T} -W_{\nu}\ge (\mu/e)(1-\beta) (\mbox{\boldmath $p$}+\mbox{\boldmath $A$})^2-(\mu/e)B -\widetilde W\ , $$ where $\widetilde W(x,e)= \widetilde G_j(x,e)+F_j(x)$ for $x \in \Gamma_j$ with $\widetilde G_j(x,e)$ defined by \begin{eqnarray*} (Z^2e/4\beta\mu) + 3Z/(2\lambda D_j) &\ \ \mbox{for }\ & |x-R_j| \le \lambda D_j\\ Z|x-R_j|^{-1} &\ \ \mbox{for }\ & |x-R_j| > \lambda D_j\ . \end{eqnarray*} Note that $\widetilde W$ depends on $e$. Again, as in part B, we can use the CLR bound on the second term in (\ref{eq:c2}) to obtain [when $1-\gamma\geq Z^2/(4\beta\mu)$] \begin{eqnarray} \sqrt 2&&L_3(1-\beta)^{-3/2} \int\biggl\{\int_{\mu}^{\infty}[B(x)-\gamma e^2/\mu]_+^{3/2} {\rm d} e \nonumber\\ &&+ {\mu}^{-3/2} \int_0^\infty[e \widetilde W(x,e) -(1-\gamma)e^2]_+^{3/2}{\rm d} e \biggr\}{\rm d} ^3x \ . \label{eq:c4} \end{eqnarray} First we compute the last integral in (\ref{eq:c4}), which is $$ \sum_{j=1}^K\int_{\Gamma_j}\int_0^\infty [e \widetilde G_j(x,e)+e F_j(x)-(1-\gamma)e^2]_+^{3/2}{\rm d} e {\rm d}^3x\ . $$ Now split the $\Gamma_j$ integral into an inner integral $|x-R_j| \le \lambda D_j$ and an outer integral $|x-R_j| > \lambda D_j$. The inner integral yields, using the definitions of $\widetilde G_j$ and $F_j$, \begin{eqnarray} \frac{3\pi^2}{32}(1-\gamma-\frac{Z^2}{4\beta \mu})^{-5/2} \int_0^{\lambda}\biggl[&&\frac{1}{2(1-r^2)}\nonumber\\ && \mbox{}+\frac{3Z}{2\lambda}\biggr]^4r^2{\rm d} r D_j^{-1}\ . \label{eq:c6} \end{eqnarray} To bound the outer integral from above we replace $\Gamma_j$ by $\mbox{\boldmath $R$}^3$ and get \begin{equation} (3 \pi^2/32)(1-\gamma)^{-5/2}(\sqrt Z +\sqrt{1/2})^8 (\lambda D_j)^{-1}\ . \label{eq:c7} \end{equation} Combining (\ref{eq:c3})--(\ref{eq:c7}) we find that the sum of the negative eigenvalues of ${\cal T} - W_{\nu}$ is bounded below by \begin{equation} -a\int W_{\nu}(x)^{5/2}{\rm d}^3x -b\int B(x)^2 {\rm d}^3x -c\sum_{j=1}^K D_j^{-1} \ . \label{eq:c8} \end{equation} Here $a=q(2 \sqrt 2/5)L_3(1-\gamma)^{-1}$, \begin{eqnarray*} b&=&q\frac{3\pi \sqrt 2}{16}L_3(1-\beta)^{-3/2} (\mu/\gamma)^{1/2}\\ c&=&q\frac{3\pi^2 \sqrt 2}{32}L_3(1-\beta)^{-3/2} \mu^{-3/2} \Bigl\{ \frac{(\sqrt Z +\sqrt{q/4})^8}{\lambda(1-\gamma)^{5/2}} \\ &&+ (1-\gamma-\frac{Z^2}{4\beta \mu})^{-5/2} \int_0^{\lambda} \Bigl[\frac{q}{4(1-r^2)}+\frac{3Z}{2\lambda}\Bigr]^4 r^2{\rm d} r \Bigr\}. \end{eqnarray*} To simplify the stability condition we have artificially increased the bounds by recalling that $q\geq 2$ and twice replacing $1/2$ by $q/4$ in the definition of $c$. We choose $\beta=1/8$, $\gamma=1/2$, $\lambda=8/9$ and $\mu$ so that $b=(8\pi\alpha^2)^{-1}$. The {\it stability condition} $c\leq Z^2/8$ [see (\ref{eq:c1})] now depends {\it only} on the 2 parameters $X=qZ\alpha^2$ and $Y=q\alpha$. A straightforward, but lengthy calculation shows that the stability condition holds if $X=X_0\equiv0.082$ and $Y=Y_0\equiv0.12$. The condition is monotone in $Y$, so it holds for $X=X_0$, $Y\leq Y_0$. Although our condition does not hold for {\it all} $X\leq X_0$, $Y\leq Y_0$ we can use the $Z$-monotonicity of $E$ to conclude stability in this range; this proves (\ref{eq:4}). With the same values of $\beta,\gamma$ and $\lambda$ and with $q=2$ the values $Z=1050$, $\alpha=1/137$ also give stability. To derive (\ref{eq:3}), note that $W(x)\leq Q|x-R_j|^{-1}$ for $x\in\Gamma_j$. Using this bound and replacing $\Gamma_j$ by $\mbox{\boldmath $R$}^3$, one easily obtains $-\sqrt{2}\pi^2L_3qKQ^3\nu^{-1/2}-N\nu$ as a lower bound on the $-a\int W_\nu^{5/2}$ term in (\ref{eq:c8}). Optimizing over $\nu$ yields (\ref{eq:3}) when $X=X_0$, $Y\leq Y_0$. In this case, $Z\geq Z_0\equiv 5.7 q$. If $X\leq X_0$, $Y\leq Y_0$ and $Z\geq Z_0$ we get a lower bound on $E$ by increasing $\alpha$ until $X=X_0$, $Y\leq Y_0$; this yields (\ref{eq:3}) with $Q=Q(Z)$. Otherwise, with $Z