Currently, I'm developing a numerical method for the nonpolynomial Gross Pitaevski equation (NPSE) [1] in one spatial dimension \begin{eqnarray} i \hbar \partial_t &= -\frac{\hbar^2 }{2 m} \psi_{zz} + V_{\text{ext}} \psi + \frac{2 \hbar \alpha_s \omega_{\perp}}{\sqrt{1 + 2 \alpha_s |\psi|^2 }} \left| \psi \right|^2 \psi +\frac{\hbar \omega_{\perp}}{2 } \left( \frac{ 1}{\sqrt{1 + 2 \alpha_s |\psi|^2 }} + \sqrt{1 + 2 \alpha_s |\psi|^2 } \right)\psi. \end{eqnarray} by using a spectral Strang splitting scheme in time. According to the lack of known reference solutions of the dynamic case I need to validate the expected convergence rates by producing a numerical reference solution at a certain time $T_{\text{sim}}$. The corresponding reference solution is constructed by hyper fine time steps towards $T_{\text{sim}}$ on a very dense grid such that the error in time and space of the method can be neglected. This reference solution is now compared to other outputs of the algorithm using much larger time steps. By plotting the the global convergence rate for the wave function $\psi$ it now turns out that it is compressed to an over all linear convergence rate which is atypical for Strang splitting methods. But the same code (with a little adaption) produces a regular second order convergence rate on the "normal" Gross Pitaevskii equation \begin{eqnarray} i \hbar \partial_t \psi =& -\frac{\hbar^2 }{2 m} \psi_{zz} +V_{\text{ext}} \psi +g \left| \psi \right|^2 \psi \end{eqnarray} where a analytical refernce solution is known.
I talked to several colleagues about this concern and it turned out they are having similar problems with this kind of saturation effects by using numerical computed reference solutions. I studied the literature on numerical analysis and found only a few notes where it was suggested to consider the ratio of halfed step sizes in the case the reference solution is unknown (https://www.csc.kth.se/utbildning/kth/kurser/DN2255/ndiff13/ConvRate.pdf). However, this turned out to work a little bit better in my case but I want to know if there is any recent publication or book covering this topic in more depth. It is still unsatisfying not to have some nice plots which describe the rate of convergence correctly.
(I also checked the convergence rate for the spatial error of my method which drops down to machine precession very fast as I'm increasing the number of grid points. Such that the linear convergence rate in time should not be able to be backtraced to the magnitude of the spatial error.)
[1] L. Salasnich, and A. Parola and L. Reatto, ”Effective wave equations for the dynamics of cigar- shaped and disk-shaped Bose condensates”, Phys. Rev. A 65 (4) 043614 (6pages), 2002
To answer your explicit question: No, as far as I know there are no recent publications that cover the relevant topic in the necessary depth.
There is a fair chance that either Zahari Zlatev at DMU in Denmark or Ole Østerby of Aarhus University in Denmark will be able to answer your question. They are both experts on the subject matter and the Zahari is an author of recent book on the stability of methods derived through Richardson extrapolation. Richardson extrapolation is the technique discussed by Oluf Rundborg.
I would like to make a few remarks after reading your question.
Please be very careful when you use a highly accurate numerical reference solution to drawn any conclusions. If you eliminate the discretization error then what you have left is rounding errors and they are much more difficult to analyze. The successful application of Richardson's techniques requires you to remain inside the so-called asymptotic range, where rounding errors are insignificant compared with with the discretization error.
For the sake of clarity consider the problem of computing a target value $T$ using a family of approximations $A_h$ that depend on a single parameter $h>0$ which we are essentially free to modify and which satisfies an asymptotic error expansion of the form $$ T - A_h = \alpha h^p + \beta h^q + O(h^r), \quad h \rightarrow 0_+ $$ where $\alpha$ and $\beta$ are constants that are independent of $h$ and $$ 0 < p < q < r.$$ In this situation it follows immediately that Richardson's fraction, namely $$ F_h = \frac{A_{2h} - A_{4h}}{A_h - A_{2h}}$$ satisfies $$ F_h \rightarrow 2^p, \quad h \rightarrow 0_+.$$ As a matter of fact, there is a constant $c$ such that $$ \frac{F_h - 2^p}{h^{q-p}} \rightarrow c, \quad h \rightarrow 0_+.$$ In principle this allows you to conduct numerical experiments that verify the existence of an asymptotic error expansion and determines the values of $p$ and $q$. Moreover, Richardson's error estimate, i.e., $$E_h^{\text{est}} = \frac{A_h - A_{2h}}{2^p-1}$$ satisfies $$T - A_h = E_h^{\text{est}} + O(h^q), \quad h \rightarrow 0.$$
There are several restrictions that must be respected here.
In your case there are (at least) two adjustable parameters namely the time step ($k$) and the spatial step size ($h$). If you are lucky then there is a more general asymptotic error expansion, say, $$ \psi(t,z) - A_{h,k}(t,z) = \alpha(t,z) h^2 + \beta(t,z) hk + \gamma(t,z) k^2 + \text{higher order terms}.$$ for suitable functions $\alpha$, $\beta$ and $\gamma$ defined throughout your domain.
In principle this can be investigated experimentally as outlined above by fixing one parameter and varying the other systematically. In particular, we want to discover if $\beta = 0$ so that there is no pesky first order contribution to the error. Be advised that there are is the possibility of localized problems near the boundary of your domain and any singularities can also obscure the picture.
Finally you may have to consider the spectral approximation that you are doing. I am not implying that you have made a mistake, I simply cannot say that is accurate enough as I have no real information about it.