jagot/HilbertSpaces.jl

HilbertSpaces.jl

A small DSL for the kind of linear algebra that shows up in quantum mechanics. Dirac notation can be approximated as

(i|C|j)

where i is a bra vector, C an operator, j a ket vector. Parentheses (although not strictly necessary, from Julia’s point of view) are used instead of 〈, 〉.

A bra is a dual vector of a ket, or a functional, which when acting on a ket, yields a scalar. An operator on the space spanned by kets (bras) will map a ket (bra) onto a ket (bra).

Ideas/ponderings:

• [ ] For different spaces, different vectors and operators can be defined, which follow the same linear algebra rules, but a spin operator might not act on a angular momentum ket (say).
• [-] Deferred evaluation: the contraction (i|C|j) should return a
  • [X] Scalar object which can subsequently be evaluated, either
  • [X] numerically (using WIGXJPF.jl for angular momentum algebra), or
  • [ ] analytically (where applicable/practical, maybe using SymPy.jl?).
  • [X] How to represent (i|j)? Maybe as (i|I|j), where I is the identity operator.
• [ ] How to represent vectors belonging to coupled Hilbert spaces (e.g. |ℓm⊗χ〉)?
• [ ] The transpose of a column vector of kets should return a row vector of bras.
• [ ] How to represent linear combinations of operators? Useful for angular momentum representation of the Cartesian basis vectors; z/r is simply C₀¹, but x/r = [C₋₁¹ - C₁¹]/√2.

Simple usage example

First we specify we want to consider the Hilbert space spanned by the spherical harmonics:

using HilbertSpaces

s = AngularSpace

We then pick two vectors (one normal, one dual), and a spherical tensor:

i = bra(s)(1, 0)
j = ket(s)(0, 0)
C = op(s)(1, 0)
i, j, C

To calculate the contraction (or matrix element), we simply write

v = (i|C|j)

〈p|C₀⁽¹⁾|s〉 ≈ 0.57735

To get the numerical value, we use

float(v)

0.5773502691896258

which is equal to 1/√3.

Here we instead loop over all angular momenta [0,10) (s,p,d,f… in spectroscopic notation) and calculate the matrix representation of z/r:

N = 10
z = zeros(N,N)
C = op(s)(1,0)
for i in 1:N
    b = bra(s)(i-1, 0)
    for j in 1:N
        k = ket(s)(j-1, 0)
        z[i,j] = round(float(b|C|k), 4)
    end
end
z

0	0.5774	0	0	0	0	0	0	0	0
0.5774	0	0.5164	0	0	0	0	0	0	0
0	0.5164	0	0.5071	0	0	0	0	0	0
0	0	0.5071	0	0.504	0	0	0	0	0
0	0	0	0.504	0	0.5025	0	0	0	0
0	0	0	0	0.5025	0	0.5017	0	0	0
0	0	0	0	0	0.5017	0	0.5013	0	0
0	0	0	0	0	0	0.5013	0	0.501	0
0	0	0	0	0	0	0	0.501	0	0.5008
0	0	0	0	0	0	0	0	0.5008	0