from ColdTOFU.Physics import *
import scipy.linalg as ln
from scipy.constants import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

$^3P_J$ state of $^{87}Sr$ intercombination line

I, L, S = SpinAngularMomenta(I=9/2, L=1, S=1) 
Jx, Jy, Jz = L[0]+S[0], L[1]+S[1], L[2]+S[2]
IdotJ = I[0]*Jx+I[1]*Jy+I[2]*Jz
LdotS = L[0]*S[0]+L[1]*S[1]+L[2]*S[2]
muB = e*hbar/(2*m_e*h*1e6*1e4) # MHz/Gauss

Hamiltonian

Full Hamiltonian including fine structure, hyperfine structure and Zeeman shift reads

$$ H = A_{hfs} \boldsymbol{I}.\boldsymbol{J} + B_{fs} \boldsymbol{L}.\boldsymbol{S} + Q \frac{3\boldsymbol{I}.\boldsymbol{J}(2\boldsymbol{I}.\boldsymbol{J}+1)-2IJ(I+1)(J+1)}{4IJ(2I-1)(2J-1)} + (g_I I_z + g_l L_z + g_s S_z)\mu_B B_{bias}$$

For $^{87}Sr$ inter-combination line, the constants $A_{hfs}$, $B_{fs}$ and $Q$ are $-260.084$ MHz, $5.57 \times 10^6$ MHz and $-35.685$ MHz, respectively. See the following reference for more details:

[1] Martin M. Boyd, Tanya Zelevinsky, Andrew D. Ludlow, Sebastian Blatt, Thomas Zanon-Willette, Seth M. Foreman, and Jun Ye. Nuclear spin effects in optical lattice clocks. Phys. Rev. A, 76:022510, Aug 2007. URL: https://link.aps.org/doi/10.1103/PhysRevA.76.022510,doi:10.1103/PhysRevA.76.022510.

def Hamiltonian(B):
    A_hfs = -260.084
    B_fs = 5.57e6
    Q = -35.685
    H = A_hfs*IdotJ + B_fs*LdotS + Q*(1.5*IdotJ*(2*IdotJ+Id(I[0].shape[0])) -99*Id(I[0].shape[0])/2)/72 + muB*(-1.3e-4*I[2]+L[2]+2*S[2])*B
    return H.toarray()

N = 3000
B = np.linspace(0,1500, N, endpoint=True)
result = np.zeros((I[0].shape[0], len(B)))
for i, bb in enumerate(B):
    result[:, i] = ln.eigh(Hamiltonian(bb), eigvals_only=True)

$^3P_0$

plt.figure()
for i in range(0,10):
    plt.plot(B, result[i,:]-result[0,0], 'C0')
plt.xlabel(r'$B_{bias}$ (Gauss)')
plt.ylabel(r'Shift (MHz)')
plt.tight_layout()

$^3P_1$

plt.figure()
for i in range(10,40):
    plt.plot(B, result[i,:]-result[23,0], 'C0')
plt.xlabel(r'$B_{bias}$ (Gauss)')
plt.ylabel(r'Shift (MHz)')
plt.tight_layout()

$^3P_2$

plt.figure()
for i in range(40,90):
    plt.plot(B[:], result[i,:]-result[68,0], 'C0')
plt.xlabel(r'$B_{bias}$ (Gauss)')
plt.ylabel(r'Shift (MHz)')
plt.tight_layout()