Chapter 5

Causality and Coupling Delay

In the analysis presented in the preceding chapters, it was tacitly assumed that the coupling was implemented using nondispersive transmission lines characterized by a phase shift of Φ generally taken to be an integral multiple of 2π (plus π in the case of series resonant oscillators). However, the theory made no provision for the transit time through the coupling line. As a result, the solutions were noncausal. That is, each oscillator in the array responded immediately upon changing the tuning of an oscillator or the phase of an injection signal no matter what the distance between the excitation and the response. This is characteristic of the diffusion equation that arises from the continuum model. Heat conduction analyzed in this manner is similarly non-causal. Following Pogorzelski [47], we propose to remedy this situation by explicitly introducing time delay in the coupling. This time delay is determined by the physical length of the line and its propagation velocity.

5.1   Coupling Delay

A nondispersive transmission line introduces a pure time delay in that the signal applied at one end of the line is duplicated at the other end after the delay time. At that point the signal is reflected if the termination is not matched to the line impedance. For our analysis we will assume a matched termination. Now, if the analysis is done via Laplace transformation of the applied signal, the transform of the delayed signal is merely the original transform multiplied by e^−sd, where d is the delay time and s is the transform variable conjugate to the time variable.

Suppose now that we envision an array of coupled oscillators and take the reference frequency to be the initial ensemble frequency of the array. We define the coupling phase delay using this reference frequency. That is,

where ℓ is the physical length of the line and ν_p is the phase velocity. The line length is chosen so that the coupling phase is a multiple of 2π (plus π in the case of series resonant oscillators). Now, using the reference frequency, we define the phase, φ, of the phasor signal voltage, V, by

where

Recall from Chapter 1 that V can be written in the form

so that

Crucial to our analysis is the fact that any function of the input signal will be delayed by the nondispersive transmission line in the same manner as the signal itself so that the Laplace transform of any function of the input signal multiplied by e^−sd will be the transform of the same function delayed. Thus, we may apply this delay factor to the Laplace transform of φ(t) given by Eq. (5.1-5) to obtain the transform of the phase delayed by the coupling line. This forms the basis of our introduction of coupling delay into the analysis of coupled oscillator arrays.

The following question regarding this treatment of time delay was posed by a particularly astute student, so we thought it appropriate to answer it here as you may be similarly puzzled. Slightly paraphrased, the student asked that we consider a linear array in which one of the oscillators is detuned upward, thus changing the ensemble frequency of the array. “Is it not then true,” he asked, “that the coupling phase produced by coupling lines of fixed length would be changed and would thus be no longer a multiple of 2π?” To clarify this point, recall that, assuming that the reference frequency is held constant, the solution for the time evolution of the oscillator phases will contain terms linear in time representing the shift in ensemble frequency (as in Eq. (3.1-51)). The slope of this linear dependence relates the time delay to an equivalent phase shift through the transmission line. So, for example, if the phase at the input end of the line is

then the phase at the output end of the line is

an effective coupling phase delay of Φ–αd. Conversely, due to the linear time dependence, the new ensemble frequency will be

and the effective coupling phase will be

So, we conclude that indeed the coupling phase has changed, but that change is embodied in the linear time dependence of the phases arising from the change in ensemble frequency and need not be explicitly imposed on the formulation by a change in the Φ parameter.

5.2   The Discrete Model with Coupling Delay

Returning to the linearized discrete model of a linear array of (2N + 1) oscillators discussed in Section 2.2, we have

Laplace transformation with respect to results in

and, introducing the coupling delay factors for delay of d inverse locking ranges, we have

Rearranging yields

These equations may be written compactly in matrix form as

in which [M] is given by

We now have two alternative approaches available for solving this system of linear equations. We can expand the solution as a sum of eigenvectors of the matrix [I]s − [M], or we can solve the system via Cramer’s rule. Following Pogorzelski [47], we choose the Cramer’s rule approach. The result is

where U is the Chebyshev polynomial of the second kind of argument a/(2b), a₀ = s + 1, a = s + 2, and b = −e^−sd. Now, U can be written in the form

and defining Q to be

U becomes

Substituting Eq. (5.2-18) into Eq. (5.2-15) yields

in which

In Eq. (5.2-19), n_> is the greater of i and j while n_< is the lesser. The form of Eq. (5.2-19) is suggestive of an image series produced by reflections at the ends of the array, where Q plays the role of a “propagator.” The series may be obtained by expanding in powers of the reflection coefficient, R. When R is set equal to zero, we obtain the solution for an infinite array,

