Hinge solitons in three-dimensional second-order topological insulators

A second-order topological insulator in three dimensions refers to a topological insulator with gapless states localized on the hinges, which is a generalization of a traditional topological insulator with gapless states localized on the surfaces. Here we theoretically demonstrate the existence of stable solitons localized on the hinges of a second-order topological insulator in three dimensions when nonlinearity is involved. By means of systematic numerical study, we find that the soliton has strong localization in real space and propagates along the hinge unidirectionally without changing its shape. We further construct an electric network to simulate the second-order topological insulator. When a nonlinear inductor is appropriately involved, we find that the system can support a bright soliton for the voltage distribution demonstrated by stable time evolution of a voltage pulse.

Traditional topological insulators in n dimensions support edge states localized at (n − 1)-dimensional boundaries, such as the Chern insulator. Recently, topological phases have been generalized to the higher-order case where gapless edge states are localized at (n − m)dimensional (with m > 1) boundaries . Such insulators are coined higher-order topological insulators, * yongxuphy@tsinghua.edu.cn or specifically, mth-order topological insulators. For instance, in two dimensions (2D), there exists type-I [40] and type-II [62] quadrupole topological insulators that support zero-energy modes localized at the corners. In three dimensions (3D), there exists a second-order topological insulator hosting chiral modes localized on the hinges, when an appropriate term breaking the timereversal symmetry is involved in a Z 2 topological insulator [50]. In the context of nonlinearity, it is natural to ask whether the second-order topological system can support solitons localized on the hinges when nonlinearity is involved.
In this paper, we theoretically predict the existence of solitons in a 3D second-order topological insulator when nonlinearity is involved. Such solitons result from the balance between nonlinearity and dispersion of the hinge modes. For the higher-order topological insulator, there are two regimes: strong regime with chiral hinge modes and weak regime with hinge modes that are not chiral. We find that the solitons can exist in both of these two regimes. Yet, we show that the soliton in the strong regime is more stable than the one in the weak regime, which gradually slows down. Furthermore, we propose a scheme to realize the second-order topological insulator in electric circuits with nonlinearity that can be realized by a voltage-controlled variable inductor. By simulating the dynamics of the circuit, we finally show that a bright soliton can stably exist in the system.

II. LINEAR MODEL HAMILTONIAN
We start by considering the model described by the following Hamiltonian in momentum space [50] where σ ν and τ ν represent Pauli matrices acting on different degrees of freedom, σ 0 and τ 0 refer to 2×2 identity matrices. J, M , ∆ 1 and ∆ 2 are real parameters. When ∆ 2 = 0, the Hamiltonian describes the 3D Z 2 topological insulator [87]. It respects both time-reversal symmetry, i.e., T H(k)T −1 = H(−k) with T = iσ y κ (κ denotes the complex conjugation operator), and the rotational symmetry, i.e., When 1 < |M/J| < 3, the Hamiltonian represents a strong topological insulator with odd number of Dirac points in the energy spectra of surface states; when |M/J| < 1, it describes a weak topological insulator with even number of Dirac points in the surface energy spectra [88,89].
To generate the higher-order topological insulating phase, Ref. [50] introduces the term involving ∆ 2 that breaks both the time-reversal symmetry and the rotational symmetry but preserves the symmetry of their product, i.e., C z 4 T . This term provides an effective mass, opening the gap of the Dirac points on the surface vertical to either x or y axis. Since the signs of the Dirac masses are opposite for neighboring surfaces, the hinges of their intersections exhibit chiral modes, similar to the chiral modes of a Chern insulator. It has been found that four chiral modes exist on the four hinges between the surfaces normal to x and y axes. The topological property of this higher-order phase can be characterized by the Chern-Simons term where A a;n,n = i ϕ n (k)| ∂ ka |ϕ n (k) with a, b, c ∈ {x, y, z} is the k a -component of Berry connection, |ϕ n (k) is the Bloch eigenstate of H L (k) and n, n are running over the occupied bands of the system. The topological invariant is quantized to 0 or π by C z 4 T symmetry [50].
For simplicity without loss of generality, we set ∆ 1 = ∆ 2 = J. In Fig. 1, we plot the energy spectra in a geometry with open boundaries along x and y and periodic boundaries along z. We see that there exist chiral modes localized on the hinges when M = 2J in the strong topological regime. When M = 0.8J in the weak topological regime, we can still observe the hinge-localized modes, but they are not chiral. We find solitons in both of these two cases, but their stabilities manifested by time evolution are different, which will be discussed in detail in the following.

