Actually, it can be theoretically derived from D'Alembert equation (that is satisfied by each component of ${\bf E}$ and ${\bf B}$ in absence of sources in view of free Maxwell equations). The idea is to compute the field (any component of ${\bf E}$ or ${\bf B}$) in $p$, when it is generated by a spherical point source localized in $q$ emitting a spherical monochromatic field with fixed (scalar) wavenumber $k$, and between $q$ and $p$ there is a screen with an aperture of known area. The mathematical tool is an integral formula, due to Kirchoff, which produces the solution in $p$ when it is known the value of the field and its normal derivative on a surface surrounding $p$. The surface is chosen to have a part adapted to the screen, including the aperture, and the remaining part is taken far away from $p$. Here, i.e., to fix the value of the field and its normal derivative, some approximations enters the computation and they usually have physical sense for $k >> d$, where $d$ is the "diameter" of the aperture. This situation is discussed in details in Jackson's textbook. The final formula obtained this way can be shown to be equivalent to directly apply Huygen's principle from scratch.

Actually, it can be theoretically derived from D'Alembert equation (that is satisfied by each component of ${\bf E}$ and ${\bf B}$ in absence of sources in view of free Maxwell equations). The idea is to compute the field (any component of ${\bf E}$ or ${\bf B}$) in $p$, when it is generated by a spherical point source localized in $q$ emitting a spherical monochromatic field with fixed (scalar) wavenumber $k$, and between $q$ and $p$ there is a screen with an aperture of known area. The mathematical tool is an integral formula, due to Kirchoff, which produces the solution in $p$ when it is known the value of the field and its normal derivative on a surface surrounding $p$. The surface is chosen to have a part adapted to the screen, including the aperture, and the remaining part is taken far away from $p$. Here, i.e., to fix the value of the field and its normal derivative on the surface, some approximations enter the computation and they usually have physical sense for $k >> d$, where $d$ is the "diameter" of the aperture. This situation is discussed in details in Jackson's textbook. The final formula obtained this way can be shown to be equivalent to directly apply Huygen's principle from scratch.