# Solving a differential eq. for the number operator in the Heisenberg picture

So I've been working from a master equation to derive a Heisenberg picture equation for the number operator on my system, using the method explained here; $$\mathcal{L}(\rho) = \dot{\rho}$$ $$\mathcal{L}^{\dagger}(a^\dagger a) = \frac{d}{dt}a^{\dagger}a$$ From which I've derived an equation of the form $$\frac{d}{dt}a^{\dagger}a = i E-A(t)B+C(t)D$$ Where $E,B,D$ are functions of the creation and annhilation operators, and $A(t),C(t)$ are fairly simple functions of $t$.

The problem is that now I don't know how to solve this. Do I need to find the equations of motion for $E,B,D$, and then solve it as a system of coupled equations? Or should I be able to do it in terms of initial values for those operators?

Edit in response to udrv's comment below:

The original master equation looks like $$\dot{\rho} = -\frac{i}{\hbar}[H,\rho(t)]-A(t)[q,[q,\rho(t)]]+C(t)[q,[p,\rho(t)]]$$ with $H = \frac{p^2}{2m}+\frac{m}{2}\omega^2q^2$. Converting all of this into ladder operators and using the method linked to, I found the above form for the eqm of the number operator, with $E,B,D$ being combinations of terms containing up to four ladder operators.

My thinking is that if I were to find the eqm's for all the operators in $E,B,D$, I would imagine that I would have to find eqm's for terms new terms generated in the process, but that this would not proliferate indefinitely. I think I wouldn't end up with any terms containing more than four ladder operators, and that there are a finite number of terms containing four or less; so ultimately there ought to be a solvable set of coupled equations. Does this sound correct?

• In the notation of the link you give, how do your $L_k$-s look like? Usually one is left with an (infinite) hierarchy of eom-s that have to be truncated by some set of approximations. But in some very simple cases it may be possible to do better. – udrv Jan 3 '17 at 8:53
• I've edited the main question to address the above comment. I can also post what $E,B,D$ are if that's useful. – Dan Goldwater Jan 3 '17 at 16:27
• That should be exactly solvable. I can try to write up something later, but you want to look up lit. on exactly solvable Lindblad models. See for instance arxiv.org/abs/1511.03347v2 (Sec.3), arxiv.org/abs/1609.07249, but there should be an even nicer algebraic solution out there. – udrv Jan 3 '17 at 17:42

This is not a solution, but some hopefully useful notes that may lead to other things useful towards a genuine solution, and do not fit into a comment.

I mentioned in a comment that the model should be exactly solvable. The issue is that your coefficients A and C are time-dependent and things get complicated. If there is any way to approximate/replace them by time-independent values, everything simplifies a lot, to the point where the problem becomes algebraically tractable. Simple oscillatory exponentials might also work to some degree.

In any case, the general idea for an exact solution is that the master equation's generator is of a form ${\mathcal L} = {\mathcal L}_0 + {\mathcal L}_1$, where ${\mathcal L}_0$ and ${\mathcal L}_1$ yield relatively simple nested commutation relations (as superoperators), which in turn allow a suitable disentangling of the time evolution (for either ${\mathcal L}$ or ${\mathcal L}^\dagger$) into simpler tractable factors. The time evolution of the number operator and its average may then be calculated exactly. An approach developed along these lines for the $C=0$ case and based on a su(1,1) algebra is given in http://lanl.arxiv.org/abs/0710.2724. Details on their su(1,1) method in http://lanl.arxiv.org/abs/quant-ph/0112090.

Anyway, here is another way to do it, possibly in more transparent notation.

Begin by defining commutators and anti-commutators of $p$ and $q$ as simple superoperators, $$Q_+({\hat u}) = \frac{1}{2}\lbrace q {\hat u} + {\hat u} q \rbrace\\ Q_-({\hat u}) = q {\hat u} - {\hat u} q\\ \Pi_+({\hat u}) = \frac{1}{2}\lbrace p {\hat u} + {\hat u} p \rbrace\\ \Pi_-({\hat u}) = p {\hat u} - {\hat u} p$$ where ${\hat u}$ is some arbitrary operator. Then observe that

