Non-Hermitian physics of levitated nanoparticle array

The ability to control levitated nanoparticles allows one to explore various fields of physics, including quantum optics, quantum metrology, and nonequilibrium physics. It has been recently demonstrated that the arrangement of two levitated nanoparticles naturally realizes the tunable nonreciprocal dipole-dipole interaction. Motivated by this development, we here propose and analyze an array of levitated nanoparticles as an ideal platform to study non-Hermitian physics in a highly controlled manner. We employ the non-Bloch band theory to determine the continuum bands of the proposed setup and investigate the non-Hermitian skin effect therein. In particular, we point out that the levitated nanoparticle array exhibits rich dynamical phases, including the dynamically unstable phase and the unconventional critical phase where the spectral singularity persists over a broad region of the controllable parameters. We also show that the long-range nature of the dipole-dipole interaction gives rise to the unique self-crossing point of the continuum band.

In this paper, we propose and analyze a 1D levitated nanoparticle array as an ideal platform to study previously unexplored regimes of non-Hermitian physics in a highly controlled manner.A prominent feature is that there exists the tunable nonreciprocal dipole-dipole interaction between particles, which is induced by the non-reciprocal interference originating from phase difference between the trapping lasers [49].The proposed system then realizes a 1D tight-binding model with arbitrarily tunable asymmetric hopping amplitudes that have possibly negative signs and long-range dependence.This high controllability allows one to explore the whole parameter region of non-Hermitian systems, thus opening the possibility to fully uncover the potential of non-Hermitian systems.The proposed setup should be contrasted to the previous non-Hermitian platforms where it remains challenging to realize long-range asymmetric and/or negative hopping amplitudes.
To determine the continuum bands and the dynamical phase diagram of the levitated nanoparticle array, we invoke the non-Bloch band theory [74,[77][78][79][80][81], a recently developed powerful tool to investigate models featuring the non-Hermitian skin effect.The non-Bloch band theory allows for calculating the asymptotic eigenvalues under open boundary conditions in the limit of a large system size.This makes contrast to the conventional Bloch band theory, where the band structure reproduces the eigenvalues under periodic boundary conditions.
On the basis of this theoretical framework, we find that the levitated nanoparticle array exhibits rich dynamical phases, including the unconventional critical phase and the dynamically unstable phase.In the former, a remarkable feature is that the non-Hermitian degeneracy of the continuum bands known as the spectral singularity appears over a broad region of the parameters.The key ingredients of the latter are negative interparticle couplings, which were difficult to realize in the existing non-Hermitian platforms.Moreover, the proposed system can naturally realize the long-range hopping amplitudes originating from the dipole-dipole interaction.We show that this long-range nature leads to the unique self-crossing point of the continuum band, which corresponds to the singularity of the generalized Brillouin zone.The distance between the nearest-neighbor particles is d0, and the mass of all the particles is m.All the trapping lasers have the power P and the wavelength λ.We set the phase of the nth trapping laser in the focal plane to be ϕ + n∆ϕ.

II. LEVITATED NANOPARTICLE ARRAY
To be concrete, we consider a 1D array of the trapped levitated nanoparticles (Fig. 1).The particles are equally spaced at the interval d 0 , and all the particles have the mass m.Let λ and P denote the wavelength and the power of all the trapping lasers, respectively.Furthermore, we assume that the motion of the particles along the plane perpendicular to the optical axis is frozen.We note that the distances between the particles are assumed to be much larger than a characteristic length scale in the collective behavior of particles [82][83][84].
The interaction between the two particles arises due to the interference between the scattered electromagnetic field and the trapping laser.Since the scattered field acquires the phase kd 0 during the propagation, the phase difference between the trapping lasers at the positions of the particles leads to the constructive and destructive interference depending on the propagation direction of the scattered field.It is this spatial asymmetry that renders the interparticle coupling nonreciprocal.We note that the effective open boundaries can be realized by arranging the tightly localized particles at the ends of the system (cf.Fig. 6 in Appendix C); the non-Hermitian skin effect then manifests itself as the large oscillation amplitudes close to the boundary regions.Due to the long-range nature of the dipole-dipole interaction, it is in general necessary to incorporate the couplings that reach up to N th nearest neighbor particles.Altogether, the linearized equation of motion of the nth particle along the z axis is given by Here, Ω is an intrinsic mechanical frequency of the particle proportional to √ P , γ is a friction coefficient, and K l and Kl are the coupling strengths given by where G has the dimension of a spring constant and is proportional to P , ∆ϕ is the optical phase difference between the nearest neighbor trapping lasers in the focal plane (Fig. 1), and k 0 (= 2π/λ) is the wavenumber of the trapping laser.We explain the detailed derivation of Eq. ( 1) in Appendix A. Importantly, the optical phase difference gives rise to the nonreciprocal couplings due to the nonzero Kl .Thus, our setup is distinct from the array of levitated nanoparticles proposed in Ref. [85], which has investigated a non-Hermitian transport phenomenon with reciprocal couplings.Furthermore, one can infer from Eq. ( 2) that the couplings are long-range because the dipole-dipole interaction is proportional to the inverse of the distance between the particles.We note that the model can be mathematically mapped to a tight-binding model with gain and loss via a similarity transformation [75].
The continuum bands of non-Hermitian tight-binding models can be obtained by invoking the non-Bloch band theory (cf.Appendix B), which reproduces the eigenvalues under open boundary conditions.Specifically, the continuum bands are calculated from the generalized Brillouin zone spanned by β ≡ e ik for the complex Bloch wavenumber k.We here apply the non-Bloch band theory to the levitated nanoparticle array; throughout this paper, we assume |K N | ̸ = KN .By substituting z n = ψ n e iωt to Eq. ( 1), the real-space eigenequation reads Importantly, the ansatz of Eq. ( 3) can be taken as where β j (= β) is the solution of the characteristic equation given by We note that Eq. ( 5) is an algebraic equation for β of 2N th degrees.The main result of the non-Bloch band theory is that the condition for the generalized Brillouin zone is obtained from the 2N solutions as follows: with The trajectories of β N and β N +1 form the generalized Brillouin zone on the complex plane, which reveals the essential features of non-Hermitian systems (see, e.g., Refs.[86][87][88][89]).Then, we can calculate the continuum bands by combining Eq. ( 5) with the generalized Brillouin zone.

III. DYNAMICAL PHASE DIAGRAM
We start our analysis from the levitated nanoparticle array with the nearest-neighbor interaction, which corresponds to N = 1 in Eq. ( 1); in the following, we assume γ > 2Ω for the sake of concreteness.From Eq. ( 6), the generalized Brillouin zone can be given by the circle with the radius r = K 1 + K1 / K 1 − K1 .The analytical form of the continuum bands reads where θ is a real number.
Since each eigenmode contributes to the dynamics through the factor e −Im(ω±) e iRe(ω±) , we can show the dynamical phase diagram depending on K 1 /m and K1 /m [Fig.2(a)].We here emphasize that the continuum bands discussed here can have direct experimental relevance.Indeed, the theoretical calculation of the eigenvalues in the coupled two levitated nanoparticles has explained well experimentally observed crossing/avoided crossing of the eigenspactra and the appearance of an exceptional point [49].In the blue-shaded regions of Fig. 2(a), all the particles oscillate with the attenuation because Re (ω ± ) ̸ = 0 and Im (ω ± ) > 0 [Fig.2(b)].In contrast, in the red-shaded regions, their motion monotonically vanishes without oscillations because Re (ω ± ) = 0 and Im (ω ± ) > 0 [Fig.2(d)].For these reasons, we term the former (latter) the dynamical phase as the underdamped (overdamped) phase.
Remarkably, we find the broad green-shaded region where the two branches coalesce at Re (ω ± ) = 0 [Fig.2(c)]; this degeneracy is called the spectral singularity.There, we find the crossover behavior where the overdamped behavior eventually sets in after the initial underdamped oscillations; we shall term this intermediate regime as the critical phase.One of its key characteristics is the presence of the particles near the boundaries which are strongly driven by the adjacent trapping lasers inducing the nonreciprocal couplings.It is worthwhile to mention that, with γ < 2Ω, the critical phase appears on the parameter region where the sign of K/m becomes negative [cf.Fig. 7(a) in Appendix C].The transient phenomenon discussed here is supported by the non-Hermitian skin effect, since the critical phase disappears under periodic boundary conditions, as discussed in Ap-pendix D. As shown in Fig. 2(f), the spectral singularity also appears along the green vertical lines in Fig. 2(a).However, one would need the fine-tuning of the parameters in this case as indicated by Figs.2(e)-2(g), where the two continuum bands are recombined across the green line.
Additionally, we also find the dynamically unstable phase as indicated by the gray-shaded region in Fig. 2(a).There, the driving forces from the adjacent trapping lasers give rise to the dynamical instability of the particles where the oscillation amplitudes diverge in the longtime limit, because negative hopping amplitudes cause the force that increasingly keeps away the particles from their equilibrium positions.We expect that nonlinear effects will eventually play a crucial role in this phase, since the amplification is eventually balanced by nonlinear suppression.The dynamically unstable phase discussed here is difficult to realize in previous non-Hermitian systems due to the lack of the ability to implement negative hopping amplitudes.It is worthwhile to mention that, in finite-size systems, the phase boundary between the overdamped and dynamically unstable phases can be slightly modified, as discussed in Appendix C.

IV. NONRECIPROCAL LONG-RANGE INTERACTION
We next investigate how the long-range nature of the couplings can affect the continuum band and the corresponding generalized Brillouin zone; in the following, we neglect the friction for the sake of simplicity.We assume that the interaction reaches up to the nextnearest-neighbor particles, which corresponds to N = 2 in Eq. ( 1).In Fig. 3, we plot the continuum bands with the positive branch of the square root and the corresponding generalized Brillouin zone at different ∆ϕ.The black dashed curves in Figs.3(d)-(f) indicate the conventional Brillouin zone formed by β ≡ e ik (k ∈ R).We note that the parameters considered in these calculations satisfy K 1 /m, K1 /m, K 2 /m, K2 /m ≪ Ω 2 and have been experimentally realized in Ref. [49] in the case of two particles.
One can see from Figs. around the self-crossing points is stronger than that of the other eigenstates away from the self-crossing point.This is because the overlap of the left and right eigenvectors becomes minimum at the self-crossing point.Thereby, such eigenstates exhibit a striking response against perturbations [91], which indicates that the excitation modes around the self-crossing point could be utilized for a highly sensitive sensor.

V. SUMMARY AND DISCUSSION
In summary, we propose and analyze the levitated nanoparticle array as an ideal platform to study new realms of non-Hermitian physics in a highly controlled manner.We show that the system exhibits the unconventional critical phase, where the spectral singularity originating from the non-Hermitian skin effect persists over a broad parameter region.We also point out that the tunable dipole-dipole interaction allows for extremely nonreciprocal hopping amplitudes with possibly negative signs, which result in the dynamical instability.We finally reveal that the long-range nature of the couplings further enriches the non-Hermitian band structures, leading to the cusps of the generalized Brillouin zone and the self-crossing point of the continuum band.
We expect that the levitated nanoparticle array will be an ideal platform to explore rich phenomena induced by nonreciprocal long-range hopping amplitudes thanks to its high tunability.For example, a nontrivial topological phase transition mediated by a topological semimetal phase with exceptional points has been proposed [90,92].Additionally, the number of branches from the self-crossing point of the continuum band increases as a coupling distance between two particles [93].Furthermore, it has been proposed that nonreciprocal long-range couplings give rise to bulk eigenstates which exhibit the crossover from a constant localization length to a systemsize-dependent localization length [94].This previous work has also indicated that, since long-range couplings suppress the non-Hermitian skin effect, the scaling of entanglement entropy can change from an area law to a subextensive law.
The continuum bands of the levitated nanoparticle array can be experimentally studied by measuring the power spectral density.It is thus feasible to directly observe the spectral singularity and the self-crossing points.Meanwhile, to access the strong coupling regime with K/mΩ 2 = O (1-10) considered in Fig. 2, it is necessary to realize a stronger dipole-dipole interaction than the one already realized in Ref. [49].We expect that this should be made possible by increasing size of the particles and using lasers with a shorter wavelength.Also, we note that the friction coefficient can be controlled by changing the pressure of the surrounding gas, while thermal noise can be suppressed by keeping its temperature sufficiently low.We discuss in detail the experimental feasibility in Appendix E.
It is interesting to further explore various aspects and potentials of a levitated nanoparticle array, since this unique platform opens a new avenue of investigating sensing, nonlinearity, and nonequilibrium quantum physics.First, the critical phase can be potentially applicable for enhanced sensing because the strong nonorthogonality associated with the spectral singularity which one can observe in the array of a few particles (cf.see Appendix C) is known to trigger the singular sensitivity to perturbation [95][96][97].Second, it is crucial to reveal the competi- tion between nonreciprocity and nonlinearity [98], where a levitated nanoparticle array is expected to exhibit rich collective phenomena such as synchronization.Third, an extension of the present analysis to quantum regimes is an intriguing open problem.For instance, it merits further study to understand how the nonreciprocal couplings affect the entanglement between particles [40], how the dynamical instability of the levitated nanoparticle array can be utilized to squeeze the mechanical mode of a levitated nanoparticle [99], and whether or not the nonorthogonality is helpful for quantum metrology [39].

Appendix A: Equation of motion
We derive the equation of motion of a single levitated nanoparticle, coupled two levitated nanoparticles, and an array of levitated nanoparticles as shown in Figs.4(a), 4(b), and 4(c), respectively.We assume that the motion of the particles along the plane perpendicular to the optical axis is frozen.

Single levitated nanoparticle
We first focus on the motion of a single levitated nanoparticle along the z axis [Fig.4(a)].Let P and λ denote the power and the wavelength of the trapping laser, respectively.In this system, the subwavelengthsized particle with the refractive index n is surrounded by a medium with the refractive index n ′ , and it can be regarded as a point dipole with the electric dipole moment given by p (r, t) = αE (r, t) . (A1) Here, E (r, t) expresses the electric field of the trapping laser, and α is the polarizability of the particle.The latter can explicitly be written as where c is the speed of light in vacuum, V is the volume of the particle, and n r ≡ n/n ′ is the relative refractive index of the particle.The electric field then acts on the point dipole, and the particle feels the gradient force [100] given by where ⟨• • • ⟩ T means time average, and I (r) is the spatial profile of the intensity of the trapping laser.We note that the trapping laser is described by the Gaussian beam, of which the electric field is given by where Here, k 0 (= 2π/λ) is the wavenumber of the trapping laser, w 0 is the beam waist, and z 0 = πw 2 0 /λ is the Rayleigh length.The combination of Eqs.(A3) and (A4) leads to the z component of the gradient force given by In the vicinity of the focal plane of the trapping laser, the linearized equation of the motion of the levitated nanoparticle is obtained as where and γ is a friction coefficient.We remark that the linearized motion of the particle exhibits the harmonic oscillation because the gradient force plays a role as the restoring force.

Coupled two levitated nanoparticles
We next consider the equation of motion of the coupled two levitated nanoparticles [Fig.4(b)].We assume that the polarization direction of the particles points to the y axis.Then, in addition to the gradient force, both particles feel the interaction caused by the dipole radiation, which is called the optical binding force [101,102].
First of all, we shall explain the mechanism of the optical binding force based on Ref. [49].Importantly, there exist two contributions to the optical binding force.The first contribution results from the combination of the electric field scattered from one particle to another particle and the electric dipole induced by the trapping laser.The second contribution originates from the acting of the trapping laser on the electric dipole induced by the scattered field.Specifically, for particle 1, the sum of these contributions leads to the form of the optical binding force as follows: Here, E 2→1 sca (r 1 , r 2 , t) is the electric field scattered from particle 2 to particle 1, and its form is written as where G (r) called the Green's tensor is the electric field propagator between the two dipoles [103], and it is given by and ∇ r1 × E (r 1 , t) = 0, we can rewrite Eq. (A9) to a brief form given by In the following, we calculate the z component of the optical binding force.In the far-field regime with k 0 d 0 ≫ 1, the dominant contribution from the Green's tensor to the optical binding force is the term proportional to 1/r.Furthermore, in the vicinity of the focal plane, the z dependence of the electric field is approximated as for j = 1, 2, where ϕ j expresses the optical phase at the focal plane.Substituting Eqs.(A11) and (A14) into Eq.(A13), we can get the optical binding force along the z axis as follows: where ∆ϕ ≡ ϕ 1 − ϕ 2 .We can similarly obtain the optical binding force for particle 2 as follows: Therefore, the explicit form of the optical binding force along the z axis is derived as Remarkably, the interaction between these two particles becomes nonreciprocal because The key ingredient of this nonreciprocity is the interference between the trapping laser and the scattered field.Let Φ j (j = 1, 2) denote the optical phase of the trapping laser at the position of the particle.The interference depends on the local phase difference ∆Φ ≡ Φ 1 − Φ 2 and the phase accumulation kd 0 which the scattered field acquires during the propagation.Specifically, while the contribution of the interference is kd 0 − ∆Φ within the propagation of the scattered field from particle 1 to particle 2, it becomes kd 0 + ∆Φ in the opposite case.As a result, the interaction originating from the interference becomes spatially asymmetric.
We can now obtain the linearized equation of motion of the coupled levitated nanoparticles as follows: where Ω is the intrinsic mechanical frequency given by Eq. (A1).The coupling constants are given by and (A20)

Levitated nanoparticle array
We finally study the arrangement of the multiple levitated nanoparticles at equal interval d 0 [Fig.4(c)].In this system, the dipole-dipole interaction among the several particles arises from the multiple scattering of the trapping lasers.Nevertheless, the dominant contribution to the dynamics of the system comes from the interaction between two particles.Thus, we neglect higher-order scattering processes.This corresponds to approximating the optical binding force up to O |p| 2 .
We shall explain how one can derive the equation of motion of the levitated nanoparticle array.We set the optical phase of the nth trapping laser in the focal plane to be ϕ + n∆ϕ.For the nth and n + lth particles, the phase difference between the trapping lasers is l∆ϕ, and the distance between the particles is ld 0 .Hence, the interaction between these particles can be obtained by the same procedure as explained above.Due to the longrange nature of the dipole-dipole interaction, it is necessary to incorporate the couplings that reach up to N th neighbor particles.Then, the equation of motion of the system is written as where the constant G is given in Eq. (A20).

Appendix B: Non-Bloch band theory
We describe the physical meaning of the non-Bloch band theory for 1D non-Hermitian tight-binding models.First of all, we show that the condition for the generalized Brillouin zone can be interpreted as the condition that the "plane waves" form the standing wave, by using the specific model.We next prove that the recombination of the "plane waves" forming the standing wave occurs at the cusps of the generalized Brillouin zone.

Condition for the generalized Brillouin zone
We consider the 1D non-Hermitian tight-binding model with asymmetric hopping amplitudes, the Hamiltonian of which reads where all the parameter are real numbers.The Schrödinger equation can be then written in the form of the real-space eigenequation given by where ψ n means the amplitude of the state at the site n.A general difference theory allows us to take the form of the linear combination which corresponds to the plane-wave expansion, as the ansatz of Eq. (A2), due to the spatial periodicity.Here, β j (= β) is the solution of the characteristic equation written as which is a quadratic equation for β.We note that it is necessary to combine the open boundary conditions obtained from Eqs. (B3) and Eq.(B4) to get the asymptotic set of the energy eigenvalues in the thermodynamic limit.Although the calculation is cumbersome for a large system size, the non-Bloch band theory allows us to avoid the procedure of calculating the continuum band [77,79].It has been shown that the values of β are restricted to lie on the closed curve so that the wavefunction satisfies the open boundary conditions.The closed curve, called the generalized Brillouin zone, is then formed by β = e ik for the Bloch wavenumber k.It is important that, for the four solutions of Eq. (B4), the condition for the generalized Brillouin zone reads with We note that the trajectories of β 2 and β 3 satisfying Eq. (B5) form the generalized Brillouin zone.Finally, the combination of Eqs.(B4) and (B5) gives the continuum band of the system.We exemplify the continuum band and the generalized Brillouin zone with the specific parameters as shown in Figs.5(a) and 5(b).The Bloch wavenumber obtained by Eq. (B5) becomes complex numbers, which indicates that the bulk eigenstates are localized at boundaries of the system due to the non-Hermitian skin effect.Meanwhile, the generalized Brillouin zone forms a unit circle, when the system becomes a Hermitian system: t +,1 = t −,1 and t +,2 = t −,2 .This means that the Bloch wavenumber takes real values in Hermitian tight-binding systems, which is consistent with the result of the conventional Bloch band theory.
It is remarkable that the condition for the generalized Brillouin zone can be interpreted in the viewpoint of the non-Hermitian skin effect.The physical meaning of Eq. (B5) is that the localization lengths of the "plane waves" corresponding to β 2 and β 3 match each other, which leads to the formation of the standing wave by the interference of the "plane waves".Furthermore, it is intriguing that the asymptotic energy eigenvalues of the open chain in the thermodynamic limit do not depend on any boundary conditions, since Eqs.(B4) and (B5) are independent of boundary conditions of the open chain.

Cusps of the generalized Brillouin zone
We nest explain the appearance mechanism of the cusp of the generalized Brillouin zone, at which it becomes indifferentiable.To this end, we investigate what happens if we impose to the system, for some i and j among the four solutions of Eq. (B4).We then show several sets of β i and Thus, the recombination of the "plane waves" forming the standing wave occurs at the cusps of the generalized Brillouin zone and the corresponding self-crossing points of the continuum band.
Appendix C: Array of finite levitated nanoparticles We investigate the array of a finite number of levitated nanoparticles as shown in Fig. 6.We assume that this system includes only the interaction between the nearestneighbor particles.Furthermore, at both boundaries of the system, we arrange the deeply trapped levitated nanoparticles, of which the motion in all the directions is frozen, and we set the optical phase of the trapping laser in the focal plane at the left and right boundaries to be ϕ − ∆ϕ and ϕ + (L + 1) ∆ϕ, respectively.In this case, these particles are coupled with the system so that they impose the fixed end boundary conditions on the system.

Eigenvalue
We first show a way to calculate the eigenvalues of the system described by with the fixed boundary condition given by z 0 = z L+1 = 0.In this equation, the coupling constants K and K are given by Eq. (A22).In the following, we suppose |K| ̸ = K .By assuming z n = ψ n e iωt , Eq. ( C1) is rewritten into From a general theory of a difference equation, we can take as an ansatz of Eq. (C2).Here, β j (= β) is the solution of the characteristic equation given by We note that Eq. (C4) is a quadratic equation for β.
The boundary conditions, ψ 0 = ψ L+1 = 0, tell us the condition that the combination coefficients ϕ (1) and ϕ (2)  take nonzero values, and it is written as Then, we can obtain the explicit form of β 1 and β 2 from Eqs. (C4) and (C5).When K + K K − K > 0, the Vieta's formula of Eq. (C4) gives where Hence, the eigenvalues of the system can be calculated as (C9) Similarly, when K + K K − K < 0, we obtain the form of β 1 and β 2 as follows: Here, and θ l is given by Eq. (C8).The eigenvalue of the system in this case is written as (C12)

Finite-size effect
In our work, we investigate the levitated nanoparticle array in the thermodynamic limit, L → ∞, and obtain the dynamical phase diagram and the continuum bands of the levitated nanoparticle array.We show the dynamical phase diagram with γ < 2Ω and γ > 2Ω in Figs.7(a) and 7(b), respectively.We note that the critical phase extends over only a narrow region in Fig. 7(a), while it clearly possesses a broad region in Fig. 7(b).We here discuss how the finite-size effects can affect the band structures in the critical phase and the boundaries between the underdamped and dynamically unstable phases.To do so, we recall that the continuum bands of the infinite-size system is given by ω± (C13) First, we consider the parameters indicated by the orange star in Fig. 7(a) and show the corresponding continuum bands ω± and the eigenvalues ω > l,± in Fig. 7(c).One can see that there is no degeneracy between ω > l,+ and ω > l,− in a strict sense, simply because θ l takes only the discrete values determined by Eq. (C9).Nevertheless, we emphasize that the eigenvectors around the spectral singularity, where ω+ touches ω− , still exhibit the strong nonorthogonality which is a hallmark of the non-Hermitian degeneracy.This result is the same as in the case of the parameter indicated by the white star in Fig. 7(b).Indeed, the eigenvalues in Fig. 7(d) are similar to those in Fig. 7(c).Next, we consider the parameters indicated by the white star in Fig. 7(a) and show the corresponding continuum band ω+ and the eigenvalues ω < l,+ in Fig. 7(c).It is found that the finite-size system does not exhibit the dynamical instability because the imaginary parts of all the discrete eigenvalues become positive.In this sense, the finitesize effect can slightly modify the boundary between the dynamically unstable phase and the other phases.One can also infer from Fig. 7(c) that a larger system size would be favorable to observe the dynamical instability.Nevertheless, we emphasize that the phase diagram of a finite-size system still remains qualitatively the same as in the result obtained in the thermodynamic limit.

Appendix D: Periodic boundary condition
We consider the levitated nanoparticle array described by Eq. (A14) with N = 1 and study the dynamical phase diagram under periodic boundary conditions to compare the result obtained in Fig. 2(a).We note that the characteristic equation of the system reads After replacing β by e ik for the real Bloch wavenumber k, we can obtain D2) which reproduces the eigenvalues of the system with periodic boundary conditions.Let us suppose that γ > 2Ω for the sake of concreteness.The condition for the appearance of the spectral singularity can be then written as The blue and gray-shaded regions are the underdamped and dynamically unstable phases, respectively.The spectral singularity appears on the green line.We set the parameters to be Ω = 1 and γ = 5.
since it appears when k = 0. Remarkably, this indicates that the critical phase obtained in Fig. 2(a) disappears in the system with periodic boundary conditions.Said differently, the appearance of the spectral singularity requires the fine-tuning of the system parameters.Indeed, one can see from the dynamical phase diagram shown in Fig. 8 that only the underdamped and dynamically unstable phases appear in this case.Thus, modifying open boundaries to periodic boundaries drastically changes the behavior of the levitated nanoparticle array.

