This is actually exactly the same question I was asking myself a while ago, and it took me quite some time to figure it out. What I ended up doing was:
- diagonalise the BdG Hamiltonian in Mathematica
- solve the expressions for the eigenvalues for $\kappa$
- neglect terms proportional to $\kappa^2$, $\kappa \Delta$ in the expressions for the eigenvectors, expand them to first order in $\kappa$, insert the expression for $\kappa$ you got from the eigenvalues
- simplify everything (using the expression for $\epsilon$ and approximate $\hbar v k \approx E_F+U$) and arrive at beenakkers simple looking expressions (tedious)
I hope that helps some, in case you haven't figured it out yourself in the meantime. If you need more info/details don't hesitate to ask, I can elaborate some more if needed!
