# Easy plane baby skyrmions

## Abstract

The baby Skyrme model is studied with a novel choice of potential, . This “easy plane” potential vanishes at the equator of the target two-sphere. Hence, in contrast to previously studied cases, the boundary value of the field breaks the residual internal symmetry of the model. Consequently, even the unit charge Skyrmion has only discrete symmetry and consists of a bound state of two half lumps. A model of long-range inter-Skyrmion forces is developed wherein a unit Skyrmion is pictured as a single scalar dipole inducing a massless scalar field tangential to the vacuum manifold. This model has the interesting feature that the two-Skyrmion interaction energy depends only on the average orientation of the dipoles relative to the line joining them. Its qualitative predictions are confirmed by numerical simulations. Global energy minimizers of charges are found numerically. Up to charge , the minimizers have half lumps positioned at the vertices of a regular -gon. For charges , rectangular or distorted rectangular arrays of half lumps are preferred, as close to square as possible.

###### pacs:

11.27.+d, 05.45.Yv, 11.10.Lm, 11.10.Jj## I Introduction

The Skyrme model, and its smaller relative, the baby Skyrme model, have received much attention during the last few decades. The models are very interesting from the mathematical physics point of view: they contain topologically nontrivial field configurations called topological solitons and the Skyrme model arises as a large- limit of QCD Witten (1983a, b). The baby Skyrme model is the two-dimensional analog of the full model and its main interests lie in its easier tractability, the possibility of using the baby Skyrme model results as guidelines in investigating the full Skyrme model and some condensed matter contexts, where the baby Skyrme model itself is of physical relevance Sondhi et al. (1993). The baby Skyrme model was investigated intensively during the 1990s but interest in the model has recently received a new boost Hen and Karliner (2008a, b); Hen and Karliner (2009); Hen and Karliner (2008c, d); Adam et al. (2009, 2010); Speight (2010), apparently due to the prevalence of baby Skyrmions in various condensed matter systems Rößler et al. (2006), such as FeCoSi Yu et al. (2010), MnSi Mühlbauer et al. (2009); Neubauer et al. (2009) and quantum Hall systems (for a review, see Ezawa and Tsitsishvili (2009) and references therein), despite the fact that these are not modeled by the baby Skyrme model. It is against this background that we investigate the baby Skyrme model.

The behavior of the solitons of the baby Skyrme model depends crucially on the choice of potential. For some potentials, there will be no smooth solutions at all; for some, like that studied in Leese et al. (1990), there will be exact solutions for the equations of motion and for the rest, such as here, numerical methods must be used. So far, only cases where the potential and boundary condition preserve the symmetry of the Lagrangian have been studied in detail. Here we study the easy plane potential, which breaks symmetry to , and the boundary condition further breaks the symmetry to . This endows the model with one massive and one massless mode and one expects significantly different behavior from symmetric cases.

We use numerical methods to demonstrate that the model has a rich static energy structure, including local minima and that multi-Skyrmions are stable in this model. We also develop an approximate way to estimate the long-range interaction of two unit baby Skyrmions and verify this numerically. We identify some open questions and present a conjecture about the energy of the global minimum, when the charge approaches infinity.

## Ii The model

The Lagrangian density of the model, defined on , is

(1) | ||||

(2) |

where the first terms are called the Dirichlet, the Skyrme and the potential terms, respectively and are free parameters (coupling constants) of the model. The subscripts are chosen to match the number of derivatives in the term and will prove useful later.

Choosing the usual metric yields the static energy density

(3) |

If one requires the solitons to have finite energy, it is necessary to impose the boundary condition

(4) | |||

and choose from the vacuum manifold of the model, i.e. such that | |||

(5) |

As usual, Eq. (4) allows us to consider the field as map , which gives the system its interesting topological properties, namely, its fields are classified topologically by the homotopy group . The homotopy invariant of the field is called the (topological) degree and is given by the expression

(6) |

Furthermore, one can expect the system to support topological solitons based on two factors. First, Derrick’s theorem of nonexistence does not apply since the terms zeroth and fourth order in derivatives behave in opposite ways upon uniform scaling (the second order term is scale invariant), thus allowing solutions with a preferred size. The second reason is the existence of a Bogomol’nyi bound for the (static) energy which depends on the degree of the field, given by Speight (2010)

