text:
# The proximity effect in normal metal - multiband superconductor hybrid structures

## Abstract

A theory of the proximity effect in normal metal - multiband superconductor hybrid structures is formulated within the quasiclassical Green’s function formalism. It is shown that the existence of multiple superconducting bands manifests itself as the occurence of additional peaks in the density of states in the structure. The interplay between the proximity effect and the interband coupling influences the magnitudes of the gaps in a superconductor in a non-trivial way. The developed theory is applied to the calculation of supercurrent in multiband Superconductor - Normal metal - Superconductor Josephson junctions with low-transparent interfaces, and the results are compared with the predictions for multiband tunnel junctions.

The proximity effect is the phenomenon that a superconducting order parameter can penetrate from a superconductor (S) into a normal metal (N), or another superconductor (S’) with a critical temperature $`T_{cS^{}}<T_{cS}`$, over a distance of the order of the coherence length, inducing a minigap in N or S. This phenomenon is well understood, both in terms of Andreev reflections as well as in terms of microscopic Green’s functions Cooper ; Werth ; deGennes ; McMillan ; KL ; Golubov ; Fominov ; Beenakker ; Belzig .

It is not known, however, how the proximity effect will manifest itself when multiple pairing potentials are present in the superconductor. This question has become relevant now that multiband superconductors are coming into practical use. The most clear example of a multiband superconductor is MgB2, for which the experimental and theoretical evidence for the coexistence of two gaps is overwhelming issue . The multiband nature of the superconductivity in MgB2 is theoretically well explained Liu  by the qualitative difference between different sheets of the Fermi surface, together with the large disparity of the electron-phonon interaction. Therefore, in this paper, the question is addressed how the multiband nature influences the proximity effect. E.g. what will the density of states look like in a SN bilayer, where S is a two-band superconductor?

Josephson and quasiparticle tunneling in hybrid structures containing multiband superconductors have been investigated theoretically in Ref. Brinkman,  and applied to the calculation of the total Josephson current in a SIS two-band Josephson tunnel junction. For all-MgB2 devices, high-quality tunnel barriers are not available yet, and realizing SNS structures is an attractive alternative, of which first systems have been realized already SNS . In this paper, the theory of the multiband proximity effect is applied to the calculation of Josephson current in SNS structures having two-band S electrodes. The practically interesting SINIS case is considered, where a non-ideal interface transparency is taken into account. Predictions are made for Josephson devices based on MgB2 and compared with those for MgB2-based tunnel junctions.

In this paper, we will use the quasiclassical Green’s function formalism in order to describe electrical transport in SS’ hybrid structures, where S’ is a single-band superconductor while S is a multiband superconductor. We will restrict ourselves to the limit of diffusive transport, which is justified if $`l_{S,S^{}}\xi _{S,S^{}}`$, where $`l_{S,S^{}}`$ and $`\xi _{S,S^{}}`$ are the electric mean free path and coherence length of the S and S’ materials respectively. In the dirty limit, the Green’s functions in the S’ metal are given by the standard Usadel equations Usadel . In the S metal in the regime of vanishing interband scattering, as is the case for MgB2 cleandirty , the Usadel equations take the following form Koshelev

$`{\displaystyle \frac{D_S^\alpha }{2\omega G_S^\alpha }}\left[(G_S^\alpha )^2\mathrm{\Phi }_S^\alpha \right]^{}\mathrm{\Phi }_S^\alpha =\mathrm{\Delta }_\alpha ,`$ (1)
$`\mathrm{\Delta }_\alpha =2\pi T{\displaystyle \underset{\beta ,\omega 0}{}}\widehat{\mathrm{\Lambda }}_{\alpha \beta }{\displaystyle \frac{G_S^\beta \mathrm{\Phi }_S^\beta }{\omega }}.`$ (2)
Here, $`\alpha `$ and $`\beta `$ are the band indices, e.g. $`\alpha ,\beta =1,2`$ in the two-band case (later we will use the band indices $`\sigma `$ and $`\pi `$ for MgB2 specifically), $`\mathrm{\Delta }_\alpha `$ is the pair potential, $`G_S^\alpha `$ and $`\mathrm{\Phi }_S^\alpha `$ are Green’s functions KL , $`\omega =\pi T(2n+1)`$ are Matsubara frequencies, $`D_S^\alpha `$ is the diffusion coefficient, and $`\widehat{\mathrm{\Lambda }}_{\alpha \beta }`$ is the matrix of effective coupling constants. The prime denotes a derivative with respect to the coordinate $`x`$ in the direction perpendicular to the S-S’ interface.

Equations (1) and (2) must in general be supplemented by boundary conditions. Zaitsev Zai1  derived boundary conditions to the quasiclassical Eilenberger equations at the S-S’ boundaries in the clean limit, which were further simplified in Ref. KL,  in the dirty limit. These boundary conditions have to be modified when S is a multiband superconductor.

In the limit of small interband scattering a multiband superconductor may be represented by separate groups of superconducting electrons which interact with each other only indirectly, via selfconsistent pair potentials in the bulk. Therefore, for the derivation of the boundary conditions for the Usadel equations, one can apply a similar procedure to that used in Ref. KL,  in the single-band case. In the multiband case, the set of interface parameters, $`\gamma ^\alpha `$ and $`\gamma _B^\alpha `$, describing the proximity effect, should be introduced for each of the bands.

The first boundary condition relates the current from the S’ metal side at the S-S’ interface, $`\sigma G^2\mathrm{\Phi }^{}`$, to that from the S side, $`_\alpha \sigma _S^\alpha (G_S^\alpha )^2(\mathrm{\Phi }_S^\alpha )^{}`$. Therefore we have

$$\xi G^2\mathrm{\Phi }^{}=\underset{\alpha }{}\frac{\xi _S^\alpha }{\gamma ^\alpha }(G_S^\alpha )^2(\mathrm{\Phi }_S^\alpha )^{},$$
(3)
with

$$\gamma ^\alpha =\frac{\rho _S^\alpha \xi _S^\alpha }{\rho \xi },(\xi _S^\alpha )^2=\frac{D_S^\alpha }{2\pi T_{cS}},\xi ^2=\frac{D}{2\pi T_{cS}}$$
(4)
Here, $`\sigma =1/\rho `$ and $`\sigma _S^\alpha =1/\rho _S^\alpha `$ are the conductivities of the S’ layer and the respective bands of the S metal, $`D`$ is the diffusion constant in S’ and $`T_{cS}`$ is the critical temperature of S. The ratio between the parameters $`\gamma ^\alpha `$ for the different bands is mainly determined by the relation between the diffusion constants, $`D_S^\alpha `$. In the case of MgB2, the $`\pi `$-band is generally considered to be more dirty than the $`\sigma `$-band cleandirty , i.e. $`D_S^\pi D_S^\sigma `$.

The second boundary condition relates the gradient of the Green’s function $`\mathrm{\Phi }`$ near the S-S’ interface to its jump at the interface due to the finite interface resistance,

