How to derive annihilation and creation operators? In the book I am studying (Zettili) the author introduces annihilation and creation operators by defining them in terms of $q$ and $p$ and showing how to write the hamiltonian of the harmonic oscilator in terms of them. I want to know if there is a way to derive them. The way they are introduced in the book seems like something that magically appears. How do they arise from the theory?
 A: The way you "derive" them in real life is that you know about them from classical mechanics.
It is absolutely astonishing that history seems to have forgotten this.
Classical harmonic oscillator in Hamiltonian form
Consider a classical undriven and undamped harmonic oscillator.
Generally speaking, it's equation of motion is
$$\ddot{q} + \omega_0^2 q = 0$$
where $q$ is the coordinate of the oscillator and $\omega_0$ is the natural resonance frequency.
If we guess that the Lagrangian is
$$ \mathcal L = \frac{\beta}{2} \dot{q}^2 - \frac{1}{2 \alpha} q^2 \, ,$$
then applying the Euler-Lagrange equation gives
\begin{align}
  \frac{d \mathcal L}{dq} - \frac{d}{dt} \frac{d \mathcal L}{d \dot{q}} &= 0 \\
  \ddot{q} + \frac{1}{\alpha \beta} q &= 0 \, ,
\end{align}
which is correct for any $\alpha$ and $\beta$ such that $1 / \alpha \beta = \omega_0^2$.
Following the usual procedure to find the Hamiltonian, we get
$$H = \frac{1}{2 \alpha} q^2 + \frac{1}{2 \beta} p^2  $$
where the momentum $p$ is defined as $p \equiv \partial \mathcal L / \partial \dot{q} = \beta \dot{q}$.
Hamilton's equations of motion are
$$
\dot p = - \frac{\partial H}{\partial q} = -q / \alpha
\qquad \text{and} \qquad
\dot q =  \frac{\partial H}{\partial p} = p / \beta
$$
or combined as a matrix equation
$$\frac{d}{dt}
\begin{bmatrix} q \\ p \end{bmatrix}
=
\left(
  \begin{array}{cc}
    0 & 1/\beta \\
    - 1 / \alpha & 0
  \end{array}
\right)
\begin{bmatrix} q \\ p \end{bmatrix} \, .
$$
Rescaled variables
If we define $x \equiv A q$ and $y \equiv B p$ with the constraint that $A/B = \sqrt{\beta / \alpha}$, then our Hamilton equations of motion become
$$\frac{d}{dt}
\begin{bmatrix} x \\ y \end{bmatrix}
= \omega_0
\left(
  \begin{array}{cc}
    0 & 1 \\
    - 1 & 0
  \end{array}
\right)
\begin{bmatrix} x \\ y \end{bmatrix} \, .
$$
This is a set of first order coupled differential equations for $x$ and $y$.
To uncouple the equations, we solve for the eigenvectors and eigenvalues of the matrix.
They are$^{[*]}$
$$a \equiv x + i y = \begin{bmatrix} 1 \\ i \end{bmatrix} \text{ with eigenvalue } i \omega_0$$
and
$$ a^* \equiv x - i y = \begin{bmatrix} 1 \\ -i \end{bmatrix} \text{ with eigenvalue } -i \omega_0 \, .$$
The $a$ and $a^*$ variables have very simple time dependence, i.e.
$$\boxed{
\begin{array}{ll}
  \dot{a} = i \omega_0 a & a(t) = a(0) e^{i \omega_0 t} \\
  \dot{a}^* = -i \omega_0 a^* & a^*(t) = a^*(0) e^{-i \omega_0 t}
\end{array}
}
$$
Discussion


*

*The equation for $a$ and $a^*$ in terms of $x$ and $y$ (or in terms of $q$ and $p$) mirror the quantum equations for $\hat a$ and $\hat a^\dagger$.

*The time evolution of $a$ and $a^*$ is the same as for $a$ and $a^\dagger$ in the Heisenberg picture of quantum mechanics.
As we can see, the $\hat a$ and $\hat a^\dagger$ operators definitely didn't come just from quantum mechanics, as they are direct analogues to the variables that diagonalize the Hamiltonian evolution matrix of the classical harmonic oscillator.
Further notes
A huge part of why $a$ and $a^*$ are useful is that they make it easy to analyze coupled and driven systems, in particular using perturbative methods.
For example, using the rotating frame and rotating wave approximation, the interaction term of two interacting oscillators
$$H_\text{interact} = g x_1 x_2 \, .$$
can be rewritten as
$$H_\text{interact} = \frac{g}{4} (a_1 a_2^* + a_1^* a_2) \, ,$$
which is exactly analogous to the quantum interaction term
$$\frac{g}{4} ( a_1 a_2^\dagger + a_1^\dagger a_2) \, .$$
Footnotes
$[*]$: I wasn't careful to make the signs come out the same way as they do in quantum mechanics, so actually we have the eigenvalues swapped compared to what you'd get with $\hat a$ and $\hat a^\dagger$. I also wasn't careful to make the normalization come out the same way as it does in quantum mechanics. This is a pretty simple detail that could be edited by someone else, if they're interested.
A: 
I want to know if there is a way to derive them

