The fact is that, in the general case
$$
\vec{E} = -\vec{\nabla}V - \frac{\partial\vec{A}}{\partial t};
$$
(signs depend on conventions used) where $\vec{A}$ is called vector potential. You can consult for example Wikipedia.
Let us consider homogeneous Maxwell equations:
$$
\begin{cases}
\vec{\nabla}\cdot\vec{B} = 0,\\
\vec{\nabla}\times\vec{E} + \frac{\partial\vec{B}}{\partial t} = 0;
\end{cases}
$$
It is well-known that every divergenceless filed on $\mathbb{R}^3$ can be written a curl of another vector field just as we know that a curless field can be written as a gradient of a scalar function on $\mathbb{R}^3$. Thus from the first equation,
$$
\vec{B} = \vec{\nabla}\times\vec{A},
$$
and substituting this in the second equation,
$$
\vec\nabla\times\left(\vec{E} + \frac{\partial\vec{A}}{\partial t}\right)=0,
$$
since one can exchange the curl with the derivative w.r.t. time, and so one can set:
$$
\vec{E} + \frac{\partial\vec{A}}{\partial t} = -\vec\nabla V,
$$
from which
$$
\vec{E} = -\vec{\nabla}V - \frac{\partial\vec{A}}{\partial t}.
$$
Note that if your magnetic field is time-independent, you recover the well-know formula
$$
\vec{E} = -\vec\nabla V.
$$
So intuitively $\nabla \cdot B = 0$ is best understood by appeal to the continuity equation $\dot \rho = -\nabla \cdot J,$ in other words if you imagine $\vec B$ as a flow field $\vec v$ for some fluid of constant density, this describes a flow of stuff which does not accumulate at any particular points: it is not "flowing into" or "flowing out of" any of those places.
The expression $\nabla\times B = 0$ is a little harder to understand when you imagine $\vec B$ as a flow field $\vec v$, but if you imagine that you insert a little pinwheel into the flow field at a point, the curl says something about how much torque there will be on that pinwheel due to the fluid trying to drag its tines, with the direction of the curl being the direction that makes the pinwheel spin the most. So $\nabla\times B=0$ means literally "this flow field does not make a pinwheel spin." That helps to understand the strange flow field $\vec v = \alpha ~ \hat \theta / r,$ which seems to "curl" around the origin but has $\nabla \times \vec v = 0$ at all points that are not the origin. If you think about a tiny little area $dr~r~d\theta,$ it has a slightly longer "outside edge" than its "inside edge" and so the flow along that outside edge drags the pinwheel slightly more; but if the flow drops like $1/r$ this slight lengthening is exactly compensated by the slight attenuation in the flow field, meaning that the pinwheel is not rotated at all.
Helmholtz's theorem says that any "nice" vector field $X$ can be
"integrated" into a sum of two fields, a curl and a gradient: $$X = \nabla \times Y - \nabla \zeta.$$Since $\nabla \cdot X = -\nabla^2 \zeta$ and $\nabla \times X = \nabla \times (\nabla \times Y)$ due to standard vector identities, whenever we know that $\nabla \cdot X = 0$ we can solve this trivially with $\zeta = 0$ and we know that $X = \nabla \times Y$ for this "vector potential" $Y.$ Similarly when we know that $\nabla\times X = 0$ we know that $Y=0$ works and that $X = -\nabla \zeta$ for this "scalar potential" $\zeta.$
How this plays out in classical E&M
The full Maxwell Equations are (SI units) $$\begin{array}{ll}\nabla \cdot E = \rho/\epsilon_0&~~\nabla \times E = -\dot B\\
\nabla \cdot B = 0&~~\nabla \times B = \mu_0 J + \mu_0\epsilon_0 \dot E.
\end{array}$$
Clearly we always have $B = \nabla \times A$ for some vector potential $\vec A$ by the third equation. More subtly, the second equation then works out to $\nabla \times (E + \dot A) = 0$ and therefore $E = -\dot A - \nabla \phi$ for some scalar potential $\phi.$ The choices that you make for these are not unique as you can add any $\nabla \psi$ to $A$ and still preserve $B,$ because the curl of a gradient is zero. If you work out the effect of this on $E$ you'll see that you have to also subtract $\dot \psi$ from $\phi$ to preserve $E$ as well. This is called the "gauge freedom" but it really reflects this freedom to choose your "constants of integration" in the Helmholtz decomposition.
The remaining equations work out to $$\begin{array}{rl}-\nabla\cdot \dot A -\nabla^2 \phi ~=& \rho/\epsilon_0\\
\nabla(\nabla \cdot A) - \nabla^2 A ~=& \mu_0 J - \mu_0\epsilon_0\ddot A - \mu_0\epsilon_0\nabla \dot \phi,\end{array}
$$ once you use this identity that $\nabla \times (\nabla \times X) = \nabla(\nabla\cdot X) - \nabla^2 X.$ These can be put into a very elegant form by recognizing the wave operator $\square X = \mu_0\epsilon_0 \ddot X - \nabla^2 X$ and then incorporating the other spare terms in the second equation into a scalar field $\lambda = \mu_0\epsilon_0\dot \phi + \nabla\cdot A.$
Note that your gauge freedom maps $\lambda \mapsto \lambda - \square \psi$ so you can basically choose any functional form you want for $\lambda$ by solving a wave equation, in theory. Note also that $\lambda$ enters into the first equation via replacing $\nabla\cdot \dot A$ with $\dot \lambda - \mu_0 \epsilon_0 \ddot \phi,$ so you get $$\begin{array}{rl}\square \phi ~=& \rho/\epsilon_0 + \dot \lambda\\
\square A ~=& \mu_0 J - \nabla\lambda.\end{array}$$The "Lorentz gauge" sets $\lambda = 0$ because $(c\rho, J)$ is a 4-vector in relativity where $c = 1/\sqrt{\mu_0\epsilon_0}$ and hence this gives you an equation $\Box(\phi/c, \vec A) = \mu_0 (c\rho, \vec J),$ and since $\Box$ is relativistically covariant you get that $(\phi/c, \vec A)$ is also a 4-vector in relativity. If that doesn't make much sense to you just accept "hey, it makes these equations look exactly alike" and that's good enough.
The "Coulomb gauge" sets $\lambda = \mu_0\epsilon_0 \dot\phi$ so that the first equation reduces to just $-\nabla^2\phi = \rho/\epsilon_0,$ and we can basically assume that the scalar potential $\phi$ updates instantaneously to all of the charges, which makes it super-simple to calculate. If you look up at the definition of $\lambda$ you find out that it also sets $\nabla\cdot A = 0$ so $A$ can be integrated in terms of a vector potential potential, $A = \nabla \times W.$
Best Answer
Charge conservation demands that $\text{div} \textbf{J} = -\frac{\partial \rho}{\partial t}$ hold always. Thus in regions where there is no electric current, $\textbf{J}=0$, the charge density $\rho(\textbf{r},t)$ must be constant in time, $\rho= \rho(\textbf{r})$ and thus contributes nothing to radiation (which is what Balanis is interested in) and can be neglected. There follows $\text{div}\textbf{D}=0$.
The purpose of introducing the "electric vector potential" is to exploit the mathematical analogy with the physically more meaningful $\textbf{A}$ vector. This is similar to what is done in magnetostatics where the concept of magnetic potential is used in analogy with the conventional electric potential used in electrostatics. Both potentials satisfy Laplace's equations $\text{divgrad}\psi=0$ in regions free of sources, and as long as one stays in simple, singly connected domains the analogy holds.
The use of the $\textbf{F}$ potential follows similar line of thinking. One solves with much pain the electric dipole radiation problem and then to solve for the radiation from a loop one replaces the loop with an equivalent magnetic dipole, a kind of dual to the electric one.