A Comprehensive Mathematical Model of Surface Electromyography and Force Generation

The purpose of this article is to provide a unified description of a comprehensive mathematical model of surface electromyographic (EMG) measurements and the corresponding force signal in skeletal muscles. The model comprises motor unit pool organization, recruitment and rate coding, intracellular action potential generation and the resulting EMG measurements, as well as the generated muscular force during voluntary isometric contractions. It consolidates and extends the results of several previous publications that proposed mathematical models for the individual model components. A parameterization of the electrical and mechanical components of the model is proposed that ensures a physiologically meaningful EMG-force relation in the simulated signals. Moreover, a novel nonlinear transformation of the excitation model input is proposed, which ensures that the model force output equals the desired target force. Finally, an alternative analytical formulation of the EMG model is proposed, which renders the physiological meaning of the model more clear and facilitates a mathematical proof that muscle fibers in this model at no point in time represent a net current source or sink. Neuromuscular physiology is a vibrant research field that has recently seen exciting advances. Many previous publications have focused on thorough analyses of particular aspects of neuromuscular physiology, yet an integration of the various novel findings into a single, comprehensive model is missing. A consistent description of a complete physiological model as presented here, including thorough justification of model component choices, will facilitate the use of these advanced models in future research. Results of a numerical simulation highlight the model’s capability to reproduce many physiological effects observed in experimental measurements, and to produce realistic synthetic data that are useful for the validation of signal processing algorithms. The model is based on recent advances in the understanding of muscular physiology and hence also applicable for analyzing the influence of various physiological and measurement setup parameters on the measured force and EMG signals.

The purpose of this article is to provide a unified description of a comprehen-

Introduction
Electromyography (EMG) denotes the measurement of the electrical fields generated by the electrophysiological processes that lead to muscle fiber contraction. EMG is highly relevant for a number of clinical and scientific application fields, since it enables monitoring and analysis of a muscle's electromechanical properties and state, both of which would otherwise remain mostly inaccessible. Surface electromyography (sEMG) denotes the noninvasive measurement of electrical muscle activity by means of electrodes placed on the skin surface, as opposed to the traditional measuring method using needle electrodes. For more background information on sEMG, its analysis and many of its applications, refer to, e.g., Merletti and Parker (30); Merletti and Farina (31).
Mathematical models of sEMG are highly useful, on the one hand to advance understanding of the underlying physiological processes, and on the other hand to analyze the sensitivity of sEMG measurements to various physiological and technical parameters, and to test and validate sEMG signal processing algorithms. Over the past decades, researchers have pursued a number of different approaches for the modelling and simulation of different aspects of sEMG measurements. Phenomenological (23; 28; 29; 50) as well as physiologically motivated (13-15; 17; 19; 20; 34) models have been proposed and analyzed. Overviews can be found in McGill (29); Rodriguez-Falces et al. (47); Stegeman et al. (51) and Merletti and Farina (31). Particular emphasis has been placed on modelling the electric signal produced by a single contraction of a single muscle fiber, the so-called single fiber action potential (SFAP). Classically, simplified dipole, tripole or quadrupole models have been employed for modelling the propagation of the action potential along a contracting muscle fiber (21; 30; 32; 41). A more general model has been proposed by Dimitrov and Dimitrova (15), and this model has been successfully employed, modified and combined with various other models for the remaining physiological processes in a number of publications (14; 17; 19; 55).
In the present article, the SFAP model originally proposed by Dimitrov and Dimitrova (15), and subsequently extended by Farina and Merletti (17), is combined with the well-known motor unit (MU) pool organization model of Fuglevand et al. (20) and the twitch force parameterization used by Raikova and Aladjov (43). Recent results regarding the modelling of MU rate coding and recruitment (10) and the variability of the inter-spike intervals (35) are incorporated. Care is taken in particular to achieve a consistent parameterization of the electrical and the mechanical components of the model, resulting in a realistic EMG-force relationship of the simulated muscle. Moreover, a novel nonlinear transformation of the rate coding model is proposed, which ensures that the desired output force is matched by the model. Finally, an alternative analytical formulation of the SFAP model of Farina and Merletti (17) is proposed, which renders the physiological meaning of the model more clear. Based on this alternative formulation, a proof is provided that in this model no fiber represents a net current source or sink at any point in time, which is a physiologically plausible property due to the quasi-static behavior of action potential generation (41). It is also in accordance with the predictions of the reknown Hodgkin-Huxley model (27).