or, using Eq. (5.2-17),

Consider now an example of a 17-element array with coupling delay of two inverse locking ranges (ILRs) and step detuning of the center oscillator by one locking range. The inverse Laplace transform of Eq. (5.2-19) may be easily obtained by expanding it in powers of −b = e^−sd. Each term of the resulting series will be of the form

which has a known inverse transform

The solution is plotted in Fig. 5-1. This solution exhibits several easily understandable features. First, the center oscillator is the only one detuned, and it is detuned at time zero. Thus, its nearest neighbors on either side do not change phase until one delay time has elapsed, giving the influence of the center oscillator detuning sufficient time to propagate to them. The center oscillator phase evolution continues exponentially and unperturbed until two delay times have elapsed. Then the influence of the phase changes of the neighboring oscillators impact the center oscillator, causing the slope change at time equal to four inverse locking ranges. This multiple reflection-like behavior continues to spread throughout the array creating the ripples in the phase visible in Fig. 5-1. For comparison, the corresponding phase behavior in the absence of coupling delay is shown in Fig. 5-2 and is noncausal.

The preceding discussion indicates that the array behavior will not begin to differ from that of the infinite array until eight delay times have elapsed (16 inverse locking ranges) and that even then the effect will begin with the outermost oscillators. The center oscillator behavior will not differ from that of the infinite array until 16 delay times have elapsed (32 inverse locking ranges). Thus, these plots do not indicate that the array is of finite size. To display finite array effects, we plot similar curves for a seven element array in Fig. 5-3 and Fig. 5-4. In this case the end effects begin to appear in the phase behavior of the end oscillators after three delay times (six inverse locking ranges). This may be seen by comparing the curves for i = ±3 in Figs. 5-1 and 5-3. Note that they differ only for time greater than six inverse locking ranges. However, the curves in Figs. 5-2 and 5-4 differ for all times because with no coupling delay the end effects begin immediately and, of course, acausally.

Fig. 5-1. Phase dynamics for a 17-element linear array with two-inverse-locking-range coupling delay. (Reproduced by permission of American Geophysical Union from Ref. 47, ©2008 American Geophysical Union.)

Fig. 5-2. Phase dynamics for a 17-element linear array with no coupling delay. (Reproduced by permission of American Geophysical Union from Ref. 47, ©2008 American Geophysical Union.)

Thus, we have shown that the introduction of coupling delay in the linearized discrete model of coupled oscillator arrays eliminates the noncausal nature of the solutions in the absence of coupling delay. We now proceed to apply this approach in the continuum model.

5.3   The Continuum Model with Coupling Delay

In this section we develop a generalization of the continuum model of Section 3.1 that accounts for coupling delay. The causality properties of this generalization will be discussed in terms of the infinite linear array which, of course, behaves identically to a corresponding finite array for times early enough to preclude end effects. Our approach will be that of Section 3.1, where we begin with Eq. (3.1-1) with the time delay, d, inserted in the arguments of the appropriate terms leading to

Fig. 5-3. Phase dynamics for a 7-element linear array with two-inverse-locking-range coupling delay. (Reproduced by permission of American Geophysical Union from Ref. 47, ©2008 American Geophysical Union.)

Fig. 5-4. Phase dynamics for a 7-element linear array with no coupling delay. (Reproduced by permission of American Geophysical Union from Ref. 47, ©2008 American Geophysical Union.)

Introducing the scaled time, τ, and the detuning function, ΔΩ_tune, as before, Laplace transformation leads to

Then, expanding in Taylor series to second order in Δx, we obtain

the analog of Eq. (3.1-4). Setting

corresponding to step detuning of the oscillator at x = y at time zero by one locking range, we obtain the Green’s function, ₁, as the differential equation solution,

At this point, a serious difficulty is encountered with respect to causality. If one were to compute numerically the inverse Laplace transform integral for Eq. (5.3-5), one would find that the influence of the nearest neighbors of the detuned oscillator begins at time d. This violates causality because, as pointed out in Section 5.2, this influence must not begin until time 2d, the round-trip transit time between the detuned oscillator and its neighbors. Following Pogorzelski [48], we begin our study of this apparent paradox by comparing the denominator of Eq. (5.3-5) with that of Eq. (5.2-22) known to be causal. That is, the denominator of Eq. (5.2-22) is

while the denominator of Eq. (5.3-5) is

