initially published version

README.md

 # Weak Decay Model Hamiltonian
 
+Finds right and left eigenstates of weak decay model hamiltonian within the span of trial neutron states.

weak_decay.py

+# Weak decay of neutrons
+
+import numpy
+
+### whether the considered states must be eigenstates of total isospin
+isospin_eigenstate = True
+
+## Model hamiltonian
+
+### modelled spatial dimensions per quark
+spatial_dimensions = 3
+
+### modelled total dimensions per quark: 2 flavors, 2 spins, and spatial dims
+quark_dimensions = 2*2*spatial_dimensions
+
+### Clebsh--Gordan coefficients for spin one and spin one-half objects
+### spin order of rows and columns is: -1,0,+1 and -1/2,+1/2, respectively
+### note: everywhere else we use [1,0] for +1/2, [0,1] for -1/2
+cg_1_1half = numpy.array(
+  [
+    [  0.,              -numpy.sqrt(2/3) ],
+    [ +numpy.sqrt(1/3), -numpy.sqrt(1/3) ],
+    [ +numpy.sqrt(2/3),  0.  ]
+  ]
+)
+
+### spin angular momentum transition from neutron to proton + W-boson
+spin_Wp_n = numpy.array(
+  [
+    [ cg_1_1half[1,1], cg_1_1half[0,1] ],
+    [ cg_1_1half[2,0], cg_1_1half[1,0] ]
+  ]
+)
+
+
+### weak decay converts flavor-down quark in flavor-up quark.
+### flavor-up quarks can be projected out, as they will have no overlap with the target state of the proton (duu)
+flavor_Wp_n = numpy.array([[0.,1.],[0.,0.]])
+
+# test action of flavor decay on down-quark
+assert((flavor_Wp_n @ [0,1]==[1,0]).all())
+
+### recursion to get n-body unit matrix, where dim is the single-body dimension
+### could also be done with numpy.identity(dim**n)
+def H0n(dim,n):
+  if n > 1:
+    return numpy.kron(numpy.identity(dim),H0n(dim,n-1))
+  else:
+    return numpy.identity(dim)
+
+### apply single-body operator h for i-th object and identiy of all others
+def H1n(h,i,n):
+  if n > 1:
+    if i > 0:
+      return numpy.kron(numpy.identity(len(h)), H1n(h,i-1,n-1))
+    else:
+      return numpy.kron(h, H0n(len(h),n-1))
+  else:
+    return h
+
+### sum application of single-body operator h for all bodies
+def Hn(h,n):
+  return numpy.sum( [H1n(h,i,3) for i in range(3)], axis=0 )
+
+### triple kronecker product
+def k3(f,g,h):
+  return numpy.kron(numpy.kron(f,g),h)
+
+### single-quark weak decay: left kronecker factor on flavor, middle on spin, right on spatial
+one_body_weak_decay = k3(flavor_Wp_n, spin_Wp_n, numpy.identity(spatial_dimensions))
+
+### 3-body quark decay model Hamiltonian: single-quark decay on each quark
+### as stated in Eq.(12) of the article
+print("bulding weak decay model hamiltonian...")
+H_beta = Hn(one_body_weak_decay,3)
+
+
+## Model states of neutron and proton in ground and excited state
+### Canonical basis states of three quarks, each with (flavor,spin,spatial)
+
+###produce i-th canonical basis state, i starts with 0
+def basis_state(i,n):
+  return numpy.concatenate(
+    (
+      numpy.zeros(i, dtype=numpy.float64),
+      numpy.ones(1, dtype=numpy.float64),
+      numpy.zeros(n-i-1, dtype=numpy.float64)
+    )
+  )
+
+assert((basis_state(3,8)==[0.,0.,0.,1.,0.,0.,0.,0.]).all())
+
+### build super index given object indices, compatible with kronecker_product:
+def super_index(indices,dim):
+  if len(indices) > 0:
+    return indices[-1] + dim*super_index(indices[0:-1],dim)
+  else:
+    return 0
+
+### conversly, get object indices from given super index:
+def object_indices(super_index,dim,n):
+  if n > 1:
+    return numpy.concatenate(
+      (object_indices(super_index//dim,dim,n-1), [super_index%dim])
+    )
+  else:
+    return [super_index]
+
+### e.g. super index of first object in state #1, second in state #0, and third in state #0, each object has 8 states
+assert(super_index([1,2,3],8)==83)
+assert((object_indices(super_index([1,2,3],8),8,3)==[1,2,3]).all())
+
+### make canonical basis states from object indices:
+def super_basis_state(indices):
+  return basis_state(
+    super_index(indices,quark_dimensions),quark_dimensions**len(indices)
+  )
+
+
+## Symmetrization and normalization of three-body states
+
+###build transposition (swapping) of first and second quark:
+transposition_12 = numpy.array(
+  [
+    basis_state(
+      super_index(object_indices(I,quark_dimensions,3)[[1,0,2]],quark_dimensions),
+      quark_dimensions**3
+    )
+    for I in range(quark_dimensions**3)
+  ],
+  dtype=numpy.float64
+)
+
+### test transposition_12 with bracket:
+assert(super_basis_state([3,1,2])@transposition_12@super_basis_state([1,3,2])==1.0)
+
+### transposition of second and thrid quark:
+transposition_23 = numpy.array(
+  [
+    basis_state(
+      super_index(object_indices(I,quark_dimensions,3)[[0,2,1]],quark_dimensions),
+      quark_dimensions**3
+    )
+    for I in range(quark_dimensions**3)
+  ],
+  dtype=numpy.float64
+)
+
+### test transposition_23
+assert(super_basis_state([3,1,2])@transposition_23@super_basis_state([3,2,1])==1.0)
+
+### finally, transposition of first and third quark:
+transposition_13 = numpy.array(
+  [
+    basis_state(
+      super_index(object_indices(I,quark_dimensions,3)[[2,1,0]],quark_dimensions),
+      quark_dimensions**3
+    )
+    for I in range(quark_dimensions**3)
+  ],
+  dtype=numpy.float64
+)
+
+assert(super_basis_state([3,2,1])@transposition_13@super_basis_state([1,2,3])==1.0)
+
+### full symmetrization can be built from identity and combinations of single transpositions:
+symmetrization = (
+  numpy.identity(quark_dimensions**3)
+  + transposition_12 + transposition_23 + transposition_13
+  + transposition_23@transposition_12 + transposition_12@transposition_23
+)
+
+### test symmetrization:
+assert(super_basis_state([1,2,3])@symmetrization@super_basis_state([1,2,3])==1.)
+assert(super_basis_state([1,3,2])@symmetrization@super_basis_state([1,2,3])==1.)
+assert(super_basis_state([2,1,3])@symmetrization@super_basis_state([1,2,3])==1.)
+assert(super_basis_state([2,3,1])@symmetrization@super_basis_state([1,2,3])==1.)
+assert(super_basis_state([3,1,2])@symmetrization@super_basis_state([1,2,3])==1.)
+assert(super_basis_state([3,2,1])@symmetrization@super_basis_state([1,2,3])==1.)
+assert(super_basis_state([1,2,2])@symmetrization@super_basis_state([1,2,2])==2.)
+assert(super_basis_state([2,1,2])@symmetrization@super_basis_state([1,2,2])==2.)
+assert(super_basis_state([2,2,1])@symmetrization@super_basis_state([1,2,2])==2.)
+assert(super_basis_state([1,1,1])@symmetrization@super_basis_state([1,1,1])==6.)
+
+### normalization of state psi
+def normalize(psi):
+  square_norm = psi @ psi
+  return psi / numpy.sqrt(square_norm) if square_norm > 1e-14 else psi
+
+
+## Neutron and proton states
+def tensor_product(As,Bs):
+  return numpy.array(
+    [
+      a*b for a in As for b in Bs
+    ]
+  )
+
+### triple tensor product
+def t3(a,b,c):
+  return tensor_product(tensor_product(a,b),c)
+
+### test compatibility with kronecker product
+assert(
+  (
+    numpy.kron(
+      numpy.identity(2), numpy.array([[0,1],[1,0]])
+    ) @ tensor_product([1,-1],[2,5]) == [5, 2, -5 ,-2]
+  ).all()
+)
+
+### flip flavor of a single quark
+one_body_flavor_flip = k3(
+  numpy.array([[0.,1.],[1.,0.]]), numpy.identity(2), numpy.identity(spatial_dimensions)
+)
+
+### flip flavor in all 3 quarks
+flavor_flip = k3(one_body_flavor_flip,one_body_flavor_flip,one_body_flavor_flip)
+
+### flip spin pf a single quark
+one_body_spin_flip = k3(
+  numpy.identity(2), numpy.array([[0.,1.],[1.,0.]]), numpy.identity(spatial_dimensions)
+)
+
+### flip spin in all 3 quarks
+spin_flip = k3(one_body_spin_flip,one_body_spin_flip,one_body_spin_flip)
+
+### returns the proton state from a given neutron state with spin up
+### that satisfies spin angular momentum conservation
+def proton(psi):
+  return flavor_flip @ (cg_1_1half[1,1]*psi + cg_1_1half[2,0]*spin_flip @ psi)
+
+### reorders the tensor product of (all flavors) x (all spins) x (all orbitals)
+### to (flavor1,spin1,orbital1) x ... x (flavor3,spin3,orbital3)
+def state(f,s,o):
+  return [
+    f[
+      4*f1+2*f2+f3
+    ] * s[
+      4*s1+2*s2+s3
+    ] * o[
+      spatial_dimensions**2*o1+spatial_dimensions*o2+o3
+    ]
+    for f1 in range(2)
+    for s1 in range(2)
+    for o1 in range(spatial_dimensions)
+    for f2 in range(2)
+    for s2 in range(2)
+    for o2 in range(spatial_dimensions)
+    for f3 in range(2)
+    for s3 in range(2)
+    for o3 in range(spatial_dimensions)
+  ]
+
+
+### convenience definitions
+### flavors
+U = [1.,0.]
+D = [0.,1.]
+### spins
+u = [1.,0.]
+d = [0.,1.]
+### ground and excited
+I = [1.,0.,0.]
+J = [0.,1.,0.]
+K = [0.,0.,1.]
+
+### test compatibility with basis states
+#assert(
+#  (state(t3(U,D,D),t3(u,d,d),t3(Px,S,S)) == super_basis_state([,6,6])).all()
+#)
+
+### test for total spin and isospin isostates
+
+print("bulding total isospin and spin operators...")
+sigma_x = numpy.array([[0.,1.],[1.,0.]], dtype=numpy.complex128)
+sigma_y = numpy.array([[0.,-1.j],[1.j,0.]], dtype=numpy.complex128)
+sigma_z = numpy.array([[1.,0.],[0.,-1.]], dtype=numpy.complex128)
+
+def total_general_spin_squared(Sx,Sy,Sz):
+  return Hn(Sx,3) @ Hn(Sx,3) + Hn(Sy,3) @ Hn(Sy,3) + Hn(Sz,3) @ Hn(Sz,3)
+
+### total spin squared
+S2 = total_general_spin_squared(
+  k3(numpy.identity(2), sigma_x, numpy.identity(spatial_dimensions)),
+  k3(numpy.identity(2), sigma_y, numpy.identity(spatial_dimensions)),
+  k3(numpy.identity(2), sigma_z, numpy.identity(spatial_dimensions))
+)
+
+### total isospin squared
+I2 = total_general_spin_squared(
+  k3(sigma_x, numpy.identity(2), numpy.identity(spatial_dimensions)),
+  k3(sigma_y, numpy.identity(2), numpy.identity(spatial_dimensions)),
+  k3(sigma_z, numpy.identity(2), numpy.identity(spatial_dimensions))
+)
+
+def cos_angle(phi,psi):
+  return phi@psi / numpy.sqrt(phi@phi * psi@psi)
+
+def test_spin(SI2, psi):
+  assert(abs(psi @ SI2 @ psi - 3.) < 1e-14)
+  assert(abs(cos_angle(psi, SI2 @ psi) - 1.) < 1e-14)
+
+
+### neutron ground states as stated in Eq.(1)
+psi_n_g = normalize(
+  symmetrization@(
+    state((2*t3(U,D,D)-t3(D,D,U)-t3(D,U,D)),(2*t3(d,u,u)-t3(u,u,d)-t3(u,d,u)),t3(I,I,I))
+  )
+)
+
+### verify that it is equivalent to Eq.(13)
+#assert(
+#  psi_n_g @ normalize(
+#    symmetrization@(
+#        2*super_basis_state([2,4,4])-(super_basis_state([0,6,4])+super_basis_state([0,4,6]))
+#    )
+#  ) == 1.
+#)
+
+### neutron excited states as stated in Eq.(14) and (15)
+psi_n_e_m1 = normalize(
+  symmetrization@(
+    state(
+      2*t3(D,D,U)-t3(D,U,D)-t3(U,D,D),
+      (2*t3(d,u,u)-t3(u,u,d)-t3(u,d,u)),
+      t3(I,J,K)-t3(I,K,J)-t3(J,I,K)+t3(J,K,I)+t3(K,I,J)-t3(K,J,I)
+    )
+  )
+)
+psi_n_e_m2 = normalize(
+  symmetrization@(
+    state(
+      2*t3(D,D,U)-t3(D,U,D)-t3(U,D,D),
+      t3(u,u,d)-t3(u,d,u),
+      2*t3(I,J,K)-t3(J,K,I)-t3(K,I,J)
+      -2*t3(I,K,J)+t3(J,I,K)+t3(K,J,I)
+    )
+  )
+)
+
+psi_n_e_p = normalize(
+  symmetrization@(
+    state(
+      2*t3(D,D,U)-t3(D,U,D)-t3(U,D,D),
+#      t3(u,u,d)-t3(u,d,u),
+      (2*t3(d,u,u)-t3(u,u,d)-t3(u,d,u)),
+      t3(I,J,K)+t3(I,K,J)+t3(J,I,K)+t3(J,K,I)+t3(K,I,J)+t3(K,J,I)
+    )
+  )
+)
+
+### systematic search for left/right eigenvectors of H_beta
+
+if isospin_eigenstate:
+  possible_flavor_states = [
+    (2*t3(U,D,D)-t3(D,D,U)-t3(D,U,D)),
+    (2*t3(D,D,U)-t3(D,U,D)-t3(U,D,D)),
+    (2*t3(D,U,D)-t3(U,D,D)-t3(D,D,U)),
+    t3(D,D,U)-t3(D,U,D),
+    t3(D,D,U)-t3(U,D,D),
+    t3(U,D,D)-t3(D,U,D)
+  ]
+else:
+  # broken isospin symmetry
+  possible_flavor_states = [
+    t3(U,D,D)
+  ]
+
+### trial vectors
+print("building trial neutron states...")
+if not isospin_eigenstate:
+  print("NOTE: isospin symmetry is turned off")
+
+excited_neutron_states = [
+  normalize(
+    symmetrization@(
+      state(
+        flavor_state,
+        spin_state,
+        spatial_state
+      )
+    )
+  )
+  for flavor_state in possible_flavor_states
+  for spin_state in [
+    (2*t3(d,u,u)-t3(u,u,d)-t3(u,d,u)),
+    (2*t3(u,u,d)-t3(u,d,u)-t3(d,u,u)),
+    (2*t3(u,d,u)-t3(d,u,u)-t3(u,u,d)),
+    t3(u,u,d)-t3(u,d,u),
+    t3(u,u,d)-t3(d,u,u),
+    t3(d,u,u)-t3(u,d,u)
+  ]
+  for spatial_state in [t3(I,J,K),t3(I,K,J),t3(J,I,K),t3(J,K,I),t3(K,I,J),t3(K,J,I)]
+]
+
+### remove zero-vectors
+excited_neutron_states = [
+  psi
+  for psi in excited_neutron_states if psi @ psi > 1e-14
+]
+
+### determine span of vectors from svd
+(
+  excited_neutron_states_left_singular_vectors,
+  excited_neutron_states_singular_values,
+  excited_neutron_states_right_singular_vectors
+) = numpy.linalg.svd(excited_neutron_states)
+
+excited_neutron_states_basis = numpy.array(
+  [
+    value_vector_pair[1]
+    for value_vector_pair in zip(
+      excited_neutron_states_singular_values,
+      excited_neutron_states_right_singular_vectors
+    ) if value_vector_pair[0] > 1e-14
+  ]
+)
+
+excited_proton_states_basis = numpy.array(
+  [
+    proton(excited_neutron_state)
+    for excited_neutron_state in excited_neutron_states_basis
+  ]
+)
+
+print("dimension of span:", len(excited_neutron_states_basis))
+
+### test if vectors are iso sin and spin eigenstates with total (iso)spin 3/4
+if isospin_eigenstate:
+  print("testing for isospin and spin eigenstates...")
+else:
+  print("testing for spin eigenstates...")
+  
+for psi in excited_neutron_states_basis:
+  if isospin_eigenstate:
+    test_spin(I2, psi)
+  test_spin(S2, psi)
+
+### project H_beta onto span of trial vectors
+H_beta_projected = numpy.array(
+  [
+    [
+      excited_proton_states_basis[i] @ H_beta @ excited_neutron_states_basis[j]
+      for j in range(len(excited_neutron_states_basis))
+    ]
+    for i in range(len(excited_proton_states_basis))
+  ]
+)
+
+### test whether it is symmetric
+assert(numpy.linalg.norm(H_beta_projected-H_beta_projected.T) < 1e-14)
+
+
+### find eigenvalues and coefficients of right eigenstates
+(
+  H_beta_eigenvalues, H_beta_right_eigenvectors
+) = numpy.linalg.eigh(H_beta_projected)
+
+print("eigenvalues(H):", H_beta_eigenvalues)
+
+### build right (neutron) eigenstates from coefficients
+excited_neutron_eigen_states = numpy.array(
+  [
+    H_beta_right_eigenvector @ excited_neutron_states_basis
+    for H_beta_right_eigenvector in H_beta_right_eigenvectors.T
+  ]
+)
+
+### left and right eigenvectors have the same basis coefficients but are different vectors
+excited_proton_eigen_states = numpy.array(
+  [
+    H_beta_left_eigenvector @ excited_proton_states_basis
+    for H_beta_left_eigenvector in H_beta_right_eigenvectors.T
+  ]
+)
+
+### verify that proton(state[i]) and state[i] are biorthogonal w.r.t. H_beta
+H_beta_proof = numpy.array(
+  [
+    [
+      proton(excited_neutron_eigen_states[i]) @ H_beta @ excited_neutron_eigen_states[j]
+      for j in range(len(excited_neutron_eigen_states))
+    ]
+    for i in range(len(excited_neutron_eigen_states))
+  ]
+)
+assert(
+  numpy.linalg.norm(H_beta_proof - numpy.diag(H_beta_eigenvalues)) < 1e-14
+)
+
+### test if the left eigenvectors are of the form proton(right eigenvectors)
+p_e_overlap = numpy.array(
+  [
+    [
+      proton(excited_neutron_eigen_states[i]) @ excited_proton_eigen_states[j]
+      for j in range(len(excited_proton_eigen_states))
+    ]
+    for i in range(len(excited_proton_eigen_states))
+  ]
+)
+assert(
+  numpy.linalg.norm(
+    p_e_overlap - numpy.identity(len(excited_neutron_states_basis))
+  ) < 1e-14
+)
+
+### overlap of found eigen states with the p_e^+ and p_e^- states from the article
+p_e_comparison = numpy.array(
+  [
+    [
+      excited_proton_eigen_states[i] @ proton(e_mp)
+      for e_mp in [psi_n_e_m1,psi_n_e_m2,psi_n_e_p]
+    ]
+    for i in range(len(excited_proton_eigen_states))
+  ]
+)
+print("<p_i|p_mp>:")
+print(p_e_comparison)
+