Tensor matrix elements

The matrix element of any irreducible tensor operator of rank $k$ is defined as

\[\begin{equation} \tag{V13.1.1} \matrixel{n'j'm'}{\tensor{T}^{(k)}_q}{njm} \defd \int\diff{\tau} \conj{\Psi}_{n'j'm'} \tensor{T}^{(k)}_q \Psi_{njm}. \end{equation}\]

For simple quantum systems, it is most easily evaluated using the Wigner–Eckart theorem; however, in most cases, the quantum system and/or the tensor is composed of multiple parts, which complicates the situation. The table below shows the most common cases:

Basis \ Tensor	Acts on entire system	Acts on subsystems
Uncoupled	Transform to coupled basis: $\eqref{eqn:coupling}$	Evaluation in subspace
Coupled	Wigner–Eckart: $\eqref{eqn:wigner-eckart}$	Uncoupling: $\eqref{eqn:uncoupling}$

Interface

There are two interfaces provided for computation of matrix element of tensor operators:

A low-level interface matrix_element((γj′, m′), 𝐓ᵏq, (γj′, m′)) (and friends), where γj′, m′, γj′, and m′ are quantum numbers.
A high-level interface dot(a, 𝐓ᵏq, b), where a and b are SpinOrbitals from AtomicLevels.jl. The high-level interface dispatches as appropriate to the low-level interface, depending on whether a and b are expressed in the coupled basis or not, and which part of the quantum system 𝐓ᵏq acts on.

High-level interface

There are two main functions in the high-level interface:

dot(a, 𝐓ᵏq, b) for one-body interactions
dot((a,b), 𝐓ᵏq, (c,d)) for two-body interactions (e.g. Coulomb interaction).

LinearAlgebra.dot — Method

dot(a::SpinOrbital,
    𝐓ᵏq::Union{TensorComponent,TensorScalarProduct},
    b::SpinOrbital)

Compute the matrix element ⟨a|𝐓ᵏq|b⟩ in the basis of spin-orbitals, dispatching to the correct low-level function matrix_element, depending on the value of system(𝐓ᵏq).

Examples with coupled orbitals

julia> a,b,c,d = (SpinOrbital(ro"2p", half(3)),
                  SpinOrbital(ro"2p", half(1)),
                  SpinOrbital(ro"2s", half(1)),
                  SpinOrbital(ro"3d", half(3)))
(2p(3/2), 2p(1/2), 2s(1/2), 3d(3/2))

julia> 𝐋, 𝐒, 𝐉 = OrbitalAngularMomentum(), SpinAngularMomentum(), TotalAngularMomentum()
(𝐋̂⁽¹⁾, 𝐒̂⁽¹⁾, 𝐉̂⁽¹⁾)

julia> 𝐋², 𝐒², 𝐉² = 𝐋⋅𝐋, 𝐒⋅𝐒, 𝐉⋅𝐉
((𝐋̂⁽¹⁾⋅𝐋̂⁽¹⁾), (𝐒̂⁽¹⁾⋅𝐒̂⁽¹⁾), (𝐉̂⁽¹⁾⋅𝐉̂⁽¹⁾))

julia> dot(a, cartesian_tensor_component(𝐉, :x), b),
           1/2*√((3/2+1/2+1)*(3/2-1/2)) # 1/2√((J+M+1)*(J-M))
(0.8660254037844386, 0.8660254037844386)

julia> dot(a, cartesian_tensor_component(𝐉, :z), a), a.m[1]
(1.5, 3/2)

julia> dot(a, TensorComponent(𝐋, 0), a)
0.9999999999999999

julia> dot(c, cartesian_tensor_component(Gradient(), :x), a)
- 0.408248(∂ᵣ + 2/r)

julia> dot(c, cartesian_tensor_component(SphericalTensor(1), :x), a)
-0.40824829046386296

julia> orbitals = rsos"2[p]"
6-element Array{SpinOrbital{RelativisticOrbital{Int64},Tuple{Half{Int64}}},1}:
 2p-(-1/2)
 2p-(1/2)
 2p(-3/2)
 2p(-1/2)
 2p(1/2)
 2p(3/2)

julia> map(o -> dot(o, 𝐉², o), orbitals)
6-element Array{Float64,1}:
 0.7499999999999998
 0.7499999999999998
 3.7500000000000004
 3.7500000000000004
 3.7500000000000004
 3.7500000000000004

julia> 1/2*(1/2+1),3/2*(3/2+1) # J(J+1)
(0.75, 3.75)

julia> dot(a, 𝐋², a), 1*(1+1)
(2.0, 2)

julia> dot(d, 𝐋², d), 2*(2+1)
(5.999999999999999, 6)

julia> dot(a, 𝐒², a), 1/2*(1/2+1)
(0.7499999999999998, 0.75)

julia> dot(a, 𝐋⋅𝐒, a)
0.4999999999999999

Examples with uncoupled orbitals

julia> a,b,c,d = (SpinOrbital(o"2p", 1, half(1)),
                  SpinOrbital(o"2p", -1, half(1)),
                  SpinOrbital(o"2s", 0, half(1)),
                  SpinOrbital(o"3d", 2, -half(1)))
(2p₁α, 2p₋₁α, 2s₀α, 3d₂β)

julia> 𝐋,𝐒,𝐉 = OrbitalAngularMomentum(),SpinAngularMomentum(),TotalAngularMomentum()
(𝐋̂⁽¹⁾, 𝐒̂⁽¹⁾, 𝐉̂⁽¹⁾)

julia> 𝐋²,𝐒²,𝐉² = 𝐋⋅𝐋,𝐒⋅𝐒,𝐉⋅𝐉
((𝐋̂⁽¹⁾⋅𝐋̂⁽¹⁾), (𝐒̂⁽¹⁾⋅𝐒̂⁽¹⁾), (𝐉̂⁽¹⁾⋅𝐉̂⁽¹⁾))

julia> dot(a, cartesian_tensor_component(𝐉, :z), a), sum(a.m)
(1.5, 3/2)

julia> dot(a, TensorComponent(𝐋, 0), a)
0.9999999999999999

julia> # Same as previous, but with spin down
       dot(a, TensorComponent(𝐋, 0), SpinOrbital(o"2p", 1, -half(1)))
0

julia> dot(d, TensorComponent(𝐋, 0), d)
1.9999999999999998

julia> dot(d, TensorComponent(𝐒, 0), d)
-0.49999999999999994

julia> dot(c, cartesian_tensor_component(Gradient(), :x), a)
- 0.408248(∂ᵣ + 2/r)

julia> dot(c, cartesian_tensor_component(SphericalTensor(1), :x), a)
-0.40824829046386285

julia> orbitals = sos"2[p]"
6-element Array{SpinOrbital{Orbital{Int64},Tuple{Int64,Half{Int64}}},1}:
 2p₋₁α
 2p₋₁β
 2p₀α
 2p₀β
 2p₁α
 2p₁β

julia> # Only 2p₋₁β and 2p₁α are pure states, with J = 3/2 => J(J + 1) = 3.75
       map(o -> dot(o, 𝐉², o), orbitals)
6-element Array{Float64,1}:
 1.7499999999999998
 3.7500000000000004
 2.75
 2.75
 3.7500000000000004
 1.7499999999999998

julia> dot(a, 𝐋², a), 1*(1+1)
(2.0, 2)

julia> dot(d, 𝐋², d), 2*(2+1)
(5.999999999999999, 6)

julia> dot(a, 𝐒², a), half(1)*(half(1)+1)
(0.7499999999999998, 0.75)

julia> dot(a, 𝐋⋅𝐒, a)
0.4999999999999999