Thus, the two solutions—the causal one Eq. (5.2-22), and the present one, Eq. (5.3-5)—have different branch points in the complex s plane. Solution Eq. (5.3-5) has branch points where s + 2 = 2e^−sd, whereas the causal solution Eq. (5.2-22) has these plus additional branch points where s + 2 = −2e^−sd. Now, computing the inverse Laplace transform via integration on the Bromwich contour will involve deformation of the contour around the branch cuts associated with these branch points. Thus, it becomes clear that the solution Eq. (5.3-5) will be missing the contribution from half of the branch cuts in the causal solution Eq. (5.2-22). As shown in Ref. 48, this is the root of the causality difficulty.

Why do we find ourselves in this situation? Our approach was successful in the absence of coupling delay, but something went wrong when delay was included. This can be understood by looking at the nature of the solutions corresponding to the two sets of branch cuts shown in Fig. 5-5 where the dots correspond to Eq. (5.3-7) and the circles to the remaining branch points of the complete set, Eq. (5.3-6).

In deriving the partial differential equation, Eq. (5.3-3), we have assumed that the solution will be smoothly varying in the interior of the array so that the inter-oscillator phase differences are small, validating the linearization of the sine functions in Adler’s formalism. Thus, in the interior of the array where the detuning is zero, the second derivative will be small and

corresponding to the dot branch points in Fig. 5-5. However, we can switch from the dots to the circles by replacing e^−sd with −e^−sd. Doing this in Eq. (5.3-2), we obtain

Fig. 5-5. Branch point locations for delay of two inverse locking ranges. (Reproduced by permission of American Geophysical Union from Ref. 48, ©2008 American Geophysical Union.)

so it is evident that φ₁, the solution associated with the circle-type branch points, alternates in sign between adjacent oscillators and is thus clearly not slowly varying. Therefore, we cannot use the Taylor expansion to advantage here. However, if we define

then is slowly varying even though is not and we may write

and expand in Taylor series to obtain

corresponding to the circle-type branch points. Thus, it becomes clear that our assumption of slowly varying phase, implicit in the use of the Taylor series, eliminated the solutions associated with the circle-type branch points. The Green’s function corresponding to these branch points is

and the causal Green’s function is a linear combination of Eq. (5.3-5) and Eq. (5.3-13); that is,

where A + B = 1 so that the proper detuning function is generated on the right side of the differential equation. It remains to determine A and B.

For large values of s, far from the origin of the s plane, we can obtain a fairly accurate estimate of the branch point locations. These locations are defined by

where the upper sign corresponds to the dots and the lower one to the circles. Inserting s = σ + jω, we obtain

For |ω| >> σ + 2, we have

Thus,

and we have

for the dots and

for the circles. Now, from Eq. (5.3-17), we obtain

so

for p = m or n. Armed with these approximate branch point locations, we are in a position to estimate g(τ, y, y)/τ, the time derivative of the phase of the detuned oscillator. This will exhibit the temporal discontinuities associated with the arrival of influence from neighboring oscillators and highlight the causal behavior. First, from Eq. (5.3-14), we have

Now, envisioning the inverse transform as a sum of branch cut integrals, we recognize that the result will be approximately

Here we have used the s ’s given by Eq. (5.3-19) through Eq. (5.3-22) and

The expression on the right side of Eq. (5.3-24) is a Fourier series except for the time dependence of the coefficients. Recall that this series was obtained using the large s approximation, so only the high-order terms are accurate. The high-order terms of this series govern the discontinuities in the time dependence. Now, looking at Eq. (5.3-24) for τ = d, we obtain

and we see that if A = B, the high-order portions of the two series will cancel term by term so that there will be no discontinuity at τ = d. However, at τ = 2d we have

and the high-order terms no longer cancel but add. Thus, there will be a discontinuity at τ = 2d. This is to be expected because it allows for one round-trip interval to the nearest neighbors from the time when the oscillator is detuned. We conclude that the discontinuities will occur at the proper times for causality to be satisfied only if A = B. From this condition and the fact that A + B = 1, we determine that both A and B are equal to , and from Eq. (5.3-14) the causal Green’s function is

As shown in Ref. 48, a better approximation to the exact discrete model solution may be obtained from the form,

with optimal selection of the constant, C. From Ref. 48, the optimal value of C is

