Symbolic manipulation of infinite-dimensional band matrices in Mathematica

There are some symbolic calculations in infinite-dimensional Hilbert spaces that can still be done entirely in Mathematica. Assume one has an operator whose infinite-dimensional matrix in some basis vanishes outside some finite number of off-diagonals, and assume there is a closed expression for the corresponding entries: then, polynomial equations in such matrices can be simplified analytically in Mathematica.

The trick is to work with finite-size matrices (sufficiently large to accomodate the finite band width) with symbolic entries, whose ‘middle’ column represents an arbitrary column of our polynomial equation. Since everything is symbolic and all nonzero entries can be accounted for by finiteness of the band width, this allows for analytic equation solving in Mathematica.

Example: generators of \(C^\infty(S^2)\).

Let’s take the generators of the smooth functions on the sphere as example. The algebra \(C^\infty(S^2)\) can be described by two generators \(a\) and \(b\), where \(b^* = b\), \(ab = ba\), \(a^* b = b a^\) and \(a^ a + b^2 = 1\). To see that this does generate \(C^\infty(S^2)\), write \(2 a = x + i y\), \(2 b = 1 - z\) and check that those must satisfy precisely our constraints.

Now, there are two irreducible representations of the algebra generated by \(a\) and \(b\), labeled by the sign \(s = \pm 1\). We want Mathematica to check that these (given) representations are indeed representations of \(C^\infty(S^2)\) by checking whether the relations still hold.

Let’s take our Hilbert space \(H\) to be the direct sum of the alleged representations. An orthonormal basis is the set of vectors \(\lvert l,m,s\rangle\), where \(l\) is a positive half-integer and \(-l <= m <= l\) is another half-integer. In fact, \(H\) is isomorphic to the space of spinors on the sphere (with the standard spin structure) and the standard Dirac operator acts as \(D \lvert l,m,s\rangle = (l+\frac12) \lvert l,m,-s\rangle\).

Now, crucially, the generators \(a\) and \(b\) in this alleged representation do not increase either \(|l|\) or \(|m|\) by more than 1 when acting on \(\lvert l,m,s\rangle\), that is, \(\langle l’,m’,s’ \rvert a \lvert l,m,s\rangle = 0\) for \(|l’-l| > 1\) or \(|m’-m|> 1\) (or \(s’ != s\)). Moreover, the relations we want to check are of order no more than 2. This means that the corresponding polynomials in \(a,b,a^*\) will be (infinite-dimensional) band matrices \(c\) with \(\langle l’,m’,s’ \rvert c \lvert l,m,s\rangle = 0\) for \(|l’-l|>2\) or \(|m’-m| > 2\) or \(s’ != s\).

Then Mathematica can represent \(a\) and \(b\) as follows.

First, setup the space of \(\lvert l,m,s\rangle\) that will be in play:

(* width of the band matrices *)
maxorder = 2;
(* the set of |l,m,s> in play *)
basis = Flatten[
  Table[{myl, mym, S},
    {myl, L - 2, L + 2}, {mym, M - 2, M + 2}]
  , 1];
(* obtain l, m, s as functions of the entry number *)
l[i_] := basis[[i, 1]];
m[i_] := basis[[i, 2]];
s[i_] := basis[[i, 3]];

Now, construct the arrays representing \(a\) and \(b\):

dims = Length[basis];
(* construct array. input: boolean function that specifies nonzero entries and coefficient function *)
bandarray[supportQ_, coeff_] := SparseArray[
  {i_, j_}/; supportQ[i,j] :> coeff[j],
  {dims, dims}];

(* coefficients for a, b *)
fa0[i_] := 2*s[i]*Sqrt[(l[i]+m[i]+1) (l[i]-m[i])] / (2 l[i] (2 l[i]+2));
faplus[i_] := Sqrt[(l[i]+m[i]+1) (l[i]+m[i]+2)] / (2 l[i]+2);
famin[i_] := -Sqrt[(l[i]-m[i]) (l[i]-m[i]-1)] / (2 l[i]);
fb0[i_] := s[i] (
     (l[i]-m[i]+1) (l[i]+m[i])-(l[i]-m[i]) (l[i]+m[i]+ 1)
     )/(2 l[i] (2 l[i]+2));
fbplus[i_] := -Sqrt[(l[i]-m[i]+1) (l[i]+m[i]+1)] / (2 l[i]+2);
fbmin[i_] := -Sqrt[(l[i]-m[i]) (l[i]+m[i])] / (2 l[i]);

