Designing non-Hermitian real spectra through electrostatics

Non-hermiticity presents a vast newly opened territory that harbors new physics and applications such as lasing and sensing. However, only non-Hermitian systems with real eigenenergies are stable, and great efforts have been devoted in designing them through enforcing parity-time (PT) symmetry. In this work, we exploit a lesser-known dynamical mechanism for enforcing real-spectra, and develop a comprehensive and versatile approach for designing new classes of parent Hamiltonians with real spectra. Our design approach is based on a novel electrostatics analogy for modified non-Hermitian bulk-boundary correspondence, where electrostatic charge corresponds to density of states and electric fields correspond to complex spectral flow. As such, Hamiltonians of any desired spectra and state localization profile can be reverse-engineered, particularly those without any guiding symmetry principles. By recasting the diagonalization of non-Hermitian Hamiltonians as a Poisson boundary value problem, our electrostatics analogy also transcends the gain/loss-induced compounding of floating-point errors in traditional numerical methods, thereby allowing access to far larger system sizes.


I. INTRODUCTION
Condensed matter physics has traditionally been studied in the Hermitian context, since real energies are necessary for observing stable quantum states. Yet, with intense recent research in non-Hermitian systems [1], it has become apparent that many of the most exciting contemporary phenomena -exceptional points [2][3][4], non-Hermitian skin localization and modified bulk-boundary correspondence [5][6][7][8][9][10][11], nontrivial spectral topology [12][13][14][15][16][17][18], negative entanglement entropy [19,20], effective non-Hermitian curved spaces [21], amplified Rabi frequencies [22] -exist only in the non-Hermitian realm. Fortunately, non-Hermitian systems are not necessarily unstable, since they can still possess real eigenenergies if appropriately designed. To guarantee real spectra, a key approach has been to enforce parity-time (PT) symmetry [23][24][25][26], such that the gains and losses conspire to give rise to eigenstates with conserved total amplitude. However, it should be emphasized that having a PT symmetric Hamiltonian is neither a necessary nor a sufficient condition for obtaining real spectra -not sufficient because we also require the PT symmetry to be unbroken [27]. This additional condition restricts candidate non-Hermitian systems to particular optical * russell.yang@chch.ox.ac.uk † phylch@nus.edu.sg media [28][29][30][31] or lattice configurations [32,33], thereby forgoing potentially rich possibilities afforded by many other platforms [34][35][36][37]. On the other hand, PT symmetry is also not a necessary condition because pseudohermiticity provides a more definitive condition for real spectra [38,39]. However, pseudo-Hermiticity exists in our context only in the a priori unknown engineered real spectrum system, and not the parent PBC Hamiltonian to be found. All in all, this work is motivated by the lesserknown observation that even without PT symmetry, non-Hermitian gain/loss can still dynamically cancel in a bounded medium, if they enter through directed unbalanced couplings. In such cases, states moving in one direction will be amplified, while those moving in another direction will be attenuated. For a bounded system, these gain/loss processes cannot continue forever, and states become stable once they are "trapped" by a boundary. This mechanism thus constitutes an alternative approach to real energy spectra (pink region in Fig. 1a), in addition to PT symmetry and/or hermiticity. Indeed, such real spectra have been found in simple bounded non-Hermitian lattices with asymmetric nearest-neighbor couplings [40,41]. However, generalizing this dynamical mechanism for real spectra to richer, more complicated lattices remains a challenge, since unbalanced couplings also destroy the cherished Bloch property associated with ordinary translation-invariant states, necessitating the additional knowledge of the generalized Brillouin zone (GBZ) [8,42]. With the form of the GBZ being difficult to obtain analytically in all but the simplest cases [8,[42][43][44][45], the comprehensive design of systems with real non-Hermitian spectra thus hinges on a fundamentally more illuminating approach for treating the effects of unbalanced non-Hermitian couplings.
In this work, we circumvent these difficulties by developing a new electrostatics analogy that intuitively reconstructs parent Hamiltonians with desired spectra and spatial eigenstate profiles (Fig. 1b), different from known electrostatics analogies that related the correlations within complex non-Hermitian spectra with a Coulomb gas. This contrasts with conventional approaches where the Hamiltonian is fixed through a combination of symmetry arguments and empirical optimization, and the spectrum (desired or otherwise) subsequently computed from it. Being a reversal of conventional approaches, our electrostatics analogy also avoids directly diagonalizing non-Hermitian Hamiltonians, which usually suffers from fatal floating-point errors in large systems [42].
Our electrostatics analogy exploits the hitherto unexploited parallel between the conformal structure of electric field lines in real space, and the complex spectral flow as non-Hermitian states accumulate along boundaries ( Fig. 1 and Table I). It identifies non-Hermitian periodic/open boundary condition (PBC/OBC) spectra respectively with the loci of grounded conductors and electrostatic charges. Their electrical potential distribution gives the corresponding extent of boundary state accumulation, a key property that stabilizes state evolution, with no Hermitian analog. As such, instead of having to solve complicated algebraic equations to determine the GBZ, one just needs to solve the equivalent electrostatics problem, i.e. the Poisson equation [46], which is more geometrically intuitive. In particular, given a set of available couplings, all possible non-Hermitian models with real OBC spectra can be obtained by solving the corresponding electrostatics problem with charges restricted to a line, and grounded conductors determined by the coupling constraints.
In the rest of this work, we first explain the mechanism of obtaining real spectra from unbalanced (asymmetric) couplings, and use that to develop our electrostatics analogy and design approach. Next, after two pedagogical warm-up examples with well-known models, we demonstrate how we can robustly discover new models with real spectra, even when unexpected from symmetry arguments. Finally, we discuss our electrostatics reconstruction beyond the context of real spectra, where it is equally mathematically valid, and still physically relevant either as an optimization avenue, or in the context of network Laplacians.   1. (a) Hermiticity and PT-symmetry are two wellknown routes towards real energy spectra, but our work provides a new approach for designing generic real-spectrum Hamiltonians (pink) satisfying neither condition. (b) Conventionally, model Hamiltonian parameters have to be repeatedly optimized to yield the desired spectrum. By contrast, our electrostatics design approach directly outputs a parent Hamiltonian H possessing almost any desired eigenspectrum and eigenstate profiles. (c) Non-Hermitian eigenstates are characterized by their complex energies E and inverse spatial decay lengths κ, which together describe a landscape κ(E). In particular, PBC eigenstates lie along κ = 0 loops, while OBC eigenstates accumulate along ridges where κ(E) is not smooth. (d) The κ(E) landscape of a non-Hermitian system is mathematically equivalent to the V (x, y) potential landscape of its electrostatics analog, with PBC and OBC spectral loci corresponding to grounded conductors (V = 0) and lines of induced charges respectively. See Tables I, Sect. II A  and the Supplementary Sections I, II, III, and IV for more  elaboration. II. RESULTS

A. Unbalanced couplings, real spectra and their electrostatics analogy
Before deriving our electrostatics analogy, we first recap how unbalanced non-Hermitian couplings can lead to real spectra and hence stable eigenstates. Non-Hermitian systems that possess real spectra experience nontrivial gain and loss, but their effects balance out such that states do not grow or decay indefinitely, whilst retaining novel non-Hermitian properties. Unlike the conventional route of attaining this balance through certain global symmetries, i.e. PT symmetry, we shall consider a less known dynamical balance approach via unbalanced couplings between lattice sites in a bounded lattice. Away from boundaries (more generally, spatial inhomogeneities like disorder or impurities), unbalanced couplings cause states to grow or decay unimpeded while moving in different directions, resulting in eigenenergies with negative/positive imaginary parts. However, since amplification/decay is tied to spatial motion, further amplifi-cation of wavepackets is inevitably resisted upon reaching a boundary, resulting in boundary-accumulated skin states. If this resistance to amplification is perfect, the spectrum will become real.
The simplest example is given by the Hatano-Nelson model [40] with unequal nearest neighbor couplings |t + | = |t − | only. Without boundaries(PBCs), amplification is unimpeded and its PBC spectrum is complex, given by E HN = t + e −ik + t − e ik , k ∈ [0, 2π). However, in the presence of a boundary, its OBC spectrumĒ HN = 2 √ t + t − cos k is real [43] due to non-Bloch boundary accumulation While it is straightforward to understand the perfect balance of gain/loss in this particular simple model, boundaries induce a very complicated deformation k → p(k) = k + iκ(k) in most other models i.e. the OBC spectrumĒ converges towards a complex momentum deformation of the PBC spectrum E(k) i.e.Ē(k) = E(k + iκ(k)) [6], as elaborated in the Supplementary Section I. As a complex extension of the Bloch momentum, κ(k) corresponds to real-space decay, and is known as the inverse skin depth of the k-eigenstate. Mathematically, κ(k) must satisfy the condition that there exist at least two different k, k such thatĒ(k) =Ē(k ) and κ(k) = κ(k ). Such conditions can be expressed as a rather complicated singular function, which is in general difficult to solve.
Alternatively, our electrostatics approach sheds geometric intuition to the rather opaque algebraic problem of finding non-Hermitian OBC spectra, and provides a more intuitive and tractable approach for engineering the most generic lattice Hamiltonians with real spectra. To understand how, note that the complex deformation k → p(k) can be understood geometrically [ Fig. 1c]: While the PBC spectrum (solid loop) traces out a loop E(k), k ∈ [0, 2π) in the complex E plane, the OBC spectrum (line of crests within the loop) is obtained by the ramping up |Im(p)| such that the PBC loop "shrinks" into its interior until it overlaps with itself i.e is degenerate everywhere [6,15,47,48]. In this OBC limit, the spectrumĒ traces out lines or curve segments connected to each other at branch points. In particular, the OBC spectrum is real if the PBC loop successfully shrinks into a segment on the real line. Our electrostatics approach below shows that this can be achieved/engineered no matter how complicated the PBC Hamiltonian is.
To connect with electrostatics ( Fig. 1d), we turn to the conformal mapping p(E), which is the inverse of the complex energy dispersion, the momentum p and energy E both regarded as complex variables. The analog of p(E) in electrostatics is the quantity φ(z), where z = x + iy represents real space, and curves of constant U = Re(φ) and V = Im(φ) represent field lines and equipotentials respectively. Corresponding to them are curves of con-stant Re(p) = k and Im(p) = κ respectively (Table I). Figure 1 compares the profile of the inverse skin depth κ(E) with its analog, the electrical potential V (z). Intuitively, we make this correspondence since any conformal map, by definition, must preserve the orthogonal nature of the field lines and equipotentials. Likewise, unique isolines in the skin depth Im(p) = κ remain orthogonal to lines of constant Re(p) = k under adiabatic OBC-PBC transformation in steps of constant κ [43]. It should be noted that this electrostatics analogy relating the non-Hermitian skin effect problem to electrostatics is unrelated to other electrostatic analogies concerning the spectra of random non-Hermitian matrices [49][50][51], where the Green's function is found to have a 1/r dependence resembling the Coulomb electrostatic potential. In these works, the eigenvalue correlations are found to resemble that of a Coulomb gas, while in our analogy, it is the complex spectral flow that is found to contain the conformal structure of an equipotential vs. electric field lines pair.
In addition, it is quite straightforward to associate grounded conductors, i.e. contours of constant V = Im(φ) = 0 with PBC eigenenergy loops E(k), k ∈ [0, 2π) (solid closed loops in Fig. 1c,d). This is because Im(p) = κ = 0 for PBC Bloch eigenstates, due to their translation invariance. However, the OBC eigenenergies E(k) = E(k + iκ(k)) do not in general correspond to any equipotential, unless κ(k) happens to be constant. Instead, we identify the OBC spectra with regions of nonzero electric charge, since they are precisely where solution surfaces κ(E) intersect (crests within the loops in Fig. 1c,d), such that the corresponding potential surfaces V (z) become piecewise continuous. Specifically, nonzero electrical charge density ε 0 ∇ 2 z V is inherited from the nonzero ∇ 2 E κ = 0 at these crests. What makes this electrostatics analogy particularly useful is that the density of states (DOS) in the non-Hermitian system corresponds to the charge density in the electrostatic problem. As such, an electrostatic solution (of the Poisson equation) can be used as a proxy for solving the equivalent non-Hermitian problem, whose numerical diagonalization may be undesirably sensitive to noise. To understand this duality between charge density and spectral density, note that the DOS along an arbitrary curve in the complex E plane of a lattice with L sites is given by since there are L eigenstates along a loop E(p) parametrized by Re(p) = k ∈ [0, 2π), with fixed Im(p) chosen such that the loop intersects . Here E and and the unit tangent vectorˆ = /| | are treated as 2D vectors. The second equality was obtained via the Cauchy-Riemann relations. According to our electrostatics analogy, κ = Im(p) correspond to the electrical potential V , and ∇ E (κ) thus corresponds to the negative electric field strength −E = ∇V . Along a path such that there is Analogues Non-Hermitian lattice Electrostatic system Complex energy E Complex position z = x + iy Inverse skin depth/ imaginary flux κ(E) Electrical potential V (z) Density of states ρ Electrostatic charge density σ OBC to PBC spectral interpolation curve Electric field lines of const.
Path of nonzero charge density TABLE I. Correspondence between a non-Hermitian lattice and its electrostatic analogy. Of primary importance are the quantities in blue, which completely define their respective systems. Analogous to the electrical potential V is the skin depth κ −1 of accumulated states arising from the unbalanced non-Hermitian hoppings. Whether the OBC spectrum E(k) = E(k + iκ) can be real depends on whether the Poisson equation on V is consistent with nonzero charge density σ only along y = 0, i.e. nonvanishing spectral density of states (DOS) ρ along Im(E) = 0. Here refers to a generic path on the complex energy (spatial) plane E (z), emphasizing that σ and ρ are in principle proportional along any curve, not just the real line.
a discontinuity in the field strength, the induced charge density can be shown via a Gaussian pillbox to be where ε 0 is the permittivity of the medium, and the sign of ∓ dependent on the direction of the discontinuity step. Evidently, then, the charge and corresponding spectral densities are proportional, σ ∝ ρ , with the proportionality constant fixed by the constraint that the total number of states in the entire spectrum is L.
Some comments are in order regarding the discreteness of the charges. So far, in formulating the electrostatics analogy, we have tacitly assumed large L, such that the DOS and hence their corresponding charge density tends towards a continuum. However, there exist scenarios of non-Bloch band collapse [52,53], where there exists couplings towards only one uncompensated direction, and the OBC spectrum shrinks into one of more isolatedĒ points, with divergent κ and hence V . Such scenarios exactly correspond to electrostatic problems with isolated point charges, which include extremely well-understood textbook examples. As elaborated in Supplementary Section V, classic approaches such as the method of Images and superposition of point charges allow for the elegant design of Hamiltonians exhibiting non-Bloch band collapse onto an arbitrary variety of point charge configurations.

B. Designing Hamiltonians with desired spectra via electrostatics
We now elaborate on how our electrostatics analogy can help us engineer realistic Hamiltonians that possess desired OBC eigenenergies, particularly real eigenenergies. Conventionally, as indicated in Fig. 1b, given a lattice Hamiltonian H(k) with known properties, i.e. topological properties, one simply computes its OBC or PBC spectrum by numerical diagonalization under the respective boundary conditions. However, even though it is easy to write down a Hamiltonian H(k) that gives rise to a particular PBC spectrum E(k), it is much more nontrivial to design a H(k) that gives a desired OBC spectrumĒ(k) = E(k + iκ(k)), since κ(k) usually takes a very complicated form. Yet, it is the OBC spectrum that holds the key to achieving real spectra beyond the symmetry constraint paradigm.
Our electrostatics analogy is tailor-fit to solve this inverse problem of finding H(k) with desired OBC spectral properties. Based on Table I and surrounding discussions, the key input specifications for the desired OBC states are exactly the crucial data for constructing the electrostatic analog. Namely, they are (i) the locus of the desired OBC spectrum in the complex E plane, which specifies the locations of the analog electrostatic charges, (ii) the desired spatial profile (skin depth κ) of its eigenstates, which specifies the electrical potentials at their analog charges, and (iii) the shape of the desired PBC spectrum.
As illustrated in Fig. 1, the OBC eigenvalues/charges manifest as ridges, and (i) and (ii) respectively give the positions and heights of these ridges. To uniquely define the electrostatic system, we additionally specify the boundary equipotentials that enclose these ridges, which correspond to (iii) the shape of the desired PBC spectrum. Combined together, these data are sufficient for determining the electrostatic potential V (z) everywhere, and ultimately recovering the parent H(k) that produces them. This can be achieved through the following workflow, whose conceptual and implementation details are elaborated in the Supplementary Section III. First, we obtain V (z) everywhere by solving the Poisson equation with respect to Dirichlet boundary conditions from stipulated input data (i) to (iii). Next, by applying Gauss's law on E = −∇V , we obtain the induced charges on the boundary conductor, which corresponds to the PBC spectral locus. Due to the proportionality between electrostatic charge σ and its corresponding spectral density ρ (Eqs. 2 and 3), we hence obtain the DOS along the PBC curve stipulated by input (iii). As illustrated in Fig. 2a, this determines the k-point spacing along the PBC curve, leading to a full reconstruction of E(k) and hence H(k) as further described in the Supplementary Section III.
With a known candidate H(k), the real-space couplings can be obtained via Fourier transformation.
For more insight into the motivation behind our approach, it is instructive to contrast generic non-Hermitian scenarios with Hermitian scenarios. In the latter, the PBC and OBC spectra are both real, and largely overlap except possibly at isolated points. Hence their potential surfaces in Fig. 1c,d would have been compressed onto the real E line. This forces inputs (i) and (iii) to the same line segment on the real line. Furthermore, Hermitian eigenstates are Bloch states with κ = 0, so input (ii) would have been moot, with the electrostatic setup reducing to a trivial equipotential line, instead of a potential surface in the 2D plane. At a deeper level, our electrostatics approach reinforces that even if two non-Hermitian lattice systems possess the same OBC spectra, they can still differ in two fundamentally independent aspects: their (ii) eigenstate skin depths κ and (iii) PBC spectra.
Before presenting new Hamiltonians designed with our approach, we shall first provide two pedagogical illustrations using familiar non-Hermitian models-the Hatano-Nelson model [40] and the non-Hermitian SSH model [5,54]. Both models are already known to possess real OBC spectra due to unbalanced directed gain/loss, not symmetry protection, since the same models under PBCs do not have real spectra. However, for the purpose of illustration, we shall assume no a priori knowledge of these models, and derive their Hamiltonians based on the input data (i) to (iii) discussed earlier: to recap, (i) range of the desired real OBC spectra, (ii) desired skin depths κ −1 of their corresponding OBC states, and (iii) complex paths traced out by their PBC spectra, which uniquely specifies the possible output models.

Warm-up I: Hatano-Nelson model
In this simplest first example, we show how the HN model (Eq. 1) can be recovered just from stipulating its real OBC spectrum, uniform κ and elliptical PBC spectral locus. From the previous subsection, its OBC eigenenergies lie along the real segmentĒ HN ∈ , and the PBC eigenenergies E HN are distributed along an ellipse with semi-major/minor axes t + ± t − in the complex plane (Fig. 2b). Note that we have not parametrized both the OBC and PBC spectra, because the DOS and hence functional forms of these spectra are a priori unknown. We also stipulate that all OBC eigenstates have similar spatial decay profiles, i.e. uniform κ in anticipation of recovering the HN model; different distributions of κ(E) can lead to very different parent Hamiltonians, albeit with identical OBC and PBC spectral loci.
In Fig. 2b, these stipulated OBC and PBC spectra are mapped onto an electrostatic system with a real line segment (blue) of elevated uniform potential V = log t − /t + enclosed by an elliptical equipotential V = 0 (black) with semimajor/minor axes t + ± t − . The Laplace equation on this boundary-value problem possesses an analytic so- ]. This gives precisely the functional form of V = Im(p) in the (complexified) energy dispersion E HN (p) = z, establishing that energy E = E HN indeed corresponds to the position z in this analogy, and that H HN is indeed the parent Hamiltonian that yields these stipulated specifications on the spectrum and κ.