One can solve for the energy eigenstates $|E_n\rangle$ of the quantum harmonic oscillator (QHO) without the use of ladder operators and find that
$$H|E_n\rangle = (n + 1/2)\,\hbar\omega\,|E_n\rangle$$
which implies that the QHO Hamiltonian can be written in the form
$$H = (N + 1/2)\hbar\omega$$
where $N$ is the number operator
$$N|E_n\rangle = n|E_n\rangle$$
Now, suppose that the (hermitian) operator $N$ is the product of a (non-hermitian) operator and its hermitian conjugate:
$$N = a^\dagger a$$
where
$$[a,a^\dagger] = c$$
which implies
$$aa^\dagger = N + c$$
An 'obvious' choice for $c$ is $c = 1$ so that $H$ can be written as
$$H = \frac{a^\dagger a + aa^\dagger}{2}\hbar\omega $$
Note that, so far, there's no mention of ladder operators.  We've simply written the Hamiltonian in terms a (hermitian) number operator and then written that operator as a product of a (non-hermitian) operator and its hermitian conjugate.
So, take the next step of asking "but what do these operators $a^\dagger$ and $a$ do to an energy eigenstate"?
Look at the commutator of $N$ with $a^\dagger$ and $a$:
$$[N,a^\dagger]|E_n\rangle = a^\dagger[a,a^\dagger]|E_n\rangle = a^\dagger|E_n\rangle$$
$$[N,a]|E_n\rangle = [a^\dagger,a]a|E_n\rangle = -a|E_n\rangle$$
thus
$$N\,a^\dagger|E_n\rangle = (n + 1)\,a^\dagger|E_n\rangle \Rightarrow a^\dagger|E_n\rangle \propto |E_{n+1}\rangle$$
$$N\,a|E_n\rangle = (n - 1)\,a|E_n\rangle  \Rightarrow a|E_n\rangle \propto |E_{n-1}\rangle$$
and see that these are in fact the ladder operators that connect energy eigenstates to their adjacent eigenstates. 
In summary, we 'derived' the ladder operators starting only with the fact that, for the QHO,
$$H|E_n\rangle = (n + 1/2)\,\hbar\omega\,|E_n\rangle$$
A: One way to do this is to search for operators $\hat a_-$ and $\hat a_+$ so that
$$
[\hat H,\hat a_\pm]=\pm \hbar\omega \hat a_\pm\, , \tag{1}
$$
where $\hat H$ is the harmonic oscillator Hamiltonian.
If can you find these, then you have
\begin{align}
\hat H\left[\hat a_+\vert n\rangle\right] &=[\hat H,\hat a_+]\vert n\rangle 
+\hat a_+\hat H \vert n\rangle\, ,\\
&=\hbar\omega \hat a_+\vert n\rangle +(n+\textstyle\frac{1}{2})\hbar\omega\hat a_+\vert n\rangle\, ,\\
&=(n+1+\textstyle\frac{1}{2})\hbar\omega\hat a_+\vert n\rangle
\end{align}
since $\hat H\vert n\rangle = (n+\frac{1}{2})\hbar\omega\vert n\rangle$.  This shows that 
$\hat a_+\vert n\rangle$ must be proportional to $\vert n+1\rangle$, with an analogous reasoning valid for $\hat a_-\vert n\rangle$.  
This will give you the usual $\hat a_\pm$ up to an overall “normalization” factor, in the sense that, if $\hat a_\pm$ satisfy (1), so will $\kappa \hat a_\pm$ for any $\kappa$.  You can fix this “normalization” by requiring that
$$
\hat a_+\hat a_-=\frac{\hat H}{\hbar \omega}-\frac{1}{2}\hat 1\, .
$$
Note that, if $H=\frac{1}{2}p^2+\frac{1}{2}x^2$ is the classical mechanical Hamiltonian for the harmonic oscillator, combinations of the type 
$\alpha = x+ ip$ and $\alpha^*=x-ip$ occur naturally as they are the combinations that have a simple time evolution, i.e. $d\alpha/dt=i\alpha$ and $d\alpha^*/dt=-i\alpha^*$; it’s not really surprising that similar combinations should occur in quantum mechanics as well.
A: Another way would be to observe that $$
\begin{align}
[\hat x^2,~\hat p]&=\hat x(\hat p \hat x+i\hbar) - \hat p\hat x\hat x\\
&= 2i\hbar \hat x
\end{align}
$$
while $[\hat p^2,~\hat x]=-2i\hbar \hat p$. So then we could see that this “flip-flopping” pattern for the harmonic oscillator Hamiltonian supports a combination,$$
\left[\hat H,~\hat x+\alpha\hat p\right] =\frac12\left[k \hat x^2 + m^{-1}\hat p^2,~ \hat x+\alpha\hat p\right] =i\hbar\left(\alpha k \hat x - m^{-1}\hat p\right).
$$
So we are looking at some sort of combination where $[\hat H,~\hat q]=\lambda \hat q$ but it only works if $\alpha=-(\alpha m k)^{-1},$ so that $\alpha=\pm i/\sqrt{mk}$ and hence $\lambda=i\hbar \alpha k=\mp \hbar\omega.$
The key here is that
$$[\hat H,~\hat q]=\lambda \hat q$$ means that given an eigenfunction $|\Psi\rangle$ with eigenvalue $v$, $\hat q|\Psi\rangle$ is an eigenfunction 
of $\hat H$ with eigenvalue $v+\lambda$, and you have the ladder operators.
