Chapter 5

The Deuteron

The Simplest Nucleus

The deuteron ($^2$H or $d$) is the only bound state of the two-nucleon system and serves as the fundamental laboratory for understanding the nucleon-nucleon interaction. With just one proton and one neutron, it is the simplest composite nucleus and its properties provide direct constraints on the nuclear force.

Experimental Properties

Binding energy: $B_d = 2.2246$ MeV

Spin-parity: $J^\pi = 1^+$

Isospin: $T = 0$

Magnetic moment: $\mu_d = 0.8574\;\mu_N$

Quadrupole moment: $Q_d = 0.2860$ fm$^2$

RMS charge radius: $r_d = 2.1421$ fm

RMS matter radius: $r_m \approx 1.97$ fm

No excited states (only one bound state)

Quantum Numbers and Selection Rules

The deuteron has $J^\pi = 1^+$. The total angular momentum is$\mathbf{J} = \mathbf{L} + \mathbf{S}$, where $\mathbf{S} = \mathbf{s}_p + \mathbf{s}_n$ is the total spin. Since $J = 1$ and the parity is $(-1)^L = +1$, the orbital angular momentum must be even: $L = 0, 2, 4, \ldots$

To couple to $J = 1$ with even $L$, we need the total spin $S = 1$ (triplet state). The allowed states are therefore:

$$^{2S+1}L_J: \quad {}^3S_1 \;(L=0) \quad \text{and} \quad {}^3D_1 \;(L=2)$$

The deuteron is predominantly in the $^3S_1$ state (S-wave, $L=0$), but the tensor force mixes in a $^3D_1$ component. The absence of a bound $^1S_0$ state (singlet, $S=0$) means the nuclear force is spin-dependent: it is stronger in the triplet channel than in the singlet.

The isospin is $T = 0$ (isoscalar). By the generalized Pauli principle for two nucleons, the total wave function must be antisymmetric. With $S = 1$ (symmetric spin) and $T = 0$ (antisymmetric isospin), the spatial wave function must be symmetric, consistent with even $L$.

Square Well Model

The simplest model treats the deuteron as a particle of reduced mass$\mu = m_N/2 \approx 469.5$ MeV/$c^2$ moving in a spherical square well potential of depth $V_0$ and radius $R$:

$$V(r) = \begin{cases} -V_0 & r < R \\ 0 & r > R \end{cases}$$

For the S-wave ($L = 0$), the radial Schrodinger equation with $u(r) = rR(r)$ gives:

$$-\frac{\hbar^2}{2\mu}\frac{d^2 u}{dr^2} + V(r)u = Eu, \qquad E = -B_d$$

The solutions are:

Interior ($r < R$):

$$u(r) = A\sin(kr), \quad k = \sqrt{\frac{2\mu(V_0 - B_d)}{\hbar^2}}$$

Exterior ($r > R$):

$$u(r) = Be^{-\kappa r}, \quad \kappa = \sqrt{\frac{2\mu B_d}{\hbar^2}} = 0.2316\;\text{fm}^{-1}$$

Matching $u$ and $u'$ at $r = R$ gives the eigenvalue condition:

$$k\cot(kR) = -\kappa$$

With $R \approx 2.1$ fm and $B_d = 2.225$ MeV, one finds $V_0 \approx 36$ MeV. The small binding energy ($B_d \ll V_0$) means the deuteron is a weakly bound system with a large spatial extent. The exponential tail $e^{-\kappa r}$ extends far beyond the range of the nuclear force, with a characteristic decay length$1/\kappa \approx 4.3$ fm.

Tensor Force and D-State Admixture

The existence of a nonzero electric quadrupole moment $Q_d = 0.286$ fm$^2$ is direct evidence that the deuteron is not a pure S-state. A state with $L = 0$ has spherical symmetry and therefore $Q = 0$. The nonzero quadrupole moment requires a D-state ($L = 2$) admixture.

The tensor force is responsible for this mixing. It has the operator form:

$$V_T(r) = V_T(r)\,S_{12}, \qquad S_{12} = 3\frac{(\boldsymbol{\sigma}_1 \cdot \hat{r})(\boldsymbol{\sigma}_2 \cdot \hat{r})}{r^2} - \boldsymbol{\sigma}_1 \cdot \boldsymbol{\sigma}_2$$

The tensor operator $S_{12}$ has the same angular structure as the magnetic dipole-dipole interaction and couples states with $\Delta L = 2$, mixing$^3S_1$ and $^3D_1$. In one-pion exchange, the tensor force arises naturally from the pseudoscalar coupling of the pion.