(7) |

where is the average value of . Unlike in some models, the bound cannot generically be saturated in this model.

In this work, we take and . This potential has the exceptional effect of allowing the breaking of the symmetry of the energy almost completely. While the energy density exhibits symmetry (from the potential), the choice of boundary condition breaks this down to . Hence, even the Skyrmion with unit topological charge will have only discrete symmetry.

We seek static solutions in various homotopy classes and investigate the forces between some of these. In order to find a static solution, we numerically minimize the static energy of an initial field configuration. Some initial configurations will be constructed from existing minima, but those that are obtained from an exact expression are taken to be

(8) | |||

(9) | |||

which for gives | |||

(10) |

where is simply a scale factor to adjust the size of the initial configuration. We also added random ripples to destroy any discrete symmetries left over from the initial configuration. It was observed that this is usually unnecessary, but some minimizers take different orientations (but keep their shape) when the random ripples are omitted.

In all that follows, we set .

## Iii The unit skyrmion

The lattice approximation of derivatives is achieved by using a simple forward differencing scheme. The discrete energy functional is then minimized using gradient based methods.

The gradients are computed from the discrete total energy in a straightforward manner. The energy minimization is then achieved with the TAOBenson et al. (2007) and PETScBalay et al. (2009, 2008, 1997) parallel numerical libraries. The libraries provide several different minimization schemes, of which we have chosen the (in this case) fastest algorithm, the limited memory quasi-Newton algorithm (also called a variable metric algorithm) with BFGSBroyden (1970); Fletcher (1970); Goldfarb (1970); Shanno (1970) formula for Hessian approximations. We consider the minimization to have converged, when the -norm of the gradient is less than . This may sound like a very strong requirement, but it turns out to be necessary. Our initial choice of works for most cases, but when the results were checked with a the stronger requirement, they experienced significant changes and decreases in energy after attaining the looser condition. Using a still stronger criterion did not change the energies or shapes of the minimizers any further.

The accuracy of the numerical scheme was tested against the known exact solution of Peyrard et al. (1992) for the potential . We find that the accuracy of the numerical scheme is heavily dependent on the size of the lattice, due to the energy density of the Dirichlet term being slowly decaying. When in a large enough lattice, the difference in energies between our minimizer and the exact solution is less than , so we expect similar accuracy in the present model as well, since it has similar decay properties.

We always use a square lattice and will refer to the number of lattice points along its side by a number . Starting with (10) on a lattice with lattice constant , the corresponding minimum energy density is depicted in Fig. (a)a. The broken rotational symmetry is immediately evident in the form of the two red peaks. A view of how the field itself behaves is shown in (b)b and an immediate observation is that the half-plane gets mapped to the upper hemisphere of the target and the half-plane gets mapped to the lower hemisphere. This, combined with the fact that each individual peak in energy density strongly resembles a model lump, leads us to call the solution on a half-plane a half lump.

Note that the energy density image is zoomed in on the central area of the numerical lattice, where almost all the energy density is concentrated. Obviously, the whole computation could have been done in a much smaller lattice, but we opted not to find the optimal lattice sizes for each and every configuration due to the fact that the amount of computer resources saved would have been small. We also confirmed that the resulting field configuration is not dependent on either the value of the lattice constant or “physical” size of the computational box by first repeating the computation with lattice constant and and then repeating again with and .

It is instructive to look at the fields themselves, in particular, their asymptotics at large . Since as , we must have

(11) |

at large , where are small. Substituting this into and keeping only leading order terms yields the Lagrangian density of an uncoupled pair of real scalar fields,

(12) |

one massless () and the other of mass (). One expects the large field of a soliton to be well approximated by a static solution of this linear model. Since , should wind once around as one traverses a circle at large . This leads one to predict that

(13) |

at large , where are unknown constants.

Numerical evidence in favor of the conjectured asymptotics for is presented in Fig. 6, which shows plots of (LABEL:sub@fig:b1masslessfit) and (LABEL:sub@fig:b1massivefit) against for an increasing sequence of values of . By fitting sinusoidal curves, we find that and . Note that, at large , decays like , so the massive field decays much more quickly than the massless field . This has important consequences for the long-range forces between Skyrmions, as we shall see in the next section.