Mathematical model
The fundamental functional unit of a skeletal muscle is the motor unit (MU), comprising a motor neuron and all muscle fibers innervated by that neuron. The following sections introduce mathematical models of the electrical and mechanical properties of MUs, as well as their organization in a muscle. Figure 1 shows a graphical summary of the main model components and their interactions and may provide a useful reference for the . . . EMG(t) Figure 1: Block diagram illustrating the main components of the proposed model of muscular force generation. The model input is the normalized desired muscle forceτ (t) with values in [0, 1], the model outputs are the total generated muscle force τ (t) and the measured EMG signal EMG(t). Firing rates of individual MUs are denoted by λ i (t), the current force contribution of each MU by τ i (t), and the current EMG contribution of each MU by MUAP i (t).
reader while following along the description of the model.

Motor unit pool structure
Every muscle consists of a number n of MUs. Each MU has various mechanical and electrical properties, most of which have been found to be closely related by means of the size principle (22): MU size as measured by the number of fibers contained in the MU is roughly proportional to force twitch amplitude, EMG twitch amplitude, and recruitment threshold. This means that larger MUs are only activated at higher levels of desired muscle force compared to smaller units, but they also add larger force and EMG contributions to the muscle output once activated. Recent results have shown that the electrical twitch conduction velocity v is linearly related to the recruitment threshold, and hence can be considered a size principle parameter aswell (11). Setting one of these parameters to a particular value for a given MU also determines the remaining parameters.
The recruitment thresholds appear to follow a continuous distribution with many MUs attaining a small recruitment threshold, and few large MUs only being recruited at high activation levels (44). This behavior is captured well by the exponential model proposed by Fuglevand et al. (20), which assigns the recruitment thresholds to MUs i = 1, . . . , n, whereτ i denotes the minimum fraction of total muscle force at which the MU is recruited, andτ * denotes the point of full recruitment, i.e., the relative level of total muscle force at which all MUs are recruited. Alternatively, following the formulation of De Luca and Contessa (8), the thresholds can also be modelled as where b denotes a scaling factor that influences the shape of the distribution. The latter model results in a more gradual slope compared to the first one and has been used by With the recruitment thresholds set, peak twitch forces are calculated as a linear function of the recruitment thresholds following as proposed by Contessa and De Luca (7), where the twitch peak range (P n /P 1 ) is typically large, e.g., P n /P 1 = 130 for the First Dorsal Interosseus (FDI) muscle (7).
Equivalently, the number η of innervated muscle fibers, which appears to be the main factor influencing MU twitch force (20), can be modelled directly proportional to the peak twitch force and hence also to the recruitment threshold, i.e., The same is true for the electrical twitch conduction velocity v, as noted above (11). Note that equation (4) determines the relation between the electrical and mechanical properties of a MU, as the number of fibers in a MU determines the amplitude of the electrical twitch response in this model. In equation (4), a linear relationship between the two has been assumed. However, it has been shown that in some muscles, the electrical twitch amplitude may rather be related to the square root of the force twitch amplitude (57); such a relation can easily be implemented into the model by modifying equation (4) to

