Multiple-Order Singularity Expansion Method

Physical systems and signals are often characterized by complex functions of frequency in the harmonic-domain. The extension of such functions to the complex frequency plane has been a topic of growing interest as it was shown that specific complex frequencies could be used to describe both ordinary and exceptional physical properties. In particular, expansions and factorized forms of the harmonic-domain functions in terms of their poles and zeros under multiple physical considerations have been used. In this work, we start from a general property of continuity and differentiability of the complex functions to derive the multiple-order singularity expansion method. We rigorously derive the common singularity and zero expansion and factorization expressions, and generalize them to the case of singularities of arbitrary order, whilst deducing the behaviour of these complex frequencies from the simple hypothesis that we are dealing with physically realistic signals.

Physical systems and signals are often characterized by complex functions of frequency in the harmonic-domain.The extension of such functions to the complex frequency plane has been a topic of growing interest as it was shown that specific complex frequencies could be used to describe both ordinary and exceptional physical properties.In particular, expansions and factorized forms of the harmonic-domain functions in terms of their poles and zeros under multiple physical considerations have been used.In this work, we start from a general property of continuity and differentiability of the complex functions to derive the multiple-order singularity expansion method.We rigorously derive the common singularity and zero expansion and factorization expressions, and generalize them to the case of singularities of arbitrary order, whilst deducing the behaviour of these complex frequencies from the simple hypothesis that we are dealing with physically realistic signals.

I. INTRODUCTION
The model of Linear and Time-Invariant Systems (LTIS) is commonly used in physics to derive the response of a medium to a excitations.Such systems are usually studied in the harmonic domain where they are associated with complex-valued transfer functions that can be used to fully describe the response of the system to an arbitrary excitation [1,2].This formalism has been mainly developed in the research field of automatic system control but it can be found in a wide variety of problems where it goes by different names.The complex impedance and admittance formalism, for instance, is used to characterize the properties of materials in electronics, acoustics and biology [3][4][5].The Scattering matrix formalism, first introduced in quantum electrodynamics [6], can be used to link outgoing waves to incoming waves, or final states to initial states in scattering problems in various situations (transport phenomena, diffraction gratings, integrated circuits, chaotic systems, ...) [7][8][9][10].
While emphasis is generally placed on real and positive frequencies with clear physical interpretations, the description of a linear transfer system together with its input and response signals often requires an analytical continuation into the complex frequency plane [11,12].Instead of considering individual responses at specific real frequencies, complex frequencies can provide us with intuition as to how a system will behave or how a signal is shaped over a large spectral width.The analysis of the electromagnetic response in terms of complex zeros and singularities has turned out to be highly efficient for several applications such as quantum waveguides [13], highly selective filters [14], plasmonic metasurfaces [15], anapoles [16], coherent perfect optical absorbers [17][18][19][20][21] and analog computing [22].In addition, studying the or-der of the zeros and singularities can help in providing better interpretations of the associated phenomena [23].
Any real temporal signal h(t) can be associated with a harmonic-domain function H(ω), which is either an input or output signal, or a transfer function linking the two and describing an LTIS as a filter acting on the input in the harmonic domain (see Fig. 1).Different methods have been developed to express H with respect to its zeros and singularities in the complex frequency plane.The most relevant methods, in our case, can be traced back to the Weierstrass and Hadamard factorization theorems for holomorphic and meromorphic functions respectively, which provide general expressions of H as opposed to local expressions such as the Laurent series expansion.
(i)-With the rapid development of electronics and automated machines, the control of the LTIS stability often favoured the factorization of the transfer function H, in the simple form of the ratio of complex polynomials [24], in order to monitor the evolution of the phase and amplitude with respect to the frequency in Bode diagrams [25].In electronics and automatic system control, singularities and zeros of multiple-order are not uncommon (n th order Butterworth filters [26] for instance).
(ii) Alongside this progress, pioneer works led to the development of the Singularity Expansion Method (SEM) [27][28][29].It was first developed to approximate the time dynamics of systems associated with arbitraryorder singularities in the harmonic-domain and helped in describing many problems in electromagnetism [29][30][31][32][33] (an exhaustive set of references dating back to the early developments of the SEM can be found in ref. [34]).The SEM has received a renewed interest in the recent years in the case of simple poles, i.e. singularities of order 1, where accurate expressions have been derived and applied to various problems [35][36][37][38][39][40].In addition, pole or singularity expansions have been increasingly used in the frame of quasi-normal modes or resonant state expansions [41][42][43][44][45].However, exact expansions have not been obtained in the case of poles of arbitrary order, despite the prospects they offer in describing more complex systems and/or their input and output signals [46][47][48].
FIG. 1: General representation of an LTIS.(a) In the temporal domain, the system is characterized by its impulse response h(t) which can be used to obtain the output y(t) after a convolution with the input x(t).(b) In the harmonic domain, the LTIS is described by its transfer function H(ω).It is a filter acting in the frequency domain on the input X(ω) to generate an output Y (ω).In this case, H(ω) is a low-pass filter which partially removes the noise present at higher frequencies, as it can be seen by comparing the temporal input x(t) to the temporal output y(t).The output Y (ω) is the product of H(ω) and X(ω).The harmonic domain functions X(ω), H(ω) and Y (ω) are obtained by bilateral Laplace transform of the temporal signals x(t), h(t) and y(t).The temporal signals can then be recovered using an inverse Laplace transform.
These two methods, i.e. the singularity expansion method for simple poles and the multiple-order pole and zero factorization used in system control and electronics, can be obtained under certain restrictive considerations.Furthermore, the derivation of one method from the other and their relationships are not obvious, in spite of the same singularities appearing in both the factorized and expanded expressions.
In this work, we first describe and explain the usual state of the art expressions of the pole expansion in the case of simple poles, and the pole and zero factorization as depicted in the fields of electronics and automatic system control.We present some of the limits of these expressions which motivated the need for the gen-eral expressions of the singularity expansion and the singularity and zero factorization which are then rigorously derived using complex analysis theorems applied to the harmonic-domain signals at play under simple hypothesis.We show that the behaviour of any function around its discrete set of singularities and/or zeros is enough to describe that function in the complete complex frequency plane.In addition, we put the emphasis on the constraints of the distribution of the singularities and zeros which arise from the general, non-restrictive consideration of physically realistic temporal signals, i.e. real and causal signals.Finally, we show how to obtain, from these harmonic-domain expressions, the temporal singularity expansion which holds information regarding the stability of the systems, the convergence of the signals, and the transient and steady-state temporal dynamics.

