In 3-space, one can interpret the 4 Maxwell equation as determining the relationship between the fields (the electric field vector and the magnetic field bivector) and all four types of possible sources.
But this is rather illusory. In relativity, the equations look quite different:
$$\begin{align*} \nabla \cdot F &= -\mu_0 J \\ \nabla \wedge F &= 0\end{align*}$$
where $F$ is the electromagnetic field bivector. The vector derivative $\nabla$ can only increase or decrease the grade of an object by 1. Since $F$ is grade 2, the divergence equation describes its relationship with a grade 1 source term (the vector four-current $J$). The curl equation describes its relationship to a grade 3 (trivector) source term (of which there is none).
The reason the 4 Maxwell equations in 3-space come out the way we do is that we ignore the timelike basis vector, which would unify the scalar charge density with the 3-current as the four-current, as well as unify the E field with the B field as a bivector. The relativistic formulation, however, is considerably more sensible, as it correctly presents the EM field as one object of a single grade (a bivector), which can only have two sources (vector or trivector). It just so happens that the EM field has no trivector source.
What if there were trivector sources? Well, as you observe, there would be magnetic charge density (monopoles), but there would also be quite a bit more. There would have to be magnetic current as well, which would add an extra term to the $\nabla \times E$ equation to fully symmetrize things.
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
Roughly speaking, yes; changing the signs of $\mu_0$ or $\epsilon_0$ (or both) can lead to "run-away" solutions.
The easiest way to see why changing the sign of $\mu_0$ and/or $\epsilon_0$ is problematic is to look at the potential energy stored in a set of electromagnetic fields: $$ U = \frac{1}{2} \iiint \left(\epsilon_0 E^2 + \frac{1}{\mu_0} B^2 \right) \, \mathrm{d} \tau $$ In the "real-world" case, both $\epsilon_0$ and $\mu_0$ are positive. This means that there is a bound on how big $|\vec{E}|$ and $|\vec{B}|$ can get (averaged over space.) There's only a certain amount of energy to go around in the system, and if I want to make a really really large magnetic field field, I am constrained by the fact that I only have a total amount of energy $U$ available to do it.
But if $\epsilon_0$ and $\mu_0$ are of different signs, then this constraint disappears. Let's suppose that $\epsilon_0 < 0$ and $\mu_0 > 0$. If that's the case, then I can make a very very strong magnetic field without violating conservation of energy, so long as I also make a large electric field! Since magnetic and electric fields would contribute to the integral with different signs, I could get arbitrarily large fields while still having the above integral being relatively small.
In fancy mathematical terms, we say that the above energy integral (with $\epsilon_0 \mu_0 < 0$) is "indefinite", since it can be either positive or negative. Models with indefinite energy often have this problem of run-away solutions. The "real-world" integral, with $\mu_0, \epsilon_0 > 0$, has a "positive definite" energy: the integral is always positive except in the case where both fields vanish, in which case the integral is zero. The positive-definiteness of the integral is important to placing bounds on the fields in the "real-world" case.
Finally, flipping both of $\epsilon_0$ and $\mu_0$ does not lead to runaway solutions in "vacuum" (i.e., due solely to the field energy.) In this case, the above integral is "negative definite" rather than "positive definite". But these fields can also (for example) do work on electric charges, giving them mechanical kinetic energy; and the kinetic energy of these charges is always positive definite. So in this case, there still exists the possibility of runaway solutions in which charged objects get very large kinetic energies (with a positive energy of great magnitude) and the fields become arbitrarily large (with a negative energy of great magnitude.)