Having found the solution, it is natural to ask whether there are multisolitons, that is localized solutions with . Given the linear dependence of the energy bound on , Eq. (6), it is conceivable that all these would consist of copies of the solution at infinite separation. Indeed, this is what happens for the potential , where the localized solutions are unstable against decay to infinitely separated solutions with Leese et al. (1990).

In that model (the so-called “holomorphic” model), the force between a pair of unit charge skyrmions turns out to be repulsive, whatever their separation and relative orientation. By contrast, we will see that in the current model (as for several other choices of potential) there is a choice of relative orientation for which solitons attract one another at long range (the so-called attractive channel), making it plausible that spatially localized solutions of every charge exist.

## Iv Interskyrmion forces

Intersoliton forces in the baby Skyrme model with an easy-axis potential were analyzed in detail by Schroers et al in Piette et al. (1995), and their methods readily adapt to the current setting. We shall assume that the dominant interaction at long range is mediated by the massless field and ignore the contribution of the faster decaying field . The key observation is that the asymptotic field

(14) |

coincides with the static solution of the linearized model

(15) |

in the presence of an appropriate point source at the origin, in this case

(16) |

where denotes the two-dimensional Dirac delta distribution. To see this, note that the elementary kernel Choquet-Bruhat et al. (1982) for Laplace’s equation on is , that is,

(17) |

Source (16) induces the field in standard orientation and position. The field rotated (spatially) through angle and translated to position is induced by the source

(18) |

where . Physically, one interprets this as the massless scalar field induced by a dipole source of moment at position . It is important to realize that the dipole moment which replicates the soliton’s asymptotic field points orthogonally to the line connecting its constituent half lumps. One should not, therefore, think of the half lumps as being scalar monopoles of opposite charge, generating a dipole moment (recall that the half lumps are structures in the field, while the long-range field is dominated by ).

Now consider the situation where two solitons with orientations are placed at positions with large. The interaction energy that they experience should, on physical grounds, coincide asymptotically with the interaction energy of the corresponding scalar dipoles in the linear theory, that is,

(19) |

where and is the field induced by dipole at , that is

(20) |

Hence

(21) | |||||

where , the unit vector directed from dipole 2 to dipole 1. It is interesting to note that this is precisely minus the interaction energy of a pair of electric dipoles in two-dimensional electromagnetism, which can be understood by noting that, in scalar field theory, like monopole charges attract, while opposite charges repel. In the case where the solitons are placed symmetrically about the origin on the axis, that is, , , one sees that

(22) |

This formula is rather remarkable. It has the interesting property that it depends only on the average orientation of the dipoles, and is independent of their phase difference . This is precisely the opposite of the situation found by Schroers et al in the easy-axis model Piette et al. (1995),

(23) |

which highlights a fundamental difference between the two models. The unit soliton in the easy-axis model is axially symmetric, and spatially rotating the soliton is identical to internally rotating it (about the axis ). Consequently, the interaction energy in the easy-axis case must be invariant under , and hence can depend only on . In the easy plane model, however, the unit soliton has no rotational symmetry. One cannot rotate it internally about (the fields so obtained do not solve the field equation), and spatially rotating the soliton produces a genuinely physically distinct solution (e.g. its energy density changes). So there is no reason to expect to be independent of . That it depends only on is more surprising, and is a basic fact about dipole interactions in two dimensions. In practice, it seems likely that higher-order corrections will break the independence of on .

The most important consequence of formula (22) is that it predicts long-range attraction between unit solitons if they are appropriately oriented relative to one another, leading one to expect there to exist stable, static, spatially localized solutions in all homotopy classes. In particular, soliton pairs initially oriented as in Fig. 7(a) or 7(b) are predicted to attract one another, while those oriented as in Fig. 7(c) or 7(d) are predicted to repel. These predictions were vindicated by the numerical simulations described in the next section.

## V Multiskyrmions

In this section, we report the numerical results on multi-Skyrmions, confirming the validity of the dipole approximation, and give examples of the very complicated landscape of local minima of higher multi-Skyrmions. We split this section in two according to two different “branches” in the shapes of the minima. Before going further, it is worth noting several general observations.