A. Simple Pole Expansion
The first singularity expansions were developed in the aforementioned SEM [27][28][29] based on the observation of a system's response to a sinusoidal input in the transient regime.It was shown to be a combination of damped sinusoidal functions associated with complex frequencies which were the singularities of the transfer function of the system.While Baum first derived the singularity expansion by taking into account the order of these singularities, most applications made use of (and later demonstrated) this expansion in the case of simple poles, i.e isolated singularities of order 1 [35,36,38].The resulting expansion is referred to as the Simple Pole Expansion (SPE), and it has the following expression: where p denotes the simple poles of H.The SPE is composed of a background term H bgrd (ω) which does not possess poles [42], and a sum of resonant terms Res(H, p)/(ω−p).This translates the idea that the shape of H(ω) in a specific frequency range is chiefly affected by the nearby singularities, with the background term compensating for the offset introduced by the singularities with respect to the average (or background) value of H in that range.When all the poles p of H are known, the background term reduces to a constant.The contribution of the poles to the shape of H is illustrated in Fig. 2 for H defined as the reflection coefficient of a thin layer of silver illuminated from one side at normal incidence: where d = 70 nm is the thickness of the silver layer, c is the speed of light in the air, n(ω) is the refractive index of silver (see the SI for the detailed expression), r(ω) = (n(ω)−1)/(n(ω)−1) is the Fresnel reflection coefficient at the air/silver interface, and t(ω) = 2/(n(ω) + 1) and t (ω) = 2n(ω)/(n(ω) + 1) are the Fresnel transmission coefficients at the air/silver and silver/air interfaces respectively.The zeros, the poles and the residues of H were determined numerically.Each pole p in the complex frequency plane, in Fig. 2 (b) is associated with a resonant term Res(H, p)(ω − p), whose moduli are plotted in Fig. 2 (a).The sum of the resonant terms associated with the poles in the plotted complex frequency window gives the red curve in Fig. 2 (a), which matches the local shape of the exact response (dashed line which corresponds to Eq. (2), i.e. the approximate position of the local minimum and maximum frequencies.The poor match between this sum of resonant terms and the exact expression is due to the missing background term H bgrd (ω) which should correct for the influence of the poles outside of that frequency window.Since most physical systems are described by poles of order 1, the SPE is well suited for the description of the transfer function of such systems, or their response to a temporal or spatial sinusoidal input [36,38,39,49,50].Other expressions must be used when the order of the poles is increased.One commonly used alternative in electrical and electrics engineering is what we will refer to as the Pole and Zero Factorization (PZF).

