Laplacian of Lorenz gauge magnetic potential My textbook, Gettys's Physics (Italian language edition), says that the Lorenz gauge choice uses the magnetic vector potential $$\mathbf{A}(\mathbf{x},t):=\frac{\mu_0}{4\pi}\int \frac{\mathbf{J}(\mathbf{y},t-c^{-1}\|\mathbf{x}-\mathbf{y}\|)}{\|\mathbf{x}-\mathbf{y}\|}d^3y $$ which is said to be such that
$$\nabla^2\mathbf{A}(\mathbf{x},t)-\varepsilon_0\mu_0\frac{\partial^2 \mathbf{A}(\mathbf{x},t)}{\partial t^2}=-\mu_0\mathbf{J}(\mathbf{x},t)$$
How is this identity derived?
If $\mathbf{J}$ does not depend on time and is assumed, as often done in physics, to be of class $C_c^2(\mathbb{R}^3)$, I know, as proved here, that $$\nabla^2\mathbf{A}(\mathbf{x})=-\mu_0\mathbf{J}(\mathbf{x}),$$ consistently with the above equation, but I am not able to derive the general case. I heartily thank any answerer.
 A: So the best way to view this is that you're going to derive $\vec A(\vec x, t)$ from the equations that you want, first and foremost. So we have to start with what equations we want.
