As I said, Heras's paper is utterly confusing. He's basically going in circles, and also his notation is unclear.
The logic in classical electrodynamics is as follows:
The fundamental laws are Maxwell's equations, which are expressed in completely in gauge-invariant terms. In (1+3)-dimensional notation in Heaviside-Lorentz units they read
$$\vec{\nabla} \times \vec{E} + \frac{1}{c} \partial_t \vec{B}=0, \quad \vec{\nabla} \cdot \vec{B}=0.$$
These homogeneous Maxwell equations are basically constraint equations for the various field components. Then there are the inhomogeneous Maxwell equations, which relate the field components to the sources:
$$\vec{\nabla} \times \vec{B} - \frac{1}{c} \partial_t \vec{E}=\frac{1}{c} \vec{j}, \quad \vec{\nabla} \cdot \vec{E}=\rho.$$
The homogeneous equations imply that the 6 field components can be expressed in terms of a scalar and a vector potential,
$$\vec{E}=-\vec{\nabla} \Phi -\frac{1}{c} \partial_t \vec{A}, \quad \vec{B}=\vec{\nabla} \times \vec{A}.$$
The potentials are defined only up to a gauge transformation, i.e., if ##\Phi## and ##\vec{A}## describe a physical situation then with any scalar field ##\chi##, also the potentials
$$\Phi'=\Phi+\frac{1}{c} \partial_t \chi, \quad \vec{A}'=\vec{A}-\vec{\nabla} \chi$$
describe the very same physical situation, indeed evaluating the fields from the new potentials you get
$$-\vec{\nabla} \Phi'-\frac{1}{c} \partial_t \vec{A}' = -\vec{\nabla} (\Phi+\frac{1}{c} \partial_t \chi) -\frac{1}{c} \partial_t (\vec{A}-\vec{\nabla} \chi)=-\vec{\nabla} \Phi-\frac{1}{c} \partial_t \vec{A} =\vec{E}$$
and
$$\vec{\nabla} \times \vec{A}'=\vec{\nabla} \times (\vec{A}-\frac{1}{c} \vec{\nabla} \chi) = \vec{\nabla} \times \vec{A}=\vec{B}.$$
Now we can use this "gauge freedom" to impose one "gauge condition" to (at least) partially fix the potentials. This can be done to make life as easy as possible.
To see what gauge choice makes life as easy as possible, we express the fields in the inhomogeneous Maxwell equations, which remain to be solved, in terms of the potentials:
$$\vec{\nabla} \times (\vec{\nabla} \times \vec{A})+\frac{1}{c} \partial_t (\vec{\nabla} \Phi+\frac{1}{c} \partial_t \vec{A}) = \frac{1}{c}\vec{j}.$$
This you can rearrange a bit further to get
$$(\frac{1}{c^2} \partial_t^2 -\Delta) \vec{A} + \vec{\nabla} (\vec{\nabla} \cdot \vec{A} + \frac{1}{c} \partial_t \Phi)=\Box \vec{A} + \vec{\nabla} (\vec{\nabla} \cdot \vec{A} + \frac{1}{c} \partial_t \Phi)=\frac{1}{c} \vec{j}.$$
The other inhomogeneous equation becomes
$$\vec{\nabla} \cdot (\vec{\nabla} \Phi+\frac{1}{c} \partial_t \vec{A})=-\rho$$
or
$$\frac{1}{c} \partial_t (\vec{\nabla} \cdot \vec{A}) +\Delta \Phi=-\rho.$$
Now there are two choices where these equations become "convenient". In some sense the most easy choice is the Lorenz gauge, defined by the contraint
$$\frac{1}{c} \partial_t \Phi + \vec{\nabla} \cdot \vec{A}=0.$$
Then in the equation with ##\vec{j}## on the right-hand side the second term vanishes, and you get simply three separate wave equations for the components,
$$\Box \vec{A}=\frac{1}{c} \vec{j}.$$
In the equation with ##\rho## on the right-hand side you use the Lorenz gauge condition to get rid of ##\vec{A}## in this equation, also yielding the simple wave equation for the scalar potential,
$$\Box \Phi=\rho.$$
Not he box operator has no unique "inverse", but there's a most convenient choice, the "retarded solutions", which give the potentials at times ##t## in terms of the sources at time arguments ##t'<t## only, i.e., you evaluate the potentials at time ##t## only from the past history of the sources' evolution. This matches the causality principle of physics: All quantities at present time depend only on the past (or the present) state of the system but never on its future. This is manifest in Lorenz gauge for the potentials, i.e.,
$$\vec{A}(t,\vec{x})=\int_{\mathbb{R}^3} \mathrm{d}^3 x' \frac{j(t-|\vec{x}-\vec{x}'|/c,\vec{x}')}{4 \pi c|\vec{x}-\vec{x}'|}, \quad \Phi(t,\vec{x}) = \int_{\mathbb{R}^3} \mathrm{d}^3 x' \frac{\rho(t-|\vec{x}-\vec{x}'|/c,\vec{x}')}{4 \pi |\vec{x}-\vec{x}'|}.$$
This makes even the somewhat stricter relativistic sense of the causality principle manifest: Nothing physical can spread faster than the speed of light. Indeed here the potentials depend on sources at times earlier by the time light needs to travel from the source point ##\vec{x}'## and the point of observation ##\vec{x}##. Since the fields are given by derivatives ##\partial_t## and ##\vec{\nabla}##, this holds true for them too, and only for them it always must hold true, because they are the only really observable fields.
This becomes clear if you choose the Coulomb gauge condition for the potentials, which often has also its advantages. There the idea is to get rid of ##\vec{A}'## (we label the potentials in the Coulomb gauge with a prime to dinstinguish them from those in the Lorenz gauge above) in the equation with ##\rho## on the right-hand side in the most simple way by demanding the Coulomb-gauge condition,
$$\vec{\nabla} \cdot \vec{A}'=0.$$
Then the said equation simply becomes
$$-\Delta \Phi'=\rho.$$
This is solved as in electrostatics with the "instantaneous integral"
$$\Phi'(t,\vec{x}) = \int_{\mathbb{R}^3} \mathrm{d}^3 x' \frac{\rho(t,\vec{x}')}{4 \pi |\vec{x}-\vec{x}'|}.$$
Now the equation with ##\vec{j}## on the right-hand side, gets somewhat more complicated, but can be written as a wave equation with the source containing not only ##\vec{j}## but also the just calculated scalar potential:
$$\Box \vec{A}'=\frac{1}{c} (\vec{j} + \partial_t \vec{\nabla} \Phi')=\frac{1}{c} \vec{j}_{\perp}.$$
This you can solve again with help of the retareded integral,
$$\vec{A}'(t,\vec{x})=\int_{\mathbb{R}^3} \mathrm{d}^3 x' \frac{\vec{j}_{\perp}(t-|\vec{x}-\vec{x}'|/c,\vec{x}')}{4 \pi c |\vec{x}-\vec{x}'|}.$$
You can easily check that the Coulomb-gauge potentials, ##\Phi'## and ##\vec{A}'## are related to the Lorenz-gauge potential ##\Phi## and ##\vec{A}## by a gauge transformation. It's a good exercise to determine the corresponding gauge field ##\chi##.
The upshot is that with the Coulomb gauge, despite its mixture of "instantaneous" and "retarded" integrals leads to precisely the same electromagnetic field which is entirely given by "retarded" integrals, as it should be since the Maxwell equations fulfill the stricter notion of causality according to relativity, because the Maxwell equations are in fact a relativistic field theory of a "massless vector field", although that's hidden in using the (1+3)-notation. Everything becomes much more lucid when writing it in the more appropriate Minkowski-space 4-vector notation, which is the natural mathematical language of everything related to relativity.
