Nuclear Physics Code

This repository contains Python notebooks covering four core topics in nuclear physics: Binding Energy, Drip Lines, Shell Model, and Spherical Harmonics. Each section includes code and plots illustrating the relevant formulas and concepts.

1. Binding Energy

Calculates the binding energy per nucleon using the semi-empirical mass formula (Weizsäcker formula). For a nucleus with mass number A = N + Z (where N = number of neutrons, Z = number of protons), the binding energy B(A, Z) is computed as:

• Volume term: E_vol = a_v · A     a_v = 15.5 MeV

• Surface term: E_surf = a_s · A^(2/3)     a_s = 16.8 MeV

• Coulomb term: E_coul = a_c · Z · (Z – 1) / A^(1/3)     a_c = 0.72 MeV

• Asymmetry term: E_asym = a_asym · (A – 2Z)^2 / A     a_asym = 23 MeV

• Pairing term: δ(A, Z) = a_p / A^(3/4)     a_p = +34 MeV for even-even nuclei     a_p = 0 for odd-A nuclei     a_p = –34 MeV for odd-odd nuclei

The total binding energy: B(A, Z) = E_vol – E_surf – E_coul – E_asym + δ(A, Z)

The code constructs lists for A (from 1 to some upper limit), computes each term, and then plots all five curves versus A. Finally, it plots the binding energy per nucleon: BE/A = B(A, Z) / A

What the code does

1. Defines constants a_v, a_s, a_c, a_asym.
2. Loops over A and chooses Z optimally (Z ≈ A/2).
3. Calculates E_vol = a_v · A.
4. Calculates E_surf = a_s · A^(2/3).
5. Calculates E_coul = a_c · Z(Z – 1)/A^(1/3).
6. Calculates E_asym = a_asym · (A – 2Z)^2 / A.
7. Determines a_p based on parity, then E_pair = a_p / A^(3/4).
8. Computes B(A, Z) and divides by A to get BE/A.
9. Plots each individual term and BE/A versus A.

2. Drip Lines

Determines the proton and neutron drip lines and plots the β-stability curve in the (N, Z) plane, using a modified semi-empirical mass formula. Drip lines are defined by separation energies:

• Proton separation: S_p(N, Z) = B(N, Z) – B(N, Z – 1) → Proton drip if S_p ≤ 0

• Neutron separation: S_n(N, Z) = B(N, Z) – B(N – 1, Z) → Neutron drip if S_n ≤ 0

Coefficients:

• a_v = 14.1 MeV
• a_s = 13.0 MeV
• a_c = 0.595 MeV
• a_asym = 19 MeV
• a_p as in Section 1

Steps:

1. Implement B(N, Z) with these coefficients.

2. For each N, Z up to A = 300, compute S_p and mark proton drip where S_p = 0.

3. Compute S_n and mark neutron drip where S_n = 0.

4. Compute β-stability: Z_stab(A) ≈ A / [2 + (a_c/a_asym) · A^(2/3)]

5. Plot:

   • Proton drip line (green)
   • Neutron drip line (brown)
   • β-stability curve (blue dots)
   • N = Z line (black)

3. Shell Model

Plots energy levels for a nuclear shell model with harmonic oscillator potential, an ℓ(ℓ + 1) term, and spin-orbit coupling. Formulas:

1. E_ho = (N + 3/2) ℓ ω, N = 2n_r + ℓ
2. E_ll = D · ℓ(ℓ + 1), D = – 0.0225 ℓ ω
3. ℓ·s = [j(j + 1) – ℓ(ℓ + 1) – 3/4]/2, C = – 0.1 ℓ ω
4. E_{ls} = C · ℓ·s

Total: E_total = E_ho + E_ll + E_{ls}

Code:

• Magic numbers: [2, 8, 20, 50, 82, 126]
• ℓ ω = 1 unit
• Loops over N = 0..6, computes energies for all (ℓ, j) combos
• Labels each level (e.g., 1s1/2, 1p3/2)
• Shows energy splitting, indicates magic numbers

4. Spherical Harmonics

Generates and plots Y_l^m(θ, φ) using SciPy. Formula:

Y_l^m(θ, φ) = sqrt[(2l + 1)/(4π) · (l – m)!/(l + m)!] · P_l^m(cos θ) · exp(i m φ)

Steps:

1. Meshgrid: θ in [0, π], φ in [0, 2π]
2. Compute Y_l^m(θ, φ) with sph_harm(m, l, φ, θ)
3. Surface radius R = |Y_l^m| or Re[Y_l^m]
4. Plot 3D surface: x = R sin θ cos φ, y = R sin θ sin φ, z = R cos θ
5. Arrange plots: (2l + 1) per row, l = 0..l_max, shared colorbar

---

How to run each section

1. Clone this repo.
2. Install: pip install numpy matplotlib scipy
3. Open Nuclear_Physics_Code.ipynb in Jupyter.
4. Run cells top to bottom. Each section begins with a header like # Binding Energy, # Drip Lines, etc.

---

Summary of Key Formulas

1. Binding Energy E_vol = a_v A, E_surf = a_s A^(2/3), E_coul = a_c Z(Z – 1)/A^(1/3) E_asym = a_asym (A – 2Z)^2/A, δ(A, Z) = ±34/A^(3/4) (even-even/odd-odd) B(A, Z) = E_vol – E_surf – E_coul – E_asym + δ

2. Drip Lines S_p = B(N, Z) – B(N, Z – 1), S_n = B(N, Z) – B(N – 1, Z) Z_stab ≈ A / [2 + (a_c/a_asym) · A^(2/3)]

3. Shell Model E_ho = (N + 3/2)ℓω, E_ll = Dℓ(ℓ + 1), E_{ls} = Cℓ·s E_total = E_ho + E_ll + E_{ls}

4. Spherical Harmonics Y_l^m = sqrt[(2l + 1)/(4π) (l – m)!/(l + m)!] · P_l^m(cos θ) · e^{i m φ}