Warm-up II: non-Hermitian SSH model
We next show how a more complicated 2-component parent Hamiltonian can be recovered via our approach without the benefit of a known analytic solution. As a warm-up demonstration, the stipulated spectral loci are based on the known non-Hermitian SSH model [5,54] H SSH , even though the model itself is assumed to be a priori unknown, and meant to be reconstructed. In this case, the real spectrumĒ SSH is specified to lie within the two real segments 1 − T < |Re(E)| < 1 + T , corresponding to the segments of constant elevated potential V in the z plane of the electrostatic analog ( Fig. 2c). They are enclosed by a grounded (V = 0) outer boundary corresponding to the PBC locus of E SSH , which we stipulate (with the benefit of hindsight, for illustration purposes) as (4) γ and t = T 2 + γ 2 parameters of the H SSH to be found. Spatial coordinates (x, y) correspond to the spectral locus (Re(E SSH ), Im(E SSH )). We emphasize that specifying the locus of E SSH does not imply knowledge of either the PBC spectrum E SSH (k) or its Hamiltonian H SSH (k), since the parametrization with k is still unknown.
For a model with 2 or more components, obtaining a physically realistic Hamiltonian H(k) from the dispersion E(k) via Det[H(k) − E(k)I] = 0 presents an additional avenue of subtlety. As there are many possible forms of H(k) which all yield the same eigenenergy E(k), it is crucial to choose a correct ansatz that allows all components of H(k) to possesses rapidly decaying Fourier coefficients, such that the solution corresponds to a local Hamiltonian. In general, the Fourier coefficients f x of a function f (k) decay like f x ∼ x −(1+β) e −ηx , where η is the distance between the real line and the nearest complex singularity of f (k) and β, the order of that singularity, which is either negative or fractional. That is, with k 0 and f 0 being constants [55,56]. In particular, local Hamiltonians should possess η > 0, which means H(k) should not contain divergences or branch points at k ∈ R.
In this case, the reconstructed PBC dispersion can be fitted to the square-root expression E SSH (k) = π 2 0 π 2 π - Exact agreement between the stipulated and reconstructed spectra in our design approach, for the two warm-up models. (a) The spectrum E(k) of the output Hamiltonian is reconstructed by setting energy intervals ∆E between equal momentum spacings ∆k to be inversely proportional to the density of states ρ , which is obtained from the induced charge density σ in the equivalent electrostatic problem. (b) Warm-up example I with stipulated OBC spectrum on the real line segment |x| ≤ 2 √ t+t− with constant eigenstate decay profile κ −1 , corresponding to a constant electrical potential V . Its elliptical PBC spectral locus corresponds to a grounded conductor, giving rise to induced charges that enable the reconstruction of the full Hamiltonian HHN(k). As a check, its PBC and OBC spectra fall exactly on the initially stipulated loci. Parameters used are t+ = 0.9, t− = 1.8 and n = 42 lattice sites. (c) Warm-up example II with two stipulated real OBC line segments, PBC locus of the form of Eq. 4 and constant κ. The corresponding electrical potential induces charges that allow the reconstruction of an ESSH(k) dispersion, corresponding to the 2-component HSSH(k) from Eq. 5 with parameters γ = 0.4, t = 0.8 and n = 50 sites. Both the PBC and OBC spectra of HSSH display perfect agreement with their initially stipulated loci. Note that we have normalized the potential to V = 1 along the real line segments in both b and c.
± 1 + 2t cos k + T 2 + 2iγ sin k, which will result in unphysical long-ranged hoppings if we use a 1-component ansatz Hamiltonian. But since E 2 SSH (k) does not contain branch cuts, a prudent choice would be a 2-component ansatz Hamiltonian with only off-diagonal elements (see Supplementary Section III), which can be shown to be Shown in Fig. 2c is the excellent match between the OBC/PBC spectra of the reconstructed H SSH (orange/yellow), compared with the initially stipulated elevated potential interval along the real line (blue) and the boundary equipotential (black). This is one of the very few models where a constant κ deformation, i.e. potential V , can recover a real OBC spectrum-specifically, In the following, we shall explore how, by varying the PBC spectral locus and the κ(E) profile for OBC states, we can engineer a much greater variety of hitherto undiscovered models with real spectra.