Appendix E: Experimental feasibility
We investigate the dynamical phase diagram of the levitated nanoparticle array with K/m, K/m ∈ −10Ω 2 , 10Ω 2 and γ = 5Ω in our work.We here discuss the experimental values to access the assumed parameter region in our setup.We note that the previous experiment realizing the coupled two levitated nanoparticles has achieved K/m, K/m ∈ −0.1Ω 2 , 0.1Ω 2 and γ ≃ 0.03Ω [49], where the power and wavelength of the trapping laser, and the gas pressure are set to be P = 400 mW, λ = 1064 nm, and P gas = 1.5 mbar, respectively.Furthermore, the previous work has utilized silica nanoparticles, the radius, the polarizability, and the mass of which are set to be r = 105 nm, α = 3.48 × 10 −32 F • m 2 , and m = 1.07 × 10 −17 kg, respectively.First, we explain a possible way to amplify the dipoledipole couplings K/m, K/m compared to Ω 2 .To this end, we consider the ration between G/m and Ω 2 , which can be obtained from Eqs. (A8) and (A20) as follows: It is important that one can achieve K/m, K/m ∈ −10Ω 2 , 10Ω 2 by using the trapping laser with the half wavelength and making the radius of the particle 2.5 times larger than the one used in the previous experiment; we recall that the larger particle results in higher polarizability.Thus, we assume that the wavelength of the trapping laser is set to be λ = 532 nm, and the radius, the polarizability and the mass of the particles are set to be r = 260 nm, α = 5.28 × 10 −31 F • m 2 , and m = 1.62 × 10 −16 kg, respectively, and the laser power is the same as in the previous experiment.
We expect that our model captures qualitative features of the motion of the levitated nanoparticle even if its size is comparable to the wavelength of the trapping laser, because the dominant contribution to the forces which the particle feels still originates from the dipole moments in the particle.Meanwhile, to develop a quantitatively accurate theory, it should be necessary to modify the description to go beyond the dipolar approximation.Specifically, one should take into account the contributions from the higher-order moments (e.g., quadrupole moments).
Second, we explain how to realize the friction coefficient γ = 5Ω.We note that, in the low-vacuum regime, the damping rate of a levitated nanoparticle is proportional to the square of the radius of the particle and the gas pressure [104].With the above parameter set, Ω approximately becomes twice larger than the value in the previous experiment, and it is possible to realize γ = 5Ω by setting the gas pressure to be P gas = 80 mbar.Thus, we expect that the levitated nanoparticle array with the above experimental value set should allow one to realize the parameters assumed in our manuscript and to observe the predicted dynamical phases accordingly.

