Why can the Lorenz gauge condition always be fullfilled? Why is the Lorenz gauge condition always possible for classical electromagnetic fields?
So far I can only understand the following:
If we perform a gauge transformation $A\mapsto A'=A+\mathrm{d}\Lambda$, then the 'physical' field $F=\mathrm{d}A$ is unchanged. 
From our definition of new potential, we have $\mathrm{d}\star A'=\mathrm{d}\star A+\mathrm{d}\star \mathrm{d}\Lambda$. 
If the Lorenz condition is always true, i.e. we can always find a gauge such that $d\star A'=0$, then from the above equation, we conclude that for any 1-form potential $A$ satisfying $F=\mathrm{d}A$, we have a function $\Lambda:\mathbb{R}^{1,3}\rightarrow\mathbb{R}$ satisfying $\mathrm{d}\star \mathrm{d}\Lambda=-\mathrm{d}\star A$ at each point in $\mathbb{R}^{1,3}$. 
This seems to be a very strong statement. Why should I not expect any singularity for the function $\Lambda$ on $\mathbb{R}^{1,3}$?
 A: As Meng Cheng said, you need to solve a wave equation with sources (see, e.g., http://physics.gmu.edu/~joe/PHYS685/Topic6.pdf , eq.(3)). The wave equation's solution can be expressed as some integral (https://en.wikipedia.org/wiki/Wave_equation#Solution_of_a_general_initial-value_problem ), so you do need some integrability condition.
A: I can (in fact!) recommend Jackson's book "Classical Electrodynamics" for a discussion about the derivation. See the first sections in Chapter 6. 
Then, in section 6.3, he arrives at the equation that has to have a solution for the Lorenz Gauge to exist, see equation (6.18):
\begin{equation*}
\label{eq:1}
\tag{1}
\nabla^2 \Lambda - \frac{1}{c^2} \frac{\partial^2 \Lambda}{\partial t^2}
=
- \left( \vec{\nabla} \cdot \vec{A} + \frac{1}{c^2} \frac{\partial \Phi}{\partial t} \right)
\end{equation*}
Unfortunately, Jackson does not then go into a discussion about existence, so here is that addition in a rather informal way, and according to my understanding.
From the construction (up until this point, see Jackson) we know that the right-hand-side is some non zero field, and we wish to find a left-hand side that is a solution for any non-zero right-hand-side. We call the right-hand-side a source term, and for convenience we can rename it as $f(t,\vec{x})$.
\begin{equation*}
\label{eq:2}
\tag{2}
\nabla^2 \Lambda - \frac{1}{c^2} \frac{\partial^2 \Lambda}{\partial t^2}
=
f(t,\vec{x})
\end{equation*}
The left-hand-side of equation (\ref{eq:2}) is the famous wave-equation, and it is linear. This means that if we can find solutions for simple source terms, and then superpose those solutions (add them) such that the corresponding superposed source-term exactly match the complicated source $f(t,\vec{x})$, then the superposition of the solutions will be the actual solution to $\Lambda$ in equation (\ref{eq:2}).
That was a bit much, but anyway, this is the approach of using Green's functions. And to read up on them I can recommend chapter 10 in Mathematical Methods for Physicists, by Arfken et al. The bottom line is that you can make use the infinite superposition of an integral, once you have figured out the so-called Green's function, $G(t,\vec{x}\,;t',\vec{x}')$ for the wave-equation. And for the wave-equation it is possible to find $G$. Then $\Lambda$ is uniquely defined by carrying out an integral.
$$
\Lambda(t,\vec{x}) =
\int\!\!\!\!\int\!\!\!\!\int\!\!\!\!\int G(t,\vec{x}\,;t',\vec{x}') \; f(t',\vec{x}') \; \text{d} t' \; \text{d}^3 x'
$$
Thus there exist a $\Lambda$ for any $f(t,\vec{x})$. $\blacksquare$
