Graph theory approach to exceptional points in wave scattering

In this paper, we use graph theory to solve wave scattering problems in the discrete dipole approximation. As a key result of this work, in the presence of active scatterers, we present a systematic method to find arbitrary large-order zero eigenvalue exceptional points (EPs). This is achieved by solving a set of non-linear equations that we interpret, in a graph theory picture, as vanishing sums of scattering events. We then show how the total field of the system responds to parameter perturbations at the EP. Finally, we investigate the sensitivity of the power output to imaginary perturbation in the design frequency. This perturbation can be employed to trade sensitivity for a different dissipation balance of the system. The purpose of the results of this paper is manifold. On the one hand, we aim to shed light on the link between graph theory and wave scattering. On the other hand, the results of this paper find application in all those settings where zero eigenvalue EPs play a unique role like in coherent perfect absorption (CPA) structures.

In this paper, we use graph theory to solve wave scattering problems in the discrete dipole approximation. As a key result of this work, in the presence of active scatterers, we present a systematic method to find arbitrary large-order zero eigenvalue exceptional points (EPs). This is achieved by solving a set of non-linear equations that we interpret, in a graph theory picture, as vanishing sums of scattering events. We then show how the total field of the system responds to parameter perturbations at the EP. Finally, we investigate the sensitivity of the power output to imaginary perturbation in the design frequency. This perturbation can be employed to trade sensitivity for a different dissipation balance of the system. The purpose of the results of this paper is manifold. On the one hand, we aim to shed light on the link between graph theory and wave scattering. On the other hand, the results of this paper find application in all those settings where zero eigenvalue EPs play a unique role like in coherent perfect absorption (CPA) structures.
Although wave scattering is an elementary process and straightforward to picture, its analysis continues to fuel developments in electromagnetic and acoustic material research. While a small object (particle) scatters as a point source with a strength proportional to the applied field, larger objects scatter the wave between their constituent parts. This multiple scattering process is an infinite chain of possible scattering events, interfering to give the total field. This complicated interaction breaks the simple relationship between the applied and scattered wave amplitudes. From this complex interaction, several fields of research emerge including metamaterials [1], photonic crystals [2], propagation and imaging through disordered media [3], and random lasing [4].
The last decade has seen a large body of research into wave scattering in non-Hermitian materials, originating from Bender's proposed parity-time symmetric extension to quantum mechanics [5]. Non-Hermitian materials differ from ordinary matter in that they are usually driven, containing regions where the wave can be amplified, in addition to regions of absorption. This absorption and re-emission of wave energy provides much more control over the wave field compared to passive structures, demonstrated in designs for invisible and reflectionless media [6][7][8], cloaking [9], one-way propagation [10], coherent perfect absorption [11,12], and disordered media without scattering [13]. Although initially an obstacle, controlled wave amplification has now been demonstrated from GHz [14] to optical frequencies [15], as well as in acoustics [16][17][18].
In this work, we investigate the problem of designing non-Hermitian arrays of particles with controllable exceptional point degeneracies. Exceptional points (EPs) are peculiar to non-Hermitian materials where two or more modes of the system have both eigenvalues and * s.scali@exeter.ac.uk eigenvectors that coalesce. They have attracted considerable interest [19] exhibiting an apparently increased sensitivity to system perturbations [20,21], with the degenerate modes transforming into one another after cycling the system parameters [22,23]. To the best of our knowledge, while extensive work has been done on higher-order exceptional points [21,24,25], no consistent method to find N th -order EPs in wave scattering systems has been presented yet. In this work, we provide a recipe based on graph theory for implementing an exceptional point of arbitrary order in a system of scattering particles. The resulting system exhibits scattering properties with an extreme sensitivity to small changes in the particles' positions.
Our graph theory approach is based on the discrete dipole approximation (DDA) [26,27]. This is an established method for calculating the field scattered from any configuration of N particles. Originally introduced by Purcell to calculate the scattering from astrophysical dust [28], this method is now commonly applied to, e.g., metamaterial design [29,30] and wave propagation in disordered media [31] thanks to its vast range of validity [32]. By treating the particles as point sources, with a strength proportional to the incident field, the scattering problem can be solved self consistently determining the field on each particle. This requires the inversion of an N × N matrix, which rapidly becomes analytically intractable as the number of particles (scatterers) increases. Here, we provide a graph theory representation of this matrix inversion. We use this to understand the requirements on the scatterer parameters for the system to exhibit an exceptional point of arbitrary order, finding a remarkably simple picture in terms of vanishing sums of graphs related to different scattering events.
The paper is organized as follows: in Sec. I, we review the discrete dipole approximation (DDA). In Sec. II, we show how to interpret DDA by means of graph theory. In Sec. III, we derive the single scattering events and define orders of interactions. By means of the graph the- Figure 1. Schematic of a source field ϕs incident onto an array of sub-wavelength size scatterers (red dots) with polarizabilities αn. The scatterers respond to the source field, producing an outgoing field ϕout = n ϕ out n .
ory interpretation, we perform and give insights on weak and strong interaction limits. In Sec. IV, we present a method to design N th -order EPs with zero eigenvalue in systems described by DDA, perhaps the most important result of this paper. To do this, we derive the conditions to find these EPs (Sec. IV A) and, consequently, we interpret these conditions in terms of graphs in a scattering setting (Sec. IV B). In this setting, we show the effects of the EPs on the system's properties (Sec. IV C), namely the total field and the power output. Finally, we show how one can exploit perturbations to the design resonant frequency to tune the dissipation balance across the array of scatterers. However, this comes at the cost of a broader power output. In Sec. V, we conclude by summarizing the results and possible next developments.

