Dynamics of Dirac solitons in networks

We study dynamics of Dirac solitons in prototypical networks modeling them by the nonlinear Dirac equation on metric graphs. Soliton solutions of the nonlinear Dirac equation on simple metric graphs are obtained. It is shown that these solutions provide reflectionless vertex transmission of the Dirac solitons under suitable conditions. The constraints for bond nonlinearity coefficients, allowing reflectionless transmission over a Y-junction are derived. The analytical results are confirmed by direct numerical simulations.


Introduction
Nonlinear evolution equations on networks have attracted much attention recently [1,2,3,4,5,6,7,8,9,10]. Such interest is caused by a broad variety of potential applications of the nonlinear wave dynamics and soliton transport in networks, such as Bose-Einstein condensates (BECs) in branched traps, Josephson junction networks, the DNA double helix, polymer chains, etc.
Despite the rapidly growing interest in wave dynamics on networks, most of the studies are mainly focused on nonrelativistic wave equations such as the nonlinear Schrödinger (NLS) equation [1,2,3,4,5,6,7,10]. Nevertheless, there is a number of works on the sine-Gordon (sG) equation in branched systems [8,9,11]. However, relativistic wave equations such as the nonlinear Klein-Gordon and Dirac equations are important in field theory and condensed matter physics and hence exploring them on metric graphs is of interest in its own right. These graphs consist of a system of bonds which are connected at one or more vertices (branching points). The connection rule is associated with the topology of a graph. When the bonds can be assigned a length, the graph is called a metric graph.
In this paper we address the problem of the nonlinear Dirac equation on simple metric graphs by focusing on conservation laws and soliton transmission at the graph vertices. Our prototypical example will be the Y-junction. Early studies of the nonlinear Dirac equation date back to Thirring [12] and Gross-Neveu [13] models of field theory.
Integrability, the nonrelativistic limit and exact solutions of the relevant models have been considered in, e.g., [14,15,16,17,18] among others. A potential application of the nonlinear Dirac equation to neutrino oscillations was discussed in [19]. Recently, the possibility of experimental realization of Dirac solitons in Bose-Einstein condensates in honeycomb optical lattices was discussed in [21] where the nonlinear Dirac equation was derived. Further analysis of this setting in Refs. [22,23,24,26] has excited a rapidly growing interest in the nonlinear Dirac equation and its soliton solutions (see, e.g. [27]- [49]). In [27,28,31,32,33,34], a detailed study of soliton solutions for different types of nonlinearity, their stability and discussions of conservation laws were presented. These developments have, in turn, had an impact also on the mathematical literature where the stability of solutions for different forms of the nonlinear Dirac equation was explored in [29,30,36]. In [37]- [41], the stability of Gross-Neveu solitons in both 1d and 3d and for other nonlinear Dirac models an important energy-based stability criterion was developed. The numerical corroboration of the stability for solitary waves and vortices in 2D nonlinear Dirac models of the Gross-Neveu type was considered quite recently in [35].
Dirac solitons in branched systems can, in principle, be experimentally realized in different systems of optics, but also envisioned elsewhere (e.g. in atomic physics etc.).
A relevant such possibility consists of the discrete waveguide arrays of [44,49] that can be formulated as a branched system to be described (in the appropriate long wavelength limit [44]) by a nonlinear Dirac equation on metric graphs. Such networks in optical systems have been proposed, e.g., earlier in [50]. More exotic possibilities can be envisioned in branched networks of honeycomb (i.e., armchair nanoribbon) optical lattices for atomic BECs [23], although we will not focus on the latter here.

