We know these:
$$
\text{In free space:}
\left\{
\begin{array}{l}
\nabla\times\mathbf{E}=0 \\
\nabla\cdot\mathbf{E}=\rho/\epsilon_0 \\
\end{array}
\right.
\ \ \ \ \ \ \text{In matter:}
\left\{
\begin{array}{l}
\nabla\times\mathbf{D}=\nabla\times\mathbf{P} \\
\nabla\cdot\mathbf{D}=\rho_{f}/\epsilon_0 \\
\end{array}
\right.
$$
Where the electric displacement is defined by $\mathbf{D}:=\epsilon_0\mathbf{E}+\mathbf{P}$.
$$\nabla\times\mathbf{E}=0 \implies \mathbf{E}=-\nabla V$$
By definition,
$$dV=\nabla V\cdot d\mathbf{l}$$
So if $d\mathbf{l}$ is perpendicular to $\nabla V$, V will not change. This means:
In 2D, If you start from any point on the plane, moving in the direction perpendicular to field lines, you have to form a closed loop. Because the surface you are in, must intersect with the equi-potential surface associated with the potential at your position. And there must exist one closed loop on that plane on which all points share the same potential with that of your position. That's because equi-potential surfaces are closed ones, no matter at infinity or not.
In 3D, If you start from any point in space, sweeping any possible point on your perpendicular-to-field-lines path, you have to from a closed surface.
This doesn't hold for $\mathbf{D}$ unless you can prove $\mathbf{P}$ is irrotational, that is, to prove it has these properties. This is often done by symmetry.
By definition,
$$\nabla\cdot\mathbf{E}=\lim_{\Delta V\to 0}\left[\frac{\oint_\mathcal{S} \mathbf{E}\cdot d\mathbf{a}}{\Delta V}\right]=\rho/\epsilon_0$$
This means if you consider a very small closed surface, The lines which come in, are exactly the same amount as the lines which get out of, the surface; if there is no charge enclosed by the small surface. This means those lines must be continuous in charge-free space.
In fact, Divergence is a measure of sources. So if you have discontinuities in some points, in which, the outgoing lines are more than incoming lines, you must be on a positive source point.
This holds for $\mathbf{D}$ too, as long as you neglect bound charges, those which are not much real ones, but almost our creation.
By application of Stokes's law on $\nabla\times\mathbf{E}=0$, we also get:
$$\oint_\mathcal{P}\mathbf{E}\cdot d\mathbf{l}=0$$
But this is not any help to our sense one field lines. Instead, we can play with the gradient formula to achieve another not-much-different-from-the-first-rule rule. :D
We know, by $\mathbf{E}=-\nabla V$, field lines must point to the steepest local decrease in the amount of potential. To if you go where $\mathbf{E}$ tells you, your potential should decrease the fastest. So if it turned out that your potential hasn't changed, or has even increased, something must be wrong about the field lines.
Some tracing will help you make sure what is drawn in front of you exists there in the real world.
But the first and the third point assume a previous sense of potential.
Looking at this integral:
$$V(\mathbf{r})=\frac1{4\pi\epsilon_0}\int\rho(\mathbf{r'})\frac{\mathbf{r}-\mathbf{r'}}{\left|\mathbf{r}-\mathbf{r'}\right|^3}d\tau'$$
It seems that potential is like a measure of positiveness inheritance which decays with getting more and more distant sharply.
If you are near a positive charge, you will inherit more positiveness. If you get farther from it, your positiveness will decrease kinda fast.
I think this is the most I could provide. And that this is probably the end of it and you can't dig more to get more sense. Instead of this, one could stick to mathematics to discover new fact!
Good luck.