Anomalous in-gap edge states in two-dimensional pseudospin-1 Dirac insulators

Quantum materials that host a flat band, such as pseudospin-1 lattices and magic-angle twisted bilayer graphene, can exhibit drastically new physical phenomena including unconventional superconductivity, orbital ferromagnetism, and Chern insulating behaviors. We report a surprising class of electronic in-gap edge states in pseudospin-1 materials without the conventional need of band-inversion topological phase transitions or introducing magnetism via an external magnetic type of interactions. In particular, we find that, in two-dimensional gapped (insulating) Dirac systems of massive spin-1 quasiparticles, in-gap edge modes can emerge through only an {\em electrostatic potential} applied to a finite domain. Associated with these unconventional edge modes are spontaneous formation of pronounced domain-wall spin textures, which exhibit the feature of out-of-plane spin-angular momentum locking on both sides of the domain boundary and are quite robust against boundary deformations and impurities despite a lack of an explicit topological origin. The in-gap modes are formally three-component evanescent wave solutions, akin to the Jackiw-Rebbi type of bound states. Such modes belong to a distinct class due to the following physical reasons: three-component spinor wave function, unusual boundary conditions, and a shifted flat band induced by the external scalar potential. Not only is the finding of fundamental importance, but it also paves the way for generating highly controllable in-gap edge states with emergent spin textures using the traditional semiconductor gate technology. Results are validated using analytic calculations of a continuum Dirac-Weyl model and tight-binding simulations of realistic materials through characterizations of local density of state spectra and resonant tunneling conductance.


