Huygens' Principle geometric derivation and elimination of the wake and backward wave

Huygens' Principle (1678) implies that every point on a wave front serves as a source of secondary wavelets, and the new wave front is the tangential surface to all the secondary wavelets. But two problems arise: portions of wavelets that exist outside of the new wave front combine to form a wake. Also there are two tangential surfaces so wave fronts are propagated in both the forward and backward directions. These problems have not previously been resolved by using a geometrical theory with impulsive wavelets that are in harmony with Huygens' geometrical description. Doing so would provide deeper understanding of and greater intuition into wave propagation, in addition to providing a new model for wave propagation analysis. The interpretation, developed here, of Huygens' geometrical construction shows Huygens' Principle to be correct: as for the wake, the Huygens' wavelets disappear when combined except where they contact their common tangent surfaces, the new propagating wave fronts. As for the backward wave, a source propagates both a forward wave and a backward wave when it is stationary, but it propagates only the forward wave front when it is advancing with a speed equal to the propagation speed of the wave fronts.

Huygens' Principle implies that every point on a wave front serves as a source of secondary wavelets. The new wave front is the tangential surface to all the secondary wavelets in the direction of propagation ( Fig. 1 presents this as a geometric construction).
An abundance of literature refers to Huygens' Principle 1-3 . References 4-32 are a very small sample which show the tremendous variation in the application of Huygens' Principle and some continuing work on it. But there has been insufficient theory developed which addresses the underlying Huygens' basic geometrical construction of how a wave propagates via Huygens' wavelets. "At present mathematical physicists are busy with the rapid developments of modern physics, and older problems are likely to be neglected" (referring to Huygens Principle in 1940) 33 . This is despite the common use of Huygens' construction to introduce wave theory, reflection, refraction, and diffraction. In his construction two problems arise: the portions of the wavelets not contacting the position of the new propagating wave must disappear ('no wake') 2,3 . Also, the portions of the wavelets propagating in the direction opposite to the new propagating wave front must vanish ('no backward wave'). There has been no clear direct geometrical theory with impulsive wavelets that are in harmony with Huygens' geometrical description which shows how this can occur.
Huygens' Principle's problems of the backward wave and the wake have been approached before by relying on physics (for example, suitable initial values of the velocity potential and condensation 3 ), or by involving doublets, dipoles, obliquity factors, sinusoidal waves, and also by using complex mathematics 2, 3 , "..but the intuitive appeal of Huygens' simple principle is lost... " 4 along with a clear direct connection to his geometrical construction and the physical insight that derives from it.
Wave propagation is linear so superposition holds: it should be possible to decompose an impulsive propagating wave front into its constituent points, then consider the impulsive wavelets radiating from each of those points at a future time, and combine those wavelets in a simple direct geometric manner to obtain the progressing wave front at that future time.
The following shows that Huygens' geometrical construction is correct as depicted in his figures when the wavelets are interpreted as Dirac Delta distributions (Delta functions) [34][35][36][37][38] , and when the propagating wave field resulting from the summation of the wavelets is differentiated to yield the wave front, and also when the speed of the source is taken into account. Given this, the backward wave front and wake both vanish.
The following development is geometric and does not involve the wave equation or any particular type of field (except in the Supplementary Information Notes). The symbol φ is used for the propagating wave fields Stationary planar source and elimination of the wake. (The infinite planar source may sometimes be referred to as 'planar source' or 'plane'). This section analyzes the waves propagating away from a motionless infinite planar source. Each point on the plane is a source of a Huygens' wavelet consisting of a spherical Dirac delta function (distribution) expanding as a function of time at a rate equal to c, the propagation speed. The strength (or amplitude) of the wavelet is α which is attenuated as the wavelet expands by spherical spreading.
It is shown that the wave field (summation of the Huygens' wavelets) observed as a function of time at some fixed point is a constant value after the wavelets first reach that point. The partial derivative with respect to time of that wave field gives the wave front which replicates the original wave (the Huygens' wavelet) at 1/2 amplitude in both the advancing and retreating directions. The derivative also shows there is no wake. See Fig. 2.
It is only necessary to derive the wave field for each point on any one line perpendicular to the source since the planar source is uniform with respect to the x and y coordinates. Only the points in the planar source intersected by the sphere with radius ct centered at P can contribute to the wave field φ(z 0 , t) when deriving the wave field at P at time t (Fig. 2a).
Initially the duration ε of the excitation pulse δ f will be finite and the pulse will have an height equal to α. Then the area (or strength) under this finite version of δ(t) is εcα which is required to be 1. Later, where needed, the limit as ε tends toward zero will be taken while still requiring α to take on a value such that εcα remains equal to 1 (i.e. α = 1/[εc]) so that δ f → δ as ε → 0.
The locus of all radiating points on the planar source which can contribute to the wave field φ(z 0 , t) during the time interval [t, t + ε] and at distance z 0 from the plane is an annular area. This annular area is bounded by two concentric spherical shells of radius ct and c[t + ε].
The area Ap of the annulus is: where and So which has the limit The wave front f in the Z + direction that is derived from the wave field which will arrive at P at t = z 0 /c. www.nature.com/scientificreports/ (and so δ f → δ). Spherical spreading attenuation specified to be of the form 1/[4πct] attenuates the wavelets and therefore the radiation from the annular area. As a result the attenuated radiation from Ap is not a function of the field point's distance, z 0 , from the planar source (except for a step at ct = z 0 as the wave front passes). This will be shown in the following: Combining the 1/[4πct] spherical spreading attenuation and the α height factor with the area Ap gives the (dimensionless) amplitude of the wave field, φ, after the impulsive spherical wavelets have all passed (i.e. for |z|≤ ct), yielding: Because cεα = 1, the result simplifies to or using a Heaviside step function, θ, which as a spatial function of z is a rectangular pulse expanding in time in both advancing (toward P) and retreating (away from P) directions (Fig. 2b).
In the advancing direction the wave field observed at P as a function of time is which has a transition at t = z 0 /c when the wave field has advanced to P. The derivative of φ yields an advancing planar wave front, f, which observed at P as a function of time is (Fig. 2c): Or for z in general This holds for every line perpendicular to the source. Note that only the single point source at the nearest point N contributes to the wave front at P: the wave field φ is a step function and its derivative is zero except at the step's location, and the step was caused by the wave field originating at N (Fig. 2a). (So there is a one to one correspondence between the points in the wave front and the points in the planar source.) Because the derivative is zero except at the step, there is no wake.
There will also be a retreating plane wave front f -(z, t) = ½ δ(t + z/c) propagating in the opposite direction, away from P, which can be derived similarly. This also propagates without a wake. (Here z is negative, and f will not be seen at z 0 , the location of P).
Therefore an infinite planar distribution of point sources radiating impulsive spherical wavelets creates advancing and retreating impulsive plane wave fronts that propagate cleanly without a wake. This shows that the Huygens' wavelets cancel each other when summed together and differentiated except on their common tangential surfaces (Fig. 1a,b).
It might be expected that the previous process can be repeated and the impulsive planar wave fronts made to progress to the next future time t + Δt, and so on. But because of the ½ factor the waves would be attenuated upon each iteration. When the source is moving at the same speed as the advancing wave front the ½ factor disappears along with the retreating or backward wave front, as will be shown: Moving planar source and elimination of the backward wave. This section analyzes the changes that occur in the results of section "Stationary planar source and elimination of the wake" when the planar source is moving with speed v toward the observation point. It is shown that the results from section "Stationary planar source and elimination of the wake" can be used if the effects of motion on the Huygens' wavelets are included. When this is done it is seen that a new additional wave is created as a function of the planar source motion, v. This additional wave reduces the amplitude of the backward wave and increases the amplitude of the forward wave. When v = c the backward wave is eliminated and the forward wave amplitude is doubled eliminating the 1/2 factor. In any case the original wave is replicated. See Fig. 3.
During the time interval ε the planar source moves a distance vε in the direction toward P. Then R2 is the radius of the base of a right circular cone of slant height c[t + ε] and R1 is the radius of the base of a right circular cone of slant height ct. These two bases define a conical frustum of height vε and slant height d. The lateral area of this frustum, π[R2 + R1]d, is the locus of all radiating points which can contribute to the field at the field point P during the time interval [t, t + ε]. This is analogous to the annular area on the infinite stationary plane previously considered. www.nature.com/scientificreports/ However, in the limit of ε → 0 the area of the frustum reduces to the area of that annulus as follows: The annulus area is π[R2 2 − R1 2 ] and the frustum area is π[R2 + R1]d. Then referring to Fig. 3, Therefore as ε → 0 the frustum area becomes π[R2 , which is also the area of the annulus.
Consequently the previously developed results and formulas for the infinite stationary plane can be used. However, the effects of the motion of the moving source on its point sources' initial radiated impulses must still be considered: See Fig. 4a-c. If the planar impulsive source is moving with positive source speed, v, toward P, the motion results in a displacement during the pulse time interval ε (an additional initial condition) which causes additional wave front pulses to propagate in both the forward and backward directions, as with the stationary planar source. (Displacement at x is approximately the product of the temporal integration interval and the displacement speed evaluated at x. See Supplementary Notes B and E) This is in addition to the previously derived wave front pulses from the stationary planar impulsive source. The displacement height is the same as the original forward wave front pulse, α/2, whether longitudinal or transverse. The displacement spatial extent is vε.
The source motion in the forward direction displaces a positive vεα/2 during the time interval ε which propagates forward as an additional wave front pulse. Since α = 1/[cε] the displacement or strength equals [1/2] [vε]/[cε] = v/2c. Do to linearity this additional wave front pulse is added to the original stationary planar source wave front pulse. Note that as ε → 0 a pulse of width ε and height 1/ε moving at speed c tends to δ(x − ct). Note also that v and c are scalars and do not have to be measured along the same direction-they effect pulse strength by different means.
Then the total resulting wave front in the forward direction is the sum of two wave front pulses which, as ε → 0, can be combined to yield: See Fig. 4d-f. Proceeding similarly with the equalizing additional backward wave where the displacement is − vεα since the propagation direction is opposite to the source motion direction. The source motion displaces an equalizing negative − vεα/2 during the pulse time interval ε which then propagates backward as an additional wave front pulse. Since α = 1/[cε] the displacement or strength equals [1/2] As before this wave front pulse is added to the original stationary planar source wave front pulse. Then the total resulting wave front in the backward direction (where z is negative) is the sum of two wave front pulses which, as ε → 0, can be combined to yield: The sign is negative for the additional wave front pulse that is propagating in the direction opposite to the motion of the source (in the backward direction). This relative direction of motion defines the terms 'forward' and 'backward' in the case of the moving planar impulsive source. (Note that the terms are not defined for a stationary planar impulsive source).
So the impulsive forward and backward planar wave fronts are Figure 3. The geometry of an infinite planar impulsive source moving toward P with speed v. (a speed away from P would be given by − v.) R1 is the radius of the circle created by the intersection at time t of a sphere of radius ct with the plane. R2 is the radius of the circle created by the intersection at time t + ε of a sphere of radius c[t + ε] with the plane. www.nature.com/scientificreports/ (At P at a positive distance along the z axis, only the forward wave front will be observed.) For the stationary planar source (where v = 0), we get f + (z,t) = ½ δ(t − z/c) for the forward wave front and f − (z,t) = ½ δ(t + z/c) for the backward wave front, as was derived before. See Fig. 4g-i. For source speed v = c the forward wave front is f + (z, t) = δ(t − z/c) or as seen at P as a function of time, f + (z 0 , t) = δ(t − z 0 /c). So the forward wave front propagates without changing and becomes the new progressing Huygens' wave front (Note the absence of the ½ factor). For the backward wave front f -(z,t) = 0. So there is no backward wave front (See SI Notes C, D).
Non expanding spherical source and elimination of the wake. This section analyzes the waves propagating away from a non expanding spherical source. Each point on the sphere is a source of a Huygens' wavelet consisting of a spherical Dirac delta function (distribution) expanding as a function of time at a rate equal to c, the propagation speed. The strength (or amplitude) of the wavelet is α which is attenuated as the wavelet expands by spherical spreading.
It is shown that the wave field (summation of the wavelets) observed as a function of time at some fixed point is a constant value after the wavelets first reach that point and remain constant for a period of time that is proportional to the sphere's diameter. The partial derivative with respect to time of that wave field gives the wave front which replicates the original wave (the Huygens' wavelet) at 1/2 amplitude in both the expanding and converging directions. That derivative also shows there is no wake. See Fig. 5.
Since the spherical source is symmetrical, it is only necessary to derive the wave field for each point on any one radial line emanating from its center. The wave field will be the same on any other radial line. As before, initially it is assumed that the temporal duration ε of the excitation pulse δ f is finite and that it has a height equal to α and an area equal to εcα, where α = 1/[εc] to maintain an area equal to 1 as ε → 0 (and δ f → δ).The surface area of a spherical cap is 2πR 0 h(t) (Fig. 5a), and the surface area of the spherical zone is As = 2πR 0 [h(t + ε) − h(t)] (Fig. 5b); Note surprisingly that the surface area of the zone depends only on its height, h(t + ε) − h(t), not on its vertical position within the sphere.
To find As the heights h(t) and h(t + ε) must be found. Two right triangles can be formed (Fig. 5a):   (The spherical shell may be referred to as a 'sphere' hereafter.) The radius is R 0 , and the sphere is centered at the coordinate origin and has impulsive excitation δ f (t). Each point on the sphere's surface represents a point source which radiates an impulsive spherical wavelet. The parameter z 0 is the radial distance from the sphere's surface to the point in space, P, at which the resultant wave field φ will be observed as a function of time. The propagation speed is c. The distance ct from P to a point on the sphere defines a spherical cap of height h(t) with radius R1. The sphere is transparent to radiation so that the whole sphere will contribute to the wave field at any one point in space. (b) The locus of all contributing radiating points is a spherical zone. These are the only points on the spherical source of radius R 0 which can contribute to the field at P during the time interval [t, t + ε] and at distance z 0 from the sphere. The zone on the spherical source is bounded by its intersection with two concentric spheres of radius ct and radius c[t − ε] centered at the field point P. The intersections defines spherical caps of heights h(t) and h(t + ε), also radii R1 and R2. www.nature.com/scientificreports/ and as ε → 0 (so δ f → δ) which is the same as the area of the stationary infinite plane annulus Ap but attenuated by the spherical spreading factor R 0 /[R 0 + z 0 ] (the radial lines represent infinitesimal solid angles). This area is also the same as the area of the plane's annulus Ap when R 0 → ∞ or when z 0 → 0. As before spherical spreading of the form 1/[4πct] attenuates the wavelets and therefore the radiation from the spherical zone. The dimensionless amplitude of φ observed at P for z 0 /c ≤ t ≤ [2R 0 + z 0 ]/c is produced by combining the area As with the attenuation 1/[4πct] and the pulse height α: Because cεα = 1 which holds for z 0 ≤ ct ≤ [2R 0 + z 0 ]. (See Supplementary Note A for a different derivation) For fixed z 0 , |φ| observed at P is a rectangular pulse in time having temporal width 2R 0 /c. Note that as z 0 → 0 or R 0 → ∞, |φ|→ 1/2, which is the same as for the planar source.
The wave field φ observed at P as a function of time is produced by using Heaviside step functions to implement the time dependence, z 0 ≤ ct ≤ [2R 0 + z 0 ]: The derivative of φ yields the wave front observed at P as a function of time (Fig. 6b): For finite z 0 and t in the equation for f(z 0 ,t), as R 0 → ∞ the argument of δ(ct − z 0 − 2R 0 ) never equals zero, so never 'activates' and δ(ct − z 0 − 2R 0 ) may be ignored. Then which is the same as the equation for the planar source advancing wave front.
Both propagating spherical wave fronts (Fig. 6a) are affected by the spherical spreading term R 0 /[R 0 + z 0 ] which supplies the additional attenuation due to the radius increasing from the initial R 0 to R 0 + z 0 (it is implied that spherical spreading attenuation is included in the initial strength of the wavelets originating on the sphere's surface). For example, [1/R 0 In − δ(t − [z 0 + 2R 0 ]/c) the inverse spreading due to going from the far surface to the center is exactly cancelled by the spreading resulting from going from the center to the near surface. So the same R 0 /[R 0 + z 0 ] spherical spreading term can ultimately be applied to both the near and far surface waves that are to be observed at P.
The wave field θ(ct − z 0 − 2R 0 ) originated as a positive amplitude propagating wave field on the far surface (in the limit as R 0 → ∞ the wave fields close to the surface of the sphere must be the same as the wave fields from the plane, and the wave fields on both sides of the plane are positive). As the initially converging spherical wave field passes through the center its amplitude changes sign (Gouy phase shift [39][40][41] : "It is well known that a spherical converging light wave undergoes a phase change of 180 degrees in passing through its focus…and is, in fact. a general property of any focused wave.", Boyd p. 877) 40 , and see Fig. 6d. www.nature.com/scientificreports/ The wave field then begins expanding. When it expands through the surface of the spherical source the now negative wave field cancels the remainder of the wave field that originated on the near surface, θ(ct − z 0 ), (Fig. 6c). The combined wave field (Fig. 6e) when observed in time at P is a rectangular pulse of time duration 2R 0 /c . Notice that only the two points N and F (Fig. 6a), contribute as sources to the wave front f(z 0 ,t). This is because φ(z 0 , t) is a rectangular pulse of temporal width 2R 0 /c for fixed z 0 and its derivative is zero except on its leading and trailing edges (which originated at the two points N and F). Because the derivative is zero there is no wake (except for an impulse corresponding to the trailing edge). Also there is a one to one correspondence between the points on the source and the points on the wave front.
Therefore a spherical distribution of point sources radiating impulsive spherical wavelets creates two impulsive spherical wave fronts (Fig. 6b). Both propagate cleanly without a wake. Again this shows that the Huygens' wavelets cancel each other when summed together and differentiated, except where they make contact with their tangent surfaces.
The same 1/2 factor occurs in f as occurred in the case of the stationary infinite planar source. Note also there are forward (expanding) and backward (converging) wave fronts. When the sphere's radius expands with a speed equal to the propagation speed, the 1/2 factor and the backward wave both disappear, as will be shown.
Expanding spherical source and elimination of the backward wave. This section analyzes the changes that occur in the results of section "Non expanding spherical source and elimination of the wake" when the spherical source is expanding with radial speed v. It is shown that the results for the non expanding sphere in section "Non expanding spherical source and elimination of the wake" can be used if the effects of motion on the Huygens' wavelets are included.
When this is done it is seen that a new additional wave is created. (A similar additional wave is discussed in section "Moving planar source and elimination of the backward wave".) The amplitude of this additional wave is a function of the radial expansion speed, v. This additional wave reduces the amplitude of the backward wave and increases the amplitude of the forward wave. When v = c the backward wave is eliminated and the forward wave amplitude is doubled. In any case the original wave is replicated and there is no wake. See Fig. 7. The expanding equivalent, As v , of the nonexpanding spherical zone area, As, must be derived. This is the lateral area of the frustum.
The frustum area is where the slant height is and the height is The height h(t) was found previously as that plus the two right triangles T1 and T2 can be solved for h(t + ε): Expanding the terms on the left side of T1-T2 gets: The impulsive spherical shell source (the 'sphere') with initial radius R 0 expanding with a positive relative radial speed v where as with the planar source, v will be positive for forward wave fronts (in the direction of the expansion) and will be proceeded by a negative sign for backward wave fronts. R1 is the radius of the circle created by the intersection at time t of a sphere centered at P of radius ct with the spherical source. R2 is the radius of the circle created by the intersection at time t + ε of a sphere centered at P of radius c[t + ε] with the spherical source. www.nature.com/scientificreports/ Expanding the two terms on the right side of T1-T2 yields for the first term: and yields for the second term: Simplifying and subtracting the two terms, yields for the right side of T1-T2: Equating the left side to the right side of T1-T2 produces Solving for h(t + ε) results in: Then and squaring both sides and taking a limit gets So the frustum area as ε → 0 is which is also the area, Ap, of an annulus with radii R1 and R2. This annulus area Ap = 2πc 2 tε was previously computed for the stationary planar source and was also used in determining the area of the nonexpanding spherical source, As = [R 0 /[R 0 + z 0 ]]2πc 2 tε. Consequently the previously developed results for the non expanding sphere can be used if the effects of the motion of the expanding spherical source on its wavelet point sources' initial impulses are included: At the wavelet point source location as z 0 → 0 the motion due to the sphere's expanding radius is the same as the motion of the moving planar source if spherical spreading is disregarded. Then the method used for the moving planar source will be used here to implement the effects of motion on the wavelet point sources initial impulses due to the sphere's radius expanding with a speed v: As with the non expanding sphere, the wavelet point sources at points N and F of the near and far surfaces of the expanding sphere are the sole contributors to the wave front along the radial line. The wavelet point sources at points N and F on the near and far surfaces have relative speeds + v and − v with respect to their originating surfaces. (The far surface is the source for a backward propagating wave front and so has relative speed − v.) Including the moving planar source's source speed effect in the equation for the wave front from a nonexpanding spherical source observed at P as a function of time yields: for the forward propagating spherical wave front from the near surface combined with the backward propagating spherical wave front from the far surface of the expanding spherical source.
When source speed equals zero, v = 0, the result is the same as for the non expanding spherical source. When source speed equals the propagation speed, v = c, the equation becomes www.nature.com/scientificreports/ the ½ factor). Since this holds for every radial line, this is the new impulsive spherical progressing Huygens' wave front which has neither a wake nor a backward wave. The spherical spreading factor R 0 R 0 +z 0 differentiates this result from the result from the moving planar source and accounts for the radial lines representing infinitesimal solid angles. As R 0 → ∞ this equation becomes the same as for the planar source moving with speed equal to c.
Huygens' Principle as depicted in his figures. It was shown that for the spherical and planar source, when the excitation is δ(t), the resulting wave field ϕ when differentiated disappears except at the locations of the propagating wave fronts f. The following shows how the intermediate step of differentiation can be eliminated.
Let L represent the linear wave propagation operator which transforms the excitation on either the plane or sphere into the wave field ϕ, If δ(t) is the excitation, the resulting wave field due to the wave propagation process is which is the response of L to an impulse. As a function of time in the case of the plane this wave field is a step function and in the case of the sphere it is a rectangular pulse.
Because L is linear (LTI) the derivative can be taken of both ϕ and of δ(t) 42 yielding But the derivative of the wave field ϕ is the wave front f, so As a result, if the doublet δ′(t) is used as the initial excitation instead of the singlet δ(t) for the wavelets originating from the source, the Huygens' wavelets when combined disappear except where they contact their common tangent surfaces. Also when the sources' motion equals the wave propagation speed the backward wave corresponding to the backward tangent surface disappears. The remaining forward tangent surface is the new propagating wave front, obtained without the intermediate step of differentiation.
In Huygens' figures the wavelets then would represent doublets [4πr] −1 δ′(t − r/c). However, notice that the resulting new wave fronts are of the form δ not δ′ so they must be differentiated before being used as the initial excitation in the next iteration of wave propagation.
Because ϕ is the response to an impulse, the analysis can be applied to other forms of excitation, say g(t), by using convolution: the convolution gives the response of L to the excitation g.

Summary
Huygens' geometrical construction has been shown to be literally correct as depicted in his figures (see Fig. 1a) when his wavelets represent the first derivative of Dirac delta distributions and the source motion is equal to the propagation speed. For the sphere and the infinite plane there is no wake and there is no backward wave, thus solving the two most long standing problems (1678).
An intuitive, relatively simple method to analyze wave propagation has been provided which provides a clear view of the propagation process starting with the very first wave propagation analysis, Huygens' Principle. This will facilitate the teaching of wave physics and provide a new model for wave propagation.