1D problem The derivation consists in the one-dimensional case of quantum tunnelling using the
WKB approximation. Considering a wave function of a particle of mass
m, we take area 1 to be where a wave is emitted, area 2 the potential barrier which has height
V and width
l (at 0), and area 3 its other side, where the wave is arriving, partly transmitted and partly reflected. For wave numbers
k [m−1] and energy
E we get: : \Psi_1=Ae^{i(kx+\alpha)}e^{-i{Et}/{\hbar}} : \Psi_2=B_1e^{-k'x}+B_2e^{k'x} : \Psi_3=(C_1e^{-i(kx+\beta)}+C_2e^{i(kx+\beta')})\cdot e^{-i{Et}/{\hbar}} where k = \sqrt{2mE/\hbar^2} and k' = \sqrt{2m(V-E)/\hbar^2}, both in [1/m]. This is solved for given
A and phase
α by taking the boundary conditions at the barrier edges, at x=0 and x=l: there \Psi_{1,3}(t) and its derivatives must be equal on both sides. For k'l \gg 1, this is easily solved by ignoring the time exponential and considering the real part alone (the imaginary part has the same behaviour). We get, up to factors • depending on the
β phases which are typically of order 1, and • of the order of {k}/{k'}=\sqrt{{E}/{(V-E)}} (assumed not very large, since
V is greater than
E (not marginally)): \Psi_1=Ae^{i(kx+\alpha)} , \Psi_3=C_1e^{-i(kx+\beta)}+C_2e^{i(kx+\beta')}, \Psi_2\approx Ae^{-k'x} +Ae^{k'x}: B_1,B_2\approx A and C_1, C_2 \approx \frac{1}{2}A\frac{k'}{k} e^{k'l}. Next, the
alpha decay can be modelled as a symmetric one-dimensional problem, with a standing wave between two symmetric potential barriers at q_0 and -(q_0+l), and emitting waves at both outer sides of the barriers. Solving this can in principle be done by taking the solution of the first problem, translating it by q_0 and gluing it to an identical solution reflected around x=0. Due to the symmetry of the problem, the emitting waves on both sides must have equal amplitudes (
A), but their phases (
α) may be different. This gives a single extra parameter; however, gluing the two solutions at x=0 requires two boundary conditions (for both the wave function and its derivative), so in general there is no solution. In particular, re-writing \Psi_3 (after translation by q_0) as a sum of a cosine and a sine of kx, each having a different factor that depends on
k and
β; the factor of the sine must vanish, so that the solution can be glued symmetrically to its reflection. Since the factor is in general complex (hence its vanishing imposes two constraints, representing the two boundary conditions), this can in general be solved by adding an imaginary part of
k, which gives the extra parameter needed. Thus
E will have an imaginary part as well. The physical meaning of this is that the
standing wave in the middle decays; the waves newly emitted have therefore smaller amplitudes, so that their amplitude decays in time but grows with distance. The
decay constant, denoted
λ [1/s], is assumed small compared to E/\hbar.
λ can be estimated without solving explicitly, by noting its effect on the
probability current conservation law. Since the probability flows from the middle to the sides, we have: : \frac {\partial}{\partial t} \int_{-(q_0+l)}^{(q_0+l)} \Psi^*\Psi\ dx = 2 \frac{\hbar}{2mi}\left(\Psi_1^* \frac{\partial \Psi_1 }{\partial x}- \Psi_1 \frac{\partial \Psi_1^* }{\partial x} \right) , note the factor of 2 is due to having two emitted waves. Taking \Psi\sim e^{-\lambda t}, this gives: : \lambda\frac{2}{4} (q_0+l)\left(A\frac{k'}{k}\right)^2 e^{2k'l}\approx2\frac{\hbar}{m}A^2k. Since the quadratic dependence on k'l is negligible relative to its exponential dependence, we may write: : \lambda\approx4\frac{\hbar k}{m(q_0+l)} \frac{k^2}{k'^2}\cdot e^{-2k'l}. Remembering the imaginary part added to
k is much smaller than the real part, we may now neglect it and get: : \lambda\approx4\frac{\hbar k}{m(q_0+l)}\cdot \frac{E}{V-E}\cdot e^{-2\sqrt{2m(V-E)}l/\hbar}. Note that \frac{\hbar k}{m}=\sqrt{2E/m} is the
particle velocity, so the first factor is the classical rate by which the particle trapped between the barriers (2q_0 apart) hits them.
3D problem Finally, moving to the three-dimensional problem, the spherically symmetric
Schrödinger equation reads (expanding the wave function \psi(r,\theta,\phi) = \chi(r)u(\theta,\phi) in
spherical harmonics and looking at the
l-th term): : \frac{\hbar^2}{2m}\left(\frac{d^2\chi}{dr^2}+\frac{2}{r}\frac{d\chi}{dr}\right)=\left(V(r)+\frac{\hbar^2}{2m}\frac{\ell(\ell+1)}{r^2}-E\right)\chi. Since \ell>0 amounts to enlarging the potential, and therefore substantially reducing the
decay rate (given its exponential dependence on \sqrt{V-E}): we focus on \ell=0, and get a very similar problem to the previous one with \chi(r) = \Psi(r)/r , except that now the potential as a function of
r is not a
step function. In short \frac{\hbar^2}{2m}\left(\ddot\chi+\frac{2}{r}\dot\chi\right)=\left(V(r)-E\right)\chi. The main effect of this on the amplitudes is that we must replace the argument in the exponent, taking an integral of 2\sqrt{2m(V-E)}/\hbar over the distance where V(r)>E rather than multiplying by width
l. We take the
Coulomb potential: : V(r) = \frac {z(Z-z) e^2}{4\pi\varepsilon_0 r} where \varepsilon_0 is the
vacuum electric permittivity,
e the
electron charge,
z = 2 is the
charge number of the alpha particle and
Z the charge number of the nucleus (
Z–
z after emitting the particle). The integration limits are then: r_2 = \frac {z(Z-z) e^2}{4\pi\varepsilon_0 E}, where we assume the nuclear
potential energy is still relatively small, and r_1, which is where the nuclear negative potential energy is large enough so that the overall potential is smaller than
E. Thus, the argument of the exponent in
λ is: : 2\frac {\sqrt{2mE}}{\hbar}\int_{r_1}^{r_2}\sqrt{\frac{V(r)}{E}-1} \,dr=2\frac{\sqrt{2mE}}{\hbar}\int_{r_1}^{r_2}\sqrt{\frac{r_2}{r}-1}\,dr. This can be solved by substituting t=\sqrt{r/r_2} and then t=\cos(\theta) and solving for θ, giving: : 2r_2\frac{\sqrt{2mE}}{\hbar}[\cos^{-1}(\sqrt{x})-\sqrt{x}\sqrt{1-x}]=2\frac{\sqrt{2m}z(Z-z)e^2}{4\pi\varepsilon_0\hbar\sqrt{E}}\left[\cos^{-1}(\sqrt{x})-\sqrt{x}\sqrt{1-x}\right] where x = r_1/r_2. Since
x is small, the
x-dependent factor is of the order 1. Assuming x\ll 1, the
x-dependent factor can be replaced by \arccos0=\pi/2, giving:\lambda\approx e^{-\sqrt{{E_{\mathrm{G}}}/{E}}} with E_{\mathrm G}=\frac{\pi^2m/2\left[z(Z-z)e^2\right]^2}{(4\pi\varepsilon_0\hbar)^2}.Which is the same as the formula given in the beginning of the article with Z_\text{a}=z, Z_\text{b}=Z-z and the fine-structure constant \alpha=\frac{e^2}{4\pi\varepsilon_0\hbar c}: \sqrt{E_{\rm G}}=\sqrt{m/2}/(4\epsilon_0\hbar)[Z_aeZ_be]. For a
radium alpha decay,
Z = 88,
z = 2 and
m ≈ 4
mp,
EG is approximately 50
GeV. Gamow calculated the slope of \log(\lambda) with respect to
E at an energy of 5
MeV to be ~ 1014 J−1, compared to the experimental value of .{{efn|1=\log\lambda\approx-100; {\operatorname{d}\log\lambda\over\operatorname{d}\!E}=\tfrac{1}{2}E_{\rm G}^{1/2}/E^{3/2} (Natural log.) Around 10 /MeV; 'kT = 0.217 fJ = 0.135 keV', 'typical core temperatures in main-sequence stars (the Sun) give kT of the order of 1 keV': 2.17\times10^{-16} joule}} == History ==