$$\xi G\mathrm{\Phi }^{}=\underset{\alpha }{}\frac{G_S^\alpha }{\gamma _B^\alpha }\left(\mathrm{\Phi }_S^\alpha \mathrm{\Phi }\right),$$
(5)
where $`\gamma _B^\alpha =R_B^\alpha /\rho \xi `$. $`R_B^\alpha `$ are the components of the specific interface resistance, describing the tunneling of an electron across the interface into the corresponding conduction band.

In order to obtain the resistances $`R_B^\alpha `$, we have to evaluate the effective junction transparency components. It was first pointed out by Mazin Mazin , that the normal state conductance $`R_\alpha ^1`$, in the limit of a specular barrier with small transparency, is proportional to the Fermi-surface average $`Nv^2_\alpha `$, where $`N`$ is the density of states and $`v`$ the Fermi velocity. In Ref. Brinkman,  it was further shown that the normal state resistance component of tunneling into band $`\alpha `$ of S is given by the contribution of the electrons in band $`\alpha `$ to the squared plasma frequency $`\left(\omega _p^\alpha \right)^2`$, which can be obtained from first principle calculations. For MgB2, the ratio $`R_B^\sigma /R_B^\pi =(\omega _p^\pi /\omega _p^\sigma )^2`$ is $`2`$ and $`100`$ for tunneling in the direction of the $`ab`$ plane and $`c`$-axis respectively Brinkman .

In the case of a SS’ bilayer, the Usadel equation (1) needs to be solved in the S as well as in the S’ layer, together with the self-consistent determination of the pair potentials in S and S’, Eq. (2). A general numerical method, using $`\mathrm{\Theta }`$-parametrization, $`\mathrm{\Phi }=\omega \mathrm{tan}\mathrm{\Theta }`$ and $`G=\mathrm{cos}\mathrm{\Theta }`$, is described for the single-band case in Ref. Golubov, . Here, we extended this method by applying the new boundary conditions, Eqs. (3) and (5). The density of states at energy $`E`$ can be obtained by applying an analytical continuation, $`\omega =iE`$, to the Usadel equations and the boundary conditions and solving the numerical scheme in the complex energy plane.

The numerically obtained dependence of the pair-potential on position is presented in Fig. 1 for the example in which the coupling constants are taken as calculated for MgB2 in Ref. const, . The parameter values are indicated in the caption. For temperatures above $`T_{cS^{}}`$ (solid lines in Fig. 1), it can be seen that the pair-potential in S’ increases towards the interface, while $`\mathrm{\Delta }_\sigma `$ decreases, as expected in analogy with the proximity effect in the single-band case. The decrease in $`\mathrm{\Delta }_\sigma `$ towards the interface can be explained by the relatively strong coupling between the $`\sigma `$ and $`\pi `$ bands. By decreasing the interband coupling constants and by decreasing the interface suppression parameters, one can obtain the opposite regime, in which $`\mathrm{\Delta }_\pi `$ increases towards the interface. For relatively large values of $`T_{cS^{}}`$, and for temperatures below $`T_{cS^{}}`$, we even predict an increase in $`\mathrm{\Delta }_\sigma `$ towards the interface, as illustrated by the dashed line in Fig. 1. This effect is further enhanced when the $`\sigma `$-band is decoupled from the S’-layer, which is the case for example when the interface normal is parallel to the crystallographic $`c`$-axis of MgB2, due to the vanishingly small ratio $`R_B^\pi /R_B^\sigma `$ in that case.

As an example, in Fig. 2 the results of a calculation of the density of states in a SN-bilayer are presented. In the considered case, the bulk energy gaps in a two-band superconductor are not too different. As is seen from the figure, the density of states in the N-layer has three peaks: the lowest energy peak corresponds to the proximity induced minigap and the two other peaks correspond to the bulk energy gaps in the two-band superconductor S. The existence of a minigap is a characteristic feature of the proximity effect in a SN-bilayer in the dirty limit, as was studied in detail in the single-band case in Ref. Golubov, . As we can see, the minigap persists in the two-band case as well and its magnitude depends on the parameters of the interface, thicknesses of the N and S and the values of the bulk gaps in the superconductor.

The next step in investigating the influence of multiband superconductivity on the proximity effect, is to study supercurrents in multiband proximized structures. We will consider double-barrier structures consisting of two S electrodes coupled by a normal metal N. As a model system we use a SINIS double-barrier hybrid structure, since in practical devices interface potential barriers are always present at the S-N interfaces, either originating from a Fermi-velocity mismatch, degradation of surface layers, or artificially deposited oxide barriers.

If the conditions of the dirty limit (electron mean free path $`\mathrm{}d,\xi `$) are fulfilled in the N interlayer, than the stationary Josephson effect in the structure can be analyzed in the framework of the Usadel equations by the method developed in Refs. KL,  and Zai1,  for the single-band case. We assume that the interface transparencies are small enough such that the condition $`1+\gamma _{B1,2}^\alpha \gamma _{1,2}^\alpha `$ holds at both NS interfaces (here and below we drop the subscript $`S`$). In this case, the suppression of superconductivity in the S layers is weak and the Green’s functions in the electrodes near the interfaces, $`G_{1,2}^\alpha `$ and $`\mathrm{\Phi }_{1,2}^\alpha `$, are equal to their bulk values. To calculate the supercurrent, it is sufficient to consider Eq. (5) at the two interfaces, giving

$$\xi G\mathrm{\Phi }^{}=\underset{\alpha }{}\frac{G_{1,2}^\alpha }{\gamma _{B1,2}^\alpha }\left(\pm \mathrm{\Phi }_{1,2}^\alpha \mathrm{\Phi }\right),x=\pm \frac{d}{2}.$$
(6)
For simplicity, we will consider symmetric junctions where $`G_{1,2}^\alpha G_S^\alpha `$ and $`\gamma _{B1,2}^\alpha \gamma _B^\alpha `$, and where the functions $`\mathrm{\Phi }_{1,2}^\alpha `$ are related to the phase shift $`\phi `$ across the junction by $`\mathrm{\Phi }_{1,2}^\alpha =\mathrm{\Delta }_\alpha \mathrm{exp}\left(\pm i\phi /2\right)`$. Further, we consider purely normal N-layer with vanishing pair potential $`\mathrm{\Delta }=0`$ and restrict ourselves to considering the limit of a small interlayer thickness, $`d\xi `$.

In the limit $`d\xi `$, there are two characteristic frequencies $`\mathrm{\Omega }_{1,2}`$ in the Usadel equations (1,2). At $`\omega \mathrm{\Omega }_1=\pi T_{cS}\frac{\xi }{d}\mathrm{\Omega }_2=\pi T_{cS}`$ we can neglect all nongradient terms in the Usadel equation. Hence, $`\left[G^2\mathrm{\Phi }^{}\right]^{}=0`$, and in the zero approximation on $`d/\xi `$ one obtaines that all $`\mathrm{\Phi }`$ functions are spatially independent constants, $`\mathrm{\Phi }=A`$. In the next approximation we have

$$\mathrm{\Phi }=A+B\frac{x}{\xi }+A\frac{x^2\beta ^2}{2\xi ^2},\beta ^2=\frac{\omega }{\pi T_{cS}G}.$$
(7)
From the boundary conditions and by taking into account that in our model $`\mathrm{\Phi }_{1,2}^\alpha =\mathrm{\Delta }_\alpha \mathrm{exp}\left(\pm i\phi /2\right)`$, we finally will have