First, the dipole approximation was confirmed for pairs of solitons. To confirm this, we used three methods: the usual BFGS method, simple gradient flow and also a Newton flow second order dynamics, where we set and . They all agree that two unit Skyrmions starting in the attractive channel coalesce into a multisoliton with an energy lower than twice that of the unit Skyrmion. The resulting two-Skyrmion is depicted in Fig. 16LABEL:sub@fig:b2. However, confirming the repulsive channel is slightly more challenging. First, the computational box needs to be large enough, otherwise the boundaries will exert pressure on the solitons and they react by rotating into the attractive channel and coalescing. Second, the solitons need to be placed exactly symmetrically on the lattice, lest they again rotate into the attractive channel and coalesce. When set up correctly, however, placing two Skyrmions in the repulsive channel results in them receding. Thus we can confirm the qualitative correctness of the dipole approximation.

Second, all the results presented below are done in a standard lattice with and , except the case , which does not fit into the smaller lattice at all. For those values of , where multiple minima were found, two additional methods were used for their initial configurations. The first method is a simple cut and paste of other solitons, for example, two -solitons were put side by side to produce a -soliton. The other method was to start with which resulted in a different minimum than . This is due to boundary pressure: if was doubled, the other result reemerged. When this happened, the solution was made sparser by a factor of 4 to put it into the standard lattice, where it was then minimized again to get a result which can be compared with the others. Where several initial states produced the same minimizer, we only mention one of the initial states.

### v.1 Charges

Minimal energy Skyrmions of charges in the interval found using our numerical scheme form regular polygons, where there are half lumps at the vertices. The energies of these minimizers are always lower than , where is the energy of the unit Skyrmion, so there is good reason to believe they are stable in accordance with the above dipole picture. It should be noted, however, that the symmetry of these minimizers is the same as that of the initial data . Therefore, one is led to suspect that the random ripples added to may not be enough to break the symmetry. It is therefore necessary to introduce a stronger symmetry break. For , the computations done to confirm the dipole approximation provide this break, and still relax to the same square final configuration. The obvious way to test the higher charges was therefore to place two lower charge solitons side by side and minimize their energy.

Indeed, for charges 5 and 6, there are minima, which were only found using this spliced initial configuration. However, they have higher energy than the respective polygons, and therefore represent a local minimum only. For the charges and , various different initial configurations were tried, but they always relaxed to the respective polygon. We present the minimizer of as an example of the polygon shaped minimizers in Fig. 16LABEL:sub@fig:b6.

### v.2 Charges

The -gon form of the global minimizers only persists up to after which, the shape of the minimizer is completely different and there is no obvious rule for the shape. We will briefly discuss these minimizers below, but first it should be noted that the symmetric (plus random ripples) initial state still produces a polygon shape local minimizer for all charges except 13 and 14, but this is no longer a global minimum. The polygon exists for charges 13 and 14, also, but was only obtained by starting the minimization process with and embedding the result in the standard lattice. The polygons are, however, quite stable and the energy differences between them and the global minima are very small. For example, we squeezed some of the polygons to half the size in the direction and they relaxed back to the polygon shape instead of the global minimum.

The solutions obtained from lower charge solitons set side by side and starting the minimization in a lattice, produce a veritable zoo of different minimizers, of which all but one have lower energies than the -gons, the exception being charge 14. Because of the amount of work involved in patching the solitons together and the total number of computations already performed, we did not attempt an exhaustive search with side-by-side-type initial configurations of all charges.

The first nonpolygon minimizer occurs at and looks like a slightly distorted rectangle of half-lumps. It is depicted in 16LABEL:sub@fig:b7 and was obtained by putting a soliton next to a one. A similar minimizer is also found for using a soliton beside a soliton as the initial configuration. These are global minima, within the limits of our numerics.

At the half lumps of the minimum energy configuration arrange themselves into a square lattice of half lumps, as seen in Fig. 16LABEL:sub@fig:b8. The initial configuration consisted of four solitons at the vertices of a square. Similar square lattice configurations are also found for and , shown in Figs. 16LABEL:sub@fig:b18 and 16LABEL:sub@fig:b32. The initial configurations were as follows. For , two and two solitons were put at the opposite vertices of a square and for , four solitons were placed at the vertices of a square. It is interesting to note that the solution can also be thought of as square lattice of half lumps, providing us with a tantalizing opportunity to conjecture that all will fall into similar square lattices. While these are all possible square lattices for , this raises the question of whether there are other symmetric minimizers and what are their symmetries, if any.