Real OBC states with non-constant skin depths κ −1
We now turn to the first nontrivial demonstration of our approach, where we design a Hamiltonian stipulated to possess an OBC spectrum occupying two real line segments, each with a different value of κ. Consider the spectra shown in Fig. 3a, which naively looks like a simple variant of the Hatano-Nelson or SSH model, with specified OBC real line segments lying symmetrically about the origin, surrounded by an also symmetric elliptical PBC spectral locus. However, what is nontrivially asymmetric is the unequal inverse skin depth κ (corresponding to V ) of the eigenstates lying on each line segment, which cannot be fulfilled via any simple extension of these models. As a property with no Hermitian analog, κ is entrenched in the complex analytic properties of the Hamiltonian, and it is nontrivial to tune it without also modifying the spectrum. This is when our electrostatics analogy becomes valuable: corresponding to κ is the electrical potential V , and the task is reduced to that of solving the Laplace equation, subject to these V = 1 and V = 1.6 segments and the V = 0 outer boundary.
Upon solving this electrostatic problem and reconstructing the corresponding dispersion E(k), we find that E(k) has a slowly decaying train of Fourier components, and cannot result from a single-component local Hamiltonian. As discussed in the previous subsection, this is generically expected, and a resolving strategy is to consider an engineered Hamiltonian with 2 or more components (bands). To accommodate the asymmetry, we generalize the previous 2-component ansatz to: where σ ± = (σ x ± iσ y )/2, which allows the reconstructed energy dispersion E(k) to be fitted to either branch of the expression where T (k), H z (k) and F(k) each contain no more than a few Fourier components. The details of the fitting and Showcased in Fig. 3a is the very good agreement between the PBC and the OBC spectra of the engineered Hamiltonian, compared to their stipulated counterparts. For this example, up to next-nearest neighbor hoppings were kept in H (see Supplementary Section VI); even better agreement and reality of spectra can obtained with further hoppings. Notice the unequal eigenenergy spacings between the left and right halves, which result "naturally" as induced charges from the asymmetric electrical potential (pale orange-purple). This asymmetric example also illustrates the distinction between the PBC dispersion E(k) and the PBC spectral locus, which is the ellipse that is initially specified. Even though the latter is known a priori, the actual dispersion and density of states, which is skewed to the right as shown in Fig. 3a, can only be known after the potential profile V (x, y), i.e. κ(E), has been solved for. Graphically, it is also evident that steeper potential gradients lead to stronger induced charges and hence DOS on the right half. In the Supplementary Section VII, it is shown that good agreement was obtained not just for the stipulated and engineered spectra, but also for the κ(Re(E)) profile, which corresponds to the V (x) potential.

