NonlinearSchrodinger: Higher-Order Algorithms and Darboux Transformations for Nonlinear Schrödinger Equations

1) calculation of analytical solutions via a numerical implementation of the Darboux transformation
2) easy-to-use program package
3) large number of propagation algorithms implemented

1) applicable only in 1+1 dimensions, which is not mentioned in the abstract or in the introduction
2) very limited in scope (in terms of physical systems)
3) lack of systematic approach
4) state-of-the-art is not properly presented and references are not appropriate

The author presents a Julia package for solving a nonlinear Schroedinger equation with a cubic attractive nonlinearity. The package is also capable of calculating analytic solutions to NLSE using the Darboux transformation.

As noted by Referee 1, the author makes quite strong claims regarding the availability of open source software packages for solving of this type of NLSEs. Actually, there are quite a lot packages that implement different algorithms and approaches, published in various journals (such as Computer Physics Communications, Journal of Computational Physics, etc.), which are capable of solving not only a very narrow type of NLSEs considered in this manuscript, but a much broader class of equations and, hence, a much broader class of physical systems. In particular, there are packages are capable of solving NLSEs with higher order polynomial nonlinearities, nonlinearities of integral or integro-differential type (e.g., packages for solving the dipolar GPE), etc. Furthermore, the software packages are available for systems in 1D, 2D, 3D (plus time), parallelized using OpenMP, MPI, CUDA, etc. The author does not comment on this, although this would be highly beneficial for the audience. Also, the fact that the presented package only works for 1D systems is not mentioned in the abstract or in the introduction.

The author also does not comment on various algorithms used to solve NLSEs, and only focuses on the spectral algorithm, emphasizing the importance of higher-order integrators. While they may be relevant and important in some use-cases, for broad classes of physical systems second-order integrators are of sufficient quality (as noted by Referee 1), in particular when experiments are modeled. On the other hand, it would be quite useful to mention other algorithms in the manuscript (e.g., various Crank-Nicolson options) and motivate the choice made in writing the presented software package.

While the scope of the current manuscript is quite limited (in terms of physical systems), it also lacks systematic presentation of the properties of employed algorithm. For instance, a large number of integrators is implemented and their performance is presented, but the behavior of numerical errors is not systematically studied. There is no analytic expression for the behavior of errors, nor it is analytically derived how the errors should scale. From a practical point of view, if the author has determined which integrator is the best, why include all others in a package intended for a broad audience?

The main strength of the presented software package, in my opinion, is its capability to obtain analytical solutions via a numerical implementation of the Darboux transformation. Notwithstanding the fact that this is only possible for quite a narrow type of NLSEs, Eqs. (29) and (30), this represents the original contribution of the manuscript. Perhaps the paper could be refocused to this topic, with possible broadening of the applicability range, Eqs. (29) and (30).

Let me also note two technical issues:
1) the correct transliteration of Schr\"odinger is Schroedinger, not Schrodinger
2) $N$ is used with different meanings in different parts of the manuscript

1) Very easy to download and use
2) Very flexible in the choice of the algorithm for time stepping
3) unique (as far as I know) implementation of Darboux transformation to obtain analytical solutions

1) Limited in scope: only a (1+1)-dimensional Non Linear Schrodinger equation with anomalous dispersion can be simulated
2) Lack of benchmarks on CPU time and scalability

The author presents a Julia package to solve the (1+1)-dimensional Non-Linear Schrodinger Equation (NLSE) with anomalous dispersion (in the language of quantum gases, in which the undersigned referee is specialized, with "attractive interaction"), by using several different schemes to various orders of accuracy, based on the very well known technique of "operator splitting" and pseudo-spectral discretization using Fast Fourier Tranforms. It also allows to form classes of analytical solutions of the equation by using the Darboux transformation.

The paper is altogether clear and well written.

To my opinion, the manuscript, upon revisions, may meet the criteria for publication in SciPost Physics Codebases.

More in detail, let me start with a discussion of the first part of the manuscript, dedicated to an introduction to the NLSE and to the typical schemes for its numerical solution.

