text:
# Quenched Spin Tunneling and Diabolical Points in Magnetic Molecules: II. Asymmetric Configurations

## I Introduction

### A The story so far

In a previous paper, hereafter cited as I, we studied the tunneling of a spin governed by the Hamiltonian

$$=k_2J_z^2+(k_1k_2)J_x^2g\mu _B𝐉𝐇,$$
(1)
where $`𝐉`$ is a dimensionless spin operator, and $`k_1>k_2>0`$. This Hamiltonian is the simplest descriptor of the magnetic properties of the molecule \[(tacn)6Fe8O2(OH)12\]8+ (or just Fe8 for short), with $`J=10`$, $`k_10.33`$ K, and $`k_20.22`$ K . Interest in this molecule arises because of its rich low temperature magnetic behavior, which include hysteresis at the level of one molecule , and more recently, the discovery  of an entire lattice of diabolical points  in its magnetic spectrum, a subset of which was predicted to exist earlier . It is the latter property that we wish to continue investigating in this paper.

In paper I, attention was confined to the case where the external magnetic field $`𝐇`$ is along the hard axis $`𝐱`$. In this case, the classical energy, which may be viewed as the expectation value $``$ of the Hamiltonian in a spin-coherent state, is symmetric about the $`xy`$ plane, and the problem is analogous to that of a massive particle in one dimension in a reflection symmetric double well potential. It is then natural to consider the tunneling between the symmetrically related states localized in the left- and right-hand wells. In the magnetic case, the analogous states are those with predominatly positive and negative values of $`J_z`$, at least as long as $`H_x`$ is not so large as to bring the classical minima very close to the $`x`$ axis. The surprise is that the tunnel splitting between the ground states oscillates as a function of $`H_x`$, vanishing exactly at a series of points. In fact, the splitting between higher pairs of levels also vanishes at just these points, as noted in I and earlier studies .

In addition to observing the quenching of ground state tunneling when $`𝐇𝐱`$, however, Wernsdorfer and Sessoli  also performed experiments with $`H_z0`$. The reflection symmetry of the classical energy is now lost, but if the value of $`H_z`$ is chosen properly, it is possible to bring an excited state in the positive $`J_z`$ well into approximate degeneracy with the ground state of the negative $`J_z`$ well. The new discovery by them is that if $`H_x`$ is now varied, the tunnel splitting between the degenerate or quasi-degenerate levels again oscillates. It is theoretically understood that if both $`H_z`$ and $`H_x`$ are properly tuned, the splitting vanishes exactly in this case too . Experimentally, of course, one can never see a perfect zero in the splitting, and Wernsdorfer and Sessoli only see a minimum in the Landau-Zener-Stückelberg transition rate between the levels in question. The minima are so deep, however, that there is little doubt that the underlying tunneling matrix element is quenched.

When $`𝐇𝐱`$, the Hamiltonian (1) is invariant under a 180 rotation about $`𝐱`$, and the quenchings can be understood from the viewpoint of the von Neumann-Wigner theorem as allowed crossings of energy levels with different parities under this rotation. A similar argument can be made when $`𝐇𝐳`$. When $`𝐇`$ has both $`x`$ and $`z`$ components, however, the Hamiltonian has no obvious symmetry, and the above theorem states that an intersection of two energy levels is infinitely unlikely as a single parameter in the Hamiltonian is varied. For a real symmetric Hamiltonian, it is known  that one must vary two parameters to obtain an intersection. Since we can choose the matrices of both $`J_z`$ and $`J_x`$ to be real in the $`J_z`$ basis, these two parameters can be taken as $`H_x`$ and $`H_z`$. The isolated points  in the $`H_x`$-$`H_z`$ plane where any two levels intersect are precisely what Herzberg and Longuet-Higgins  call conical intersections and what Berry and Wilkinson  call diabolical points. The latter terminology originates in the resemblance of the energy surface—a double cone with a common vertex —to an Italian toy called the diavolo.

### B Content and plan of this paper; the perfect lattice hypothesis

In this paper we shall allow $`𝐇`$ to lie in the $`xz`$ plane, with a view to studying the tunneling in the asymmetrical well, and locating the diabolical points. As in paper I, our analysis is based on the discrete phase integral (DPI), or Wentzel-Kramers-Brillouin method . This method is semiclassical in character, with $`1/J`$ playing the same role as $`\mathrm{}`$ in the continuum phase integral method. To introduce this method, let us write the Schrödinger equation $`|\psi =E|\psi `$ in the $`J_z`$ eigenbasis $`\{|m\}`$. With $`J_z|m=m|m`$, $`m|\psi =C_m`$, $`m||m=w_m`$, and $`m||m^{}=t_{m,m^{}}`$ ($`mm^{}`$), we have

$$\stackrel{}{}_{n=m2}^{m+2}t_{m,n}C_n+w_mC_m=EC_m,$$
(2)
where the prime on the sum indicates omission of the $`n=m`$ term. A vivid picture of the approximation can be obtained if we think of Eq. (2) as the tight-binding model for an electron in a one-dimensional lattice with sites labelled by $`m`$, and on-site ($`w_m`$) and hopping ($`t_{m,m\pm 1}`$, $`t_{m,m\pm 2}`$) energies. If $`J1`$, these matrix elements vary slowly with $`m`$, on a length scale of order $`J`$ in fact. We may then use the approximation of semiclassical dynamics by working entirely in terms of wavepackets whose spatial extent is much less than the length scale over which the properties of the lattice vary, i.e., $`J`$, and whose spread in Bloch vectors is much less than the width of the Brillouin zone, i.e., $`2\pi `$. These ideas have close counterparts in the continuum quasiclassical method, and the DPI method is nothing but the discrete analog.

When $`𝐇𝐱`$, the problem can also be approached using instantons—indeed the oscillations in the splitting were discovered in this way. When $`𝐇`$ has other components besides $`H_x`$, however, the DPI method is, to our knowledge, the only successful one to date. Villain and Fort’s approach  is also an approximate version of this method that makes artful use of some special features of the Fe8 problem and works for small values of the field. Our analysis is more prosaic, and almost self-evident once one has understood how to deal with the new feature in Eq. (2)—the presence of second neighbour hopping. This gives rise to novel turning points with no continuum analogues. Connection formulas for these turning points are given in Ref. . We have quoted the results of our analysis before (b), but the details are only being presented here.

Our main result for the specific Hamiltonian (1) is for the locations of the diabolical points. We find that the $`\mathrm{}^{}`$th level in the negative $`J_z`$ well (where $`\mathrm{}^{}=0`$ denotes the lowest level) and the $`\mathrm{}^{\prime \prime }`$th level in the positive $`J_z`$ well are degenerate when (see Fig. 1)

$`{\displaystyle \frac{H_z(\mathrm{}^{},\mathrm{}^{\prime \prime })}{H_c}}`$ $`=`$ $`{\displaystyle \frac{\sqrt{\lambda }(\mathrm{}^{\prime \prime }\mathrm{}^{})}{2J}}`$ (3)
$`{\displaystyle \frac{H_x(\mathrm{}^{},\mathrm{}^{\prime \prime })}{H_c}}`$ $`=`$ $`{\displaystyle \frac{\sqrt{1\lambda }}{J}}\left[Jn\frac{1}{2}(\mathrm{}^{}+\mathrm{}^{\prime \prime }+1)\right],`$ (4)
with $`n=0,1,\mathrm{},2J(\mathrm{}^{}+\mathrm{}^{\prime \prime }+1)`$. Here, $`\lambda =k_2/k_1`$, and $`H_c=2k_1J/g\mu _B`$.

It should be stressed that Eqs. (3) and (4) only represent the first terms of a series in $`1/J`$, and that our DPI calculations give no reason to believe that the higher order terms are absent. Yet a large amount of empirical evidence suggests just this, i.e., that the formulas are exact as written! That the diabolical points lie on a perfect centered rectangular lattice, and that many pairs of levels are simultaneously degenerate, we refer to as the perfect lattice hypothesis. It has previously been made for the ground state on the $`H_z=0`$ line by Villain and Fort , and extended to include simultaneous degeneracy of the higher states by us (b). This was shown by perturbation theory in $`\lambda `$ to $`O(\lambda ^3)`$ for all $`J`$, and analytically for $`J2`$. We have since found analytically that for $`J=3/2`$, the additional diabolical point at $`(H_x,H_z)=(\sqrt{1\lambda },\sqrt{\lambda })H_c/3`$ is exact. Of course, Kramers theorem guarantees double degeneracy of all energy eigenvalues at the point $`H_x=H_z=0`$ for all half-integral $`J`$. We have also checked this hypothesis numerically for a variety of values of $`\lambda `$, and $`J`$ up to 10. We present our results for $`J=5`$ in Table I. The deviations in the locations of the diabolical points from the perfect lattice hypothesis predictions are never more than $`10^{10}`$, certainly much less than order $`1/J`$. In fact, the values listed are below the numerical tolerance that we prescribed.

Further support for the perfect lattice hypothesis comes from the fact that it is consistent with the following duality property of the Hamiltonian (1). If we set $`k_1=1`$, $``$ may be written as $`(\lambda ,H_x,H_z)`$, showing its dependence on the three parameters, $`\lambda `$, $`H_x`$ and $`H_z`$. By a $`90^{}`$ rotation about the axis $`(\widehat{𝐱}+\widehat{𝐳})/\sqrt{\mathrm{𝟐}}`$, we obtain the transformation

$$(\lambda ,H_x,H_z)(1\lambda ,H_z,H_x).$$
(5)
In particular, the spectra of the two Hamiltonians are so related, and ranking the levels is order of increasing energy, we see that if the levels with ordinal numbers $`k`$ and $`k+1`$ are degenerate when $`H_x=f_x(\lambda )`$ and $`H_y=f_y(\lambda )`$, then level numbers $`2J+2k`$ and $`2J+1k`$ are degenerate when $`H_x=f_y(1\lambda )`$ and $`H_y=f_x(1\lambda )`$. These conditions do not constrain the functions $`f_x`$ and $`f_y`$ in any real way, however, so the discovery of formulas (3) and (4) must presently be put down to serendipity.

The plan of our paper is as follows. In Sec. II, we analyze the asymmetric double well problem in completely general terms. We briefly review the DPI method and the asymmetric double-well problem in one dimensional quantum mechanics. In contrast to the situation that prevails in that case, we must deal with four independent DPI wavefunctions at every point. Consequently the matching problem is harder and more subtle. Its study forms the bulk (subsections C to G) of the analysis. The main result of this analysis is a formula for the diabolo \[Eq. (73)\], and a two level system Hamiltonian \[Eq. (74)\]. The tunneling amplitude appearing in these formulas is found in a general form \[Eq. (70)\] involving a small number of action integrals running between various turning points. Readers who are not interested in the details of the analysis should skip to subsection H and read on from Eq. (73). They may also find the last paragraph of Sec. II interesting.

In Sec. III, we apply the general analysis to the model (1) for Fe8. To keep the analysis tractable, we will assume that $`H_z/H_c`$ is small, although there is no reason why one could not apply the results of Sec. II numerically for arbitrary values of $`H_z`$. It may in fact be interesting to do this to further investigate the perfect lattice hypothesis. Further, as discussed in Sec. III, the duality property of $``$ immediately extends our work to large $`H_z`$ and small $`H_x`$.

## II General Formula for Diabolo in terms of Action Integrals

### A Summary of DPI method, critical curves, etc.

Very briefly, the DPI analysis proceeds as follows. (See Ref.  and paper I for details.) Let the energy of an electron wavepacket in the equivalent tight-binding model be given by $`_{\mathrm{sc}}(q,m)`$ where $`q`$ and $`m`$ are the mean wavevector and position of the wavepacket, respectively. Holding $`m`$ fixed, we may think of $`_{\mathrm{sc}}(q,m)`$ as a dispersion relation for the local energy band. The mean velocity of the wavepacket is then given by $`v(q,m)=_{\mathrm{sc}}(q,m)/q`$. If we regard $`m`$ as a continuous variable, and approximate $`w_m`$, $`t_{m,m\pm 1}`$, and $`t_{m,m\pm 2}`$ by smooth functions $`w(m)`$, $`t_1(m)`$, and $`t_2(m)`$, then