C. Real spectra without any symmetry
Even though our electrostatics approach can in principle generate a parent Hamiltonian corresponding to any desired spectral loci and skin depth profiles, the usefulness of the reconstructed Hamiltonian hinges on its experimental feasibility. As such, it is oftentimes preferable to restrict the reconstructed Hamiltonians to a small number of relatively local real-space hoppings. With further hoppings truncated, the compromise is that the reconstructed spectra may no longer be exactly real, even though the stipulated desired spectra lie along the real line. Even then, our approach allows for a systematic exploration of the region in parameter space where a class of local hopping models possesses almost real spectra, unaided by any specific symmetry.
As a demonstration, consider the rather demanding scenario where the desired PBC spectrum is stipulated along a so-called Reuleaux triangle (Fig. 3b), which is parametrically defined by the expression E Reuleaux (t) = √ 3e . Such a shape has 3-fold rotational symmetry about the complex origin, which generically do not encourage real OBC spectra. Meanwhile, the stipulated OBC spectral locus is set to be a segment of length d along the real line x ∈ (−1, 2), although it remains to be seen how faithfully the reconstructed OBC spectrum can reproduce this. As an additional complication, κ(Re(E)) (which corresponds to potential V (x)) is set to with offset x 0 = −1 for definiteness, and A, d parameters to be tuned. We proceed as in the previous examples, and is able to obtain local and simple Hamiltonians that produce almost exactly real OBC spectra while respecting the stipulated inputs of this problem. For instance, the spectra showcased in Fig. 3b arose from the engineered Hamiltonian (9) whose real spectrum possess eigenstate localization profiles closely obeying V (x), despite hoppings with no apparent symmetry whatsoever that even suggests of the possibility of a real spectrum. Containing only up to next-nearest neighbor hoppings, it is simple enough to realize in photonic systems given recent developments in the capability to control long range hoppings [57][58][59]. Mechanical and ultracold atomic systems also provide lee ways into implementing such systems [60,61]. But perhaps the most versatile means is mechanical or electrical systems where hoppings of arbitrary order can be generated with appropriate circuit configurations [35,62,63].
The propensity of real spectra can be increased by optimizing parameters A and d. Shown in Fig. 3c is the value of OBC Max|Im(E)| in the parameter space of A and d, with hoppings across up to 5 sites. In general, although the reconstructed OBC spectra are not always perfectly real, they are still mostly real with "branches" pointing away from the real line (See Fig. S3). Notably, there exists a large parameter region (purple) where the OBC spectrum is almost perfectly real, with Im(E) at least two orders of magnitude smaller than Re(E). Examples of almost-real spectra at different points in the parameter space can be found in the Supplementary Section VIII.
It is remarkable that this Hamiltonian can possess any real OBC spectrum at all, given that it not only lacks suitable conventional symmetry such as hermiticity or PT-symmetry, but also lacks any other ostensible symmetry that can encourage PBC→OBC spectral flow towards the real line.

D. Design of generic spectra from electrostatics
While most of this work is focused on the engineering of new classes of model Hamiltonians with real spectra, our approach is equally useful for obtaining those with desired non-real spectra of any shape, since the electrostatics analogy remains valid for PBC or OBC spectra, i.e. boundary potentials of any shape. We hereby describe two additional applications, where only certain branches of the OBC spectra need to be real.

Hamiltonians with real spectral branch
When desired states can be selectively filled, it suffices that the spectrum possesses a real branch which can Designing beyond real spectra. (a) Almost perfect agreement between the stipulated and engineered spectra in the right lobe was achieved, with only up to next-nearest neighbor hoppings, by discarding solutions from the two other extraneous lobes. The stipulated PBC spectrum is given by the parametrization Eq. 10. (b) Excellent reconstruction can also be achieved by enlarging the number of components of the ansatz Hamiltonian and its solution branches. Here, the 3-component ansatz H 6-fold (Eq. 13) with ∆ = 0 allows for an excellent local fitting of E(k) 3 with minimal Fourier components. Stipulated OBC and PBC spectra are given by Eqs. 11 and 12, with rmin = 0.9835, rmax = 2.9532. We used 30 and 120 lattice sites for (a) and (b) respectively, the large number in the latter necessary for demonstrating excellent convergence in all the branches.
be occupied, such that the reality of other branches become immaterial. Relaxing the requirement for the entire spectrum to be real greatly enlarges the space of candidate Hamiltonians. Shown in Fig. 4a, for instance, is a simple example where the stipulated real OBC branch x ∈ (0, 1.9) lies within a lobe-shaped stipulated PBC loop parametrized by Considered in isolation, the reconstructed Hamiltonian would have required many non-local terms due to the cusp in the PBC loop. However, if this PBC lobe and real OBC branch were part of the full 3-fold symmetric spectrum (faded in Fig. 4a), the solution to E(k) would be simply given by the expression of E(t) above.

Multi-component Hamiltonians
With this example in mind, we can repeat our methodology for more sophisticated multi-component models. Here, the entire stipulated OBC spectra is not restricted to the real line. For definiteness, suppose we desire to have both the OBC and PBC spectra exhibit approximately straight line segments, such that we can get various real spectra sectors upon suitable energy translation or rotation. In Fig. 4b, the stipulated OBC spectrum is given 6 straight line segments related by π/3 rotations, together with a zero mode, while the stipulated PBC spectrum is parametrized by which is a sixfold rotationally symmetric figure with approximately straight sides when a = 10. It should be noted that this t parametrization will not correspond to the expected PBC dispersion relation E(k) of the engineered Hamiltonian, since only the locus of points needs to be matched. From the geometric symmetry of the spectra, we chose a 3-component ansatz for the engineered Hamiltonian of the form with dispersion E(k) = ∆ + e 2πin/3 3 f (k), n = 0, 1, 2, such that f (k) can be found from the truncated Fourier transform of the reconstructed (E(k)−∆) 3 obtained from our induced electrostatic charge solutions. For Fig. 4b, we have set the energy shift ∆ to be zero; we could have set ∆ = 3i had we wanted to have a segment of approximately real PBC spectrum. As shown, keeping only hoppings within the same unit cell and across 3 unit cells in f (k) (see Supplementary Section VI), the reconstructed spectra agree with the stipulated spectra almost perfectly. This would not have been possible had we chosen another ansatz that involves some root of f (k) other than 3 f (k).

