Two-oscillator system coupled to reservoir I am implementing the Monte Carlo wave-function approach to dissipation problems.  So far, I have simulated the quantum harmonic oscillator coupled to a finite temperature reservoir given in section 5A of the above paper.  I must now implement the two-oscillator system, whose Hamiltonian I am told is of the form
$$
H=\hbar\omega_aa^\dagger a+\hbar\omega_bb^\dagger b-\hbar ga^\dagger a(b^\dagger+b).
$$
My question is this:


*

*Can I form my wave-function vector in such a way that, for example, every odd entry represents the first oscillator in the corresponding state and every even entry represents the second?  As an example, if both oscillators were in the second state, the wave-function would be
$$
\left|\psi\right\rangle=\frac{1}{\sqrt{2}}\begin{pmatrix}0\\0\\1\\1\\0\\0\\\vdots\end{pmatrix}
$$

*If this is possible, what sort of form would the ladder operators $a,a^\dagger,b,b^\dagger$ take and would the corresponding number operators for each oscillator still be $a^\dagger a$ and $b^\dagger b$?


Sorry if this is a dumb question, I have not yet taken a real QM course.  Thanks for the help.
 A: I'll add to Trimok's answer in the light of your further question as to "what does it look like in higher dimensional bases. The operators going to be higher order tensors now?". Since you're "implementing" as you say, I presume you are wanting to do calculations in software.
You can, if you like, represent the state vector 
$\psi = \sum_{m,n} \psi_{m,n}|m,\,n \rangle$
in the "product" basis (I will steer clear of the word "tensorial" for now) as a two-dimensional array $\psi_{n,m}$. But it is important to recall that this is now NOT a tensor in the sense of being a multilinear homogeneous functional of two vectors and it is NOT a matrix in the sense of being homogeneous linear operator on a space of vectors. It is still essentially a one-dimensional, discrete vector, normalized so that:
$\sum_{m,n} \left|\psi_{m,n}\right|^2 = 1$
and the coefficient $\psi_{m,n}$ is simply the probability amplitude that the first oscillator is in its m-photon state and the second in its n-photon state (I'm using the words "n-photon-state" simply to mean a harmonic oscillator raised to the $n^{th}$ energy eigenstate with energy $n\,h\,\nu_0$ above the ground state: here $\nu_0$ is the oscillator's natural frequency). 
Let's discuss the product basis of two, finite-dimensional quantum systems. Making it really simple, let's think of two, three-dimensional systems. Suppose a quantum measurement $\hat{\mathrm{M}}$ can take the values $-1,\,0,\,1$, and we choose a basis wherein $\hat{\mathrm{M}}$ is diagonal. Our system state now has the form:
$\left(\begin{array}{c}\psi_{-1}\\\psi_0\\\psi_1\end{array}\right)$
where $\psi_m$ is the probability amplitude that we will measure value $m$ with $\hat{\mathrm{M}}$ and:
$\hat{\mathrm{M}} = \left(\begin{array}{ccc}-1&0&0\\0&0&0\\0&0&1\end{array}\right)$
Now let's do a quantum system comprising two quantum systems of the kind just discussed, perhaps coupled. Now the basis states are:
$\begin{array}{l}
\left.\left|-1,\,-1\right.\right>\\
\left.\left|-1,\,0\right.\right>\\
\left.\left|-1,\,1\right.\right>\\
\\
\left.\left|0,\,-1\right.\right>\\
\left.\left|0,\,0\right.\right>\\
\left.\left|0,\,1\right.\right>\\
\\
\left.\left|1,\,-1\right.\right>\\
\left.\left|1,\,0\right.\right>\\
\left.\left|1,\,1\right.\right>\end{array}$
There are now $9 = 3\times3$ basis vectors and the general state vector will be:
$\left(\begin{array}{c}\psi_{-1,-1}\\\psi_{0,-1}\\\psi_{1,-1}\\\psi_{-1,0}\\\psi_{0,0}\\\psi_{1,0}\\\psi_{-1,1}\\\psi_{0,1}\\\psi_{1,1}\end{array}\right)$
Operators and observables (i.e. operators like $\hat{\mathrm{M}}$ above) are now $9\times9$ matrices. In a general, coupled system, they are general, $9\times9$ Hermitian matrices. The matrix for the observable $\hat{\mathrm{M}}$ applied to the second quantum system alone is simply a block diagonal $9\times9$ matrix with three identical $3\times3$ copies of the $3\times3$ version of $\hat{\mathrm{M}}$ above along the leading diagonal. The matrix for observable $\hat{\mathrm{M}}$ applied to the first quantum system alone is the Kronecker Product of the $3\times3$ version of $\hat{\mathrm{M}}$ with the $3\times3$ identity matrix , $i.e.$
$\hat{\mathrm{M}}_{9\times9} = \hat{\mathrm{M}}_{3\times3} \otimes \mathrm{I}_{3\times3}$
Slightly more generally, the observables representing the linear combination $\alpha \hat{\mathrm{A}} + \beta \hat{\mathrm{B}}$ where $\hat{\mathrm{A}}$ is a $3\times3$ observable applying to the first system alone and  $\hat{\mathrm{B}}$ is a $3\times3$ observable applying to the second system alone is:
$\alpha\,\hat{\mathrm{A}} \otimes \mathrm{I}_{3\times3} + \beta\, \mathrm{I}_{3\times3}\otimes\hat{\mathrm{B}}$
General observables and operators for this combined system do not in general have this product structure: they are general Hermitian $9\times9$ matrices.
So now I've reduced the system state to a column vector, but it is a little different to your notation in your question. Your notation in your question will only work for factorisible or unentangled states

