Let be
$$\frac{2a}{Q}V(\theta,\varphi)=f(\theta,\varphi)=2\sin\theta\cos\varphi+\cos^2\theta.\tag 1$$
The Laplace spherical harmonics form a complete set of orthonormal functions and thus form an orthonormal basis of the Hilbert space of square-integrable functions. On the unit sphere, any square-integrable function can thus be expanded as a linear combination of these:
$$
f(\theta,\varphi)=\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell f_\ell^m \, Y_\ell^m(\theta,\varphi)\tag 2
$$
where $Y_\ell^m( \theta , \varphi )$ are the Laplace spherical harmonics defined as
$$
Y_\ell^m( \theta , \varphi ) = \sqrt{{(2\ell+1)\over 4\pi}{(\ell-m)!\over (\ell+m)!}} \, P_\ell^m ( \cos{\theta} ) \operatorname{e}^{i m \varphi } =N_{\ell}^m P_\ell^m ( \cos{\theta} ) \operatorname{e}^{i m \varphi }\tag 3
$$
and where $N_{\ell}^m$ denotes the normalization constant
$ N_{\ell}^m \equiv \sqrt{{(2\ell+1)\over 4\pi}{(\ell-m)!\over (\ell+m)!}},$ and $P_\ell^n(\cos\theta)$ are the associated Legendre polynomials.
The Laplace spherical harmonics are orthonormal
$$
\int_{\theta=0}^\pi\int_{\varphi=0}^{2\pi}Y_\ell^m \, Y_{\ell'}^{m'*} \, d\Omega=\delta_{\ell\ell'}\, \delta_{mm'},
$$
where $δ_{ij}$ is the Kronecker delta and $\operatorname{d}\Omega = \sin\theta \operatorname{d}\varphi\operatorname{d}\theta$.
The expansion coefficients are the analogs of Fourier coefficients, and can be obtained by multiplying the above equation by the complex conjugate of a spherical harmonic, integrating over the solid angle $Ω$, and utilizing the orthogonality relationships. This is justified rigorously by basic Hilbert space theory. For the case of orthonormalized harmonics, this gives:
$$
f_\ell^m=\int_{\Omega} f(\theta,\varphi)\, Y_\ell^{m*}(\theta,\varphi)\operatorname{d}\Omega = \int_0^{2\pi}\operatorname{d}\varphi\int_0^\pi \operatorname{d}\theta\,\sin\theta f(\theta,\varphi)Y_\ell^{m*} (\theta,\varphi). \tag 4
$$
where $ Y_\ell^{m*} (\theta, \varphi) = (-1)^m Y_\ell^{-m} (\theta, \varphi)$.
The evaluation of the expansion $f_\ell^m$ may be very long in this way...
We can use some tricks in your case for $f(\theta,\varphi)=2\sin\theta\cos\varphi+\cos^2\theta$ observing that
$$
\sin\theta=-P_1^1(\cos\theta)
$$
and
$$
Y_{1}^{-1}(\theta,\varphi) - Y_{1}^{1}(\theta,\varphi)= {1\over 2}\sqrt{3\over 2\pi}\cdot e^{-i\varphi}\cdot\sin\theta-{-1\over 2}\sqrt{3\over 2\pi}\cdot e^{i\varphi}\cdot\sin\theta =\sqrt{3\over 2\pi} \sin\theta\cos\phi
$$
so that
$$
2\sin\theta\cos\varphi=2\sqrt{\frac{2\pi}{3}}\left(Y_1^{-1}(\theta,\varphi)-Y_1^{1}(\theta,\varphi)\right)=-2P_1^1(\cos\theta)\cos\varphi\tag 5
$$
and observing that
$$
\cos^2\theta=\frac{1}{3}P_0^0+\frac{2}{3}P_2^0
$$
and using the relation $Y_\ell^0(\theta,\varphi)=\sqrt{\frac{2\ell+1}{4\pi}}P_\ell^0(\cos\theta)$ where $P_\ell^0(\cos\theta)$ are the ordinary Legendre's polynomials $P_\ell(\cos\theta)$, we obtain
$$
\cos^2\theta=\frac{1}{3}P_0^0(\cos\theta)+\frac{2}{3}P_2^0(\cos\theta)=2\sqrt{\pi}Y_0^0(\theta,\varphi)+\frac{4}{3}\sqrt{\frac{\pi}{5}}Y_2^0(\theta,\varphi).\tag 6
$$
Finally, putting together (5) and (6) in (1) we obtain
$$
f(\theta,\varphi)=2\sqrt{\frac{2\pi}{3}}Y_1^{-1}(\theta,\varphi)-2\sqrt{\frac{2\pi}{3}}Y_1^{1}(\theta,\varphi)+2\sqrt{\pi}Y_0^0(\theta,\varphi)+\frac{4}{3}\sqrt{\frac{\pi}{5}}Y_2^0(\theta,\varphi)\tag 7
$$
so that, comparing (7) and (2), the coefficients $f_\ell^m$ are
$$
f_1^{-1}=2\sqrt{\frac{2\pi}{3}}\qquad
f_1^{1}=-2\sqrt{\frac{2\pi}{3}}\qquad
f_0^{0}=2\sqrt{\pi}\qquad
f_2^{0}=\frac{4}{3}\sqrt{\frac{\pi}{5}}\tag 8
$$
So you have
$$\small
V(\theta,\varphi)=\frac{Q}{a}\left[\sqrt{\frac{2\pi}{3}}Y_1^{-1}(\theta,\varphi)-\sqrt{\frac{2\pi}{3}}Y_1^{1}(\theta,\varphi)+\sqrt{\pi}Y_0^0(\theta,\varphi)+\frac{2}{3}\sqrt{\frac{\pi}{5}}Y_2^0(\theta,\varphi)\right]
$$
The general solution to the Laplace equation outside the sphere is
$$
\Phi(r,\theta,\varphi)=\sum_{\ell=0}^\infty \sum_{m=-\ell}^\ell \frac{B_\ell^m}{r^{\ell+1}} \, Y_\ell^m(\theta,\varphi)\tag 9
$$
with $B_\ell^m=\frac{Q}{a}f_\ell^m$,
that is
$$\small
\Phi(r,\theta,\varphi)=\frac{Q}{a}\left[\frac{1}{r^2}\left(\sqrt{\frac{2\pi}{3}}Y_1^{-1}(\theta,\varphi)-\sqrt{\frac{2\pi}{3}}Y_1^{1}(\theta,\varphi)\right)+\frac{\sqrt{\pi}}{r}Y_0^0(\theta,\varphi)+\frac{1}{r^3}\left(\frac{2}{3}\sqrt{\frac{\pi}{5}}\right)Y_2^0(\theta,\varphi)\right]
$$
I think the handiest trick to dealing with these situations is by noting that, if you have a solution $\varphi$ to the Poisson equation $\nabla^2 \varphi = -\rho/\varepsilon_0$, then you can get new ones by applying the cartesian derivatives $\partial/\partial x$ and so on. (This is because $\partial/\partial x$ commutes with the laplacian; the trick fails in curvilinear coordinates where that's no longer the case.)
Many standard solutions can be seen as arising in this way: for example, the point-dipole field is obtained by taking the point-charge field,
$$
\varphi(\mathbf r) = \frac{q}{4\pi\varepsilon_0}\frac{1}{r},
\qquad\text{obeying}\qquad
\nabla^2\varphi = q\delta(\mathbf r)/\varepsilon_0,
$$
and differentiating it with respect to $x$, giving you
$$
\partial_x\varphi(\mathbf r) = -\frac{q}{4\pi\varepsilon_0}\frac{x}{r^3},
\qquad\text{obeying}\qquad
\nabla^2\partial_x\varphi = q\partial_x\delta(\mathbf r)/\varepsilon_0,
$$
where $\partial_x\delta(\mathbf r) = \delta'(x)\delta(y)\delta(z)$ is a point dipole's representation as a distribution-valued charge density. (The higher multipoles can also be obtained in this way, though if you just use cartesian derivatives then you'll need to be careful with how you handle their linear combinations.)
Similarly, the two-displaced-spheres calculation has a natural home in this formalism, where you take as the starting point the field of a uniform ball of charge:
$$
\varphi(\mathbf r) = \begin{cases}\frac{Q}{4\pi\varepsilon_0}\frac{1}{r} & r>a \\ \frac{Q}{4\pi\varepsilon_0} \frac{r^2}{a^3} & r<a\end{cases},
\qquad\text{obeying}\qquad
\nabla^2\varphi = \rho_0\theta(a-r)/\varepsilon_0,
$$
where you now need to differentiate the Heaviside step function to give you
$$
\frac{\partial}{\partial x} \theta(a-r) = \frac{\partial r}{\partial x} \theta'(a-r) = \frac{x}{r} \delta(a-r),
$$
and the internal derivative $\partial_x r^2 = 2x$ gives you a linear potential and therefore a constant electric field. (I'm possibly fudging some constants. Work this stuff through in detail to check it's right.)
Moreover, it should be clear why this derivative trick works for this situation: if you want to model a uniform ball of charge centered at $(\Delta x,0,0)$, you can say that it's the result of a uniform ball of charge at the origin plus the $\sigma\propto\cos(\theta)$ double-shell thing, or you can do a Taylor expansion in $\Delta x$ to get
$$
\varphi_{\Delta x}(\mathbf r) = \varphi_{0}(\mathbf r) + \Delta x \left[ \frac{\partial}{\partial \Delta x}\varphi_{\Delta x}(\mathbf r) \right]_{\Delta x=0} + \mathcal O((\Delta x)^2),
$$
where the derivative along $\Delta x$ is equivalent to a derivative along $x$ (in the same way that active and passive frame transformations are equivalent) and the neglected terms in $\mathcal O((\Delta x)^2)$ are also thrown away in the hand-wavy example.
OK, so that's all nice, but how do you apply it to your cylinder? Well, you start with a uniform cylinder of charge, and you take the derivative of the charge distribution and its solutions with respect to a coordinate orthogonal to the axis. Happy calculating!
Best Answer
Consider first a charged sphere of radius $a$ with a uniform density $\rho=\sigma_0/d$. Since $\sigma_0$ is surface density, one has to divide by a length to get a volume density. Write the potentiel outside the sphere ($r>a$): $$\varphi(r)={4\over 3}\pi a^3\rho {1\over 4\pi\varepsilon_0r}$$ Now consider two spheres with uniform charge densities $\rho=\pm \sigma_0/d$ respectively whose centers are located at $\vec r_\pm=\pm {d\over 2}\vec u_z$. Write the total potential $$\varphi(r)={a^3\sigma_0/d\over 3\varepsilon_0} \left[{1\over ||\vec r-d/2\vec u_z||}-{1\over ||\vec r+d/2\vec u_z||} \right]$$ Perform a Taylor expansion to lowest order (same calculation as the dipole).
The last step is to convince yourself that the two spheres are equivalent to a single sphere with a surface density $\sigma_0\cos\theta$ (draw a figure for example). When you substract the two spheres, you end up with a thin spherical layer whose thickness is $d\cos\theta$. Consider a surface element $dS$ on this layer. The charge carried is $\rho dV=\sigma_0/d\times dS\times d\cos\theta=\sigma_0\cos\theta dS$. Therefore, it can be interpreted as a sphere carrying a surface density $\sigma=\sigma_0\cos\theta$.