Finding the matrix representation of a superoperator I am trying to express superoperator (e.g. the Liouvillian) as matrices and am having a hard time finding a way to do this.
For instance, given the Pauli matrix $\sigma_y$, how do I find the matrix elements of the commutator superoperator? So far I have been trying to figure this out by trial and error (making sure the superoperator acting on the operator vector still gives $[\sigma_y, \rho]$). In the end I want to find superoperators in larger bases, so I am looking for a systematic method to find the matrix elements.
 A: This is exactly analogous to the procedure for finding matrix elements of normal operators. Let's first recall how this works in the familiar case. You choose an orthonormal basis of vectors, say $|n\rangle$, with $n = 1,2,\ldots D$, where $D$ is the dimension of the Hilbert space, such that $\langle n\rvert m\rangle = \delta_{mn}$. Now the matrix elements of an operator $A$ are given by 
$$A_{mn} = \langle m\rvert A\lvert n\rangle.$$
The procedure for superoperators is the same, but the inner product is different. Here it is convenient to use the Hilbert-Schmidt product of two operators:
$$(A,B) = \mathrm{Tr}\{A^{\dagger}B\}.$$
You now must find a complete orthonormal basis with respect to this Hilbert-Schmidt product, i.e. a set of matrices $M_\mu$, with $\mu = 1,2,\ldots,D^2$, such that $(M_\mu,M_\nu) = \delta_{\mu\nu}$. A convenient choice for $D=2$ is the Pauli basis: $$M_\mu \in \left\lbrace\frac{1}{\sqrt{2}}\mathbb{1},\frac{1}{\sqrt{2}}\sigma^x,\frac{1}{\sqrt{2}}\sigma^y,\frac{1}{\sqrt{2}}\sigma^z\right\rbrace.$$
Another simple choice of basis that is easy to generalise is the set of $D^2$ matrices which have one element with value $1$, and all other elements are $0$.
Now if you have a superoperator $\mathcal{L}$, you find its matrix elements through the formula
$$\mathcal{L}_{\mu\nu} = (M_\mu, \mathcal{L}[M_\nu] ).$$
For example, if you have a Hamiltonian $H$ generating a Liouvillian $\mathcal{L}[\bullet] = -i[H,\bullet]$, one of its matrix elements in the Pauli basis would be found from
$$ \mathcal{L}_{xy} = \frac{-i}{2}\mathrm{Tr}\left\lbrace\sigma^x [H,\sigma^y]\right\rbrace. $$
A: If you want to write a super-operator representing left- or right-multiplication, there is a distinct method which is simpler and more elegant. Let us define the left-multiplication superoperator by
$$ \mathcal{L}(A)[\rho] = A\rho,$$
and the right-multiplication superoperator by
$$  \mathcal{R}(A)[\rho] = \rho A.$$
It should be clear that these operations commute, i.e. $\mathcal{L}(A)\mathcal{R}(B) = \mathcal{R}(B)\mathcal{L}(A)$. Many common superoperators can be represented as a sum of these elementary components, for example the commutator: $$ [H,\rho] = \mathcal{L}(H)[\rho] - \mathcal{R}(H)[\rho].$$ Actually I believe all superoperators can be represented in terms of these elementary operations, although I have never proven it: it seems rather obvious.
Now, in order to represent these operations as matrices, you need to flatten your target operator into a vector. One way of performing this mapping is the following
$$ \tag{1}\rho = \sum_{i,j} \rho_{ij}\;\lvert i\rangle\langle j\rvert \to \sum_{i,j} \rho_{ij}\;\lvert i\rangle\otimes\lvert j \rangle . $$ 
In this flattened representation we find
$$\mathcal{L}(A)[\rho] =  \sum_{i,j} \rho_{ij}\;A\lvert i\rangle\langle j\rvert \to \sum_{i,j} \rho_{ij}\;(A\lvert i\rangle)\otimes\lvert j \rangle =   \sum_{i,j} \rho_{ij}\;(A\otimes \mathbb{1})\lvert i\rangle\otimes\lvert j \rangle. $$
Therefore the left-multiplication superoperator is represented by the matrix $\mathcal{L}(A)= (A\otimes\mathbb{1})$. Similarly, you should be able to show that $\mathcal{R}(A) = (\mathbb{1}\otimes A^T)$.
Be warned: many standard computer linear algebra packages do not automatically perform the flattening map according to Eq. (1). For example, the MATLAB reshape() function uses a different convention, meaning that these formulas must be adapted. 