$`_{\mathrm{sc}}(q,m)`$ $`=`$ $`w(m)+2t_1(m)\mathrm{cos}q+2t_2(m)\mathrm{cos}(2q),`$ (6)
$`v(q,m)`$ $`=`$ $`2\mathrm{sin}q(m)\left(t_1(m)+4t_2(m)\mathrm{cos}q(m)\right).`$ (7)
In close analogy with the continuum method, one can show that the quasiclassical solution to the wavefunction for a given energy $`E`$ is given by linear combinations of solutions of the form

$$C_m\frac{1}{\sqrt{v(m)}}\mathrm{exp}\left(i^mq(m^{})𝑑m^{}\right),$$
(8)
where $`q(m)`$ and $`v(m)`$ are given by

$$E=_{\mathrm{sc}}(q(m),m),v(m)=v(q(m),m).$$
(9)
It follows from Eq. (8) that the quasiclassical solution breaks down when $`v(m)=0`$—which we call a turning point, as in the continuum case. At these points, the solutions must be augmented by connection formulas. In addition to the values $`q=0`$ and $`q=\pi `$, $`v(m)`$ also vanishes at $`q=q^{}(m)=\mathrm{cos}^1\left(t_1(m)/4t_2(m)\right)`$. We will get turning points whenever $`E`$ equals $`U_0(m)`$, $`U_\pi (m)`$, or $`U_{}(m)`$, where these three functions are $`_{\mathrm{sc}}(0,m)`$, $`_{\mathrm{sc}}(\pi ,m)`$, and $`_{\mathrm{sc}}(q^{},m)`$ respectively. These three curves, which we call critical curves, collectively play the same role as the potential energy in the continuum case.

For Fe8, in the same reduced variables as in I \[$`\mu =m/\overline{J}`$, energies in units of $`k_1\overline{J}^2`$, $`\overline{J}=(J+\frac{1}{2})`$\], the on-site and hopping energies are given by

$`w(m)`$ $`=`$ $`\frac{1}{2}(1+\lambda )(1\mu ^2)2h_z\mu ,`$ (10)
$`t_1(m)`$ $`=`$ $`h_x(1\mu ^2)^{1/2},`$ (11)
$`t_2(m)`$ $`=`$ $`\frac{1}{4}(1\lambda )(1\mu ^2).`$ (12)
Here, $`\lambda =k_2/k_1`$, $`h_x=JH_x/\overline{J}H_c`$, $`h_z=JH_z/\overline{J}H_c`$, and $`H_c=2k_1J/g\mu _B`$. Thus the critical curves are given by

$`U_0(\mu )`$ $`=`$ $`1\mu ^22h_x(1\mu ^2)^{1/2}2h_z\mu ,`$ (13)
$`U_\pi (\mu )`$ $`=`$ $`1\mu ^2+2h_x(1\mu ^2)^{1/2}2h_z\mu ,`$ (14)
$`U_{}(\mu )`$ $`=`$ $`\lambda (1\mu ^2){\displaystyle \frac{h_x^2}{1\lambda }}2h_z\mu .`$ (15)
These curves are shown in Fig. 2. The energy $`U_\pi (m)`$ is the upper edge of the band $`_{\mathrm{sc}}(q,m)`$ for all $`m`$. The lower band edge is given by $`U_0(m)`$ for values of $`m`$ close to $`\pm J`$, and by $`U_{}(m)`$ for $`m`$ in the central region. In the central region, the second neighbour hopping element $`t_2(m)`$ is sufficiently large that the local energy band $`_{\mathrm{sc}}(q,m)`$ has its global minimum not at $`q=0`$, but at $`q^{}`$. Since the energy $`U_0`$ lies in the band, $`U_\pi >U_0>U_{}`$. By contrast, in the outer $`m`$ regions, $`_{\mathrm{sc}}(q,m)`$ has only one minimum (at $`q=0`$) and only one maximum (at $`q=\pi `$) for real $`q`$. The energy $`U_{}`$ lies outside the band, so that once again $`U_\pi >U_0>U_{}`$. The curves $`U_0(m)`$ and $`U_{}(m)`$ touch one another with a common tangent at $`m=\pm m^{}`$, where

$$m^{}=\overline{J}\left[1\frac{h_x^2}{(1\lambda )^2}\right]^{1/2}.$$
(16)
The turning point where $`E=U_{}`$ is special if it happens at a value of $`m`$ where $`U_{}`$ lies below the lower band edge. The wavevector $`q^{}`$ at the turning point is then complex, and the wavefunction $`C_m`$ changes from a decaying (or growing) exponential on one side to a decaying (or growing) exponential with an oscillatory envelope on the other side. These turning points are the new feature caused by second neightbor hopping that we referred to in Sec. I.

### B Nature of asymmetric double well wavefunctions

To understand how the DPI solutions are to be used to find the eigenstates, it is useful to discuss the corresponding problem for an asymmetric double well in the continuum case. Suppose the potential is as drawn in Fig. 3. The potential minima are at $`x_{0\pm }`$, and for the energy $`E`$ drawn, the turning points are at $`x_a^{}`$ and $`x_b^{}`$ in the left well, and $`x_b^{\prime \prime }`$ and $`x_a^{\prime \prime }`$ in the right well. We use the quasiclassical approximation to find a wavefunction $`\psi ^{}(x)`$ on the left hand side as follows. First, we choose the solution which decays exponentially as $`xx_a^{}\mathrm{}`$. This solution is matched via connection formulas at $`x_a^{}`$ to an oscillatory solution in the region $`x_a^{}<x<x_b^{}`$. We then use connection formulas at $`x_b^{}`$ to find the quasiclassical solution in the region $`x>x_b^{}`$. We repeat this procedure on the right hand side to find a wavefunction $`\psi ^{\prime \prime }(x)`$ that decays to zero as $`xx_a^{\prime \prime }\mathrm{}`$. The last step is to demand that the wavefunctions $`\psi ^{}(x)`$ and $`\psi ^{\prime \prime }(x)`$ be the same in the central region, i.e., in the vicinity of $`x=0`$. This demand will be unsatisfiable for an arbitrarily chosen energy $`E`$, and will provide one with the eignevalue condition.

There are two remarks that we wish make about the above procedure. The first remark concerns the basic nature of the solution. In general, in the central region, the left solution $`\psi ^{}(x)`$ will be a linear combination of two parts, $`\psi _d^{}(x)`$ and $`\psi _g^{}(x)`$, that are exponentially decaying and growing as $`xx_b^{}`$ increases, respectively. Likewise the right solution, $`\psi ^{\prime \prime }`$, will be a sum of parts that decay ($`\psi _d^{\prime \prime }`$) and grow ($`\psi _g^{\prime \prime }`$) as $`x_b^{\prime \prime }x`$ increases. The key point is that the growing parts $`\psi _g^{}`$ and $`\psi _g^{\prime \prime }`$ must be present in an eigenstate, for without them, there is no way that the values and slopes of $`\psi ^{}`$ and $`\psi ^{\prime \prime }`$ could be made to agree at $`x=0`$, say.

The second remark is technical. If the potential well is reasonably parabolic near the minimum at $`x_0`$, the Schrödinger equation can be solved directly for any choice of $`E`$ in terms of parabolic cylinder functions, and we can always find a linear combination that will decay to zero as $`xx_a^{}\mathrm{}`$. This linear combination will have both growing and decaying pieces as $`xx_b^{}`$ grows. In this way we can obtain the wavefucnction $`\psi ^{}(x)`$ on the entire left hand side without using connection formulas at $`x_a^{}`$ or $`x_b^{}`$. Once one has found $`\psi ^{}(x)`$ sufficiently far to the right of $`x_b^{}`$ in this way, one can write the parabolic cylinder functions in quasiclassical form which can then be extended in this form all the way to $`x=0`$. The right-hand wavefunction can be treated in the same way. This device leads to considerable savings in labor.

We now apply these ideas to our problem. In what follows, we will denote quantities pertaining to the left hand solution or the left hand side of the well by either a single prime or a suffix $``$, while analogous right-hand quantities will carry a double prime or a $`+`$ suffix. We consider an energy $`E`$ as drawn in Fig. 4, which leads to turning points at $`m_a^{}`$, $`m_b^{}`$, $`m_c^{}`$ on the left hand side, and $`m_c^{\prime \prime }`$, $`m_b^{\prime \prime }`$, and $`m_a^{\prime \prime }`$ on the right. The wavefunction $`C_m^{}`$ will decay away from the well bottom as $`mm_a^{}`$ decreases, oscillate in the classically allowed region $`m_a^{}<m<m_b^{}`$ . In the region just to the right of $`m_b^{}`$ it will consist of a decaying part and a growing part. The new feature will be encountered at $`m_c^{}`$ where $`E=U_{}`$. For $`m>m_c^{}`$, both the growing and decaying parts will acquire oscillatory envelopes. Similar remarks apply to the right side wavefunction $`C_m^{\prime \prime }`$.

### C DPI wavefunction in the leftmost forbidden region

We are now ready to find the wavefunction explicitly. Let us start constructing $`C_m^{}`$ from the left. The Hamilton-Jacobi equation in Eq. (9) has the general solution

$`\mathrm{cos}q(m)`$ $`=`$ $`{\displaystyle \frac{t_1(m)\pm [t_1^2(m)4t_2(m)f(m)]^{1/2}}{4t_2(m)}};`$ (17)
$`f(m)`$ $`=`$ $`w(m)2t_2(m)E.`$ (18)
This leads to four values of $`q(m)`$ for any $`E`$, since if $`q`$ is a solution, so is $`q`$. In the region $`mm_a^{}`$, since $`E<U_0`$, all four solutions are pure imaginary. We write the two which lead to decaying wavefunctions as $`\mu \mu _a^{}`$ becomes large and negative as

$$q_{1,2}=i\kappa _{1,2}(\mu ),$$
(19)
where $`\kappa _2>\kappa _1>0`$, and the corresponding DPI solutions as

$`C_{m,1}^{}`$ $`=`$ $`A^{}|v_1(m)|^{1/2}\mathrm{exp}\left(i{\displaystyle ^m}q_1(m^{})𝑑m^{}\right),`$ (20)
$`C_{m,2}^{}`$ $`=`$ $`B^{}|v_2(m)|^{1/2}\mathrm{exp}\left(i{\displaystyle _{m_c^{}}^m}q_2(m^{})𝑑m^{}\right),`$ (21)
with $`v_i(m)=v(q_i(m),m)`$. We take $`A^{}`$ and $`B^{}`$ to be real without any loss of generality. Note that we have left the lower limit of the phase integral for $`C_{m,1}^{}`$ unspecified and written it as $`m_c^{}`$ for $`C_{m,2}^{}`$. The reasons for this will become clear shortly.

To see how these two solutions behave as $`m`$ approaches $`m_a^{}`$, and continues beyond this point, let us note that

$$\mathrm{cosh}\kappa _{1,2}=\frac{|t_1|[t_1^24t_2f]^{1/2}}{4t_2}.$$
(22)
As $`mm_a^{}`$, $`\mathrm{cosh}\kappa _11`$, i.e., $`\kappa _10`$, while $`\mathrm{cosh}\kappa _21+2|t_1|/4t_2>1`$. As we cross the point $`m_a^{}`$, $`q_1`$ will become real, while $`q_2`$ will continue to be pure imaginary and large. Thus the solution $`C_{m,2}^{}`$ continues to hold at $`m_a^{}`$, while $`C_{m,1}^{}`$ breaks down at $`m_a^{}`$ , and must be related to a the solution for $`m>m_a^{}`$ by a connection formula. It is clear that the wavevector(s) for $`C_{m,2}^{}`$ will continue to be given by Eq. (22) as $`mm_b^{}`$, while those for $`C_{m,1}^{}`$ will again approach zero, necessitating the use of connection formulas to go on to $`m>m_b^{}`$. It is here that the technical remark about sidestepping the use of connection formulas that was made in connection with the continuum antisymmetric double-well is relevant. The solution $`C_{m,1}^{}`$ can be approximated by a harmonic oscillator wavefunction provided the energy $`E`$ is not very far from the minima of $`U_0`$. The asymptotic forms of this wavefunction give us the quasiclassical wavefunction in the regions $`m<m_a^{}`$ and $`m_b^{}`$ more simply. We therefore turn to this subproblem.