B. Pole and Zero Factorization
The PZF is obtained by making the assumption that the harmonic-domain function H is the ratio of two polynomials of the complex frequency variable ω [1,2,24], and that any system responds instantaneously to an input signal: where iz and ip are the zeros of N and D respectively, with potentially identical zeros and poles (in which case their order is higher).Let us point out that we consider, in this case, the variable iω instead of ω for N and D in order to remain consistent with the usual definition of the Laplace transform formalism.Using a partial fraction expansion, it is possible to write an expression similar to the SPE: where H bgrd = α ,0 and the poles p are considered with their order or multiplicity ν .The constants α ,1 can be identified as the residues, and they can be calculated with the zeros and poles using Eq. ( 3) [27,36]: As stated in the introduction, the PZF and the resulting pole expansion are widely used in electronics and system control where they provide highly accurate results when studying systems at time scales larger than  2), as well as its resonant terms Res(p )/(ω − p ) associated with the poles p l and the modulus of the sum of these resonant terms.The frequencies range from 1.040 × 10 1 6 Hz to 1.093 × 10 1 6 Hz, which is equivalent to ultraviolet wavelengths between 172.9 nm and 181.3 nm.(b) Log-amplitude of H(ω) in the complex frequency plane, in a limited complex frequency window.The poles in (a) are highlighted in (b) as red points.The sum over the poles gives an approximation of the shape of H(ω) on the real axis (the minimum and maximum frequencies), but with the background term of Eq. ( 1) omitted, the reconstructed red curve poorly matches the exact expression.
the transient time scale.However, this is not the case for other fields such as wave physics, in which case the PZF is missing a complex exponential factor depending on the frequency ω which can be interpreted as the result of a time delay required by a system to respond to different frequencies (i.e. a dispersive system) [36,38].This term, which is particularly important for systems which interact with signals with respect to both time and space cannot be obtained from the restrictive hypothesis that H is a ratio of two polynomials.In the next section, we present a singularity expansion taking into account the potentially infinite order of the singularities.From this expression, a more general singularity and zero factorization will then be derived.

A. Multiple-Order Singularity Expansion
We now consider a meromorphic function H, i.e. holomorphic everywhere on C except for a set of points P which is the set of its singularities which are all assumed to be poles or isolated essential singularities.Using the Cauchy integration theorem and the residue theorem, the function H can be expressed as the sum of an integral term and an expansion on its set of singularities (see Eqs. (S1 -S5) in the SI): where γ is a closed curve around the singularities {p} of H, and Res(F, p) is the residue of F at the singularity p.
The residues of F can be analytically determined from the Laurent series coefficients α(H, p, n) of H around the poles p, which provides local information regarding the behaviour of H in the vicinity of its singularities.This leads to the following expression (details of the calculations can be found in Eqs.(S6 -S17) of the SI): where ν p is the potentially infinite order of the singularity p. Eq. 7 shows two contributions to the behaviour of H at the frequency ω: (i) the nearby singularities {p}, within the closed curve γ; (ii) the set of values of H on γ, represented by the integral term, which accounts for the contribution of all the singularities outside of γ.
The bigger the closed curve gets, the less the integral affects the value of H since the singularities outside the closed curve become too far from ω, provided that H does not grow faster than |ω| in the complex plane (which is usually the case for physically realistic signals).When this hypothesis holds, we show that the integral term can be replaced by a constant by replacing the curve by a circle of infinite radius (Eqs.(S18-S24) of the SI): where H NR is a constant non-resonant term, and H R (ω) is a resonant shaping term which depends on the frequency ω.The non-resonant term H NR is defined as: where a an arbitrary complex frequency which is not a pole.Let us stress that the choice of a has no influence over the value of the constant term H NR .The resonant term H R (ω) possesses the following expression: If 0 is not a pole, we usually set a = 0 and the nonresonant term H NR can be expressed using the static response H(0): We refer to Eqs. ( 8) to (9) as the Multiple-Order Singularity Expansion Method (MOSEM).The accuracy of MOSEM is shown in the case of poles of order 2 in a purely theoretical example in the SI.When the order ν p of all the poles is 1, the expression becomes: Let us stress that the residue of H associated with the pole p is defined as Res(H, p) = α(H, p, −1).Therefore by identifying H NR as a constant background term H bgrd when all the singularities are taken into account in the resonant terms, we obtain the SPE expression defined in Eq. ( 1): If H is reconstructed using only a finite set of singularities, two strategies can be adopted to still obtain a good match between the resulting truncated SPE or MOSEM and the exact expression; (i) as mentioned before, H bgrd can be considered as a holomorphic function, i.e. with no singularities, which models the influence of the singularities missing from the resonant term H R (ω); (ii) the influence of the poles outside of the region of interest can be accounted for by considering an additional virtual singularity within the resonant term as was done in ref [39].