FIG. 1 .
FIG. 1. Schematic figure of the levitated nanoparticle array.The distance between the nearest-neighbor particles is d0, and the mass of all the particles is m.All the trapping lasers have the power P and the wavelength λ.We set the phase of the nth trapping laser in the focal plane to be ϕ + n∆ϕ.
Figures 2(b)-2(d) and 2(e)-2(g) plot the evolutions of the continuum bands along the black and white arrows indicated in Fig. 2(a), respectively.

FIG. 2 .
FIG. 2. Dynamical phase diagram and continuum bands of the levitated nanoparticle array.(a) Dynamical phase diagram exhibiting the underdamped, critical, overdamped, and dynamically unstable phases shown in the blue, green, red, and gray-shaded regions, respectively.The spectral singularity (SS) appears in the green-shaded region and on the green lines.|K1| = K1 is satisfied on the black dashed lines.We set the parameters to be Ω = 1 and γ = 5. (b-g) Evolutions of the continuum bands along the arrows in panel (a).The magenta and cyan express ω+ and ω−, respectively.The numerical values in each panel specify K1, K1 .

FIG. 4 .
FIG. 4. Schematic figure of a single levitated nanoparticle, coupled two levitated nanoparticles, and an array of leviated nanoparticles.In panels (a), (b), and (c), the mass of the particle is m, and the trapping lasers have the power P and the wavelength λ.In panels (b) and (c), the distance between the particles is d0.In panel (c), we set the phase of the nth trapping laser in the focal plane to be ϕ + n∆ϕ.
ACKNOWLEDGMENTS K.Y. was supported by JSPS KAKENHI through Grant No. JP21J01409.Y.A. acknowledges support from the Japan Society for the Promotion of Science through Grant No. JP19K23424 and from JST FOREST Program (Grant No. JPMJFR222U, Japan).

