High-order exceptional points in supersymmetric arrays

We employ the intertwining operator technique to synthesize a supersymmetric (SUSY) array of arbitrary size $N$. The synthesized SUSY system is equivalent to a spin-$(N-1)/2$ under an effective magnetic field. By considering an additional imaginary magnetic field, we obtain a generalized parity-time-symmetric non-Hermitian Hamiltonian that describes a SUSY array of coupled resonators or waveguides under a gradient gain and loss; all the $N$ energy levels coalesce at an exceptional point (EP), forming the isotropic high-order EP with $N$ states coalescence (EPN). Near the EPN, the scaling exponent of phase rigidity for each eigenstate is $(N-1)/2$; the eigen frequency response to the perturbation $\epsilon$ acting on the resonator or waveguide couplings is $\epsilon^{1/N}$. Our findings reveal the importance of the intertwining operator technique for the spectral engineering and exemplify the practical application in non-Hermitian physics.

Different types of energy level coalescence exist in non-Hermitian systems. The most common types of EPs are the two-state coalescence (EP2) that exhibits a squareroot dependence on the system parameters [6][7][8]30] and * jinliang@nankai.edu.cn the three-state coalescence (EP3) that exhibits a cubicroot dependence on the system parameters [31,[45][46][47][48]. Even four-state coalescence (EP4) are accessible in the coupled resonators [49,50]. Recently, the high-order EP of arbitrary order is realized in coupled resonators [51]. And the scaling law for the eigenvalue and eigenstate confirm the sensitive property of the high-order EP [52][53][54][55][56][57]; the dynamics near a high-order EP exhibits a power law dependence on the order of EPs for the maximal amplification [58].
The intertwining operator technique is a useful method for the spectral engineering [59,60]. Such technique is capable of eliminating a target energy level in the spectral to create the isospectral SUSY partner [61][62][63][64], which has an identical spectrum except for the eliminated target level. In parallel, the intertwining operator technique can add target energy level or realize Hamiltonian with spectrum fully constituted by the desirable energy levels [60]. Exact solvable models with desirable energy levels can be synthesized employing the intertwining operator technique. Thus, the intertwining operator technique is beneficial for proposing non-Hermitian Hamiltonian with multiple energy level coalescence. The concept of supersymmetry, originated from quantum field theory [65], has boomed during recent years in the research fields of optics and photonics. It is possible to create designed spectrum and propose intriguing applications using the synthesized system. The SUSY array synthesized through intertwining operator technique can be utilized for optical sensing [30,31], single mode lasing [61][62][63], and optical mode converting [64]. The integrated lasing array usually has multiple mode emission. To acquire single mode lasing, one can design an isospectral partner SUSY array using the intertwining operator technique. The spectrum of the partner SUSY array is engineered to be constituted by all the excited-state mode except for the ground-state mode of the lasing array. Coupling the partner SUSY array to the lasing array and intentionally inducing loss in the partner SUSY array can enable the ground-state single mode lasing [61][62][63]. The SUSY array can remove the ground-state mode of a multimode light field and manipulate the modal content of the light field through a hierarchical sequence of partner SUSY arrays [64]. The synthetic SUSY is a hypercube useful for quantum information science [66], and the proposed SUSY array has equally spaced energy levels; thus the SUSY array is capable of realizing a perfect state transfer that the initial state is exactly mapped from one side of the SUSY array to the other side of it [53].
In this paper, we introduce the intertwining operator technique to propose a non-Hermitian SUSY array of arbitrary size. The energy levels of the proposed SUSY array are equally spaced square-root branches. The non-Hermitian phase transition in the proposed SUSY array is associated with an isotropic high-order EP. In contrast to the anisotropic high-order EP [55], the isotropic highorder EP has identical parameter dependence for independent system parameters in the parameter space. The Hamiltonian of the SUSY array can be understood either as many uncorrelated spin-1/2 particles in a magnetic field or the noninteracting bosonic many-particle system in a two-site model. For the arbitrarily high-order EP, topological properties and the frequency response to perturbation on the resonator couplings are investigated. In the SUSY array of N sites, the isotropic EPN has the phase rigidity scaling exponent (N − 1) /2. The eigen frequency sharply responses to the coupling perturbation with the form 1/N near the EPN. The results are in accord with topological features of the EPN. The SUSY is proper for the perfect state transfer in quantum information science.
The remainder of the paper is organized as follows. In Sec. II, we introduce the intertwining operator technique to synthesize the PT -symmetric SUSY array. In Sec. III, we investigate the topological properties of the arbitrarily high-order EP through the phase rigidity. In Sec. IV, we focus on the eigen frequency response to the coupling perturbation to reflect features of the arbitrarily highorder EP. In Sec. V, we summarize the results.

II. PT -SYMMETRIC SUSY ARRAY
In this section, we introduce the intertwining operator technique to propose a synthetic SUSY array with all the energy levels equally spaced. The synthetic SUSY is a hypercube, and the hypercube of arbitrary dimension can be synthesized. The schematic of a synthetic sixsite SUSY array of coupled resonators (upper panel) or coupled waveguides (lower panel) is shown in Fig. 1(a). The frequency of each resonator (waveguide) is ω 0 . The resonators (waveguides) are coupled through evanescent tunneling between neighboring ones. The coupling amplitudes of the SUSY array are determined from the intertwining operator technique. In Fig. 1(b), schematics of energy levels of the non-Hermitian SUSY arrays with different site numbers are shown; all the levels are equally spaced and the energy difference between each pair of neighboring energy levels is 2J. To constitute the SUSY 1. (a) The SUSY array has a linear gradient on the gain and loss distribution. The resonant frequency is ω0 for each resonator or waveguide. (b) Proposal of engineering the SUSY array.
array, the couplings are required to be engineered at the proper amplitudes.
In the framework of quantum mechanics, the procedure is as follows. The Hamiltonian h N has N energy levels and we aim to remove an energy level ε N from h N to construct a target superpartner h N . The Hamiltonian h N is factorized into where The target superpartner Hamiltonian is obtained as To synthesize the SUSY array, we consider an inverse process and gradually increase the size of the target Hamiltonian by adding target energy levels one by one. Details of synthesizing the SUSY array are displayed as follows. We start with a single-level system h 1 , shift its energy level by 2J and add a zero energy (ε 2 = 0) to obtain energy spectrum {0, 2J} of the target Hamiltonian h 2 . We factorize the Hamiltonian The factorization is in the form of Then, following the intertwining operator technique, we interchange R 1 and Q 1 in the matrix product to obtain the target Hamiltonian The target Hamiltonian is a 2 × 2 matrix where σ x is the Pauli matrix. Next, we offset the energy by −J to get h 2 − JI 2 = Jσ x and gather the Hermitian SUSY array of N = 2. The matrix form of h 2 after subtracting out the term JI 2 is equivalent to the single-particle Hamiltonian for the two-site model. Jσ x also describes a spin-1/2 particle in an effective magnetic field J along the x direction, that is (B x , B y , B z ) = (J, 0, 0). Furthermore, in order to construct the non-Hermitian Hamiltonian with the EPs, we consider an additional effective imaginary magnetic field iγ applied along the z direction; the PT -symmetric non-Hermitian dimer is obtained with the expression Jσ x + iγσ z and the corresponding matrix form is given by where ω 0 is the on-resonator frequency of each resonator or waveguide. H 2 has a pair of EP2s at γ = ±J, where two eigenstates coalesce to one (±i, 1) T / √ 2. We repeat the above procedure to construct the non-Hermitian SUSY array with high-order EPs. The spectrum of h 2 is shifted by 2J to {2J, 4J}; and then ε 3 = 0 is added to obtain the spectrum {0, 2J, 4J} of the target Hamiltonian h 3 [ Fig. 1 where δ means the Kronecker delta function. The factorization gives Following the intertwining operator technique, we obtain a 3 × 3 Hamiltonian after interchanging R 2 and Q 2 in the matrix product, that is The target Hamiltonian h 3 can be expressed in the form is the angular momentum operator for spin-1 in the x (z) direction. We remove the overall energy background 2JI 3 from h 3 to obtain the Hermitian SUSY array of N = 3. h 3 − 2JI 3 = 2JS x describes a spin-1 particle in an effective real magnetic field (J, 0, 0). By considering an effective imaginary field iγ applied in the z direction (0, 0, iγ), we obtain the PT -symmetric non-Hermitian trimer JS x + iγS z . The corresponding non-Hermitian SUSY array is The generalized SUSY array here is the two-particle Hamiltonian for the two-site model in non-Hermitian cases, with the matrix form corresponding to the Hamiltonian described by Eq. (8). The quantum theory of spin (angular momentum) tells us that the energy spectrum of the Hamiltonian JS x +iγS z is restricted to n J 2 − γ 2 with n = −2, 0, 2. Each EP of H 3 belongs to EP3 with the eigen frequency ω 0 and coalescence eigenstate (−1, ±i √ 2, 1) T /2 at the critical condition γ = ±J. To synthesize the 4 × 4 SUSY array, we shift all the energy levels of h 3 by 2J to obtain {2J, 4J, 6J} and add a zero level to compose the spectrum {0, 2J, 4J, 6J} of the target Hamiltonian h 4 . We factorize h 3 in the form of Following the intertwining operator technique, we obtain the 4 × 4 target Hamiltonian We offset energy 3J from h 4 to obtain the Hermitian SUSY array of N = 4, which describes a spin-3/2 particle in an effective real magnetic field (J, 0, 0). The non- is the angular momentum operator for spin 3/2 in the x (z) direction. The matrix form of the SUSY array is The J related terms are the couplings between neighbor resonators induced by the evanescent fields, and depend on the distance between them. The distribution of gain and loss in the SUSY array has a gradient. The PT -symmetric non-Hermitian 4×4 SUSY array is the three-particle Hamiltonian for the two-site model; the spectrum of which is restricted to n J 2 − γ 2 with n = −3, −1, 1, 3, with the appearance of EP4 at γ = ±J. We can synthesize the SUSY array of arbitrary size via repeating the same procedure. The PT -symmetric non-Hermitian 5 × 5 SUSY array has the form The PT -symmetric non-Hermitian 6 × 6 SUSY array illustrated in Fig. 1(a) has the form The SUSY arrays H 5 and H 6 are the four-particle and five-particle Hamiltonians for the two-site model [57], respectively. In both cases, the EP occurs at γ = ±J, but being EP5 in H 5 and EP6 in H 6 . Currently, experimental investigation on the high-order EP has sprung up rapidly. For example, a fifth-order EP can be designed via tuning parameters in nitrogen-vacancy centers [67]; a PT -symmetric electronic circuit has been proposed to study sensing at a sixth-order EP [68].
In general, we shift all the energy levels of h N −1 by 2J to obtain a spectrum {2J, 4J, · · · , 2(N − 1)J} and add ε N = 0 by the intertwining operator technique to obtain h N with the spectrum {0, 2J, · · · , 2(N −1)J}. We obtain is the angular momentum operator for spin (N − 1) /2 in the x (z) direction. Removing the offset energy (N − 1) J and introducing an imaginary magnetic field in the z direction, we obtain the PT -symmetric non-Hermitian N × N model JS x + iγS z [53]. H N describes a N -site PT -symmetric non-Hermitian SUSY array [62][63][64]. The imaginary magnetic field corresponds to tilted on-site imaginary potentials in the form of gain and loss, which linearly depends on the site number. The concise form of H N is given by The non-Hermitian generalized SUSY array is synthesized by a recursive bosonic quantization technique in the coupled resonators or waveguides [54].
The first line in H N is the Hermitian SUSY array [66]. The energy of JS x relates to the (quantized) possible value of spin angular momentum in the x direction, being nJ for the integer n = − (N − 1) , − (N − 3) , · · · , (N − 3) , (N − 1). H N describes a spin-(N − 1)/2 particle in an effective magnetic field (B x , B y , B z ) = (J, 0, iγ). This indicates that the SUSY array H N has the frequency Notably, the EPs of H N are exactly EPNs at γ = ±J, where all the levels are square-root branches and coalesce to the resonant frequency ω 0 . Furthermore, introducing the angular momentum op- ≡ H two−site can be alternatively understood as the noninteracting bosonic many-particle system in a PT -symmetric non-Hermitian two-site model [57], where a † 1(2) and a 1 (2) are the creation and annihilation operators of the first (second) site, respectively. [S a , S b ] = 2i abc S c , where abc is the Levi-Civita symbol and a, b, c are x, y, z. The commutation relation is equivalent to that of Pauli matrices with spin 1/2.
The basis set for the single particle system is chosen as (5) is H two−site in the single-particle basis. If we consider the two-particle problem, the basis set is The factor 1/ √ 2 in the basis is not only the normalization factor to ensure 2 m|n 2 = δ mn (m, n = 1, 2, 3); but also the normalization factor in the Fock representation for the occupation number of two. H 3 in Eq. (8) is H two−site in the twoparticle basis. Moreover, the basis set for three particle is (11) is H two−site in the three-particle basis. H 2 [30], H 3 [31], and H 4 [52] are experimentally realized in different physical setups. In general, the basis for (N − 1)-particle system is |l N −1 = (a † 1 ) N −1−l (a † 2 ) l / (N − 1 − l)!l! |vac , where l = 0, 1, · · · , N − 1 and the subscript in the basis stands for the number of particles. H two−site in the basis {|l N −1 } gives H N in Eq. (14).
The notion of SUSY plays an important role for plenty of intriguing optical properties and functionalities as well as for a number of practical applications of optical metamaterials [61][62][63][64]. The non-Hermitian SUSY array also provides a promising platform for the study of the topology of arbitrary high-order EP. The geometric topological properties reflect the order of EPs and are the essential features of different EPs, which are captured by the phase rigidity scaling exponents.
The basis and eigenstates of the (N − 1)-particle Hamiltonian are chosen under the Fock representation. Considering the direct product representation that employs the single-particle eigenstates as the basis, the expression of eigenstates of the (N − 1)-particle system is the direct product of the N −1 numbers of single-particle eigenstates |u 1 and |u 2 . The general expression of the normalized eigenstate of a spin-(N − 1) /2 system in the direct product representation is (|u ; the scaling exponent ν for EPN is ν = log 10 |r l 1 r N −1−l 2 |/ log 10 |γ − γ EP | = (N − 1) /2 (see Appendix A). The analysis on the scaling exponents of higher-order EPs is numerically verified in Fig. 2, where the phase rigidities and scaling exponents of EP3, EP4, EP5 and EP6 under the Fock representation Hamiltonian [Eq. (14)] are shown for N = 3, 4, 5, 6. The scaling exponents are 1.0, 1.5, 2.0 and 2.5 for the EP3, EP4, EP5, and EP6, respectively. Notably, there exists a zero-energy flat band for odd N , which participates in the coalescence of eigenstates (see Appendix B). Thus it exhibits an identical scaling behavior to other levels. We emphasis that the high-order EPs in the SUSY array are isotropic. Replacing iγ by ∆ + iγ, we observe the scaling exponent ν = (N − 1) /2 for tuning the detuning ∆ to approach the EPNs |r| = |∆ − ∆ EP | ν . The scaling exponent of the phase rigidity is robust to the perturbation [55,71]. For an anisotropic EPN in the system with asymmetric couplings, the scaling exponents of the phase rigidity are (N − 1)/2 and N − 1 when approaching EPN from two independent parameters, respectively [55].

IV. EIGEN FREQUENCY RESPONSE TO PERTURBATION
The non-Hermitian SUSY array at the high-order EP enhances the susceptibility in optical sensing, the frequency response near the high-order EP is greatly increased [30,31,44,71]. Near the high-order EP in the non-Hermitian SUSY array, a remarkable point is the enhanced frequency response to the detuning as well as the coupling when the array is subjected to the perturbation . The SUSY array is a hypercube with high symmetry, the response to acting on the coupling appears similar response to that acting on the detuning when approaching an EPN. And the eigen frequency response is ∼ 1/2 for EP2 [30,39] and ∼ 1/3 for EP3 [31], which is distinct from the linear response to the perturbation strength ∼ near the degeneracy point in Hermitian systems. The sharp response is a typical feature of the EP that paves the way of the application of sensors. Moreover, the SUSY array we constructed holds an arbitrary EPN (N ≥ 2) at γ = ±J, which has striking features. According to the Newton-Puiseux series expansion [72,73], the frequency splitting ω l is a function of the perturbation c 1 1/N + c 2 2/N + · · · [71], where the corresponding coefficients c 1 , c 2 , · · · are complex numbers. In the perturbation theory, the unperturbed Hamiltonian is H N , and the perturbation Hamiltonian is H Notably, in general cases, ω l ≈ c 1 1/N for comparable c 1 ∼ c 2 . We also show the eigen frequency response to the coupling perturbation characterized by the order of magnitude ∼ 1/N near an EPN in the SUSY array is similar to that of the detuning perturbation.
The coupling is determined by the distance between resonators or waveguides. In practice, the resonator frequency can be accurately fabricated, and the coupling J may have imperfections. Thus, it is reasonable to consider the resonator coupling perturbations for the EPN of the SUSY array H N . We take the example in Fig. 3 that all the coupling terms have perturbations, that is, the J related terms are J m (N − m)+ (m = 1, 2, · · · N −1). Figure 3 depicts the energy levels, frequency splitting, as well as the logarithmic plots of the frequency splitting as a function of coupling perturbation . The coupling perturbation presents in each coupling J m (N − m). The whole spectrum is symmetric about zero energy due to the equal amount of chosen without breaking the chiral symmetry of H N . The frequency splitting response proportional to 1/N is shown in Fig. 3; this differs from the response to the coupling perturbation J + instead of J in Eq. (15), which leads to the square-root dependence ∼ 1/2 . The frequency splitting dependence on the coupling perturbation between the first two cavities exhibits slopes of 1/3 and 1/4 for H 3 and H 4 in the logarithmic plots; this reveals the cubic-root and quartic-root dependence near the EP3 and EP4 (see Appendix C). In general cases, the coupling perturbation in an arbitrary resonator of the SUSY array leads to the similar response ω l ≈ c 1 1/N ; besides, we can observe the response is in the order of ω l ≈ c 2 2/N . The frequency splitting response as a function of coupling perturbation is identical to the detuning perturbation [56,57], which reflects the high symmetry feature of the hypercube (SUSY array).

V. CONCLUSION
The intertwining operator technique is an important approach for the spectral engineering. We employ the intertwining operator technique to propose the non-Hermitian SUSY array with arbitrarily high-order exceptional points. The Hamiltonian of the proposed array is a non-Hermitian generalized SUSY lattice chain for perfect state transfer in quantum information science, which is equivalent to a noninteracting many-particle Hamiltonian of the two-site non-Hermitian PT -symmetric dimer. At the EPN of the SUSY array with N coupled resonators or waveguides, all the energy levels are equally spaced, being square-root branches and coalescing at the EPN. The phase rigidity of each eigenstate reaches zero and the scaling exponent is (N − 1) /2 for the EPN; the eigen frequency response to perturbation is 1/N for coupling amplitude perturbation in certain resonators or waveguides of the SUSY array. The intertwining operator technique provides a promising method for synthesizing artificial optical metamaterial.
In general, the eigenstate for the N × N Hamiltonian H N with eigenvalue lε 1 + (N − 1 − l) ε 2 is given by (|u 1 ) l ⊗ (|u 2 ) N −1−l , where l = 0, 1, · · · , N − 1. Therefore, the phase rigidity vanishes with r EP = 0 at EPN γ EP = J and the corresponding phase rigidity r is The scaling exponent ν for EPN is ν = log 10 r l 1 r N −1−l 2 log 10 |γ − γ EP | = log 10 |r 1 | l log 10 |γ − γ EP | + log 10 |r 2 | APPENDIX B: ENERGY BANDS FOR ODD N A zero-energy flat band exists in H N when N is odd because H N has the chiral symmetry; therefore, the spectrum of H N is symmetric about the zero energy, and there is a zero-energy flat band if the system has an odd number of energy levels. At the EP γ = ±J, all energy levels coalesce, H N is nondiagonalizable and reduces into a Jordan block   Imposing perturbation on the couplings of the first two resonators, we can get the matrix form for Hamiltonians H 3×3 , H 4×4 and H 5×5 , Figure 5(c) represents the frequency splittings of bands with cyan solid and yellow dashed lines in Fig. 5(a) and Fig. 5 (b) for H 3×3 , the real and imaginary parts of which are depicted in cyan solid and orange dashed lines. The logarithmic relationship between the frequency splitting and disturbance is shown in Fig. 5(d), where the slope 1/3 indicates the order of EP3.
Two energy bands with the maximum and minimum real parts in Fig. 6(a) hold the identical positive imaginary parts in Fig. 6(b); the middle two bands in 6(a) correspond to the same negative imaginary parts in Fig. 6(b). Figure 6(c) represents the frequency splittings of energy bands with cyan solid and orange dash-dot lines in Fig. 6(a) and Fig. 6(b) for H 4×4 , the real (imaginary) part of which is depicted with cyan solid (orange dashed) line. The corresponding logarithmic relationship with is shown in Fig. 6(d), where the slope 1/4 indicates the order of EP4. Figure 7(c) represents the frequency splittings of energy bands with cyan solid and green dashed lines in Fig. 7(a) and Fig. 7(b) for H 5×5 , the real and imaginary parts of which are depicted in cyan solid and orange dashed lines. The logarithmic relationship between the frequency splitting and perturbation is shown in Fig. 7(d), where the slope 1/5 indicates the order of EP5.