The temporal behavior of the phase of each oscillator in the array is most easily seen by plotting the time derivative of the phase because this makes more obvious the times at which the influences from the neighboring oscillators arrive. Thus, in Figs. 5-6 through 5-10, we compare the result of the approximate continuum formula Eq. (5.3-29) in solid lines with that of the discrete model Eq. (5.2-22) in dashed lines considered to be the exact result. The coupling delay in this example is two inverse locking ranges (ILRs). The fine scale wiggles shown in the inset of Fig. 5-6 arise from the truncation of the series of branch cut integrals to a finite number of terms.

Fig. 5-6. Time derivative of the phase of the detuned oscillator. (Reproduced by permission of American Geophysical Union from Ref. 48, ©2008 American Geophysical Union.)

Fig. 5-7. Time derivative of the phase of the nearest neighbors of the detuned oscillator. (Reproduced by permission of American Geophysical Union from Ref. 48, ©2008 American Geophysical Union.)

Fig. 5-8. Time derivative of the phase of the second nearest neighbors of the detuned oscillator. (Reproduced by permission of American Geophysical Union from Ref. 48, ©2008 American Geophysical Union.)

Fig. 5-9. Time derivative of the phase of the third nearest neighbors of the detuned oscillator. (Reproduced by permission of American Geophysical Union from Ref. 48, ©2008 American Geophysical Union.)

Fig. 5-10. Time derivative of the phase of the fourth nearest neighbors of the detuned oscillator. (Reproduced by permission of American Geophysical Union from Ref. 48, ©2008 American Geophysical Union.)

Notice that the more distant the oscillator from the detuned one, the later the response by exactly two inverse locking ranges (one delay time) per oscillator. Moreover, the influence of the nearest neighbors of the detuned oscillator does not impact that oscillator until four inverse locking ranges (two delay times) have elapsed. Similar delays of two delay times are visible in all of the curves corresponding to round-trip delays between the oscillators. All of these behaviors are consistent with a causal solution.

We return now to the previous analysis of the location of the branch points to highlight two properties that may not have been obvious in the earlier discussion. First, as the delay time is decreased, there is a critical value at which the distribution of the branch points changes character. If the delay time is equal to 0.139232271 inverse locking ranges, the smallest circle-type branch points merge at σ = −9.18224297. For delays less than that, say for a delay of 0.12 inverse locking ranges, the branch point locations are as shown in Fig. 5-11.

Second, as the delay approaches zero, all of the branch points move to infinity except two: one at the origin and one at −4. Thus, in this zero-delay limit we have from Eq. (5.2-22), taken to be the exact solution, that

which, perhaps surprisingly, does not agree with Eq. (3.1-6). It does agree in the limit of small s, so one can expect that the time functions will agree for late times, but there will be a difference at early times. When x = y; that is, for the detuned oscillator, the inverse Laplace transforms of Eqs. (3.1-6) and (5.3-31) can be computed analytically, and we thus obtain from Eq. (5.3-31)

where I_n is the Bessel function of imaginary argument and from Eq. (3.1-6), we get

Fig. 5-11. Branch point locations for d = 0.12 inverse locking ranges. (Reproduced by permission of American Geophysical Union from Ref. 48, ©2008 American Geophysical Union.)

For comparison, these two functions are plotted in Fig. 5-12. The solid curve is Eq. (5.3-32), and the short dashed curve is Eq. (5.3-33) while the long dashed curve is the difference. Note that, although Eq. (5.3-33) neglects the alternating sign solution, it is nevertheless a very good approximation to the exact solution Eq. (5.3-32).

5.4   Beam Steering in the Continuum Model with Coupling Delay

In this section we apply what we have learned so far regarding the analytical treatment of coupling delay to the analysis of beam-steering of oscillator arrays embodying such delay. We begin with the continuum generalization of the linearized discrete model solution for one detuned oscillator in an infinite array given by Eq. (5.2-22); that is,

Fig. 5-12. Comparison of the single-branch-point solution (dashed) with the two-branch-point solution (solid). The difference is shown in long dashes. (Reproduced by permission of American Geophysical Union from Ref. 48, ©2008 American Geophysical Union.)

Our approach will be to devise a differential equation having Eq. (5.4-1) as its Green’s function. When x is not equal to y, this solution satisfies the differential equation,

However, at x = y there will be a discontinuity in the slope of the phase that gives rise to a delta function. Evaluating the magnitude of this slope discontinuity, we determine that the Green’s function Eq. (5.4-1) satisfies