B. Singularity and Zero Factorization
The Weierstrass and Hadamard factorization theorems state that any meromorphic function can be written as the ratio of two complex polynomials multiplied by a complex exponential.Starting from this consideration, expressions such as the PZF (in which the complex exponential is usually removed) can be obtained (although not straightforwardly) to study the phase and amplitude variations at real frequencies.The poles and zeros can be defined in the resulting rational fractions as the zeros of the denominator and numerator respectively.Alternatively, a factorization involving the singularities and zeros can be obtained from MOSEM, providing a clearer link between the expanded and factorized forms as was done in Refs.[15,36] in the case of poles of order 1.
Here, we aim at deriving a generalized factorized expression from MOSEM expression of the function H in Eqs.(8) to (9).H can always be written as: with m ≥ 0 the order of the zero 0 of H, and G a meromorphic function which does not possess 0 as a zero.If m = 0, we have G = H.Let us consider F the logderivative of G, with G the complex derivative of G: G is calculated by taking the derivative of MOSEM applied to G (Eqs. ( 8) to ( 9)).The poles of F are the zeros z and the singularities p of G, and they are all of order 1 ((see Eqs.(S33 -S36) in the SI): In addition, their associated residues are the order ν z and ν p of the zeros and poles of G respectively : We apply MOSEM to F taking these values into account: with p ∈ {z , p }, ν p = 1, and α(F, p, −1) Finally, we obtain the following expression of F : By integrating F on a curve from the arbitrary complex frequency a to the frequency of interest ω, and applying the exponential function to the result, we derive the following expression (Eqs.(S37 -S41) in the SI): where G(a) is defined as: and G (a) is defined as: We refer to Eq. ( 21) as the Singularity and Zero Factorization (SZF).Let us point out the presence of the aforementioned complex exponential e −iτ ω missing from the PZF, but also of the known response G(a) (the static response in the case of a = 0).Let us also stress that in the case of an LTIS, the phase introduced by the complex exponential can be set to an arbitrary position by changing the time origin or the space-origin (for a system with coupled space and time variables).Shifting the time origin by τ 0 in the temporal domain results in a multiplication by e iωτ0 of the input X(ω) in the harmonic domain.This is tantamount to considering the transfer function H(ω)e iωτ0 with the same input X(ω) as depicted in Fig. 3 for H defined as: In Fig. 3 (a), the phase of the transfer function H(ω) is shown in the complex plane.In Fig. 3 (b), it is plotted for H(ω) multiplied by e iωτ0 with τ 0 arbitrarily set to π/6.We show that the position of the zeros and singularities remains the same after switching from H(ω) to H(ω)e iωτ0 , but a non-constant phase-shift is introduced in the complete complex frequency plane by the phasor e iωτ0 .It is thus possible to set τ 0 in such a way that τ + τ 0 = 0, where τ is the time constant naturally appearing in the SZF in Eq. (21).
By setting τ = 0, a = 0, and appropriately defining two constants N 0 and D 0 appropriately, we recover the PZF described in the previous section by considering only N s poles and N z zeros: The MOSEM and SZF expressions, which are reminded in Fig. 4, are two equivalent means to characterize a function: (i) MOSEM expression, which relies on the behaviour in the vicinity of the singularities only and thus depend on the singularities and Laurent series coefficients associated with them, is useful to get a fast and accurate approximation of a function in the harmonicdomain; (ii) the SZF is more convenient to look at the phase and better understand the behaviour on the real frequency axis.It is easy to obtain the Laurent series coefficients, and thus MOSEM expression, from the SZF using the definition of the Laurent series coefficients: where p 0 is a pole of order ν p0 of H, and H(ω) is given by the SZF in Eq. (21).It is however more difficult to obtain the zeros from MOSEM expression, although they can be approximated in a specific frequency range by writing MOSEM as a rational function involving only the singularities in the selected complex frequency window and solving for the zeros of the numerator.

IV. HERMITIAN SYMMETRY AND CONSTRAINTS IN THE HARMONIC DOMAIN
We now wish to take advantage of the physical nature of the signals to derive some constraints on the parameters of MOSEM and SZF.

Let us consider a real-valued function h(t) in the temporal domain, associated with a complex function H(ω).
Since h is real, H possesses a hermitian symmetry in the complex plane [36,51]: This property leads to constraints on the distribution of the poles and singularities in the complex frequency plane.By evaluating the complex conjugate of the singularity expansion in Eq. ( 8) evaluated at −ω * , it can be shown that for any pole p of order ν p , −p * is also a pole of order ν p .We can also obtain the Laurent coefficients associated with −p * to those of p via: and we show similarly that the non-resonant term H NR must be real.The Laurent series coefficients of H at −p * are the opposite of the complex conjugate of those at p.This leads to the following MOSEM expression: Let us point out that Eq. ( 28) restricts α(H, p, −m) to iR if p ∈ iR.The poles and Laurent series coefficients thus always come in pairs, as shown in Fig. 5 for the transfer function H defined in Eq. ( 2).In (a) and (b), the logamplitude in a complex frequency window as well as its symmetric window is plotted, highlighting the symmetry of the amplitude and thus of the singularities.In (c) and (d), the same frequency windows were chosen for the phase plots.They show the antisymmetry of the phase relative to the imaginary axis in the complex frequency plane, which is tantamount to an antisymmetry of the Laurent series coefficients.Let us stress that the hermitian symmetry only arises from the fact that we have real-valued signals in the temporal domain.Considering a causal plane wave e −iω0t u(t), with u the Heaviside step function and ω 0 > 0 is not equivalent to considering that time flows backward.It only means that the variations of the phase are opposed to those of e iω0t u(t).

B. Singularity and Zero Factorization
If we now look at the complex conjugate of the SZF (Eq.( 21)) evaluated at −ω * , we can show that if z is a zero of order ν z of H, then −z * is also a zero of order ν z .In addition, let us show that the Hermitian symmetry forces the time constant τ to be real-valued.If τ were complex-valued, it could be written τ = τ R + iτ I , with τ R and τ I the real and imaginary parts of τ respectively.Since the opposite of the complex conjugate of zeros and poles are also zeros and poles of the same order, the Hermitian symmetry would thus lead to: This condition would not be satisfied as ω tends towards +∞ unless τ I = 0. Therefore, τ ∈ R. If we set a = 0, the SZF thus yields the following expression of τ : with G (0)/G(0) ∈ iR.The time constant τ is thus the sum of two contributions: (1) the phase-shift introduced by the imaginary part of the singularities and the nonnull zeros, (2) a constant term depending on the static response of the derivatives of H.In physical systems for which the space and time variables are coupled, the phase shift can be modified by moving the spatial or temporal origin.In this case, the zeros and poles are identical, and only the constant term linked to the static response is changed.

C. Stability and Causality
Stability and causality are linked but distinct principles which can both be expressed in terms of the position of the singularities in the complex plane depending on the convention used to perform a Fourier transform [52].Causality states that any signal h must be generated at a certain time t h , and that it cannot depend on its future values.If h is a causal signal, it can therefore be written, using the Heaviside step function u, as: As long as h does not diverge faster than an exponential function for a long time t, it can thus always be regularized using a function h γ , γ > 0, which converges: Let us point out that we set t h = 0 in u, and we can do so without losing in generality.By construction, h γ possesses a Fourier transform H γ from which the harmonic domain function H associated with h can be determined: h can be retrieved by integrating H (multiplied by a complex exponential) over a horizontal line within the region of convergence of H, lower-bounded by the amplitude of the smallest diverging exponential function diverging faster than h (see Fig. 6).Using the residue theorem, it can thus be shown that any pole possessing a positive imaginary part in the complex frequency plane is associated with a causal diverging exponential function in the temporal domain, i.e. an unstable signal.Therefore, the poles of the harmonic-domain function associated with any stable signal must have a negative imaginary part.
The stability of the signals can be interpreted with energy considerations.As an input signal interacts with a system, it exchanges energy with it.This leads to a modification of the input signal which result in the output signal.For a passive system, the energy is transferred from the input to the system.The output thus has a lower energy and cannot diverge if the input is stable.For the output to diverge or become unstable, it must result from a sufficient energy transfer from the system to the input, or the interaction with an already diverging input with the system.Therefore, the only way to obtain an unstable output is through a high-energy, unstable input, or an active system.In terms of singularities, this means that the singularities of the transfer function of a passive system always possess a negative imaginary part (using our Fourier transform convention).If we inject energy into the system, i.e. the system is active, we move the singularities closer and closer to the real axis until the system is unstable and at least one singularity possess a positive imaginary part.

V. TEMPORAL EXPRESSIONS WITH THE MULTIPLE-ORDER SINGULARITY EXPANSION METHOD
The expansion obtained with MOSEM can be used to derive an analytical expression in the temporal domain [38,39].For physical systems, unstable behaviours might appear, which prevent the use of the inverse Fourier transform [53].It is therefore preferable to use the more general inverse Laplace transform to retrieve the temporal dynamics of a system or its response (and/or input) [52].In this section, we derive a generalized expression of the temporal domain signal h knowing its singularity expansion.Similarly to the previous section, if h is physically realistic, then causality implies that : Let us point out that if h is the response of an LTIS, the time constant t h corresponds to the time it takes for the system to interact with the input signal and produce the response.We shift the time origin to set t h = 0.The Laplace transform L and the Fourier transform F are equivalent in this case: Using these conventions, more properties regarding the inverse Laplace and Fourier transforms can be deduced.Let us consider the temporally diverging (or unstable) signal h(t) increasing slower than an exponential function in Fig. 6 (a).The Laplace or Fourier transform of h(t) has at least one pole in the upper half of the complex plane.This pole, which is called ip I in Fig. 6 (b), imposes the Region Of Convergence (ROC) of the Laplace transform of the signal (green band at the top).The inverse Laplace transform is defined as an integral over a horizontal in the complex plane.If we choose that horizontal line above the imaginary part of all the poles, within the ROC, we retrieve the original signal h(t) (from (b) to (a)).Otherwise, since only the contribution of the singularities below are taken into accounts and we obtain another temporal function g(t).The inverse Fourier transform is a special case of the inverse Laplace transform in which the horizontal line is the real frequency axis.Therefore, performing an inverse Fourier transform on the harmonic-domain function of an unstable signal does not allow the retrieval of the temporal signal (from (b) to (c) in Fig. 6).It is therefore necessary to perform, in general, an inverse Laplace transform over a horizontal line above the real axis and any potentially unstable pole.
Let us now calculate the temporal-domain function h(t) associated with the harmonic-domain function H(ω).We apply MOSEM to H (see Eq. 8) and derive the inverse Laplace transform of every term: A. Response to a Sinusoidal Input In many problems, the main focus is the global temporal response of the LTIS at a specific frequency.In these situations, the characterization of the system in the harmonic domain is sufficient since it is almost equivalent to a study of the temporal permanent-regime.We expect the response of a stable system to an excitation to be composed of that same input, scaled and phase-shifted in the permanent-regime, as well as exponentially decaying functions in the transient-regime [39].It is worth demonstrating this known result by deriving the response y(t) to a sinusoidal input x(t) = Re[e −iω0t u(t − t e )] of an LTIS of impulse response h(t).Since y and x are real, h must be real, which leads to: The response to a sinusoidal input is the real part of the response to a causal plane wave t → e −iω0t u(t).In the harmonic domain, Y and X can be written as: Therefore, the set of isolated singularities of the harmonic response includes the singularities of the transfer function with the same order, as well as ω 0 of order 1 introduced by the harmonic input signal X.Applying MOSEM to Y gives: with  24) to which a pole ip I with residue 8 + 0.1i was added.(b) The harmonic-domain function H(ω), which can be obtained by Fourier or Laplace transform of h(t).(c) A stable temporal function g(t) obtained via inverse Fourier transform of H.The Laplace and Fourier transformation are equivalent since the signal is causal, i.e. h(t) = h(t)u(t).The Region of convergence (ROC) corresponds to the part of the complex frequency plane where the usual Laplace transform is defined, and within which the inverse Laplace transform should be performed to retrieve the signal h(t) in (a).Below the ROC, another function would be obtained.In particular, the inverse Fourier transform is performed on the real axis, and leads to the stable function g(t) in (c) The calculation of the Laurent series coefficients of Y leads to: for the singularity p of order ν p .In addition, the residue of Y at ω 0 is proportional to the transfer function evaluated at the plane wave frequency ω 0 : By replacing the corresponding terms in Eqs.(40,41), we obtain: It follows that the inverse Laplace transform of Y can be written as the sum of a Dirac's delta function multiplied by Y NR , an expansion on the dynamic states similarly to the resonant term in Eq. ( 37) after replacing α(H, p, n) by α(Y, p, n), and the scaled and phase-shifted sinusoidal function t → e −iω0t : with The initial value theorem imposes that Y NR must be null if x is a continuous input.Since lim t→0 + y(t) = 0, i.e. the system cannot immediately respond to a physical signal, the initial value theorem applied to y gives: It follows that the constant, non-resonant term must be equal to 0, and that the sum of the residues of Y is null: We obtain the aforementioned expected result: the response y of a stable system to a causal sinusoidal input is only composed of the scaled and phase-shifted input in the permanent regime.In the transient regime, y must be expanded on the set of dynamic states f p,n of the system: This expansion remains valid in the case of an unstable system, although the concept of transient and permanent regime would no longer hold.In this case, the dynamic states would hold information regarding the diverging speed of y.

