The Analogy Is With Non-Abelian Gauge Fields
The relevant analogy is with non-Abelian gauge fields. Using the conventions I wrote out in my reply here
Yang-Mills vs Einstein-Hilbert Action
the relevant analogy is with the non-linear equations:
$$dA + A A = F, \hspace 1em dF + AF - FA = 0.$$
For electromagnetism and gauge fields, one has
$$A = A^a Y_a = A_μ dx^μ = A^a_μ Y_a dx^μ, \hspace 1em
F = F^c Y_c = ½ F_{μν} dx^μ dx^ν = ½ F^c_{μν} Y_c dx^μ dx^ν,$$
where I'm using juxtaposition $dx^μ dx^ν$ for wedge products $dx^μ∧dx^ν$ and adopting the convention of allowing the gauge basis $Y_a$ to freely intermix with the form, e.g. $Y_a dx^μ = dx^μ Y_a$.
For the Cartan structure equations
$$dθ^a + ω^a_c θ^c = Θ^a, \hspace 1em dω^a_b + ω^a_c ω^c_b = Ω^a_b,$$
one has
$$A = θ^a Y_a + ω^a_b Y_a^b, \hspace 1em F = Θ^a Y_a + Ω^a_b Y_a^b,$$
with the frame one-forms $θ^a = h^a_μ dx^μ$, connection one-forms $ω^a_b = Γ^a_{μb} dx^μ$, torsion two-forms $Θ^a = ½ T^a_{μν} dx^μ dx^ν$ and curvature two-forms $Ω^a_b = ½ R^a_{bμν} dx^μ dx^ν$.
This plays the role of $dA + A A = F$, for a suitably-defined set of Lie brackets defined for the basis elements $Y_a$ and $Y_a^b$. (The Lie brackets are those for the Lie algebra of the gauge group $GA(4)$ - the 4D general affine group).
The role of the "Bianchi identity" equations $dA + A F - FA = 0$ is then played by the equations derived from the Cartan structure equations:
$$dΘ^a + ω^a_c Θ^c = Ω^a_c θ^c, \hspace 1em dΩ^a_b + ω^a_c Ω^c_b = Ω^a_c ω^c_b.$$
Because this is subject to the constraint that the (constant) frame metric $g_{ab}$ has zero-covariant derivative,
$$∇g_{ab} ≡ dg_{ab} - ω^c_a g_{cb} - ω^c_b g_{ac} = 0,$$
then along with the equations derived from it:
$$d∇g_{ab} - ω^c_a ∇g_{cb} - ω^c_b ∇g_{ac} = Ω^c_a g_{cb} - Ω^c_b g_{ac},$$
this entails that the connection one-form and curvature two-form are anti-symmetric in their frame indices, when raised or lowered by the metric. The preferred position that dovetails into the above-mentioned reply is with them up:
$$ω^a_c g^{cb} = ω^{ab} = -ω^{ba}, \hspace 1em Ω^a_c g^{cb} = Ω^{ab} = -Ω^{ba}.$$
In that case, the constraint can be removed, reducing the one-form $A$ and two-form $F$ to:
$$A = θ^a P_a + ½ ω^{ab} S_{ab}, \hspace 1em F = Θ^a P_a + ½ Ω^{ab} S_{ab},$$
the metric $g_{ab}$ now being woven into the structure coefficients for the basis elements $P_a$ and $S_{ab}$. (The relevant gauge group is the inhomogeneous extension $ISO(3,1)$ of the Lorentz group $SO(3,1)$.)
The actual reduction - with more detail on the correspondence laid out here - is provided in my earlier reply here
What are the analogues of $F_{\mu\nu}$ in General Relativity?
The space-time metric, itself, is the frame metric converted to space-time indices: $g_{μν} = h^a_μ h^b_ν g_{ab}$.
In my later reply, cited on top above, which starts out with this reduced form, I made a minor change to the $[P,P]$ Lie brackets to produce this as the field strength
$$F = Θ^aP_a + ½\left(Ω^{ab} + λθ^aθ^b\right)S_{ab}$$
so as to directly weave in the cosmological coefficient, in the process also enabling the action to be expressed equivalently as a quadratic function of the gauge field strengths. Depending on what the coefficient $λ$ used with the modification is, the Lie algebra can be that for either $SO(3,2)$ or $SO(4,1)$.
What Are The "Second Set" Of Fields?
First, you should identify what actually are the "second set" of fields; and that is obscured by the form you wrote the equations in. I will rewrite your equations, plus the continuity equation (which you left out), in the following form:
$$
∇·𝐞 = 4πr, \hspace 1em ∇×𝐛 = \frac{4π}{c}𝐣 + \frac{1}{c}\frac{∂𝐞}{∂t}
\hspace 1em⇒\hspace 1em
∇·𝐣 + \frac{1}{c}\frac{∂r}{∂t} = 0, \\
∇×𝐞 = -\frac{1}{c}\frac{∂𝐛}{∂t}, \hspace 1em ∇·𝐛 = 0 \hspace 1em⇐\hspace 1em
∇·𝐞 = -\frac{1}{c}\frac{∂𝐚}{∂t} - ∇f, \hspace 1em 𝐛 = ∇×𝐚.
$$
The first set of equations arise from an action principle given by the Lagrangian density
$$𝔏 = \frac{|𝐞|^2 - |𝐛|^2}{8π} + 𝔏_0(𝐚,f),$$
the additional "interaction" part $𝔏_0(𝐚,f)$ arises from all other fields that have non-trivial gauge transformations.
The extra fields arise as derivatives of the Lagrangian density, whose total variational can be written in terms of them as:
$$δ𝔏 = δ𝐚·𝐣 - δf r + δ𝐞·𝐝 - δ𝐛·𝐡,$$
with the Lagrangian density, itself, serving as a device to generate a set of constitutive relations that connect those fields with the gauge fields:
$$
𝐝 = \frac{∂𝔏}{∂𝐞} = \frac{𝐞}{4π}, \hspace 1em 𝐡 = -\frac{∂𝔏}{∂𝐛} = \frac{𝐛}{4π}, \hspace 1em
r = -\frac{∂𝔏}{∂f} = -\frac{∂𝔏_0}{∂f}, \hspace 1em 𝐣 = \frac{∂𝔏}{∂𝐚} = \frac{∂𝔏_0}{∂𝐚}.
$$
From this, your first set of equations arise as Euler-Lagrange equations
$$∇·𝐝 = r, \hspace 1em ∇×𝐡 = \frac{𝐣}{c} + \frac{1}{c}\frac{∂𝐝}{∂t},$$
for the response fields $(𝐝,𝐡)$ and sources $(r,𝐣)$.
This is best seen by restoring the constitutive coefficients, making the change:
$$(f,𝐚,𝐛,𝐞) = \sqrt{4πε_0} (φ,c𝐀,c𝐁,𝐄), \hspace 1em (r,𝐣,𝐝,𝐡) = \frac{(ρ,𝐉,𝐃,𝐇/c)}{\sqrt{4πε_0}},$$
thereby giving us back the Maxwell (and SI) form of the equations:
$$
F = dA \hspace 1em⇒\hspace 1em 𝐁 = ∇×𝐀, \hspace 1em 𝐄 = -∇φ - \frac{∂𝐀}{∂t}, \\
dF = 0 \hspace 1em⇒\hspace 1em ∇·𝐁 = 0, \hspace 1em ∇×𝐄 + \frac{∂𝐁}{∂t} = 𝟎,
$$
with the one-form $A$ and two-form $F$ being written
$$A = 𝐀·d𝐫 - φ dt, \hspace 1em F = 𝐁·d𝐒 + 𝐄·d𝐫 dt,$$
where $d𝐫 = (dx, dy, dz)$ and $d𝐒 = (dy dz, dz dx, dx dy)$. The non-Abelian analogue of this would be:
$$
𝐁 = ∇×𝐀 + 𝐀×𝐀, \hspace 1em 𝐄 = -∇φ - \frac{∂𝐀}{∂t} + φ𝐀 - 𝐀φ, \\
∇·𝐁 + 𝐀·𝐁 - 𝐁·𝐀 = 0, \hspace 1em ∇×𝐄 + \frac{∂𝐁}{∂t} + 𝐀×𝐄 + 𝐄×𝐀 + 𝐁φ - φ𝐁 = 𝟎.
$$
There are no $c$'s in any of this, because there is no connection with the Minkowski - or any - metric in these equations.
The second set of fields are the response fields and sources that are generated by the Lagrangian density, arising as its derivatives, with the variational being:
$$δ𝔏 = δ𝐀·𝐉 - δφ ρ + δ𝐄·𝐃 - δ𝐁·𝐇.$$
With the Lagrangian density now expressed as
$$𝔏 = \frac{ε_0|𝐄|^2}{2} - \frac{|𝐁|^2}{2μ_0} + 𝔏_0(φ,𝐀),$$
where $μ_0 = 1/(ε_0c^2)$ this yields the constitutive relations
$$𝐃 = \frac{∂𝔏}{∂𝐄} = ε_0 𝐄, \hspace 1em 𝐇 = -\frac{∂𝔏}{∂𝐁} = \frac{𝐁}{μ_0}, \hspace 1em ρ = -\frac{∂𝔏}{∂φ} = -\frac{∂𝔏_0}{∂φ}, \hspace 1em 𝐉 = \frac{∂𝔏}{∂𝐀} = \frac{∂𝔏_0}{∂𝐀},$$
along with the following as the Euler-Lagrange equations and (derived from it) the continuity equation
$$∇·𝐃 = ρ, \hspace 1em ∇×𝐇 - \frac{∂𝐃}{∂t} = 𝐉, \hspace 1em ∇·𝐉 + \frac{∂ρ}{∂t} = 0.$$
In terms of differential forms, the response fields and sources are:
$$G = ½ 𝔊^{μν} ∂_ν ˩ ∂_μ ˩ d^4x = 𝐃·d𝐒 - 𝐇·d𝐫 dt, \hspace 1em J = 𝔍^μ ∂_μ ˩ d^4x = ρ dV - 𝐉·d𝐒 dt,$$
where
$$dV = dx dy dz, \hspace 1em d^4 x = dt dV = dt dx dy dz.$$
The contraction operator $∂_μ ˩ α$ is defined recursively on differential forms by $∂_μ ˩ (dx^ν ∧ α) = δ_μ^ν α - dx^ν ∧ (∂_μ ˩ α)$ and $∂_μ ˩ f = 0$ for scalars (a.k.a. 0-forms) $f$. For instance,
$$
\frac{∂}{∂t} ˩ d^4 x = \frac{∂}{∂t} ˩ dt dV = dV, \\
\frac{∂}{∂x} ˩ d^4 x = \frac{∂}{∂x} ˩ dt dx dy dz = -dt \left(\frac{∂}{∂x} ˩ dx dy dz\right) = -dt dy dz = -dy dz dt.
$$
The non-Abelian version of these equations are:
$$
∇·𝐃 + 𝐀·𝐃 - 𝐃·𝐀 = ρ, \hspace 1em ∇×𝐇 - \frac{∂𝐃}{∂t} + 𝐀×𝐇 + 𝐇×𝐀 + φ𝐃 - 𝐃φ = 𝐉, \\
∇·𝐉 + \frac{∂ρ}{∂t} + 𝐀·𝐉 - 𝐉·𝐀 + ρφ - φρ = 𝐃·𝐄 - 𝐄·𝐃 + 𝐁·𝐇 - 𝐇·𝐁.
$$
The response fields and sources are now:
$$G = G_c y^c = ½ 𝔊_c^{μν} y^c ∂_ν ˩ ∂_μ ˩ d^4x, \hspace 1em J = J_c y^c = 𝔍^μ ∂_μ ˩ d^4x = 𝔍_c^μ y^c ∂_μ ˩ d^4x.$$
They are referenced to the dual Lie basis $y^a$ and are not even the same type of objects as the gauge fields. They Lie co-vectors, while the gauge fields are Lie vectors.
That's what the "second set" of fields actually are.
In terms of the Lagrangian 4-form $L = 𝔏 d^4x$, the variational may be written:
$$ΔL = \left(ΔA^a\right) J_a - \left(ΔF^a\right) G_a.$$
For Yang-Mills actions and Lagrangians, $FG - GF = 0$ and the terms $𝐃·𝐄 - 𝐄·𝐃 = 0$ and $𝐁·𝐇 - 𝐇·𝐁 = 0$ drop out from the continuity equation.
All of this is obscured by the form your equations were written in, which essentially equated $𝐝$ with $𝐞$ and $𝐛$ with $𝐡$.
Where Are The Second Set Of Fields?
The Einstein-Hilbert Lagrangian with a cosmological coefficient - as a 4-form - can be written, along with the corresponding action $S$, as:
$$S = \int L, \hspace 1em L = ε_{abcd} \sqrt{|g|} \left(k Ω^{ab} + l θ^aθ^b\right)θ^cθ^d,$$
for suitably-defined coefficients $k$ and $l$.
Technically, if the torsion is not constrained to be 0, it is the Einstein-Cartan action (with cosmological term), while it is the Einstein-Hilbert action (with cosmological term) expressed on a Riemann-Cartan geometry as a Palatini action, if the torsion is constrained to be 0. But with matter-free Lagrangians, or outside of matter of a matter term is included, the torsion reduces to 0 as a result of the field equations.
If $l ≠ 0$, the action can be equivalently written as
$$S = \int L', \hspace 1em L' = m ε_{abcd} \sqrt{|g|} \left(Ω^{ab} + λθ^aθ^b\right)\left(Ω^{cd} + λθ^cθ^d\right),$$
where
$$m = \frac{k^2}{4l}, \hspace 1em λ = \frac{2l}{k},$$
because $ε_{abcd} \sqrt{|g|} Ω^{ab} Ω^{cd}$ drops out from the action as a boundary term.
Up to proportionality $ε_{abcd} \sqrt{|g|} Ω^{ab} θ^cθ^d$ is just $R \sqrt{|g|} d^4x$, itself, where $R$ is the curvature scalar.
In both cases, you can write down the variational in terms of a set of response fields and sources:
$$ΔL = \left(Δθ^a\right) J_a + ½ \left(Δω^{ab}\right) J_{ab} - \left(ΔΘ^a\right) G_a - ½ \left(ΔΩ^{ab}\right) G_{ab}.$$
But the $G$'s and $J$'s are the bona fide response and source fields only for the field strength $F$ without the $λ$ adjustment.
With the $λ$ adjustment, the fields would be:
$$
ΔL = \left(Δθ^a\right) \bar{J}_a + ½ \left(Δω^{ab}\right) J_{ab} - \left(ΔΘ^a\right) G_a - ½ \left(Δ\bar{Ω}^{ab}\right) G_{ab},\\
\bar{Ω}^{ab} = Ω^{ab} + λθ^aθ^b, \hspace 1em \bar{J}_a = J_a + λG_{ab}θ^b.
$$
The 3-currents associated with the source fields
$$J_a = 𝔍^μ_a ∂_μ ˩ d^4 x, \hspace 1em J_{ab} = 𝔍^μ_{ab} ∂_μ ˩ d^4 x,$$
are naturally identified as the 3-currents, respectively, for momentum $p_a$ and internal angular momentum $s_{ab}$.
The corresponding Euler-Lagrange equations are
$$
dG_a - G_c ω^c_a = J_a, \\
dG_{ab} - G_{cb} ω^c_a - G_{ac} ω^c_b = J_{ab} + θ_a G_b - θ_b G_a,
$$
the resulting "continuity" equations being:
$$
dJ_a + J_c ω^c_a = -G_c Ω^c_a, \\
dJ_{ab} + J_{cb} ω^c_a + J_{ac} ω^c_b - θ_b J_a + θ_a J_b = -G_{ac} Ω^c_b - G_{cb} Ω^c_a - Θ_a G_b + Θ_b G_a.
$$
So ... what are the extra fields? With the unadjusted Lagrangian, they are:
$$
G_a = 0, \hspace 1em G_{ab} = 2k \sqrt{|g|} ε_{abcd} θ^c θ^d, \hspace 1em
J_c = \sqrt{|g|} ε_{abcd} \left(2k Ω^{ab} + 4 l θ^a θ^b\right) θ^d, \hspace 1em J_{cd} = 0.
$$
Up to a constant factor the curvature term $ε_{abcd} Ω^{ab} θ^d$ in the source $J_c$ is the Einstein tensor, itself, expressed as a 3-current density: $\sqrt{|g|} e^ν_c \left(R^μ_ν - ½ δ^μ_ν R\right) ∂_μ ˩ d^4 x$, where $e^ν_c$ is the inverse to the frame matrix $h^μ_a$.
The $λ$-adjusted response fields and sources, in terms of the $λ$-adjusted action and Lagrangian, take on a much simpler form:
$$G_a = 0, \hspace 1em G_{ab} = 2m \sqrt{|g|} ε_{abcd} \bar{Ω}^{cd}, \hspace 1em \bar{J}_a = 0, \hspace 1em J_{ab} = 0.$$
The response field $G_{ab}$ picks up an extra curvature term $2m \sqrt{|g|} ε_{abcd} Ω^{cd}$, compared to the unadjusted response field. It arises from the boundary term $ε_{abcd} \sqrt{|g|} Ω^{ab} Ω^{cd}$ and drops out from the Euler-Lagrange equations (as a result of the "Bianchi identities" that came from $dF + AF - FA = 0$). The unadjusted source term $J_a$ picks up a contribution proportional to $λ G_{ab} θ^b$ that yields the 3-current density associated with the Einstein tensor.