### D Jumping across the potential well

The assumption that $`EU_0(\mu _{0\pm })`$ is not very large, means that $`q_1(m)`$ is never very far from zero, and we may expand $`_{\mathrm{sc}}(q,m)`$ in powers of $`q`$ and $`m+m_0^{}`$. As in I, we write

$$_{\mathrm{sc}}(q,m)=E_{}+\frac{1}{2M_{}}q^2+\frac{1}{2}M_{}\omega _{}^2(mm_{0,})^2+\mathrm{}.$$
(23)
Since $`q`$ and $`m`$ are conjugate variables, we can write $`C_{m,1}^{}`$ as the solution to the differential equation

$$_{\mathrm{sc}}(i_m,m)C_{m,1}^{}=EC_{m,1}^{}.$$
(24)
Introducing two new variables $`z`$ and $`\nu ^{}`$ by the equations

$`m`$ $`=`$ $`m_0+(2M_{}\omega _0)^{1/2}z,`$ (25)
$`E`$ $`=`$ $`E_{}+(\nu ^{}+\frac{1}{2})\omega _0,`$ (26)
within the approximation (23), the differential equation becomes that for the parabolic cylinder functions:

$$\left[\frac{d^2}{dz^2}+\left(\nu ^{}+\frac{1}{2}\frac{z^2}{4}\right)\right]C_{m,1}^{}=0.$$
(27)
If we take as the two linearly independent solutions the standard forms $`D_\nu ^{}(z)`$ and $`D_\nu ^{}(z)`$ , the former must be rejected as it diverges for $`z\mathrm{}`$. We accordingly write

$$C_{m,1}^{}=A^{}(1)^{\mathrm{}^{}}D_\nu ^{}(z),$$
(28)
where $`A^{}`$ is the constant in Eq. (20), and the additional factor $`(1)^{\mathrm{}^{}}`$, where we shall define $`\mathrm{}^{}`$ shortly, is another constant introduced for later convenience.

As $`z\mathrm{}`$, $`D_\nu ^{}(z)(z)^\nu ^{}e^{z^2/4}`$, and one can show with a little work that \[modulo the factor $`(1)^{\mathrm{}^{}}`$\] this is indeed the DPI form for $`C_{m,1}^{}`$ with the approximation (23) . For $`z+\mathrm{}`$, on the other hand,

$`D_\nu ^{}(z)`$ $``$ $`\mathrm{cos}(\pi \nu ^{})z^\nu ^{}e^{z^2/4}\left(1{\displaystyle \frac{\nu ^{}(\nu ^{}1)}{2z^2}}+\mathrm{}\right)`$ (30)
$`+{\displaystyle \frac{\sqrt{2\pi }}{\mathrm{\Gamma }(\nu ^{})}}z^{\nu ^{}1}e^{z^2/4}\left(1+{\displaystyle \frac{(\nu ^{}+1)(\nu ^{}+2)}{2z^2}}+\mathrm{}\right).`$
Note that this form has both decaying and growing components. In fact, the latter component vanishes only if $`\nu ^{}`$ is a positive a integer. As explained earlier, it is essential for our DPI solution $`C_m^{}`$ to contain a growing component. We therefore allow for its presence by writing

$$\nu ^{}=\mathrm{}^{}+\frac{ϵ^{}}{\omega _0},$$
(31)
where $`\mathrm{}^{}`$ is a positive integer, and $`ϵ^{}`$ is a shift defined to lie in the interval $`(1/2,1/2)\omega _0`$. In fact, we expect $`ϵ^{}`$ to be very close to zero for any state in which there is large probaility of finding the particle in the left well. We can make the vanishing of the growing component in $`C_{m,1}^{}`$ more manifest by writing

$$\frac{1}{\mathrm{\Gamma }(\nu ^{})}=\frac{\mathrm{sin}(\pi \nu ^{})}{\pi }\mathrm{\Gamma }(1+\nu ^{})(1)^{\mathrm{}^{}+1}(\mathrm{}^{}!)\frac{ϵ^{}}{\omega _0}.$$
(32)
Combining Eqs. (28)–(32), we thus find that for $`m`$ beyond $`m_b^{}`$,

$$C_{m,1}^{}A^{}\left(\mathrm{cos}\frac{\pi ϵ^{}}{\omega _0}z^\nu ^{}e^{z^2/4}\sqrt{2\pi }(\mathrm{}^{}!)\frac{ϵ^{}}{\omega _0}z^{\nu ^{}1}e^{z^2/4}\right).$$
(33)
### E DPI form in ordinary forbidden region

The next step is to write the solution for $`C_m^{}`$ in such a way that it holds in the entire region $`m_b^{}<m<m_c^{}`$ . For $`C_{m,2}^{}`$, Eq. (21) already meets this demand, since just as at the turning point $`m_a^{}`$, $`q_2`$ stays imaginary and negative as $`m`$ passes through $`m_c^{}`$. For $`C_{m,1}^{}`$, on the other hand, Eq. (33) only holds in a region where the parabolic approximation to $`U_0`$ is good, and may not hold near $`m_c^{}`$. It is clear, however, that the wavevector $`q_1`$ associated with $`C_{m,1}^{}`$ is again given by $`i\kappa _1`$ with $`\kappa _1`$ given by Eq. (22). Hence it must be possible to write $`C_{m,1}^{}`$ in the DPI form

$$C_{m,1}^{}=|v_1(m)|^{1/2}\left[Q^{}\mathrm{exp}\left(i_{m_b^{}}^mq_1(m^{})𝑑m^{}\right)+R^{}\mathrm{exp}\left(i_{m_b^{}}^mq_1(m^{})𝑑m^{}\right)\right],$$
(34)
where $`Q^{}`$ and $`R^{}`$ are coefficients which we expect to be proportional to $`A^{}`$. To find these, let us first calculate the phase integral in the parabolic approximation. Since $`q_1=0`$ at $`m_b^{}`$, we have

$$E=_{\mathrm{sc}}(0,m_b^{})=_{\mathrm{sc}}(i\kappa _1(m),m),(m>m_b^{}).$$
(35)
Using Eq. (23), we get

$$\kappa _1(m)M_{}\omega _0[(mm_0)^2(m_b^{}m_0)^2)^{1/2}].$$
(36)
Therefore, if $`(mm_b^{})(m_b^{}m_0)`$, we obtain

$`\mathrm{exp}\left({\displaystyle _{m_b^{}}^m}\kappa _1(m^{})𝑑m^{}\right)`$ $`=`$ $`\left(2{\displaystyle \frac{mm_0}{m_b^{}m_0}}\right)^{\nu ^{}+{\scriptscriptstyle \frac{1}{2}}}e^{{\scriptscriptstyle \frac{1}{2}}\left(\nu ^{}+{\scriptscriptstyle \frac{1}{2}}\right)}`$ (38)
$`\times \mathrm{exp}\left({\displaystyle \frac{1}{2}}M_{}\omega _0(mm_0)^2\right),`$
where we have used the fact that

$$\frac{1}{2}M_{}\omega _0^2(m_b^{}m_0)^2=(\nu ^{}+\frac{1}{2})\omega _0.$$
(39)
Next we note that (a) from the definition of $`z`$, Eq. (25), the last exponential in Eq. (38) is nothing but $`\mathrm{exp}(z^2/4)`$, and that (b) $`|v_1(m)|\kappa _1(m)/M_{}\omega _0(mm_0)`$. Since $`mm_0z`$, it follows that the firat term in Eq. (34) varies as $`z^\nu ^{}\mathrm{exp}(z^2/4)`$, i.e., precisely as the first term in Eq. (33). Thus $`Q^{}`$ is indeed proportional to $`A^{}`$, and a little algebra plus the use of Eq. (39) gives

$$Q^{}=\left(\frac{\omega _0}{2M_{}}\right)^{1/4}e^{{\scriptscriptstyle \frac{1}{2}}\left(\nu ^{}+{\scriptscriptstyle \frac{1}{2}}\right)}(\nu ^{}+\frac{1}{2})^{{\scriptscriptstyle \frac{1}{2}}\left(\nu ^{}+{\scriptscriptstyle \frac{1}{2}}\right)}\mathrm{cos}\frac{\pi ϵ^{}}{\omega _0}A^{}\alpha ^{}A^{}.$$
(40)
In the same way, we can show that

$$R^{}=\left(\frac{\omega _0}{2M_{}}\right)^{1/4}e^{{\scriptscriptstyle \frac{1}{2}}\left(\nu ^{}+{\scriptscriptstyle \frac{1}{2}}\right)}(\nu ^{}+\frac{1}{2})^{{\scriptscriptstyle \frac{1}{2}}\left(\nu ^{}+{\scriptscriptstyle \frac{1}{2}}\right)}\sqrt{2\pi }(\mathrm{}^{}!)\frac{ϵ^{}}{\omega _0}A^{}\beta ^{}\frac{ϵ^{}}{\omega _0}A^{}.$$
(41)
The definitions of the factors $`\alpha ^{}`$ and $`\beta ^{}`$ which we have introduced for later convenience can be read off these equations. Note that since we took $`A^{}`$ to be real, $`Q^{}`$ and $`R^{}`$ are also real.

To summarize where we are, the complete DPI solution for $`C_m^{}`$ in the region $`m_b^{}<m<m_c^{}`$ is given by the sum of Eqs. (21) and (34). The terms in $`B^{}`$ and $`R^{}`$ are exponentially growing with increasing $`m`$, while the term in $`Q^{}`$ is exponentially decreasing. The next step is to connect this solution to the DPI form in the region $`m_c^{}<m`$. As already stated, the turning point $`m_c^{}`$ is the irregular one under the barrier, where exponentially growing and decaying wavefunctions acquire oscillatory envelopes as it is crossed.

Before we use the connection formulas at $`m_c^{}`$, it is useful to see how the quasiclassical wavevector behaves near this point. Since $`\mathrm{cos}q(m_c^{})=\mathrm{cos}q^{}=t_1(m_c^{})/4t_2(m_c^{})`$, it follows that the discriminant in Eq. (17) vanishes at $`m_c^{}`$, and both $`q_1`$ and $`q_2`$ tend to the same value $`i\kappa _c^{}`$, where

$$\mathrm{cosh}\kappa _c^{}=(|t_1|/4t_2)_{m=m_c^{}}.$$
(42)
As we cross $`m_c^{}`$, the discriminant in Eq. (17) becomes negative and $`\mathrm{cos}q(m)`$ \[and therefore $`q(m)`$\] becomes complex. We separate $`q(m)`$ into its real and imaginary parts, and write two distinct solutions as

$$q_{d,g}(m)=\pm i\kappa (m)+\chi (m),(m>m_c^{})$$
(43)
where both $`\kappa `$ and $`\chi `$ are real and positive. The subscripts $`d`$ and $`g`$ stand for ‘decaying’ and ‘growing’.

### F DPI form in oscillatory forbidden region

The connection formulas to be used at $`m_c^{}`$ were derived in Sec. IV of Ref. . The parts of $`C_m^{}`$ multiplying $`B^{}`$, $`Q^{}`$, and $`R^{}`$ correspond respectively to the cases there labelled $`(\sigma _1,\sigma _2)=(+1,1)`$, $`(1,+1)`$, and $`(1,1)`$. The $`B^{}`$ part is given by

$$C_{m,2}^{}=B^{}[\frac{1}{\sqrt{s_g(m)}}\mathrm{exp}(i_{m_c^{}}^mq_g(m^{})dm^{}\frac{\pi }{2})+\mathrm{c}.\mathrm{c}.],$$
(44)
where

$`s_g(m)`$ $`=`$ $`8t_2(m)\mathrm{sinh}\kappa (m)\mathrm{sin}\chi (m)\mathrm{sin}q_g(m)`$ (45)
$`=`$ $`8t_2\mathrm{sinh}\kappa \mathrm{sin}\chi (\mathrm{sin}\chi \mathrm{cosh}\kappa i\mathrm{cos}\chi \mathrm{sinh}\kappa ).`$ (46)
Likewise, the coefficient of $`Q^{}`$, which we call part 1a, connects to