I. DISCRETE DIPOLE APPROXIMATION
For simplicity, we restrict our theory to scalar waves of amplitude ϕ (e.g., the pressure of an acoustic wave in a fluid or, in two dimensions, the fundamental mode of a waveguide), although there is no obstacle to adapting our theory to vector waves. A model of the system presented in the following is shown in Fig. 1. We take N scattering particles of polarizability α n , with n = 1, 2, · · · , N . Subject to an incoming wave of amplitude ϕ inc , each of these particles will act as a point source s n of strength Note that the incoming field ϕ inc (x n ) is defined as the total field at position x n (position of the scatterer α n ) minus the self-field of the scatterer. The total field ϕ(x) obeys the three dimensional Helmholtz equation, including the sources of scattered waves given in Eq. (1), where k 0 = ω 0 /c is the wavenumber with ω 0 the resonant frequency, and s(x) is the externally driven source of waves in the system. Throughout this paper, we assume c = 1. The solution to the Helmholtz equation (2) can be written in terms of the 3D Green's function Integrating the Green function against the right hand side of Eq. (2) we have the solution to Eq. (2), which takes the form where ϕ s (x) is the integral of the Green's function over the source s(x). To determine the unknowns ϕ inc (x n ), Eq. (3) is evaluated on each of the N scatterers, excluding the infinite self-field, and demanding self-consistency, To write the problem in a more convenient form, we scale our field amplitudes by the polarizability, defining the new set of unknownsφ inc (x n ) = α n ϕ inc (x n ). Writing Eqs. (4) in matrix form, the solution is where the interaction matrix M is given by with the source field vector ϕ s = (ϕ s (x 1 ), ϕ s (x 2 ), · · · , ϕ s (x N )) T , and the incident field vectorφ inc = (φ inc (x 1 ),φ inc (x 2 ), · · · ,φ inc (x N )) T . Note that, in general, the matrix M is non-Hermitian, being both complex and symmetric. In non-reciprocal systems [33], the interaction matrix is both complex and asymmetric. From Eq. (5), we can therefore find a solution for the incident fieldsφ inc and consequently the total field ϕ(x) using Eq. (3). This is the discrete dipole approximation (DDA) method for solving scattering problems [26,28,34], reducing the entire problem to the matrix inversion M −1 . This must be done numerically even for a small number of scatterers [27].

