$\newcommand{\bra}[1]{\left\langle#1\right|}\newcommand{\ket}[1]{\left|#1\right\rangle}\newcommand{\proj}[1]{|#1\rangle\langle#1|}\newcommand{\half}{\frac12}$In answer to your first question, I wrote myself some notes some time ago about my understanding of how it worked. The notation is probably a bit different (I've tried to bring it more into line, but it's easy to miss bits), but attempts to explain that choice of the state $|\Psi_0\rangle$. There also seem to be some factors of $\frac12$ floating around in places.
When we first study phase estimation, we're usually thinking about it in respect to use in some particular algorithm, such as Shor's algorithm. This has a specific goal: getting the best $t$-bit approximation to the eigenvalue. You either do, or you don't, and the description of phase estimation is specifically tuned to give as high a success probability as possible.
In HHL, we are trying to produce some state
$$
\ket{\phi}=\sum_j\frac{\beta_j}{\lambda_j}\ket{\lambda_j},
$$
where $\ket{b}=\sum_j\beta_j\ket{\lambda_j}$, making use of phase estimation. The accuracy of the approximation of this will depend far more critically on an accurate estimation of the eigenvalues that are close to 0 rather than those that are far from 0. An obvious step therefore, is to attempt to modify the phase estimation protocol so that rather than using `bins' of fixed width $2\pi/T$ for approximating the phases of $e^{-iAt}$ ($T=2^t$ and $t$ is number of qubits in phase estimation register), we might rather specify a set of $\phi_y$ for $y\in\{0,1\}^t$ to act as the centre of each bin so that we can have vastly increased accuracy close to 0 phase. More generally, you might specify a trade-off function for how tolerant you might be of errors as a function of the phase $\phi$. The precise nature of this function can then be tuned to a given application, and the particular figure of merit which you will use to determine success. In the case of Shor's algorithm, our figure of merit was simply this binning protocol -- we were successful if the answer was in the correct bin, and unsuccessful outside it. This is not going to be the case in HHL, whose success is more reasonably captured by a continuous measure such as the fidelity. So, for the general case, we shall designate a cost function $C(\phi,\phi')$ which specifies a penalty for answers $\phi'$ if the true phase is $\phi$.
Recall that the standard phase estimation protocol worked by producing an input state that was the uniform superposition of all basis states $\ket{x}$ for $x\in\{0,1\}^t$. This state was used to control the sequential application of multiple controlled-$U$ gates, which are followed up by an inverse Fourier transform. Imagine we could replace the input state with some other state
$$
\ket{\Psi_0}=\sum_{x\in\{0,1\}^t}\alpha_x\ket{x},
$$
and then the rest of the protocol could work as before. For now, we will ignore the question of how hard it is to produce the new state $\ket{\Psi_0}$, as we are just trying to convey the basic concept. Starting from this state, the use of the controlled-$U$ gates (targeting an eigenvector of $U$ of eigenvalue $\phi$), produces the state
$$
\sum_{x\in\{0,1\}^t}\alpha_xe^{i\phi x}\ket{x}.
$$
Applying the inverse Fourier transform yields
$$
\frac{1}{\sqrt{T}}\sum_{x,y\in\{0,1\}^t}e^{ix\left(\phi-\frac{2\pi y}{M}\right)}\alpha_x\ket{y}.
$$
The probability of getting an answer $y$ (i.e. $\phi'=2\pi y/T$) is
$$
\frac{1}{T}\left|\sum_{x\in\{0,1\}^t}e^{ix\left(\phi-\frac{2\pi y}{T}\right)}\alpha_x\right|^2
$$
so the expected value of the cost function, assuming a random distribution of the $\phi$, is
$$
\bar C=\frac{1}{2\pi T}\int_0^{2\pi}d\phi\sum_{y\in\{0,1\}^t}\left|\sum_{x\in\{0,1\}^t}e^{ix\left(\phi-\frac{2\pi y}{T}\right)}\alpha_x\right|^2C(\phi,2\pi y/T),
$$
and our task is to select the amplitudes $\alpha_x$ that minimise this for any specific realisation of $C(\phi,\phi')$. If we make the simplifying assumption that $C(\phi,\phi')$ is only a function of $\phi-\phi'$, then we can make a change of variable in the integration to give
$$
\bar C=\frac{1}{2\pi}\int_0^{2\pi}d\phi\left|\sum_{x\in\{0,1\}^t}e^{ix\phi}\alpha_x\right|^2C(\phi),
$$
As we noted, the most useful measure is likely to be a fidelity measure. Consider we have a state $\ket{+}$ and we wish to implement the unitary $U_\phi=\proj{0}+e^{i\phi}\proj{1}$, but instead we implement $U_{\phi'}=\proj{0}+e^{i\phi'}\proj{1}$. The fidelity measures how well this achieves the desired task,
$$
F=\left|\bra{+}U_{\phi'}^\dagger U\ket{+}\right|^2=\cos^2\left(\frac{\phi-\phi'}{2}\right),
$$
so we take
$$
C(\phi-\phi')=\sin^2\left(\frac{\phi-\phi'}{2}\right),
$$
since in the ideal case $F=1$, so the error, which is what we want to minimise, can be taken as $1-F$.
This will certainly be the correct function for evaluating any $U^t$, but for the more general task of modifying the amplitudes, not just the phases, the effects of inaccuracies propagate through the protocol in a less trivial manner, so it is difficult to prove optimality, although the function $C(\phi-\phi')$ will already provide some improvement over the uniform superposition of states. Proceeding with this form, we have
$$
\bar C=\frac{1}{2\pi}\int_0^{2\pi}d\phi\left|\sum_{x\in\{0,1\}^t}e^{ix\phi}\alpha_x\right|^2\sin^2\left(\half\phi\right),
$$
The integral over $\phi$ can now be performed, so we want to minimise the function
$$
\half\sum_{x,y=0}^{T-1}\alpha_x\alpha_y^\star(\delta_{x,y}-\half\delta_{x,y-1}-\half\delta_{x,y+1}).
$$
This can be succinctly expressed as
$$
\min\bra{\Psi_0}H\ket{\Psi_0}
$$
where
$$
H=\half\sum_{x,y=0}^{T-1}(\delta_{x,y}-\half\delta_{x,y-1}-\half\delta_{x,y+1})\ket{x}\bra{y}.
$$
The optimal choice of $\ket{\Psi_0}$ is the minimum eigenvector of the matrix $H$,
$$
\alpha_x=\sqrt{\frac{2}{T+1}}\sin\left(\frac{(x+1)\pi}{T+1}\right),
$$
and $\bar C$ is the minimum eigenvalue
$$
\bar C=\half-\half\cos\left(\frac{\pi}{T+1}\right).
$$
Crucially, for large $T$, $\bar C$ scales as $1/T^2$ rather than the $1/T$ that we would have got from the uniform coupling choice $\alpha_x=1/\sqrt{T}$. This yields a significant benefit for the error analysis.
If you want to get the same $|\Psi_0\rangle$ as reported in the HHL paper, I believe you have to add the terms $-\frac14\left(\ket{0}\bra{T-1}+\ket{T-1}\bra{0}\right)$ to the Hamiltonian. I have no justification for doing so, however, but this is probably my failing.