$$C_{m,1a}^{}=e^\mathrm{\Gamma }^{}Q^{}[\frac{1}{\sqrt{s_d(m)}}\mathrm{exp}\left(i_{m_c^{}}^mq_d(m^{})dm^{}\right)+\mathrm{c}.\mathrm{c}.],$$
(47)
where

$$s_d(m)=8t_2(m)\mathrm{sinh}\kappa (m)\mathrm{sin}\chi (m)\mathrm{sin}q_d(m)=[s_g(m)]^{},$$
(48)
and $`\mathrm{\Gamma }^{}`$ is the phase integral,

$$\mathrm{\Gamma }^{}=_{m_b^{}}^{m_c^{}}\kappa (m^{})𝑑m^{},$$
(49)
which we acquire in changing the lower limits of the $`m`$ integrals from $`m_b^{}`$ to $`m_c^{}`$. Lastly, the term in $`R^{}`$, which we call part 1b, connects to

$$C_{m,1b}^{}=\frac{1}{2}e^\mathrm{\Gamma }^{}R^{}[\frac{1}{\sqrt{s_g(m)}}\mathrm{exp}\left(i_{m_c^{}}^mq_g(m^{})dm^{}\right)+\mathrm{c}.\mathrm{c}.].$$
(50)
Equations (44), (47), and (50), give us the complete wavefunction $`C_m^{}`$ in the central region near $`m=0`$. To simplify the writing, we denote

$$\mathrm{\Phi }_{\lambda _1\lambda _2}^{}(m)=_{m_c^{}}^m[\lambda _1\kappa (m^{})+i\lambda _2\chi (m^{})]𝑑m^{},$$
(51)
where $`\lambda _1`$ and $`\lambda _2`$ can be $`\pm 1`$ independently. In other words, the subscripts on $`\mathrm{\Phi }^{}`$ give the signs of the real and imaginary parts \[$`\mathrm{\Phi }_{++}^{}`$ is the integral of $`(\kappa +i\chi )`$, $`\mathrm{\Phi }_+^{}`$ that of $`(\kappa +i\chi )`$, etc.\]. The complete DPI solution for $`C_m^{}`$ can then be written as

$$C_m^{}=\left[e^\mathrm{\Gamma }^{}Q^{}\frac{e^{\mathrm{\Phi }_+^{}(m)}}{\sqrt{s_g^{}(m)}}+\left(\frac{1}{2}e^\mathrm{\Gamma }^{}R^{}iB^{}\right)\frac{e^{\mathrm{\Phi }_{++}^{}(m)}}{\sqrt{s_g(m)}}\right]+\mathrm{c}.\mathrm{c}.$$
(52)
The wavefunction from the right, $`C_m^{\prime \prime }`$ can now be written down at once. We define quantities with double primes in exact correspondence with those for $`C_m^{}`$. The analog of Eq. (52) is then

$$C_m^{\prime \prime }=\left[e^{\mathrm{\Gamma }^{\prime \prime }}Q^{\prime \prime }\frac{e^{\mathrm{\Phi }_+^{\prime \prime }(m)}}{\sqrt{s_g^{}(m)}}+\left(\frac{1}{2}e^{\mathrm{\Gamma }^{\prime \prime }}R^{\prime \prime }iB^{\prime \prime }\right)\frac{e^{\mathrm{\Phi }_{++}^{\prime \prime }(m)}}{\sqrt{s_g(m)}}\right]+\mathrm{c}.\mathrm{c}.$$
(53)
The only issue requiring any thought is what the sign suffixes in $`\mathrm{\Phi }^{\prime \prime }`$ should mean. By defining a new variable $`n=m`$, so that the problem for $`C_m^{\prime \prime }`$ becomes completely isomorphic to that for $`C_m^{}`$, and then transforming back to $`m`$, one can show that

$`\mathrm{\Gamma }^{\prime \prime }`$ $`=`$ $`{\displaystyle _{m_c^{\prime \prime }}^{m_b^{\prime \prime }}}\kappa (m^{})𝑑m^{},`$ (54)
$`\mathrm{\Phi }_+^{\prime \prime }(m)`$ $`=`$ $`{\displaystyle _m^{m_c^{\prime \prime }}}[\kappa (m^{})+i\chi (m^{})]𝑑m^{},`$ (55)
etc. Note that since $`m<m_c^{\prime \prime }`$ in the center, and $`m_c^{\prime \prime }<m_b^{\prime \prime }`$, these integrals are written so that the lower limit is less than the upper limit. Thus the suffixes on $`\mathrm{\Phi }^{\prime \prime }`$ give the true signs of its real and imaginary parts.

### G Matching of left and right wavefunctions

It remains to see if Eqs. (52) and (53) are the same function. We note that $`\mathrm{\Phi }_+^{\prime \prime }(m)`$ has the same integrand as $`\mathrm{\Phi }_+^{}(m)`$ if $`m`$ is taken to be the upper limit for both integrals. We further note that the remaining $`m`$ dependence in both terms is $`(s_g^{})^{1/2}`$. Similar remarks apply to the $`\mathrm{\Phi }_+^{}`$ and $`\mathrm{\Phi }_+^{\prime \prime }`$ terms. Thus, we conclude that $`C_m^{}`$ will equal $`C_m^{\prime \prime }`$ if the following conditions are obeyed:

$`e^{\mathrm{\Gamma }^{\prime \prime }}Q^{\prime \prime }e^{\mathrm{\Phi }_+^{\prime \prime }(m)}`$ $`=`$ $`\left(\frac{1}{2}e^\mathrm{\Gamma }^{}R^{}+iB^{}\right)e^{\mathrm{\Phi }_+^{}(m)},`$ (56)
$`e^\mathrm{\Gamma }^{}Q^{}e^{\mathrm{\Phi }_+^{}(m)}`$ $`=`$ $`\left(\frac{1}{2}e^{\mathrm{\Gamma }^{\prime \prime }}R^{\prime \prime }+iB^{\prime \prime }\right)e^{\mathrm{\Phi }_+^{\prime \prime }(m)}.`$ (57)
To simplify these conditions we note that

$`\mathrm{\Phi }_+^{}(m)\mathrm{\Phi }_+^{\prime \prime }(m)`$ $`=`$ $`{\displaystyle _{m_c^{}}^{m_c^{\prime \prime }}}[\kappa (m)i\chi (m)]𝑑m,`$ (58)
$``$ $`\mathrm{\Gamma }_ci\mathrm{\Lambda }_c,`$ (59)
where $`\mathrm{\Gamma }_c`$ and $`\mathrm{\Lambda }_c`$ are the real and imaginary parts of the integral, and the subscript ‘c’ indicates that the integrals extend over the central region of $`m`$. Equations (56) and (57) can now be written as

$`Q^{\prime \prime }`$ $`=`$ $`\left(\frac{1}{2}e^\mathrm{\Gamma }^{}R^{}+iB^{}\right)e^{\mathrm{\Gamma }^{\prime \prime }}e^{\mathrm{\Gamma }_ci\mathrm{\Lambda }_c},`$ (60)
$`Q^{}`$ $`=`$ $`\left(\frac{1}{2}e^{\mathrm{\Gamma }^{\prime \prime }}R^{\prime \prime }+iB^{\prime \prime }\right)e^\mathrm{\Gamma }^{}e^{\mathrm{\Gamma }_ci\mathrm{\Lambda }_c}.`$ (61)
If we recall that \[see Eqs. (40) and (41)\] $`Q^{}`$ and $`R^{}`$ are proportional to $`A^{}`$, and likewise for $`Q^{\prime \prime }`$ and $`R^{\prime \prime }`$, these equations are two complex equations in the four real quantities $`A^{}`$, $`B^{}`$, $`A^{\prime \prime }`$ and $`B^{\prime \prime }`$. To solve them we first note that the imaginary parts on the right hand sides must vanish. This yields

$`B^{}`$ $`=`$ $`\frac{1}{2}R^{}e^\mathrm{\Gamma }^{}\mathrm{tan}\mathrm{\Lambda }_c,`$ (62)
$`B^{\prime \prime }`$ $`=`$ $`\frac{1}{2}R^{\prime \prime }e^{\mathrm{\Gamma }^{\prime \prime }}\mathrm{tan}\mathrm{\Lambda }_c.`$ (63)
Substituting these in Eqs. (60) and (61), we obtain after some simplification

$`R^{}=2e^{\mathrm{\Gamma }_G}\mathrm{cos}\mathrm{\Lambda }_cQ^{\prime \prime },`$ (64)
$`R^{\prime \prime }=2e^{\mathrm{\Gamma }_G}\mathrm{cos}\mathrm{\Lambda }_cQ^{},`$ (65)
where $`\mathrm{\Gamma }_G`$ is the total Gamow factor

$$\mathrm{\Gamma }_G=_{m_b^{}}^{m_b^{\prime \prime }}\kappa _1(m)𝑑m,$$
(66)
and the subscript ‘1’ on $`\kappa `$ is to remind us that we must use the imaginary part of the wavevector that goes to zero at the turning points $`m_b^{}`$ and $`m_b^{\prime \prime }`$.

### H The eigenvalue condition, and the diabolo

The simplest way of solving Eqs. (64) and (65) is terms of the ratios $`\alpha ^{}`$ and $`\beta ^{}`$ defined in Eqs. (40) and (41), and the analogous ratios $`\alpha ^{\prime \prime }`$ and $`\beta ^{\prime \prime }`$. Equating the products of the left hand and right hand sides, and simplifying a little, we get

$$ϵ^{}ϵ^{\prime \prime }=4\frac{\alpha ^{}\alpha ^{\prime \prime }}{\beta ^{}\beta ^{\prime \prime }}\omega _0\omega _{0+}\mathrm{cos}^2\mathrm{\Lambda }_ce^{2\mathrm{\Gamma }_G}.$$
(67)
This is our eigenvalue condition. To understand it better, we first note that the right hand side is exponentially small on account of the square of the Gamow factor $`e^{\mathrm{\Gamma }_G}`$. Thus, ignoring for the moment the possibility that $`\mathrm{cos}\mathrm{\Lambda }_c`$ may vanish, either both $`ϵ^{}`$ and $`ϵ^{\prime \prime }`$ must be of order $`e^{\mathrm{\Gamma }_G}`$, or at the other extreme, one must be of order unity, and the other of order $`e^{2\mathrm{\Gamma }_G}`$. Suppose that $`ϵ^{}=O(e^{2\mathrm{\Gamma }_G})`$, and $`ϵ^{\prime \prime }=O(1)`$. Let us take $`A^{}=1`$. Then from Eqs. (40) and (41), we see that $`Q^{}=O(1)`$, while $`R^{}=O(e^{2\mathrm{\Gamma }_G})`$. It then follows from Eqs. (64) and (65) that $`Q^{\prime \prime }R^{\prime \prime }=O(e^{\mathrm{\Gamma }_G})`$, and in turn from the double primed analogs of Eqs. (40) and (41) that $`A^{\prime \prime }=O(e^{\mathrm{\Gamma }_G})`$. Lastly, since $`\mathrm{\Gamma }^{\prime \prime }<\mathrm{\Gamma }_G`$, Eq. (63) implies that $`B^{\prime \prime }`$ is also exponentially small . Hence, it follows that the entire wavefunction on the right hand side of the well, $`C_m^{\prime \prime }`$ is exponentially small compared to the left hand part $`C_m^{}`$. In other words, there is negligible mixing of the states in the left and right hand well. This is exactly what we expect when the energies of the two states in the absence of tunneling differ by much more than the tunneling matrix element itself.

The more interesting case, therefore, is that in which both $`ϵ^{}`$ and $`ϵ^{\prime \prime }`$ are of order $`e^{\mathrm{\Gamma }_G}`$. In the defining equations for the $`\alpha `$’s and $`\beta `$’s, Eqs. (40) and (41), we can then neglect the $`ϵ`$’s to very good approximation. This yields

