No explicit complexification is needed to derive this breakdown of Maxwell's equations. This can be understood wholly through the real vector space of special relativity.
Let's start with Maxwell's equations for the EM field, in the clifford algebra language called STA: the spacetime algebra. Maxwell's equations take the form
$$\nabla F = -J$$
where $\nabla F = \nabla \cdot F + \nabla \wedge F$, $F = e_0 E + B \epsilon_3$, in the $(-, +, +, +)$ sign convention.
Let $x$ be the spacetime position vector. It's generally true that, for a vector $v$ and a constant bivector $C$,
$$\nabla (C \cdot x) = -2 C, \quad \nabla (C \wedge x) = 2C \implies \nabla (Cx) = 0$$
One can then evaluate the expression
$$\nabla (Fx) = (\nabla F)x + \dot \nabla (F \dot x)$$
where the overdot means that only $x$ is differentiated in the second term; using the product rule, $F$ is "held constant" and so the above formulas apply. We just argued that the second term is zero, so we get $\nabla (Fx) = (\nabla F) x$. Thus, we arrive at the following transformation of Maxwell's equations:
$$\nabla (Fx) = -Jx$$
Now, we could always write $F$ as a "complex bivector" in the sense that, using $\epsilon = e_0 \epsilon_3$, and $\epsilon \epsilon = -1$, we have
$$F = e_0 E - B \epsilon_3 e_0 \epsilon_3 \epsilon = e_0 (E + \epsilon B)$$
It's crucial to note that $\epsilon$ does not commute with any vector.
What are the components of $Fx$? Write $x = t e_0 + r$ and we can write them as
$$Fx = e_0 (Ex + \epsilon Bx) = e_0 (E \cdot r + E \wedge r - e_0 Et + \epsilon B \cdot r - e_0 B \times r + \epsilon B t e_0)$$
This too can be written in a "complex" form:
$$Fx = (e_0 E \cdot r + Et + B \times r) + \epsilon (E \times r + e_0 B\cdot r + Bt)$$
We seem to differ on some signs, but this is recognizably the same quantity you have called $G$.
Now, to talk about how these equations break down, let's write $G = G_1 + G_3$, where $G_1 = (e_0 E \cdot r + \ldots)$ and $G_3 = \epsilon (E \times r + \ldots)$. Let's also write for $R = Jx = R_0 + R_2$.
Maxwell's equations then become
$$\nabla \cdot G_1= R_0, \quad \nabla \wedge G_1 + \nabla \cdot G_3 = R_2, \quad \nabla \wedge G_3 = 0$$
The first and third equations are the components of the Gauss dipole; the second equation is the Ampere-Faraday dipole equation.
Now, what does it all mean? The expression for $G = Fx$ includes both rotational moments of the EM field as well as some dot products, so it measures both how much the spacetime position is in the same plane as the EM field as well as how much the spacetime position is out of the plane.
It's probably more instructive to look at the source term $-Jx$. This tells us both about the moments of the four-current as well as how it goes toward or away from the coordinate origin. The description for the moments is wholly in the Ampere-Faraday dipole equation. What kinds of moments would this describe? A pair of two opposite point charges at rest, separated by a spatial vector $2 \hat v$ and centered on the origin, each with current at rest $j_0$, would create a $R = Jx = + j_0 e_t \hat v - j_0 e_t (-\hat v) = 2 j_0 e_t \hat v$, so this would be described wholly by the A-F dipole equation.
That's at time zero, however. At later times, $R$ will pick up these weird time terms. Say we're at time $\tau$. Then $R = 2 j_0 e_t \hat v + j_0 e_t (\tau e_t) - j_0 e_t (\tau e_t)$. So for this case, there's no problem: the extra stuff will just cancel. A single charge, however, would start picking up this term.
In a few words, these equations are weird.
The Biot-Savart law says that under magnetostatic conditions ($\frac{\partial}{\partial t}\rightarrow 0$),
$$\mathbf B(\mathbf r) = \frac{\mu_0}{4\pi}\int \frac{\mathbf J(\mathbf r') \times (\mathbf r - \mathbf r')}{|\mathbf r - \mathbf r'|^3} dV'$$
Noting that
$$ \frac{\mathbf r - \mathbf r'}{|\mathbf r - \mathbf r'|^3} = \nabla \left(\frac{1}{|\mathbf r - \mathbf r'|}\right)$$
where $\nabla$ refers to differentiation by the unprimed coordinates, this can be written
$$\mathbf B(\mathbf r) = \nabla \times\frac{\mu_0}{4\pi} \int \frac{\mathbf J(\mathbf r')}{|\mathbf r - \mathbf r'|}dV'$$
Taking the curl of this and using the fact that $\nabla \times (\nabla \times \mathbf F) = \nabla(\nabla \cdot \mathbf F) - \nabla^2 \mathbf F$,
$$\nabla \times \mathbf B(\mathbf r) = \nabla\left(\int J(\mathbf r')\cdot \nabla\left[\frac{1}{|\mathbf r - \mathbf r'|}\right]dV'\right) - \nabla^2 \frac{\mu_0}{4\pi} \int \frac{\mathbf J(\mathbf r')}{|\mathbf r - \mathbf r'|}dV'$$
Noting that
$$\nabla\left[\frac{1}{|\mathbf r - \mathbf r'|}\right] = -\nabla'\left[\frac{1}{|\mathbf r - \mathbf r'|}\right]$$
we can integrate the first term by parts to obtain
$$\nabla\left(\int \nabla' \cdot \left[\frac{\mathbf J(\mathbf r')}{|\mathbf r - \mathbf r'|}\right]dV' - \int\frac{\nabla' \cdot \mathbf J(\mathbf r')}{|\mathbf r - \mathbf r'|} dV'\right)$$
The first term is a surface term, and vanishes if we assume that $\mathbf J(\mathbf r') \rightarrow 0$ as $|\mathbf r'| \rightarrow \infty$. The second term vanishes because according to the continuity equation, $\nabla\cdot\mathbf J = -\frac{\partial \rho}{\partial t} = 0$ in magnetostatics. This leaves us with
$$\nabla \times \mathbf B(\mathbf r) = \nabla^2 \frac{\mu_0}{4\pi} \int \frac{\mathbf J(\mathbf r')}{|\mathbf r - \mathbf r'|}dV'$$
and since
$$\nabla^2 \frac{1}{|\mathbf r - \mathbf r'|} = 4\pi \delta^{(3)}(\mathbf r - \mathbf r')$$
we have
$$\nabla \times \mathbf B = \mu_0 \mathbf J$$.
Again, Biot-Savart is valid only under magnetostatic conditions, and therefore so is this version of Ampere's law. It would be nice to relax these conditions and re-do this derivation more generally, but we don't yet know what to replace Biot-Savart with.
Instead, let's see how this version of Ampere's law fails when we go to general electrodynamics. Clearly since $\mathbf J \propto \nabla \times \mathbf B$, we have that $\nabla \cdot \mathbf J = 0$. However, according to the general continuity equation, $\nabla \cdot \mathbf J = -\frac{\partial \rho}{\partial t}$.
To fix this, let's assume that we need a new term, so
$$\nabla \times \mathbf B = \mu_0 \mathbf J + \mathbf G$$
for some vector field $\mathbf G$. Taking the divergence of both sides yields
$$ 0 = - \mu_0 \frac{\partial \rho}{\partial t} + \nabla \cdot \mathbf G$$
$$\implies \nabla \cdot \mathbf G = \mu_0 \frac{\partial \rho}{\partial t}$$
From Gauss' law for electric fields, we know that $\rho = \epsilon_0 \nabla \cdot \mathbf E$, and so
$$\nabla \cdot \mathbf G = \epsilon_0 \mu_0 \frac{\partial}{\partial t} \nabla \cdot \mathbf E = \nabla \cdot \left(\epsilon_0 \mu_0 \frac{\partial}{\partial t}\mathbf E\right)$$
and so we can simply postulate that
$$\mathbf G = \epsilon_0\mu_0 \frac{\partial}{\partial t} \mathbf E$$
so
$$\nabla \times \mathbf B = \mu_0 \mathbf J + \epsilon_0\mu_0 \frac{\partial}{\partial t} \mathbf E$$
This was Maxwell's correction to Ampere's law, and it has been validated over and over by experiment.
In summary, magnetostatics + Biot-Savart gives us $\nabla \times \mathbf B = \mu_0 \mathbf J$. Predictably, this fails when we leave the domain of magnetostatics, and in particular is inconsistent with the continuity equation. We don't know how to generalize Biot-Savart, but patching up the problem with the continuity equation in the simplest possible way yields the correct Ampere's law, $\nabla \times \mathbf B = \mu_0 \mathbf J + \epsilon_0 \mu_0 \frac{\partial}{\partial t}\mathbf E$.
From this, we can work backward to find the correct generalization of Biot-Savart; this is one of Jefimenko's equations.
EDIT:
Returning to the original derivation after eliminating the surface term (but before sending $\nabla'\cdot \mathbf J(\mathbf r')\rightarrow 0$), we have
$$\nabla \times \mathbf B(\mathbf r) = \mu_0 \mathbf J(\mathbf r) - \frac{\mu_0}{4\pi}\int \frac{\nabla '\cdot \mathbf J(\mathbf r')}{|\mathbf r - \mathbf r'|}$$
Under the conditions for which Biot-Savart holds, the latter term is equal to zero; however, we can be bold and throw those restrictions aside just to see what happens. Under general conditions, $\nabla \cdot \mathbf J = -\frac{\partial \rho}{\partial t}$, so that term becomes
$$\frac{\mu_0}{4\pi} \nabla \int \frac{\partial \rho}{\partial t} \frac{1}{|\mathbf r - \mathbf r'|} dV'= \frac{\partial}{\partial t} \frac{\mu_0}{4\pi} \nabla \int\frac{\rho(\mathbf r')}{|\mathbf r - \mathbf r'|}$$
Defining $$ \phi(\mathbf r) = \int \frac{\rho(\mathbf r')}{4\pi \epsilon_0 |\mathbf r-\mathbf r'|}$$
and letting $\mathbf E = -\nabla\phi$, this becomes
$$\nabla \times \mathbf B = \mu_0 \mathbf J + \epsilon_0 \mu_0 \frac{\partial}{\partial t} \mathbf E$$
What we did here - simply disregarding the conditions under which Biot-Savart is applicable and plugging in the more general continuity equation - is morally the same as Maxwell's addition of the extra term to compensate for the nonzero divergence of $\mathbf J$.
Note also that we have glossed over how to go from magnetostatics to electrodynamics $\big(\mathbf J(\mathbf r) \rightarrow \mathbf J(\mathbf r,t), \rho(\mathbf r)\rightarrow \mathbf \rho(\mathbf r,t)\big)$. Simply plugging in a $t$ to Biot-Savart and letting it "go along for the ride" is insufficient; working backwards from the full Maxwell's equations demonstrates the need to introduce the retarded time $t_r = t - \frac{|\mathbf r - \mathbf r'|}{c}$, indicating that Biot-Savart is genuinely wrong for electrodynamics.
Best Answer
Starting from this equation for the components: $$ \frac{dP_i}{dt}=\sum q_jE_i+\epsilon_0 c^2\int \left((\nabla \times B)\times B\right)_i d^3r-\epsilon_0 \int \left(E\times (\nabla \times E)\right)_i d^3r $$
Use the identity: $$ \epsilon_{ijk}\epsilon_{ilm} = \delta_{jl}\delta_{km} - \delta_{jm}\delta_{kl} $$ to rewrite the starting equation $$ \frac{dP_i}{dt}=\sum q_jE_i+\epsilon_0 c^2\int \left(\epsilon_{ijk}(\nabla_l B_m)\epsilon_{jlm}B_k\right) d^3r - \epsilon_0\int \left(\epsilon_{ijk}E_j(\nabla_l E_m)\epsilon_{klm}\right) d^3r $$ as: $$ \frac{dP_i}{dt}=\sum q_jE_i+\epsilon_0 c^2\int \left((\delta_{kl}\delta_{im} - \delta_{km}\delta_{il})(\nabla_l B_m)B_k\right) d^3r-\epsilon_0 \int \left((\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl})E_j(\nabla_l E_m)\right) d^3r $$
Now, we can use the fact the the fields vanish at infinity to integrate by parts and we have: $$ \frac{dP_i}{dt}=\sum q_jE_i-\epsilon_0 c^2\int \left((\delta_{kl}\delta_{im} - \delta_{km}\delta_{il})B_m(\nabla_l B_k)\right) d^3r+\epsilon_0 \int \left((\delta_{il}\delta_{jm}-\delta_{im}\delta_{jl})E_m(\nabla_l E_j)\right) d^3r $$
Then, Maxwell's equations give us: $$ \frac{dP_i}{dt}=\epsilon_0 c^2\int \delta_{km}\delta_{il}B_m\nabla_l B_k d^3r+\epsilon_0 \int \delta_{il}\delta_{jm}E_m\nabla_l E_j d^3r $$
Or: $$ \frac{dP_i}{dt}=\epsilon_0 c^2\int B_m\nabla_i B_m d^3r+\epsilon_0 \int E_m\nabla_i E_m d^3r\;, $$ where the repeated indices (m) are summed over.
But it turns out that each of the two terms on the RHS are individually equal to zero. This is true because, as we have already assumed, the fields vanish at infinity. So we can integrate by parts to our heart's content.
For example: $$ \int B_m(\nabla_i B_m) d^3r = -\int (\nabla_iB_m)(B_m) d^3r\;, $$ which says that this integral is equal to the negative of itself--that is, zero! (You can also see that it is zero by the reasoning of the other answer, namely, that it is half the total derivative of the square of the field.)
And similarly for the electric field integral.
Therefore the time derivative of the total momentum is zero (total momentum is conserved).