given I perform a gauge transformation:
$(\vec{A},\phi)->(\vec{A'},\phi')$
Such that,
$\vec{A'} = \vec{A} + \nabla f$
$\phi' = \phi - \frac{\partial f}{\partial t}$
The new potentials, $\vec{A'},\phi'$, Leave the field invariant.
Performing the actual gauge transformation would change the potentials in the potential formulation of maxwell equation to be $\vec{A'}$ and $\phi'$
Suppose further, that aswell as leaving the field invariant, we would also like the new potentials $\vec{A'},\phi'$
To satisfy the condition that :
$\nabla \cdot \vec{A'} + \mu_0 \epsilon_0 \frac{\partial \phi'}{\partial t} = 0$
Thus, by substituting the definition of $\vec{A'},\phi'$ into the lorenz gauge condition. We are essentially saying this: Given the potentials $\vec{A'},\phi'$ satisfy the lorenz gauge condition, AND they are in a form that leaves the field invariant. What should the function "f" be?
If we can prove than an "f" exists given the previous 2 statements, then we have shown that:
Given the potentials $\vec{A'},\phi'$ satisfy the lorentz gauge condition, they can be written in a form that leaves the field invariant. Aka, we can prescribe the divergence of $\vec{A'}$ to follow the lorenz gauge condition. And the potentials still gives the correct E,B field.
Let's solve for the function F, to prove that F exists!
$\nabla \cdot \vec{A'} + \mu_0 \epsilon_0 \frac{\partial \phi'}{\partial t} = 0$
Sub in:
$(\vec{A'} = \vec{A} + \nabla f$),
$(\phi' = \phi - \frac{\partial f}{\partial t})$
$\nabla \cdot ( \vec{A} + \nabla f ) + \mu_0 \epsilon_0 \frac{\partial }{\partial t}(\phi - \frac{\partial f}{\partial t}) = 0$
$\nabla \cdot \vec{A} + \nabla^2 f + \mu_0 \epsilon_0 \frac{\partial \phi }{\partial t}-\mu_0\epsilon_0 \frac{\partial^2 f}{\partial t^2} = 0$
$\nabla^2 f -\mu_0\epsilon_0 \frac{\partial^2 f}{\partial t^2} = -(\nabla \cdot \vec{A} + \mu_0 \epsilon_0 \frac{\partial \phi}{\partial t})$
note the function on the right isn't the same as what you have gotten.
This IS a solvable equation( "inhomogenous wave equation"). meaning F exists. Meaning the potentials can be written in a form that leaves the field invariant whilst also satisfying the lorenz gauge condition.
To answer your question more directly, this is the conditions on F such that the new potentials satisfy the lorenz gauge
Let's set $c=1$. There really isn't anything special about retarded potentials - they are particular solutions of Maxwell's equations for the potential in the Lorenz gauge
$$ \partial_\mu\partial^\mu A^\nu = J^\nu\tag{1}$$
(where $A^0 = \phi$ and $J^0=\rho$ - since the Lorenz gauge conditions is Lorentz-invariant, we can use covariant language here). Since eq. (1) is the inhomogeneous version of
$$\partial_\mu\partial^\mu A^\nu = 0,\tag{2}$$
the general theory of inhomogeneous equations applies - if $A_\text{ret}^\mu$ is a particular of eq. (1) (for instance the usual retarded potentials) and $A_\text{hom}$ is any solution of eq. (2), then $A_\text{ret}^\mu + A_\text{hom}^\mu$ is a solution of eq. (1), too. The Lorenz gauge condition is linear as well: If $\partial_\mu A_\text{ret}^\mu = 0$ and $\partial_\mu A_\text{hom}^\mu = 0$, then $\partial_\mu(A_\text{ret}^\mu + A_\text{hom}^\mu) = 0$, too. In particular you can choose $A_\text{hom}^\mu$ to be the usual plane waves for EM waves in all of space that don't vanish at infinity.
This is exactly analogous to us solving the equations in terms of $\vec E$ and $\vec B$ and being able to add homogeneous solutions to one particular inhomogeneous solution - it's the same mathematical principle, and it's completely irrelevant that the equations are in terms of potentials and not fields now.
Best Answer
The gauge $\phi = A_0 = 0$ is called Weyl gauge or temporal gauge.
This gauge is incomplete, as one can see from the definition of a gauge transformation, $$A_\mu \to A_\mu + \partial_\mu \alpha(x).$$ We can still perform any gauge transformation with gauge parameter $\alpha$ independent of $t$, as this keeps $A_0$ the same. To remove some of the residual gauge freedom we could, e.g. impose $$A_z|_{t = 0} = 0.$$ The proof this gauge can be reached is just the exact same as the proof that Weyl gauge can be reached, except with effectively one less dimension since nothing depends on $t$. At this point we are still not done, because we can still preserve both gauge conditions using any $\alpha$ independent of both $t$ and $z$. So we impose the further condition $$A_y|_{t = z = 0} = 0$$ leaving only $\alpha$ dependent on $x$, which are removed by imposing $$A_x|_{t = z = y = 0} = 0.$$ These four conditions together are a complete gauge fixing. It's quite a mouthful, which is why you won't see it written out in textbooks too often.
Whether or not you want to perform a complete gauge fixing is up to taste. For example, in the standard presentation of the QCD $\theta$-vacua, one takes the incomplete gauge fixing $A_0 = 0$ and then argues there are multiple vacua $|n \rangle$. But there is a completely equivalent presentation where one takes the complete gauge fixing I gave above (also mentioned here), and finds a unique vacuum but the exact same physical effects. This is related to whether one chooses to regard large gauge transformations as "do-nothing" transformations.