Geometrical distribution of motor units and muscle fibers
Muscle fibers belonging to the same motor unit are distributed in a territory that may span a large portion of the muscle cross-section (3; 4; 53). The territories of the different MUs overlap, leading to fibers belonging to multiple MUs being intermingled (3; 53).
Motor unit territories have been found to attain an irregular round shape (3), whence we propose the use of an elliptic model for the MU cross-sections. With the elliptic axis ratio being fixed, the MU cross-sectional area-and thus the axes lengths-is calculated by dividing the number of innervated fibres η by the desired MU fiber density ρ (fibers/area): The midpoints of all MUs are then distributed uniformly over the muscle cross-section.
Note that without further assumptions, the above model directly leads to overlapping regions between MUs, which is a desirable feature, as noted above. Finally, fibers belonging to the MU are then again distributed uniformly inside the elliptic MU crosssection. The model is equivalent to the propositions of Fuglevand et al. (20), except for the division of the muscle into multiple parts to reduce fiber density variability, and the fact that they used circular MU territory shapes as opposed to the more flexible elliptic shape proposed here.
A key decision when modelling random MU placement concerns the treatment of muscle boundaries. Those parts of MU territories that exceed the muscle territory must be cut off and the question remains how to account for this loss in MU territory. Several approaches to solving the problem are conceivable (47): 1. All fibers belonging to the MU are placed in the remaining parts of the MU territory.
This approach leads to an increase in the fiber density of boundary MUs, and hence also in the overall fiber density towards the muscle boundaries.
2. The number of fibers innervated by the MU is reduced proportionally. This approach keeps the assigned fiber density constant but reduces the number of innervated fibers and thus the size of boundary MUs.
3. The axes lengths of the elliptic MU region are adjusted in such a way as to keep the MU area at the desired value in spite of the cut-off. This keeps the number of fibers and the fiber density at the desired values but likely leads to strongly increased overall fiber density towards the muscle centre, due to many adjusted MU regions overlapping there.
There are advantages and disadvantages to each approach, and it does not yet seem to be clear if one of the proposed approaches is generally superior to the others, or which approximates reality best (47). However, muscle fiber diameters appear to be approximately constant throughout a muscle (24; 49), whence a constant fiber density throughout the muscle cross-section seems desirable. To this end, the second of the above approaches has been pursued here. In order to avoid a high variability of the fiber density due to the random MU placement, it is advisable to divide the muscle cross-section into M parts of equal size, and then distribute n/M MUs uniformly in each part.
Reducing only the number of fibers in MUs close to the muscle boundary without adjusting the other MU parameters as well would disturb the relation between the electrical and mechanical properties of these MUs. To account for this disturbance, the recruitment thresholdτ i , the peak twitch force P i and the electrical twitch conduction velocity v have been recomputed following equations (3) and (4) for these MUs. Note that this, in turn, distorts the exponential distribution of MU parameters described by equations (1) and (2). This, however, was considered less grave than a disturbed electromechanical relationship in a significant number of MUs.
An alternative model of MU placement has been proposed by Navallas et al. (36).
They explicitly considered the optimization problem of minimizing the variability of muscle fiber density throughout the muscle, regardless of variances in MU fiber density and while maintaining the exponential relationship (1). Their model has recently been extended to account for regionalized MU placement (46). While all of these are desirable properties of a MU placement algorithm, the variability of the fiber density may also be reduced by the much simpler division of the muscle region into distinct parts and random placement proposed above; MU regionalization could be implemented using this same division into distinct parts; and the distortion to the exponential relationship (1) when using this algorithm has been found to be rather small in practice, see section 3.
In summary, both the placement algorithms of Navallas et al. (36) and the one proposed here represent viable modelling choices, with the former probably being preferrable if the particular influence of different aspects of MU geometry is of interest, and the latter being a much simpler algorithm.

Firing rates
Following the experimental results of Milner-Brown et al. (33), firing rates have traditionally been modelled as a linear function of excitatory drive (16; 20). This simplified model has recently been superseeded by the linear-exponential firing rate model proposed by De Luca and Hostage (10). In this model, the firing rate of a MU is described as whereτ = τ /τ max ∈ [0, 1] denotes the relative level of total desired muscle force τ , and the remaining constants are shape parameters. This model reproduces many phenomena observed experimentally, such as the onion-skin phenomenon: MUs recruited first appear to attain a higher firing rate throughout the whole contraction than those MUs recruited later (9). Interestingly, however, the model also misses an essential feature of the model proposed earlier by Erim et al. (16), namely, the (experimentally observed (9)) increased slope of the firing rate characteristics at excitation levels exceeding the pointτ * of full recruitment. We will come back to this issue in section 2.8. For example values of the shape parameters for different muscles, refer to De Luca and Hostage (10).
Note that although the model is defined as a function of the desired total muscle force, it does not guarantee in any way that this force level is actually attained by the muscle. This inconsistency will be addressed in section 2.8 of this article by introducing a suitably defined nonlinear input transformation act(τ ) that is then used as the input to the rate coding model presented above, instead of using the desired muscle force levelτ directly. One benefit of the nonlinear input transformation proposed in section 2.8 is that it effectively results in an increased slope of the firing rate characteristics at excitation levels exceedingτ * , as has been observed experimentally (9).