1) $Q_\pm$ and $P_\pm$ satisfy very nice commutation relations of their own: $$\left[ Q_-, \; Q_+\right] = \left[ \Pi_-, \; \Pi_+\right] = \left[ Q_-, \; \Pi_+\right] = \left[ Q_+, \; \Pi_- \right] = 0\\ \left[ Q_\pm, \; \Pi_\pm\right] = i$$

2) In terms of $Q_\pm$ and $P_\pm$, your generator simplifies to $${\mathcal L} = {\mathcal L}_0 + {\mathcal L}_1(t)\\ {\mathcal L}_0 = \Pi_+ \Pi_- + \omega^2 Q_+ Q_- \\ {\mathcal L}_1(t) = -A(t) Q_-^2 + C(t) \Pi_+ Q_-$$ Note: There is nothing new here, and in principle one may dispense entirely with points (1) and (2), but the compact notation and the simple algebra do make superoperator calculus much more straightforward.

3) As anticipated, it is now easy to see that ${\mathcal L}_0$ and ${\mathcal L}_1$ also have simple nested commutation relations: $$\left[ {\mathcal L}_0, \;{\mathcal L}_1(t)\right] = 2iA(t) \Pi_+ Q_- - iC(t)\left(\Pi_+^2 - \omega Q_-^2 \right) \equiv {\mathcal K}_1(t)\\ \left[ {\mathcal L}_0, \left[ {\mathcal L}_0, \;{\mathcal L}_1(t) \right] \right] = \left[ {\mathcal L}_0, \;{\mathcal K}_1(t)\right] = 4 \left[ A(t) \left(\Pi_+^2 - \omega Q_-^2 \right) + \omega^2 C(t)\; \Pi_+ Q_- \right] \equiv {\mathcal K}_2(t)\\ \left[ {\mathcal L}_0, \left[ {\mathcal L}_0, \left[ {\mathcal L}_0, \;{\mathcal L}_1(t) \right] \right] \right] = \left[ {\mathcal L}_0, \;{\mathcal K}_2(t)\right] = 2\omega^2 {\mathcal K}_1(t)\;,\;\;\;\text{etc.}$$ where $$\left[ {\mathcal K}_1(t), \; {\mathcal K}_2(t)\right] = 0\;,$$ and also $$\left[ {\mathcal L}_1(t), \left[ {\mathcal L}_0, \;{\mathcal L}_1(t) \right] \right] = \left[ {\mathcal L}_1(t), \;{\mathcal K}_1(t)\right] = 0\\ \left[ {\mathcal L}_1(t), \left[ {\mathcal L}_0, \left[ {\mathcal L}_0, \;{\mathcal L}_1(t) \right] \right] \right] = \left[ {\mathcal L}_1(t), \;{\mathcal K}_2(t)\right] = 0\;, \;\;\;\text{etc.}$$

If $A$ and $C$ were time-independent, one could apply the Zassenhaus formula to the evolution $e^{{\mathcal L} t} = e^{\left({\mathcal L}_0 + {\mathcal L}_1\right) t}$, $$e^{\left({\mathcal L}_0 + {\mathcal L}_1\right) t} = e^{{\mathcal L}_0 t} e^{{\mathcal L}_1 t} e^{-(t^2/2!) \left[ {\mathcal L}_0, \;{\mathcal L}_1\right]} e^{(t^3/3!) \left(\left[ {\mathcal L}_0, \left[ {\mathcal L}_0, \;{\mathcal L}_1 \right] \right] + 2 \left[ {\mathcal L}_1, \left[ {\mathcal L}_0, \;{\mathcal L}_1 \right] \right] \right)} \cdot \dots$$ and observe that all factors after the first two, $e^{{\mathcal L}_0 t} e^{{\mathcal L}_1 t}$, mutually commute and that the exponents reduce to a single term each and can be eventually summed up.

