Warning: not an expert, so could be major mistakes in this.
The multiplicity is one for every element in the global packet, in the non-CM case. In the CM case,
half of the packet has multiplicity one and the other half has multiplicity zero.
For the general story, I'd suggest looking at Arthur's conjectures, which at least conjecturally give a very nice picture. In general one should talk of Arthur packets, not Langlands packets. They coincide in your case.
For CM representations on $SL_2$, the associated global parameter has dihedral image inside $PGL_2(\mathbb{C})$. Its centralizer $S$ has size $2$. According to Arthur, the obstruction to an element of the global packet being automorphic is valued in the dual of $S$; thus, "one-half" is automorphic. More precisely, the set of irreducible representations in a global $L$-packet is a principal homogeneous space for a certain product of $Z/2Z$s (in the obvious way: one $Z/2Z$ for each $p$ where $a_p = 0$; for simplicity, suppose that local multiplicity $4$ doesn't occur). The automorphic ones correspond exactly to one of the fibers of the summation map to Z/2Z. Thus, if you take an automorphic representation in this packet, and switch it at one place to the other element of the local packet, it won't be automorphic any more.
Concretely: Start with $\Pi$ an automorphic representation for $GL_2$, it maps
by restriction to the space of automorphic forms for $SL_2$. The above remarks
suggest that this restriction map must have a huge kernel when $\Pi$ is CM. But one can almost see this by hand: $\Pi$ is isomorphic
to its twist by a certain quadratic character $\omega$. This, and multiplicity one for $GL_2$,
means that, for $f \in \Pi$, the function $g \mapsto f(g) \omega(\det g)$ also belongs to $\Pi$. This forces extra vanishing, although I didn't work out the details: For instance, if $\omega$ is everywhere unramified and $f$ the spherical vector, then $f(g) \omega (\det g)$ must be proportional to $f$, so $f$ vanishes whenever $\omega(\det g)$ is not $1$, i.e. various translates of $f$ do vanish when restricted to $SL_2$.
Remark/warning: If you are concerned with multiplicity one, there is a further totally distinct phenomenon that causes it to fail for $\mathrm{SL}_n, n \geq 3$: Two non-conjugate homomorphisms of a finite group into $\mathrm{PGL}_n$ can be conjugate element by element.
There's a paper of Blasius about this. It has nothing to do with packets.
Excellent question indeed. The quick answer is that $E_2(z)$ is an almost holomorphic modular form of weight $2$ and level $1$, so the automorphic representation generated by it is not irreducible. For more details (and my thought process), read below.
Consider the Maass raising operator in weight $0$,
$$ R:=y\left(i\frac{\partial}{\partial x}+\frac{\partial}{\partial y}\right).$$
Let $(m,n)\in\mathbb{Z}^2$ be a nonzero pair of integers. Then a small calculation gives that, for $z=x+iy$ and $s\in\mathbb{C}$,
$$ R\left(\frac{y^s}{|mz+n|^{2s}}\right) =\frac{sy^s}{(mz+n)^2|mz+n|^{2s-2}}.$$
Now let us introduce the usual weight $0$ level $1$ (nonholomorphic) Eisenstein series
$$ E(z,s):=\sum_{\substack{m, n \in \mathbb{Z} \\ (m, n) \ne (0,0)}} \frac{\operatorname{Im}(z)^s}{|mz + n|^{2s}},\qquad \operatorname{Re}(s)>1$$
then we see that
$$ R\,E(z,s+1) = (s+1)\,yE_2(z,s),\qquad \operatorname{Re}(s)>0.\tag{1}$$
On the right hand side, $yE_2(z,s)$ is the canonical weight $2$ level $1$ (nonholomorphic) Eisenstein series, the one which transforms as a weight $2$ Maass form. It is worthwhile to recall here that weight $k$ holomorphic forms embed into the weight $k$ Maass spectrum by multiplying each weight $k$ holomorphic form by $y^{k/2}$. In our case $k=2$, which explains why we multiply by $y$.
So your Eisenstein series, after inserting the factor $y$ to make it into a canonical weight $2$ form, and also inserting the scaling factor $s+1$, equals the Maass raising shift of $E(z,s+1)$. It belongs to the same automorphic representation as $E(z,s+1)$, hence it has the same Langlands parameters as $E(z,s+1)$ at every place. In particular, the archimedean Langlands parameters are
$$ (s+1)-\frac{1}{2}=s+\frac{1}{2}\qquad\text{and}\qquad \frac{1}{2}-(s+1)=-s-\frac{1}{2}.$$
Added and revised. Well, we still need to specify all this to $s=0$, but in this case the above argument breaks down, because $E(z,s+1)$ has a pole at $s=0$. So we need to be more careful. Let us use the results and notation of Section 4 of Duke-Friedlander-Iwaniec: The subconvexity problem for Artin L-functions (Invent. Math. 149 (2002), 489-577). Then for $\operatorname{Re}(s)>1$ we have the Fourier decomposition
\begin{align*} \frac{1}{2}E(z,s)&=\ \zeta(2s)y^s+\pi^{2s-1}\frac{\Gamma(1-s)}{\Gamma(s)}\zeta(2-2s)y^{1-s}\\&+\ \frac{\pi^s}{\Gamma(s)}\sum_{n=1}^\infty\frac{\sigma_{2s-1}(n)}{n^s}\bigl\{f_0^+(nz,s)+f_0^-(nz,s)\bigr\}.\end{align*}
Let us replace $s$ by $s+1$ here, and then apply the raising operator $R$ along with $(1)$. Then for $\operatorname{Re}(s)>0$ we obtain the Fourier decomposition
\begin{align*} \frac{1}{2}E_2(z,s)&=\ \zeta(2s+2)y^s+\pi^{2s+1}\frac{\Gamma(1-s)}{\Gamma(2+s)}\zeta(-2s)y^{-s-1}\\&-\ \frac{\pi^{s+1}}{y\Gamma(2+s)}\sum_{n=1}^\infty\frac{\sigma_{2s+1}(n)}{n^{s+1}}\bigl\{f_2^+(nz,s+1)+s(s+1)f_2^-(nz,s+1)\bigr\}.\end{align*}
The right hand side is indeed holomorphic at $s=0$, and at this value it specifies to
\begin{align} \frac{1}{2}E_2(z)&=\ \frac{\pi^2}{6}-\frac{\pi}{2y}- \frac{\pi}{y}\sum_{n=1}^\infty\frac{\sigma_1(n)}{n}f_2^+(nz,1)\\
&=\ \frac{\pi^2}{6}-\frac{\pi}{2y}-4\pi^2\sum_{n=1}^\infty\sigma_1(n)e(nz).\tag{2}\end{align}
It is clear now that the $L$-function of $E_2(z)$ is $\zeta(s-1/2)\zeta(s+1/2)$, and $E_2(z)$ should belong to the holomorphic discrete series of weight $2$ even though its constant term is not holomorphic. I think this paradox arises from the fact that $E_2(z)$ is not a true automorphic form. (More precisely, $E_2(z)$ is not a vector from an irreducible automorphic representation, see the Added 3 section below.)
Added 2. Indeed, $E_2(z)$ is an almost holomorphic modular form of weight $2$ and level $1$: it equals $2G_2^*(z)$ in the notation of Section 2.3 of Bruinier-v.d.Geer-Harder-Zagier's book "The 1-2-3 of modular forms". In particular, (19) and (21) in this book reveal that $E_2(z)$ transforms precisely like a holomorphic modular form of weight $2$ and level $1$, even though it is not holomorphic. Of course, the same also follows from the fact that $E_2(z)=\lim_{s->0+}E_2(z,s)$, where $yE_2(z,s)$ for $\operatorname{Re}(s)>0$ transforms precisely like a Maass form of weight $2$ and level $1$. One can read more about almost holomorphic modular forms in Section 5.3 of the book.
Added 3. We can learn more by applying the Maass lowering operator in weight $2$,
$$L:=1+y\left(i\frac{\partial}{\partial x}-\frac{\partial}{\partial y}\right).$$
Using the Fourier decomposition (2), we can see directly that
$$ L(yE_2(z))=-\pi, $$
which harmonizes with the facts that $ L(yE_2(z,s)) = -s E(z,s+1)$ for $\operatorname{Re}(s)>0$ and $\operatorname{res}_{s=0}E(z,s+1)=\pi$. So we see that the automorphic representation generated by the weight $2$ vector $yE_2(z)$ is reducible: it contains the trivial representation $\mathbb{C}$ as a subrepresentation, while its quotient by $\mathbb{C}$ is irreducible and belongs to the holomorphic discrete series of weight $2$. See also Theorem 2.5.2 in Bump: Automorphic forms and representations, specified to $k=2$ and $\lambda=0$, which helped me understand what is going on.
Best Answer
Jared Weinstein and I wrote a paper on how to compute $\pi_p$: see here.
As Olivier says, $a_p$ will often be zero, and in fact if the central character is trivial (or has conductor coprime to $p$) this is always the case when $p^2$ divides the level of $f$. One can get a bit futher by twisting: you can always twist a newform by Dirichlet characters, and Atkin and Li have shown that $\pi_p$ is principal series or Steinberg at $p$ if and only if there is some Dirichlet character $\chi$ such that the twist of $f$ by $\chi$ is a newform with nonzero Fourier coefficient at $p$ (or an oldform attached to such a newform).
So that leaves the supercuspidal cases, and here Hecke theory won't help you at all: no matter how you twist your form, the Hecke eigenvalues are all zero. One can actually show (the "local converse theorem") that $\pi_p$ is uniquely determined by the Atkin-Lehner pseudo-eigenvalues of all of the twists of $f$; but it is not so easy to calculate these, or to explicitly identify $\pi_p$ from a list of them once you've done so. In my paper with Jared you can find details of a different approach, using Bushnell and Kutzko's theory of "types", which seems to work quite well.
These algorithms are implemented in recent versions of Magma (and should be in Sage fairly shortly as well, once someone gets around to reviewing my code).