The Huygens-Fresnel principle comes from looking at the Green's function for the wave equation.
If we look at the response of free space to a point disturbance, it turns out to be a spherical propagating wave, hence Huygens principle.
Mathematically, we solve for the Green's function $G$ in the following equation.
$$\left(\nabla^2-\frac{1}{c^2}\frac{\partial^2}{\partial t^2}\right)G(\mathbf{x},t;\mathbf{x}_0,t_0) = \delta(\mathbf{x}-\mathbf{x}_0) \delta(t-t_0)$$
With some math, this can be shown to be a spherically propagating wave.
The idea behind Green's functions is that once you know $G$, you can write any wave in roughly the following way:
$$u(\mathbf{x},t) \sim \int d\mathbf{x}^{\prime} G(\mathbf{x},t;\mathbf{x}^{\prime},0)u(\mathbf{x}^{\prime},0)$$
This is why, given the initial distribution of the wave $u(\mathbf{x}^{\prime},0)$, you can find the wave at a later time by propagating each point of the wave out with the spherical wave Green's function. For the exact form for the wave, see here https://math.dartmouth.edu/~ahb/notes/waveequation.pdf
Edit: As pointed out in the comments, there is a bit of subtlety going on with getting the homogenous solution from the inhomogenous Green's function, but looking at OP's question, that is more of a technicality. Using the inhomogenous Green's function for the inhomogenous case and vice-versa is justified by Duhamel's principle for Partial Differential Equations