On the other hand, the time-dependence of $A$ and $C$ turns the overall evolution into $$e^{{\mathcal L} t} = {\mathcal T} e^{\int_0^t{d\tau\;\left({\mathcal L}_0 + {\mathcal L}_1 \right)} }$$ where $\mathcal T$ stands for time ordering, and the procedure complicates as usual. But perhaps in simpler cases it may still provide some good guidelines. For the current problem the final target is of course an explicit algebraic expression for the time-dependence of the number operator.

• Very elegant formulation of the problem. I don't know if it would work or not, but the Magnus representation of the propagator might be useful here since the commutators are "nice". – Mark Mitchison Jan 9 '17 at 16:03
• Thanks, half-baked unfortunately. The Magnus representation is indeed what one defaults to eventually, but it takes some work. – udrv Jan 9 '17 at 22:39

I've worked in a different direction to that posted by @urdv, and arrived at a sperate possible answer; though I'm not entirely certain it is a proper answer. My working is this:

We start again with $$\mathcal{L}(\rho) = -\frac{i}{\hbar}[H,\rho]-A(t)[q,[q,\rho]]+C(t)[q,[p,\rho]].$$ We can expand the above in terms of the ladder operators in order to find the adjoint of $\mathcal{L}$, $$\mathcal{L}^{\dagger}(O) = i k_1 \mathcal{L}_1^\dagger(O)-k_2A(t)\mathcal{L}_2^\dagger+k_3C(t)\mathcal{L}_3^\dagger(O)$$

Where $$\mathcal{L}_1^\dagger (O) = (\alpha+\beta)((a^\dagger a^\dagger+aa)O-O(a^\dagger a^\dagger +aa))+(\alpha-\beta)((a^\dagger a+a a^\dagger )O-O(a^\dagger a +a a^\dagger ))$$ using $$\alpha =\frac{-1}{2m}, \beta = \frac{m\omega^2}{2}$$ $$\mathcal{L}_2^\dagger = (a a + a^\dagger a + a a^\dagger + a^\dagger a^\dagger )O-2(a O a +a^\dagger O a +a O a^\dagger + a^\dagger O a^\dagger )+O(aa + a a^\dagger + a^\dagger a + a^\dagger a^\dagger )$$ and $$\mathcal{L}_3^\dagger (O) = i((aa+a^\dagger a-a a^\dagger -a^\dagger a^\dagger )O+2(a^\dagger Oa^\dagger -2aOa)+O(aa+a a^\dagger -a^\dagger a-a^\dagger a^\dagger )),$$ We can then use these to find the equation of motion for the number operator: $$\mathcal{L}^\dagger(a^\dagger a)= \frac{d}{dt}a^\dagger a = 2ik_1(\alpha+\beta)(aa-a^\dagger a^\dagger )-2k_2A(t)$$ We can see that this is contingent upon knowing the values of $aa$ and $a a^\dagger$, both functions of time as well. We can find eqms for these via the same method as above, giving; $$\frac{d}{dt}aa = -2ik_1((\alpha+\beta)(a^\dagger a+aa^\dagger )+2(\beta-\alpha)aa)-2k_2A(t)-2ik_3C(t)$$ and the adjoint equations

$$\frac{d}{dt}aa^\dagger =-2ik_1(\alpha+\beta)(a^\dagger a^\dagger -aa)-2k_2A(t)$$ $$\frac{d}{dt}a^\dagger a^\dagger =ik_1((\alpha+\beta)(a^\dagger a+aa^\dagger )+2(\beta-\alpha)a^\dagger a^\dagger )-2k_2A(t)+2ik_3C(t)$$