VI. CONCLUSION
We extended in this work the singularity expansion method to the general case of multiple order singularities in the complex frequency plane.Starting from simple considerations regarding the physical nature of the signals, we detailed the derivation of a more general singularity expansion, from which we deduced the singularity and zero factorization of a function.We calculated the exact temporal-domain expression of the response or the impulse response using the inverse Laplace transform of the generalized singularity expansion.By considering the case of the response to a sinusoidal input, we show that the singularities form a natural basis for the expansion of temporal responses in the transient regime, while they only influence the amplitude of the sinusoidal signal in the permanent regime.Finally, we inferred the constraints put on the poles and zeros in the complex plane for physically realistic signals possessing a Hermitian symmetry in the harmonic domain.Furthermore, causality was assumed to discriminate between stable and unstable poles based on the sign of their imaginary part.
We believe that the multiple-order singularity expansion method will find applications as a means to unveil specific properties of linear systems and their response by linking physical phenomenon to the distribution of the singularities and zeros, but also as a numerical tool to obtain highly accurate approximations of functions knowing only some of their singularities.Further works will explore the richness of these properties, in both the harmonic and time domains.
We wish to derive a global expansion of a complex function in terms of its singularities.The emphasis is put on the harmonic-domain function H associated with the temporal signal h.In what follows, H is assumed to be meromorphic, i.e. holomorphic everywhere on C except for a set of points P which is the set of its isolated singularities.Let us consider a domain of C containing ω, but no singularity of H, bounded by a curve γ C .By definition, H is holomorphic in this domain, and the Cauchy Integral Theorem gives: It is therefore possible to obtain H at the frequency ω from the knowledge of H around this frequency, in a small domain.As the size of the domain increases, more and more singularities are contained within it, and we have to make use of the residue theorem to fully characterize H at ω.
Let us now consider another closed curve γ (black curve in figure S1) around a domain D containing a subset of the singularities of H, and ω ∈ D. The integral over γ can be expressed as the integral over γ P , the closed curve surrounding the singularities and excluding ω (red curve in S1), plus the integral over γ C , a closed curve around ω containing no singularities (green curve in figure S1).
Therefore, the generalized version of the Cauchy's integral formula [2] leads to: The integral over γ C is 2iπH(ω) as stated in Eq.S4.The integral over γ P can be calculated using the residue theorem for the function R defined as: R is holomorphic everywhere on C except on its isolated singularities: shown in green in figure S1.The following expression is obtained: The residue of R for the singularity p can be expressed as a combination of the coefficients of negative order in the Laurent expansion of the function H around p, and the coefficients of the Taylor or Laurent series expansion of D : z → 1 z−w .The Taylor series coefficients of D can be explicitly calculated by definition, leading to an expression of Res(R, p) where only the Laurent series coefficients of H are unknown.The Laurent series expansion of the transfer function H for z in a small punctured disk around p is written as [3]: where ν p is the order of the singularity p and α(H, p, ) is the Laurent series coefficient of H of order for the singularity p.Similarly, the expression of the Taylor expansion of D : z → 1 z−ω is: β(D, p, m) being the Laurent series coefficient of D of order m around the complex frequency p (which is not a pole of D).R is meromorphic and its singularities, other than ω, are the same as H with the same order.Therefore, its Laurent series expansion near the singularity p is: We show the accuracy of MOSEM with M = 5 in Fig. S2, where we compare it to the expression of H from Eq. (S32).