a0 = bandarray[(l[#1]==l[#2] && m[#1]==m[#2] + 1 && s[#1]==s[#2]) &, 
  fa0[#1] &];
aplus = bandarray[(l[#1]==l[#2] + 1 && m[#1]==m[#2] + 1 && s[#1]==s[#2]) &,
  faplus[#1] &];
amin = bandarray[(l[#1]==l[#2] - 1 && m[#1]==m[#2] + 1 && s[#1]==s[#2]) &,
  famin[#1] &];
b0 = bandarray[(l[#1]==l[#2] && m[#1]==m[#2] && s[#1]==s[#2]) &, 
  fb0[#1] &];
bplus = bandarray[(l[#1]==l[#2] + 1 && m[#1]==m[#2] && s[#1]==s[#2]) &, 
  fbplus[#1] &];
bmin = bandarray[(l[#1]==l[#2] - 1 && m[#1]==m[#2] && s[#1]==s[#2]) &, 
  fbmin[#1] &];

a = (a0 + aplus + amin);
b = (b0 + bplus + bmin);

Due to the usage of sparse arrays and the symbolic nature of \(a,b\), the resulting computational object is still very easy to handle, and of course matrix multiplication automatically does the right thing with respect to the composition of linear operators.

In order to use Mathematica’s built-in simplification (the main point of the exercise, after all), we’ll set up some assumptions:

(* assume we don't take (nonzero) square roots of negative numbers *)
mypairs = Flatten[
   Table[
    {(myl+mym+1) (myl-mym),
     (myl+mym+1) (myl+mym+2),
     (myl-mym) (myl-mym-1),
     (myl-mym+1) (myl+mym+1),
     (myl-mym) (myl+mym)} /. {myl -> L+x, mym -> M+y}
    , {x, -maxorder, maxorder}, {y, -maxorder, maxorder}]];
possqrts = Map[# >= 0 &, mypairs];
assmptns = 
  FullSimplify[
   Element[S,Integers] && S^2==1
     && Element[L + 1/2, Integers]
     && Element[M + 1/2, Integers]
     && L > 0 && M + L >= 0 && L - M >= 0
     && Apply[And, possqrts]
   ];

We can recover the \(\lvert l, m, s \rangle\) notation using the following helper:

(* ArrayRules on a column returns {i} -> v[i], ..., so the argument should be {index_} *)
stringform[{index_}] := "|" <>
  ToString[l[index]] <> "," <>
  ToString[m[index]] <> "," <>
  ToString[s[index]] <> "⟩";

Let’s set up the simplification procedure. Given an operator \(A\), we want to obtain (symbolically) the columns \(A \lvert l,m,s\rangle \).

lmsindex = Position[basis, {L,M,S}][[1,1]];

(* turn matrix into dictionary {entry -> value} *)
associationform[mat_] := Association[ ArrayRules[mat] ];
(* get the column A|e_n> *)
roughlmscolumn[mat_] := mat[[All,lmsindex]];

(* simplified column association *)
simpleassoc[mat_] := Map[
    FullSimplify[#,Assumptions->assmptns, TimeConstraint->10]&,
    associationform[mat]
  ]
nonzeropart[assoc_] := Select[#=!=0&][assoc];

(* simplified columns, only the nonzero part *)
simplecolumn[mat_] := nonzeropart[simpleassoc[roughlmscolumn[mat]]];

(* print simplified |l,m,s> column in our basis *)
lmscolumn[mat_] := KeyMap[stringform, simplecolumn[mat]];

Now, let’s implement the generating equations of \(C^\infty(S^2)\).

id = bandarray[(#1 == #2)&, 1&];
abcomm = b.a - b.a;
bminbst = b-ConjugateTranspose[b];
absphere = ConjugateTranspose[a].a + b.b - id;

Test the result:

lmscolumn[abcomm]

lmscolumn[bminbst]

lmscolumn[abspere]

Although abcomm et al, even simplified, are certainly not zero (due to boundary effects), the output is three empty associations. Since the column was selected symbolically and we know maxorder was large enough (i.e. all other coefficients in the column must be zero), we’ve now analytically checked the generating equations in our representations for all \(\lvert l, m, s \rangle\). That is, our \(a,b\) operators on \(H\) are indeed a representation of \(C^\infty(S^2)\). Hooray!