$`A`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}\eta }{\stackrel{~}{G}}},B={\displaystyle \frac{i\stackrel{~}{\mathrm{\Delta }}}{G}}{\displaystyle \frac{d}{\xi }}\mathrm{sin}(\phi /2),`$ (8)
$`G`$ $`=`$ $`{\displaystyle \frac{\omega }{\sqrt{\omega ^2+A^2}}}={\displaystyle \frac{\omega \stackrel{~}{G}}{\sqrt{\omega ^2\stackrel{~}{G}^2+\stackrel{~}{\mathrm{\Delta }}^2\eta ^2}}},`$ (9)
where $`\gamma _{BM}^\alpha =\gamma _B^\alpha d/\xi `$, $`\eta ^2=\mathrm{cos}^2(\phi /2)`$, and

$$\stackrel{~}{G}=\underset{\alpha }{}\frac{G_S^\alpha }{\gamma _{BM}^\alpha }+\frac{\omega }{2\pi T_{cS}},\stackrel{~}{\mathrm{\Delta }}=\underset{\alpha }{}\frac{G_S^\alpha \mathrm{\Delta }_\alpha }{\gamma _{BM}^\alpha }.$$
(10)
The density of states $`N(E)=\mathrm{Re}(\mathrm{G})`$ in the interlayer of the double-barrier junction can now be found from an analytical continuation of Eq. (9) to real energies, $`\omega =iE`$. The results for the two-band case are plotted in Fig. 3. The known density of states for a single-band SINIS junction PRB  is shown in the inset. For $`\gamma _{BM}1`$, the single-band results show a peak in the density of states at $`\mathrm{\Delta }\mathrm{cos}(\varphi /2)`$, while the density of states in the two-band junction in this regime is predicted to have a peak at a value that is even lower than $`\mathrm{\Delta }_\pi \mathrm{cos}(\varphi /2)`$. For larger values of $`\gamma _{BM}`$, the density of states shows three peaks: at the minigap and at $`\mathrm{\Delta }_\pi `$ and $`\mathrm{\Delta }_\sigma `$, in analogy with the two peaks in the density of states of a single-band SINIS junction.

Substituting Eq. (8) into the supercurrent expression,

$$I=\sigma 2\pi T\mathrm{Im}\underset{\omega 0}{}\frac{1}{\omega ^2}\mathrm{G}^2\mathrm{\Phi }^{}\mathrm{\Phi }^{},$$
(11)
and taking $`\mathrm{\Phi }^{}\mathrm{\Phi }^{}`$ in lowest order equal to $`A^{}B`$, we obtain

$$I=\frac{\pi T}{\xi \rho }\underset{\omega 0}{}\frac{d\stackrel{~}{\mathrm{\Delta }}^2\mathrm{sin}(\phi )}{\omega \xi \sqrt{\omega ^2\stackrel{~}{G}^2+\stackrel{~}{\mathrm{\Delta }}^2\eta ^2}}.$$
(12)
A generalization to take boundary asymmetry and a finite $`\mathrm{\Delta }`$ in the interlayer into account can be made straightforwardly.

In the two-band case in the limit $`\gamma _{BM}^\sigma \mathrm{}`$, which is for example the case for tunneling in the MgB2 $`c`$-axis direction, the normal metal is only proximized by the $`\pi `$-gap of electrode S and Eq. (12) gives

$$I=\frac{\pi T}{\xi \rho \gamma _B^\pi }\underset{\omega 0}{}\frac{G_S^\pi \mathrm{\Delta }_\pi ^2\mathrm{sin}(\phi )}{\omega \sqrt{\left[\omega +\frac{\omega ^2\gamma _{BM}^\pi }{2\pi T_{cS}G_S^\pi }\right]^2+\mathrm{\Delta }_\pi ^2\mathrm{cos}^2(\frac{\phi }{2})}},$$
(13)
which has been previously obtained KL ; Kupr  for SINIS junctions with single-band superconductivity in S. It also follows from Eqs. (12) and (13) that the critical current of a single-band SINIS junction always exceeds the critical current of a two-band SINIS junction with a vanishing second gap, $`\mathrm{\Delta }_\pi =0`$.

The temperature dependence of the critical current can now be calculated for SINIS Josephson structures for different orientations of the crystallographical axis with respect to the interface normal. The gap functions, $`\mathrm{\Delta }_{\pi ,\sigma }(T)`$ and the ratio $`\gamma _B^\sigma /\gamma _B^\pi `$, follow from band structure calculations Brinkman . The results are shown in Fig. 4 for vanishingly small $`\gamma _{BM}`$, and compared to the calculation results for SIS junctions Brinkman . The full specific interface resistance of a SINIS junction, $`R_N=R_B^\sigma R_B^\pi /[R_B^\sigma +R_B^\pi ]`$. It is clearly seen that the critical current of SINIS junctions is larger than in SIS structures, practically in the whole temperature region, as is the case for single band superconductors KL ; Kupr . At low temperatures the $`I_cR_N`$ product can be as large as 5.2 mV when only the $`\pi `$-band contributes to the current and close to 7.3 mV when the sum over different band contributions can be taken into account, as is the case for tunneling in the direction of the $`ab`$ plane. The negative curvature of $`I_cR_N(T)`$ is a direct consequence of the two-band nature of superconductivity and is absent in $`I_cR_N(T)`$ of single-band SINIS junctions in the regime of small $`\gamma _{BM}`$ Kupr .

In summary, we have formulated a microscopic theory of the proximity effect in hybrid structures based on multiband superconductors in the diffusive limit. We have shown that the existence of multiple superconducting bands manifests itself in the proximity effect between a normal metal and a superconductor as the occurence of additional peaks in the density of states at the normal metal side. The interplay between the proximity effect and interband coupling determines the gap magnitudes at the interfaces. The supercurrent in multiband proximized Josephson junctions was calculated and compared to known single-band results and predictions for multiband tunnel junctions.

We acknowledge useful discussions with O.V. Dolgov and A.E. Koshelev and support from INTAS project 2001-0617. M.Yu.K. acknowledges support from the Russian Ministry of Industry, Science and Technology.


format: text

text:
# Velocity excitations and impulse responses of strings — Aspects of continuous and discrete models

## I Introduction

Our purpose here is to discuss aspects of the relationship of the solution of the one-dimensional second order wave equation to two discrete models thereof. The first discrete model is the digital waveguide model in one spatial dimension. The second discrete model is a finite difference model in the time domain. In particular we will also discuss how this relationship explains different behavior between the discrete models. This relationship has drawn much attention recentlyBBKS03 ; Bilbao01 ; Bilbao03 ; BS03 ; EK02 ; Karjalainen03 ; KE03 ; KS03 ; Smith03

It is shown that the finite difference model can account for solutions of the wave equation and that these solutions are physically meaningful.

This allows for a direct interpretation of recent results by Karjalainen and Erkut KE03  from the fundamental solution of the wave equation. Karjalainen and Erkut gave the restricting conditions necessary to make finite difference models in the time domain and digital waveguide models connect.

Regarding the stability behavior of the discrete models, the continuous stability is discussed and it is shown that physically stable responses of the string may appear unstable in a discrete model or signal-processing sense.

The paper is structured as follows. First known derivations of the solution of the wave equation is given, both classically and via the theory of fundamental solutions of its partial differential operators. A discrete comparison of the finite difference Ames92  and the digital waveguide model222For a friendly introductory exposition see CookCook02 . follows. Smith’s textSmith03  provides the authoritative summary in of digital waveguides with respect to the continuous derivations given earlier. We conclude with implications of these observations.

## II Solution of the Wave Equation in One Dimensions

The results in this section are well-known. They are repeated here to facilitate arguments in the following sections. Of concern is the general digital simulation of a string under force. The free ideal string is well-described by the the $`1+1`$ dimensional d’Alembertian operator on a scalar field333We limit our discussion to this particular form. A note regarding the relationship of this equation to pairs of first order equations can be found in Appendix B.:

$`\mathrm{}y(x,t)\stackrel{\text{def}}{=}({\displaystyle \frac{^2}{x^2}}c^2{\displaystyle \frac{^2}{t^2}})y(x,t)`$ (1)
An external force leads to the inhomogeneous case of the same equation:

$`\mathrm{}y(x,t)=f(x,t)`$ (2)
Without loss of generality we assume that $`c=1`$ for our discussion and hence one gets the factored form of the d’Alembertian:

$`\mathrm{}={\displaystyle \frac{^2}{x^2}}{\displaystyle \frac{^2}{t^2}}=\left({\displaystyle \frac{}{x}}+{\displaystyle \frac{}{t}}\right)\left({\displaystyle \frac{}{x}}{\displaystyle \frac{}{t}}\right)`$ (3)
The factored form suggests the substitution $`\xi =xt`$ and $`\eta =x+t`$ we obtain the canonical form of the wave equation:

$`{\displaystyle \frac{^2u}{\xi \eta }}={\displaystyle \frac{1}{4}}f({\displaystyle \frac{\xi +\eta }{2}},{\displaystyle \frac{\eta \xi }{2}})=\varphi (\xi ,\eta )`$ (4)
This can be directly integrated, yielding CH68

$`u={\displaystyle _{\eta _0}^\eta }{\displaystyle _{\xi _0}^\xi }\varphi (u,w)𝑑u𝑑w+h_1(\xi )+h_2(\eta )`$ (5)
where $`h_1`$, and $`h_2`$ are “constants of integration”. However, this notation hides that either one of these constants has been integrated over with regards to its parameter. Also these are not uniquely defined functions from a yet undefined functional space but rather any function from suitable family of functions CH68 . The derivation procedure suggests that $`h_1()`$ and $`h_2()`$ are one and twice differentiable everywhere in the solution space, but which one is twice differentiable depends in what order the solution has been integrated over.

The solution of the homogeneous case (1) with initial conditions $`y(x,0)=f(x)`$ and $`y_t(x,0)=g(x)`$ is well known Kreyszig99 ; SS03  to correspond to d’Alembert’s solution444It is noteworthy to mention that the nomenclature for “d’Alembert’s solution” is not universally agreed upon. Some sources use it to describe the the general form of the integration of the homogeneous equation in the absence of initial value or force (also called Cauchy) data (like Kreyszig99 ) whereas in other sources, for exampleEKS99 , it refers to the solution including Cauchy (initial, boundary, external force) data. We will use it to denote the solution of the initial value problem.:

$`y(x,t)={\displaystyle \frac{1}{2}}\left(f(x+ct)+f(xct)\right)+{\displaystyle \frac{1}{2c}}{\displaystyle _{xct}^{x+ct}}`$ $`g(s)ds`$ (6)
Hence initial displacements travel left and right. Initial velocity, however smears over a widening domain of influence.

Alternatively this results can be derived using the theory of partial differential operators. Writing $`L=\mathrm{}`$ and $`X=(x,t)`$ one arrives at the generic partial differential operator equation:

$`Lu(X)=k(X)`$ (7)
Here it is important to note that it is no longer required that $`u()`$ is in the class of twice differentiable functions $`𝒞^2`$, but rather that $`u()`$ and $`k()`$ are distributions555Laurent Schwartz, the originator of the theory of distributions recently wrote an autobiography Schwartz01  that includes a readable account of the historical development of his discovery. We quote two statements from his book that are of interest relating to this article. He writes, when describing the prehistory of his discovery (Schwartz01, , p. 212): “One of the most important precursors of distributions was the electrical engineer Heaviside.” In a section headed Vibrating strings, harmonic functions he writes (Schwartz01, , p. 218): “We had learned that the general solution of this equation is of the form $`u(t,x)=f(x+vt)+g(xvt)`$, where $`f`$ and $`g`$ are arbitrary functions of one variable. Naturally, this presupposes that $`f`$ and $`g`$ are $`C^2`$, so as to be able to differentiate them. What should one think of a function $`u`$ which would be analogous except that $`f`$ and $`g`$ would be merely $`C^1`$ (continuous), or not even continuous? Is it a wave or not? I was obsessed by this question for some time, then I stopped thinking about it and relegated the question to a corner of my mind for future reflection.” or generalized functions. What this means in detail we will have to defer to expositions available elsewhere Edwards95 ; Hulshof93 ; EKS99 . For our purpose it is interest to note that in the theory of distributions jumps and discontinuity are gracefully and meaningfully included in the formalism including the definition derivatives of entities like the Dirac-delta $`\delta ()`$ over the space of suitable functions testing for this property. In addition it makes continuous convolutions a central operation to the calculation of continuous solution, which make it’s treatment very similar to the study of discrete models OS89  which is regularly used for digital waveguide modelsSmith03 .

In our treatment here we will closely follow Hulshof’sHulshof93  and Joshi and Wassermann’s lecture notesJW03  and also Edwards’ textEdwards95  which are more accessible than, for example the technical survey by Egorov, Komech and ShubinEKS99 . The interested reader is pointed to the latter for statements of necessary theorems as well as proofs or detailed references to original proofs.

The fundamental solution of the equation (7) the effect of the operator $`L`$ on the distribution $`u`$ when it sees as input a Dirac-delta (as noted earlier this is not a function, but a distribution). In the digital signal processing literatureOS89  $`u`$ is called the impulse response of $`L`$. In a dynamical sense it is the response to the inhomogeneous equation where the external force distribution is a Dirac-delta $`\delta `$. Of immediate interest are solution forward in time the discussion is restricted to fundamental solutions in the positive half-plane with respect to time. This will be indicated by the superscript $`+`$ to the symbol $``$ for the fundamental solution:

$`L^+(X)=\delta (X)`$ (8)
It can be proved that the solution of both the homogeneous equation with initial value data and the inhomogeneous equation with some external force distribution can be recovered from the fundamental solution EKS99 . This is done, in analogy to the impulse response convolution in digital signal processing OS89  by convolution of the fundamental solution with the force distribution and the initial value data $`u(x,0)=f()`$ and $`u_t(x,0)=g()`$. In the literature these are call Cauchy data.

The fundamental solutionHulshof93 ; Edwards95  of the one-dimensional wave equation can be derived to be:

$`^+(x,t)={\displaystyle \frac{1}{2}}H(t)\left[H(x+t)H(xt)\right]`$ (9)
Here $`H()`$ is the Heaviside distribution, which is the distributional integral of the Dirac-delta distribution $`\delta ()`$. In conventional functional form the Heaviside step-“function” can be written asAS72 :

$`H(x)=\{\begin{array}{cc}0\hfill & x<0,\hfill \\ \frac{1}{2}\hfill & x=0,\hfill \\ 1\hfill & x>0.\hfill \end{array}`$ (10)
The interpretation of this equation is indeed important because it indicates, that the “system response” of a wave operator to an input impulse are not isolated traveling wave pulses but traveling step-distributions.

It may be convenient to think of the fundamental solution as the “distributional continuous impulse response”. This makes sense because the solution of equation can be recovered by convolution of the fundamental solution with the Cauchy data. The continuous convolution has the familiar form Hulshof93 ; EKS99 :

$`u(x)=(^+f)(x)\stackrel{\text{def}}{=}{\displaystyle _{\mathrm{}}^{\mathrm{}}}^+(xs)f(s)𝑑s`$ (11)
where $`f(x)`$ is Cauchy data in one variable.

If both data and the fundamental solution are in two dimensions (as is possibly the case for external force distributions), one need to convolve over both variables.

Specifically it can be shown (see EKS99 ; Hulshof93 ) that for the set of Cauchy data $`u(x,0)=f(x)`$, $`u_t(x,0)=g(x)`$ and $`Lu(x,t)=k(x,t)`$ one gets the complete solution:

$`u(x,t)=u_f(x,t)+u_g(x,t)+u_k(x,t)`$ (12)
with the convolutions:

$`u_f(x,t)`$ $`=_t^+(,t)f(x)`$ (13)
$`u_g(x,t)`$ $`=^+(,t)g(x)`$ (14)
$`u_k(x,t)`$ $`=^+k(x,t)`$ (15)
Performing the convolutions yields the the solution which is equivalent to the conventional inhomogeneous initial value solution of the wave equation Hulshof93 :

$`\begin{array}{cc}\hfill y(x,t)=& {\displaystyle \frac{1}{2}}\left(f(x+t,0)+f(xt,0)\right)\hfill \\ & +{\displaystyle \frac{1}{2}}{\displaystyle _{xt}^{x+t}}g(s,0)𝑑s\hfill \\ & +{\displaystyle \frac{1}{2}}{\displaystyle _0^t}{\displaystyle _{x(t\tau )}^{x+(t\tau )}}k(s,\tau )𝑑s𝑑\tau \hfill \end{array}`$ (16)
For proofs of uniqueness see Egorov, Komech and ShubinEKS99 . In the absence of external forces this reduces to the initial value solution (6).

The theory of generalized functions for partial differential operators explains why equations (5), derived for the forced case, and (6), derived for the homogeneous initial value case have similar structure. The inhomogeneous case can in a generalized sense be made to include the homogeneous initial value problem (also called Cauchy problem). For example, if we symbolically write $`\stackrel{~}{k}(x,t)=\delta (t)k(x,t)`$ one sees that the external force response matches the response to the initial velocity.

$`u_g(x,t)={\displaystyle \frac{1}{2}}{\displaystyle _{xt}^{x+t}}g(s,0)𝑑s`$ (17)
$`u_{\stackrel{~}{k}}(x,t)={\displaystyle \frac{1}{2}}{\displaystyle _{xt}^{x+t}}\stackrel{~}{k}(s,0)𝑑s`$ (18)
Hence an external force distribution which is impulsive in time $`\stackrel{~}{k}(x,)`$ is indistinguishable from the equivalent initial velocity distribution $`g(x)`$. Conversely it is noteworthy that initial values do not prescribe a state of the string alone, but also prescribe a sudden onset of such state. Hence initial values are not necessarily a free-field solution of the wave equation, that is the state of a string in the absence of force. Rather impulsive onset states act like external forces.

In particular one can write the solution of the wave equation for the inhomogeneous initial value problem in the following simple form JW03 :

$`u(x,t)=^+g(x)+_t^+f(x)+{\displaystyle _0^t}_{ts}^+k(x,s)𝑑s`$ (19)
Next we note the differentiation of distributions on step-functions (see also the Appendix A):

$`H^{},\varphi =\delta ,\varphi `$ (20)
observe, using this relationship, that the derivative of the fundamental solution (9) are two propagating Dirac-deltas, here again in symbolic notation:

$`_t^+={\displaystyle \frac{1}{2}}H(t)\left(\delta (x+t)+\delta (xt)\right)`$ (21)
For the current discussion it remains only to point out that the first contribution of (19) look like propagating impulses under differentiation.

From (19) we see that the solution of the wave equation will only stay on the characteristic lines $`\xi `$ and $`\eta `$ for a restricted class solutions of the wave equation. Only for very exceptional cases of initial velocities $`g`$ and external forces $`k`$ will the solution not integrate into the domain. Rather generically one ought to expect them to integrate into the inside of the forward characteristic cone $`x\pm t0`$ as depicted in Figure 1.

The condition in which integration into the interior of the characteristic cone does not occur will be discussed in the discrete case and in this context has been discovered by Karjalainen and Erkut KE03 . For this to be appropriate, care needs to be taken when assuming the impulse response of the system to have a particular form, otherwise one ought to expect both contributions at the same time. Additionally, in a physically realistic situation, when it cannot be guaranteed that an excitation is of purely displacement-type, one ought to expect that the dynamics of the system to integrate and persist over the whole inside of the characteristic cone. It also is worthwhile to point out that this is the mathematically consistent solution of the wave equation EKS99 . Hence, assuming that the wave equation is a reasonable model for a given physical situation, one ought to expect such a behavior to exist and be observable.

As a note, we state, that it is in fact well known that the Huygens’ Principle, the isolated propagation of wave fronts, only holds for d’Alembertians of odd spatial dimensions greater or equal to $`3`$ EKS99 ; CH68  meaning that only in this case does the fundamental solution concentrate on the characteristic cone $`t^2=|x|^2`$. Other cases, including one spatial dimension666The case of two spatial dimensions is relevant, but is not in the main thrust of this papers discussion. An illustration of the one and two-dimensional wake of the wave-equation can be found in GraffGraff91 , p. 220, Fig 4.4., have a wave influence inside the characteristic cone. The idea that the solution would concentrate on the characteristic cone for all dimensions was held for a long time by mathematicians working on the wave-equation until Hadamard opened the development and the situation is since well understood Garding98 .

Regarding Huygens’ principle, a short popular exposition can be found in Veselov02  whereas a comprehensive technical exposition can be found in Zubelli97 .

## III Comparison of Leapfrog and Waveguide Solvers

Next the leapfrog finite difference scheme Ames92  will be compared with the digital waveguide method Smith03 . The full treatment of derivations will not be repeated here and the reader is referred to these sources for details.

Now let us discuss the properties of the so-called leap-frog finite difference molecule for the wave equation. The explicit time-stepping equation reads Ames92 ; Smith03 :

$`\begin{array}{cc}\hfill y(n+1,m)=& y(n,m1)+y(n,m+1)\hfill \\ & y(n1,m)\hfill \end{array}`$ (22)
where the relationship between the discrete time step $`T`$ and the spatial discretization $`X`$ is chosen to satisfy $`c=X/T`$. In this case the leapfrog molecule can be shown to be consistent at sampling points with the wave equation Ames92 . It can also be shown that waveguide solutions are solutions of the leapfrog Smith03 . As Karjalainen observes, the converse does not hold Karjalainen03 .

For future discussion we will use the following symbols:

$`y_+`$ $`=y(n+1,m)`$ (23)
$`y_>`$ $`=y(n,m+1)`$ (24)
$`y_<`$ $`=y(n,m1)`$ (25)
$`y_{}`$ $`=y(n1,m)`$ (26)
The condition that an impulse at the root of the molecule will only create responses along the characteristics of the wave can be expressed by the condition $`y_+=0`$, i.e. there is no data within the characteristic domain of the molecule.

Hence the non-integrating molecule condition reads:

$`y_<+y_>y_{}=0`$ (27)
From this we get the relationship of waves on the characteristics to their sum:

$`y_<+y_>=y_{}`$ (28)
The updating rules for waveguides777Strictly speaking digital waveguide synthesis can be formulated in various ways. We will not use any arguments from transmission-line theory here. For treatment of those aspects we refer to Smith03 . are:

$`y_l(n,m)=y_l(n1,m+1)`$ (29)
$`y_r(n,m)=y_r(n1,m1)`$ (30)
with the external force rule:

$`y_l(n,m)={\displaystyle \frac{1}{2}}f(n,m)`$ (31)
$`y_r(n,m)={\displaystyle \frac{1}{2}}f(n,m)`$ (32)
The wave reconstruction rule is:

$`y(n,m)=y_l(n,m)+y_r(n,m)`$ (33)
in response to an external force function $`f(n,m)`$. If we take an impulse of height $`y_{}`$ at time $`n1`$ we get:

$`y_l(n,m1)=y_l(n1,m)=f(n1,m)={\displaystyle \frac{1}{2}}y_{}`$ (34)
$`y_r(n,m+1)=y_r(n1,m)=f(n1,m)={\displaystyle \frac{1}{2}}y_{}`$ (35)
The reconstructed wave using (33) is zero everywhere except at:

$`y(n,m1)=y_l(n,m1)={\displaystyle \frac{1}{2}}y_{}`$ (36)
$`y(n,m+1)=y_r(n,m+1)={\displaystyle \frac{1}{2}}y_{}`$ (37)
$`y(n,m)=y_l(n1,n)+y_r(n1,m)=y_{}`$ (38)
and we see that the non-integrating case of the leapfrog (28) is satisfied with $`y_<=1/2y_{}`$ and $`y_>=1/2y_{}`$.

The leapfrog will “integrate” whenever condition (28) is not satisfied. To study the behavior within the characteristic domain it is first assumed that the elements on the characteristic of the molecule $`y_<`$ and $`y_>`$ are unaltered. That is, the same waves as before travel outward in the molecule. This leaves us to study an altered relationship between $`y_+`$ and $`y_{}`$.

Let $`y_{}`$ be the difference of $`y_{}^0`$, the molecule value for the non-integrating case (28), and $`\stackrel{~}{y}_{}`$, an assumed contribution to the interior of the characteristic domain. Then we get:

$`y_+`$ $`=y_<+y_>(y_{}^0+\stackrel{~}{y}_{})`$ (39)
$`0`$ $`=y_<+y_>y_{}^0`$ (40)
Subtracting (40) from (39) we get:

$`y_+`$ $`=\stackrel{~}{y}_{}`$ (41)
Hence the response at at the interior point of the characteristic domain is constant with regards to the contribution of the incoming wave that violates the non-integration condition (28).

To study how this behavior, one can illustrate the response of the leapfrog to an initial $`1`$:

$$\begin{array}{ccccccccc}1& & 1& & 1& & 1& & 1\\ & 1& & 1& & 1& & 1& \\ & & 1& & 1& & 1& & \\ & & & 1& & 1& & & \\ & & & & 1& & & & \end{array}$$
(42)
and compare it to an excitation which observes (28). $`y_<=y_>=1`$ and $`y_{}=2`$:

$$\begin{array}{ccccccccc}1& & & & & & & & 1\\ & 1& & & & & & 1& \\ & & 1& & & & 1& & \\ & & & 1& & 1& & & \\ & & & & 2& & & & \end{array}$$
(43)
Observe that (43) appears to be the sum of traveling histories and they are time-symmetric around the intersection point. A time-symmetric solution is an equal contribution to the solution traveling forward and backward in time and their sum yielding the complete solution.

D’Alembert’s solution (6) can be used to investigate this observation when writing it in the form following Alpert, Greengard and HagstromAGH00 :

$`\begin{array}{cc}\hfill y(x,t)+y(x,t)=& {\displaystyle \frac{1}{2}}\left(f(x+ct)+f(xct)\right)\hfill \\ & +{\displaystyle \frac{1}{2c}}{\displaystyle _{xct}^{x+ct}}g(s)𝑑s\hfill \\ & +{\displaystyle \frac{1}{2}}\left(f(xct)+f(x+ct)\right)\hfill \\ & +{\displaystyle \frac{1}{2c}}{\displaystyle _{x+ct}^{xct}}g(s)𝑑s\hfill \end{array}`$ (44)
Taking the time-symmetric sum we get:

$`y(x,t)+y(x,t)=`$ $`f(x+ct)+f(xct)`$ (45)
Similarly, by taking the difference, one finds the time-asymmetric case:

$`y(x,t)y(x,t)={\displaystyle \frac{1}{c}}{\displaystyle _{x+ct}^{xct}}`$ $`g(s)ds`$ (46)
In the discrete case it is easy to see this property preserved in the leapfrog case:

$$\begin{array}{ccccccccc}& & 1& & 1& & 1& & \\ & & & 1& & 1& & & \\ & & & & 1& & & & \\ & & & 0& & 0& & & \\ & & & & 1& & & & \\ & & & 1& & 1& & & \\ & & 1& & 1& & 1& & \end{array}$$
(47)
$$\begin{array}{ccccccccc}& & 1& & & & 1& & \\ & & & 1& & 1& & & \\ & & & & 2& & & & \\ & & & 1& & 1& & & \\ & & 1& & & & 1& & \end{array}$$
(48)
While (48) nicely illustrates the time-symmetry and the “interference” of waves at the interaction point, (47) is insightful as it clearly shows the properties of a velocity excitation. The displacement vanishes at the interaction moment, while the temporal slope is maximal. It should be made clear, that vanishing of data at one time-step in the leap-frog simulation does not imply velocity solutions. This can be seen if a positive and a negative impulsive wave cross, creating a time-symmetric situation that will not integrate inside the domain:

$$\begin{array}{ccccccccc}& & 1& & & & 1& & \\ & & & 1& & 1& & & \\ & & & & 0& & & & \\ & & & 1& & 1& & & \\ & & 1& & & & 1& & \end{array}$$
(49)
We observe that this situation does satisfy the time-symmetric equation (45).

Alternatively similar results can be derived by discretizing the initial velocity $`g(x)`$ directly using a matching center difference scheme Kreyszig99 :

$`y^+y^{}=g`$ (50)
and including an arbitrary background field one gets:

$`y^+={\displaystyle \frac{1}{2}}(f^<+f^>)+g`$ (51)
where we use the notation $`f^+`$ and $`f^{}`$ to denote initial displacement wave contributions aligned with the left-right branch of the leapfrog-molecule. Observe that (51) does satisfy the same integration (41) and non-integration (28) conditions. Hence a velocity contribution $`g`$ can be interpreted as any violation of the rule of the sum of incoming traveling waves.

This behavior has been observed earlier. Karjalainen observed that an asymmetric pair of impulses need to be fed into a leapfrog motivated junction formulations to avoid integration behavior Karjalainen03b . The subsequent physical interpretation is derived in KE03  from a center-difference time-discrete velocity excitation. An interpretation of this result follows next.

## IV Singular Propagation from Integration

The non-integrating condition can be algorithmically enforced by a method discovered by Karjalainen and Erkut KE03 ; Karjalainen03 . Hence we will call this the Karjalainen-Erkut condition. The rule is to present the excitation through a feed-forward filter Karjalainen03 :

$`H(z)=1z^2`$ (52)
which can be derived from physical conditions by using a center difference velocity term KE03 .

We observe that the Karjalainen-Erkut condition (52) creates two impulses from one and those impulses are center symmetric and sign-inverted. If we calculate those two impulse responses separately and then create the sum, we see that the impulsive propagating solution comes about as the difference of two Heaviside distributions. The pulses are represented by their sign only as the amplitudes are assumed to be matched:

$$\begin{array}{ccccccccc}+& & +& & +& & +& & +\\ & +& & +& & +& & +& \\ & & +& & +& & +& & \\ & & & +& & +& & & \\ & & & & +& & & & \end{array}$$
(53)
$$\begin{array}{ccccccccc}& & & & +& & & & \end{array}$$

$$\begin{array}{ccccccccc}0& & & & & & & & 0\\ & 0& & & & & & 0& \\ & & 0& & & & 0& & \\ & & & 0& & 0& & & \\ & & & & 0& & & & \end{array}$$
(54)
$$\begin{array}{ccccccccc}& & & & =& & & & \end{array}$$

$$\begin{array}{ccccccccc}+& & 0& & 0& & 0& & +\\ & +& & 0& & 0& & +& \\ & & +& & 0& & +& & \\ & & & +& & +& & & \\ & & & & +& & & & \end{array}$$
(55)
Hence we see that in an impulse-response interpretation of the leap-frog, a Heaviside integration over the characteristic cone of influence is sensible and the Karjalainen-Erkut condition ensures that each Heaviside integration is matched with a delayed sign-inverted response that cancels all interior integration of the first impulse to leave unaltered the traveling impulse solution.

## V Effects of the Boundary

Next the effect of imposing spatial boundary conditions is studied. For this it is assumed that the solution of the wave equation is only meaningful and defined for a compact domain $`\mathrm{\Omega }`$. The length of the domain is denoted by $`L=|\mathrm{\Omega }|`$. For the boundary of the domain we write $`\mathrm{\Omega }`$ and the interior of the domain is defined by the quotient $`\mathrm{\Omega }\mathrm{\Omega }`$. On each distinct point of the boundary $`\mathrm{\Omega }`$ we impose one boundary condition. Fixed ends $`u(\mathrm{\Omega })=0`$ we call Dirichlet boundary conditions, whereas open ends $`u_t(\mathrm{\Omega })=0`$ we call Neumann boundary conditions. Note that a circular domain $`u(|\mathrm{\Omega }|)=u(0)`$ is a periodic unbounded domain.

The behavior at the boundary can be conveniently studied by extension of the domain. If the boundary is of Dirichlet type, the value of $`u`$ needs to vanish at the boundary and hence the extension needs to be odd. In the case of Neumann conditions $`u`$ needs to be even. As the resulting infinite extension is periodic in $`2L`$ this extension can be interpreted as a periodic unbounded domain of this length SS03 . We will denote the original domain by subscript $`0`$ and extended domains by indicies $`n0`$. The periodicity implies $`\mathrm{\Omega }_m+2=\mathrm{\Omega }_m`$ for all $`m`$.

The following discussion will be restricted to the behavior in response to the velocity term $`g`$ in (6). Observe that with periodicity we can write the integral as the sum of contributions of the periodic domains:

$`{\displaystyle _{\xi =xt}^{\eta =x+t}}g(s)𝑑s={\displaystyle \underset{m}{}}{\displaystyle _{\mathrm{\Omega }_m\xi ,\eta }}g(s)𝑑s`$ (56)
Hence we integrate over all contributions above the characteristic lines from an excitation point of the periodic domain.

With Dirichlet conditions one gets the odd extension Graff91 :

$`g(x)=g(2m|\mathrm{\Omega }|x)`$ (57)
and for Neumann conditions we get the even extension:

$`g(x)=g(2m|\mathrm{\Omega }|x)`$ (58)
with $`x\mathrm{\Omega }_0`$ and $`m0`$.

Integrating up to the point where the characteristic lines have reached $`2|\mathrm{\Omega }|`$ one gets for Dirichlet boundary conditions:

$`u(\mathrm{\Omega })={\displaystyle _{\mathrm{\Omega }_0}}g(s)𝑑s{\displaystyle _{\frac{1}{2}\mathrm{\Omega }_{\pm 1}}}g(s)𝑑s=0`$ (59)
Hence integral contributions cancel every $`2|\mathrm{\Omega }|`$ and the maximum amplitude is bounded by the integral of $`g()`$ over the original domain $`\mathrm{\Omega }_0`$.

The same procedure for Neumann boundary conditions leads to:

$`\begin{array}{cc}\hfill u(\mathrm{\Omega })& ={\displaystyle _{\mathrm{\Omega }_0}}g(s)𝑑s+{\displaystyle _{\frac{1}{2}\mathrm{\Omega }_{\pm 1}}}g(s)𝑑s\hfill \\ & =2{\displaystyle _{\mathrm{\Omega }_0}}g(s)𝑑s\hfill \end{array}`$ (60)
Figure 3 shows the behavior for a string tied at the ends (Dirichlet conditions) after an initial impulsive distribution. It reveals many properties of the effect of the boundary on the integration of velocities under Dirichlet boundary conditions. It shows the odd-periodic extension of the domain $`\mathrm{\Omega }_0`$ to $`\mathrm{\Omega }_n`$, $`n`$, it also shows the cancellation and constructive interference effect of overlapping integration regions. It also shows the $`2|\mathrm{\Omega }|`$ cancellation of waves. Erkut and KarjalainenEK02b  (compare their Figure 7) reported numerical simulations using the leapfrog molecule with comparable results, which hence matches the situation of the continuous model.

### V.1 Linear Growth of Displacement

Observe that the Neumann condition leads to a linear increase in the displacement as a response to velocity or force data being present.

The difference between the Dirichlet condition and the Neumann condition can be interpreted as the difference between an alternating sum and an accumulative sum.

In the Dirichlet case the sign of the area integrated over alternates with periodicity $`|\mathrm{\Omega }|`$ and hence any finite bounded signal $`g()`$ will produce an alternating sum which is bounded similarly but infinitely periodic.

In the Neumann the signs match and hence the area integrated over increases with every iteration over the domain by the integral over the finite bounded signal $`g()`$. Once the support of $`g()`$ has been exhausted, this obviously corresponds to a linear increase.

This is however, not an unphysical situation. This corresponds to constant kinetic energy being present in the string and hence implies that the energy is bound. This can easily be understood as linear increasing displacement implies constant velocity, which in turn implies constant energy. Hence energy is conservedSS03 . It can be interpreted as a string moving at constant velocity, which is meaningful as the Neumann conditions imply that the string is not tied down at the boundaries. Hence linear buildup in a displacement-like wave variable is energy-conserving888Evidently this argument is valid for any amplitude, also small amplitude oscillations, for which the wave-equation is valid, as constant displacement does not alter curvature..

Numerically this is still an undesirable situation because even in the absence of numerical imprecision, the dynamic range of numbers are bounded and hence an infinite increase cannot be represented.

The case of Neumann boundary conditions is interesting because it highlights the difference between notions of stability as customary in discrete signal literatureOS89  and stability in physical situations. The Neumann displacement response to any bounded input will be unbounded and hence is evidently not bounded-input bounded-output (BIBO) stable, see Oppenheim and SchaferOS89 , p. 20. We suggest that this BIBO-unstable but energy-conserving system be called physically stable. The BIBO-instability is a discrete computational problem and not one of the physical situation999Note that the velocity response is in fact BIBO-stable and hence treatment of the problem in a velocity variable will not suffer this problem..

## VI Implications

This paper discussed the linear lossless wave equation and its relationship to discrete models, namely a finite difference scheme called leapfrog, and the digital waveguide method. It is shown that the waveguide model corresponds to the finite difference scheme in the absence of integration. In the continuous case, integration can be expected to occur when initial velocities or external forces are present. In this light the observed integrating behavior of finite difference discretization in the time domain using the leapfrog molecule displays results consistent with the continuous model. Here we assume that the wave equation is at least in principle physically meaningful for the modeled situation. If this is the case one should expect consistent behavior of the related discrete models.

In relation to this argument, a use of waveguide discretization that does not include contributions inside the characteristic cone, does not include the integrating behavior of the model equation. In general both integrating and non-integrating responses are to be expected and hence should be present unless they can be explicitly excluded for physical reasons.

This also implies that the impulse response in just one variable in general will not carry the full dynamics of the string. Hence any assumption of the non-integrating impulse response in one variable in the construction of physical models might contain deviations for it only covers a reduced set of the solution space.

###### Acknowledgements.

Much thanks to Sile O’Modhrain for her support and input. The author also has much gratitude to send to Matti Karjalainen and Cumhur Erkut for stimulating discussions relating to this topic. I am also grateful for their sending of reprints and the graceful sharing of novel unpublished manuscripts. This work was made possible by the kindness of employment of Media Lab Europe and access to its academic resources.

## Appendix A Properties of Distributions

Let $`f`$ be a distribution on a real open interval $`\mathrm{\Omega }`$ and let $`\varphi `$ be in the the set of test functions $`𝒟(\mathrm{\Omega })`$ then one has Edwards95 :

$`\varphi ,f^{}={\displaystyle _\mathrm{\Omega }}f\varphi ^{}𝑑\mu =\varphi ^{},f`$ (61)
and for arbitrary derivatives:

$`\varphi ,^pf=(1)^{|p|}^p\varphi ,f`$ (62)
the Dirac delta $`\delta `$ has the property:

$`\varphi ,\delta =\varphi (0)`$ (63)
hence returns the value of $`\varphi `$ at $`0`$. By the differentiation rule the higher order derivatives of the Dirac delta returns the higher order derivatives at $`0`$ with alternating sign:

$`\varphi ,^p\delta =(1)^{|p|}^p\varphi (0)`$ (64)
Let $`H`$ be the Heaviside distribution. It is defined as Hulshof93 :

$`\varphi ,H={\displaystyle _{\mathrm{}}^{\mathrm{}}}H(s)\varphi (s)𝑑s={\displaystyle _0^{\mathrm{}}}\varphi (s)𝑑s`$ (65)
It hence permits the positive part of $`\varphi `$ over the domain. The derivative of the Heaviside distribution $`H`$ yields (using (61) and (64)) the Dirac-delta:

$`\varphi ,H^{}=\varphi ^{},H=\varphi (0)=\varphi ,\delta `$ (66)
## Appendix B The Wave Equation and First Order Systems

In order to derive the relationship between the wave equation to first order systems, we discuss two forms of such systems, namely, two transport equations in one variable and two transport equations in a mixed pair of variables.

A generic version of a system of inhomogeneous first order hyperbolic equations reads:

$`a{\displaystyle \frac{y}{x}}+b{\displaystyle \frac{y}{t}}=h_1(x)+h_2(t)`$ (67)
$`c{\displaystyle \frac{y}{x}}+d{\displaystyle \frac{y}{t}}=h_3(x)+h_4(t)`$ (68)
For simplicity assume that the force terms are separated in the independent dimensions. Then a second order version is usually derived taking the derivative of one equation with respect to $`t`$ and the other one with respect to $`x`$. The cross-term $`y_{xt}`$ can be eliminated and one gets two equations:

$`{\displaystyle \frac{b}{d}}{\displaystyle \frac{^2y}{t^2}}{\displaystyle \frac{c}{a}}{\displaystyle \frac{^2y}{x^2}}={\displaystyle \frac{1}{d}}{\displaystyle \frac{}{t}}h_2(t){\displaystyle \frac{1}{a}}{\displaystyle \frac{}{x}}h_3(x)`$ (69)
$`{\displaystyle \frac{a}{c}}{\displaystyle \frac{^2y}{t^2}}+{\displaystyle \frac{d}{b}}{\displaystyle \frac{^2y}{x^2}}={\displaystyle \frac{1}{c}}{\displaystyle \frac{}{x}}h_1(x){\displaystyle \frac{1}{b}}{\displaystyle \frac{}{t}}h_4(t)`$ (70)
The key observation is that one second order equation (69) or (70) is not strictly equal to the system of first order equations (67) and (68). It is only equal up to two functions (whichever got eliminated, $`h_1,h_4`$ or $`h_2,h_3`$). They are equivalent up to two ”constants of integration”.

For systems of first order linear equations in two independent variables a related proof holds. An intuitive interpretation is that in fact for first order equations of the type:

$`u_x+w_t`$ $`=g_1(t)`$ (71)
$`w_x+u_t`$ $`=g_2(x)`$ (72)
one sees that the reduction to second order equations in $`u`$ by differentiating (71) with respect to $`x`$ and (72) with respect to $`t`$ one gets:

$`u_{xx}u_{tt}`$ $`=0`$ (73)
$`w_{xx}+w_{tt}`$ $`=\dot{g_1}(t)g_2^{}(x)`$ (74)
Note that differentiation eliminated $`g_1`$ and $`g_2`$ in one case and hence the homogeneous wave equation is again indistinguishable for both the homogeneous and a class of inhomogeneous systems of first order equations and in this sense they are equivalent only up to two functions.