Conservation laws and vertex boundary conditions
The nonlinear Dirac system that we are going to explore, i.e., specifically the Gross-Neveu model, follows from the Lagrangian (in the unitsh = c = 1) [27,28,31] and describes interacting Dirac fields with 0 and 1 corresponding to time and coordinate variables, respectively. Here g is the nonlinearity coefficient characterizing the strength of the nonlinear interaction. The field equations for this Lagrangian lead to the nonlinear Dirac model of the form where and For a metric star graph consisting of three semi-infinite bonds (see, e.g., Fig. 1), the Lagrangian density of the nonlinear Dirac field for each bond can be written as L j where j parametrizes the bond and L for each is of the form of Eq. (1).
The field equation following from this Lagrangian density can be written as Eq. (2) before, where the spatial coordinates are defined as x 1 ∈ (−∞, 0] and x 2,3 ∈ [0, ∞), while 0 coincides with the graph vertex.
The formulation of an evolution set of equations on metric graphs requires imposing vertex boundary conditions which provide "gluing" of the graph bonds at the graph vertices. For the linear Dirac equation on a metric graph studied earlier by Bolte and Harrison in [51], the general vertex boundary conditions have been derived from the self-adjointness of the Dirac operator on a graph. In that case such conditions led to Kirchhoff rules and continuity of the wave function at the graph vertex [51]. However, for the nonlinear problem in order to derive vertex BCs, it is arguably more natural to use suitable conservation laws of the nonlinear flow, which give rise to appropriate BCs. For the nonlinear Dirac system fundamental conservation laws can be presented in terms of the momentum-energy tensor given by [27,28] The energy on each bond of a star graph can be written as where integration is performed along the bond b j . Here, to derive the vertex boundary conditions (VBC) we use conservations of charge and energy. The total charge for the star graph is defined as where the charge for each bond is given by The conservation of charge (Kirchhoff rule for charge current),Q = 0 yields the following vertex boundary condition: Here we used the asymptotic conditions We note that the boundary conditions (6) can be derived also from the conservation of the current density given by j = ϕχ * + χϕ * .
The energy conservation,Ė = 0 leads to As per the above analysis, the VBC for NLDE on metric graphs have been obtained from the energy and charge conservation. However, the VBC given by Eqs. (6) can be fulfilled, if the following linear relations at the vertices are imposed (see Appendix for details): where α 1 , α 2 , α 3 are the real constants which will be determined below. In the following we will use Eqs. (9) as the vertex boundary conditions for Eq. (2) on a metric star graph. When applied to the exact soliton solutions of the NLDE, the VBC will lead to algebraic conditions connecting α j and g j , as will be derived in the next section. We note that in the linear limit (g j → 0) the vertex boundary conditions given by Eq. (9) preserve the self-adjointness of the (linear) Dirac equation on metric graph, i.e. belong to the class of general boundary conditions derived by Bolte and Harrison in [51].

Soliton dynamics and vertex transmission
In [27,28,31]  We look for the soliton solution of NLDE on a metric star graph in the form Then, from Eqs. (2) and (10) we have The vertex boundary conditions for the functions A j and B j can be written as using the standing wave soliton solutions given by Eqs. (14) and (15)(x 0 = 0) for A prototypical static solution of the system (11), (12) vanishing at x → ±∞ can be written as [28] A where x 0 is the position of the soliton's center of mass and β = √ m 2 − ω 2 . In order for these solutions of Eqs. (11) and (12) to solve the problem on the metric graph, they need to also satisfy the vertex boundary conditions (13). This can be achieved if the constants α j and coupling parameters g j fulfill the following conditions: which stems from the first of Eqs. (9) and in accordance with the second one of conditions (9). The combination of the two leads to the "sum rule": It is important to note that this sum rule is derived by assuming that incoming wave comes from the first bond, b 1 , while if it comes from the bond b 2 (or b 3 ) one should replace g 1 in Eq. (18) with g 2 (or g 3 ). In Fig. 2 plots of the soliton solutions of NLDE on metric star graph corresponding to Eqs. (14) and (15), satisfying the vertex conditions are given.
We note that soliton solutions given by Eqs. (14) and (15) describe the standing waves in metric graphs. The traveling wave (soliton) solutions of NLDE can be obtained by considering the case of the moving frame and taking into account the Lorentz invariance [28], in the case of the homogeneous domain (i.e., without a Y-junction).
Assuming implicitly that the solitary wave is centered around a position well to the left of the origin, the Lorentz transformations between the frames moving with relative velocity, v can be written as [28] x where Using these transformations, the traveling wave (soliton) solutions of the nonlinear Dirac equation in the moving frame determined by the constraints (16) can be written where x ′ , t ′ and η are given by Eqs. (19) and (20).
It is important to note that this transformation does not affect the scaling dependence on the coupling parameters g j . For that reason, we expect the form of Eqs. (16) to still be valid in the case of traveling/moving wave solutions.
To analyze the dynamics of the Dirac solitons on a metric graph, we solve numerically the time-dependent NLDE given by Eq. (2) for the vertex boundary conditions given by Eq. (9) by considering the an initially traveling waveform localized in bond 1, far from the vertex. We do this for the case when the constraint given by Eq. (18) is fulfilled, as well as for the case of arbitrary g j which do not respect the relevant constraints. Fig. 3(a) presents the plots of the solution of Eqs. (2) and (9) obtained numerically for the values of g j obeying the constraint (18). For vertex boundary conditions given by Eq.(13) parameters α j are chosen as α 2 = 1/|g 2 |, α 3 = 1/|g 3 |, One can observe the absence of the vertex reflection in this case. Fig. 3 conditions (9) obtained numerically for those values of g j which do not fulfill Eq. (18).
Here, reflection can be clearly discerned. Fig. 4(a) shows the conservation of the energy and reflectionless transmission through the Y-junction during the propagation of Dirac soliton in graph. The emergence of the vertex reflection when the constraints (18) are not fulfilled, can also be seen from the Fig. 4(b). Here the vertex reflection coefficient, which is determined according to the definition R = E 1 (t = 7)/E(t = 7) (with E(t = 7) = E 1 (t = 7) + E 2 (t = 7) + E 3 (t = 7)) is plotted as a function of g 3 for fixed values of g 1 and g 2 . This systematic analysis clearly illustrates the necessity of the symmetric scenario of g 3 = √ 2, consonant with our vertex sum rule, for reflection to be absent. A similar phenomenon was observed in the case of other nonlinear PDEs on metric graphs, such as the NLS and sG equations, considered earlier in the Refs. [1] and [11], respectively. For the initial condition corresponding to a traveling wave of