Firing instants
Given the instant t i(j−1) of the last firing of MU i and the time course of the MU's firing rate λ i , the next firing instant t ij can be calculated. To model the stochastic distribution of the inter-spike intervals, these are assumed to follow a normal distribution with a coefficient of variation that decreases with increasing activation, following as proposed by Moritz et al. (35). The j th inter-spike interval ISI j then is drawn from where the mean t * ij of the time of the next firing event and the corresponding mean inter-spike interval ISI * j are obtained by solving for t * ij as proposed by Fuglevand et al. (20), withτ (t) denoting the normalized requested muscle force level at time t. The time t ij of the j th firing event is then calculated as Note that several distributions other than the normal distribution have been proposed for modelling the distribution of the inter-spike intervals (1; 25). For reasons of simplicity, and as the influence of the shape of the distribution of the inter-spike intervals on the overall EMG and force signals found by Barry et al. (1) appeared to be rather negligible, a normal distribution is used in the simulation described in section 3.

Intracellular action potential propagation
The propagation of an intracellular action potential (IAP) from the neuromuscular junction (NMJ) of a muscle fiber along both directions towards the two fiber ends can be modelled by representing the actively firing fiber by a distributed current source and sink.
In the model originally proposed by Dimitrov and Dimitrova (15), this distributed fiber membrane current sourceî(z, t) is composed of two propagating wave fronts and localized contributions at the NMJ and the two fiber ends. These localized contributions model the IAP generation and extinction process. In the formulation of Farina and Merletti (17), the model readŝ Here, z denotes the spatial variable along the muscle fiber, z i the location of the NMJ, L 1 and L 2 are the distances between the innervation zone and the right and left tendon, respectively, and v denotes the IAP's propagation velocity. Moreover, denotes the voltage gradient across the fiber membrane along the fiber axis, where the function V m (z) prescribes a model for the trans-fiber membrane voltage wave shape and can be chosen arbitrarily to match simulated or measured data. Refer to Plonsey and Barr (41) for details on the significance of V m (z). Here, the analytical model function with D 1 = 96 mV mm −3 and D 2 = −90 mV will be used, as originally proposed by Rosenfalck (48) and as has been done by Farina and Merletti (17).
The IAP model in equation (12) can be shown to be equivalent to choosinĝ with the Dirac distribution δ, the end-of-fiber components and and the potential generation component Here, again, the end-of-fiber components describe the IAP extinction process at the fiber ends, and the potential generation component models the influence of IAP generation at the innervation zone on the membrane current.
This formulation renders -to the authors' opinion -the structure of the model more obvious, by clearly distinguishing between propagating and non-propagating signal components, and by revealing the non-smoothness of the resulting distributed current source, the latter following from the presence of the stationary Dirac distributions at the two fiber ends and the location of the innervation zone. The equivalence of the two formulations of the model is summarized in the following lemma, the proof of which is given in the appendix. (12) and (15) to (18) are equivalent,

Lemma 1. The expressions given in equations
Interestingly, one can also show that this IAP model ensures which implies that the fiber does not represent a net current source or sink at any point in time. This property is well motivated by physiology, considering that the invoked electrodynamical processes can be considered quasi-static (41), and it is also in accordance with the predictions of the reknown Hodgkin-Huxley model for action potential propagation (27). It is precisely a result of the presence of the three Dirac distributions in the model equation (15), as these have the combined effect of collecting all remaining currents exerted by the intermediate fiber sections due to their higher potential. This result is the subject of the following lemma, the proof of which again is deferred to the appendix.