$`{\displaystyle \frac{\alpha ^{}}{\beta ^{}}}`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2\pi }\mathrm{}^{}!}}\left(\mathrm{}^{}+\frac{1}{2}\right)^{\mathrm{}^{}+{\scriptscriptstyle \frac{1}{2}}}e^{(\mathrm{}^{}+{\scriptscriptstyle \frac{1}{2}})},`$ (68)
$`=`$ $`{\displaystyle \frac{g_{\mathrm{}^{}}}{2\pi }},`$ (69)
where $`g_n`$ is the standard curvature correction  in the phase integral expression for the tunnel splitting (see e.g., Eq. (4.10) of paper I). This quantity tends to 1 rapidly as $`n`$ gets large: $`g_0=(\pi /e)^{1/2}1.075`$, $`g_11.028`$, $`g_21.017`$, $`\mathrm{}`$. The right hand side of Eq. (67) can thus be written as $`\mathrm{\Delta }^2(\mathrm{}^{},\mathrm{}^{\prime \prime })/4`$, where

$$\mathrm{\Delta }(\mathrm{}^{},\mathrm{}^{\prime \prime })=\frac{2}{\pi }(g_{\mathrm{}^{}}g_{\mathrm{}^{\prime \prime }})^{1/2}(\omega _0\omega _{0+})^{1/2}e^{\mathrm{\Gamma }_G}\mathrm{cos}\mathrm{\Lambda }_c.$$
(70)
We further define

$`ϵ`$ $`=`$ $`\frac{1}{2}(ϵ^{}+ϵ^{\prime \prime })=E\frac{1}{2}\left(E_{}+E_++(\mathrm{}^{}+\frac{1}{2})\omega _0+(\mathrm{}^{\prime \prime }+\frac{1}{2})\omega _{0+}\right),`$ (71)
$`\delta `$ $`=`$ $`ϵ^{\prime \prime }ϵ^{}=\left(E_{}E_++(\mathrm{}^{}+\frac{1}{2})\omega _0(\mathrm{}^{\prime \prime }+\frac{1}{2})\omega _{0+}\right).`$ (72)
With these definitions, Eq. (67) can be rewritten as

$$ϵ=\pm \frac{1}{2}[\delta ^2+\mathrm{\Delta }^2(\mathrm{}^{},\mathrm{}^{\prime \prime })]^{1/2}.$$
(73)
Equation (73) is the complete formal solution to the problem of tunneling in an asymmetric double well. Along with Eqs. (70)–(72), it is the analog of the general phase integral formula for the tunnel splitting in a symmetric double well that we found in I \[see Eq. (4.38) there\]. Since there is no great need for having the final answer in simple closed form for the specfic problem (1), and since the general procedure is fully explained in Sec. IV.E of I, we do not bother to extract the singular $`\mathrm{ln}J`$ parts of the $`\mathrm{\Gamma }_G`$ integral.

It is immediately obvious that the eigenvalues in Eq. (73) are what we would get from a two level system Hamiltonian

$$_{\mathrm{TLS}}=\frac{1}{2}\left(\begin{array}{cc}\delta & \mathrm{\Delta }(\mathrm{}^{},\mathrm{}^{\prime \prime })\\ \mathrm{\Delta }(\mathrm{}^{},\mathrm{}^{\prime \prime })& \delta \end{array}\right),$$
(74)
which is of course, just what we would expect. The quantity $`ϵ`$ is the energy measured from a convenient reference point, while $`\delta `$, which depends on the fields $`h_x`$, $`h_z`$, and the quantum numbers $`\mathrm{}^{}`$ and $`\mathrm{}^{\prime \prime }`$ of the states whose mixing is being examined, is the offset between these energy levels in the absence of tunneling. Equation (70) gives the tunneling amplitude between these levels when the offset is small, i.e., when the two levels are in approximate degeneracy. Note that although this amplitude is defined even for relatively large offsets—offsets comparable to the intrawell spacings $`\omega _{0\pm }`$—and indeed is not very sensitive to the value of the offset, the concept of tunneling is physically sensible and useful only when the offset is comparable to or less than the amplitude $`\mathrm{\Delta }`$. If $`\delta \mathrm{\Delta }`$, we get $`ϵ\pm \delta /2`$, i.e., $`ϵ^{\prime \prime }\delta `$, $`ϵ^{}\mathrm{\Delta }^2/\delta `$, or the other way around. Then by the argument given after Eq. (67), the mixing between the wells is negligible.

One may wonder if the above conclusion does not invalidate the entire calculation. After all, we defined the shifts $`ϵ^{}`$ and $`ϵ^{\prime \prime }`$ assuming the wells were parabolic. Surely the corrections to the energies from cubic and higher order corrections to the potential are far larger than $`\mathrm{\Delta }`$. Since the offset between the levels must be tuned with exponential sensitivity, should not we know the energies of the levels before tunneling to the same sensitivity? The answer is no. The reason can be seen from Eq. (33). As we have seen, the amplitude of the growing part of the wavefunction plays a key role in the tunneling. From Eq. (33), we see that this amplitude is only linearly dependent on $`ϵ^{}`$. Hence, a small error in locating the absolute position of the level has little effect on the computed value of $`\mathrm{\Delta }`$. To say it another way, even though we know that the bottoms of the wells must be tuned to exponential accuracy to get significant mixing between the two wells, we cannot and need not determine the center of gravity of the two levels, $`ϵ`$, to the same accuracy to determine the tunneling amolitude itself. This feature is also present in the symmetric case studied in I, although there it is not so apparent, since Herring’s formula gives $`\mathrm{\Delta }`$ directly without making reference to the absolute energy level.

Secondly, it should be noted that Eq. (73) is nothing but the equation for the diabolo. The splitting vanishes only when $`\delta `$ and $`\mathrm{\Delta }`$ both vanish, which furnish the two conditions required to determine the diabolical point. Since both $`\delta `$ and $`\mathrm{\Delta }`$ will in general have linear terms in the deviation from the diabolical point, the energy surface is a double cone as asserted earlier.

Thirdly, let us ask if we recover the results of paper I in the symmetric case, i.e., when $`h_z=0`$. Then $`E_+=E_{}`$, $`M_+=M_{}`$, and $`\omega _0=\omega _{0+}`$. The only sensible case is $`\mathrm{}^{}=\mathrm{}^{\prime \prime }n`$, so that $`ϵ^{}=ϵ^{\prime \prime }`$, and $`\delta =0`$. The splitting is (up to an irrelevant sign) $`\mathrm{\Delta }(n,n)`$, which is precisely the tunnel splitting $`\mathrm{\Delta }_n`$ computed in paper I. In addition, however, we now have more explicit information about the wavefunction. Proceeding as before, we see that $`Q^{}A^{}1`$, $`R^{}e^{\mathrm{\Gamma }_G}`$, and $`B^{}e^{(\mathrm{\Gamma }_G\mathrm{\Gamma }^{})}`$. (The double primed quantities are equal to their single primed counterparts.) The conclusion about $`B^{}`$ is totally consistent with the approach in paper I, which was based on Herring’s formula. There one takes the wavefunctions as $`(C_{m,d}\pm C_{m,d})/\sqrt{2}`$, where $`C_{m,d}`$ is the wavefunction of a state localized in the left well, and which decays away from that well in both directions. If we equate $`C_{m,d}`$ with $`C_{m,1a}^{}`$ (and, therefore, $`C_{m,d}`$ with $`C_{m,1a}^{\prime \prime }`$) in the central region, then we see that $`B^{}`$, which by Eq. (21) gives the magnitude of the growing part of $`C_m^{}`$ at $`m=m_c^{}`$, is also the order of magnitude of $`C_{m,1a}^{\prime \prime }`$, the decaying part of $`C_m^{\prime \prime }`$ at $`m=m_c^{}`$.

### I What is the origin of the diabolical points?

Lastly, it is extremely instructive to examine the problem when $`\mathrm{\Delta }=0`$, i.e., $`\mathrm{cos}\mathrm{\Lambda }_c=0`$, without necessarily imposing the condition $`\delta =0`$, for this gives insight into what causes the quenching of the tunnel splitting. Taking the imaginary parts of Eqs. (60) and (61) we see directly that we must have $`R^{}=R^{\prime \prime }=0`$, and that

$$Q^{\prime \prime }=\pm B^{}e^{\mathrm{\Gamma }_c+\mathrm{\Gamma }^{\prime \prime }},Q^{}=\pm B^{\prime \prime }e^{\mathrm{\Gamma }_c+\mathrm{\Gamma }^{}}.$$
(75)
Going back to Eq. (41), we see that $`R^{}=0`$ requires either $`ϵ^{}=0`$ or $`A^{}=0`$, and likewise for $`R^{\prime \prime }`$, $`ϵ^{\prime \prime }`$, and $`A^{\prime \prime }`$. If $`\delta 0`$, then both $`ϵ^{}`$ and $`ϵ^{\prime \prime }`$ can not vanish, and the only solution is $`ϵ^{}0`$, $`ϵ^{\prime \prime }=0`$, $`A^{\prime \prime }=Q^{\prime \prime }=B^{}=0`$, (or the one obtained by interchanging single and double primes.) The only non-zero coefficients are $`A^{}`$, $`Q^{}`$, and $`B^{\prime \prime }`$. From Eqs. (40) and (75), we see that only one of these coefficients is independent, which can only be fixed by normalization. Thus, we see that the part $`C_{m,2}^{\prime \prime }`$ proportional to $`B^{\prime \prime }`$ should really be regarded as the extreme right-hand tail of a state localized in the left well. If we denote this state by $`|L`$, and the wavefunction $`m|L`$ by $`C_L(m)`$, then ,

$$C_L(m)=\{\begin{array}{cc}A^{}|v_1(m)|^{1/2}\mathrm{exp}\left(i^mq_1(m^{})𝑑m^{}\right),\hfill & m<m_a^{}\text{,}\hfill \\ A^{}(1)^{\mathrm{}^{}}D_{\mathrm{}^{}}(z),\hfill & mm_a^{}\mathrm{to}mm_b^{}\text{,}\hfill \\ Q^{}|v_1(m)|^{1/2}\mathrm{exp}\left(i_{m_b^{}}^mq_1(m^{})𝑑m^{}\right),\hfill & m_b^{}<m<m_b^{\prime \prime }\text{,}\hfill \\ B^{\prime \prime }|v_1(m)|^{1/2}\mathrm{exp}\left(i_{m_b^{\prime \prime }}^mq_1(m^{})𝑑m^{}\right),\hfill & m_b^{\prime \prime }<m\text{.}\hfill \end{array}$$
(76)
\[Note that in the second line, we wrote $`D_{\mathrm{}^{}}`$, not $`D_\nu ^{}`$, and that in the third line, we have not bothered to write the oscillatory exponential continuation in the region $`m_c^{}<m<m_c^{\prime \prime }`$ correctly—see Eq. (47)—since our aim now is merely to indicate the general structure.\] We can define a right-hand function $`C_R(m)`$ analogously. Indeed it is now clear that the two state Hamiltonian (74) is a truncation of the full Hamiltonian (1) in the $`|L`$, $`|R`$ basis.

The above argument shows that for any $`H_z`$ ($`\delta 0`$), we can tune $`H_x`$ so that $`\mathrm{\Delta }`$ vanishes, at which point, the energy eigenfunction is like $`C_L(m)`$ \[or $`C_R(m)`$\], which is localized in one well and does not “see” the other well at all! In ordinary one dimensional quantum mechanics, this is of course impossible, since a wavefunction like $`C_L(m)`$ which continues decaying with increasing $`m`$ in the classically allowed region of the right hand well has the wrong sign of the curvature in that well. In fact, this argument does not depend on having $`\delta 0`$. If $`\delta =0`$ in addition to $`\mathrm{\Delta }=0`$, $`C_L(m)`$ and $`C_R(m)`$ are independent solutions of Schrödingers equation, as is any linear combination, since they are degenerate.

The above point of view helps elucidate the origin of the quenching more clearly. Indeed, it is better to think about the non-symmetric situation ($`H_z0`$) than the symmetric one ($`H_z=0`$). In continuum problems with a symmetric double well, Herring’s formula gives the splitting as proportional to $`[\psi _d(x)(d\psi _d^{}/dx)]_{x=0}`$, where $`\psi _d(x)`$ is a left-well localized wavefunction . It is tempting to think that the splitting in the spin problem vanishes because the oscillatory envelope in $`C_{m,d}`$ in the central region allows the discrete analog of $`\psi _d(x)`$ or $`\psi _d^{}(x)`$ to vanish at the midpoint. This reasoning is false, as one can see from a close examination of the symmetric case wavefunction, or even more clearly, by looking at the situation when $`\mathrm{\Delta }=0`$ but $`\delta 0`$. The condition $`\mathrm{\Delta }=0`$ can not be reduced to a local property of the wavefunction such as its value or its slope at a particular point. Rather it is the global property that the phase integral $`\mathrm{\Lambda }_c`$ be an odd integer times $`\pi /2`$. From this perspective, the quenching is perhaps better visualized as a manifestation of interfering Feynman trajectories and the Berry phase, even though the value of this phase is more easily found using the DPI method.

## III Application to Fe8

Let us now apply our general formalism to Fe8. The problem of greatest interest is the location of the diabolical points, and for that we need only solve the conditions $`\delta =\mathrm{\Delta }=0`$. We have already given formulas for the matrix elements and the critical curves in Eqs. (10)–(15). The problem that remains is to use these formulas to find $`\delta `$ and $`\mathrm{\Delta }(\mathrm{}^{},\mathrm{}^{\prime \prime })`$. To keep the problem tractable and obtain answers in closed form, we will assume that $`h_z`$ is small. Specifically, we will assume that the reduced field $`h_z`$ defined above Eq. (13) is formally of order $`1/\overline{J}`$. This enables us to evaluate the turning points and action integrals as expansions in powers of $`h_z`$. Also, it is convenient to carry out all calculations in terms of the reduced variable $`\mu `$.

The first step is to obtain $`\delta `$. For this, we need to analyze the critical curve $`U_0(\mu )`$. Its minima $`\mu _{0\pm }`$ are found to be located at

$$\mu _{0\pm }=\pm \mu _0+\frac{h_zh_x^2}{1h_x^2}+O(h_z^2),$$
(77)
where

$$\mu _0=(1h_x^2)^{1/2}.$$
(78)
The quantities $`E_\pm `$, $`\omega _{0\pm }`$, and $`M_\pm `$ defined through Eq. (23) are given by

$`E_\pm `$ $``$ $`U_0(\mu _{0\pm })=\left[h_x^2\pm 2\mu _0h_z+{\displaystyle \frac{h_z^2h_x^2}{1h_x^2}}\right],`$ (79)
$`\omega _{0\pm }`$ $`=`$ $`{\displaystyle \frac{2\lambda ^{1/2}\mu _0}{\overline{J}}}\left[1\pm {\displaystyle \frac{h_z}{2\mu _0}}\left({\displaystyle \frac{1}{\lambda }}+{\displaystyle \frac{1+2h_x^2}{1h_x^2}}\right)+O(h_z^2)\right],`$ (80)
$`M_\pm `$ $`=`$ $`{\displaystyle \frac{1}{2\lambda h_x^2}}\left[1{\displaystyle \frac{h_z}{\mu _0}}\left({\displaystyle \frac{1}{\lambda }}2\right)+O(h_z^2)\right].`$ (81)
Substituting these results in Eq. (72), we obtain

$$\delta (h_z,\mathrm{}^{},\mathrm{}^{\prime \prime })=4\mu _0h_z+\frac{2\sqrt{\lambda }\mu _0}{\overline{J}}(\mathrm{}^{}\mathrm{}^{\prime \prime })\frac{\sqrt{\lambda }h_z}{\overline{J}}(\mathrm{}^{}+\mathrm{}^{\prime \prime }+1)c_1(h_x)+O(\overline{J}^3),$$
(82)
where

$$c_1(h_x)=\frac{1h_x^2+\lambda (1+2h_x^2)}{\lambda (1h_x^2)}.$$
(83)
The next step is to evaluate $`\mathrm{\Delta }(\mathrm{}^{},\mathrm{}^{\prime \prime })`$, or more precisely $`\mathrm{\Lambda }_c`$, the imaginary part of the phase integral defined in Eq. (59), since that by itself locates the diabolical points. This in turn requires expressions for $`\chi (m)`$ and the points $`m_c^{}`$ and $`m_c^{\prime \prime }`$. To obtain $`\chi (m)`$, we return to the Hamilton-Jacobi equation Eq. (9), write $`q=\kappa +i\chi `$, and separate the equation into its real and imaginary parts. Eliminating $`\kappa (m)`$ from the two equations that result, we obtain a single equation for $`\chi (m)`$, which can be written as

$`4t_2(m)X^2g(m)X+{\displaystyle \frac{t_1^2(m)}{4t_2(m)}}=0;`$ (84)
$`g(m)=w(m)+2t_2(m)E,`$ (85)
where $`X\mathrm{cos}\chi (m)`$. \[We can find the equation obeyed by $`\kappa (m)`$ similarly, and we discover that it is again Eq. (84) with $`X=\mathrm{cosh}\kappa (m)`$.\]

What value of $`E`$ should we use in Eq. (84)? As stated earlier, the value of $`\mathrm{\Delta }`$ is relatively insensitive to small changes in the absolute position of $`E`$. Thus we certainly needn’t incorporate the exponentially small shifts caused by the tunneling itself. Secondly, the tunneling is relevant only when $`\delta \mathrm{\Delta }`$. Thus it suffices to set both $`ϵ^{}`$ and $`ϵ^{\prime \prime }`$ to 0 in this part of the calculation. Then, Eqs. (71), (79), and (80), we get

$`E`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(E_{}+E_++(\mathrm{}^{}+\frac{1}{2})\omega _0+(\mathrm{}^{\prime \prime }+\frac{1}{2})\omega _{0+}\right)`$ (86)
$`=`$ $`h_x^2+{\displaystyle \frac{\sqrt{\lambda }\mu _0}{\overline{J}}}(\mathrm{}^{}+\mathrm{}^{\prime \prime }+1)+O(\overline{J}^2).`$ (87)
Substituting this along with the formulas for $`w(m)`$, $`t_1(m)`$, and $`t_2(m)`$ in Eqs. (84) and (85), yields the following equation for $`X`$:

