Let the eigenvectors of the phase estimation be $|\phi_i\rangle$. There are lots of these (an entire orthonormal basis). You know that what the phase estimation achieves is
$$
|\phi_i\rangle \longrightarrow |\phi_i\rangle|\theta_i\rangle
$$
where $\theta_i$ is the vector that contains the information about the eigenvalue.
So, if your input state is $|b\rangle$, you can always decompose this in terms of the eigenvector basis, i.e. there exist some numbers $\alpha_i$ such that
$$
|b\rangle=\sum_i\alpha_i|\phi_i\rangle.
$$
Now, everything that we are doing is linear, so if we know the action on the $|\phi_i\rangle$, we know the action on $|b\rangle$:
$$
|b\rangle\longrightarrow\sum_i\alpha_i|\phi_i\rangle|\theta_i\rangle.
$$
About $|\Psi_0\rangle$ (at least a little bit of intuition): Usually when you perform phase estimation for something like Shor's algorithm, your initial state is $|000\ldots 0\rangle$ and is first acted upon by $H$ on every qubit. So, really, we could write the state at that point as the input, and it would be
$$
|\Psi_{\text{original}}\rangle=\frac{1}{\sqrt{T}}\sum_{\tau=0}^{T-1}|\tau\rangle
$$
where $\tau$ is just the decimal representation of a $t$ bit number, if you're doing the phase estimation with a register of $t$ bits. This choice of initial state does a particular job -- whatever the eigenvalue that's being measured, if it isn't exactly one of the phases that's detected by the phase estimation, you are still guaranteed that with high probability you get the best approximation to the eigenvalue. In HHL, this is not exactly what you care about. For eigenvalues $\lambda$, you're going to be computing $\frac{1}{\lambda}$ in the inversion step. So we care far more about making sure we get the small $\lambda$ terms accurate than we do the large $\lambda$ terms if we're trying to keep the overall error down. Thus, the idea of using a slightly different initial state is to try and minimise this error.