In the introduction, the author expresses some concern about the rarity of open source packages devoted to the study of the NLSE. Moreover, a strong opinion is expressed by the author on the validity of works employing second order split-operator methods, and whether such methods should be used or not.
Now, contrary to what the author states, on one side, there are already different packages devoted at least to the solution of the so-called Gross-Pitaevskii equation (GPE), which is basically a NLSE used for the description of Bose-Einstein condensates (see, e.g., https://hal.archives-ouvertes.fr/hal-02957115/file/PaperBEC2HPC.pdf and references therein). Moreover, most, if not all, the works in which the numerical solutions of the GPE are obtained, a second-order split-operator method is used, typically in combination with pseudo-spectral discretization based on FFTs, just like the author of this manuscript does. This method has been used now for more than two decades for the study of Bose-Einstein condensates by hundreds of researchers, including the undersigned.
Even if the author of this manuscript thinks the opposite, the use of such methods (with some care, of course) allows to obtain fast, accurate and robust codes for the solution of the NLSE, leading to high quality research, also for very complex simulations. The assertion that such methods should not be used for anything but simple simulations is by far over ambitious and should be avoided.

Let me now discuss some aspect of the part devoted to the benchmark of the codes. Here, the author shows a comparison of the numerical accuracy of the methods implemented in the package with respect to the conservation laws of the system, as a function of the number of points with which the x-axis is sampled (i.e., the number of "time"-steps, for a fixed length of the simulation). While one can appreciate the accuracy achieved by using higher order integrators for the same number of "time"-steps, I think that for a good benchmark it is fundamental to also include details about CPU-time and scalability. Moreover, the author states that it is necessary to use higher order integrators in order to be sure that numerical artifacts are not produced. This is simply not true: as the author also states later, being careful with the finess of the mesh is enough to avoid bad surprises, even without extremely accurate integrators and quadruple precision. Checking the stability of solutions after changing the time step is a very basic procedure, performed even by the most inexperienced programmers, so that I cannot see the necessity of making such a recommendation to people using lower order integrators.

I appreciate the fact that installing and using the package is really easy even for people with a very little experience in Julia, at least for the limited scope for which the package is designed. However, as stated formerly, let me underline again that the package can handle only the NLSE in (1+1) dimensions with anomalous dispersion, while adding other terms to the equation (necessary for handle a broad spectrum of physical systems modeled by the NLSE) seems not possible in current implementation.

The author finally presents another functionality of the package, i.e. the generation of analytical solutions of the NLSE using the Darboux transformation. This is, in my opinion, the main point of strength and of interest of the paper, since to my knowledge there are no other packages implementing such feature, but also because this part of the package can handle some more general cases. However, also in this case one must notice that the defocusing case ("repulsive" interaction) cannot be handled, as well as the presence of an eventual external potential. This again limits the scope of the package. Finally, a benchmark of the package based on CPU time and a discussion on scalability are, in my opinion, necessary.

In conclusion, the paper presents a Julia package through which the 1+1 dimensional NLSE with anomalous dispersion can be solved using a wide variety of approaches. The package is simple to install and use, and may be useful for researchers in the field of fiber optics, as well as for newcomers. Nonetheless, to my opinion, before being worth publication on SciPost Physics Codebases, the author must address the following points:
1) loosen certain claims on the validity of some well established and long used methods in the field;
2) provide benchmarks on CPU times;
3) discuss the scalability of the package to multicore CPUs;
4) discuss possible extensions that include, at least, the defocusing case, the presence of external potentials, higher dimensional equations (2+1 and 3+1). If this is not possible or not on schedule, it should be stated clearly at least in the abstract.

Upon considerations of these comments, I would recommend publication in SciPost Physics Codebases

As stated above:

1) loosen certain claims on the validity of some well established and long used methods in the field;
2) provide benchmarks on CPU times;
3) discuss the scalability of the package to multicore CPUs;
4) discuss possible extensions that include, at least, the defocusing case, the presence of external potentials, higher dimensional equations (2+1 and 3+1). If this is not possible or not schedule, it should be stated clearly at least in the abstract.