Soliton transport in tubular networks: transmission at vertices in the shrinking limit
Abstract
Soliton transport in tubelike networks is studied by solving the nonlinear Schrödinger equation (NLSE) on finite thickness (”fat”) graphs. The dependence of the solution and of the reflection at vertices on the graph thickness and on the angle between its bonds is studied and related to a special case considered in our previous work, in the limit when the thickness of the graph goes to zero. It is found that both the wave function and reflection coefficient reproduce the regime of reflectionless vertex transmission studied in our previous work.
I Introduction
Particle and wave transport in branched structures is of importance for different topics of contemporary physics such as optics, cold atom physics, fluid dynamics and acoustics. For instance, such problems as light propagation in optical fiber networks, BEC in network type traps and acoustic waves in discrete structures deal with wave transport in branched systems. In most of the practically important cases such transport is described by linear and nonlinear Schrödinger equations (NLSE) on graphs. The latter has become the topic of extensive study during past few years Sobirov et al. (2010); Cascaval and Hunter (2010); Banica and Ignat (2011); Adami et al. (2011, 2012, 2012, 2013); Noja (2014); Gnutzmann et al. (2011); Sabirov et al. (2013) and is still rapidly progressing. Such interest in the NLSE on networks is mainly caused by possible topologydependent tuning of soliton transport in branched structures which is relevant to many technologically important problems such as BEC in network type traps Leboeuf and Pavloff (2001); Paul et al. (2007); de Oliveira et al. (2010), information and charge transport in DNA double helix Yomosa (1983); Yakushevich et al. (2002), light propagation in waveguide networks Burioni et al. (2006) etc.
Soliton solutions of the NLSE on simplest graphs and connection formulae are derived in Sobirov et al. (2010), showing that for certain relations between the nonlinearity coefficients of the bonds soliton transmission through the graph vertex can be reflectionless (ballistic). Dispersion relations for linear and nonlinear Schrödinger equations on networks are discussed in Banica and Ignat (2011). The problem of fast solitons on star graphs is treated in Adami et al. (2011) where estimates for the transmission and reflection coefficients are obtained in the limit of high velocities. The problem of soliton transmission and reflection is studied in Cascaval and Hunter (2010) by solving numerically the stationary NLSE on graphs. More recent progress in the study of the NLSE on graphs can be found in Adami et al. (2012, 2012, 2013); Noja (2014). Scattering solutions of the stationary NLSE on graphs are obtained in Gnutzmann et al. (2011), and analytical solutions of the stationary NLSE on simplest graphs are derived in Sabirov et al. (2013).
In metric graphs the bonds and vertices are one and zero dimensional, respectively. However, in realistic systems such as electromagnetic waveguides and tubelike optical fibers, the wave (particle) motion may occur along both longitudinal and transverse directions Markowsky and Schopohl (2014); Demidov et al. (2009); Biggs (2012). Therefore it is important to study below which (critical) thickness the transverse motions become negligible and the wave(particle) motion can be treated as onedimensional. In other words, studying the regime of motion when wave dynamics in such tubelike network can be considered the same as that in metric graph is of importance.
In this paper we study the NLSE on socalled fat graphs, i.e. on twodimensional networks having finite thickness. The geometry will be explained in more detail below, but see Fig.1 for a sketch. In particular, we consider the same relations between the bond nonlinearity coefficients as those in Sobirov et al. (2010) and study the shrinking of the fat graph into the metric graph keeping such relations. Initial conditions for the NLSE on fat graph are taken as quasi 1D solitons. By solving the NLSE on fat graphs we find that in the shrinking limit such fat graphs reproduce the reflectionless regime of transport studied in Sobirov et al. (2010), i.e., the vertex transmission becomes ballistic.
The linear Schrödinger equation on fat graphs was the subject of extensive study during the past decade (see, e.g. Rubinstein and Schatzman (2001); Kuchment and Zeng (2001); Exner and Post (2005); Post (2006); Exner and Post (2007, 2009); Exner (2011); Exner and Post (2013); Cacciapuoti and Exner (2007); Post (2012); Ruedenberg and Scherr (1953); Molchanov and Vainberg (2007, 2010); Dell’Antonio and Costa (2010); Kosugi (2002); Grieser (2008a)). The first treatment of particle transport on fat graphs dates back to Rudenberg and Scherr Ruedenberg and Scherr (1953), who used a Green function based heuristic approach. A pioneering study of particle transport in fat networks comes from the paper by Mehran Mehran (1978) on particle scattering in microstrip bends and junctions, where theoretical results on reflection and transmission are compared with experimental data. However, the dependence of the scattering on the bond thickness and the shrinking limit is not considered in Mehran (1978).
The main problem to be solved in the treatment of the Schrödinger equation on fat graphs is reproducing of vertex coupling rules in the shrinking limit, i.e., when the fat graph shrinks to the metric graph. In case of metric graphs, ”gluing” conditions, or vertex coupling rules, are needed to ensure selfadjointness of the Schrödinger equation. The most important example of a vertex coupling is the Kirchhoff condition. For fat graphs there are no such coupling rules; they only appear in the shrinking limit, and their form depends on specifics of the fat graph, for example on the boundary conditions imposed at the lateral boundary. For Neumann boundary conditions, the resulting vertex coupling is the Kirchhoff condition, as was shown in Rubinstein and Schatzman (2001); Kuchment and Zeng (2001), who study convergence of the eigenvalue spectrum of the Schrödinger equation, and in a series of papers by Exner and Post Exner and Post (2005)Exner and Post (2013), who study various aspects of the Schrödinger equation with Neumann boundary conditions (including transport, resonances and magnetic field effects). The vertex couplings obtained in the shrinking limit of the Schrödinger equation on the fat graph with Dirichlet and other boundary conditions were obtained in Molchanov and Vainberg (2007); Grieser (2008a). Recent studies of the linear Schrödinger equation on fat graphs focused on the inverse problem of finding a suitable fat graph problem which reproduces a given coupling rule in the shrinking limit Cacciapuoti and Exner (2007). Further references on linear Schrodinger equation on fat graphs are Exner et al. (1988); Exner and Seba (1989); Kostrykin and Schrader (1999); Kottos and Smilansky (1999); Gnutzmann and Smilansky (2006); Gnutzmann et al. (2010); Dell’Antonio and Costa (2010); Molchanov and Vainberg (2010); Exner (2011); Exner and Post (2013), and the reviews Grieser (2008b); Post (2012). All the above results have been limited to linear and stationary cases, and spectral results. Related problems also have a long history in (nonlinear) PDEs, see Raugel (1995) and the references therein, where however the focus is on dissipative systems, and on damped wave equations.
The case of the NLSE on fat graphs is much more complicated than the linear case. Therefore one may expect that the treatment of the NLSE with the same success as for the linear problem is not possible. To our knowledge, the only work dealing with nonlinearities on fat graphs is by Kosugi Kosugi (2002), who considers semilinear elliptic problems and shows convergence of solutions towards solutions of the metric graph problem. However, for problems such as soliton transport, scattering and interaction with external potentials which are described by timedependent evolution equations on fat graphs, we have to rely to a large extent on numerics.
In this paper, using the numerical solution of the NLSE on fat graph we explore the dependence of soliton transmission and reflection at the fat graph vertex on the bond thickness and the angle between the bonds. It is organized as follows. In the next section we give detailed formulations of the problems both for fat and metric graphs. Section III presents numerical (soliton) solutions of the NLSE on fat graphs, and analysis of the soliton reflection at the graph vertex in the shrinking limit, including the dependence of reflection coefficient on the angle between the graph bonds. Section IV gives conclusions, while the Appendix contains some details of the numerics.
Ii The NLSE on metric and fat graphs
Consider the nonlinear Schrödinger equation
(1) 
on a metric star graph with edges , and nonlinearity coefficients . The graph is assumed to have semiinfinite bonds , , but the main part of our analysis will be numerical, for which we assume finite lengths of bonds, with coordinates , , and homogeneous Dirichlet boundary conditions at . Furthermore, we assume that the solutions, obey the vertex (at ) conditions
(2) 
with parameters , where it is understood that () denote the derivatives from the left (right). In the following we call Eqs.(1) and (2) problem (P).
Soliton solutions of the problem (P) that propagate without reflection (i.e., ballistically) were obtained analytically in Sobirov et al. (2010) for the special case when the nonlinearity coefficients satisfy the relation
(3) 
These solutions have, after properly identifying with on the form
(4) 
with free parameters amplitude , speed (wavenumber ), and reference position . Fig.2 presents amplitudes, for Kirchhoff boundary conditions () and for the boundary conditions given by Eq.(2).
The vertex boundary conditions given by (2) are one possibility to make the linear part of (1) skewadjoint. The problem (P) conserves the norm and the Hamiltonian given by
(5)  
(6) 
It is a question of normalization to set
(7) 
which leaves 4 parameters for (P), and, of course, the choice of the initial conditions.
Our goal is to compare exact and numerical solutions of (P) with the numerical solutions of an associated NLSE on a fat graph presented in Fig. 1, i.e.,
(8) 
where , and consists of a “vertex–region” of diameter , and tubes around , see Fig. 1. In the following Eq.(8) will be called the problem (P). We also use the notation for .
It is clear that different versions of are possible. Here we choose to give the following 5 parameters to not a priori present in Eq.(1):