The deuteron wave function is therefore:

$$|\Psi_d\rangle = \frac{u(r)}{r}|{^3S_1}\rangle + \frac{w(r)}{r}|{^3D_1}\rangle$$

where $u(r)$ and $w(r)$ are the S-wave and D-wave radial functions. The D-state probability is:

$$P_D = \int_0^\infty w(r)^2\,dr \approx 4\text{--}7\%$$

Evidence for the Tensor Force

1. Quadrupole moment: $Q_d > 0$ requires $L = 2$ admixture
2. Magnetic moment: Pure S-state predicts $\mu_d = \mu_p + \mu_n = 0.8798\;\mu_N$, but the experimental value is $0.8574\;\mu_N$. The 2.5% discrepancy is explained by D-state mixing
3. Photo-disintegration: The angular distribution of $\gamma + d \to p + n$ requires D-state contributions
4. Asymptotic D/S ratio: $\eta = w(r)/u(r)|_{r\to\infty} = 0.0256 \pm 0.0004$, a precisely measured quantity

Magnetic and Quadrupole Moments

The magnetic moment operator for the deuteron is:

$$\hat{\mu} = \mu_p \hat{\sigma}_p^z + \mu_n \hat{\sigma}_n^z + \frac{\mu_N}{2}\hat{L}_z$$

where $\mu_p = 2.7928\;\mu_N$ and $\mu_n = -1.9130\;\mu_N$. For a mixed S-D state with $M_J = J = 1$:

$$\mu_d = (1 - P_D)(\mu_p + \mu_n) + P_D\left[\frac{3}{4}(\mu_p + \mu_n) - \frac{1}{4}\mu_N\right]$$

Setting $\mu_d = 0.8574\;\mu_N$ gives $P_D \approx 4\%$, but exchange currents (meson exchange corrections to the electromagnetic current operator) modify this estimate to$P_D \approx 5\text{--}7\%$.

The quadrupole moment provides a complementary constraint. For the mixed state:

$$Q_d = \frac{1}{\sqrt{50}} \int_0^\infty w(r)\left[u(r) - \frac{w(r)}{\sqrt{8}}\right] r^2\,dr$$

Effective Range Theory

The low-energy np scattering in the triplet ($^3S_1$) channel is intimately connected to the deuteron bound state through effective range theory. The S-wave phase shift is parameterized as:

$$k\cot\delta_0(k) = -\frac{1}{a_t} + \frac{1}{2}r_t k^2 + \cdots$$

where $a_t = 5.424$ fm is the triplet scattering length and $r_t = 1.759$ fm is the triplet effective range. For a bound state at $k = i\kappa$:

$$-\kappa = -\frac{1}{a_t} - \frac{1}{2}r_t\kappa^2 + \cdots$$

At leading order, $\kappa \approx 1/a_t$, which gives $\kappa \approx 0.184$ fm$^{-1}$, compared to the exact value $\kappa = 0.2316$ fm$^{-1}$. The 20% discrepancy shows that the effective range correction is significant: the deuteron is large but not infinitely large compared to the range of the force.

The fact that $a_t \gg r_t$ (and $a_t \gg R$ where $R$ is the force range) means the deuteron is close to the unitarity limit, a regime where the scattering length diverges and universal properties emerge, independent of the details of the short-range interaction.

Python Simulation: Deuteron Wave Function

Solves for the deuteron bound state in a square well potential via eigenvalue matching, computes the radial wave function, probability density, and magnetic moment with D-state corrections.

Deuteron Bound State

Python

S-wave wave function, probability density, and potential for the np bound state

script.py162 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
# scipy removed - not needed for this simulation

# ──────────────────────────────────────────────
# Deuteron Bound State: Radial Wave Functions
# Square well + D-state admixture
# ──────────────────────────────────────────────

hbar2_over_m = 41.47  # hbar^2/(2*mu) in MeV*fm^2, mu = m_N/2

# Deuteron properties
B_d = 2.2246  # Binding energy in MeV
E = -B_d
kappa = np.sqrt(-E * 2 / (2 * hbar2_over_m))  # decay constant outside well

# Square well parameters (fit to reproduce binding energy)
R = 2.1  # well radius in fm
V0 = 35.0  # well depth in MeV (approximate)

# Solve for V0 that gives correct binding energy
# Inside: u'' + (2*mu/hbar^2)(V0 + E)*u = 0  =>  u'' + k_in^2 * u = 0
# Outside: u'' - kappa^2 * u = 0