It seems that using various spliced initial configurations, the square symmetry can be relaxed to produce other rectangular shapes. We have already seen two of these at and , which are not global minima and at and , which are the global minima for these charges. For all four cases, the short edge was the shortest possible, but for a more squarelike rectangular lattice of , , and half lumps, respectively, exists. These are all global minima; they were all obtained using spliced initial configurations.

The global minima for and have the shape of a distorted rectangle, where the distortion consists of an extra two half lumps in one corner. These were both achieved by using initial states (as always, these were later coarsened and embedded into the standard lattice). Of these, we show the minimizer in Fig. 16LABEL:sub@fig:b13 as an example.

As the charge increases, finding the global minimum – or what can reasonably be believed to be one – becomes increasingly difficult. It is likely that the (possibly distorted) rectangular shapes prevail as global minima for charges above , but we have not been able to ascertain this. As an example of how difficult it is to find the global minimum, consider charge . There is a local minimum of the distorted rectangle type with two extra half lumps in one corner. To converge into this configuration, an initial configuration of charges and is needed, but it is not enough to just put these side by side. They must be placed at very precise locations in the lattice. Although this arrangement is natural, place the solution outside the solution along one of its diagonals, it is very difficult to guess the correct arrangement for higher charges. Moreover, at higher charges, there is an ever increasing number of different combinations to choose from.

In addition to the global and polygon-shaped minima, we found some other local minima. Most of these have energies between the global minimum and the polygon and have various (or no obvious) symmetries, like the distorted rectangle of mentioned above. The one notable exception is the lucky cloverleaf at , whose energy exceeds that of polygon and has an obvious symmetry. It is depicted in Fig. 16LABEL:sub@fig:b14.

We now briefly return to the case . Since it is the last one where the polygon seems to be the global minimum, it deserves special attention: how can we be certain some other configuration does not give a lower energy? Since it seems conceivable that a rectangular lattice of half lumps could exist, we tried several methods of producing it: we split the rectangular minimizer in two, keeping just one half; spliced together all possible combinations of two lower charge solutions side by side and even spliced together two and two solitons in a crosslike configuration. They all relax either to the polygon or the rectangle. Of course, it is possible, that the rectangles reported above for the higher charge solitons might not be the true global minima, but we believe they are: any other rectangle would necessarily have a bigger difference between its long and short sides. Therefore, we did not attempt to find any other rectangles.

Finally, we summarize the charges, normalized energies per and shapes of our global minimizers in Table 1. We also conjecture an exponential decay law for of the global minimizer which is approaching the value when increases without limit. This behavior, along with a least squares fit of the form , is displayed in Fig.17.

B | Global minimum E/B | Polygon E/B | Shape |
---|---|---|---|

1 | 1.707 | 1.707 | pair of half lumps |

2 | 1.591 | 1.591 | square |

3 | 1.565 | 1.565 | regular polygon |

4 | 1.556 | 1.556 | regular polygon |

5 | 1.551 | 1.551 | regular polygon |

6 | 1.548 | 1.548 | regular polygon |

7 | 1.547 | 1.547 | rectangle |

8 | 1.544 | 1.546 | square |

9 | 1.543 | 1.545 | rectangle |

10 | 1.541 | 1.545 | rectangle |

11 | 1.540 | 1.544 | distorted rectangle |

12 | 1.537 | 1.544 | rectangle |

13 | 1.537 | 1.544 | distorted rectangle |

14 | 1.535 | 1.544 | rectangle |

18 | 1.532 | 1.543 | square |

32 | 1.526 | 1.543 | square |

## Vi Concluding remarks

We have used numerical minimization to show that multi-Skyrmions in the easy plane baby Skyrme model exist and are stable. For each charge above 4, we have found at least two minima, sometimes three or even four. Of these, we have most likely identified the global one for most charges in the range , , and . These minima exhibit various different types of symmetries and it is conceivable that a method exists for determining the shape of the minima without finding the minima themselves. However, this remains an open question. The prevalence of rectangles or slightly distorted rectangles in this work might provide a way to answer the question. We also present a conjecture that .

Like baby Skyrme models with some other potentials, we have identified an attractive and repulsive channel for two well-separated unit Skyrmions using a dipole approximation and confirmed its validity by numerical simulations.

