The most important statement in this answer to your question is: Yes, you can superimpose a constant magnetic field. The combined field remains a solution of Maxwell's equations. $\def\vB{{\vec{B}}}$
$\def\vBp{{\vec{B}}_{\rm p}}$
$\def\vBq{{\vec{B}}_{\rm h}}$
$\def\vE{{\vec{E}}}$
$\def\vr{{\vec{r}}}$
$\def\vk{{\vec{k}}}$
$\def\om{\omega}$
$\def\rot{\operatorname{rot}}$
$\def\grad{\operatorname{grad}}$
$\def\div{\operatorname{div}}$
$\def\l{\left}\def\r{\right}$
$\def\pd{\partial}$
$\def\eps{\varepsilon}$
$\def\ph{\varphi}$
Since you are using plane waves you even cannot enforce the fields to decay sufficiently fast with growing distance to the origin. That would make the solution of Maxwell's equations unique for given space properties (like $\mu,\varepsilon,\kappa$, and maybe space charge $\rho$ and an imprinted current density $\vec{J}$). But, in your case you would not have a generator for the field. Your setup is just the empty space. If you enforce the field to decay sufficiently fast with growing distance you just get zero amplitudes $\vec{E}_0=\vec{0}$, $\vec{B}_0=\vec{0}$ for your waves. Which is certainly a solution of Maxwell's equations but also certainly not what you want to have.
For my point of view you are a bit too fast with the integration constants. You loose some generality by neglecting that these constants can really depend on the space coordinates.
Let us look what really can be deduced for $\vB(\vr,t)$ from Maxwell's equations for a given
$\vE(\vr,t)=\vE_0 \cos(\vk\vr-\om t)$ in free space.
At first some recapitulation:
We calculate a particular B-field $\vBp$ that satisfies Maxwell's equations:
$$
\begin{array}{rl}
\nabla\times\l(\vE_0\cos(\vk\vr-\om t)\r)&=-\pd_t \vBp(\vr,t)\\
\l(\nabla\cos(\vk\vr-\om t)\r)\times\vE_0&=-\pd_t\vBp(\vr,t)\\
-\vk\times\vE_0\sin(\vk\vr-\om t) = -\pd_t \vBp(\vr,t)
\end{array}
$$
This leads us with $\pd_t \cos(\vk\vr-\om t) = \om \sin(\vk\vr-\om t)$ to the ansatz
$$
\vBp(\vr,t) = -\vk\times\vE_0 \cos(\vk\vr-\om t)/\om.
$$
The divergence equation $\div\vBp(\vr,t)=-\vk\cdot(\vk\times\vE_0)\cos(\vk\vr-\om t)/\om=0$ is satisfied and the space-charge freeness
$0=\div\vE(\vr,t) = \vk\cdot\vE_0\sin(\vk\vr-\om t)$ delivers that $\vk$ and $\vE_0$ are orthogonal. The last thing to check is Ampere's law
$$
\begin{array}{rl}
\rot\vBp&=\mu_0 \eps_0 \pd_t\vE\\
\vk\times(\vk\times\vE_0)\sin(\vk\vr-\om t)/\om &= -\mu_0\eps_0 \vE_0 \sin(\vk\vr-\om t) \om\\
\biggl(\vk \underbrace{(\vk\cdot\vE_0)}_0-\vE_0\vk^2\biggr)\sin(\vk\vr-\om t)/\om&= -\mu_0\eps_0 \vE_0 \sin(\vk\vr-\om t) \om
\end{array}
$$
which is satisfied for $\frac{\omega}{|\vk|} = \frac1{\sqrt{\mu_0\eps_0}}=c_0$ (the speed of light).
Now, we look which modifications $\vB(\vr,t)=\vBp(\vr,t)+\vBq(\vr,t)$ satisfy Maxwell's laws.
$$
\begin{array}{rl}
\nabla\times\vE(\vr,t) &= -\pd_t\l(\vBp(\vr,t)+\vBq(\vr,t)\r)\\
\nabla\times\vE(\vr,t) &= -\pd_t\vBp(\vr,t)-\pd_t\vBq(\vr,t)\\
0 &= -\pd_t\vBq(\vr,t)
\end{array}
$$
That means, the modification $\vBq$ is independent of time. We just write $\vBq(\vr)$ instead of $\vBq(\vr,t)$.
The divergence equation for the modified B-field is $0=\div\l(\vBp(\vr,t)+\vBq(\vr)\r)=\underbrace{\div\l(\vBp(\vr,t)\r)}_{=0} + \div\l(\vBq(\vr)\r)$ telling us that the modification $\vBq(\vr)$ must also be source free:
$$
\div\vBq(\vr) = 0
$$
Ampere's law is
$$
\begin{array}{rl}
\nabla\times(\vBp(\vr,t)+\vBq(\vr)) &= \mu_0\eps_0\pd_t \vE,\\
\rot(\vBq(\vr))&=0.
\end{array}
$$
Free space is simply path connected. Thus, $\rot(\vBq(\vr))=0$ implies that every admissible $\vBq$ can be represented as gradient of a scalar potential $\vBq(\vr)=-\grad\ph(\vr)$.
From $\div\vBq(\vr) = 0$ there follows that this potential must satisfy Laplace's equation
$$
0=-\div(\vBq(\vr)) = \div\grad\ph = \Delta\ph
$$
That is all what Maxwell's equations for the free space tell us with a predefined E-field and without boundary conditions:
The B-field can be modified through the gradient of any harmonic potential.
The thing is that with problems in infinite space one is often approximating some configuration with finite extent which is sufficiently far away from stuff that could influence the measurement significantly.
How are plane electromagnetic waves produced?
One relatively simple generator for electromagnetic waves is a dipole antenna. These do not generate plane waves but spherical curved waves as shown in the following nice picture from the Wikipedia page http://en.wikipedia.org/wiki/Antenna_%28radio%29.
Nevertheless, if you are far away from the sender dipol and there are no reflecting surfaces around you then in your close neighborhood the electromagnetic wave will look like a plane wave and you can treat it as such with sufficiently exact results for your practical purpose.
In this important application the plane wave is an approximation where the superposition with some constant electromagnetic field is not really appropriate.
We just keep in mind if in some special application we need to superimpose a constant field we are allowed to do it.
Best Answer
Two reasons:
Maxwell's equations are linear. Also the operations ${\rm Re}$, ${\rm Im}$ and $z\mapsto z^*$ are linear in the sense that a sum mapped by the respective operator is the sum of the mapping of the addends by that operator.
With a time-varying sinusoidal quantity the mapping between the quantities $|a| \cos(\omega\,t+\arg(a));\,a\in\mathbb{C}, \omega\,t\in\mathbb{R}$ (or $i\,|a| \,\sin(\omega\,t+\arg(a))$) on the one hand and $a\,e^{-i\,\omega\,t};\,a\in\mathbb{C}$ on the other is one-to-one and onto. So for every entity of the form $a\,e^{-i\,\omega\,t};\,a\in\mathbb{C}$ there is a unique $|a| \cos(\omega\,t+\arg(a))$ and contrawise. Explicitly:
$$|a| \cos(\omega\,t+\arg(a)) = {\rm Re}(a\,e^{-i\,\omega\,t})$$
and, because $\omega\,t$ takes on every value in $[0,2\pi]$, one can uniquely infer $\arg(a)$, $|a|$ and $\omega$ from the values of $f(t) = |a| \cos(\omega\,t+\arg(a))$ as $t$ varies. Likewise for the inversion of ${\rm Im}$.
It's important to keep in mind that this seeming "trickery" works because $\omega\,t$ varies, so we see the variation of the real and imaginary parts with time. In contrast, taking the real or imaginary part of a lone complex is of course irreversible: the imaginary (or real) part is lost and one cannot get the original complex number back from only its real (or imaginary) part!
You can think of the above pithily as: the phasor signal is the "single sideband" (if you recall this archaic modulation scheme) version of the real valued signal. That is, the negative frequency component of a real valued signals is simply the complex conjugate of the positive frequency component. So, for a linear calculation, there is no need to "process" both positive and negative components: one calculates with the positive frequency component and then recovers the real signal by taking the complex conjugate of the calculation outcome, thus finding the negative frequency component, and adding it back to the "single-sideband" positive frequency component.
Phasors are simply a tool of convenience: calculations with $e^{-i\,\omega\,t}$ are easier than those with $\cos$ and $\sin$. However, in the special case of Maxwell's equations, one can interpret the complex quantities as more than simply phasors (although the technique turns out to be the same). See my answer here where I show that the complex quantities are intimately linked to the unique decomposition of the electromagnetic field into its left and right hand circularly polarised components.