def find_V0(V0_trial, R_well, E_bind):
    """Find the well depth that gives correct binding energy"""
    k_in = np.sqrt((V0_trial + E_bind) / hbar2_over_m)
    kappa_out = np.sqrt(-E_bind / hbar2_over_m)
    # Matching condition: k_in * cot(k_in * R) = -kappa
    lhs = k_in / np.tan(k_in * R_well)
    rhs = -kappa_out
    return lhs - rhs

# Bisection to find V0
V0_low, V0_high = 20.0, 80.0
for _ in range(100):
    V0_mid = 0.5 * (V0_low + V0_high)
    if find_V0(V0_mid, R, E) * find_V0(V0_low, R, E) < 0:
        V0_high = V0_mid
    else:
        V0_low = V0_mid
V0 = V0_mid

k_in = np.sqrt((V0 + E) / hbar2_over_m)
kappa_out = np.sqrt(-E / hbar2_over_m)

print(f"Deuteron Bound State Properties")
print(f"{'='*50}")
print(f"Binding energy: B_d = {B_d:.4f} MeV")
print(f"Well radius:    R   = {R:.2f} fm")
print(f"Well depth:     V0  = {V0:.3f} MeV")
print(f"Interior k:     k   = {k_in:.4f} fm^-1")
print(f"Exterior kappa: k   = {kappa_out:.4f} fm^-1")
print()

# Build wave function
r_in = np.linspace(0.001, R, 500)
r_out = np.linspace(R, 15.0, 1000)
r_full = np.concatenate([r_in, r_out])

# u(r) = A*sin(k*r) for r < R
# u(r) = B*exp(-kappa*r) for r > R
A = 1.0
u_in = A * np.sin(k_in * r_in)
B_coeff = A * np.sin(k_in * R) / np.exp(-kappa_out * R)
u_out = B_coeff * np.exp(-kappa_out * r_out)
u_full = np.concatenate([u_in, u_out])

# Normalize: integral of |u(r)|^2 dr = 1
norm = np.trapezoid(u_full**2, r_full)
u_full /= np.sqrt(norm)
u_in = u_full[:len(r_in)]
u_out = u_full[len(r_in):]

# Probability outside well
P_out = np.trapezoid(u_out**2, r_out)
print(f"Probability outside well: {P_out*100:.1f}%")

# Hulth\'en wave function for comparison (analytic)
alpha = kappa_out
beta = 6.0 * alpha  # typical value
N_h = np.sqrt(2 * alpha * beta * (alpha + beta) / (alpha - beta)**2)
r_h = r_full
u_hulthen = N_h * (np.exp(-alpha * r_h) - np.exp(-beta * r_h))
norm_h = np.trapezoid(u_hulthen**2, r_h)
u_hulthen /= np.sqrt(norm_h)

# RMS matter radius
r2_avg = np.trapezoid(u_full**2 * r_full**2, r_full)
r_rms = np.sqrt(r2_avg)
print(f"RMS matter radius:   r_rms = {r_rms:.3f} fm")
print(f"Experimental value:  r_rms ~ 1.97 fm")

# D-state admixture model (simplified)
# P_D ~ 4-7% from tensor force
P_D = 0.058  # 5.8% D-state probability
P_S = 1.0 - P_D
print(f"\nD-state probability: P_D = {P_D*100:.1f}%")
print(f"S-state probability: P_S = {P_S*100:.1f}%")

# Magnetic moment calculation
mu_p = 2.7928  # nuclear magnetons
mu_n = -1.9130
mu_S = mu_p + mu_n  # S-state only prediction
mu_D = 0.75 * (mu_p + mu_n) - 0.25  # with D-state correction
mu_exp = 0.8574  # experimental
mu_theory = P_S * (mu_p + mu_n) + P_D * (0.75*(mu_p + mu_n) - 0.25)
print(f"\nMagnetic moment analysis:")
print(f"  mu_p + mu_n (S-state only): {mu_S:.4f} mu_N")
print(f"  With D-state ({P_D*100:.1f}%):      {mu_theory:.4f} mu_N")
print(f"  Experimental:               {mu_exp:.4f} mu_N")

fig, axes = plt.subplots(1, 3, figsize=(18, 5.5))
fig.patch.set_facecolor('#0a0a1a')

for ax in axes:
    ax.set_facecolor('#0a0a1a')
    for spine in ax.spines.values():
        spine.set_edgecolor('#334155')