I. INTRODUCTION
The physics of quantum materials hosting a flat band, such as the magic-angle twisted bilayer graphene, has become a forefront area of research. These materials can generate surprising physical phenomena such as unconventional superconductivity [1,2], orbital ferromagnetism [3,4], and the Chern insulating behavior with topological edge states. The purpose of this paper is to report the surprising emergence of a class of in-gap edge states in two-dimensional Dirac/Weyl pseudospin-1 materials, which cannot be fit into any of the known scenarios for producing such states. The uncovered states, at their birth, exhibit topologically nontrivial domain-wall like pseudospin textures.
In modern physics, the emergence of low-dissipation or dissipationless topological surface or edge states in condensed matter systems is a fascinating phenomenon [5][6][7], as exemplified by topological insulators (TIs) [8][9][10][11][12][13][14][15][16]. A TI has a bulk band gap so its interior is insulating but there are gapless surface states within the bulk band gap, which are protected by the time-reversal symmetry that renders the states robust against backscattering from nonmagnetic impurities. These topologically protected surface or edge states possess a perfect spin-momentum locking characterized by the invariance of spin orientation with respect to the direction of the momentum. Quite recently, high-order TIs hosting, e.g., robust in- * Ying-Cheng.Lai@asu.edu gap excitations of zero-dimensional corner modes have been uncovered [17][18][19]. Topological states of matter, in addition to their importance to fundamental physics, have potential applications in electronics and spintronics [20]. For electronic systems, current understanding of the physical mechanisms behind the topological edge states requires a discontinuous change in the associated bulk topological invariants across the interface/edge rendered by, e.g., a strong external magnetic field in a twodimensional electron gas [6], band inversion driven by spin-orbit coupling [8,21,22], introduction of ferromagnetism in topological insulators [23], presetting domain walls in gapped Dirac materials [24][25][26][27], stacking order in layered two-dimensional materials [28], and particular spatial crystalline symmetries [29].
Pseudospin-1 type of low-energy excitations beyond the Dirac-Weyl-Majorana paradigm have recently been realized in electronic lattice systems [30][31][32][33][34][35]. In a broader context, two-dimensional massive spin-1 bulk excitations can arise in classical nonlinear physical systems such as rotating shallow water in a horizontally unbounded plane [36] and the wave system of magnetoplasmons, where the corresponding Hamiltonian representations [37] can be effectively reduced to the Dirac-like equation for spin-1 particles. Applying the sign-changing Dirac mass scenario to the systems leads to an extension of the Jackiw-Rebbi mechanism that serves to ascertain the topological origin of, e.g., the equatorial waves [36], as well as rich topological phenomena in bosonic and classical systems [38][39][40][41][42].
The subject of our study is pseudospin-1 relativistic quantum systems described by the generalized Dirac-Weyl equation which are fundamentally linear. Specifically, low-energy excitations in condensed matter systems such as graphene [43] and topological insulators [8][9][10][11][12][13][14][15][16] and in analogous physical systems of molecules, cold atoms, cavity polaritons, light and even mechanical waves in judiciously designed lattices [38][39][40][41][42] are described by the Dirac-Weyl equation. In those circumstances, if the corresponding quasiparticles are massless, the energy band structure contains a pair of Dirac cones characteristic of the relativistic energy-momentum dispersion relation. A finite mass leads to a band gap, giving rise to unconventional topological phases in Dirac material systems [44] with unusual physical properties associated with tunneling, confinement and transport, which have no analogies in quantum systems described by the Schrödinger equation. Among those, the physics of edge states and robust in-gap excitations are of fundamental interest. Jackiw and Rebbi [45] predicted a surprising zero-energy bound state solution of the Dirac equation in the presence of a kink-shaped mass profile that generates a domain wall separating regions with sign-changing Dirac mass. The realization in polyacetylene [46,47] and the theoretical studies of narrow-gap semiconductors [48,49] led to the discovery of the phenomenon of band-gap inversion enabling topologically protected conducting interface states and localized sub-gap excitations in TIs [8][9][10][11][12][13][14][15][16][17][18][19]. In the description based on the massive Dirac equation, band-gap inversion is equivalent to a sign change in the mass. The topological edge states give rise to appealing physical properties and phenomena such as robust low-power-dissipation wave transport [26], electrically tunable magnetism [50], and quasiparticles analogous to elementary fermionic particles in high-energy physics [51].
Our main finding is that, in pseudospin-1 systems with an energy gap, a surprising class of in-gap edge bound states can arise without band or mass inversion based domain-walls that separate the regions with different kinds of bulk band topology and any external magnetic interaction, but these states are remarkably robust against geometric deformations and impurities. In fact, they are generated through only a local electrostatic potential barrier of the repulsive type in the underlying insulating spin-1 systems. We uncover a number of remarkable, quite unusual spectral properties of these modes. Unlike the topological edge states previously discovered and studied, the states reported here require no established topological restrictions such as interfacing domains/systems of different bulk topological invariants and any particular type of discrete symmetries. In fact, through self-inducing topological spin textures, the uncovered states possess the degree of robustness enjoyed by conventional topological states but they belong to a distinct class due to the following physical reasons: three-component spinor wave function, unusual boundary conditions, and a shifted flat band induced by the external electrical potential. Experimentally, these states can be generated readily through routine electro-static gating within the same material (or within a single device), rendering them promising in applications, e.g., a gate-controlled spin-1 Dirac electron transistor of high on/off ratio.
Schematic illustration of the system setting and main finding. (A) A side view of the setting leading to ingap edge modes without magnetism and the conventional band inversion topological phase transition, where a gapped two-dimensional system (thick gray line) hosting Dirac-like low-energy excitations of massive spin-1 is subject to a locally applied electrostatic potential. The energy band diagram in the absence of the potential is shown on left side of the bottom panel [below (A)]. (B) In-gap edge bound modes with an emergent domain-wall like spin ordering/texture (top panel) arise in the presence of a repulsive type of potential, which defines an antidot profile as shown in the bottom panel. The criterion for the stable emergence of the in-gap states is |Vg − ∆| ∆/2. Figure 1A presents a schematic illustration of the system setting, whose effective Hamiltonian is H ef f = v FŜ ·p + ∆Ŝ z + U (r), where the first term describes the bulk low-energy excitation of a massive spin-1 particle with quasi-momentump = (p x , p y ), the second term represents the generalization of the Dirac mass withŜ z being a component of the spin-1 matrix vectorŜ, and the last term is the locally applied electrostatic potential of height V g which defines a closed interface at the boundary. As we will establish, this magnetism-free configuration permits in-gap edge states, and the states with higher angular momenta possess highly organized domain-wall like spin textures, as illustrated in Fig. 1B. In general, for the in-gap states to emerge and be stable, the perturbation in the form of the applied gate potential V g cannot be negligibly small in comparison with the pristine band gap ∆. Neither can the perturbation be too large to result in a substantially reduced effective band gap size. In fact, the inequality V g < 2∆ is required and the reduced band gap (2∆ − V g ) should be comparable to the pristine one. In terms of ∆ and V g as defined in Fig. 1, the criterion for the stable emergence of the in-gap states is |V g − ∆| ∆/2. The case shown in Fig. 1 is for V g ∆.