The rectangular lattices formed by the half lumps are reminiscent of condensed matter systems, where (baby) Skyrmions have been experimentally observed (see, for example, Mühlbauer et al. (2009); Yu et al. (2010)). These are not, however, modeled by the baby Skyrme model, so it would be interesting to investigate whether the easy plane baby Skyrme model could be used as an effective model for some condensed matter systems.

###### Acknowledgements.

This work was supported by the UK Engineering and Physical Sciences Research Council. JJ also wants to thank Lisandro Dalcin for help with the tao4py library.### References

- E. Witten, Nucl. Phys. B223, 422 (1983a).
- E. Witten, Nucl. Phys. B223, 433 (1983b).
- S. L. Sondhi, A. Karlhede, S. A. Kivelson, and E. H. Rezayi, Phys. Rev. B 47, 16419 (1993).
- I. Hen and M. Karliner, Phys. Rev. E 77, 036612 (2008a), eprint 0711.1974.
- I. Hen and M. Karliner, Phys. Rev. D 77, 116002 (2008b), eprint 0802.2348.
- I. Hen and M. Karliner, Theor. Math. Phys. 160, 933 (2009).
- I. Hen and M. Karliner, Phys. Rev. D 77, 054009 (2008c), eprint 0711.2387.
- I. Hen and M. Karliner, Nonlinearity 21, 399 (2008d), eprint 0710.3939.
- C. Adam, P. Klimas, J. Sánchez-Guillén, and A. Wereszczyński, Phys. Rev. D 80, 105013 (2009), eprint 0909.2505.
- C. Adam, T. Romańczukiewicz, J. Sánchez-Guillén, and A. Wereszczyński, Phys. Rev. D 81, 085007 (2010), eprint 1002.0851.
- J. M. Speight, J. Phys. A 43, 405201 (2010), eprint 1006.3754.
- U. K. Rößler, A. N. Bogdanov, and C. Pfleiderer, Nature 442, 797 (2006), eprint arXiv:cond-mat/0603103 and arXiv:cond-mat/0603104.
- X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, Nature 465, 901 (2010).
- S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Science 323, 915 (2009), eprint 0902.1968.
- A. Neubauer, C. Pfleiderer, B. Binz, A. Rosch, R. Ritz, P. G. Niklowitz, and P. Böni, Phys. Rev. Lett. 102, 186602 (2009), eprint 0902.1933.
- Z. F. Ezawa and G. Tsitsishvili, Rep. Prog. Phys. 72, 086502 (2009).
- R. A. Leese, M. Peyrard, and W. J. Zakrzewski, Nonlinearity 3, 773 (1990), URL http://stacks.iop.org/0951-7715/3/773.
- S. Benson, L. C. McInnes, J. Moré, T. Munson, and J. Sarich, Tech. Rep. ANL/MCS-TM-242, Mathematics and Computer Science Division, Argonne National Laboratory (2007), http://www.mcs.anl.gov/tao.
- S. Balay, K. Buschelman, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith, and H. Zhang, PETSc Web page (2009), http://www.mcs.anl.gov/petsc.
- S. Balay, K. Buschelman, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith, and H. Zhang, Tech. Rep. ANL-95/11 - Revision 3.0.0, Argonne National Laboratory (2008).
- S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith, in Modern Software Tools in Scientific Computing, edited by E. Arge, A. M. Bruaset, and H. P. Langtangen (Birkhäuser Press, Boston, 1997), pp. 163–202.
- C. G. Broyden, J. Inst. Math. Appl. 6, 76 (1970).
- R. Fletcher, Comput. J. 13, 317 (1970).
- D. Goldfarb, Math. Comput. 24, 23 (1970).
- D. F. Shanno, Math. Comput. 24, 647 (1970).
- M. Peyrard, M. B. Piette, and W. J. Zakrzewski, Nonlinearity 5, 563 (1992).
- B. M. A. G. Piette, B. J. Schroers, and W. J. Zakrzewski, Z. Phys. C 65, 165 (1995), eprint arXiv:hep-th/9406160.
- Y. Choquet-Bruhat, C. DeWitt-Morette, and M. Dillard-Bleick, Analysis, Manifolds and Physics, Part I (North-Holland, Amsterdam, The Netherlands, 1982), p. 495.