# Left: wave function u(r)
axes[0].plot(r_full, u_full, color='#ef4444', linewidth=2.5, label='Square well $u_S(r)$')
axes[0].plot(r_h, u_hulthen, color='#3b82f6', linewidth=2, linestyle='--', label='Hulthen')
axes[0].axvline(x=R, color='#fbbf24', linewidth=1, linestyle=':', alpha=0.7, label=f'R = {R} fm')
axes[0].set_xlabel('r [fm]', color='#94a3b8', fontsize=12)
axes[0].set_ylabel('u(r) [fm$^{-1/2}$]', color='#94a3b8', fontsize=12)
axes[0].set_title('Deuteron S-wave Function', color='white', fontsize=13)
axes[0].set_xlim(0, 12)
axes[0].legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[0].tick_params(colors='#94a3b8')

# Middle: probability density
axes[1].fill_between(r_full, u_full**2, color='#ef4444', alpha=0.3)
axes[1].plot(r_full, u_full**2, color='#ef4444', linewidth=2, label='$|u(r)|^2$')
axes[1].axvline(x=R, color='#fbbf24', linewidth=1, linestyle=':', alpha=0.7)
axes[1].set_xlabel('r [fm]', color='#94a3b8', fontsize=12)
axes[1].set_ylabel('$|u(r)|^2$ [fm$^{-1}$]', color='#94a3b8', fontsize=12)
axes[1].set_title('Probability Density', color='white', fontsize=13)
axes[1].set_xlim(0, 12)
axes[1].tick_params(colors='#94a3b8')
axes[1].legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')

# Right: potential
r_pot = np.linspace(0.001, 10, 500)
V_pot = np.where(r_pot < R, -V0, 0.0)
axes[2].plot(r_pot, V_pot, color='#22c55e', linewidth=2.5, label='Square well')
axes[2].axhline(y=E, color='#ef4444', linewidth=2, linestyle='--', label=f'E = {E:.3f} MeV')
axes[2].axhline(y=0, color='#64748b', linewidth=0.5, linestyle='-')
axes[2].set_xlabel('r [fm]', color='#94a3b8', fontsize=12)
axes[2].set_ylabel('V(r) [MeV]', color='#94a3b8', fontsize=12)
axes[2].set_title('Deuteron Potential', color='white', fontsize=13)
axes[2].set_xlim(0, 8)
axes[2].set_ylim(-V0 - 5, 10)
axes[2].legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[2].tick_params(colors='#94a3b8')

plt.tight_layout()
plt.savefig('deuteron_wavefunction.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation: Deuteron Properties

Solves the eigenvalue condition for the square well, computes wave function properties, and performs a D-state admixture scan for the magnetic moment.

Deuteron Bound State Solver

Fortran

Finds well depth, wave function, and scans D-state effects on magnetic moment

deuteron_solver.f90127 lines

program deuteron_solver
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), parameter :: pi = 3.14159265358979323846_dp
  real(dp), parameter :: hbar2_2mu = 41.47_dp   ! hbar^2/(2*mu) MeV*fm^2
  real(dp), parameter :: B_exp = 2.2246_dp       ! Binding energy MeV

real(dp) :: V0, R_well, k_in, kappa, E_bind
  real(dp) :: V0_low, V0_high, V0_mid, f_low, f_mid
  real(dp) :: mu_p, mu_n, mu_S, mu_D, P_D, mu_theory
  real(dp) :: Q_exp, eta_D
  integer :: i, n_iter
  integer, parameter :: nr = 5000
  real(dp) :: r(nr), u(nr), dr, norm, r2_avg, r_rms
  real(dp) :: P_outside

R_well = 2.1_dp   ! fm
  E_bind = -B_exp

! Find V0 by bisection: matching condition k*cot(k*R) = -kappa
  V0_low = 20.0_dp
  V0_high = 80.0_dp

do n_iter = 1, 200
    V0_mid = 0.5_dp * (V0_low + V0_high)
    f_low = match_func(V0_low, R_well, E_bind)
    f_mid = match_func(V0_mid, R_well, E_bind)
    if (f_low * f_mid < 0.0_dp) then
      V0_high = V0_mid
    else
      V0_low = V0_mid
    end if
  end do
  V0 = V0_mid

k_in = sqrt((V0 + E_bind) / hbar2_2mu)
  kappa = sqrt(-E_bind / hbar2_2mu)

