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.


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)$.


  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?

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

    KroneckerProduct[A, B]
  2. 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?

  • $\begingroup$ Regarding question (3), what would you like to do with "the product of the whole space"? Perhaps that can help to better define what you want. $\endgroup$
    – jjc385
    Commented Jun 20, 2017 at 20:07
  • $\begingroup$ I guess I would like first to generate a basis for this product space. If I have a matrix which stores as its coulmns the the the basis of choice (states) of each subspace, How would I multiply the matrices to produce a basis for the product space? How do I then identify the quantum numbers which belong to each state in the product space. Trivially what I wrote applies only to spaces with finite dimension. Can this somehow be applied to spaces of infinite dimension such as multiplying momentum degree of freedom with spin? Can I in general do linear algebra of infinite spaces in MMA? $\endgroup$
    – Yair M
    Commented Jun 21, 2017 at 4:26
  • $\begingroup$ Perhaps you could further explain why the linked Q&A (which I answered) is not sufficient? It seems to me to be basically completely general, but if you cab explain why it doesn't do what you want, I can expand on it in an answer here. $\endgroup$
    – march
    Commented Jun 23, 2017 at 1:08
  • $\begingroup$ Obviously, this only works for finite-dimensional spaces. For infinite-dimensional spaces, we'd have to generate some symbolic algebra, and that would be specific to the problem at hand. For instance, there are a couple of questions here about implementing harmonic oscillator algebras. $\endgroup$
    – march
    Commented Jun 23, 2017 at 1:14
  • $\begingroup$ @march, For instance let us assume that I have a Hamiltonian decrying electron states, with this Hamiltonian have a basis of dimension $k$, and I wish do add spin degrees of freedom. How would I add multiply the two degrees of freedom when they don't share the same dimensionality? $\endgroup$
    – Yair M
    Commented Jun 25, 2017 at 12:22

1 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

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

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]


(* 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$:

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

enter image description here

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

(* {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}]


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

enter image description here

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]
(* {ψ[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];


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

enter image description here

  • $\begingroup$ @YairM. I plan on adding interaction operators, but I have other things to work on right now. This is a little different from the way you're implementing it, but perhaps you can adapt some of these ideas. What I like about this approach is that we can take advantage of the fact that things like KroneckerProduct, Tuples, Flatten, Partition, etc. all respect the same ordering, so we can build operators and states in tensor-product spaces where the ordering of the basis is consistent throughout (which I've attempted to show). $\endgroup$
    – march
    Commented Jun 27, 2017 at 17:39
  • 1
    $\begingroup$ Thanks again. I just wanted to comment that this can also be done in the ordering which I use in the linked file (in the question). I just did it, and personally I prefer it for one reason. The total particle number in this ordering in monotonically increasing along the diagonal (unlike the ordering you provide). As a result every Hamiltonian which one might consider which conserves the overall particle number will be block diagonal according to my ordering. $\endgroup$
    – Yair M
    Commented Jun 28, 2017 at 0:34
  • 1
    $\begingroup$ Just a small suggestion (I upvoted this of course): using ann -> SparseArray[{{0, 1}, {0, 0}}] makes the annihilation operators come out as SpareseArrays which may help conserve memory (I found it very useful when dealing with product spaces here. $\endgroup$
    – Jens
    Commented Jun 28, 2017 at 1:32
  • 1
    $\begingroup$ a note before I forget: the first definition could be made more concise using ToeplitzMatrix @@ (PadRight[{ann}, 2^n, #] & /@ {id, sZ}) $\endgroup$
    – Mr.Wizard
    Commented Jun 28, 2017 at 15:25
  • 1
    $\begingroup$ @YairM.Yeah: the overlaps are already in $U$, so if you have an interaction written in terms of number operators, just use those. As for the basis ordering: it' s relatively straight-forward to re-order the basis after everything is constructed in order to block-diagonalize in particle-number. I've added a section at the end on how to do this. $\endgroup$
    – march
    Commented Jun 28, 2017 at 18:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.