multiplication of vector spaces

Question

Original Question

Disclaimer: The question motivated by the product of Hilbert spaces which arises from separate degrees of freedom which a system might contain in quantum technical system.

Examples are for instance momentum and spin of an electron.

The question is related to this one: How to create simple (tensor) product spaces? but I think the example there is rather specific, and I'm not sure how generalize it, nor do I think that it gives a full answer to my question.

Intro

Assuming I have two independent degrees of freedom (random variables), the Probability space describing each of them is an Hilbert space, with the dimensions of the relevant sample space.

Denoting these spaces, by $V$ and $W$, we may multiply these spaces to get a vector space which incorporates the two: $V\otimes W$ (joint probability space).

We may denote two vectors as $v\in V$, $w\in W$, and the corresponding vector in the product space would be $v\otimes w$.

Lastly we may apply some linear operator to the vector. Denoting $C=A\otimes B$, we'll get: $C(v\otimes w)=(Av)\otimes (Bw)$.

Questions

1. How can I achieve the multiplication of vectors? So far I've done:

   Flatten[c*(KroneckerProduct[v, w] + KroneckerProduct[w, v])]
   

where c is just a normalization factor so that I can preserve the meaning of probability. Is there a cleaner way to do this. Is it generally correct/robust?

2. How can I achieve multiplication of linear operators? will the following do?

   KroneckerProduct[A, B]
   

3. Is there a way to take the product of a whole space? I admit this question isn't well defined, i.e. I'm not sure what would this mean in MMA

Edit according to comments

As suggested in the comments, I narrow down my question as it is too general, and solutions may differ greatly per system.

I consider a model of an elecronic system with two sites, and spin degree of freedom. Since I have two separate degrees of freedom, each admitting to 2 values of there quantum numbers, my Fock space is finite and has a dimension $dim=16$.

What I have so far

I have already implemented a few functions that build for me a Fock space with the dried algebra, for a space of of arbitrary dimension. These can be seen here: notebook. This file implements function that eventually gets as input the numberOfSites in the desired space and returns:

• the number numberOfSites
• the dimension of the space (2^numberOfSites)
• a list of states written as occupation vectors. Each state written as a list (vector), of 0 and 1 according to the occupation of each site in the specific state.
• a list of states in the standard basis, i.e. a list of unit orthogonal unit vectors of the proper dimension, i.e. Identity matrix of the proper dimension.
• a list of creation operators (matrices) obeying the anti-commutation algebra, which may be applied to states in the standard basis representation.

Issues / What I need

Even though I can write the basis for my Hilbert space in both representations, and even though I have I have all ladder operators I still think my solution falls a bit short. My issue is with the fact that since my space is so abstract I cannot attribute values of quantum numbers to my states.

For instance I would like each of states to have as an attribute a definite spin so that I could apply an $S_z$ operator to it and get the correct eigenvalue.

The way I see it there are two options to improve what I have:

1. build the Hilbert space as an Associative Array. With the help of the key-value structure, my states could be more than just vectors, they can also hold other attributes. However, this may prove problematic at a some stage. Details on that at the bottom.
2. Perhaps I could build product spaces, in some sense, and then multiply them. Building spin operators for only spinors for instance is trivial. These could then be multiplied perhaps to obtain the product space version for each operator.

A note about the problem with Associative Array

Let us assume for instance, that i describe a system as discussed above, and build all states with having a definite location (site index), and spin. Let us assume now the I consider a Hamiltonian in which I have hopping between the two sites. Obviously the site index isn't a good quantum number, i.e. if I diagonalize my Hamiltonian the eigenstates wouldn't have a definite site index. As I would want to be able to rotate my states between bases, it means that I'll need to build the mechanism which resolves the the attributes attributed to a state after rotation. For instance in the example I gave, one could use intersection to determine the the site indices are different, and therefore this attribute should be dropped. Can I generalize this?

Answer 1

Here's a take that allows one to keep track of the order of things carefully. Note that this is similar in nature to the answer here.

Annihilation operators