EMG measurements
Biological tissues can be considered volume conductors (41). The existence of an electric field implies the existence of electric currents travelling through the tissue, and vice versa. 3 Due to the comparably low rate of change of physiological systems, it is justified (41) to assume these time-varying electric fields to behave as if they were static at each instant of time, whence they are called quasi-static. This assumption amounts to a neglection of the capacitive properties of the tissues. Accordingly, as for static fields, the electric field in a physiological volume conductor is considered equal to the negative gradient of a scalar potential ϕ, namely, By Ohm's law, the current density (current per unit of cross-sectional area) in a volume conductor is proportional to the electric field, that is, compactly supported ψ(z) ∈ C ∞ , but equation (14) could trivially be modified to comply with this requirement, e.g. by convolution with a smooth and compactly supported mollifier. As this is a purely theoretical operation with no practical relevance at all, this has not been pursued here. 2 See previous footnote. 3 Note that in biological tissues, the charge carriors are ions, as opposed to electrons in electric wires (41, p. 25).
where σ denotes the conductivity of the medium. Defining a distributed current density source I throughout the region of interest, the divergence of the current density is constrained by Combining equations (21) and (22) and assuming a homogeneous, isotropic medium yields Poisson's equation for the diffusion of the potential, namely, In the following, the electric field generated by point sources in planar tissue layers will be considered as a model for flat and large muscles, such as the recti abdominis simulated in section 3. The muscle layer is assumed to be infinitely extended and planar, and to be covered by an infinitely extended planar layer of fat and an infinitely extended planar layer of skin. Muscle tissue is considered anisotropic in order to reflect the difference in conductivity between currents along the muscle fiber axis and currents across the muscle fiber axis, whereas fat and skin tissue are considered isotropic. Muscle fibers are assumed to run along the z direction, with the x and z dimensions spanning the skin plane, and the y dimension being orthogonal to the skin plane, positive vectors pointing outwards.
The geometrical set-up described above has been analyzed by Farina and Rainoldi (18). For a point source of strengthÎ located at (0, y 0 , 0), the authors derive the 2-D spatial Fourier transform of the resulting potential distribution at the skin surface to be with the abbreviations and where ω x = 2πf x and ω z = 2πf z denote the spatial angular frequencies in the x and z directions, respectively. The coefficients EMG measurements are usually taken differentially between a set of electrodes. Consider a regular grid of R×S electrodes with interelectrode distances d x and d z , respectively, where R = R a + R b + 1 and S = S a + S b + 1. The variables subscripted by a and b denote the number of electrodes on the two sides of an arbitrarily chosen reference electrode.
The grid is assumed to be aligned parallel to the z axis. Assigning weights ζ k to the electrodes and assuming all electrodes to attain the same transfer function, the (spatial) transfer function from a given surface potential distribution to the potential measured by such an electrode configuration at each point on the surface is given by (17) H For the transfer function of a single electrode, arbitrary model assumptions can be made.
Concatenating the spatial transfer functions H vc of the volume conductor, H ec of the electrode configuration and H ele of the electrodes themselves, the global transfer function of the combined system emerges as From this, the 2-D potential distribution on the skin surface can generally be calculated where i(ω x , y, ω z , t) = F xz {i(x, y, z, t)} is the 2-D Fourier transform of the current density source i(x, y, z, t), and * (x,z) denotes 2-dimensional convolution in the x and z variables.
For a particular electrode location on the skin surface and a muscle fiber following a straight line parallel to the skin surface, equation (31) simplifies, and the resulting single-fiber action potential (SFAP) ϕ(t) can be calculated numerically (17; 38). One can prove that in this case the integration kernel only has removable singularities, which ensures the convergence of a numerical integration scheme (38). Figure 2 shows exemplary SFAPs resulting from the evaluation of equation (31) for such fibers using nested numerical integration schemes. Note that while the above derivation has been performed for the case of planar volume conductors, similar models have been derived for cylindrical volume conductors -much more appropriate for the simulation of limb muscles -as well (19).
The numerical solution of the (simplified version of) equation (31) where n denotes the number of MUs, N i the number of firing events of that MU, and the motor unit action potential, which is obtained by summing over the contributions, i.e., the SFAPs, of all muscle fibers belonging to MU i. Figure 2 shows exemplary SFAPs simulated using a numerical implementation of the above model. The simulated SFAPs correctly reproduce the dependency of the SFAP shape on the relative position of the recording electrode and the depth of the muscle fiber, as well as the distinction between propagating and localized signal components at the NMJ and the two fiber ends. In particular, the experimentally observed end-of-fiber extinction signals are included.