$`(1\lambda )(1\mu ^2)X^2[1+h_x^2\mu ^2\zeta (\mu )]X+{\displaystyle \frac{h_x^2}{1\lambda }}=0;`$ (88)
$`\zeta (\mu )={\displaystyle \frac{\sqrt{\lambda }}{\overline{J}}}[\mu (\mathrm{}^{\prime \prime }\mathrm{}^{})+\mu _0(\mathrm{}^{}+\mathrm{}^{\prime \prime }+1)].`$ (89)
We have separated out the term $`\zeta (\mu )`$ in Eq. (88) as it is of order $`\overline{J}^1`$ relative to the other terms. We can solve the quadratic equation and expand the result as a power series in $`\zeta `$. Recalling that $`X=\mathrm{cos}^2\chi `$, we obtain

$$\mathrm{cos}\chi (\mu )=\frac{h_x}{\sqrt{(1\lambda )(1\mu ^2)}}\left[1+\frac{\zeta (\mu )}{2(1h_x^2\mu ^2)}+O(\zeta ^2)\right].$$
(90)
By setting $`\mathrm{cos}\chi =1`$, we can determine the points $`\mu _c^{}`$ and $`\mu _c^{\prime \prime }`$, for which it is now more convenient to write $`\mu _{c\pm }`$ instead. If the correction $`\zeta `$ were absent, it is easy to see that these points would be at $`\pm \mu _{c0}`$, where

$$\mu _{c0}=[(1\lambda h_x^2)/(1\lambda )]^{1/2}.$$
(91)
Now, with $`\zeta 0`$, we can find $`\mu _{c\pm }`$ as a series in $`\overline{J}^1`$. Up to leading corrections, we get

$$\mu _{c\pm }=\pm \mu _{c0}\frac{1}{2\sqrt{\lambda }\mu _{c0}\overline{J}}[\mu _{c0}(\mathrm{}^{\prime \prime }\mathrm{}^{})\pm \mu _0(\mathrm{}^{}+\mathrm{}^{\prime \prime }+1)].$$
(92)
The phase integral $`\mathrm{\Lambda }_c`$ is given by

$$\mathrm{\Lambda }_c=\overline{J}_{\mu _c}^{\mu _{c+}}\chi (\mu )𝑑\mu .$$
(93)
We will only evaluate this accurate to terms of order $`\overline{J}^0`$. Let us first consider the corrections entailed by replacing the limits by $`\pm \mu _{c0}`$. It is not difficult to show that without the $`\zeta `$ term in Eq. (90),

$$\chi (\mu )\frac{\sqrt{1\lambda }}{h_x}(\mu _{c0}^2\mu ^2)^{1/2}\mathrm{as}\mu \mu _{c0}.$$
(94)
Inclusion of the $`\zeta `$ correction does not change the square root approach to zero, and can not change the coefficient to leading order in $`\overline{J}^1`$. From Eq. (92), $`\mu _{c\pm }`$ differs from $`\pm \mu _{c0}`$ by terms of order $`\overline{J}^1`$, so ignoring these shifts in the limits of the integral causes an error of order $`\overline{J}^{3/2}`$ in the integral. Hence, we have

$$\mathrm{\Lambda }_c\overline{J}_{\mu _{c0}}^{\mu _{c0}}\chi (\mu )𝑑\mu +O(\overline{J}^{1/2}).$$
(95)
The error is smaller than $`O(\overline{J}^0)`$, and so may be ignored. The remaining integral may be done exactly as in paper I. We write the solution to Eq. (90) by $`\chi =\chi _0+\mathrm{\Delta }\chi `$, where $`\chi _0`$ is the solution when the $`\zeta (\mu )`$ correction is ignored, and $`\mathrm{\Delta }\chi `$ is the $`O(\zeta )`$ correction. The $`\chi _0`$ term can be integrated by parts, and yields

$`\mathrm{\Lambda }_{c0}`$ $`=`$ $`2\overline{J}{\displaystyle _0^{\mu _{c0}}}\chi _0(\mu )𝑑\mu `$ (96)
$`=`$ $`\pi \overline{J}[1h_x(1\lambda )^{1/2}].`$ (97)
The $`\mathrm{\Delta }\chi `$ term yields

$`\mathrm{\Delta }\mathrm{\Lambda }_c`$ $`=`$ $`\overline{J}{\displaystyle _{\mu _{c0}}^{\mu _{c0}}}{\displaystyle \frac{\sqrt{1\mu _{c0}^2}}{\sqrt{\mu _{c0}^2\mu ^2}}}{\displaystyle \frac{\zeta (\mu )}{2(1h_x^2\mu ^2)}}𝑑\mu `$ (98)
$`=`$ $`(\mathrm{}^{}+\mathrm{}^{\prime \prime }+1){\displaystyle \frac{\pi }{2}},`$ (99)
where the result follows after an elementary integration by substitution. Note that the first term in $`\zeta (\mu )`$ integrates to 0 as it is odd. The total result for $`\mathrm{\Lambda }_c`$ is (reverting to unscaled variables)

$`\mathrm{\Lambda }_c`$ $`=`$ $`\mathrm{\Lambda }_{c0}+\mathrm{\Delta }\mathrm{\Lambda }_c`$ (100)
$`=`$ $`{\displaystyle \frac{\pi }{2}}\left[2J(\mathrm{}^{}+\mathrm{}^{\prime \prime })2J{\displaystyle \frac{H_x}{H_c\sqrt{1\lambda }}}\right].`$ (101)
It is now elementary to find the conditions for a diabolical point. The offset $`\delta `$ becomes zero at a certain value of $`h_z`$, and the tunneling amplitude $`\mathrm{\Delta }`$ becomes zero when $`\mathrm{\Lambda }_c`$ is an odd multiple of $`\pi /2`$. Adding arguments $`\mathrm{}^{}`$ and $`\mathrm{}^{\prime \prime }`$ to indicate the level numbers becoming degenerate, the diabolicity conditions are,