E. Discussion and conclusion
In a handful of memorable pages on electrostatic analogs and the unity of nature [64], Richard Feynman pointed out how many different phenomena in physics can all be explained using the equations of electrostatics. These include the heat flow between plates held at different temperatures, the vibrations of a drumhead, the diffusion of neutrons and the flow of a fluid past a sphere. Here we unravelled a new electrostatics analogy from non-Hermitian condensed matter physics and inverse quantum engineering. Specifically, we rigorously prove that the seemingly hard problem to synthesize non-Hermitian spectra in tight-binding lattices with modified bulk-boundary correspondence can be mapped onto a simple electrostatic problem. The working principle behind this analogy is the intimate relation between the conformal structure of electrostatics and the complex spectral flow in non-Hermitian systems. Our work opens up a new paradigm for engineering non-Hermitian spectra, particularly real spectra, in various settings, such as cold atoms [57,58], photonics [24,28], metamaterials [65,66], mechanical and acoustic systems [60]. While real spectra are important for state stability in the majority of experiments, we point out that non-real spectra present further possibilities in terms of topological sophistication [43,67], and are just as physically relevant in the form of the Laplacian spectra of steady-state networks such as electrical circuits [35,63,[68][69][70].
Our approach will impact the design of non-Hermitian sensors [71] in optics and electronics. As finite-sized devices, they are OBC systems, and the requisite sensing properties would specify the non-Hermitian skin state profiles, which correspond to the potential profile of the analog electrostatic charges. Such sensors may find applications in measuring glucose concentrations [72] and wireless sensors [73].
Besides its versatility in discovering new classes of stable non-Hermitian Hamiltonians, our electrostatics analogy also allows numerical access to non-Hermitian lattices of far larger system sizes. Till now, the numerical computation of OBC skin spectra have been limited by rapidly compounding floating point errors, which limit accurate diagonalization results to systems smaller than O(10 2 ) sites at standard machine precision. Analytic results only exists for a small subset of systems where the generalized Brillouin zone is not excessively complicated. However, our electrostatics approach trades this non-Hermitian diagonalization problem with a boundaryvalue partial differential equation problem, whose numerical solution do not suffer from non-Hermitian sensitivity at all, and can be extended even towards the thermodynamic limit.

F. Author Contributions
CHL developed the electrostatics analogy and initiated the project. CHL and RY led the project. JWT performed the spectral reconstruction. RY, JWT, TT and JMK performed the numerical computations. SL generalized the electrostatics analogy to point charges. SL, LL and CHL also took on advisory roles. All authors contributed to the writing of the manuscript.

G. Acknowledgements
This work is supported by Singapore's MOE Tier I grant WBS no. A-800022-00-00.

I. OBTAINING REAL OBC SPECTRA THROUGH THE DYNAMIC CANCELLATION OF GAIN AND LOSS
In the Hatano-Nelson model discussed in the main text, we have seen the simplest instance of real spectra originating not from symmetry constraints, but from the cancellation of unbalanced gain and loss as the state is pumped against a boundary. To rigorously understand this mechanism, we review the physics of state pumping due to unbalanced gain/loss, formally known as the non-Hermitian skin effect (NHSE) [5,6]. Consider a 1D system with L sites governed by a Hamiltonian H = k H(k)c † k c k . Its PBC eigenenergies E are simply given by the momentum eigenstate solutions of the characteristic polynomial (energy dispersion equation) Det[H(k) − E I] = 0, where k = 2nπ/L, n ∈ Z. Under OBCs and other scenarios with broken translation symmetry, however, eigensolutions must vanish at both boundaries, and must thus consist of a linear combination of at least two of such "momentum" eigenstates -which agree with the true momentum eigenstates only between the open boundaries, but vanish elsewhere. Yet for NHSE systems with unbalanced couplings, such degenerate eigensolutions for the same OBC eigenenergȳ E are spatially decaying, and are thus characterized by momenta p(k) = k + iκ(k) with nonvanishing imaginary parts which we shall call Im(p(k)) = κ(k). For their linear combination to exist in the thermodynamic limit, without one solution dominating the other, we thus require degeneracy also in the exponent of ψ k (x) ∼ e −κ(k)x , i.e. degeneracy of κ(k). Such OBC solutions are exponentially decaying away from the boundaries with localization lengths κ −1 (k), and are thus known as skin states. Specifically, we are concerned with engineering arbitrarily complicated H(k) such that the OBC eigenenergiesĒ ofH(k) = H(k + iκ(k)) are real.

II. OVERVIEW OF SPECTRAL RECONSTRUCTION WORKFLOW
In this section, we present a more detailed run-through of our workflow in engineering the parent Hamiltonians of desired spectra. The inputs and outputs of our engineering scheme are summarized in Table S I. As introduced in the main text, the inputs represent the desired features of the Hamiltonian to be engineered, namely the complete specification of the OBC spectrumĒ and eigenstate decay rate (skin depth κ −1 ), as well as the shape (locus) of the PBC spectrum E ∈ PBC in the complex energy plane. Note that specifying the PBC spectral locus does not entail specifying the PBC energy dispersion E(k) ∈ PBC -in the latter, the parametric dependence on the momentum k also contains information on the density of states |∇ k E| −1 , and allows the construction of possible parent Hamiltonian H(k) consistent with the dispersion. Likewise, the OBC Hamiltonian is also not completed specified with the input data, since its DOS is not given. The inverse skin depth κ(E) associated with each eigenenergy E on the OBC spectral locus E ∈ OBC represents the complexification of the Bloch momentum, and only exists in non-Hermitian systems with broken translation symmetry, i.e. OBCs.
Our engineering workflow is designed to leverage readily available tools in solving the Poisson equation in electrostatics. Based on the three desired aspects of the non-Hermitian lattice system that are a priori specified, we can unambiguously specify an analogous electrostatic system. As summarized in Table I of the main text, the PBC spectral locus PBC in the complex E plane corresponds to the shape of the grounded equipotential, i.e. conductor in the electrostatic analog. The OBC spectral data does not correspond to any equipotential in general, since OBC eigenstates are characterized by the inverse skin depth κ −1 (E), which correspond to non-constant potential V (z) on OBC . In essence, our input data of the desired non-Hermitian model features completely specifies the corresponding electrostatic problem by specifying all the boundary potentials (at PBC and OBC ).
This boundary value problem can be directly treated via a numerical Laplace equation solver. In the equivalent electrostatic problem, the effects of the charges are encapsulated in the boundary potentials V (z). Solving it, we recover the (output) potential V (z) everywhere. While the non-Hermitian analog to V (z) is κ(E), which is not a very useful quantity by itself, what is useful is the induced charge on the conductors σ , = PBC due to the output V (z). In particular, the spectral density of states ρ is proportional to electrostatic charge σ. Together with the shape of PBC on the complex E plane, we can reconstruct the complex energy dispersion E(k) and hence its parent Hamiltonian through techniques elaborated below. The profile of V (z) similarly allows us to recover the DOS of OBC states, but that cannot be directly used for reconstructing the parent Hamiltonian.

III. WORKFLOW DETAILS
Below, we do an explicit walk-through of our workflow, starting from specifying the problem to obtaining the final reconstructed Hamiltonian and its spectra. For illustration purposes, we shall base this demonstration on the non-Hermitian SSH model (warm-up example II in the main text), with stipulated input data generated from parameters γ = 0.4 and t = 0.8. Step Step 5: Step 4: (4.a) Compute Re(E ) and Im(E ) against implicit indexed variable k Without Trace With Trace Nonherm. lattice Electrostatic system Input OBC spectral loci OBC Location of given charges

Input
OBC skin depth κ −1 (E) on OBC Potential at given charges Input PBC spectral loci PBC Grounded conductor shape Output DOS of PBC states Ind. charge on conductor Output DOS of OBC states Values of given charges Output κ(E) profile everywhere Potential V (z) everywhere TABLE S I. Summary of input and output quantities of our spectral engineering workflow, together with their electrostatic analogs. Blue quantities are directly involved in recovering the engineered parent Hamiltonian. The inputs specified are desired features of engineered model, which include a complete specification of the OBC states and target PBC eigenenergies. The latter, combined with their corresponding density of states (DOS) which is output from our workflow, completely specifies the PBC and hence parent Hamiltonian. As a bonus, our solution recovers the potential V for the electrostatic analog, as well as its sources/sinks (charges) everywhere.
A. Defining the input by stipulating the desired spectra i.e. band structure/boundary conditions In the thermodynamic limit L → ∞, the PBC and OBC spectra trace out loci given respectively by a continuous PBC loop and an OBC graph. While an actual lattice possess only a finite L number of points, a continuum approximation is crucial in defining the electrostatic problem, such that the charge density σ , and by extension the density of states (DOS) ρ is properly defined everywhere.
By definition, the PBC loop corresponds to κ = V = 0, and is thus represented by a grounded conductor that defines the outer boundary of the electrostatic problem (outer loop in Step 1 of Fig. S1). The OBC spectral locus lies within the PBC loop, and generically corresponds to non-constant potential V (x, y) (schematic blue branching lines in Step 1), as stipulated by the OBC inverse skin depths, even though it is constant for the non-Hermitian SSH model. This OBC potential forms another set of boundary potentials.

B. Solving the Electrostatic Problem
After the electrostatic problem is completely specified via its boundary potentials, the Laplace equation is solved within the region between the boundaries (PBC and OBC spectral loci), via a finite-element method scheme. This is a unique Dirichelet problem with open curves inside a closed loop (OBC loci) [46], and there is no charge away from the OBC and PBC loci, since all charges only reside along the OBC loci, where the potential is allowed to be discontinuous. Importantly, as a fixed-size problem treated by efficient finite-element solvers, our approach is applicable for arbitrarily large systems, even in the thermodynamic limit, unlike the explicit computation of non-Hermitian spectra through diagonalization, which can run into convergence issues due to the NHSE even with only O(10 2 ) sites. After obtaining the potential V (x, y) everywhere along the interior of the PBC loop, one may then calculate the electric field E = −∇V . By Gauss's law, the charge density σ along the PBC boundary is given by wheren is the vector normal to the PBC loop (grounded conductor), pointing inwards toward the OBC loci. For the purpose of obtaining the DOS ρ ∝ σ along the PBC loop, we compute and plot σ along the arc length s of the PBC loop. In 2b of Fig. S1, we observe two peaks and troughs for σ over a period of L = 80 sites, consistent with the 2-fold rotational symmetry of the setup. In our implementation, this step is integrated with the following steps 3 and 4 into an automated routine. The engineering of Hamiltonians given desired specifications is thus made simple, without the need for manual intervention. Details on the code are available upon request.

C. Discretization of Points as specified by density of states/charges
The next step in obtaining the dispersion E(k) is to infer the rate of change of E on k along the PBC loop, from density of states (DOS) data. In our electrostatics analogy, this DOS ρ is conveniently proportional to the previously computed charge density σ .
We partition the PBC loop with L points, such that the separation between the points are inversely proportional to the DOS ρ = ∆k/∆E. In k-space, each point is thus separated by a fixed momentum interval ∆k = L/2π, and the set of all points therefore gives the PBC spectrum evaluated E(k) at L equally spaced k points. In step 3 of Fig. S1, regions of low/high DOS, i.e. charge correspond to points that are further/closer together. It should also be noted that it is irrelevant where the k = 0 reference point is placed, and simply starting at any point in the loop and picking a consistent traversal direction is sufficient.

D. Generating the momentum dispersion
Through simple interpolation, the functional forms of Re(E(k)) and Im(E(k)) are computed, from the numerical data at the L equally-spaced momentum points. If necessary, the higher powers of the reconstructed E can be similarly computed.

E. Reconstructing the Hamiltonian through Fourier series approximations
In this final step, we engineer a reasonably local Hamiltonian that yields the spectrum E(k) reconstructed from the DOS data. First, an ansatz Hamiltonian H(k) needs to be selected, which can be any n × n component (band) matrix. For practicality, however, we would like n to be as small as possible, subject to the constraint that the components of H(k) contain a small number of Fourier components. Higher-order Fourier components correspond to distant couplings, thus allowing too many of them ruins the locality of the Hamiltonian, and is experimentally undesirable.
While a few scenarios may already be satisfactorily solved with a 1-component model, usually at least 2 components are needed. This is to prevent unnecessary branching of the OBC locus caused by slow power-series convergence under a 1-component model. When the reconstructed dispersion E(k) cannot be expressed as the sum of a few Fourier terms, better convergence can often be achieved by taking higher powers of E(k) or its truncations. As such, by using a multi-component ansatz Hamiltonian, one arrives at matrix elements that are functionally dependent on the dispersion in ways that could result in far enhanced convergence. Two examples of 2-component ansatzes are given in step 5a of Fig. S1. Here, we discuss the simpler one with H(k) = ((0, F(k)), (1, 0)), where there is only one undetermined matrix element F(k) to be fitted, and F(k) = E(k) 2 . This demonstrative ansatz has been chosen for simplicity -in well-known examples like the SSH model, both off-diagonal elements are present, which usually makes for easier physical realization. For the particular example featured in Fig. S1, this ansatz with a single functional degree of freedom F proves sufficient, since the numerically fitted Fourier series (step 5b) contains only a few powers of e ik with non-negligible coefficients, and corresponds to the sum of a few local couplings.
After obtaining the engineered (output) Hamiltonian, one may verify (step 5c) that its PBC and OBC spectra indeed agrees well with the originally stipulated (input) spectra. One can then be confident that numerical and truncation errors accumulated in the solution are accept- FIG. S2. Example use of a trace term in drastically improving convergence to the stipulated real spectrum. (a) With a traceless ansatz H(k) = F(k)σ+ + σ− with F(k) given by Eq. 13, there is poor agreement between stipulated (blue) and engineered (orange) OBC spectra, particularly at the discontinuity centered at the origin. (b) Upon implementation of an optimized ansatz with trace term as detailed in Eqns. 6 and 13, the spectral flow from the PBC to purely real OBC line becomes what is desired, with the two solution branches morphing into a pair of purely real spectra segments (orange), accompanied by unwanted eigenenergy bands shown in green and red for PBCs and OBCs respectively. (c) After discarding the unwanted eigenenergies, we reach an almost exact agreement between the stipulated and reconstructed real spectra.
ably small. In Fig. S1, we observe that, for our illustrative system, the agreement is virtually perfect.

A. Exploiting additional degrees of freedom in constructing Hamiltonians
In many cases, the dispersion E(k) does not correspond to a rapidly decaying Fourier series, thereby precluding a local single-component parent Hamiltonian. An arbitrarily high number of Fourier coefficients may be required to achieve the desired accuracy.
However, a multi-component Hamiltonian with local matrix elements provide the requisite degrees of freedom that satisfies the same dispersion. As previously discussed, a simplest 2-component extension goes by the ansatz H(k) = F(k)σ + + σ − , which works when E(k) 2 = F(k) instead of E(k) is local. However, this form is still too restrictive, applicable to only a special subset of models.
By introducing a trace term T (k) = Tr(H(k)) = 0, we can expand the possible functional forms of characteristic polynomials f (E, e ik ) that can be achieved i.e. to that of the more general form E 2 − T (k)E + det(H(k)) = 0, with two independent functional degrees of freedom T (k) and det(H(k)) in deciding the form of the Hamiltonian. With added algebraic complexity, this can analogously be extended to Hamiltonians of higher bands.
While a multi-component Hamiltonian gives rise to multiple energy bands, only one band can correspond to the desired i.e. stipulated spectrum. In many experimental realizations, the other bands can be left unoccupied and thus ignored. In implementing T (k), one finds that typically only terms up to the next-nearest hoppings are needed to converge closely to the desired Hamiltonian. In Fig. S2, we show how this can be per-formed to achieve a branch of almost perfectly real OBC spectrum, with the Hamiltonian ansatz being given by H(k) = H z (k)σ z + T (k)I + F(k)σ + + σ − (Eqn. 6) with: Ultimately, the choice of ansatz is also heavily influenced by practical physical considerations, such as the structure and nature of the physical lattice used to realize the engineered Hamiltonian. For instance, asymmetric couplings are easily realized in classical electrical circuits, but less straightforwardly in quantum circuits [74,75].

B. Exploiting n-fold rotational symmetry
In cases where the spectra exhibit n-fold rotational symmetry in both the PBC loop and OBC loci, one may expand the scope of candidate Hamiltonians by partitioning the energy spectrum into equally spaced intervals of k = 2π/n where n is the order of symmetry by rotation, such that one interval coincides with the stipulated spectrum. We demonstrate this method with the Hamiltonian H = z + 1/z 2 previously presented in Fig 4. Visually, the PBC spectra consists of 3 equally sized "lobes" offset by a rotation 2π/3. One may break each lobe into separate momentum intervals (0, 2π/3), (2π/3, 4π/3), (4π/3, 2π). The benefit of doing so is illustrated and elaborated in Fig. S3. By segmenting the fourier series into piecewise separate piecewise functions, one avoids the Gibbs phenomenon-where arbitrarily large partial sums fail to converge point-wise to the target function.  At the sharp jump, the Gibbs's phenomenon can be observed where the difference converges to a finite value instead of 0. This is contrasted by i and k where only 2 out of 20 Fourier terms are used are used and j and l show that the difference in the approximation is 1-2 orders smaller than the preceding case.

V. NON-BLOCH BAND COLLAPSE AND THE METHOD OF IMAGES
An interesting limiting case is the one where the OBC energy spectrum collapses to one point E = E 1 (or to a set of points) in the interior of the PBC energy spectrum, corresponding to so-called non-Bloch band collapse or flat band limit. In the electrostatics analogy, the energy E = E 1 is a singularity of the potential V (vanishing of the skin depth or κ = V = ∞) due to a point charge placed at the point z = z 1 = E 1 . A useful analytical tool to solve certain electrostatic problems involving point charges in the presence of conductors is the method of images, where insertion of additional virtual point charges makes the conductor surface isopotential. This method solves, for example, the problem of a point charge or an electrostatic dipole placed in front of a plane (flat) conductor or in the interior of a sphere. In the spirit of our analogy, the method of images can be harnessed to synthesize a NH Hamiltonian displaying a non-Bloch band collapse (flat band) at some target real energy E 1 . For illustrative purposes, let us synthesize using the method of images a NH Hamiltonian E = H(p), with complex momentum p = k + iκ, such that the PBC energy spectrum, i.e. E = H(p) with κ = 0 and p = k, describes the circle of radius R, centered at the origin of the complex energy plane, whereas the OBC energy spectrum collapses to the real energy E = E 1 in the interior of the circle, i.e. −R < E 1 < R (see Fig. S4). In the electrostatics analogy, the complex momentum p corresponds to the complex potential φ, i.e. p = φ, whereas energy E corresponds to the complex spatial position z = x + iy, i.e. E = z. We are thus faced with the problem of finding a complex potential φ(z) = U (z) + iV (z) satisfying the following requirements (i) φ(z) shows a singularity at z = E 1 ; (ii) ∇ 2 z φ = 0 in the interior |z| < R and z = E 1 , (iii) on the circle |z| = R the potential V (z), i.e. imaginary part of φ(z), must vanish (grounded conductor). From the wellknown form of the complex potential of a point charge in two dimensions, we can satisfy the above requirements by looking for a solution of the form where Q is the point charge at the singularity z = E 1 , internal to the conductor, q is the image charge, external to the conductor and placed at a distance E 1 , and φ 0 is a constant potential term (see Fig. S4). Clearly, ∇ 2 z φ(z) = for |z| < R and z = E 1 . The image point charge q, its distance E 1 from the center of the conductor and the constant potential term φ 0 should be determined by imposing conditions (ii) and (iii). Clearly, on the conductor z = R exp(iϕ) one has are the squares of the distances ρ 1 = |z−E 1 |, ρ 2 = |z−E 1 | of the two point charges from the point z = R exp(iϕ) on the conductor. Eliminating ϕ from Eqs.S6 and S7, one obtains S4. Method of images and the synthesis of a NH lattice displaying non-Bloch band collapse. A point charge Q is placed in the hollow region of a grounded (zeropotential) 2D conductor surface (the circle of radius R). The charge is placed at a distance E1 from the center. The electrostatic problem is solved by considering a virtual (image) point charge q, in the exterior of the circle, at a distance E 1 .
In order for V (ϕ) to be independent of ϕ, we set Under such conditions, from Eq.S5 one obtains so that condition (iii) is met by letting φ 0 = (Q/2π) log(R/E 1 ). The complex potential φ(z) thus reads The Hamiltonian E = H(p) is finally obtained by solving Eq.S12 with respect to z, and letting z = E = H(p), φ = p. This yields The point charge Q should be chosen such that, for p = k real, the exponential function exp(−2πip/Q) is invariant for k → k + 2π. A possible choice is Q = 2π, so as one finally obtains In physical space, this Hamiltonian is realized in a lattice with hopping amplitudes t n defined by the Fourier series H(k) = n t n exp(ikn), i.e.
In the limiting case E 1 = 0, i.e. when the OBC energy spectrum collapses at the center of the PBC energy spectrum, one has H(p) = −R exp(−ip), which basically describes the Hatano-Nelson model in the limit of a unidirectional nearest-neighbor hopping R. For E 1 = 0, there are non-nearest hopping. To compute t n , one can calculate the integral on the right hand side of Eq.S15 by the residue theorem. For n ≥ 0, after letting ξ = exp(−ik) one has S16) For n < 0, after letting ξ = exp(ik) one obtains Interestingly, for E 1 = 0 the hopping in the lattice remains unidirectional, since t n = 0 for n > 0 while t n = 0 for n ≤ 0, with |t n | exponentially decaying toward zero as n → −∞. The matrix Hamiltonian in physical space is thus an upper triangular matrix with E 1 on the main diagonal. Under OBC, the L × L matrix is defective with all eigenvalues collapsing to the value E 1 , which is a Lorder exceptional point, L being the number of lattice sites. Finally, we mention that, owing to the superposition principle, the previous result can be readily extended to the case of an arbitrary number N of point charges in the interior of the circle, placed at the complex energies (positions) E 1 , E 2 ,..., E N . The corresponding synthesized lattice shows a non-Bloch collapse under OBC to the energies E 1 , E 2 ,..., E N . In this case, assuming again Q = 2π, the complex potential φ(z) reads and the Hamiltonian H = H(p) is implicitly defined as a root of the algebraic equation (S19)