Force twitches
Each MUAP generates a corresponding force contribution, denoted as a force twitch f i (t). Most previously proposed muscle models (13; 20) have employed the force twitch parameterization of Milner-Brown et al. (33), where a force twitch is completely described by its twitch rise time T ri and its peak twitch force P . This model, however, does not allow setting the half relaxation time T hr independently of T ri , which is essential for modelling different muscle fiber types (refer to Fig. 1.2 in (30) and which results in the following force twitch model: with and κ = log 2 This model satisfies f i (T ri,i ) = P i , and f i (T hr,i ) = P i /2. The parameter T EMD denotes the electromechanical delay between the onset of electrical and mechanical activity of the motor unit. The three parameters T ri , T hr and T EMD are sampled from a Weibull distribution for each MU, as proposed (for the former two) by Contessa and De Luca (7).
with γ 2 and γ 1 constant muscle parameters, as proposed by Contessa and De Luca (7). Here, g ij denotes the gain factor assigned to the jth firing of MU i, andλ ij is the normalized instantaneous firing rate at that firing event: Finally, the total force generated by a muscle is calculated as the superposition  Figure 1 illustrates the model of muscular force generation described so far, introducing a muscle activation signal act ∈ [0, 1]. Motor Unit firing rate models such as the model of De Luca and Hostage (10) described in section 2.3 usually define the firing rate as a function of the desired normalized muscle forceτ , i.e., they choose act ≡τ . This is a consequence of the fact that these models are derived from experimental measurements of MU firing rates at different muscle force levels. Now, defining the firing rates λ i of all

Excitation-force relationship
MUs (e.g. as described in section 2.3) and also defining the force generating properties of these MUs (e.g. as described in section 2.7) uniquely determines the generated muscle force τ , see fig. 1. It is, however, by no means guaranteed that the generated normalized muscle forceτ will be equal to the force levelτ that has been used as an input to the firing rate model. In other words, the generated muscle force output does not match the desired muscle force. include some remarks on physiological feedback processes, neither model was meant to replicate properties of actual physiological feedback control, but rather to account for the input-output inconsistency of the respective model. Although there certainly is some feedback element in physiological force control, it would appear reasonable that this feedback is mainly necessary to deal with external disturbances and changing muscle properties, not to account for a static input-output inconsistency. For these reasons, we propose a novel solution to this problem which does not require the introduction of a feedback loop. Our approach is based on a static input nonlinearity.
Denoting the desired normalized muscle force byτ ∈ [0, 1], our desire is to choose whereτ (act) is the normalized generated force output. This can be achieved by calculating τ = ξ(act), i.e., the normalized generated muscle force as a (nonlinear) function ξ of the activation input. Utilizing act = ξ −1 (τ ) (42) as the activation input to the firing rate model then yields a simulation model that satisfies condition (41). Averaging over individual firing events, the mean generated muscle force is given by where with Γ(x) the Gamma function, denotes the total impulse generated by a single force twitch f i (t) of MU i, and g i denotes the nonlinear force gain factor defined in equation (38). By evaluatingτ for different values of the activation level act, one can determineτ = ξ(act). Employing act = ξ −1 (τ ) as the activation input to the firing rate model then yields a simulation model that satisfies τ (act(τ )) τ max =τ ∀τ ∈ [0, 1] in the mean over time. Condition (46) ensures that while there may be differences between desired and generated muscle force at individual time instants due to the stochastic nature of the firing instants (see section 2.4), the two forces agree on average. Finally, note that the introduction of a static input nonlinearity leads to a distortion of the firing rate model, which now no longer receivesτ as an input, but rather its nonlinear transformation ξ −1 (τ ). The influence of this distortion will be exemplified in the following section, and its plausibility will be discussed.
R programming language (42). In total, 300 MUs have been simulated, each consisting of between 30 and 150 muscle fibers, organized into three separate muscle bellies of each of the two recti (left and right). Single differential detection was assumed for the EMG signals, with the two electrodes placed on the linea alba between the second and third belly. Where available, parameter values available in the literature have been chosen (12; 45; 52). A series of constant-force isometric contractions at progressively increasing force levels has been simulated, each lasting for 4 s. In order to assess the effect of the the inherent randomness of the MU placement algorithm, the muscle generation procedure and the subsequent simulation each have been executed seven times. The simulated signals are available at the Dryad data repository (39).

Discussion and conclusion
In this article, a comprehensive mathematical model of surface EMG and force generation during voluntary isometric contractions in skeletal muscles has been described. The

Conflict of interest statement
Part of the research leading to the presented mathematical model was financially supported by Drägerwerk AG&Co. KGaA, Lübeck, Germany.

Author contributions
Eike Petersen und Philipp Rostalski contributed to the conception of the work. The resulting mathematical model was derived and implemented by Eike Petersen under the supervision of Philipp Rostalski. Eike Petersen wrote the first draft of the manuscript.
All authors contributed to manuscript revision and approved the submitted version.
Proof. For compactly supported ψ(z), generally holds. Furthermore, and equivalently for EOF 2 . Combining everything yields which concludes the proof.

Appendix B: Distributions
A distribution T is a continuous linear mapping T : D → R, where D is a given set of so-called test functions. In the following, the value of a distribution T acting on a test function ζ shall be denoted by the duality T, ζ , with test functions chosen from the set D of compactly supported, smooth functions. Furthermore, two basic properties of distributions are ηT, ζ = T, ηζ ∀ η ∈ C ∞ , ζ ∈ D and T , ζ = − T, ζ ∀ ζ ∈ D.