B. Emergence of in-gap edge states
For a circular domain of radius R, the electrostatic potential is given by U (r) = V g Θ(R − r), where Θ is the Heaviside function. The system as governed by H ef f ψ = Eψ can be solved analytically in the polar coordinates r = (r, θ) to yield closed-form solutions of the form where µ = I, O labels the inner and outer regions as de- are the Bessel and the Hankel functions of the first kind with j being the integer angular momentum quantum number. As explained in Appendices A and B, the in-gap modes uncovered take the form of a three-component evanescent edge state. In comparison with known edge states, either topological or nontopological, the states uncovered here belong to a distinct class due to the following physical reasons: threecomponent spinor wave function, unusual boundary conditions, and a shifted flat band induced by the external scalar potential. Particularly, for |j| 1, we calculate the eigenenergy E ≈ V g /2 and the resulting spin textures where with ξ = (V g + 2∆)r/2 v F . Concretely, for a representative parameter setting, e.g., V g = ∆ = 6 v F /R, we calculate the resulting energy spectra as a function of the angular momentum quantum number j, as shown in Fig. 2A. We see that additional bounded eigenstates arise in the gap, i.e., those in the shaded area in Fig. 2A. The striking feature is that these states emerge for E µ < ∆, where the system is an insulator. In this case, without any change in the band topology (e.g., due to band inversion), conventional understanding of TIs stipulates that such states are impossible.
A feature of the spin textures is worth mentioning. If we calculate the topological number defined as where n = S/|S|, we get N = −sign(j)/2, signifying vortex-like spin textures that can arise from in-gap excitations of meron-like skyrmions [52]. Similar features have been predicted in chiral p-wave superconductors [53,54] that have the same symmetry class as the spin-1 Dirac Hamiltonian studied in this paper [Eq. (A1)].
To gain further insights, we characterize the energy spectra using two experimentally relevant quantities: the local density of state (LDOS) and spin-LDOS defined as D(E, r) = ν ν|ν δ(E − E ν ) and D s (E, r) = ν ν|S z |ν δ(E − E ν ), respectively, where ν is the eigenstate label. As shown in Fig. 2C, the in-gap modes are localized at the boundary and exhibit distinct domain wall spin textures, where the energy broadening effect (e.g., caused by measurement) has been taken into account by approximating the delta function as Figure 2D shows the spatial distributions of the corresponding wave density and spin texture for a representative state (indicated by the red arrow in Fig. 2A). A calculation of the associated spin projection S z versus j in the inner and outer regions reveals that the domain wall spin ordering is more pronounced for states with higher angular momenta, as shown in Fig. 2B. Associated with the strengthening of the spin ordering, the energy flow tends to decrease, as revealed by a nearly dispersionless dependence of the energy level on the angular momentum quantum number, as shown in Fig. 2A. Figure 2B demonstrates the emergence of spin-angular momentum locking that depends on the side of the interface in which the state is located, suggesting that these states are robust.