$`{\displaystyle \frac{H_z(\mathrm{}^{},\mathrm{}^{\prime \prime })}{H_c}}`$ $`=`$ $`{\displaystyle \frac{\sqrt{\lambda }(\mathrm{}^{\prime \prime }\mathrm{}^{})}{2J}}\left[1+{\displaystyle \frac{\sqrt{\lambda }}{4\mu _0J}}(\mathrm{}^{}+\mathrm{}^{\prime \prime }+1)c_1(h_x)+O(J^2)\right],`$ (102)
$`{\displaystyle \frac{H_x(\mathrm{}^{},\mathrm{}^{\prime \prime })}{H_c}}`$ $`=`$ $`{\displaystyle \frac{\sqrt{1\lambda }}{J}}\left[Jn\frac{1}{2}(\mathrm{}^{}+\mathrm{}^{\prime \prime }+1)\right],n=0,1,\mathrm{},2J(\mathrm{}^{}+\mathrm{}^{\prime \prime }+1).`$ (103)
Note that what appears in these equations is $`J`$, not $`\overline{J}`$. Secondly, the restrictions on the integer $`n`$ in Eq. (103) are obtained, as explained in I, by demanding that $`\mathrm{\Lambda }_c`$ be positive.

We obtain the perfect centered rectangular lattice of diabolical points introduced in Sec. I if we ignore the $`c_1`$ term in Eq. (102), and also the restriction $`H_z/H_c1`$ used to perform the calculation. From the viewpoint of our DPI calculations, this exactness is somewhat mysterious. However, by the duality argument of Sec. I, if formulas (102) and (102) are correct for small $`H_z`$ and large $`H_x`$, then they also yield diabolical points for large $`H_z`$ and small $`H_x`$. Of course if the former set of points corresponds to low lying levels, i.e., small values of $`\mathrm{}^{}`$ and $`\mathrm{}^{\prime \prime }`$, to which our analysis applies, the latter set of points corresponds to rather highly excited levels, to which the analysis does not apply prima facie. It is nevertheless surprising that the formulas should fit together so neatly. It is also somewhat surprising that the experimentally determined Fe8 diabolical points \[for the cases $`(\mathrm{}^{},\mathrm{}^{\prime \prime })=(0,0)`$, $`(0,1)`$, and $`(0,2)`$\] should lie on a centered rectangular structure so closely, since, as is known, the value of the measured $`H_x`$ period on the $`H_z=0`$ line is almost 50% different from that predicted by Eq. (103). This feature is at present understood only on the basis of numerical diagonalization of the spin Hamiltonian including fourth order terms in $`𝐉`$. An analytic approach to this aspect of the problem remains for the future.

###### Acknowledgements.

This work is supported by the NSF via grant number DMR-9616749. I am indebted to Wolfgang Wernsdorfer and Jacques Villain for useful discussions and correspondence about Fe8.


format: text

text:
# Cutoff dependence of the Casimir effect

## Abstract

The problem of calculating the Casimir force on two conducting planes by means of the stress tensor is examined. The evaluation of this quantity is carried out using an explicit regularization procedure which has its origin in the underlying (2+1) dimensional Poincaré invariance of the system. The force between the planes is found to depend on the ratio of two independent cutoff parameters, thereby rendering any prediction for the Casimir effect an explicit function of the particular calculational scheme employed. Similar results are shown to obtain in the case of the conducting sphere.

In 1948 Casimir  first predicted that two infinite parallel plates in vacuum would attract each other. This remarkable result has its origin in the zero point energy of the electromagnetic field. While the latter is highly divergent, the change associated with this quantity for specific plate configurations has been found in numerous calculations to be finite as well as cutoff dependent and thus in principle observable. Early work to detect this small effect  was characterized by relatively large experimental uncertainties which left the issue in some doubt. More recent efforts  have provided quite remarkable data, but are based on a different geometry from that of Casimir. Since a rigorous theoretical calculation has never been carried out for the latter configuration, there remains room for skepticism as to whether the Casimir effect is as well established as is frequently asserted.

The most elementary calculation of the Casimir effect between two parallel conducting planes located at $`z=0`$ and $`z=a`$ employs a mode summation in the framework of a regularization which depends only on the frequency $`\omega _k=[𝐤^2+(\frac{n\pi }{a})^2]^{\frac{1}{2}}`$ where $`n=0,1,2,\mathrm{}`$. Upon combining the result obtained with the corresponding result for for the interval $`azL`$ where $`La`$ is the $`z`$-coordinate of a third conducting plane, a finite cutoff independent result $`F/A=\pi ^2/240a^4`$  is obtained for the Casimir pressure on the plate at $`z=a`$.

A considerably more elegant approach to this problem is that of Brown and Maclay  who employ an image method to calculate $`0|T^{\mu \nu }(x)|0`$. Thus they showed that the photon propagator in the presence of conducting planes at $`z=0`$ and $`z=a`$ could be expressed in terms of an infinite sum over the usual (i.e., $`\mathrm{}<z<\mathrm{}`$) photon propagator with the $`z`$-coordinate of each term in the sum displaced by an even multiple of $`a`$. Since the stress tensor for the electromagnetic case is given by

$$T^{\mu \nu }(x)=F^{\mu \alpha }F_\alpha ^\nu \frac{1}{4}g^{\mu \nu }F^{\alpha \beta }F_{\alpha \beta }$$
(1)
where

$$F^{\mu \nu }(x)=^\mu A^\nu (x)^\nu A^\mu (x),$$

it follows that upon taking appropriate derivatives with respect to the propagator arguments $`x`$ and $`x^{}`$ and invoking the limit $`xx^{}`$ a formal expression can be obtained for the vacuum expectation value of the stress tensor. On the basis of covariance arguments together with the divergence and trace free property of $`T^{\mu \nu }(x)`$ it was then found in ref. 5 that

$$0|T^{\mu \nu }(x)|0=(\frac{1}{4}g^{\mu \nu }\widehat{z}^\mu \widehat{z}^\nu )(\frac{1}{2\pi ^2a^4})\underset{n=1}{\overset{\mathrm{}}{}}n^4$$
(2)
where $`\widehat{z}^\mu `$ is the unit vector (0,0,1,0) in the $`z`$-direction normal to the conducting planes.

However, there is some reason to question whether this approach has adequately dealt with the divergences which invariably occur in Casimir calculations. One notes in particular that the result (2) is obtained only after an obviously singular $`n=0`$ term has been dropped from the sum which occurs in that equation. While one can argue as in  that such an $`a`$-independent term can be freely omitted since it is merely the usual subtraction of the large $`a`$ result, it is well to note that the entire sum over $`n`$ is required for a demonstration that the propagator satisfies correct boundary conditions at $`z=0,a`$. Moreover, as is shown in this work, an appropriately regularized form of (2) does not necessarily allow a separation into cutoff dependent terms and $`a`$-dependent terms, in contrast with the result found in . Of still greater import is the fact that more general regularizations than those usually considered in this calculation lead to an explicit cutoff dependence of the Casimir stress, a circumstance which would seem to deny its physical significance.

To establish the above claims one reverts from the image approach to one based on expansion of the Green’s function in terms of orthogonal functions . To this end one notes that the free field propagator in the radiation gauge can be written as

$`G^{ij}(𝐱𝐱^{},z,z^{},tt^{})={\displaystyle \underset{n\lambda }{}}{\displaystyle \frac{d𝐤d\omega }{(2\pi )^3}e^{i\omega (tt^{})}}`$ (3)
$`\times {\displaystyle \frac{A_{n\lambda }^i(𝐤,z)A_{n\lambda }^j(𝐤,z^{})}{k^2\omega ^2+(n\pi /a)^2iϵ}}e^{i𝐤(𝐱𝐱^{})}`$ (4)
where $`\lambda =1,2`$ refers to the polarization, and spatial coordinates orthogonal to the $`z`$-direction are denoted by a boldface notation. The eigenfunctions $`A_{n\lambda }^i(𝐤,z)`$ satisfy the equation

$$\left[\frac{^2}{z^2}+(n\pi /a)^2\right]A_{n\lambda }^i(𝐤,z)=0$$

and are given explicitly by

$$A_{n1}^i(𝐤,z)=\frac{\overline{k}_i}{|𝐤|}(\frac{2}{a})^{\frac{1}{2}}sin(n\pi z/a)$$
(5)
and

$$A_{n2}^i(𝐤,z)=\frac{1}{|𝐤|\omega _k}\left(\widehat{𝐳}^i\omega _k^2+\widehat{𝐳}^i\right)(\frac{2}{a})^{\frac{1}{2}}cos(n\pi z/a)$$
(6)
where $`\overline{k}_iϵ^{ij}k_j`$ with $`ϵ^{ij}`$ being the usual alternating symbol. It is important to note that each eigenfunction $`A_{n\lambda }^i(𝐤,z)`$ satisfies the boundary conditions $`\widehat{𝐳}\times 𝐄=\widehat{𝐳}𝐁=0`$ at $`z=0,a`$. This means that it is possible to introduce a regularization such that contributions from large values of $`|𝐤|`$ and/or $`n`$ are reduced without destroying the validity of the boundary conditions. This stands in marked contrast with the image method which has no mechanism for the consistent suppression of the contributions of higher order reflections.

In order to determine the regularization appropriate to this calculation one should ideally make reference to the underlying symmetry. Since the latter consists of the reflection $`zaz`$ and the (2+1) dimensional Poincaré group, it is natural to seek to classify regularization schemes according to representations of the latter. The usual cutoff method for this problem invokes a parameter which damps out the large $`\omega _k`$ contributions, an approach which makes no reference to the underlying Lorentz invariance. A far more appropriate technique is to generalize this to a cutoff based on a vector $`\sigma ^\mu `$ in (2+1) dimensions as well as a scalar cutoff $`\mathrm{\Sigma }`$ which can be used to suppress large values of the (2+1) dimensional invariant $`E^2𝐏^2`$ where $`E`$ and $`𝐏`$ are respectively the energy and momentum operators associated with this (2+1) dimensional subspace. Clearly, the credibility of the Casimir effect requires that the result be independent of the relative importance of these two competing cutoffs.

The calculation proceeds by noting that since the limit $`xx^{}`$ is to be taken symmetrically at some point, it is appropriate to use only the imaginary part of the propagator. An appropriately regularized version of this function can be inferred from Eq.(3) to be

$`I`$ $`mG_{\sigma ,\mathrm{\Sigma }}^{ij}(x,x^{})=\pi {\displaystyle \underset{n\lambda }{}}{\displaystyle \frac{d^3k}{(2\pi )^3}\delta (k^2+(n\pi /a)^2)}`$ (8)
$`\times A_{n\lambda }^i(𝐤,z)A_{n\lambda }^j(𝐤,z^{})e^{ik^\mu (xx^{})_\mu }e^{\sigma _\mu k^\mu ϵ(k^0)}e^{\mathrm{\Sigma }(k^2)^{\frac{1}{2}}}`$
where $`ϵ(k^0)`$ is the alternating function and a summation convention convention has been introduced in the Lorentz invariant subspace. Note that since both $`\sigma ^\mu `$ and $`k^\mu `$ are three vectors in that space, they satisfy the orthogonality conditions $`\widehat{z}^\mu \sigma _\mu =\widehat{z}^\mu k_\mu =0`$ . In addition it is clearly necessary to impose $`\overline{\sigma }^2\sigma ^\mu \sigma _\mu >0`$ and $`\sigma ^0>0`$ in order that this propagator exist. It will subsequently be found that its existence also requires $`\mathrm{\Sigma }<\overline{\sigma }`$.

To proceed one uses the regularization (6) and the form of the stress tensor (1). When used in conjunction with the eigenfunctions (4) and (5) the vacuum expectation value of the regularized stress tensor can be determined. With some effort this is found by straightforward calculation to yield the coordinate independent result

$`0`$ $`|`$ $`T^{\mu \nu }|0={\displaystyle \frac{2\pi }{a}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{d^3k}{(2\pi )^3}\delta (k^2+(n\pi /a)^2)e^{\sigma _\mu k^\mu ϵ(k^0)}}`$
$`\times e^{\mathrm{\Sigma }n\pi /a}\left[k^\mu k^\nu +\widehat{z}^\mu \widehat{z}^\nu (n\pi /a)^2\right]`$
which is manifestly both symmetric and traceless. It can be more usefully written as