My hope was to find a neat solution by taking the second derivative of $a^\dagger a$, and plugging in the first order derivatives of the other variables as appropriate. However, as far as I can see this does not lead to a tidy solution.

Alternatively, I can plug the four first order coupled ODE's into Mathematica, and ask it to solve the system. This does yield an answer, but it's an extremely long one, far too long to write out and certainly far too long to be able to glean any meaningful understanding from.

Can anyone see a way to get a neat solution from the coupled ODE's? Either by taking the second order, or some other method?

• Using $aa^\dagger = 1+ a^\dagger a$ I think your problem actually reduces to just three coupled (first-order inhomogeneous) differential equations. – Mark Mitchison Jan 9 '17 at 16:09
• @DanGoldwater You seem to apply the Heisenberg picture eom-s with a Lindblad generator and an irreversible dynamics the way you'd apply them with a Hamiltonian and a reversible dynamics. But under a Lindblad generator the evolved of a product of operators is no longer the product of the evolved operators. For instance $e^{\mathcal{L}^\dagger t}\left(a^\dagger a\right) \neq e^{\mathcal{L}^\dagger t}\left(a^\dagger\right) e^{\mathcal{L}^\dagger t}\left(a\right)$, and similarly for commutators and time derivatives. See for instance pg.253 in uni-due.de/~hp0198/pubs/lnp2.pdf. – udrv Jan 9 '17 at 22:28
• Did you check that the evolution you work with does indeed preserve the commutation relations you casually assume and apply? – udrv Jan 9 '17 at 22:29
• @urdv, I didn't check that the commutation relations were preserved, and after reading what you've linked to it's clear that my approach isn't correct. I'm currently trying my hand at solving it using your formulation of the problem, and I'll update when I have some progress. – Dan Goldwater Jan 10 '17 at 15:31
• Don't give up on a direct calculation just yet either. What you wrote in the answer actually does hold in terms of averages, provided you interpret things correctly. This is because, by hermiticity as applied to superoperators, $$\frac{d}{dt}\langle O \rangle = Tr\left[ O \mathcal{L}(\rho(t))\right] = \left( Tr\left[ \rho(t) \mathcal{L}^\dagger(O) \right]\right)^* = Tr\left[ \mathcal{L}^\dagger(O) \rho(t) \right]$$The last equality holds because both $\rho(t)$ and $\mathcal{L}^\dagger(O)$ are hermitian operators. – udrv Jan 10 '17 at 18:05

Attempting the direct solution again, I think I made a very basic error in my first attempt, and that possibly the whole problem has an almost trivial solution. In the above, I think I expanded $H$ wrong, and that instead it can be replaced in the standard form $H = \hbar \omega (a^\dagger a+\frac{1}{2})$. Since the adjoint of this commutes with the number operator, the whole thing simplifies and we are left with no operator terms. From the top: $$\mathcal{L}(\rho) = -\frac{i}{\hbar}[H,\rho]-A(t)[q,[q,\rho]]+C(t)[q,[p,\rho]]$$ \begin{align} \frac{d}{dt}\langle O \rangle &= \mathcal{L}^\dagger( O)\\ &=i \omega \mathcal{L}_1^\dagger(O)-A(t)\mathcal{L}_2^\dagger (O)+C(t)\mathcal{L}_3^\dagger(O ) \end{align} with $\mathcal{L}_1^\dagger(O) = [a a^\dagger,O]$, and $\mathcal{L}_2^\dagger,\mathcal{L}_3^\dagger$ as in the above post. Since $\mathcal{L}_1^\dagger(a^\dagger a) =0$ and $\mathcal{L}_3^\dagger ( a^\dagger a ) = 0$, we are left with $$\mathcal{L}^\dagger(a^\dagger a ) = \frac{1}{2}A(t).$$ So $$\frac{d}{dt}\langle a^\dagger a\rangle = \frac{1}{2}A(t).$$