What equations do we want?
The Maxwell equations are:$$\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}$$And our idea is that we see $\nabla\cdot B = 0$ and we know that this means $B = \nabla\times A$ for some $A$. Now looking at $\nabla \times E = -\nabla \times \dot A $ we see that actually $\nabla \times (E + \dot A) = 0$ which means that there must also be a scalar potential $E + \dot A = -\nabla \phi$ for some $\phi$.
Now we ask, "to what extent was that choice of $A, \phi$ free, and to what extent were we constrained?" And the answer is that we have to preserve the fields $E$ and $B$. We know that we preserve $B$ if we add any $\nabla \psi$ to $A$ because the curl of a gradient is zero. But what does that do to our equation for $E$? It gives us $E + \dot A + \nabla \dot\psi = -\nabla \phi.$ Conclusion: we can add any $\nabla \psi$ to $A$ but only if we also subtract $\dot \psi$ to $\phi$, so that we preserve $E$ as well. 
This ability to add $\psi$ is called the gauge freedom and it is analogous to having a freedom to choose a constant of integration; $A$ and $\phi$ are in some sense integrations of $E$ and $B$. 
Now we have two more equations above that we haven't automatically guaranteed by the above construction. Using the identity $\nabla \times (\nabla \times X) = \nabla (\nabla \cdot X) - \nabla^2 X$ we can simplify these considerably.$$
-\nabla^2 \phi - \nabla \cdot \dot A = \rho/\epsilon_0, \\
\nabla (\nabla \cdot A) - \nabla^2 A =\mu_0 J - \mu_0\epsilon_0 \big(\nabla \dot \phi + \ddot A\big).
$$The first thing we see is the emergence of some operator $\Box X = \mu_0\epsilon_0 \ddot X - \nabla^2 X,$ and the next thing is some sort of potential field $\lambda = \nabla \cdot A + \mu_0 \epsilon_0 \dot \phi.$ Combining these into the above two equations gives:$$
\Box \phi = \rho/\epsilon_0 + \dot \lambda, \\
\Box A =\mu_0 J - \nabla \lambda.
$$
Okay now remember what I said about the gauge freedom? We can add $\nabla \psi$ to $A$ if we also subtract $\dot\psi$ from $\phi$? What does that do to $\lambda$ as defined above? It changes $\lambda\to\lambda - \Box\psi.$ It turns out that this exactly guarantees that the equations don't change (Exercise: prove this, then reflect on why that has to be true.) 
We can effectively use this to replace $\lambda$ with whatever we want. For example in the Coulomb gauge, we can first imagine that we solve for some arbitrary $\phi(\vec x,~ t)$ and some arbitrary $\vec A(\vec x,~ t),$ so $\lambda(\vec x,~t)$ is some complicated mess. But if we first solve $\Box \psi = \lambda - \mu_0 \epsilon_0 \dot \phi$ the resulting gauge transform maps $\lambda \to \mu_0 \epsilon_0 \dot\phi$ and the above equation becomes $-\nabla^2 \phi = \rho/\epsilon_0,$ solved with the standard Coulomb potential (hence this choice is known as the Coulomb gauge). You might object that this creates instantaneous action-at-a-distance, but remember that $E$ is not $-\nabla \phi$ in general; it is $-\nabla\phi - \dot A.$ The instantaneous action-at-a-distance in $\phi$ is balanced out by instantaneous action-at-a-distance in $A$ to keep the fields all right.
So the "equations that we want" explicitly contain this no-action-at-a-distance and come when we solve $\Box \psi = \lambda$ to gauge-transform $\lambda\to 0.$
In fact the two equations combine together into a "four-vector" expression: $\Box (\phi/c,~ A) = \mu_0 (c \rho, ~ J).$ And those are the equations that we want. (In college it took me a long time to grapple with the question "what happens if the equations of motion take the system to a different $\lambda$?" which turns out to be a pointless question. Recall that we're solving for the fields at all times.) 
Note that this last step of choosing the gauge basically says "you can always do this," but it doesn't quite construct it for you: our procedure right now is very good for analytical theory but very clumsy for practical theory: "find the fields, guess some $A$, figure out some $\phi$, calculate $\lambda$, solve for the $\psi$ of your dreams, use that to correct $A$ and $\phi$ whose equation no longer will have $\lambda:$ hooray, you're finally in a mathematically pretty place." We want to turn this process around: "take your $\rho, J,$ calculate the right $\phi, A$ in this mathematically pretty place, now use it to get the right $E, B$." Then we become a lean, mean theory machine!
Now how do we solve these equations?
Okay, we know that in 3D space we could solve expressions of the form $\nabla^2 \alpha = -\beta$ by the solution you provided, $$\alpha(\vec r)  = \frac{1}{4\pi} ~ \int_V d^3 r' ~\frac{\beta(\vec r')}{|\vec r - \vec r'|} = \int d^3 r' ~\eta(\vec r, \vec r').$$ There is a nice way to interpret this in terms of Dirac $\delta$-functions called the "Green's functions" approach; this says that our source equation involves a linear differential operator and therefore obeys the principle of superposition, so that if you know its response to a stimulus at any one point, $\beta = \delta_3(\vec r - \vec r')$ is the so-called "Green's function" $\alpha = f(\vec r, \vec r')$, then the actual driver is just a superposition $\beta(\vec r) = \int d^3r' ~\delta_3(\vec r - \vec r') \beta(\vec r'),$ and therefore by linearity the general solution is the weighted superposition of these "Green's functions," $\int d^3r~\beta(\vec r')~f(\vec r, \vec r').$
So the given function $-1/({4\pi|\vec r - \vec r'|})$ is the $\alpha$ you calculate when $\beta(\vec r) = \delta_3(\vec r - \vec r')$ is the 3D Dirac $\delta$-function. You can see this by recognizing the Laplacian as the divergence of a gradient; this describes a field whose divergence is 0 except at one point $\vec r'$ where the divergence suddenly spikes to infinity in such a way so that the surface integral bounding the point is 1. We know that that requires an inverse-square field which goes like $(\vec r - \vec r')/|\vec r - \vec r'|^3$ and we see that that sort of thing would come out as the gradient of $1/|\vec r - \vec r'|,$ and the $4 \pi$ comes from the fact that the surface area of the sphere is naturally $4 \pi r^2$.
So now we come to the idea of $\Box = c^{-2} \partial_t^2 - \nabla^2$ and we want to do the same thing, but we know that $\Box \alpha = 0$ is solved by a superposition of traveling waves, $\alpha(\vec r, t) = \sum_i f_i(\vec r - \vec u_i~t)$ where for all $i$ we know that $|\vec u_i| = c.$ If we imagine a $\delta$-function in space and time, we imagine that it would have to produce a spherical shell propagating in all directions, attenuating as it goes, and so we guess something like $\alpha(\vec r, t) = \delta(t - t' - |\vec r - \vec r'|/c) / (4 \pi |\vec r - \vec r'|).$
This guess turns out to be exactly correct, hence the general solution is the superposition over both coordinates: $$ \alpha(\vec r, t) = \int dt' ~\int d^3 r'~ \delta\left(t - t' - \frac{|\vec r - \vec r'|}c\right)~\frac{\beta(\vec r', t')}{4\pi|\vec r - \vec r'|}.$$Performing the integral over $t'$ you can just use the definition of the $\delta$-function to replace t' with $t - |\vec r - \vec r'|/c$, your expression above.
Alternatively, just apply $\Box$ to your expression. This is the "guess and check" method of solving an equation; I've given you the Green's function reason why you should guess this expression as the solution to the equation; now all you have to do is check it and if it's correct then it's correct!
For this you're simply doing a divergence of a gradient, where it's useful to know that $\nabla$ obeys the normal product rule for derivatives with $\nabla f\big(|\vec r - \vec r'|\big) = f'\big(|\vec r - \vec r'|\big) ~ (\vec r - \vec r')/|\vec r - \vec r'|.$ We have to say very carefully that $\dot \beta$ is the derivative of $\beta$ with respect to its second argument overall which we'll call $\tau$ for short. So you start with:
$$\nabla \left(\frac{\beta(\vec r', \tau)}{4\pi|\vec r - \vec r'|}\right) = \frac{\vec r - \vec r'}{|\vec r - \vec r'|} \left(\frac{\dot \beta(\vec r', \tau)}{4\pi|\vec r - \vec r'|} - \frac{\beta(\vec r', \tau)}{4\pi|\vec r - \vec r'|^2}\right),$$ and then do the same with the divergence. 