We may now follow the procedure of Section 3.1 to express the finite-array Green’s function as a sum of the eigenfunctions of the differential operator in this equation. In order to do this, we will need the boundary conditions at the ends of the array. Recall that the reflection coefficient at the array ends was given by Eq. (5.2-20), which is a fairly complicated function of s. However, following Pogorzelski [49], we may simplify matters by assuming the addition of half-length coupling lines at the ends of the array. If this is done, the reflection coefficient becomes unity because the array boundary then becomes an image plane. (See Pogorzelski [47].) A reflection coefficient of unity corresponds to the familiar Neumann condition of zero phase slope. Using this boundary condition, the even and odd normalized eigenfunctions are seen to be

Choosing the detuning time dependence to be a unit step at time zero and following he approach of Section 3.1, the Laplace transform of the phase distribution may be written in terms of the eigenfunctions as

We will obtain the inverse Laplace transform via residue calculus. The poles are determined by

Taking the hyperbolic secant of both sides and then the reciprocal, we obtain the equivalent condition

This equation can be solved in terms of the Lambert W function defined by

In terms of this function, the solution of Eq. (5.4-7) is

These pole locations are plotted in Fig. 5-13.

Fig. 5-13. Pole locations for delay of two inverse locking ranges. Black dots denote odd values of ℓ and gray dots denote even values of ℓ. (Reproduced by permission of American Geophysical Union from Ref. 49, ©2008 American Geophysical Union.)

The overall array time constant is determined by the poles closest to the origin. We therefore set about solving Eq. (5.4-9) approximately for small s. To do this we expand the Lambert W function in a Taylor series about 2de^2d.

Now, W satisfies the differential equation

Thus, the first two terms of the Taylor series yield

and, using Eq. (5.4-9), we obtain

Recall that z is the argument of the W function, so from Eq. (5.4-9) we have

Substituting Eq. (5.4-14) into Eq. (5.4-13) and setting m = 0, we obtain

The pole at ℓ = 0 together with the denominator s from the step detuning function produce the double pole at the origin leading to the linear time dependence or shift in ensemble frequency due to the detuning. For the antisymmetric detuning used in beam-steering, the even ℓ poles do not contribute, so the dominant pole is the one for ℓ = 1 lying on the real axis at

so the time constant of the array is

or just (2d + 1) times the time constant without coupling delay. (Compare with Eq. (2.2-40).)

Returning now to Eq. (5.4-5), we form the solution for beam-steering by combining two solutions of the form Eq. (5.4-5), one for detuning of the oscillator at −N and one for detuning of the oscillator at N, at each end of the array.

To obtain the residues, we define q(s) to be the denominator,

and expand in the Taylor series,

We thus obtain the aperture phase a function of time in the form

Recall that at each oscillator, x is an integer and that the phase only has physical meaning at these integral values of x. As a result, the sum on n need only extend from 0 to N − 1 because for integral x, these terms are equal to those for n = N + 1 through 2N with the order reversed. Higher-order terms in n only affect the phase values between the oscillators and thus are not relevant.

Typically, the time at which the time function becomes nonzero is determined by when the Bromwich contour used in the inverse Laplace transform integral can be closed in the left half-plane, thus enclosing the poles. Prior to that time, the contour may only be closed in the right half-plane; and since it encloses no poles there, the solution is zero. This is the usual way in which causality enters such analysis. In this case, however, Eq. (5.4-18) reveals that the contour may be closed in the left half-plane beginning at τ = −d. Causality dictates that the solution remain zero until τ = 0. Thus, it turns out that the residue sum remains zero even though the contour is closed in the left half-plane and only becomes nonzero after τ = 0. This is illustrated in Eq. (5.4-21), in which the contour was closed in the left half-plane beginning at τ = −d/4 leading to the term. Figure 5-14 shows the resulting solution for each oscillator of a 21-element array with coupling delay of two inverse locking ranges. Causality is obviously satisfied regardless of this unusual closing of the contour. Figure 5-15 shows the same solution extending to later times showing that in steady state the phase increments between oscillators become equal, implying a linear phase progression as needed for beam-steering.

The data in Figs. 5-14 and 5-15 are re-plotted in Figs. 5-16 and 5-17, respectively. Here one may view the aperture phase distribution at all values of time simultaneously.

Finally Fig. 5-18 shows a particular range of time specifically for comparison with Fig. 5-19, which is the same case but with no coupling delay. Pay particular attention to the time scales in these plots.

The point made by comparing Fig. 5-18 with Fig. 5-19 is that the coupling delay of two inverse locking ranges has slowed the response of the array by 2d + 1 or a factor of five, just as predicted by Eq. (5.4-17).

