Let's quickly review the quantum harmonic oscillator. We have a single particle moving in one dimension, so the Hilbert space is $L^2(\mathbb{R})$: the set of square-integrable complex functions on $\mathbb{R}$. The harmonic oscillator Hamiltonian is given by
$$H= \frac{P^2}{2m} + \frac{m\omega^2}{2}X^2$$
where $X$ and $P$ are the usual position and momentum operators: acting on a wavefunction $\psi(x)$ they are $X \psi(x) = x\psi(x)$ and $P \psi(x) = -i\hbar\ \partial \psi / \partial x$. Of course, we can also think of them as acting on an abstract vector $|\psi\rangle$.
By letting $P \to -i\hbar\ \partial/\partial x$ we could solve the time independent Schrödinger equation $H \psi = E \psi$, but this is a bit of a drag. So instead we define operators $a$ and $a^\dagger$ as in your post. Notice that the definition of $a$ and $a^\dagger$ has nothing whatsoever to do with our Hamiltonian. It just so happen that these definitions are convenient because the Hamiltonian turns out to be $\hbar \omega (a^\dagger a + 1/2)$.
For convenience we define the number operator $N = a^\dagger a$; at this stage number is just a name with no physical interpretation. Using the commutation relation $[a,a^\dagger] = 1$ and some algebra we notice that $N$ has a nondegenerate spectrum given by the natural numbers. In other words, the eigenvalues of $N$ are $\{0,1,2,\dots\}$, and to each eigenvalue $n$ there corresponds a single state $|n\rangle$ with $N|n\rangle = n |n\rangle$. Notice that, again, $N$ is independent of our Hamiltonian. However, because the Hamiltonian turns out to be $\hbar \omega (N+1/2)$ we immediately know that the states $|n\rangle$ are its eigenvectors, with energies $\hbar \omega (n + 1/2)$.
Now you are given a different Hamiltonian. The Hilbert space is still exactly the same, and so are $a$, $a^\dagger$ and $N$, because their definition had nothing to do with the original Hamiltonian. You can still use their properties to find energies, eigenvectors, and so on. The states $|n\rangle$ are still the eigenstates of $N$, though a priori they might not be eigenstates of the new $H$ (exercise 31 asks you to prove that they in fact are eigenstates of the new $H$). The important point here is that operators are (usually) defined independently of the Hamiltonian. They characterize the physical system. After all, you know that there are operators $X$ and $P$, and you have no qualms about using them with different Hamiltonians. The Hamiltonian gives the energy and the time evolution, but the observables and related operators are independent of your choice of Hamiltonian.
About the physical interpretation... exercise 31 asks you to prove that $H=\hbar\omega_0 N + \hbar \omega_1 (N^2-N)$; notice that we have gotten rid of $\hbar\omega_0 /2$ since it is just a constant. I would usually expect $\omega_1$ to be smaller than $\omega_0$ so this is a small perturbation (for small $n$ at least), but we don't really care about that right now. You can see that $|n\rangle$ are still the eigenstates of the Hamiltonian; all we did is shift the energies by an amount $\hbar \omega_1 (N^2-N)$.
$\renewcommand\mean[1]{\langle #1\rangle}$
$\renewcommand\norm[1]{||#1||}$
$\renewcommand\h{\hbar}$
$\renewcommand\ket[1]{|#1\rangle}$
$\renewcommand\expval[2]{\langle #1|#2|#1\rangle}$
$\renewcommand\braket[2]{\langle #1|#2\rangle}$
$\renewcommand\Braket[3]{\langle #1|#2|#3\rangle}$
OP's original question
The step that really confuses me here is that we seemingly use
\begin{equation} \Braket{x}{p}{0}=p\braket{x}{0}\end{equation} where I fail to see why this holds true.
has already been answered in the replies of @Cosmas and @Zack. I won't go over these again. Instead, I propose here an alternative method of finding the ground state of the Harmonic Oscillator, using this time the concept of coherent states. I am going to assume the reader knows nothing about them, so I will build up everything here to avoid losing anyone.
Introduction to coherent states
There are three different - equivalent - ways to define a coherent state. For the sake of answering your question, I will only need one of them (sadly the less intuitive of them all) so I won't go all the way to introduce them all. The plan is thus the following: defining the so-called Heisenberg Uncertainty Principle and using it to define a coherent state. From there, I will be able to build a mathematical structure that will allow us to find the ground state of the Quantum Harmonic Oscillator.
Theorem 1 (Heisenberg Uncertainty Principle) For any operators $\hat{A}$ and $\hat{B}$, we have that $\Delta\hat{A}\hat{B}\geq\frac{1}{2}\norm{\mean{[\hat{A},\hat{B}]}}$.
In the case of the position and momentum operators, one has that $[\hat{X},\hat{P}] = i\h$. Applying this to the uncertainty principle, we find the famous relation $\Delta\hat{X}\Delta\hat{P}\geq \frac{\h}{2}$.
Definition 2 (Coherent states) A coherent state is defined as a state for which the Heisenberg Uncertainty Principle (applied to the operators $\hat{X}$ and $\hat{P}$ in the context of your question) is saturated, such that $\Delta\hat{X}\Delta\hat{P}=\frac{\h}{2}$.
Wave function of a coherent state
My goal here is to construct a wave function that depicts how a coherent state behaves. Let us introduce $\mathcal{O}(\lambda)=f^\dagger(\lambda)f(\lambda)$, where we define $f(\lambda) = \hat{A'}+i\lambda\hat{B'}, \hat{A'} = \hat{A}-\mean{\hat{A}}$ and $\hat{B'} = \hat{B}-\mean{\hat{B}}$.
Let us observe that the mean of $\mathcal{O}$ is positive. Indeed, for any $\ket{\psi}\in\mathcal{H}$, one has that
\begin{equation}
\expval{\psi}{f^\dagger(\lambda)f(\lambda)} = \norm{f(\lambda)\ket{\psi}}^2 \geq 0
\end{equation}
However, one also has that
\begin{equation}\tag{1}\label{mean value O}
\mean{\mathcal{O}(\lambda)} = \lambda^2\mean{\hat{B'}}+\lambda\mean{i[\hat{A'},\hat{B'}]}+\mean{\hat{A'}^2}.
\end{equation}
Notice that this is a quadratic equation in $\lambda$. We can compute its minimum through derivation. One finds that $\lambda_{min} = -\frac{\mean{i[\hat{A'},\hat{B'}]}}{2\left(\Delta\hat{B'}\right)^2}$. Notice that it corresponds indeed to a minimum, for the coefficient in front of $\lambda^2$ is positive in \eqref{mean value O} (which is basically a parabola).
We now apply our findings to $\hat{A}=\hat{X}$ and $\hat{B}=\hat{P}$. The minimum is then $\frac{\h}{2\left(\Delta\hat{P}\right)^2}$. In particular, one has that $\mean{\mathcal{O}(\lambda)} = 0$ if and only if
\begin{align}
\left(\hat{X'}+i\lambda_{min}\hat{P'}\right)\ket{\psi} &= 0\\
\left(X-\mean{X}\right)\ket{\psi} &= -\frac{i\h}{2\left(\Delta\hat{P}\right)^2}\left(\hat{P}-\mean{\hat{P}}\right)\ket{\psi}
\end{align}
Let us now write $l = \Delta\hat{X}$. Because we define a coherent state as the saturation of the Heisenberg uncertainty principle, we can write $2\left(\Delta\hat{P}\right)^2 = \frac{\h^2}{l^2}$. Hence, one has that
\begin{equation}
\left(\hat{X}-\mean{\hat{X}}\right)\ket{\psi} = -i\frac{l^2}{\h}\left(\hat{P}-\mean{\hat{P}}\right)\ket{\psi}\tag{2}\label{b wave function}
\end{equation}
Using the x-representation expression of $\hat{P}$ (see answers above for more details on this), we find the differential equation
\begin{align}
-x\psi(x)+\expval{\hat{X}}\psi(x)&=\frac{il^2}{\h}\left(\frac{\h}{i}\psi'(x)-\expval{\hat{P}}\psi(x)\right)\notag\\
\psi'(x) &= \left[\frac{\mean{\hat{X}}-x}{l^2}+\frac{i}{\h}\mean{\hat{P}}\right]\psi(x)\tag{3}\label{wave function}
\end{align}
The solution to this is given by $\psi(x) = \psi_0\exp(-\frac{\left(x-\mean{\hat{X}}\right)^2}{2l^2}+\frac{i}{\h}\mean{\hat{P}}x)$ where $\psi_0$ is a normalization constant.
Through an appropriate Gaussian integral, one finds that the wave function of a coherent state is given by
\begin{equation}
\psi(x) = \frac{\exp(-\frac{i}{2\pi}\mean{\hat{X}}\mean{\hat{P}})}{\left(\pi l^2\right)^{1/4}}\exp\left(-\frac{\left(x-\mean{\hat{X}}\right)^2}{2l^2}+\frac{i\mean{\hat{P}}x}{\hbar}\right)\tag{4}\label{wave function coherent states}
\end{equation}
\eqref{wave function coherent states} shows that the sought after wave function is a gaussian.
Notice we added a global phase of $\exp(-\frac{i}{2\pi}\mean{\hat{X}}\mean{\hat{P}})$. Global phases are arbitrary and do not change the physics of the problem. An appropriate phase might be chosen in such a way as to simplify future considerations, and that is what is done here.
Interlude : other formulations of coherent states
I am going to develop here the mathematics behind coherent states, so as to find two other formulations of them that are easier to use. In particular, we're going to be able to use algebra instead of calculus - which is arguably (in my experience) simpler to manipulate. The reader only interested in the proposed alternative solution to the problem arised by the OP may skip to the next section.
Eigenvectors of the annihilation operator
Let us introduce $\hat{a} = \frac{1}{2}\left(\frac{\hat{X}}{l}+\frac{il\hat{P}}{\h}\right)$. One easily checks that $[\hat{a},\hat{a}^\dagger] = \mathbb{I}$. We call $\hat{a}$ the annihilation operator, whereas $\hat{a}^\dagger$ is known as the creation operator. The reasons behind these names are obvious when one studies the eigenvalues of $\hat{N} = \hat{a}^\dagger\hat{a}$.
Let $\ket{\psi}$ be an element of an Hilbert space $\mathcal{H}$. What are the solutions of the eigen-value problem $\hat{a}\ket{\psi} = \alpha\ket{\psi}$? Using the defintion of $\hat{a}$ provided, one has that $\frac{1}{2}\left(\frac{\hat{X}}{l}+\frac{il\hat{P}}{\h}\right)\ket{\psi} = \alpha\ket{\psi}$. One steadily finds that $\alpha = \frac{1}{2}\left(\frac{\mean{\hat{X}}}{l}+\frac{il\mean{\hat{P}}}{\h}\right)$. Putting this in the eigen-value problem, one finds (after a few appropriate manipulations) the equation \eqref{b wave function} again.
This means that the wave function defined by \eqref{wave function coherent states} corresponds to the eigenvectors of $\hat{a}$. Per convention, we write $\ket{\psi} = \ket{\alpha}$.
Theorem 3 $\ket{\alpha}$ is a coherent state (as defined in Definition 2) if and only if $\ket{\alpha}$ verifies
\begin{equation}
\hat{a}\ket{\alpha}=\alpha\ket{\alpha}\tag{5}\label{coherent in eigen}
\end{equation}
for any $\alpha\in\mathbb{C}$.
Displacement operator
We define the displacement operator as $\hat{D}(\alpha) = e^{\alpha\hat{a}^\dagger-\alpha^*\hat{a}}$.
Property 4 (Glauber) Let $\hat{A}$ and $\hat{B}$ be two operators. The relation $e^\hat{A}e^\hat{B} = e^{\hat{A}+\hat{B}}e^{-\frac{1}{2}[\hat{A},\hat{B}]}$ holds true if and only if ($\hat{A},\hat{B}$) commutes with [$\hat{A},\hat{B}$].
Through an appropriate use of Property 4, one shows the displacement operator defined above can equivalently be defined as $\hat{D}(\alpha) = e^{-\frac{1}{2}\mathbb{I}}e^{\alpha\hat{a}^\dagger}e^{-\alpha^*\hat{a}}$.
Property 5 One can show that a coherent state can be defined as the Displacement operator applied to the ground state $\ket{0}$:
\begin{equation}
\ket{\alpha} = \hat{D}(\alpha)\ket{0}\tag{6}\label{eq:coherent in displacement}
\end{equation}
Here you go! Three equivalent formulations of coherent states. Of course, this was a very quick treatment and the reader interested in further developments in the topic may consult other sources. I suggest reading Prof. Reinhold's lecture notes on the matter.
Alternative solution to the question
Without any further mathematical developments or proofs, let us bring to the reader's attention (I may edit the post, later on, to add some maths to back up the following claim) that coherent states closely imitates the classical solutions of the harmonic oscillator. That means that if applied to the Hamiltonian of an harmonic oscillator, the means values of the position and impulsion operators follow the classical solution of the Harmonic Oscillator. The reader may read more about this here.
The OP did not specify, but I reckon he used the ground state in the Fock space $\ket{n}$. Although coherent states and Fock states are not the same, it turns out their ground states are equivalent. Moreover, one may prove the following relation
\begin{equation}
\ket{\alpha} = e^{-\frac{1}{2}\norm{\alpha}^2}\sum_n \frac{\alpha^n}{\sqrt{n!}}\ket{n}
\end{equation}
which would allow one to define a coherent state as a linear combination of Fock states. Using this, one may see that the ground state of a quantum harmonic oscillator is $\braket{x}{0} = \psi_0(x)$. Taking $\ket{0}$ to be the coherent void, one only has to put $\alpha = 0$ in \eqref{wave function coherent states}. This gives us an alternative way of computing of the ground state of the quantum harmonic oscillator!
PS: I don't have time to read everything I wrote now. There is probably lots of grammatical mistakes, and I hope as little as possible mathematical ones. Feel free to edit if it feels necessary.
Best Answer
You could have equally well used $$ \hat{a}^\dagger= \frac{1}{\sqrt{2}}\int\!\!d\xi ~ |\xi\rangle \left(\xi-\frac{\partial}{\partial\xi}\right)\langle \xi|, $$ instead, since the hermitian conjugate of your starting expression is $$ \langle 0| \hat{a}^\dagger=0. $$ You then have $$ 0=\langle 0| \hat{a}^\dagger = \frac{1}{\sqrt{2}}\int\!\!d\xi ~ \langle 0|\xi\rangle \left(\xi-\frac{\partial}{\partial\xi}\right)\langle \xi| \\= \frac{1}{\sqrt{2}}\int\!\!d\xi ~ \left(\xi \psi^*_0(\xi) +\frac{\partial \psi^*_0(\xi) }{\partial\xi}\right)\langle \xi| ~~, $$ the last step involving integration by parts and use of the definition $\langle 0|\xi\rangle\equiv \psi^*_0(\xi) $.
Consequently, you get the same equation you had before, $$ \xi \psi^*_0(\xi) +\frac{\partial \psi^*_0(\xi) }{\partial\xi}=0, $$ with real solution, $$ \psi^*_0(\xi) \propto e^{-\xi^2/2} ~, $$ so you may complex conjugate and suitably normalize, etc.
It should be evident that the crucial step connecting states to functions is identical.