C. Robustness
The robustness of the QSH and QAH edge states are known to be protected by the presetting discontinuous change in the associated bulk topological invariants across the interface, such as the Z 2 index and Chern number, all requiring some sort of magnetic interaction. However, for the edge modes demonstrated in Fig. 2, there is no such a priori topological origin/restriction. The question is whether the modes are protected or stable against irregular perturbations. To address this question, we consider a general type of perturbation: geometric deformation of the potential domain. A significant challenge is to obtain accurate eigensolutions of the massive spin-1 Dirac equation as, with an irregular domain, analytic solutions are no longer possible. We have developed an accurate and efficient numerical method to find solutions for arbitrarily shaped domain interfaces (Appendix C). As an illustration, we create deformed domains via the superformula in botany that can generate a great diversity of natural shapes with only a few parameters [55]. Figure 3A shows, for thirteen deformed boundary shapes (insets), the corresponding energy spectra resolved by the total density of states (DOS). With respect to the eigenstates of the circular geometry, there are considerable shifts (up or down) in some eigenenergies of the strongly deformed domains but, importantly, there are stable states with virtually no changes in their energies in spite of the severe deformations.
To ascertain the nontrivial feature of the in-gap states, we examine the associated spin properties. In particular, we introduce an effective exchange energy penalty: to identify a domain wall like spin ordering structure between the inner and outer regions. It can be seen from Fig. 3B that the stable modes insensitive to deformation attain large energy penalties, a strong indication of the emergence of domain wall spin ordering, while the modes with small values of E w are sensitive to deformations. Figure 3C shows the real-space wave density and the corresponding spin texture patterns of three representative states as indicated in Fig. 3B. The wave density topography associated with the strong domain wall spin texture is mainly contributed by the high angular momentum states (those with distinctly more angular nodes -c.f., middle panel of Fig. 3C). This agrees with the prediction of the continuum theory that a nearly perfect out-of-plane spin-angular momentum locking should emerge for the high orbital angular momentum states, as shown in Fig. 2B, providing the physical reason for the robustness. (Intuitively, this behavior can be understood that a faster spinning egg is able to stand upright in a more stable manner.) The unambiguous signature of spin-angular momentum locking can greatly circumvent mode coupling due to backscattering caused by the deformation. For those modes, the conventionally anticipated level repulsion/shifting effect due to geometric deformation is greatly suppressed, an unequivocal indication that the modes with spin-angular momentum locking are robust with self-induced protection.
As in most studies of TIs [8,56], we have employed a sharp potential boundary to demonstrate the findings. However, by performing calculations using a finite difference method for realistic and smoothly varying potential profiles, we find that the topological states as exemplified in Figs. 2 and 3 persist (Appendix C). We also find that these states can tolerate strong disorders. The in-gap excitations predicted have the striking physical properties of dispersionless spectral flow and spontaneous domain wall spin ordering. They manifest themselves as distinct real-space topographies of LDOS and spin-LDOS, which can be experimentally mapped out using the low-temperature scanning tunneling spectroscopy technique [57,58]. With advances in Dirac materials in recent years, realizing the spin-1 generalization of ordinary Dirac/Weyl fermions in the form of low-energy collective states or quasiparticles is experimentally possible in condensed matter systems [31,33,34,59,60], photonic crystals [60,61], and even classical systems [37].
Our theoretical prediction is general for gapped systems of massive spin-1 particles subject to an electrostatic potential applied to a finite domain. The bandgap associated Dirac-like mass generation can be implemented in alternative ways. For example, for a twodimensional lattice with three sublattices [33,[59][60][61][62], such as a Lieb or a dice lattice, the generalized mass term can be induced via a staggered sublattice potential that breaks the inversion symmetry, which is an extension of the standard Dirac mass term in, e.g., graphene. As a way of example, we consider the case of a dice lattice model as illustrated in Fig. 4(a), which are relevant to emerging 2D Dirac materials such as transition metal dichalcogenide/dihalide monolayers [63], monolayer Mg 2 C (MXene) [64], decorated graphene [65] etc. Its tight-binding Hamiltonian in real space is given by where c † νi (c νi ) with ν = A, B, C are creation (annihilation) operators of the localized states |νi at site i belonging to the sublattice ν, i, j denotes pairs of nearestneighbor sites with the tunneling strength (hopping energy) of t. The last term represents a staggered sublattice potential that is responsible for the Dirac-type mass based gap opening. In the absence of any external field, we obtain the bulk energy band structure and corresponding LDOS spectra, as shown in Fig. 4(b). We see that, near the K point, the system behaves as a band insulator hosting Dirac-like quasiparticles of massive spin-1. Notably, the flat band leads to a sharp peak in the LDOS, but has a vanishing group velocity as well as a vanishing out-of-plane pseudospin polarization/orientation [c.f. right panel of Fig. 4 An electrostatic potential of height V 0 /t is locally applied to a small region of an undoped dice lattice sheet to realize the gate controlled quantum dot structures. Concretely, for ∆/t = 0.439 and V 0 /t = ∆/t(< 2∆/t), we calculate the sLDOS measured at the boundary of the gated region for three different domain shapes with a characteristic size parameter R = 5nm as depicted in insets of Figs. 5(a-c). The results are displayed by red curves, while those for the (ungated) case of V 0 /t = 0 (black curves) are also shown for comparison. Signified by dramatic changes in the sLDOS spectra with large amplitudes, a number of in-gap states emerges. As displayed in the middle panel of Fig. 5, they are highly localized edge modes. This result agrees with that obtained from the analytic continuum spin-1 Dirac model in Sec. II.
We also consider a lead-contacted dice lattice flake with a circular gate-defined quantum dot as schematically illustrated in the top panel of Fig. 5(d) for a possible experimental detection via transport measurements. One typical simulation result is given in bottom panel of Fig. 5(d).
Remarkably, the emerging in-gap modes acting as "doorway" states can actuate resonant tunneling through the device with large conductance. Because of the east of realizing control with an electrostatic gate potential, the setup can act as a novel quantum switch or transistor of high on/off ratio with spin-1 Dirac electrons.
Alternatively, associated with triple point semimetals of bulk massless spin-1 excitations described by a three-band extension of the Weyl Hamiltonian [34], i.e., H 3 ∝ k x S x + k y S y + k z S z , a thin film structure of thickness L in the z direction can host the two-dimensional spin-1 quasiparticles with an analogous finite mass ∝ π/L due to the confinement effect. This provides another potential experimental platform. In addition, the massive spin-1 physics turns out to be accessible in a dimerized quantum magnet [31] and is even relevant to classical systems of two-dimensional magnetoplasmon [37], where the mass term is induced by an applied magnetic field.
We have also solved a gapped Dice lattice in a semiinfinite geometry in the presence or absence of a locally applied gate potential, which represents a trivial bulk band insulator with low energy, massive, pseudospin-1 excitations. For comparison, we have also included the known case of gapped graphene. The results are shown in Fig. 6. It can be seen that, in contrast to the well studied graphene case [(c) and (d)], in the Dice lattice with massive pseudospin-1 quasiparticles, the in-gap states emerge as the result of simply applying an electrostatic potential to a trivial bulk band insulator. As shown in (e), they are localized edge states that are distinct from the dispersionless flat band states [top panel in (e)] and from the typical quantum well bound states [bottom panel in (e)] as well. These results agree with the prediction from the general continuum model.

IV. CONCLUSION AND DISCUSSION
To summarize, we have predicted a class of in-gap edge excitations with spontaneous domain-wall spin textures in insulating Dirac-type systems of massive spin-1 particles with only a locally applied electrostatic potential. Despite the absence of magnetism and any a priori topological origin, these states are extremely robust against boundary deformation and disorders. The remarkable property of these states is the self-induced emergence of domain-wall spin ordering that renders distinct spin-angular momentum locking on different sides of the domain interface. Consequently, the states are stable against impurities and/or geometric deformation. The in-gap modes are formally three-component evanescent wave solutions, bearing certain resemblance with the Jackiw-Rebbi type of bound states. The modes belong to a distinct class due to the following physical reasons: three-component spinor wave function, unusual boundary conditions, and a shifted flat band induced by the external scalar potential. Our findings provide a fully electrostatic based route to generating protected, robust spin ordering edge states without requiring any sort of magnetism, extrinsic or intrinsic. The states can be exploited for spintronics and quantum information processing applications, e.g., realization of a gatecontrolled spin-1 Dirac electron transistor or quantum switch. With rapid advances in generalized Dirac materials, especially those hosting the spin-1 generalization of ordinary Dirac/Weyl fermions, and with the state-ofthe-art measurement technologies, experimental confirmation of the states discovered here is possible.
We note a distinct feature of the system studied: the inherent mid-gap flat band hosting macroscopically degenerate states. Without the applied electrostatic potential (V g = 0), we obtain the flat band states, i.e., E(p) = 0, given by (nonnormalized) with the wavevector k = (k x , k y ) ≡ p/ making an angle ζ = arctan(k y /k x ) with the x axis. The states result in a vanishing current and a trivial spin distribution over the space as well as a vanishing Chern number [37,66,67]. Our finding is that a locally applied potential shifts the flat band relative to the surrounding and surprisingly leads to a class of exotic edge excitations that inherit the (quasi)flat dispersionlessness but attain a nontrivial feature associated with the emerging domain-wall like spin ordering. Due to the vanishing Chern number of the flat band, in the configuration in Fig. 1B, the regions with different applied potential V g possess the same Chern number. This indicates that the uncovered in-gap states do not have a topological origin. It has been known that flat bands can lead to exotic physical phenomena such as zero-refractive index, unconventional Anderson localization [68,69], itinerant ferromagnetism [70], and unconventional superconductivity [1,71,72]. Moreover, the finite gap opening makes it possible to categorize the unperturbed bulk system into the phase of class-D with a particle-hole symmetry and a broken time-reversal symmetry, which also arises in p + ip superconductors [37]. In this regard, the two-dimensional gapped pseudospin-1 system represents a paradigm to investigate high-spin topological phases with exotic edge excitations and flatband physics. With enriched pseudospin degrees of freedom, graphene-based heterostructures, such as graphene-In 2 Te 2 bilayer [65] and twisted bilayer graphene superlattice [73], can also be exploited for possible experimental realization of the topological edge states uncovered in this paper.
Taken together, the main contributions of this paper are: (1) in-gap edge modes can arise in a topologically trivial spin-1 Dirac insulators with local electrical gating or nonmagnetic doping, (2) the in-gap edge modes possess pseudospin polarized textures akin to localized domain walls of either the hedgehog or the vortex type without requiring any external pseudospin resolved field, (3) the edge modes are robust against boundary deformations and disordered scalar impurities, (4) the edge modes are nearly dispersionless in energy and intrinsically possess the capability of strong charge and spin confinement/localization, and (5) all these features of the ingap edge modes can be electrically controlled within the same material setting. We note that, the existing mechanisms for in-gap bound modes or excitations can be either topological or nontopological. Examples are the extensively studied topological in-gap edge modes [15,16], the nontopological Yu-Shiba-Rusinov bound states associated with magnetic impurities in superconductors [74][75][76], vacancy defects or particular lattice terminations induced bound states in crystalline lattice systems [77,78], and modes induced by nonmagnetic impurities in topologically nontrivial band insulators [79]. There was also a recent work [80] hinting that the multicomponent character of the Dirac-Bloch wavefunction and the associated boundary conditions would enable nontopological Dirac materials, through proper engineering of the graphene lattice boundaries, to potentially host robust surface states. The system of pseudospin-1 Dirac insulators that we have studied does not require any special lattice engineering, does not involve any magnetic-type of perturbations or defects either, nor does it have a nontrivial band topology. Yet, robust in-gap edge modes can arise.
Our system thus does not fall into any known category of systems in which in-gap bound modes can arise, and the edge modes uncovered belong to a distinct class due to the three-component spinor wave function and the unusual boundary conditions as well as an electrically induced shift of the flat band. In the position representation r = (x, y), the Hamiltonian for a massive spin-1 generalization of Dirac/Weyl fermion readsĤ where v F is the Fermi velocity,p is the momentum operator,Ŝ = (S x , S y ) andŜ z are spin-1 matrices, ∆ denotes a Dirac-type mass, and U (r) is a scalar type of perturbation (e.g., an electrostatic potential). The energy eigenstates Ψ(r) = [ψ 1 (r), ψ 2 (r), ψ 3 (r)] T can be determined by the generalized Dirac-Weyl equation For a spatially homogeneous/constant potential, e.g., U (r) = V 0 , the eigenenergies are E = V 0 and V 0 + s ∆ 2 + v F |k| 2 with s = ± being the dispersion band index. The corresponding plane wave solutions can be written as where the wavevector k = (k x , k y ) has length k = √ 2 − δ 2 with = (E − V 0 )/ v F , δ = ∆/ v F , which makes an angle ζ = arctan(k y /k x ) with the x axis. Other factors are α = k/( − δ) and β = k/( + δ). The current operator is defined based on Eq. (A1) aŝ The local current associated with state Ψ(r) = [ψ 1 , ψ 2 , ψ 3 ] T can be calculated from the local expecta-tion value ofû as By definition, the local current is the local probability density of spin vector (S x , S y ). Using the plane wave solution (A3), we obtain The effects of the applied scalar potential are to shift the Dirac point (k = 0) in the energy domain, to tune the kinetic energy = (E − V 0 )/ v F , and to alter the particle attributes from hole-to electron-type, and vice versa.
The time-reversal symmetry operator is where K is the operator for complex conjugation. Due to the Dirac-like mass term, the time-reversal symmetry is broken.

Appendix B: Eigensolutions of type-II quantum dots of massive spin-1 particles
We obtain the eigensolutions of the spin-1 massive Dirac system where an electrostatic potential is applied to a circular domain: U (r) = V 0 Θ(r − R). This is effectively a type-II quantum (anti-)dot configuration for Dirac-like massive spin-1 particles. Because of the rotational symmetry, it is convenient to use polar coordinates r = (r, θ), where the eigenequation iŝ Because the total angular momentum operatorĴ z = −i∂ θ +Ŝ z commutes with the HamiltonianĤ, the common set of eigenstates has the general form with l ∈ Z. For the dispersive bands, we have where the index µ = I, O labels the inner and outer regions of the circular domain boundary, m (x) are the Bessel and the Hankel functions of the first kind, respectively. Matching the spinor wavefunctions Ψ I l and Ψ O l at the domain boundary (interface) r = R yields the following transcendental equation which can be calculated numerically to yield the eigenenergies and eigenstates with high accuracy. Figure 7 shows some representative results. For reference, we have also included the corresponding results for the standard massive, spin-1/2 Dirac fermion system. We see that, for the massive spin-1 system, apart from the conventional quantum dot bound states, an additional group of modes emerge in the gap. While edge states can arise in the band gap as in conventional topological insulators, some kind of magnetic perturbations are required [8][9][10][11][12][13][14][15][16]. As there is no magnetic perturbation of any sort in our quantum dot system for massive spin-1 Dirac particles, the emergence of the states in the band gap is quite counterintuitive and striking.
We show analytically that the modes in the band gap possess a unique spectral peculiarity and are in fact edge states with domain-wall like, topologically nontrivial spin textures. In particular, in the gap |E µ | < |∆|, the radial wavenumbers are purely imaginary, which can be redefined as With the substitutions we rewrite the eigenvalue equation Eq. (B4) as with the associated eigenstates given by where ρ = r/R, I l (x) and K l (x) are modified Bessel functions. Making use of asymptotic expansions of high order Bessel functions [81], i.e., l 1: Using the identity lim n→∞ (1 + 1/n) n = e, we arrive at an equation that can be solved to yield the asymptotic eigenenergies: The eigenenergies are independent of the angular momentum and are thus in-gap (energy) dispersionless excitations. The associated eigenstates are approximately given by So, inside the domain ρ < 1, we have Outside of the domain ρ > 1, we have We thus have that the in-gap excitations are localized edge modes and exhibit domain-wall like spin textures for high angular momentum values. Note that, for a given value of l, in the semiclassical limit p, q 1, we have, approximately, From the eigenvalue equation, we have which leads to the same in-gap spectral properties as those from the large l regime. The associated semiclassical eigenstates are where κ = p ≈ q ∼ δ 2 − v 2 0 /4 1. We obtain the resulting spin textures as electrostatic potential. Semiclassically, the in-gap states are thus exponentially localized edge modes with spontaneously topological spin textures, which are reminiscent of the interfacial Jackiw-Rebbi modes but here the modes have a distinct spectral features and an unconventional physical origin.
Appendix C: Effects of smoothly varying electrostatic potential profiles and impurities on in-gap modes in massive spin-1 Dirac systems Realistically, the applied electrostatic potential will not be infinitely sharp at the domain boundary but, rather, the potential file varies smoothly across the boundary. From an experimental standpoint, it is necessary to investigate if the in-gap states can persist when the domain boundary is "smeared." The test would provide further support for the robustness and topological origin of those states. To be concrete, we use the following smoothly varying potential profile: where d (1/d) characterizes the boundary smoothness (sharpness) with d = 0 corresponding to the ideal case of an infinitely sharp boundary. Generally, for a finite value of d, it is not feasible to write down explicit solutions of the spin-1 Dirac equation. We thus exploit the finite difference method (FDM) recently developed for massless spin-1/2 Dirac fermions [82][83][84][85] and generalize it to massive spin-1 particles. In particular, taking advantage of the rotational symmetry of U (r) and using the polar decomposition ansatz we obtain the corresponding radial eigenvalue equation of the three-component spinor wherê When discretizing this equation on a finite lattice/grid, we need to judiciously specify the difference scheme and the boundary conditions at the ends of the lattice so as to preserve the Hermiticity of the Hamiltonian. A feasible procedure is to use the backward-forward-backward (BFB) difference scheme to approximate the derivatives of the three components in Eq. (C3): where h = L/(N + 1) is the discretization step size for the system in the range 0 < r < L with N + 2 lattice points. The boundary conditions can be deduced from the Hermitian constraint ofĤ r : which can be explicitly written as (C5) The specific boundary conditions on R(0) and R(L) then become R 1 (0)+R 3 (0) = 0 and R 2 (L) = 0. Implementing this procedure results in an eigenvalue problem for a 3N × 3N Hermitian matrix H 3N ×3N = [H µν ] with entries given by for n = 1, · · · , N . For n < N , the matrix elements are We use the typical experimental values of the local density of states (DOS) [82,84,85] to measure the spectral features and study the effects of the smooth potential profile and impurity on the in-gap states, where the DOS is defined as with ν labeling the obtained radial eigenstates for fixed l, and represents a spatial average of the wave function centered at r = r 0 with a Gaussian weight λ. We approximate the delta function by a Lorentzian with the broadening parameter Γ. In our simulations, we use a system of size L/R = 10 and discretize it with a uniform lattice of N = 600 sites. Other parameters are chosen as Γ/E * = 0.2 and λ = 0.01R. Representative results are shown in Figs. 8, 9 and 10, which provide strong support for the persistence of the in-gap modes in massive spin-1 Dirac systems in realistic systems with a smooth potential profile and impurities.
Appendix D: Multiple multipoles method: calculation of eigenenergies and eigenstates of massive spin-1 Dirac particle in arbitrary domains To test the robustness and to establish the topological origin of the in-gap states for massive spin-1 Dirac particles analytically predicted from the setting of a circular potential domain, we seek to search for such states in systems with a deformed domain. A difficulty that must be overcome is to calculate the eigenenergies and eigenstates of massive spin-1 Dirac particle in deformed domains of an arbitrarily geometric shape. We have succeeded in generalizing the multiple multipole expansion method originally developed in optics [86][87][88][89][90] to massive spin-1 Dirac particles. The end result of this nontrivial generalization is a systematic, reliable, accurate, and efficient computational paradigm incorporating the evanescent waves to detect and ascertain the existence of in-gap excitations/modes for arbitrarily shaped electrostatic potential domains.

FIG. 11.
Schematic illustration of the setting of multiple multipole expansion method. The domain in which an electrostatic potential is applied has boundary Γ separating regions I and II. The basis functions originated at rm I (blue circular dots) are used to determine the wavefunction in region II, while those at rm II (red circles) determine the wavefunction in region I. The boundary conditions for the massive spin-1 Dirac wavefunctions are imposed at the collocation points rj ∈ Γ.

Method implementation
A concrete setting of a single potential domain of arbitrary shape is illustrated in Fig. 11, where the exact shape of the geometric boundary is specified according to the superformula in botany [55], a simple but powerful prescription that can generate a vast variety of complex geometric shapes. In polar coordinates, the superformula is where the parameters (m 1 , m 2 , n 1 , n 2 , n 3 ; a, b) control the shape. The boundary defines two sub-regions, one exterior another interior, denoted by I and II, respectively, as shown in Fig. 11. The three-component spinor wave equation for a massive spin-1 Dirac particle in each subregion τ ∈ {I, II} reads where δ = ∆/ v F and τ = (E − V τ )/ v F . In polar coordinates r = (r, θ), the spinor cylindrical wave basis of the solutions with angular momentum l is where being the Hankel function of the first kind), we have that the Dirac-type expansion basis wavefunctions originated at r mτ for the specific region τ are given by where τ denotes the complement of τ , d mτ ≡ |d mτ | = |r − r mτ | and θ mτ = Angle(r − r mτ ) with r ∈ τ . Carrying out the expansion in region II, we obtain the wavefunction as The wavefunction in region I has the form where Ψ in (r) = 1 2 denotes the input source triggered by an applied external excitation outside of the domain [c.f., top panel of Fig. 11]. Imposing the relevant boundary conditions parameterized by the angle φ between the outward normal at any boundary point r j and the x-axis: j χ in = 1 2 α I e iφ + β I e −iφ e ik I (|rj | cos θj −x0) .
In principle, the set consists of an infinite number of equations with an infinite number of undetermined expansion coefficients C m II l and C m I l . To solve the system numerically, a finite truncation is necessary, which turns out to be feasible in practice by discretizing the boundary to a finite number of points J and setting the number of basis functions M τ in the specific region τ and l ∈ [−L, L] for all the functions. Carrying out the discretization procedure, we arrive at the following finite dimensional matrix equation where N = (2L + 1) × (M I + M II ) = N I + N II and the compact substitutions are As the expansions are generally nonorthogonal, more equations are required than the number of unknowns to enable the deduction of an over-determined matrix system with 2J N , which can be solved by the standard pseudo-inverse algorithm: C = −pinv(M) * Y . In particular, we use the residual error evaluated at the boundary Error = ||M * C + Y || ||Y || as the criterion to test convergence. We adjust the number, the order and/or positions of the multipoles to ensure Error < tolerance. After the unknown coefficients C have been obtained, the associated wavefunctions and hence the local density of states in the specific region can be calculated accordingly.

Method validation
To validate the method, we exploit the analytically solvable case of circular geometry. Figure 12 shows a comparison of the eigenenergy spectra obtained analytically and calculated from the multiple multipole method. The agreement is remarkable.