II. GRAPH THEORY INTERPRETATION OF WAVE SCATTERING
Graph theory is a branch of mathematics rooted in Euler's solution to the problem of the seven bridges of Königsberg [35]. From here, graph theory stemmed and evolved, finding applications to many problems in science and engineering [36].
The interaction matrix M in Eq. (6) can be represented as a graph (e.g., in panel (a) of Fig. 2), where the diagonal elements (the particles' self-interaction 1/α i ) are represented as vertices, and their interaction (−G(x i , x j )) as edges. Multiple scattering events between the particles can thus be represented as a path on this graph, known as a Coates digraph. This representation links interactions and objects to edges and vertices respectively, fundamental constituents of any graph.
For example, take a 4-scatterer system whose matrix M 4 is the 4×4 equivalent of Eq. (6). In panel (a) of Fig. 2, we represent the matrix M 4 as the complete Coates digraph D * (M 4 ). Following convention [37,38], we refer to the Coates digraph using a star superscript. The Coates digraph is constructed as follows: the scatterers are represented by vertices, the Green's function interactions take the role of the edges, and the intrinsic (inverse) polarizabilities of the single scatterers are identified by the vertices' self-loops. This graph earns the technical name of vertex-labeled directed weighted simple graph permitting loops [39,40]. From now on, we will shorten and refer to this type of graphs as digraphs or simply graphs. This interpretation of the interaction matrix allows us to calculate the inversion of the matrix in Eq. (5) using graph theory. To do this, we consider the usual formula for the inversion of a matrix [41], where adj(M ) and det(M ) are the adjugate (transpose of the cofactor matrix) and the determinant of M , respectively. The i, j-th element of the adjugate matrix is defined as adj(M ) i,j = (−1) i+j det M (j,i) , where M (j,i) is the minor 1 built by removing the j th row and the i th column from the matrix M . Therefore, both terms on the right hand side of Eq. (7) depend on determinant evaluations. This form of inversion has a distinct interpretation in graph theory. It is thanks to this graph interpretation that we will be able to distinguish and identify different scattering events, ultimately solving for the total field of the system. In addition, using the same interpretation, we will illustrate a new visual way to build the condition to find zero eigenvalue EPs in scattering systems.
The determinant of a generic matrix A can be calculated using the Coates' determinant formula [37,38,42], where N is the number of vertices of the Coates digraph D * (A) and L is an element in the set L(A) of all the possible linear subdigraphs of the Coates digraph D * (A) [37]. A linear subdigraph of the Coates digraph D * (A) is a subdigraph of D * (A) in which exactly one edge enters and exactly one edge leaves each vertex [38,42]. The term γ(L) is the product of the weights of the edges of L, and c(L) is the number of cycles contained in L, i.e., the number of closed loops of the specific graph.
In panel (b) of Fig. (2), we show an example of a linear subdigraph L of the Coates digraph D * (M 4 ) (with N = 4). Following the just mentioned definition, note that exactly one edge enters and leaves each vertex. The number of cycles of this graph is c(L) = 2, while its weight is γ(L) = −α −1 1 G 2,3 G 3,4 G 4,2 . Following the same procedure applied in this example, we obtain the determinant of the matrix A by simply adding, according to Eq. (8), the appropriately-signed weights of the linear subdigraphs of D * (A).
Using a similar construction, the expression for the adjugate of a generic matrix A is [38], where the sum runs over all the possible 1-connections is obtained from a linear subdigraph (containing the edge j → i) by simply removing the edge j → i. Note that, in the case i = j, this corresponds to removing the self-loop at vertex i.
An example 1-connection is shown in panel (c) of Fig. 2. Starting by considering the linear subdigraph L in panel (b), we remove the edge j → i, that is, the self-loop 1 → 1. In this way, we obtain the corresponding 1-connection having number of cycles c(D * [1 → 1]) = 1 and weight γ(D * [1 → 1]) = −G 2,3 G 3,4 G 4,2 . Following the same procedure applied in this example, we obtain the adjugate element i, j of the matrix A by simply adding, according to Eq. (9), the appropriately-signed weights of the 1-connections of D * ([i → j]). See appendix A for further examples and more formal definitions of Coates digraphs, linear subdigraphs, and 1connections.
As a result, we can graphically represent Eqs. (9) and (8) for the matrix inversion (7), key for the evaluation of the total field of the system (3). These graph theory constructions, namely 1-connections and linear subdigraphs, give us a visual and systematic way of computing the elements of the inverse matrix M −1 . I.e., each element (M −1 ) i,j = adj(M ) i,j / det(M ) is evaluated by dividing the weighted sum of the 1-connections from vertex i to j by the weighted sum of the linear subdigraphs of D * (M ). As seen in section I, this inverse allows us to solve for the total field of the system (3). Although graph theory doesn't reduce the number of calculations required to perform this inversion, it provides a intuitive represen- tation of any scattering process in terms of a sequence of multiple scattering events. As we shall see, this allows us to give a graphical recipe for finding exceptional points in resonant scatterer arrays.

III. IDENTIFICATION OF DIFFERENT SCATTERING ORDERS
Before treating the problem of exceptional points in these scatterer arrays, we show how we can use Eqs. (8) and (9) for the construction of the elements of the inverse matrix M −1 in the case of weak and strong interaction limits of the system. These limits are taken by controlling the order of magnitude of the distance between the scatterers relative to the magnitude of the wavenumber used to probe the system. This results in a change of the interaction terms in the form of Green's functions G. To show how to evaluate these limits, we firstly demonstrate how 1-connections and linear subdigraphs capture all the possible interaction paths of the signal in the system. This allow us to identify scattering events of different orders to build approximations.
As a simple example, we consider a system of two scatterers characterized by polarizabilities α 1 and α 2 , symmetrically interacting via the Green's function G 1,2 . Now, we constructively build all the possible paths (or scattering events) of the system. To do this, we evaluate the incident field on the first scatterer, ϕ inc (x 1 ), while analogous considerations can be done for the second scatterer. The field ϕ inc (x 1 ) is the sum of all the possible paths starting from the different scatterers of the system and ending in scatterer 1. All these signals are scaled by the polarizability of the scatterer itself, α 1 . We start adding the contribution of a signal generated in scatterer 1, where the first term on the RHS is given by the source field. Proceeding in the same way, a signal propagating from the second scatterer is scaled by the polarizability of the scatterer itself, α 2 , then weighted by the interaction G 1,2 connecting the two scatterers, obtaining While these contributions account for the "one-round trips", the signals can propagate back and forth in the systems. Considering "multiple-round trips", we obtain 1 + α 1 α 2 G 2 1,2 + (α 1 α 2 G 2 1,2 ) 2 + · · · , where the term in the second square bracket accounts for the paths of different orders and extend to an infinite number of interactions. In the case of |α 1 α 2 G 2 1,2 | < 1, this last term can be written using the closed form of the geometric series as This is the analytical solution to Eq. (5) for the incident field ϕ inc (x 1 ) in the case of a symmetric 2-scatterer system. Note that, in Eq. (11), the terms in the numerator (i.e., the adjugate terms or 1-connections) represent the single scattering events, while the denominator (i.e., the determinant or linear subdigraphs) represent the possible multiple repetitions of the single scattering events. We identify the single scatter events and multiple repetitions by their order in the interaction G. For example, in Eq. (11), the numerator is made of 0 th and 1 st -order scattering events. In the same way, the denominator is made of 0 th -and 2 nd -order multiple repetitions. Proceeding in the same way for an arbitrary number of scatterers, we can build single scattering events and identify paths containing i th -order interactions. Now, we translate this interpretation of scattering events into the graph theory picture of Sec. II and we define different regimes of approximation. To do this, as a second example, we consider again the system described by and strong (panel (b)) coupling limits. By means of the construction shown above, in the case of weakly interacting scatterers, we restrict the sums in Eqs. (9) and (8) to those 1-connections/linear subdigraphs carrying weights γ up to second order in the interactions G (i.e., up to G 2 ), similar to the truncation of the Born series to second order [43]. With this approximation, we account for all those scattering processes whose graphs include no more  Unlike the Born series, which typically diverges in the limit of strong scattering, we can also take the limit of very strongly coupled particles, isolating those graphs with the largest number of edges (i.e., the highest nontrivial power of the inter-particle interaction G). Thus, we keep only the highest-order interaction terms of the sum in the adjugate terms and in the full determinant. In Fig. 3 panel (b), we see how these correspond to 1connections of order N − 1 for the adjugate and linear subdigraphs of order N for the determinant. Consequently, the most significant scattering event in the case of strongly interacting scatterers is represented by a signal traveling across the entire system and interacting with the highest number of scatterers 2 . Therefore, graph theory allows for a systematic way to calculate the total field ϕ(x) to any order in the interaction.
This graph interpretation results in a very efficient way of getting a good approximation of the total field ϕ(x) while only including the dominant scattering events in the weak (0 th , 1 st , 2 nd -order) and strong (N th , (N −1) thorder) cases. We show this in Fig. 4, where we evaluate the average percentage error of the absolute value of the approximated fields |ϕ(x) weak | (in panel (a)) and |ϕ(x) strong | (in panel (b)) against the absolute value of the corresponding non-approximated field |ϕ(x)|. The percentage error is averaged over 100 random values of scatterers' polarizations.

IV. N-TH ORDER EXCEPTIONAL POINTS
An exceptional point (EP) of a system is a non-Hermitian degeneracy in parameter space that emerges whenever two or more eigenvectors coalesce. The order of the EP is determined by the number of coalescing eigenvectors. At the EP, the matrix of the system is not diagonalizable but still admits a Jordan form [44]. In such form, the dimension of the Jordan blocks correspond to the order of the eigenvectors' coalescence, e.g., a 2×2 Jordan block corresponds to a 2 nd -order coalescence and so on. Finding these non-Hermitian singularities in small-dimensional systems is straightforward and an analytical solution can be quickly determined. Both 2 ndorder and limited higher-order EPs have been thoroughly studied [45][46][47] and experimentally realized [21,48,49]. However, no consistent method to find N th -order EPs in wave scattering systems has been presented yet. Note that we focus on those EPs with degenerate zero eigenvalue due to their clear physical implications on the total field of the system. In fact, since the total field depends on the inverse of the determinant, these eigenvalues are the cause to its highly degenerate responsiveness to parameter perturbation.
In the following, we use the transpose Frobenius companion matrix and its characteristic polynomial to explore N th -order zero eigenvalue EPs [25] and we interpret the result from a graph theory perspective. Note that, in a similar fashion, companion matrices and N -th order EPs have been recently studied in a tropical geometric framework [50]. We then design an EP in a scattering setting and probe the system's response against parameter perturbations.

A. EPs conditions
We now consider a system of N scatterers and impose the condition that, at some desired resonant frequency ω 0 , the interaction matrix (6) exhibits an N th -order EP whose eigenvalues coalesce to zero. As the outgoing field from the system depends on the inverse of the interaction matrix, this ought to yield a system whose power output diverges at the design frequency, and yet is also very sensitive to small perturbations (as in [21]), e.g., the scatterer positions.
We first consider the transpose Frobenius companion matrix M Frob associated with the matrix M of Eq. (6) [51]. The companion matrix is defined such that it generates the same polynomial for the eigenvalues λ of M , and is given by where the c i are the coefficients of the powers of λ in the characteristic polynomial, The form of the companion matrix is useful to us as it is closely related to the single N × N Jordan block matrix, The two matrices (12) and (14) take the same form once all the c i in (12) are zero. We assume that the interaction matrix M in Eq. (6) differs from (12) by a similarity transformation, an assumption which holds for the cases considered below. It is, in fact, sufficient for the interaction matrix M to have N distinct roots (in regime of no EPs) for the transformation M Frob = T −1 M T to exist [52]. The transformation matrix T = P Q −1 is derived as the product of the non-singular matrix P whose columns are the eigenvectors of M and Q whose columns are made of the set of N eigenvectors of M Frob , q i = (1, λ i , λ 2 i , · · · , λ N −1 i ) T relative to its eigenvalues λ i [53]. 3 For details on the derivation of such transformation T see App. B. With this assumption, there is an N th -order non-Hermitian degeneracy in the spectrum of M when all the c i are zero. By means of this simple requirement, we can engineer an zero eigenvalue EP of desired order by solving the set of non-linear equations given by the conditions c i = 0 for i = 0, 1, · · · , N − 1.
These coefficients c i can be evaluated relying on the expansion of the determinant in terms of its minors. Our system of equations for an N th -order EP with zero eigenvalue thus becomes where I m is the set of indices I m = {i 1 , i 2 , · · · , i m } defining the minor and S m ([n]) is the collection of size-m combinations within the set [n] = {1, 2, · · · , n}. Therefore, M (i,i) is the first minor obtained by removing the i-th row and column, M (i,i),(j,j) is the second minor obtained by removing i-th and j-th rows and columns, and so on.
Using this form to construct the coefficients c i , we numerically evaluate the solution to the non-linear system, identifying the parameters for an N th -order EP. Importantly, the EP conditions (15) are given in terms of sums of minors of the interaction matrix, which we have given a graph theoretic interpretation for in Eq. (8) and Eq. (9). For instance, satisfying the final condition in Eq. (15) requires a vanishing sum of the 1 × 1 minors, which equals the trace of the interaction matrix. From the identification shown in Fig. 2, this condition requires the vanishing sum of the self-interactions in the system. Thus, at a zero eigenvalue N th -order exceptional point, we require (among others) the condition that the inverse polarizabilities α −1 i sum to zero. Since the polarizabilities are complex, both the real and imaginary parts of the α will have to sum to zero, which is only possible in the presence of active scatterers, i.e., scatterers that exhibit gain. Moving up through the conditions (15), from c N −1 to c N −2 and so on, we see that all the second order interactions within the 2×2 minors must also sum to zero (equivalent to considering the 2×2 interaction matrix for every pair of particles in the system), as must the third order interactions defined within the 3 × 3 minors and so on. We thus reach the conclusion that an N th -order exceptional point can be associated with N conditions, each requiring the vanishing sum of sub-scattering events between a fixed number of particles. Note that the latter zero trace and determinant conditions found in the scattering matrix are reminiscent of the ones found in the case of systems described by a Hamiltonian with pseudochiral symmetry [25].
In addition to the maximal N th -order EP, we can also find n th -order singularities with n < N by requiring only the first n coefficients c 0 , c 1 , · · · , c n−1 to vanish. This generates a smaller non-linear system whose solution identifies an n th -order EP. This is only possible if n coefficients vanish in ascending order, starting from c 0 . In fact, this condition allows one to collect a factor λ n in the polynomial in Eq. (13), producing an n th -order λ = 0 solution. This solution corresponds to the n × n Jordan block relative to the n th -order EP. Any other combination of vanishing coefficients results in a diagonalizable system, without non-Hermitian singularities. Finally, note that the construction of EPs is inevitably dependent on the presence of interaction G in the system. In fact, in the case of no interaction, we would be left with a diagonalizable system.

B. Graph theory conditions for EPs
As an example, we now design a scattering configuration exhibiting a 4 th -order exceptional point and we interpret the condition of non-Hermitian degeneracy in terms of graphs. In the next subsection, we show how the scattered total field depends on a chosen parameter, in our case, the position of the first scatterer.
For the purpose of simplicity and readability, we now find the parameters (in our case, the polarizabilities α) that satisfy the EP conditions in a system in which the scatterers' positions are fixed. As sketched in Fig. 5, we equidistantly inscribe our scatterer array in a circle of radius r EP , simplifying the interaction matrix such that it contains only N/2 different Green's functions G when N is even, and (N − 1)/2 when N is odd. Given the limited number of Green's functions, this configuration is particularly convenient for an efficient search of the EPs. The interaction matrix associated with these cyclic polygons of scattering particles is, where G 1 represents the nearest-neighbor interactions, Figure 5. Example construction of a scattering system with an N th -order exceptional point. The system consists of N scatterers (here N = 4) with polarizabilities αn forming a cyclic polygon on a circle with radius rEP. Since the scatterers are equidistantly spaced, the angle θ is uniquely determined by the number of scatterers N , θ = 2π/N . The scatterers interact with the nearest neighbors via the Green's function G1 and with the next-to-nearest neighbors via G2. When probing the total field and the power output, we use the radial distance of the first scatterer r as the tunable parameter to scan through the exceptional point in parameter space.  , and power output (panel (c)) of the system in response to a change in the tuning parameter, that is, the radial distance of the first scatterer r. The latter ranges in r ∈ [rEP − ϵk −1 0 , rEP + ϵk −1 0 ], where ϵ defines a small deviation from the exceptional point. In panel (a), we show the absolute value of the total field normalized against the source field. In the scan from left to right (indicated by the white arrow), the tunable radial distance r is shifted by the amounts ϵ ∈ [−10 −3 , −10 −5 , +10 −3 ]. Note how the total field experiences a sudden peak in the proximity of the EP (middle plot). In panel (b), we show the global Euclidean distance of the right eigenvectors (see Eq. (17)) as a function of the tuning parameter r. This measure goes to zero when r = rEP. At this point, all the right eigenvectors (and corresponding left eigenvectors) merge into a single one. In panel (c), we show the power output (see Eq. (18)) with respect to the tunable parameter r for different purely imaginary shifts of the resonant frequency, Im{ω0} ∈ {0, 5 · 10 −4 , 1 · 10 −3 , 2 · 10 −3 }. For increasing imaginary shifts, the power response of the system broadens while the peak power at the EP reduces. Note that, given the general high gain of the system determined by the polarizabilities, the baseline power output remains of order 10 7 even for significant shifts from the ideal EP condition. G 2 represents the next-to-nearest-neighbor interactions, and so on. The angle between two consecutive scatterers is θ = 2π/N . In the figure, we also represent the tunable parameter, that is, the radial distance of the first scatterer r. While this parameter is not used to find the EP condition of Eq. (16) (it would indeed change the periodic-chain-like structure of the matrix in Eq. (16)), it will be needed later for the numerical analysis on the system's sensitivity to parameter perturbations. Our graph theory description previously introduced illustrates the meaning of this set of vanishing sums. For example, in Fig. 6, we show the condition c 1 = 0 which requires all the 3 rd -order scattering events to sum to zero. It is worth recalling that the zero condition of the i th -order coefficient is entirely independent of scattering events of any other order. This means that asking for the single coefficient c i to be zero is equivalent to asking for all the scattering events of order N − i to sum to zero. Thus, to find a 4 th -order EP, we need the condition c i = 0 to be satisfied by the scattering events of every order, that is, c i = 0 for i = 0, 1, 2, 3.

Our system is described by the 4 × 4 matrix M 4,sym
We finally note that, while a graph can be associated to the matrix of eigenvectors of the system, we could not find any particular interpretation to the coalescence of multiple eigenvectors in terms of graphs. Moreover, in the case of EPs of non-trivial order, a mathematical expression for the eigenvectors becomes highly cumbersome and strongly dependent on the system described. The non-trivial problem of finding a general expression for the eigenvectors of high-order EPs and an associated graph theoretic interpretation is left for further studies.

C. Trading sensitivity for dissipation balance
In Sec. IV B, we gave an example of a convenient system to find a 4 th -order EP. On this system, we interpreted the condition to find such EPs from a graph theory perspective. We now show how the presence of this high-order EP affects the total field of the system with respect to perturbations to the chosen parameter. In our case, this parameter is the position of the first scatterer r as depicted in Fig. 5.
In Fig. 7 panel (b), we show the coalescence of the eigenvectors in the range of parameter r ∈ [r EP − ϵk −1 0 , r EP + ϵk −1 0 ] with ϵ = 0.02 by means of the vanishing total Euclidean distance. This distance is defined as where v i and v j are the right eigenvectors of the matrix M 4,sym and the sum takes care of not double-counting terms. This quantity vanishes when r = r EP , signaling the coalescence of all the N eigenvectors relative to the degenerate eigenvalue 0. This is the N th -order exceptional point. Note that, given the high-order nature of the exceptional point, known EP measures like the phase rigidity of the eigenvectors and the condition number of the eigenvector matrix do not entirely capture the features of the singularity [54]. Note also that while the distance in Eq. 17 serves as an intuitive quantity to witness full eigenvector degeneracy, it is unable to give insight on the eigenvector scaling around the EPs. To do so, one can still access the phase rigidity's critical exponent [55,56]. The immediate effects of the EP on the total field are shown in Fig. 7 panel (a). In this figure, we scan, from left to right, through the EP with the tunable parameter r. In proximity of the EP, the absolute value of the total field |ϕ(x)| rapidly increases before attenuating again, once the singularity is passed.
In the same way, we can probe the EP just obtained by measuring the power output of our system of scatterers (see Fig. 5), which we define as Eq. (18) is derived, after little manipulation, by integrating the LHS of Eq. (2) (multiplied from the left by the complex conjugate field ϕ(x) * ) in a volume surrounding all the scatterers. We obtain the closed surface integral in Eq. (18) by means of the divergence theorem. The power output, as written in Eq. (18), depends on the sum of the incident fields on the different scatterers of the system weighted by the imaginary parts of the polarizabilities. In our case, the entire dependence of the power response on the tuning parameter r is contained in the incident field. This is uniquely determined by the matrix M sym . It is common to express the sensitivity (in our case, in the form of power output) of the system at the EPs in terms of a perturbation to the system matrix [21,57]. Thus, to express the power output in Puiseux series, one would need to rederive the scattering matrix in terms of a perturbation around the EP as for example M sym = J +εM ′ where J is the full Jordan matrix (14) and M ′ is a non trivial perturbation matrix [25]. Doing so, if the perturbation around the EP lifts the coefficient c N −1 such that c N −1 ̸ = 0, the Puiseux series λ = λ 0 + ∞ i=1 ε i/N λ i exists and refers to the N th -order EP. However, in case the perturbation leaves c N −1 = 0, the perturbed eigenvalues split in k different cycles of order n k of the form λ k = λ 0 + ∞ i=1 ε i/n k λ k,i with the various n k < N summing to N as k n k = N [25,45,58]. Note that, in the case of the system described in Eq. (16), a perturbation in the radial distance r indeed lifts the coefficient c N −1 such that c N −1 ̸ = 0.
Given the high order of the exceptional point, the power output of the system shows extreme sensitivity to perturbations in parameter space. In Fig. 7 panel (c), we show the power output of Eq. (18) versus the tunable parameter r for different imaginary offsets of the resonant frequency ω 0 at which the EP is found.
The introduction of an imaginary part in the design frequency has multiple functions. On the one hand, it helps to understand how possible experimental inaccuracies can affect peak and shape of the power output of the system. On the other hand, it shows how "ad-hoc" imaginary shifts in the design frequency of the system can help to adjust the distribution of gain/loss across the scatterers. Since the system is then probed with real frequencies ω 0 , introducing an imaginary shift in the design frequency results in a quasi-coalescence of the eigenvectors causing a drop in the system responsiveness to the singularity. This is shown in the figure by means of the amplitude reduction and broadening of the power output curves when increasing the imaginary shift of ω 0 . Note that the curve with Im{ω 0 } = 0 (solid blue curve in the figure), which is set to cross the EP, is re-scaled by a factor 10 −26 in order to fit into the graph and give some insight of the power output behavior.
We now show how we can tune the distribution of the gain/loss of the system across the scatterers in order to finely adjust possible experimental setups, where it is preferred to have a set of scatterers with the least possible gain. The presence of exceptional points inevitably depends on the scatterers' structure and, in particular, on their active nature. The condition for a scatterer j to be . This has been done in a similar fashion to Ref. [18]. The "cross" marker indicates Im{ω0} = 0 while the "left-caret" and "right-caret" indicate the end of the imaginary ranges, Im{ω0} = −1 and Im{ω0} = 1, respectively. The scattering elements are passive when the α lay on the positive semi-plane (red semi-plane), therefore they satisfy the inequality. In the figure, all the solutions α of the system considered are active, therefore laying on the negative semi-plane (green semi-plane). The figure shows how imaginary shifts in the resonant frequency ω0 used to design the EP allows one to tune the distribution of the gain/loss of the system across the scatterers.
passive is that its polarizability α j satisfies the inequality [59] Im{α j } > k 0 4π obtained by asking for a negative divergence of the power output in the case of passive scatterers. Eq. (19) is derived for the case of 3D Green's function as considered in this paper. As a reminder, k 0 = ω 0 /c with c = 1 in this paper. If the polarizabilty of a scatterer satisfies this inequality, the scatterer acts as a passive, lossy medium. As we have seen in Eq. (15), one requirement to obtain an N th -order EP is Tr(M ) = 0, i.e., all the scattering events of 0-th order have to sum to 0 while individually being non-vanishing. This implies having active elements in the system, i.e., scatterers with Im{α} < 0 which cannot satisfy Eq. (19). On the other hand, elements with Im{α} > 0 do not necessarily satisfy Eq. (19), thus, are not necessarily passive. By means of this inequality, we define a polarizability regime in which energy has to be injected into the system to obtain these N th -order EPs.
In Fig. 8, we show this inequality test for the polarizabilities of the system described by M 4,sym . In this case, none of the polarizabilities satisfy the inequality (no polarizabilities lie on the positive half of the plane), indicating that no passive scatterers are found in the system 4 . The test consists of a scan in the imaginary shift range Im{ω 0 } ∈ [−1, 1], where ω 0 is the resonant frequency at which the EP is evaluated. The "cross" marker indicates Im{ω 0 } = 0 while the "left-caret" and "right-caret" indicate the end of the imaginary ranges, Im{ω 0 } = −1 and Im{ω 0 } = 1, respectively. The semi-transparent lines represent all the intermediate α's solutions found in this range.
Note that, we already implemented this imaginaryshifted resonant frequency in order to control the spectral width of the scattering resonance of the system (see panel (c) of Fig. 7). However, in this case, one can use the imaginary shift to move the gain/loss bias on different scatterers. Therefore, an imaginary shift in the design resonant frequency allows one to fine tune the dissipation balance of the system in exchange of a broadening of the power output with respect to the EP parameter. This fine tuning capability becomes crucial in experimental setups that aim for the least possible gain in their set of scatterers.

V. CONCLUSION
In this paper, we used graph theory to solve wave scattering problems within the discrete dipole approximation (DDA).
Firstly, we showed how to use graph theory to develop a diagrammatic method for understanding multiple scattering processes. These processes are encoded in the inverse of the interaction matrix used to find the analytical total field of the system. We interpreted single scattering events in terms of 1-connections and linear subdigraphs and used these to approximate weakly and strongly coupled systems. This is a convenient machinery to calculate the total field ϕ(x) when the dimensionality of the system makes finding a full analytical solution impractical.
Secondly, by exploiting the Frobenius companion matrix associated with the system, we developed a systematic procedure to find N th -order zero eigenvalue exceptional points (EPs). The EPs are found by making vanish the sum of the 1-connections associated with scattering events of the same order. At a zero eigenvalue EP, the scattering becomes singular, causing the divergence of the emitted power. In our example, the perturbation coincided with a single-particle displacement from the EP configuration of the order of 1/100 of a wavelength. Although such a sharp sensitivity is achieved in position basis, one could describe the system in terms of the directions of input and output waves. Note that, as shown in this paper, one can also generate n th -order zero eigenvalue EPs where n < N . This might be useful to trade part of the scattered field sensitivity with a reduced number of conditions in the non-linear system. This further reduces the requirement for gain, crucial in certain experimental settings. The generation of N th -order zero eigenvalue EPs can be of particular interest for coherent perfect absorption (CPA) structures [60,61]. Here, the signature of the zero eigenvalue EPs (referred to as CPA EPs) is a quartic behavior of the absorption line shape in the perfectly absorbed channel. In addition, we believe the graph theoretical approach to be a promising tool to describe EPs associated with PT symmetry breaking in scattering systems [62,63] and the non-Hermitian skin effect in the case of non-reciprocal 1D chains of scatterers [64][65][66].
Finally, to control the spectral width of the exceptional points, we explored the consequences of displacing the design resonant frequency into the complex plane. We found that it is possible to trade the required gain/loss of the single scatterers with a broadened response. This would allow one to choose the preferred dissipation balance throughout the array of elements at the expenses of a reduction in the power output of the system. It might be possible to explore this trade-off as well as the entirety of multiple scattering physics in programmable metamaterials such as those demonstrated by Cho et al. [18].

SOFTWARE PACKAGE
The Julia package developed for solving the wave scattering problems found in this paper is available at https://github.com/mekise/graph-theory-dda. Note that, while the code should be easily readable for the user, it is not documented. Reasonable requests may be addressed to SS. ACKNOWLEDGMENTS SS thanks Federico Cerisola for stimulating discussions. JA and SARH thank the Royal Society for support. SS is supported by a DTP grant from EPSRC (EP/R513210/1). JA acknowledges funding from EP-SRC (EP/R045577/1). SARH acknowledges the Royal Society and TATA for financial support through the grant URF\R\211033.
The authors declare no conflicts of interest.
Appendix A: Graph theory fundamentals "A graph G is an ordered pair of disjoint sets (V, E), such that E is a subset of unordered pairs of V" [40]. The set V defines the vertices of the graph, i.e., the interacting elements of a structure we consider. The interactions between these elements are defined by the edges in the set E. In the case of interacting discrete scatterers, the set of vertices V represents the set of scatterers and the set of edges E correspond to the set of interactions between the Figure 9. Example of linear subdigraphs associated with the graph K3. For the example graph K3, there are 6 linear subdigraphs in total. Note that, given that we deal with directed graphs, we distinguish between subdigraphs with edges linking the same nodes but in opposite directions as in the case of the last and second-to-last subdigraphs in the figure. scatterers. Note that, in general, these interactions are not symmetric. By means of these fundamental blocks, we can translate every matrix M of the form Eq. (6) into a graph of the form 5. The resulting graph will represent the polarizabilities α as self-loops (or self-edges) and the Green's functions G(x i , x j ) as edges starting from the vertex i and ending in the vertex j. This directed edges, from i to j, promote the graph to a directed-graph or digraph. As mentioned in the main text, this graph is the Coates digraph D * (M ) associated with the matrix M . Note that the asterisk superscript takes care of the historical definition of the Coates digraph, i.e., the digraph associated with the transpose of the matrix we intend to represent [37,38]. In the main text, we refer to this kind of graphs as vertex-labeled directed weighted simple graph permitting loops. "Vertex-labeled", as the name suggests, indicates that the scatterers are distinguishable, "directed" means that interactions between scatterers are not necessarily symmetric, "weighted" indicates a non-unit interaction, "simple" indicates the presence of a single directional interaction between edges, while "permitting loops" identifies a graph that allows for selfinteraction, in our case, the polarizabilities.  A subdigraph is a digraph with V ′ ⊂ V vertices and E ′ ⊂ E edges. In addition, to earn the name of linear subdigraph, the vertices in V ′ must have in-degree and out-degree equal to 1, i.e., every vertex must have exactly one edge entering and one edge leaving. In Fig. 9, we report the entire set of linear subdigraphs for an example digraph K 3 . In the main text, we use these set of graphs to construct the determinants in the adjugate inversion formula. We now show how we use these linear subdigraph constructions for the determinant evaluation. Consider the sparse matrix G, We can work out the digraph associated with the matrix G and its linear subdigraphs to evaluate the determinant. To do this, we use Eq. (8), i.e., we search for all the subdigraphs whose vertices have in-degree and out-degree equal to 1. We show the results in Fig. 10, where on the LHS we find the digraph D * (G) associated with the matrix G and on the RHS we find the determinant of G, consisting of the only linear subdigraph of the graph D * (G). Summing the weights of the edges of the subdigraph, we obtain the determinant, det(G) = (−1) 4 (−1) 2 (1·3·2·2) = 12, where the first term Figure 11. Construction of the adjugate element adj(G)1,2, built using the off-diagonal 1-connections from vertex 1 to vertex 2. On the left, again the Coates digraph of the matrix G. On the right, the 1-connections of the matrix G which define the adjugate term adj(G)1,2 as per Eq. (9). The latter is built by the corresponding linear subdigraphs by removing the edge 2 → 1 as described in the text.
accounts for the factor (−1) N , the second term accounts for the number of cycles (−1) c(L) , and the last accounts for the weights of the subdigraphs γ(L).

1-connections
Consider the adjugate expression in Eq. (9), expression for the construction of adjugate terms by means of graphs. We report the expression here for convenience, (A3) In this expression, the terms D * [i → j] are the 1connections from vertex i to vertex j while all the other elements of the equation have an analogous meaning as in the determinant expression. The 1-connection D * [i → j] is obtained from the corresponding linear subdigraph L ∋i→j (linear subdigraph that includes the edge i → j) by simply removing the edge j → i. Note that, in the case i = j, this corresponds to removing the self-loop at vertex i. This definition leads to the following relation between the number of cycles in a linear subdigraph L ∋i→j and the relative 1-connection D * [i → j] [38], which justifies the "+1" in the adjugate expression. More formally, following the definition of a 1-connection reported in Ref. [38], we call 1-connection from vertex i to vertex j, the spanning subdigraph D * [i → j] such that, • for i ̸ = j, all vertices k with k ̸ = i, j must have in-degree and out-degree equal to 1, vertex i must have in-degree equal to 0 but out-degree equal to 1 and vertex j must have in-degree equal to 1 but out-degree equal to 0. The resulting spanning subdigraph therefore has a path from vertex i to vertex j, • for i = j, all vertices must have in-degree and outdegree equal to 1, while vertex i = j must have in-degree and out-degree equal to 0.
As mentioned in the main text, the 1-connections are closely related to the linear subdigraphs. In fact, the 1-connections D * [i → j] obtained using the definition above are equivalent to those obtained from the corresponding linear subdigraph L ∋i→j simply by removing the edge j → i. By means of this definition, we now show the construction of an off-diagonal adjugate term. Consider again the matrix G, we now build the adjugate term adj(G) 1,2 consisting of the 1-connections D * [1 → 2]. To do this, we consider all the linear subdigraphs that include the edge 1 → 2 (one single subdigraph in our example) and remove the edge from vertex 2 → 1, as shown in Fig. 11. Summing the weights of the edges of the 1-connections, we obtain the adjugate term, adj(G) 1,2 = (−1) 4 (−1) 1+1 (1 · 2 · 2) = 4, where the first term accounts for the factor (−1) N , the second term accounts for the number of cycles (−1) c(D * [i→j])+1 , and the last accounts for the weights of the 1-connections γ(D * [i → j]).

Appendix B: Similarity transformation between a matrix and its Frobenius companion form
Consider a matrix A ∈ C N ×N and its Frobenius companion matrix (see Eq. 12 in the main text), where c i are the coefficients of the characteristic polynomial of A. If there exists a row vector b ∈ C 1×N such that the matrix is non-singular, then the matrix A is similar to its Frobenius companion matrix A Frob [52], Here, the matrix T is the Vandermonde matrix, whose entries are thus given by where A i−1 :,j denotes the j-th column of A i−1 . Note that, in what follows, we do not compute the vector b but rather we infer the transformation matrix T while keeping b implicit. The Vandermonde determinant can be expressed as det(T ) = 1≤i<j≤N (bA j−1 − bA i−1 ). (B5) It is immediate to see that T is non-singular if and only if det(T ) ̸ = 0, thus the N rows b, bA, · · · , bA N −2 , bA N −1 are distinct. The rows of the Vandermonde matrix are generated by the different powers of A multiplied by the same vector b. Therefore, asking for T to be non-singular is equivalent to ask that the characteristic polynomial of A has N distinct roots, which requires A to be diagonalizable. Thus, we can write A = P DP −1 , where D is a diagonal matrix whose entries are the eigenvalues of A, and P is a non-singular matrix whose columns are the eigenvectors of A. Now, if A has N distinct roots, the Frobenius companion matrix can be diagonalized by the matrix Q whose columns are made of the set of N eigenvectors of A Frob [53] q i = (1, λ i , λ 2 i , · · · , λ N −1 i ) T (B6) relative to its eigenvalues λ i . Note that a consistent order of the eigenvalues must be kept throughout the calculations. We have that