Extension for other graphs
The above treatment of the NLDE on the metric star graph can be extended to the tree graph presented in Fig. 5.  Similarly to that for star graph, one can obtain soliton solutions of the nonlinear Dirac equation on the tree graph for the moving frame.
Another graph for which a soliton solution of the nonlinear Dirac equation can be obtained, is presented in Fig. 6. It has the form of a triangle whose vertices are connected to outgoing semi-infinite leads. The vertex boundary conditions for Eqs. (11) and (12) which follow from the conservation of current and energy can be written as In short, the conditions matching the of charge and energy conservations must be applied to each node connecting the multiple vertices of the graph. The soliton solution obeying these boundary conditions can be written as where b is considered to be an index running from 1 to 6, s 1 = s 2 = s 3 = l, s 4 = l + L 2 , s 5 = l + L 2 , s 6 = l + L 3 , L 3 = L 2 + L 4 (this condition is necessary for the waves traveling along b 3 and those split along b 2 and b 4 to arrive at the vertex between b 3 and b 4 concurrently) and coefficients g j fulfill the following constraints: Again, the traveling wave (soliton) solutions can be obtained in case of the moving frame using the Lorentz transformations. This approach for designing soliton solutions of the nonlinear Dirac equation can be applied to other graph topologies with multiple junctions, provided a graph consists of finite parts connected to outgoing semi-infinite bonds.
In Fig. 7 the solution of NLDE on each bond of this triangle graph is plotted for the cases when the sum rule given by Eqs. (24)-(26) is fulfilled (a) and is broken (b). Absence of the vertex reflection in Fig. 7(a) is clearly seen, while in Fig. 7(b) the transmission of Dirac solitons through the graph vertices is clearly impeded by the mismatch at the VBC. Thus, reflection events are clearly discernible in bonds such as 1, 3 and 4.

Conclusions & Future Challenges
In this paper we studied dynamics of Dirac solitons in networks by considering the case of metric graphs. We obtained soliton solutions of the nonlinear Dirac equation on some of the simplest metric graphs such as the star, tree and triangle graphs. Constraints Our computations clearly suggest that the considered vertex boundary conditions are necessary (yet not necessarily sufficient) conditions for reflectionless transmission through the junctions of interest. This is also in line with the recent homogenization calculation of [54]. It is thus natural to formulate the following conjecture. Assume that the initial data correspond to a traveling wave of the NLDE in the form of Eq. (21) centered at x 0 , with v > 0. Then for t → +∞, the solution can consist asymptotically of a superposition of solitary waves in branches 2 and 3, only if the condition of Eq. (18) holds. This conjecture poses a challenging mathematical question for further rigorous study in the problem at hand, both at the level of the NLDE, but also by direct analogy in the context of the NLS and other similar models.
It should be noted that an important, additional set of challenges emerges as regards the study of stationary solutions involving the vertex. In the language of the very recent work of [54], the states considered herein are the so-called half-soliton states. Yet in the NLS framework of the latter work, additional so-called shifted (non-monotone in the different branches) states also exist, yet they may be spectrally unstable, depending on whether they are monotonic or non-monotonic in the outgoing edges of the metric graph. The half-soliton considered here is especially interesting as at the NLS level it turns out to be spectrally stable, yet nonlinearly unstable. This issue of the (more involved at the NLDE level [55]) spectral stability of these stationary states is also an especially interesting direction for future study.
Naturally, the above study can be extended for other simple topologies, such as a graph with one or multiple loops (e.g. the dumbbell graph), or other combinations of star, tree and loop graphs.
It will be extremely interesting if these ideas induce experimental interest in systems such as suitably tailored discrete optical waveguides [52,53] or in atomic settings, in the same way as they have for instance for the propagation of traveling waves through networks of granular crystals [56].

Acknowledgement
This work is partially supported by the grant of the Ministry for Innovation Developments of Uzbekistan (Ref. No.BF-2-022).

Appendix A
Here we will show that the linear vertex boundary conditions given by Eqs. (9) lead to the ones given by (6)- (8). Consider the following (linear) relations given at the vertex