We first construct the annihilation operators as follows. This uses a trick I picked up somewhere during grad school that builds in the anti-commutation rules automatically by taking tensor products with the Pauli matrix $\sigma_z$. That is, if there are $n$ total single-particle states (in this case $n=2^m$, where $m$ is the number of sites), then we construct the annihilation operators as $$ \hat{a}_j\to\hat{I}\otimes\cdots\otimes\hat{I}\otimes\hat{\sigma}_+\otimes\hat\sigma_z\otimes\cdots\otimes\hat\sigma_z, $$ where $$ \hat{\sigma}_z= \left[\begin{array}% 1 & 0 \\ 0 & -1 \end{array}\right] $$ and $$\hat{\sigma}_+= \left[\begin{array}% 0 & 1 \\ 0 & 0 \end{array}\right]$$ and it is in the $j$'th spot out of $n$ spots. Translating this construction literally, we can a list of annihilation operators as

Clear[annihilationOperators]
annihilationOperators[numSites_ /; numSites > 1] := 
  KroneckerProduct @@@ Table[
    Join[Table[id, {n}], {ann}, Table[sZ, {2^numSites - n - 1}]],
    {n, 0, 2^numSites - 1}
   ] /. {id -> IdentityMatrix[2], sZ -> PauliMatrix[3], ann -> {{0, 1}, {0, 0}}}

although Mr. Wizard points out that this can be done more succinctly as