write(*,'(A)') '=== Deuteron Bound State Solver ==='
  write(*,'(A)') ''
  write(*,'(A,F10.4,A)') 'Binding energy:  B_d = ', B_exp, ' MeV'
  write(*,'(A,F10.4,A)') 'Well radius:     R   = ', R_well, ' fm'
  write(*,'(A,F10.4,A)') 'Well depth:      V0  = ', V0, ' MeV'
  write(*,'(A,F10.4,A)') 'Interior k:      k   = ', k_in, ' fm^-1'
  write(*,'(A,F10.4,A)') 'Exterior kappa:  k   = ', kappa, ' fm^-1'

! Build wave function
  dr = 20.0_dp / real(nr, dp)
  norm = 0.0_dp
  do i = 1, nr
    r(i) = (real(i, dp) - 0.5_dp) * dr
    if (r(i) < R_well) then
      u(i) = sin(k_in * r(i))
    else
      u(i) = sin(k_in * R_well) * exp(-kappa * (r(i) - R_well))
    end if
    norm = norm + u(i)**2 * dr
  end do

! Normalize
  do i = 1, nr
    u(i) = u(i) / sqrt(norm)
  end do

! RMS radius and probability outside
  r2_avg = 0.0_dp
  P_outside = 0.0_dp
  do i = 1, nr
    r2_avg = r2_avg + u(i)**2 * r(i)**2 * dr
    if (r(i) > R_well) then
      P_outside = P_outside + u(i)**2 * dr
    end if
  end do
  r_rms = sqrt(r2_avg)

write(*,'(A)') ''
  write(*,'(A,F10.3,A)') 'RMS matter radius:       ', r_rms, ' fm'
  write(*,'(A,F10.1,A)') 'Prob. outside well:      ', P_outside*100.0_dp, ' %'

! Magnetic moment analysis
  mu_p = 2.7928_dp
  mu_n = -1.9130_dp
  mu_S = mu_p + mu_n

write(*,'(A)') ''
  write(*,'(A)') '--- Magnetic Moment Analysis ---'
  write(*,'(A,F10.4,A)') 'Pure S-state:  mu = mu_p + mu_n = ', mu_S, ' mu_N'
  write(*,'(A,F10.4,A)') 'Experimental:  mu = ', 0.8574_dp, ' mu_N'

! D-state admixture scan
  write(*,'(A)') ''
  write(*,'(A)') 'D-state admixture scan:'
  write(*,'(A)') '  P_D (%)   mu_theory    mu_exp     diff'
  write(*,'(A)') '  -------   ---------    ------     ----'
  do i = 0, 10
    P_D = real(i, dp) / 100.0_dp
    mu_theory = (1.0_dp - P_D) * mu_S + P_D * (0.75_dp * mu_S - 0.25_dp)
    write(*,'(F8.1,F12.4,F12.4,F10.4)') P_D*100.0_dp, mu_theory, 0.8574_dp, &
         mu_theory - 0.8574_dp
  end do

! Quadrupole moment
  Q_exp = 0.2859_dp  ! fm^2
  write(*,'(A)') ''
  write(*,'(A,F10.4,A)') 'Experimental quadrupole moment: Q_d = ', Q_exp, ' fm^2'
  write(*,'(A)') 'A nonzero Q requires D-state admixture (tensor force).'
  write(*,'(A)') 'Pure S-state (l=0) has Q = 0 by spherical symmetry.'

! Effective range parameters
  write(*,'(A)') ''
  write(*,'(A)') '--- Effective Range Parameters (triplet channel) ---'
  write(*,'(A)') 'Scattering length:  a_t =  5.424 fm'
  write(*,'(A)') 'Effective range:    r_t =  1.759 fm'
  write(*,'(A,F10.4,A)') 'Bound state:  kappa = ', kappa, ' fm^-1'
  write(*,'(A,F10.4,A)') '1/a_t = ', 1.0_dp/5.424_dp, ' fm^-1 (should ~ kappa)'

contains

function match_func(V0_val, R_val, E_val) result(f)
    real(dp), intent(in) :: V0_val, R_val, E_val
    real(dp) :: f, k_loc, kappa_loc
    k_loc = sqrt((V0_val + E_val) / hbar2_2mu)
    kappa_loc = sqrt(-E_val / hbar2_2mu)
    f = k_loc / tan(k_loc * R_val) + kappa_loc
  end function match_func

end program deuteron_solver

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Nuclear Sizes & Densities Nucleon-Nucleon Scattering →