We now compute the far-zone radiated field when the oscillators in this 21-element array with coupling delay are used to excite the elements of a phased array with half-wavelength element spacing. The result is shown in Fig. 5-20 where we see beam behavior very similar to that of arrays without delay but slower by (2d + 1).

Fig. 5-14. Phase evolution of the oscillators in a 21-element linear array with coupling delay of two inverse locking ranges. (Reproduced by permission of American Geophysical Union from Ref. 49, ©2008 American Geophysical Union.)

Fig. 5-15. Phase evolution of the oscillators in a 21-element linear array with coupling delay of two inverse locking ranges over a longer duration. (Reproduced by permission of American Geophysical Union from Ref. 49, ©2008 American Geophysical Union.)

Fig. 5-16. Three-dimensional representation of the phase evolution in a 21-element array at early times. (Reproduced by permission of American Geophysical Union from Ref. 49, ©2008 American Geophysical Union.)

Fig. 5-17. Three-dimensional representation of the phase evolution in a 21-element array at later times. (Reproduced by permission of American Geophysical Union from Ref. 49, ©2008 American Geophysical Union.)

Fig. 5-18. Phase evolution over 250 inverse locking ranges for a 21-element array with coupling delay of two inverse locking ranges. (The vertical scale is from −10 to 10 radians as in Fig. 5-19.) (Reproduced by permission of American Geophysical Union from Ref. 49, ©2008 American Geophysical Union.)

Fig. 5-19. Phase evolution over 250 inverse locking ranges for a 21-element array with no coupling delay. (Reprinted from Ref. 38 with permission, ©2000 IEEE.)

Fig. 5-20. Dynamic behavior of the far-zone radiated field for a 21-element array with coupling delay of two inverse locking ranges. (Reproduced by permission of American Geophysical Union from Ref. 49, ©2008 American Geophysical Union.)

The final value theorem applied to Eq. (5.4-18) gives the steady-state phase distribution as

and for half-wavelength spacing of the radiating elements, we find that the steady-state beam position is at

This one-radian inter-oscillator phase difference stresses the linear approximation a bit in that the error in linearizing the sine function is about 19 percent. However, the dynamic behavior is still qualitatively approximated.

In the above analysis a large number of residues are required for early times and very few are required for late times. However, returning to the discrete model, an alternative formulation is available that provides for more efficient computation for early times. Returning to Eq. (5.2-19) and specializing to the present case of a 21-element array, we have that

where, as before, a₀ = s + 1, a = s + 2, and b = −e^−sd. We now expand this expression in powers of −b, and, as before, the inverse Laplace transform of each term in the expansion can be computed analytically. The number of terms required is determined by the time interval over which the response is desired because each term has a delay factor e^−psd, where p is the power of −b in the term in the expansion and d is the coupling delay. So, for sufficiently large p, the term will be zero for the interval in question. Thus, in contrast with the eigenfunction expansion, for early times very few terms are required.

This approach was applied to the 21-element array with coupling delay treated earlier, and the results are plotted in Fig. 5-21 for comparison with Fig. 5-14. Interestingly, this power series approach is a bit more flexible in terms of boundary conditions. Recall that without the added half-length coupling lines at each end of the array, the previous method was complicated. Here, however, the use of Eq. (5.2-20) in Eq. (5.2-19) to model an array without the added lines poses no difficulty. The expansion in powers of b proceeds as before and the result is plotted in Fig. 5-22. Notice the difference in the early time ripples due to this alternative boundary condition when compared with Fig. 5-21.

Fig. 5-21. Phase evolution in a 21-element array with coupling delay of two inverse locking ranges via expansion in powers of b. (Reproduced by permission of American Geophysical Union from Ref. 49, ©2008 American Geophysical Union.)

Fig. 5-22. Phase evolution via expansion in powers of b for a 21-element array with coupling delay of two inverse locking ranges but without the added half-length coupling lines at the ends. (Reproduced by permission of American Geophysical Union from Ref. 49, ©2008 American Geophysical Union.)

5.5   Conclusion

The primary motivation for this chapter was the issue of causality in coupled oscillator arrays. Because the continuum model leads to a diffusion equation, the response to an excitation always begins immediately regardless of the physical separation of the two. Here, by appropriately introducing a delay factor in the Laplace transforms, we render the solutions causal in that there appears a finite “propagation delay” between the excitation and the response. The result is a more realistic representation of the array response, not to mention some rather interesting inverse Laplace transforms encountered along the way.