FIG. 6 .
FIG.6.Schematic figure of the array of the L levitated nanoparticles.The distance between the nearest-neighbor particles is d0, and the mass of all the particles is m.The trapping lasers have the power P and the wavelength λ.We set the phase of the nth trapping laser in the focal plane to be ϕ + n∆ϕ.At both boundaries of the system, we arrange the two levitated nanoparticles frozen in the motion in all the directions.The trapping lasers at the left and right boundaries have the optical phase ϕ − ∆ϕ and ϕ + L∆ϕ in the focal plane, respectively.

FIG. 7 .
FIG. 7. Dynamical phase diagram, continuum bands, and eigenvalues of the levitated nanoparticle array.(a,b) Dynamical phase diagrams with Ω = 3 and γ = 2 and Ω = 1 and γ = 2. respectively.Panel (b) is the same as Fig. 2(a).The blue, green, red, and gray-shaded regions indicate the underdamped, critical, overdamped, and dynamically unstable phase, respectively.The spectral singularity (SS) appears in the green-shaded region and on the green lines.(c) and (d) Continuum bands ω+ (magenta) and ω− (cyan), and eigenvalues ω > l,+ (red) and ω > l,− (blue).(e) Continuum band ω+ (black), and eigenvalues ω < l,+ (red).We set L = 8 in panels (c-e) and choose the values of K, K at the orange star for panel (c) and at the white and black stars for panels (d) and (e), respectively.

FIG. 8 .
FIG. 8. Dynamical phase diagram of the levitated nanoparticle array obtained from the conventional Bloch band theory.The blue and gray-shaded regions are the underdamped and dynamically unstable phases, respectively.The spectral singularity appears on the green line.We set the parameters to be Ω = 1 and γ = 5.