Heisenberg picture: harmonic oscillator operators I took an introductory quantum mechanics course in my university's fall 2019 quarter, so I am quite familiar with the Schrodinger picture we were taught. I was introduced to the creation and annihilation operators one uses to "intelligently" quantize a harmonic oscillator. I recently heard about the Heisneberg picture, which is apparently similar to the Schrodinger picture but includes time-dependent operators and time-independent wave functions for energy eigenfunctions (i.e. they do not contain an exponential rotating phase factor like $e^{\omega t}$). I am trying to solve the quantum harmonic oscillator problem in the Heisenberg picture. Specifically, I want to get the gaussian-polynomial type expressions I found for stationary states when I use the Schroedinger picture.
Here's what my attempts at googling my problem gave me:

We first start with analyzing the evolution of the operators in
  the Heisenberg picture. We have
  $$\frac{da}{dt} = i[\mathcal{H}, a] = i[ω(a^†a+\frac{1}{2}), a] = −iωa → a(t) = a(0)e^{−iωt}$$
  Similarly:
  $$\frac{da^†}{dt}=i[H, a^†] = i[ω(a^†a+\frac{1}{2}), a^†] = iωa^† → a^†(t) = a^†(0)e^{iωt}$$
  Notice that we could have found this last relationship just by taking the hermitian conjugate of the first one.
  Using these results, we can also find the time evolution of the position and momentum operators:
  $$x(t) = x(0) \cos(ωt) + \frac{p(0)}{m\omega}\sin(ωt)$$
$$p(t) = p(0) \cos(ωt) − mωx(0) \sin(ωt)$$
  and the corresponding expectation values, e.g.
  $$\left<x(t)\right> = \left<x(0)\right> \cos(ωt) + \frac{\left<p(0)\right>}{m\omega} \sin(ωt)$$
  (https://ocw.mit.edu/courses/nuclear-engineering/22-51-quantum-theory-of-radiation-interactions-fall-2012/lecture-notes/MIT22_51F12_Ch9.pdf)

The wikipedia page on the Heisenberg picture says something similar. That sounds helpful and it makes sense. I can substitute in the definitions of the $x$ and $p$ operators from the Schroedinger picture, since I was given to understand that operators in the two pictures are identical at time $t=0$.
1) What are the tools from the harmonic oscillator procedure in the Schrodinger picture that can be used here? Does the canonical commutator relationship still hold? Do the commutators of $a$ and $a^†$ become useful?
2) What is the ground state in the Heisenberg picture? Can I still say that the lowermost state is the one which, when acted upon by $a$, gives 0? I'm not sure because the operator itself is evolving over time.
3) What equation can give me the actual wave function, assuming it is possible to find the discrete energy values and find allowed energy eigenvalues?
4) Any books or articles with a complete derivation? I appreciate it if an answer can give a full derivation, but this question may get closed for that so it's not necessary. 
 A: Just unfold the recipe. Set $\hbar=1$ and $\omega=1$ for simplicity, while you can reinstate them uniquely by dimensional analysis if you need to. Use plain symbols a for constant Schroedinger operators, and time-dependence parentheses a(t) for Heisenberg time-dependent operators, where you already saw that a(0)=a. 
You already found 
$$
H= 1/2 + a^\dagger a , \leadsto \\
a(t) = e^{iHt} a e^{-iHt} = e^{-it} a , \leadsto \\
H(t)=H,  \qquad [a(t),a^\dagger(t)] =1.
$$
The Fock space, unsurprisingly, runs on constant states, so defined through Schroedinger operators, as in WP,
$$
a|0\rangle =0 , \qquad H| 0\rangle = 1/2 |0\rangle \\
a^\dagger |0\rangle \equiv |1\rangle  , \qquad \leadsto H|1\rangle = 3/2|1\rangle, ... 
$$
These are the Heisenberg picture states, constant, as required, as well. This should answer your question (2), and (1): no you definitely don't use Heisenberg oscillators to define time-independent Fock states. 
The Schroedinger states,  by contrast, are time-dependent solutions of the TDSE,
$$
e^{-it/2} |0\rangle, \qquad e^{-it~3/2} |1\rangle, ...  
$$
(Somewhat perversely, they could be defined by piling up Heisenberg oscillators. Books wisely avoid dwelling on such to not invite confusion.)
So now you assume you've run the ladder and determined the above eigenstates of the hamiltonian with the eigenvalues spaced equally, etc. To find Schroedinger wavefunctions, (3), you dot constant $\langle x|$ kets on the time-dependent Schroedinger states, so, as per the WP discussion,
$$
\psi_0 (x,t) = \langle x|e^{-it/2} |0\rangle= (m/\pi)^{1/4} e^{-it/2   -mx^2/2}, \\  \psi_1 (x,t) = e^{-it~3/2} \langle x |1\rangle= 
  (m/\pi)^{1/4} \sqrt{m/2} (x-\partial_x/m) e^{-mx^2/2} \\ =(m/\pi)^{1/4} \sqrt{m/2} ~ 2x e^{-it~3/2 -mx^2/2}  ,...