III. SINGULARITY AND ZERO FACTORIZATION IN THE HARMONIC DOMAIN A. General Expression with the Singularity Expansion
We now use MOSEM to derive a factorized form of H using its zeros z m of order ν zm , and its singularities p of order ν p .If 0 is a zero of order m of H, i.e. ∀m ≤ m − 1, ∂ n ω [H](0) = 0, then the Hadamard factorization theorem tells us that H can be expressed in terms of another meromorphic function G which does not have 0 as a zero: The singularities and zeros (other than 0) of H and G, as well as their orders, are identical.Applying MOSEM (Eq.(S27)) to G yields: (S36) The complex log-derivative of G is the meromorphic function F defined as F = G G , with G = ∂ z (G) the first order derivative of G relative to the complex variable ω, which is also meromorphic.Using the Laurent and Taylor expansions of G and G , it can be shown that the singularities and of F are the zeros and singularities of H. Their p is a singularity of order 1 of F , and its associated residue is its order ν p .Similarly, for ω in the vicinity of a zero z of G:

FIG. 2 :
FIG. 2: (a) Modulus of the transfer function H(ω) defined in Eq. (2), as well as its resonant terms Res(p )/(ω − p ) associated with the poles p l and the modulus of the sum of these resonant terms.The frequencies range from 1.040 × 10 1 6 Hz to 1.093 × 10 1 6 Hz, which is equivalent to ultraviolet wavelengths between 172.9 nm and 181.3 nm.(b) Log-amplitude of H(ω) in the complex frequency plane, in a limited complex frequency window.The poles in (a) are highlighted in (b) as red points.The sum over the poles gives an approximation of the shape of H(ω) on the real axis (the minimum and maximum frequencies), but with the background term of Eq. (1) omitted, the reconstructed red curve poorly matches the exact expression.