III. HINGE SOLITONS
To create the hinge solitons, we add nonlinear terms into the higher-order topological system so that the dynamics is governed by the following nonlinear Schrödinger equation, where H N (ψ) (r,λ),(r ,λ ) = δ rr δ λλ g λ |ψ λ (r)| 2 represents the short-range nonlinear interactions with g λ denoting their strength for each component, H L is the real space counterpart of the Hamiltonian (1) and ψ(r, t) is the wave function. Here we consider a cubic local nonlinear interaction, which widely exists in ultracold atoms and nonlinear optics [2,90]. For simplicity, we set g λ = g > 0.  Following Ref. [91], we construct a wave packet using the hinge modes of H L to determine a soliton solution, where λ = 1, 2, 3, 4 is the component index, ϕ kz denotes the hinge mode of H L (k z ) corresponding to the eigenenergy ε kz and a(κ, t) denotes the weights for the corresponding hinge mode over time.
We now apply the Taylor-series expansion in κ for ϕ kz+κ,λ (x, y) in the above integral to simplify Eq. (4) to where a(z, t) = π −π a(κ, t)e iκz dκ is the envelope function of the corresponding soliton. Applying H L to the wave function ψ in Eq. (5) and then performing the Taylorseries expansion for both ε kz+κ and ϕ kz+κ,λ (x, y) yields Suppose that ϕ kz varies much more slowly with respect to k z than ε kz , we can make the approximation that ∂ j (εϕ)/∂k j z ≈ ϕ∂ j ε/∂k j z . By further assuming that a(z, t) varies slowly in time and space, we obtain an approximate nonlinear Schrödinger equation for the envelop where ε = ∂ε kz /∂k z , ε = ∂ 2 ε kz /∂k 2 z and g eff = g x,y,λ |ϕ kz,λ (x, y)| 4 . Evidently, this equation is exactly the same as the Gross-Pitaevskii equation for a BEC in a free space in a reference moving with the velocity of −ε with 1/ε playing the role of an effective mass. If g eff = 0, this equation does not support a localized solution. Because of the chiral property of the hinge modes in the strong regime, the solitons can only move along one direction on a hinge, which is fundamentally different from solitons in a one-dimensional system whose velocities can be either positive or negative.
When ε < 0 for g eff > 0, this equation has a bright soliton solution with a density peak, where δ > 0 describes the energy difference between the nonlinearity-induced energy eigenvalue µ and the linear one ε kz at k z . While the Eq. (8) supports solutions for any positive value of δ, δ cannot be too large for a hinge soliton solution. On the one hand, the width of the soliton decreases as δ increases, meaning that if δ is very large, the width is very small so that we can not apply the approximation that the wave packet varies slowly in space. On the other hand, for large δ, µ = ε kz + δ can be moved into the bulk energy spectra, resulting in the scattering of the soliton into the bulk states. For a BEC, the total atom number N is given by λ |ψ λ (r)| 2 dr = N . For clarity, if we consider a zero-order approximation based on Eq. (5) so that ψ λ (r, t) ≈ e −iε kz t+ikzz ϕ kz,λ a B (z, t), the atom number constraint yields |a B (z, t)| 2 dz = N . We thus obtain , showing that the bright soliton becomes narrower with its peak value being increased as we increase g with N being fixed.
To show that a bright soliton localized on a hinge can exist in the second-order topological insulator, we now study the dynamics of the soliton by solving the nonlinear Schrödinger equation numerically. We first construct an initial wave packet using the hinge modes in the strong topological insulator based on Eqs. (5) and (8) and then observe the dynamics of the wave packet as time evolves. Fig. 2(a) illustrates that the wave packet moves unidirectionally along the z axis without changing its shape with the velocity determined by ε , implying the existence of a travelling soliton. In comparison, we also calculate the dynamics of the same initial wave packet in the linear case with g = 0 (without nonlinearities) and plot the results in Fig. 2(b). The figure shows that the wave packet diffuses significantly as time progresses. The diffusion of the wave packet happens in the linear scenario because the hinge spectrum is not perfectly linear. To show the difference quantitatively, we further compare the time evolution of the maximum value and the full width at half maximum of the wave packet for these two cases in Fig. 2(c). We can see clearly that in the nonlinear case, the peak value becomes stable after some slight oscillations, whereas in the linear case, it suffers a continuous decline. The propagation of the hinge soliton over time is also visually displayed in Fig. 2(d). The above analysis demonstrates the existence of a stable bright hinge soliton in the second-order topological insulator.
When ε > 0 for g eff > 0, Eq. (7) allows for the existence of the following dark soliton solution with a density dip, a D (z, t) = (δ/g eff ) 1/2 tanh[(δ/ε ) 1/2 (z − ε t)]e −iδt . (9) For a BEC with a fixed atom number, similar to the case of the bright soliton, if we consider atoms trapped in a cubic box with size L 1, then we obtain an expression of δ in terms of g ef f and N , i.e., δ = [ It indicates that the dark soliton becomes narrower when g is increased. To numerically demonstrate the dynamics of the dark soliton, we first construct an initial wave packet with a density dip based on Eqs. (5) and (9) and then compute the dynamics of the wave packet. Fig. 3 illustrates that the wave packet travels along z with its dip shape remaining almost unchanged during the time evolution. This demonstrates that a dark soliton can stably exist on the hinges in the second-order topological insulator.
In the strong topological regime when 1 < |M/J| < 3, there are chiral hinge modes in the energy gap. In contrast, in the weak topological regime when |M/J| < 1, instead of being chiral for the hinge modes, there are two states with opposite group velocities on each hinge for each energy in the energy gap [see Fig. 1(b)]. We now apply the same method to explore whether a stable bright soliton can exist on a hinge of a weak second-order topological insulator. Similar to the strong topological case, we find that the wave packet travels steadily along z without noticeable change of its shape, suggesting the existence of a stable bright soliton. But the soliton gradually slows down over time as shown in Fig. 4, possibly due to the scattering of the state to other modes with the same energy and opposite group velocity. In comparison, we also plot the velocity of the soliton with respect to time for the former case with M = 2J in Fig. 4, showing that the velocity remains stable. While the figure shows a very small decline of the velocity, it may be caused by the numerical error, as the decline becomes smaller when we use a smaller step size of time in the computation. Note that this numerical error in the latter case with M = 0.8J also exists but is not as noticeable as in the former one.

IV. HINGE SOLITONS IN ELECTRIC CIRCUITS
In this section, we will introduce a practical scheme to realize the second-order topological insulator in 3D and demonstrate that a bright soliton for the voltage distribution can exist in the presence of nonlinear inductors.
Let us now consider an electric network with a number of nodes labelled by m. Suppose that an electric current I m is externally flowing into node m and the current then flows into other nodes, which are connected with node m. According to Kirchhoff's law, we have where I mn denotes the current flowing from node m to node n, Y mn and Y m represent the admittance between node m and node n and the admittance between node m and the ground, respectively, and V m denotes the voltages measured at node m against ground. The relation above can be written in a matrix form as where L is the Laplacian with L mn = −Y mn + δ mn (Y m + n Y mn ) and I and V are the column vectors consisting of electric currents and potentials at each node. Consider the alternating current with the frequency of ω 0 (i.e., I = Ie iω0t and V = Ve iω0t ), we have To simulate the Hamiltonian in Eq.
(1), we need to engineer electric connections so that L = H L where H L is the real space Hamiltonian of H L (k) in Eq. (1). Since there are four degrees of freedom in each unit cell, we consider four nodes in each unit cell as shown in Fig. 5(a). Evidently, inductors and capacitors between two neighboring nodes contribute terms corresponding to the hopping with positive and negative amplitude, respectively. For the hopping with imaginary amplitudes, we utilize negative impedance converters with current inversion (INICs) [81], whose structure is shown in Fig. 5(a). Its node voltage equation is where I N = I 1 I 2 T , V N = V 1 V 2 T and υ = R b /R a (consider υ = 1 here) as shown in Fig. 5(a).
To simulate the dynamics of a soliton, we solve the following equation of motion of the circuit [81] d dt where C, Σ, and Π are the real-valued matrices characterizing the capacitance of capacitors, the conductance of the INICs, and the inductance of inductors, respectively.. These matrices in momentum space are given by where A C +M τ z and A L are contributed by the grounded capacitors and inductors, respectively. When J C = ∆ 1C and J L = ∆ 1L without M C and M L , we have A C = 8J C , M = 2J C , and A L = 8/∆ 1L using the grounding configurations listed in Table I. For alternating input currents with the frequency of ω [I(t) = Ie iωt and V (t) = Ve iωt ], this differential equation gives rise to the circuit Laplacian so that I = L(ω)V, which is the same as Eq. (12).
Let us now consider a specific case without any input currents, i.e., I(t) = 0. In this case, we reduce the Eq. (15) to the following first-order differential equation similar to the Schrödinger equation for the higher-order topological insulator circuit (HTIC), where Ψ = (V (t), V (t)) T , and This differential equation has the solution of Ψ(t) = e iωnt φ n , where φ n is an eigenvector of H HTIC corresponding to the eigenvalue ω n . Clearly, the eigenvalues of H HTIC appear in pairs as (ω n , −ω * n ). While H HTIC is non-Hermitian, its eigenvalues are real due to the positive semidefinite properties of C and Π [81], implying that the eigenvalues emerge in pairs as (ω n , −ω n ).
Since ω n is the alternating current frequency allowed by Eq. (20), we must have L(ω n )V = 0 given that the input currents are zero. To have nontrivial solutions, we require det(L(ω n )) = 0. In other words, there exist eigenvectors V n of L(ω n ) with zero eigenvalue, i.e., One can change the equation into the following form where It also tells us that if L(ω n ) has zero eigenvalues, then ω n is an eigenvalue of H HTIC . The intimate connection suggests that the band structures of H HTIC and H E (ω 0 ) = L(ω 0 ) share some common properties. For example, if H E (ω 0 ) is in a higher-order topological phase with chiral hinge modes, then there exist four zero-energy hinge modes at k z = π [see Fig. 1(a)]. Based on the discussion above, ω 0 should be an eigenvalue of H HTIC . In addition, according to Eq. (24), we have four modes of H HTIC with frequency of ω 0 localized on the hinges, suggesting that H HTIC has chiral hinge modes. Indeed, Fig. 6(b) illustrates that H HTIC exhibits the chiral hinge modes across the gap, which do not exist for periodic boundaries along x and y [see Fig. 6(a)]. The voltage distribution of these states also shows that they are localized on the hinges [see Fig. 6(c)]. We note that gapless nontrivial states in circuits have also been found in other 2D systems [67,81].
In stark contrast to a soliton for the atomic density in BECs, in our case, a soliton is for the voltage distribution in the electric circuit. To generate the soliton, we add nonlinear inductors in the circuit, which may be realized by a voltage-controlled variable inductor [92,93]. Each nonlinear inductor connects each node with the ground (see Fig. 5), contributing a nonlinear inductance described by Π N (r,λ),(r ,λ ) = g c δ rr δ λλ |V r,λ | 2 , where V r,λ denotes the voltage at the node labelled by r, λ.
To create a bright soliton, we initialize the state at t = 0 as Ψ(r, t = 0) = a 0 sech(λ 0 z)e ikzz φ kz (x, y), and then observe its dynamics. Here a 0 and λ 0 are real parameters and φ kz (x, y) is a hinge mode of H HTIC for a quasimomentum of k z in a geometry with open boundaries along x and y and periodic boundary along z. To simulate the dynamics, we impose open boundary conditions along all directions, in which case the parameters of the grounding elements at the corners and on the surface perpendicular to the z axis need to be adjusted to cancel the diagonal parts contributed by the electric elements connecting distinct nodes. For a specific example with g c = 10mH −1 V −2 , k z = 6π/7, ∆ 1C = ∆ 2C = J C = 50 µF, ∆ 1L = ∆ 2L = J L = 200 µH, ∆ 2R = 2 Ω, and M C = 100 µF, we take a 0 = 1 and λ 0 = 0.5. In Fig. 6(d) and (e), we plot the dynamics of the voltage peak with and without nonlinearity, respectively. It can be seen that the voltage distribution with nonlinearity is stable without conspicuous dispersion as time evolves. In stark contrast, in the case without nonlinearity, the wave packet for the voltage distribution becomes wider as it evolves over time. This shows the existence of a hinge bright soliton in the nonlinear circuit system.

V. CONCLUSION
In summary, we have theoretically studied solitons in a 3D second-order topological insulator when nonlinearity is involved. We find that solitons can exist in both the strong and weak regimes, but the soliton in the former regime is more stable than that in the latter regime given that the latter one gradually slows down as time progresses. In addition, the soliton on a hinge in the strong regime can only move along one direction due to the chiral property of hinge modes. To realize the soliton in a practical experimental system, we further introduce an electric network to simulate the 3D second-order topological insulator. By calculating the dynamics of a voltage pulse in the electric circuit, we show the existence of a soliton for the voltage distribution when nonlinear inductors are involved. Such solitons can also be experimentally observed in other nonlinear systems, such as ultracold atomic gases and nonlinear optics.

ACKNOWLEDGMENTS
This work is supported by the start-up fund from Tsinghua University, the National Thousand-Young-Talents Program and the National Natural Science Foundation of China (11974201).