VI. DETAILS OF ENGINEERED HAMILTONIANS
This section gives explicit expressions for the various engineered Hamiltonians designed in the main text, namely, those in Fig. 3a, 3b and 4b. These Hamiltonians were engineered via our electrostatics approach, based on stipulated spectral loci and inverse state localization lengths κ that were used as inputs.
For the OBC and PBC spectra recovered in Fig. 3a, the engineered Hamiltonian takes the following ansatz with trace and off-diagonal terms truncated to Due to the degrees of freedom afforded by the presence of both T (k) and F(k), good agreement between the stipulated and engineered spectra is achieved with hoppings of up to next-nearest neighbors only. The engineered Hamiltonian featured in Fig. 3b takes the form and gives a good fit with the stipulated real OBC spectrum despite being a simple single-component model with up to next-nearest neighbor couplings. This example one of the several further presented below, in Fig. S2, where different parameters for the real OBC spectral segment was stipulated.
Finally, the engineered Hamiltonian for Fig. 4b is S23) which exhibits excellent stipulated vs. engineered spectra agreement, and contains hoppings only within the same unit cell and across three unit cells.

VII. VERIFICATION OF κ-PROFILE
In this section, we seek to verify that the stipulated engineered skin depth distribution κ(E) along the OBC line indeed agree. It is known from the non-Hermitian Skin Effect that each state will decay exponentially from the system edge. In other words, pumping towards the edge causes an exponential accumulation of the edge state. By expressing the Hamiltonian in position space, one may obtain the eigenvalues and just as importantly the eigenstates in real space, where their spatial profiles can be compared with the stipulated κ decay.
We systematically log-plot the eigenvectors arranged by site and obtain the exponential decay factor κ. For the Hatano-Nelson warm-up model, the decay factor remains constant across all sites, which corresponds to the original stipulated potential V = 1 in the electrostatic model up to normalization. This is expected since the Hatano-Nelson model is known to have a constant κ. Doing the same for the model (Fig. 3a of  with two disjointed OBC lines with spatially inhomogeneous potential i.e. κ profile, we find that for V = 1 and V = 1.6 respectively that the average of κ across each line also agrees up to normalization. However, the exact profile need not agree entirely and is not expected to do so given the system does not respect 2-fold rotational symmetry.

VIII. DETAILS OF ENGINEERED SPECTRA IN THE PHASE DIAGRAM IN FIG. 3C OF THE MAIN TEXT
In Fig. 3c of the main text, we studied how reliably our approach can engineer local Hamiltonians with real spectra, in a setting with no manifest symmetry that encourages real spectra. We used a Reuleaux triangle as the PBC spectral locus, and stipulated the target OBC eigenstates to have a maximum κ(E) amplitude given by A and OBC locus to be a real line segment of length d, as given in Eq. 8 of the main text.
Due to the requirement of sufficient locality, some minute discrepancy from a purely real spectrum is unavoiable. Plotted in Fig. 3c of the main text, which is reproduced here as Fig. S5a, is the maximum value of Im(E) in the parameter space of A and d, for hoppings truncated beyond 5 sites. The plot is helpful in searching for a parameter region in which the energies were purely real and thereby admitting experimentally realizable states. 9 We do this by computationally sweeping a grid consisting of A ∈ (−1, 4) and d ∈ (0, 5.5) while setting the left most point of the OBC loci to be at (−1, 0), the choice of which was somewhat arbitrary. In Fig. S5a we see a region in which the maximum magnitude of Im(E) is very close to 0. Each point within the parameter space corresponds to a specific spectral configuration. A desirable configuration is already given in Fig. 3b of the main text ; some examples in which the imaginary eigenenergies are slightly non-vanishing are listed in Fig  S5b-f. Here we have arbitrarily cut off the number of Fourier modes (coupling distance) at 5 and it is expected that as we increase the Fourier modes that the PBC spectrum will come closer to the stipulated PBC spectral loci.  Labelled site numbers corresponding to the eigenvalues marked in red with their respective eigenvectors plotted over every site shown below. The log of each eigenvalue was fitted to a straight line whose gradient is κ. Within the bandstructure, the value for κ at each site was then plotted across the Re(E). b A similar procedure was performed for the case of two disjointed OBC lines (Fig. 3a of the main text) with potentials set at V = 1 for the leftward line and V = 1.6 for the rightward line. The lack of 2-fold rotational symmetry gives rise to the κ profile observed. However, averaged across the line, the value of κ across the rightward line is twice the value as the leftward line.

IX. CIRCUIT REALIZATION OF NON-HERMITIAN LATTICES
Topolectrical circuits [76][77][78][79][80][81][82] have provided a fresh avenue into realizing new phenomena. Their extreme versatility in the implementation of distant couplings, as well as precise control over linearity and/or non-linearity, has brought to reality many models that are otherwise hard to replicate in conventional materials or metamaterials. In this section, we present a generalized circuit implementation for realizing models containing up to next-nearest neighbor hoppings, including many of our designed models.
Any electrical circuit can be represented as a graph. Such a representation can be described by Kirchhoff's laws in nodal form [83,84].
Here, c ij represent the admittances of the circuit elements between nodes i and j, V i and V j are the electric potentials of nodes i and j, and I i is the current flowing into the node from the sources connected directly to node i. In compact matrix notation, Kirchhoff's law for a list of N nodes can be expressed as where L is the circuit Laplacian. Written explicitly in terms of the matrix elements, the equation reads Notice that if one strategically matches the admittances of each node in the Circuit Laplacian with numerical values that correspond to the original Hamiltonian, that we have effectively reproduced the Hamiltonian, albeit with a different physical interpretation as a currentvoltage Laplacian. Most often, such circuits are implemented with capacitors and inductors with admittances (impedances) that depend on the frequency of a driving alternating current. Furthermore, non-hermiticity can be achieved with non-reciprocal elements that asymmetrically control the current flow in one direction differently from the other. In practice, this can be implemented with externally powered directional operational amplifiers whose amplification ratio can be finely tuned by the (a) (c) S7. a. Generalized circuit realization of any 2-band, nearest neighbour hopping and b. next-nearest neighbour circuit. c. Schematic representation of the non-reciprocal element that leads to non-hermiticity using Negative Impedance Converters (INICs). d. Grounded nodes with each unit cell attached to active INIC elements effectively shift the total energy of the system, serving as an enabler for the realization of trace terms introduced in Sec III E.
user. With reference to [85], the circuit Laplacian reads In order to realize for instance the Hamiltonian in Eqn. 9, one may simply tune the experimental values of the electrical components in Eqn. S27 to match the coefficients of the engineered Hamiltonian circuit with experimental components. In an engineered Hamiltonian where a trace ansatz was used, a grounded node such as the one in Fig. S7c can be used to set a trace term of choice To measure the Laplacian eigenspectrum, one can do so from the impedance readings between all pairs of nodes, effectively measuring the circuit Laplacian node by node. The impedance Z a0 between two nodes a and a ground is given by where ψ a,i and φ a,i are the left and right eigenvectors of the Laplacian respectively and λ i the eigenvalues of the Laplacian, which can then be solved by inverting the data obtained.

X. ADDITIONAL EXMAPLES
Our electrostatics formalism goes beyond examples with real OBC spectra, and can readily inverse engineer a complex spectra. For concreteness, we consider the onedimensional nonreciprocal superlattice model proposed by Zeng et al [86], where the nonreciprocal hopping is modulated by a cosinusoidal function. By tuning the modulation of the hoppings, they demonstrate robust real energy OBC spectra. Here, we would like to reverse engineer the Hamiltonian parameters given the spectra and the form of the Hamiltonian. We choose the parameters α = 1/3, γ = 2, δ = 0.6π, t = 1, which yields a periodic chain with 3 sites in each unit cell, i.e. a 3-band model. The OBC spectrum consists of a finite real line and two parabolas, while the shape of the PBC spectrum is more sophisticated than our warm-up examples. The OBC and PBC spectra are given in Fig. S8.
To reiterate, suppose the parameters γ and δ are not known a priori, we seek to find these 2 parameters given only the spectra and the form of the bulk momentum space Hamiltonian, given by where γ j = γ cos(2παj + δ) is the modulation. By setting t = 1 without loss of generality, we yield the characteristic equation where c = 3−γ 2 1 −γ 2 2 −γ 2 3 , g = −(1−γ 1 )(1−γ 2 )(1−γ 3 ) and h = −(1 + γ 1 )(1 + γ 2 )(1 + γ 3 ). Given the numerical values of g and h, we may obtain γ and δ via a root finding algorithm. Thus, we will focus only on determining g and h. Note that by taking the real (imaginary) parts on both sides of the characteristic equation, we end up with a sinusoid of amplitude g + h (|g − h|) on the right-hand side. Thus, to find the values of g and h, we must find E(k) for one energy band (any suffice) and determine the numerical value of the constant c such that −E 3 + cE is purely sinusoidal. After this value of c is determined, the maximum value f Re[−E 3 + cE] can be computed; this is g + h. The maximum value of Im[−E 3 + cE] is similarly g − h.
To find the energy dispersion of a single band, we follow the same procedure as described in our paper (Fig. S1 Steps 1 to 4). In our formalism, we prescribe that the entire OBC line has a constant skin depth which we can set to V = 1 without loss of generality. The PBC line is grounded. The resulting solution of the Laplace Equation is Fig. S8b. As mentioned, the density of states (Fig S8c) is proportional to the charge density. The charge density is much higher on the y = 0 line due to the close proximity to the OBC line (and hence a higher ∂E ∂x ). From the density of states, we can then assign a region of 1/3 the length of the PBC line to be a single energy band, from which we can compute the values of the constants g and h. To clarify, Step 5 in Fig S1 is only done to construct a Hamiltonian that recovers the same OBC and PBC spectra, which is typically not unique. In this example, our ansatz is the form of the Hamiltonian itself H(k) and we seek to recover the parameters γ and δ.
In general, our method could be readily applied to high-dimensional systems where the Hamiltonian is separable. One example is the 2D Hatano-Nelson model on a square lattice with asymmetric left/right and up/down hopping amplitudes. For generic higher-dimensional systems, we could formulate the system as an array of 1D chains, again with asymmetric hopping amplitudes.