$`0`$ $`|`$ $`T^{\mu \nu }|0={\displaystyle \frac{1}{a}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}e^{\mathrm{\Sigma }n\pi /a}`$ (10)
$`\times \left({\displaystyle \frac{}{\sigma _\mu }}{\displaystyle \frac{}{\sigma _\nu }}\widehat{z}^\mu \widehat{z}^\nu {\displaystyle \frac{^2}{\sigma ^\alpha \sigma _\alpha }}\right)\mathrm{\Delta }^{n\pi /a}(i\sigma )`$
where $`\mathrm{\Delta }^{n\pi /a}(x)`$ is the (2+1) dimensional function

$$\mathrm{\Delta }^{n\pi /a}(x)=2\pi \frac{d^3k}{(2\pi )^3}e^{ikxϵ(k^0)}\delta (k^2+(n\pi /a)^2)$$

for a particle of mass $`n\pi /a`$. Since this is an $`O(2,1)`$ scalar, $`\mathrm{\Delta }^{n\pi /a}(i\sigma )`$ is a function of only the invariant $`\overline{\sigma }`$ which has the explicit form

$$\mathrm{\Delta }^{n\pi /a}(i\sigma )=\frac{1}{2\pi \overline{\sigma }}e^{\overline{\sigma }n\pi /a}.$$

The insertion of this result into Eq.(7) clearly implies that the sum over $`n`$ exists only for the case that $`\mathrm{\Sigma }<\overline{\sigma }`$ as previously stated. Upon performing the summation over $`n`$ it follows that

$$0|T^{\mu \nu }|0=(\frac{}{\sigma _\mu \sigma _\nu }\widehat{z}^\mu \widehat{z}^\nu \frac{^2}{\sigma ^\alpha \sigma _\alpha })F(\overline{\sigma },\mathrm{\Sigma })$$

where

$$F(\overline{\sigma },\mathrm{\Sigma })=\frac{1}{2\pi a\overline{\sigma }}\frac{1}{1e^{(\mathrm{\Sigma }\overline{\sigma })\pi /a}}.$$

One now performs the usual expansion of the denominator of this expression, discarding terms which give no contribution in the limit of vanishing cutoff, thereby obtaining

$$F(\overline{\sigma },\mathrm{\Sigma })[\frac{1}{2\pi ^2}\frac{1}{\overline{\sigma }}\frac{1}{\overline{\sigma }\mathrm{\Sigma }}+\frac{1}{4\pi a\overline{\sigma }}(1\frac{\mathrm{\Sigma }\pi }{6a})\frac{(\overline{\sigma }\mathrm{\Sigma })^3\pi ^2}{1440\overline{\sigma }a^4}].$$

Upon performing the derivatives and rearranging terms there finally results

$`0`$ $`|`$ $`T^{\mu \nu }|0=[g^{\mu \nu }+3{\displaystyle \frac{\sigma ^\mu \sigma ^\nu }{\overline{\sigma }^2}}\widehat{z}^\mu \widehat{z}^\nu \left]\right\{{\displaystyle \frac{1}{4\pi a\overline{\sigma }^3}}(1{\displaystyle \frac{\mathrm{\Sigma }\pi }{6a}})+`$
$`{\displaystyle \frac{(2\overline{\sigma }\mathrm{\Sigma })(\overline{\sigma }\mathrm{\Sigma })+\frac{2}{3}\overline{\sigma }^2}{2\pi ^2\overline{\sigma }^3(\overline{\sigma }\mathrm{\Sigma })^3}}+{\displaystyle \frac{\pi ^2}{1440a^4}}{\displaystyle \frac{\mathrm{\Sigma }}{\overline{\sigma }}}({\displaystyle \frac{\mathrm{\Sigma }^2}{\overline{\sigma ^2}}}1)\}`$
$`+({\displaystyle \frac{1}{4}}g^{\mu \nu }\widehat{z}^\mu \widehat{z}^\nu )\left[(1{\displaystyle \frac{\mathrm{\Sigma }}{\overline{\sigma }}}){\displaystyle \frac{\pi ^2}{180a^4}}{\displaystyle \frac{4}{3\pi ^2}}{\displaystyle \frac{1}{\overline{\sigma }}}{\displaystyle \frac{1}{(\overline{\sigma }\mathrm{\Sigma })^3}}\right].`$
If (following ) one subtracts the $`a\mathrm{}`$ result, this reduces to the more tractable form

$`0`$ $`|`$ $`\overline{T}^{\mu \nu }|0=\left[g^{\mu \nu }+3{\displaystyle \frac{\sigma ^\mu \sigma ^\nu }{\overline{\sigma }^2}}\widehat{z}^\mu \widehat{z}^\nu \right]`$
$`\times \left\{{\displaystyle \frac{1}{4\pi a\overline{\sigma }^3}}\left(1{\displaystyle \frac{\mathrm{\Sigma }\pi }{6a}}\right)+{\displaystyle \frac{\pi ^2}{1440a^4}}{\displaystyle \frac{\mathrm{\Sigma }}{\overline{\sigma }}}\left({\displaystyle \frac{\mathrm{\Sigma }^2}{\overline{\sigma }^2}}1\right)\right\}`$
$`+\left({\displaystyle \frac{1}{4}}g^{\mu \nu }\widehat{z}^\mu \widehat{z}^\nu \right)\left(1{\displaystyle \frac{\mathrm{\Sigma }}{\overline{\sigma }}}\right){\displaystyle \frac{\pi ^2}{180a^4}}`$
where an overbar notation has been used to denote this subtraction. It is noteworthy that even this removal of the large $`a`$ result does not lead to regularization independent results, a fact which has been remarked upon earlier.

Of particular interest to Casimir calculations are the stress components $`0|\overline{T}^{33}|0`$ and the energy density per unit area $`a0|\overline{T}^{00}|0`$ which are given by

$$0|\overline{T}^{33}|0=\frac{\pi ^2}{240a^4}(1\frac{\mathrm{\Sigma }}{\overline{\sigma }})$$
(11)
and

$``$ $`=`$ $`{\displaystyle \frac{\pi ^2}{720a^3}}\{(1{\displaystyle \frac{\mathrm{\Sigma }}{\overline{\sigma }}}){\displaystyle \frac{3\sigma _0^2\overline{\sigma }^2}{2\overline{\sigma }^2}}`$ (14)
$`\times [{\displaystyle \frac{\mathrm{\Sigma }}{\overline{\sigma }}}({\displaystyle \frac{\mathrm{\Sigma }^2}{\overline{\sigma }^2}}1){\displaystyle \frac{30\mathrm{\Sigma }a^2}{\pi ^2\overline{\sigma }^3}}]\}+`$
$`(a\mathrm{independent}\mathrm{terms})`$
respectively. It is striking that each of these terms retains a significant dependence on the cutoff details. In addition the usual relation assumed (as in ) to hold between $``$ and the stress components, namely

$$0|\overline{T}^{33}|0=\frac{}{a},$$
(15)
is manifestly contradicted by Eqs.(8) and (9) in agreement with results found earlier in the context of the Casimir energy of a sphere . It is significant that the relation (10) asserts an equality of the vacuum stress $`0|\overline{T}^{33}|0`$ which transforms under $`O(2,1)`$ as a scalar while the right hand side transforms as the $`\mu =\nu =0`$ component of a symmetric tensor under this group. Finally, note should be made of the fact that Eq.(9) predicts an additional Casimir force proportional to the divergent indeterminate form $`\mathrm{\Sigma }/a^2\overline{\sigma }^3`$.

To reinforce the conclusions reached here in the case of parallel plates it is useful to consider also the case of the conducting sphere, the only other geometry in three dimensions which has proved amenable to exact calculation . This case was first solved by Boyer  and subsequently verified by a number of authors \[12-15\]. Following reasoning similar to that of the parallel plate case note is made of the fact that the unbroken symmetry in this case consists of time translation and rotational invariance. Thus the natural cutoff parameters in this problem should refer to the energy and angular momentum. The former is the standard one and is well known to give cutoff independent results. It will be the goal here to examine the situation which occurs when a combination of these two is considered.

This is most economically achieved by reference to  which provides a useful separation of the Casimir energy into a finite part and one which requires regularization. Thus one writes for a sphere of radius $`a`$

$$E_c=E_{fin}+E_\sigma $$

where $`E_c`$, $`E_{fin}`$, and $`E_\sigma `$ are respectively the total, the regularization independent part, and the formally divergent parts of the Casimir energy. The quantity $`E_\sigma `$ is given by

$`E_\sigma `$ $`=`$ $`{\displaystyle \frac{1}{4\pi a}}{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}Ree^{i\varphi }{\displaystyle _0^{\mathrm{}}}𝑑yexp(i\nu \sigma ye^{i\varphi })y{\displaystyle \frac{d}{dy}}`$
$`\times (1+y^2e^{2i\varphi })^3`$
where $`\nu =l+\frac{1}{2}`$, $`0<\varphi <\frac{\pi }{2}`$, and $`\sigma `$ is a dimensionless cutoff used to suppress the high frequency modes. Upon choosing a secondary cutoff of the form $`e^{\mathrm{\Sigma }\nu }`$ it is readily found that $`\mathrm{\Delta }E_\sigma `$ (the change induced in $`E_\sigma `$ in the limit of small cutoff) is given by

$$\mathrm{\Delta }E_\sigma =\frac{3\mathrm{\Sigma }}{2\pi a\sigma ^2}_0^{\mathrm{}}𝑑y\frac{y^2}{(1+y^2)^4}\frac{1}{y^2+(\mathrm{\Sigma }^2/\sigma ^2)}.$$

This is evaluated to yield

$$\mathrm{\Delta }E_\sigma =\frac{3}{64a\sigma }\frac{\mathrm{\Sigma }}{(\mathrm{\Sigma }+\sigma )^4}[\mathrm{\Sigma }^2+4\sigma \mathrm{\Sigma }+5\sigma ^2],$$

a result which displays yet again the cutoff dependence of the Casimir effect for a more general choice of regularization. It may be noted that aside from confirming the vanishing of $`\mathrm{\Delta }E_\sigma `$ for $`\mathrm{\Sigma }=0`$, this result shows that $`\mathrm{\Delta }E_\sigma `$ diverges for $`\sigma 0`$ with all intermediate values obtained for finite $`\mathrm{\Sigma }/\sigma `$.

In this work it has been shown that the Casimir effect is, prevailing opinion notwithstanding, highly dependent on the particular form of regularization employed for the extraction of the force. As remarked earlier as well as in ref. the recent experiments which have seemed to many to provide the long awaited precision verification of this highly subtle effect are not based upon rigorous mathematical calculation. While the parallel plate Casimir experiment is fraught with difficulties beyond the ken of this author, it would seem that the successful completion of such experiments would be invaluable for purposes of setting to rest some of the issues which have been raised in this work.

Finally, it would be remiss not to mention in some way the very extensive work on the calculation of Casimir forces using the technique of zeta function regularization . Historically, the successes of the Casimir approach in dealing with the parallel plate geometry and the sphere were obtained using conventional field theoretical subtraction procedures. Specifically, it was noted that only changes relative to the vacuum could be considered observable and it was therefore totally consistent to perform subtractions relative to the $`a\mathrm{}`$ vacuum. However, this step did not succeed in allowing one to obtain finite and observable results in more general applications. Eventually it was realized, however, that the application of zeta function regularization to such problems could yield finite results for some fairly general cases while at the same time agreeing with those obtained in the few instances in which more conventional subtractions could be successfully applied. This work makes no claim to having established any inconsistencies in the derivation of finite results for the Casimir effect when those efforts are based on the twin axioms of vacuum energy and zeta function regularization. Rather, the calculations presented here establish that the Casimir effect is generally cutoff dependent and hence incapable of being reliably determined whenever such calculations are performed using conventional (i.e., physically plausible) subtraction procedures.

###### Acknowledgements.

This work is supported in part by the U.S. Department of Energy Grant No.DE-FG02-91ER40685.