Planck, BICEP, et al are all detecting electromagnetic radiation, but the "E-modes" and "B-modes" refer to polarization characteristics of this radiation, not the actual electric and magnetic fields. As you surmised, the names derive from an analogy to the decomposition of a vector field into curl-less (here "E" for electric or "G" for gradient) and divergence-less ("B" for magnetic or "C" for curl) components, as follows...
The first step is the measurement of the standard Stokes parameters $Q$ and $U$. In general, the polarization of monochromatic light is completely described via four Stokes parameters, which form a (non-orthonormal) vector space when the various waves are incoherent. For light propagating in the $z$ direction, with electric field:
$$ E_x = a_x(t) \cos(\omega_0 t - \theta_x (t)) \, \, , \quad E_y = a_y(t) \cos(\omega_0 t - \theta_y (t)) $$
the Stokes parameters are:
- $ I = \langle a_x^2 \rangle + \langle a_y^2 \rangle $ , intensity
- $ Q = \langle a_x^2 \rangle - \langle a_y^2 \rangle $ , polarization along $x$ (Q>0) or $y$ (Q<0) axes
- $ U = \langle 2 a_x a_y \cos(\theta_x - \theta_y) \rangle $ , polarization at $\pm 45$ degrees
- $ V = \langle 2 a_x a_y \sin(\theta_x - \theta_y) \rangle $ , left- or right-hand circular polarization
In cosmology, no circular polarization is expected, so $V$ is not considered. In addition, normalization of $Q$ and $U$ is traditionally with respect to the mean temperature $T_0$ instead of intensity $I$.
The definitions of $Q$ and $U$ imply that they transform under a rotation $\alpha$ around the $z$-axis according to:
$$ Q' = Q \cos (2 \alpha) + U \sin (2 \alpha) $$
$$ U' = -Q \sin (2 \alpha) + U \cos (2 \alpha) $$
These parameters transform, not like a vector, but like a two-dimensional, second rank symmetric trace-free (STF) polarization tensor $\mathcal{P}_{ab}$. In spherical polar coordinates $(\theta, \phi)$, the metric tensor $g$ and polarization tensor are:
$$ g_{ab} = \left( \begin{array}{cc} 1 & 0 \\ 0 & \sin^2 \theta \end{array} \right) $$
$$ \mathcal{P}_{ab}(\mathbf{\hat{n}}) =\frac{1}{2}
\left( \begin{array}{cc} Q(\mathbf{\hat{n}}) & -U(\mathbf{\hat{n}}) \sin \theta \\
-U(\mathbf{\hat{n}}) \sin \theta & -Q(\mathbf{\hat{n}})\sin^2 \theta \end{array} \right) $$
As advertised, this matrix is symmetric and trace-free (recall the trace is $g^{ab} \mathcal{P}_{ab}$).
Now, just as a scalar function can be expanded in terms of spherical harmonics $Y_{lm}(\mathbf{\hat{n}})$, the polarization tensor (with its two independent parameters $Q$ and $U$) can be expanded in terms of two sets of orthonormal tensor harmonics:
$$ \frac{\mathcal{P}_{ab}(\mathbf{\hat{n}})}{T_0} = \sum_{l=2}^{\infty} \sum_{m=-l}^{l} \left[ a_{(lm)}^G Y_{(lm)ab}^G(\mathbf{\hat{n}}) +
a_{(lm)}^C Y_{(lm)ab}^C(\mathbf{\hat{n}}) \right]$$
where it turns out that:
$$ Y_{(lm)ab}^G = N_l \left( Y_{(lm):ab} - \frac{1}{2} g_{ab} {Y_{(lm):c}}^c\right) $$
$$ Y_{(lm)ab}^C = \frac{N_l}{2} \left( Y_{(lm):ac} {\epsilon^c}_b + Y_{(lm):bc} {\epsilon^c}_a \right)$$
where $\epsilon_{ab}$ is the completely antisymmetric tensor, "$:$" denotes covariant differentiation on the 2-sphere, and
$$ N_l = \sqrt{\frac{2(l-2)!}{(l+2)!}} $$
The "G" ("E") basis tensors are "like" gradients, and the "C" ("B") like curls.
It appears that cosmological perturbations are either scalar (e.g. energy density perturbations) or tensor (gravitational waves). Crucially, scalar perturbations produce only E-mode (G-type) polarization, so evidence of a cosmological B-mode is (Nobel-worthy) evidence of gravitational waves. (Note, however, that Milky-Way "dust" polarization (the "foreground" to cosmologists) can produce B-modes, so it must be well-understood and subtracted to obtain the cosmological signal.)
An excellent reference is Kamionkowski. See also Hu.