Clear[annihilationOperators]
annihilationOperators[numSites_ /; numSites > 1] := 
  KroneckerProduct @@@ ToeplitzMatrix @@ (PadRight[{ann}, 2^numSites, #] & /@ {id, sZ})
    /. {id -> IdentityMatrix[2], sZ -> PauliMatrix[3], ann -> {{0, 1}, {0, 0}}}

It also might be worth making this a SparseArray by instead using the replacement ann -> SparseArray[{{0, 1}, {0, 0}}]; as mentioned by Jens, this can result in significant speed-up, especially since these are very sparse arrays.

A site-number operator can be formed as

Transpose[#].# &@annihilationOperators[2][[1]] // MatrixForm

and the total number operator can be formed as

Transpose[#].# & /@ annihilationOperators[2] // Total // MatrixForm

We can verify that the anti-commutation relations are satisfied like so:

verifyAnticommutationRelations[numSites_ /; numSites > 1] := And @@ Flatten@Outer[
  Array[0 &, {2^numSites, 2^numSites}] === Simplify[#1.#2 + #2.#1] &,
  #, #, 1
 ] &@ annihilationOperators[numSites]

Then:

verifyAnticommutationRelations[2]
(* True *)

Single-particle operators

Now, we construct the single-particle operators. Here, we carefully think about the order of the basis. I prefer something like the following. First, make a test wave function:

wf = ψ @@@ Tuples[{{"up", "down"}, {1, 2}}]
(* {ψ["up", 1], ψ["up", 2], ψ["down", 1], ψ["down", 2]} *)

This fixes the order of our basis. Now, let's create a single-particle operator, say, $\sigma_z$:

Clear[spinZ]
spinZ[numSites_] := KroneckerProduct[PauliMatrix[3], IdentityMatrix[numSites]]/2; 
spinZ[2] // MatrixForm

This makes it clear what the order of the basis is: note that

spinZ[2].wf
(* {1/2 ψ["up", 1], 1/2 ψ["up", 2], -(1/2) ψ["down", 1], -(1/2) ψ["down", 2]} *)

which matches what it should.

We then construct the corresponding Fock-space operator in the normal way, i.e. using the exact definition, $$ \hat{O}\to\sum_{ij}\hat{a}^{\dagger}_j\langle j|\hat{O}|i\rangle \hat{a}_i. $$ We do:

makeSingleParticleFockOperator[op_, numSites_] /; Length@op == 2 numSites := Sum[
  ConjugateTranspose[annihilationOperators[numSites][[j]]].annihilationOperators[numSites][[i]] op[[j, i]],
  {j, 4}, {i, 4}]

Then:

makeSingleParticleFockOperator[spinZ[2], 2] // MatrixForm

Full basis test

To understand how the Fock basis is ordered, consider the following. Create a test Fock-space wave function as

fock[numSites_] := ψ @@@ Tuples[{0, 1}, 2^numSites]
fock[2]
(* {ψ[0, 0, 0, 0], ψ[0, 0, 0, 1], ψ[0, 0, 1, 0], ψ[0,0, 1, 1],
    ψ[0, 1, 0, 0], ψ[0, 1, 0, 1], ψ[0, 1, 1, 0], ψ[0, 1, 1, 1],
    ψ[1, 0, 0, 0], ψ[1, 0, 0, 1], ψ[1, 0, 1, 0], ψ[1, 0, 1, 1],
    ψ[1, 1, 0, 0], ψ[1, 1, 0, 1], ψ[1, 1, 1, 0], ψ[1, 1, 1, 1]} *)

These are explicitly the occupation numbers of the four states, ordered correctly. That is, ψ[0, 1, 0, 1] means that there are two particles in the system occupying single-particle states ψ["up", 2] and ψ["down", 2], so that there are two particles occupying site 2, one with spin-up and one with spin-down. (Recall the ordering of the single-particle basis: ψ["up", 1], ψ["up", 2], ψ["down", 1], ψ["down", 2].)

To verify that this is indeed the correct order, we can test this by acting on this test wave function with various Fock-space operators. For instance,

makeSingleParticleFockOperator[spinZ[2], 2].fock[2]
(* {0, -(1/2) ψ[0, 0, 0, 1], -(1/2) ψ[0, 0, 1, 0], -ψ[0, 0, 1, 1],
    1/2 ψ[0, 1, 0, 0], 0, 0, -(1/2) ψ[0, 1, 1, 1],
    1/2 ψ[1, 0, 0, 0], 0, 0, -(1/2) ψ[1, 0, 1, 1],
    ψ[1, 1, 0, 0], 1/2 ψ[1, 1, 0, 1], 1/2 ψ[1, 1, 1, 0], 0} *)

By comparing to fock[2] above, we can see that things work out correctly. For instance, the spot ψ[1, 1, 0, 1] gets multiplied by 1/2 because there are three particles, two occupying the up state and one occupying the down state. Or, notice there is a zero where the ψ[0, 1, 0, 1] spot is, because the spins cancel in that case.

We can also make sure that the single-particle states are ordered correctly by considering the number operators. For instance:

ConjugateTranspose[#].#.wf &@annihilationOperators[2][[1]]
(* {0, 0, 0, 0, 0, 0, 0, 0,
    ψ[1, 0, 0, 0], ψ[1, 0, 0, 1], ψ[1, 0, 1, 0], ψ[1, 0, 1, 1],
    ψ[1, 1, 0, 0], ψ[1, 1, 0, 1], ψ[1, 1, 1, 0], ψ[1, 1, 1, 1]} *)

Notice that this is zero when we have ψ[0, ___] and 1 when we have ψ[1, ___], i.e. if state 1 is occupied, we get on particle in state 1. Or,

ConjugateTranspose[#].#.wf &@annihilationOperators[2][[2]]
(* {0, 0, 0, 0,
    ψ[0, 1, 0, 0], ψ[0, 1, 0, 1], ψ[0, 1, 1, 0], ψ[0, 1, 1, 1],
    0, 0, 0, 0,
    ψ[1, 1, 0, 0], ψ[1, 1, 0, 1], ψ[1, 1, 1, 0], ψ[1, 1, 1, 1]} *)

which is the same thing except for state 2.

Interaction operators

Finally, we can build an interaction operator. To do this, we again use the explicit exact form of a Fock-space interaction operator, i.e. $$ \hat{V} \to \sum_{ijkl}\hat{a}^{\dagger}_i\hat{a}^{\dagger}_j \langle ij|\hat{V}|kl\rangle\hat{a}_l\hat{a}_k $$ Remember that ordering is important!

Alternatively, interaction operators can often be written in terms of the number operators, as in the Fermi-Hubbard model: $$ \hat{V} = U\sum_i\hat{n}_{\uparrow,i}\hat{n}_{\downarrow,i} $$ This is easier to implement.

Block-diagonalization in particle-number

Typically, the number of particles is conserved, which means that the Hamiltonian is then block-diagonal in particle-number. In order to take advantage of this, we re-order the Fock-basis as follows:

basisReordering = ConjugateTranspose[#].# & /@ annihilationOperators[2] // Total // Diagonal // Ordering
(* {1, 2, 3, 5, 9, 4, 6, 7, 10, 11, 13, 8, 12, 14, 15, 16} *)

Then, transform the annihilation operators as

anns = #[[basisReordering, basisReordering]] & /@ annihilationOperators[2];

Then:

ConjugateTranspose[#].# & /@ anns // Total // MatrixForm