Edit: 
To keep things simple, let's assume your oscillators as isolated systems are still finite (say $N$)-dimensional quantum systems (this is the way it's going to be in software!). Hereafter, let the "First" or $\alpha$ oscillator be the one with Hamiltonian and ladder operators  $\hbar\omega_a a^\dagger a$, $a^\dagger$ and $a$ when in isolation and let the "Second" or $\beta$ be the one with $\hbar\omega_b b^\dagger b$, $b^\dagger$ and $b$ when in isolation. Let $\sum_{m,n} \psi_{m,n} |m,\,n \rangle$ be the general state where $\psi_{m,n}$ is the probability amplitude that the First or $\alpha$ oscillator is in the $m^{th}$ raised state and that the second or $\beta$ is in the $n^{th}$. Assume now that we stack the probability amplitudes into the $N^2$-element column vector as:
$\Psi = \left(\begin{array}{c}\psi_{0,0}\\\psi_{1,0}\\\psi_{2,0}\\\vdots\\\psi_{N,0}\\\psi_{0,1}\\\psi_{1,1}\\\psi_{2,1}\\\vdots\\\psi_{N,1}\\\psi_{0,2}\\\psi_{1,2}\\\psi_{2,2}\\\vdots\\\psi_{N,2}\\\vdots\end{array}\right)$
Then, to promote any $N\times N$ operator or observable $\hat{\mathrm{A}}_{N}$ acting on the $N$-dimensional state space of the First or $\alpha$ oscillator in isolation to a corresponding operator  $\hat{\mathrm{A}}_{N^2}$ in the $N^2$-dimensional product space, we form the Kronecker product:
$\hat{\mathrm{A}}_{N^2} = \hat{\mathrm{A}}_{N} \otimes \mathrm{I}_{N\times N}$
where $\mathrm{I}_{N\times N}$ is the $N\times N$ identity. To promote any operator or observable $\hat{\mathrm{B}}_{N}$ acting on the $N$-dimensional state space of the Second or $\beta$ oscillator in isolation to a corresponding operator  $\hat{\mathrm{B}}_{N^2}$ in the $N^2$-dimensional product space, we form the Kronecker product:
$\hat{\mathrm{B}}_{N^2} = \mathrm{I}_{N\times N} \otimes \hat{\mathrm{B}}_{N}$
(take heed that this Kronecker product is the other way around). Note that now any promoted operator acting on the $\alpha$ oscillator commutes with any promoted operator acting on the $\beta$ oscillator:
$\hat{\mathrm{A}}_{N^2} \hat{\mathrm{B}}_{N^2} = \left(\hat{\mathrm{A}}_{N} \otimes \mathrm{I}_{N\times N}\right) \left(\mathrm{I}_{N\times N} \otimes \hat{\mathrm{B}}_{N}\right) = \left(\hat{\mathrm{A}}_{N} \mathrm{I}_{N\times N}\right)\otimes\left(\mathrm{I}_{N\times N} \hat{\mathrm{B}}_{N}\right) = \hat{\mathrm{A}}_{N}\otimes\hat{\mathrm{B}}_{N}$
$\hat{\mathrm{B}}_{N^2} \hat{\mathrm{A}}_{N^2} = \left(\mathrm{I}_{N\times N} \otimes \hat{\mathrm{B}}_{N}\right) \left(\hat{\mathrm{A}}_{N} \otimes \mathrm{I}_{N\times N}\right) = \left(\mathrm{I}_{N\times N} \hat{\mathrm{A}}_{N}\right)\otimes\left(\hat{\mathrm{B}}_{N}\mathrm{I}_{N\times N}\right) = \hat{\mathrm{A}}_{N}\otimes\hat{\mathrm{B}}_{N}$
which you can check with the algebraic rules given on the Wiki page for the Kronecker product. There is a subtlety here when coding, but first let me complete my description of how the Hamiltonian looks. Applying the above promotion rules it is:
$\hat{\mathrm{H}}=\hbar\omega_a \left(a^\dagger a\right) \otimes \mathrm{I}_{N\times N} +\hbar\omega_b \mathrm{I}_{N\times N}\otimes\left(b^\dagger b\right)-\hbar g \left(a^\dagger a\right)\otimes(b^\dagger+b)$
Lastly, take good heed of the subtlety that tripped me up in software once. Going back to our basic commutativity law: 
$\hat{\mathrm{A}}_{N^2} \hat{\mathrm{B}}_{N^2} = \hat{\mathrm{B}}_{N^2} \hat{\mathrm{A}}_{N^2} = \hat{\mathrm{A}}_{N}\otimes\hat{\mathrm{B}}_{N}$
i.e. both ways around, the commuting product equates to $\hat{\mathrm{A}}_{N}\otimes\hat{\mathrm{B}}_{N}$, and the Kronecker product order in this latter quantity is set by how you define your state vector, not by anything else. So, if we would have stacked the probability amplitudes the "other way around" so that now our state vector were:
$\Psi = \left(\begin{array}{c}\psi_{0,0}\\\psi_{0,1}\\\psi_{0,2}\\\vdots\\\psi_{0,N}\\\psi_{1,0}\\\psi_{1,1}\\\psi_{1,2}\\\vdots\\\psi_{1,N}\\\psi_{2,0}\\\psi_{2,1}\\\psi_{2,2}\\\vdots\\\psi_{2,N}\\\vdots\end{array}\right)$
then the Kronecker products in the promotion rules would be the other way around:
$\hat{\mathrm{A}}_{N^2} =  \mathrm{I}_{N\times N} \otimes \hat{\mathrm{A}}_{N}$
$\hat{\mathrm{B}}_{N^2} =  \hat{\mathrm{B}}_{N} \otimes \mathrm{I}_{N\times N} $
and we would have 
$\hat{\mathrm{A}}_{N^2} \hat{\mathrm{B}}_{N^2} = \hat{\mathrm{B}}_{N^2} \hat{\mathrm{A}}_{N^2} = \hat{\mathrm{B}}_{N}\otimes\hat{\mathrm{A}}_{N}$
i.e. same commutativity, but different end product. This kind of thing is not very apparent in the notation often used in the product space where one just gives the commutation rule that ladder operators from different oscillators commute, nor should it need to be. It's just something you need to be a little bit careful of in coding. Trimok's edited answer shows the explicit index mapping rules which will be useful to you in software, but it is the same as the Kronecker product notation above. Be careful that, as with my Kronecker product notation, Trimok's rules assume the probability amplitudes in the product space are stacked one particular way into the state column vector. Incidentally, the Hamiltonian in this second stacking of the state column vector:
$\hat{\mathrm{H}}=\hbar\omega_a \mathrm{I}_{N\times N}\otimes \left(a^\dagger a\right)   +\hbar\omega_b \left(b^\dagger b\right)\otimes\mathrm{I}_{N\times N}-\hbar g \left(b^\dagger+b\right)\otimes \left(a^\dagger a\right)$
so that all the Kronecker products have swapped order.
If you get to countably infinite dimensions instead of $N$ for true harmonic oscillators (of course you won't in software :) ), then I have seen both "infinite" Kronecker product notation, which is intuitive but also a stacking of the state vector so that two countably infinite sequences are flattened into one countably infinite sequence in the same way that fractions as ordered pairs (represent $\frac{3}{4}$ by $(3,4)$) can be put in one-to-one correspondence with $\mathbb{N}$ to prove the countability of $\mathbb{Q}$ (see here for example). The promotion rules are then much more tangled.
A: You may use a tensorial harmonic oscillator energy basis $|n,m \rangle = |n \rangle |m \rangle = |n \rangle \otimes |m \rangle$, ($n \ge 0, m \ge0$, $n$ and $m$ integers) with operators $a, a^+$ acting on the states $n$, and operators $b, b^+$ acting on states $m$: 
$a^+|n,m \rangle = \sqrt{n+1}|n+1,m \rangle$
$a|n,m \rangle = \sqrt{n}|n-1,m \rangle$
$a^+a|n,m \rangle = n|n,m \rangle$
$b^+|n,m \rangle = \sqrt{m+1}|n,m+1 \rangle$
$b|n,m \rangle = \sqrt{m}|n,m-1 \rangle$
$b^+b|n,m \rangle = m|n,m \rangle$
$a^+a(b^+ +b)|n,m \rangle = n(\sqrt{m+1}|n,m+1 \rangle + \sqrt{m}|n,m-1 \rangle )\rangle$
The states $|n,m \rangle $ are normed and orthogonal : $\langle n' m'|nm \rangle = \delta_{nn'} \delta_{mm'}$
An element of the hamiltonian matrix can be expressed as : 
$$H_{(n'm') (nm)} = \langle n' m'|H|nm \rangle $$
A energy eigenfunction is a state $\psi = \sum_{n,m} \psi_{nm}|nm \rangle$, such that $H|\psi \rangle = E \psi \rangle$, where $E$ is the energy eigenvalue.
[EDIT]
The link between the usual matrix representation for H and the use of the tensorial basis was given above : $H_{(n'm') (nm)} = \langle n' m'|H|nm \rangle $
So, how it works: we have to find a correspondance betwenn $H_{I'I}$ and $H_{(n'm') (nm)}$
We are going to limit us to $0 \le n < N,  0 \le m < N,$.
So the number of possible different $I$ or $I'$ is $N^2$ $(0 \le I < N^2,  0 \le I' < N^2)$, so $H_{I'I}$ is a $N^2 *N^2$ matrix.
This $N^2 *N^2$ matrix apply on a $N^2$ dimensional vector $\psi_I$. So, we have to find a correspondance between $\psi_I$ and $\psi_{nm}$ too.
To be sure, you have, for instance, $H_{I'I} \psi_I = H_{(n'm') (nm)}\psi_{nm}$
The correspondance can be written in 2 ways.
Starting from $n,m$, you construct $I=(nm)$ as : $$I = Nn+m$$
Starting from $I$, you can obtain $n,m$  as :
$$m = I\mod N, \quad  n = I/N$$
That is, $n$ is the result of the integer division of $I$ by $N$, and $m$ is the 
remainder of this division, that is $I$ modulo $N$.
With this correspondance ,you can write your vector $\psi_I$ as : 
$$\psi_0,\psi_1,\psi_2....,\psi_{N^2-1} = \psi_{0~0},\psi_{0~1}...\psi_{0~N-1} \psi_{1~0},\psi_{1~1}...\psi_{1~N-1}...\psi_{2~0}...\psi_{N-1~0},\psi_{N-1~1}...\psi_{N-1~N-1}$$
Let's take an example with $N=5$.  The $H_{II'}$ is a $25*25$ matrix. $\Psi_I$ is a 25 dimensional-vector.
What is $\psi_{24}$? 
From $I=24$, we get $n=24/5 = 4$, and $m = 24 ~mod ~5 = 4$, so the correspondance is : 
$$\psi_{24} = \psi_{(4~4)}$$
What is $H_{17~~9}$? 
From $I'=17$, we get $n=17/5 = 3$, and $m = 17 ~mod ~5 = 2$, so the correspondance is :
From $I=24$, we get $n=9/5 = 1$, and $m = 9 ~mod ~5 = 4$, so the correspondance is :
$$H_{17~~9} = H_{(3 ~2) (1~4)}$$

Now, how it works, for instance, for finding a representation $I'I$ for practical operators. For instance, take $a^+$. We have :  
$$a^+|n,m \rangle = \sqrt{n+1}|n+1,m \rangle$$
This means that the only non-null terms of the operator are : 
$$a^+_{(n+1~m)(n~m)} = \sqrt{n+1}$$
With the correspondance law, this could be written ;
$$a^+_{I'~~~I} = a^+_{N(n+1)+m~~~Nn+m} = \sqrt{n+1}$$
Or, in a more practical manner : 
$$a^+_{I+N~~~I} = \sqrt{(I/N)+1}$$
where, as usual, $I/N$ is a integer division.
So, with these rules, you can build the entire hamiltonian $N^2*N^2$ matrix $H_{IJ}$, and then use some computer program to find eigenvalues, for instance.