FIG. 4 :
FIG. 4: Any physically realistic input X(ω), output Y (ω) or transfer function H(ω) can be described using either MOSEM or the SZF.The two expressions are equivalent.

FIG. 5 :
FIG. 5: Complex Bode diagram of the function H defined in Eq. (2).(a),(b) Log-amplitude of H for a complex frequency window in (a) and its symmetric window with respect to the imaginary axis in (b).(c),(d) The phase of H in the same frequency window (c) and its symmetric window (d).The amplitude is symmetric with respect to the imaginary axis and results in a symmetric distribution of the poles (red points).The phase is antisymmetric, as shown with the red and blue arrows indicating a clockwise 2π phase-shift around the poles with a positive and negative real part respectively.

FIG. 6 :
FIG.6:(a)An unstable temporal signal h(t) equivalent to e p I t for long times, with p I = 0.9.It is obtained via inverse Laplace transform of the function H(ω) from Eq. (24) to which a pole ip I with residue 8 + 0.1i was added.(b) The harmonic-domain function H(ω), which can be obtained by Fourier or Laplace transform of h(t).(c) A stable temporal function g(t) obtained via inverse Fourier transform of H.The Laplace and Fourier transformation are equivalent since the signal is causal, i.e. h(t) = h(t)u(t).The Region of convergence (ROC) corresponds to the part of the complex frequency plane where the usual Laplace transform is defined, and within which the inverse Laplace transform should be performed to retrieve the signal h(t) in (a).Below the ROC, another function would be obtained.In particular, the inverse Fourier transform is performed on the real axis, and leads to the stable function g(t) in (c) FIG.S1: Separation of the integration curve γ (black curve) into a circle γ C (red curve) around ω (blue cross) and a curve γ P (green curve) around the singularities.The zoom in is used to stress the complementary of the red and green curves to produce the black one.

∞
FIG. S2: Real part of the function H defined in Eq. (S32).The poles of H are the complex frequencies ω (n) p = 2nπ + π/2.The exact expression in black is compared to the truncated MOSEM in red: −M ≤ n ≤ M , with M = 5.

order is 1 ,F
and the associated residues are the orders of the zeros and singularities of G. Let us first consider a singularity p of G, of order ν p , and ω in the vicinity of p:G(ω) ≈ α(H, p, −ν p ) (ω − p) νp G (ω) ≈ −ν p α(H, p, −ν p ) (ω − p) νp+1