the angles between the bonds and and the –axis,

the widths of the different bonds.
In the numerical calculations we impose homogeneous Dirichlet boundary conditions (DBC) for both, (P) and (P), at the “ends” of bonds, and for (P) homogeneous Neumann boundary conditions (NBC) everywhere else. As our simulations will run on time–scales where the solitons will be well separated from the ends of the bonds, we could as well pose NBC there. Also note that strictly speaking (4) is not a solution over the finite graph, but it is exponentially small at the ends of the bonds.
We take constant on bond and with suitable jumps near . Furthermore, we set
, and  (9) 
and write for fixed , For definiteness we choose
(10) 
and thus . Motivated by as , corresponding to on we define the scaled norms
(11)  
(12) 
Then is conserved for (8), and the indicate how much “mass” is in the different bonds.
For the linear problem it is known, Exner and Post (2009), that under the scaling
(13) 
the vertex conditions (2) appear in the limit . Then, at least formally, we can expect (P) as a “limit” of (P) if
(14) 
If (or ), then the boundary conditions (2) gives jumps from to (resp. ) at the vertex. This, however, is merely a question of scaling. For instance, setting (cf. (13)), we obtain
(15) 
i.e., continuity at the vertex, where , as in (14). The scaling given by Eqs.(1),(2) is more custom Sobirov et al. (2010); Exner and Post (2009) than (15), and therefore we stick to (1),(2) as the “limit problem”. Note that the angles of the fat graph do not appear in (P).
We expect that for solutions of (P) behave like with being the solutions of (P), i.e., are constant in transverse direction on each bond , with width . Therefore, from Eqs. (12) and (13) we expect
(16) 
In the numerical calculations, in addition to we explore the following functions (dropping the dependence on parameters and ):
(17)  
(18)  
(maximal 
Here is the extension of to , constant in transverse direction, and for we either use the explicit formula (4) if (3) holds, or numerics for (P) if not. Note that (18) ignores phase differences between and , as these are less important from the viewpoint of applications.
Iii Soliton transport in fat graphs
The main practically important problem in the context of wave propagation in branched systems is energy and information transport via solitary waves. The dependence of the soliton dynamics on the topology of a network makes such systems attractive from the viewpoint of tunable transport in low dimensional optical, thermal and electronic devices. Therefore, the treatment of the problems (P) and (P) from the viewpoint of vertex soliton transmission is of importance. Our main purpose is to compare propagation of solitons in with that in in particular we are interested to study the “lift” the earlier results Sobirov et al. (2010) from to . Transition from two to onedimensional wave motion in the shrinking limit is of special importance for this analysis.
In a typical simulation, for (P) we use solitontype initial condition given as
(19) 
where and are chosen in such a way that is very close to . Similarly, for (P) we choose
(20) 
i.e., we extend the initial conditions (19) trivially in transverse direction. We then run both, (P) and (P) until some final time such that the solitons launched by (19) and (20), respectively, have interacted with the vertex, and have been reflected or transmitted sufficiently far into the bonds (see the appendix for the numerical methods used). Our main solution diagnostics will be the time dependent norms , the amplitudes , the distances , and the reflection coefficients defined below.
For definiteness, we consider as the “incoming” bond and as “outgoing” ones. In Fig. 3 solutions of the problem (P) for the Kirchhoff boundary conditions are presented for the case of a “relatively fat” graph (), while Fig. 4 shows the plots of the corresponding norms and amplitudes for the simulation for (P) in Fig. 3 (Kirchhoff case), together with the respective quantities for (P). At this relatively large there is a significant difference between (P) and (P).
In the following we focus on soliton reflection and transmission in the shrinking limit, for the ”ballistic” boundary conditions given by (3) on (P). In Fig.5 and 6 we plot the diagnostics defined above for different on an otherwise fixed graph fulfilling (3), i.e., for the ballistic case. As , the amplitudes and masses in the different bonds get close to the metric graph case, and also the (numerical) wave functions as a whole converge to the ones on the metric graph, with one small qualification: While the main mismatches between (P) and (P) result from reflection and position shifts of the incoming soliton during interaction with the vertex around , already for , i.e., before interaction of the soliton with the vertex, there is a small linear growth of , i.e., of the amplitude mismatch in the incoming bond. This is not a property of the fat graph itself, but related to the fact that it is difficult to accurately resolve the speed of the soliton numerically. In other words, for small , a significant part of mismatch between our (numerical) fat graph solution and the (analytical) metric graph solution from Sobirov et al. (2010) is not due to the behaviour at the vertex, but due to an error in (numerical) soliton speed, which results in a position mismatch growing in time (see the Appendix for further discussion). However, noting the different scales in panels (a4),(b3) and (c3) strongly indicates the convergence of the (P) wave function to the (P) wave function in (modulo phases), uniformly on bounded time intervals.
From the viewpoint of practical applications, probably the most important question is how much of an incoming soliton is reflected resp. transmitted in the vertex region of a fat graph. To display this in a concise way, for (P) we define the reflection and transmission coefficients
(21)  
(22) 
where again we dropped the dependence on parameters and here, but will plot as functions of some parameters below. Thus, e.g., (and thus also ) means zero reflection of an incoming soliton at the vertex, while, e.g., means that all of the “mass” was transmitted to bond two. These extreme cases of course do not occur, but the goal is, e.g., to tune . The corresponding quantities for (P) are defined as
(23) 
and the transmission formula (3) means, that in the limit of infinite bonds and of .
Fig.7 presents plots of the reflections as a function of different parameters such as bond thickness, , angle between the bonds, and coefficient, . Of course, the ballistic regime is an idealization, and in order to quantify the reflections when we perturb it, in Fig. 7(a) we study the dependence of and on the mismatch . As we have chosen , , in (a) we fix and let vary between and , such that again corresponds to the minimal reflection case. Clearly, the minima of and are attained close to , and these functions are somewhat steep. Similar graphs were obtained in other geometries (e.g., angles), and thus in applications it appears desireable to move as close as possible to the ballistic case by, e.g., fine tuning the widths of the bonds. In Fig. 7(b) the vertex reflection coefficients (both for norm and amplitude) are plotted as functions of the graph thickness for the case described by (3). The limit again shows a rather smooth transition from the “scattering” to ballistic regime.Figure 7(c) presents the dependence of on the graph thickness and the angles . Even though the angles do not appear in the limit (P), at finite they of course play an important role. Ballistic transport through the vertex occurs in the shrinking limit as well as in the limit of small angles.
Besides the equal angle case , we checked a variety of other configurations with , for various between and . The results remain qualitatively similar to Figs. 5–7, i.e., in the ballistic case the reflection coefficients vanish as , and as above the (P) wave functions converge to the independent wave functions of (P). As the convergence for is clearly linear, an interesting question is how to choose a first order in correction of the fat graph geometry or NLSE coefficients that minimizes also for finite .
An important issue for particle and wave transport in fat graphs is the dependence of the scattering on initial soliton velocity and amplitude. In Fig.8 reflection coefficients are plotted as functions of bond thickness for different initial velocities (a) and amplitudes (b). The dependence of reflection on initial data is significant for fat graphs, with, e.g., less reflection for slower and longer waves, as should be expected. However, in the shrinking limit the reflections vanish in all cases considered.
Iv Conclusions
We studied soliton transport in tube like networks modeled by the timedependent NLSE on fat graphs, i.e. graphs with finite bond thickness. We numerically solved the NLSE on fat graphs for different values of thickness, and studied behavior of solutions and vertex reflection coefficients in the shrinking limit. It is found that in the shrinking limit solutions of the NLSE on fat graphs converge to those on the associated metric graphs, and hence the conditions (3) for reflectionless transport also work on fat graphs with small . The dependence of the vertex reflection coefficient on the bond thickness and on the angle between the bonds of fat graph is also studied.
At this point it is not clear in which norms we can expect or analytically show convergence of solutions of (P) to solutions of (P), as . First, following Kosugi (2002) this will be discussed for the stationary case, including some potentials at the vertex in order to have nontrivial stationary solutions for the fat graph and the metric graph, cf. Adami et al. (2012); Sabirov et al. (2013). An important point in the study of wave(particle) dynamics in fat graphs is the definition of the fat graph thickness at which one can neglect transverse motion and consider the system as onedimensional. The above treatment allows us to define such a regime. However, the transition from two to one dimensional motion is rather smooth and there is no critical value of the bond thickness at which a ”jump” from the fat to the metric graph occurs. In any case, we believe that our numerical results should be considered as a first step in the way for the study of particle and wave transport described by nonlinear evolution equations on fat graphs. In addition, can be useful for further analytical studies of the NLSE on such graphs.
Acknowledgement
We thanks Pavel Exner and Riccardo Adami for discussions and useful comments. This work is supported by a grant of the Volkswagen Foundation.
Appendix: Details of the numerical approach
We discretize (P) by second order spatial finite differences (FD) and denote , , , , such that, e.g., . Moreover, we set
, , ,  
and , , and . 
The vertex conditions then are , Using onesided FD for and we have
(24) 
which expresses and hence in terms of . The resulting can be best expressed by a matrix vector multiplication . The scheme differs from the one in Sobirov et al. (2010), where the PDE is extended up to and including the vertex from the left, which works well to discretize the reflectionsless solutions (4) in case of (3), but it introduces an asymmetry between the bonds not present in (P).
To integrate the resulting ODEs , where with obvious meaning, we use an explicit scheme with stepsize in , namely
(25) 
where and similar for and . For this conserves with high accuracy, and also .
To simulate (P) we write it as a 2component real system for where . We set up and discretize the domain using routines from pde2path Uecker et al. (2014) which are based on the FEM from the Matlab pdetoolbox. For efficiency it is quite useful to apply some local mesh refinement near the vertex. We typically work with meshes of 500020000 triangles, requiring a maximal mesh size of before local refinement. Eq. (8) then translates into the system of ODEs
(26) 
where is the mass matrix, is the stiffness matrix, and is the FEM nonlinearity. For the time integration of (26) we use a semilinear trapez rule, i.e., setting , with typically to ,
(27) 
Over relevant timescales (27) conserves (the discretized version of) from (12) reasonably well, see, e.g., Fig. 4, 5(a2), and 6(b1), (c1), but, as already indicated, dependent on the discretizations there are some slightly more significant errors in the numerical soliton speed. To quantify this, we used (27) to propagate a soliton of amplitude and speed on a straight bond from to position at for various time steps and meshsizes , and calculated the error
where and denote the numerical and the exact solution, respectively. In reasonable ranges turns out to be roughly a linear function of . See Fig. 9(a) for as an example, while (b) shows for and thus illustrates that (27) propagates the soliton too quickly.
We also tried the relaxation scheme from Besse (2004) which conserves slightly better, but becomes computationally much slower, mainly since one can no longer prefactorize . On the other hand, the stability requirements for explicit schemes like (25) become prohibitive for fine meshes near the vertex. For (27), typical computation times for the propagation of a solitary wave through the network are on the order of 1 minute.
References
 Z. Sobirov, D. Matrasulov, K. Sabirov, S. Sawada, and K. Nakamura, Phys. Rev. E 81, 066602 (2010).
 R. Cascaval and C. Hunter, Libertas Math. 30, 8598 (2010).
 V. Banica and L. Ignat, J. Math. Phys. 52, 083703 (2011).
 R. Adami, C. Cacciapuoti, D. Finco, and D. Noja, Rev. Math. Phys. 23, 409 (2011).
 R. Adami, C. Cacciapuoti, D. Finco, and D. Noja, Europhys. Letters (2012).
 R. Adami, C. Cacciapuoti, D. Finco, and D. Noja, J. Phys. A: Math. Theor. 45, 192001 (2012).
 R. Adami, D. Noja, and C. Ortoleva, J.Math. Phys. 54, 013501 (2013).
 D. Noja, Phil. Trans. R. Soc. A 372 (2014).
 S. Gnutzmann, U. Smilansky, and S. Derevyanko, Phys. Rev. A 83, 033831 (2011).
 K. Sabirov, Z. Sobirov, D. Babajanov, and D. Matrasulov, Phys. Letters A 377, 860865 (2013).
 P. Leboeuf and N. Pavloff, Phys. Rev. A 64, 033602 (2001).
 T. Paul, M. Hartung, K. Richter, and P. Schlagheck, Phys. Rev. A p. 063605 (2007).
 I. N. de Oliveira, F. A. B. F. de Moura, M. L. Lyra, J. S. Andrade, Jr., and E. L. Albuquerque, Phys. Rev. E 81, 030104(R) (2010).
 S. Yomosa, Phys. Rev. A 27, 2120 (1983).
 L. V. Yakushevich, A. V. Savin, and L. I. Manevitch, Phys. Rev. E 66, 016614 (2002).
 R. Burioni, D. Cassi, P. Sodano, A. Trombettoni, and A. Vezzani, Phys. Rev. E 73, 066624 (2006).
 A. Markowsky and N. Schopohl, Phys. Rev. A 89, 013622 (2014).
 V. E. Demidov, J. Jersch, S. O. Demokritov, K. Rott, P. Krzysteczko, and G. Reiss, Phys. Rev. B 79, 054417 (2009).
 N. R. Biggs, Wave Motion 49, 24–33 (2012).
 J. Rubinstein and M. Schatzman, Arch. Rat. Mech. Anal. 160, 271 (2001).
 P. Kuchment and H. Zeng, J. Math. Anal. and Appl. 258, 671 (2001).
 P. Exner and O. Post, J. Geom. Phys. 54, 77 (2005).
 O. Post, Ann. Henri Poincaré 7, 933 (2006).
 P. Exner and O. Post, J. Math. Phys. 48, 092104+43 (2007).
 P. Exner and O. Post, J. Phys. A, Math. Theor. 42, 22 (2009).
 P. Exner, in Mathematics in science and technology. Mathematical methods, models and algorithms in science and technology. Proceedings of the satellite conference of the ICM 2010, New Delhi, India, August 14–17, 2010. (2011), pp. 71–92.
 P. Exner and O. Post, Commun. Math. Phys. 322, 207 (2013).
 C. Cacciapuoti and P. Exner, J. Phys. A 40, F511 (2007).
 O. Post, Spectral analysis on graphlike spaces. (Berlin: Springer, 2012).
 K. Ruedenberg and C. W. Scherr, J. Chem. Phys. 21, 1565 (1953).
 S. Molchanov and B. Vainberg, Commun. Math. Phys. 273, 533 (2007).
 S. Molchanov and B. Vainberg, in Integral methods in science and engineering. IMSE 2008, Santander, Spain (Basel: Birkhäuser, 2010), pp. 255–278.
 G. Dell’Antonio and E. Costa, J. Phys. A, Math. Theor. 43, 23 (2010).
 S. Kosugi, Journal of Differential Equations 183, 165 (2002).
 D. Grieser, Proc. Lond. Math. Soc. 97, 718 (2008a).
 R. Mehran, IEEE Transactions on Microwave Theory and Techniques 26, 4005 (1978).
 P. Exner, P. Seba, and P. Stovicek, J. Phys. A 21, 4009 (1988).
 P. Exner and P. Seba, Rep. Math. Phys. 28, 7 (1989).
 V. Kostrykin and R. Schrader, Rev. Math. Phys. 11, 187 (1999).
 T. Kottos and U. Smilansky, Ann. Phys. 274, 76 (1999).
 S. Gnutzmann and U. Smilansky, Adv. Phys. 55, 527 (2006).
 S. Gnutzmann, J. Keating, and F. Piotet, Ann. Phys. 325, 2595–2640 (2010).
 D. Grieser, in Analysis on graphs and its applications. Isaac Newton Institute, Cambridge, UK, January 8–June 29, 2007 (AMS, 2008b), pp. 565–593.
 G. Raugel, in Dynamical systems (Montecatini Terme, 1994) (Springer, Berlin, 1995), vol. 1609 of Lecture Notes in Math., pp. 208–315.
 H. Uecker, D. Wetzel, and J. Rademacher, Numerical Mathematics: Theory Methods Applications 7, 58106 (2014).
 C. Besse, SIAM J. Num. Anal. 42, 934 (2004).