$$


*

*But the whole point of  Dirac's superior ladder method is to not have to solve the Schroedinger equation, and so consider these at all! They are truly useless! The above formulas are for formal reassurance--and of course trump direct solution of the TDSE: 19th century physics.


Thus, you never need such space wavefunctions in the  Heisenberg picture, since you  don't solve the space Schrödinger equation in matrix mechanics; still, if you wanted them, in order to connect to the sinful past, you just look at the zero-time behavior of these Schroedinger ones to get the Heisenberg ones. Any expectation value you'd need evaluate with them will involve the time-evolution in the propagator-sandwich of the operators detailed in the first link above.
A: *

*The entire procedure is almost exactly like that used in the Schroedinger picture. You just need to be clear about the definition of the Heisenberg picture and the properties of operators and kets therein. To answer question (1), yes, the canonical commutator between $\hat{x}$ and $\hat{p}$ holds in the Heisenberg picture, as mentioned on wikipedia (see section Commutation relations) and discussed in this Physics SE post. Similarly, we acquire similar forms of the commutator for $a$ and $a^\dagger$ in the Heisenberg picture.
Let's make the following definitions:
$$a=\sqrt{\frac{m\omega}{2\hbar}}\hat{x}+\mathrm{i}\sqrt{\frac{1}{2m\omega\hbar}}\hat{p};\quad a^\dagger=\sqrt{\frac{m\omega}{2\hbar}}\hat{x}-\mathrm{i}\sqrt{\frac{1}{2m\omega\hbar}}\hat{p};\quad a=a_0\mathrm{e}^{\mathrm{i}\omega t}$$
$a$ and $a^\dagger$ are time-dependent, as are $\hat{x}$ and $\hat{p}$. 
Using the earlier premise that the canonical commutation relation is availible for use, we find that another familiar commutator relationship holds:
$$[a,a^\dagger]=1$$
Similarly, you should have no trouble finding $$[H,a]=-\hbar\omega a.$$
The time-dependent operators are not treated differently when you're just doing algebraic stuff like this.

*The ground state is exactly the same, without the factor of $\mathrm{e}^{\mathrm{i}E_0t/\hbar}$. You find it in a perfectly standard way: for eigenvalue $E_n$ corresponding to ket $|n\rangle$, $$\langle n|\hbar\omega aa^\dagger|n\rangle\geq0\Rightarrow\left(E_n-\frac{1}{2}\hbar\omega\right)\langle n|n\rangle\geq 0.$$ For a ground state, we have equality for both of those expressions, which means $$E_0=\frac{1}{2}\hbar\omega;\quad a|0\rangle=0$$

*To recreate the states from the Schroedinger picture, use the definition of $a$ and plug in familiar operators for $\hat{x}$ and $\hat{p}$. $$\left(\hat{x}+\mathrm{i}\frac{\hat{p}}{m\omega}\right)|0\rangle=0\Rightarrow x\langle x|0\rangle+\frac{\hbar}{m\omega}\frac{\partial}{\partial x}\langle x|0\rangle=0$$ The result for $\langle x|0\rangle$ is the ground state wave function you've calculated in the Schroedinger picture.

*Paul Dirac's Principles of Quantum Mechanics contains a similar procedure, and doesn't go into a discussion of directly solving differential equations at all.
