# High Order Rogue Waves on the Salerno Lattice

## Context
We're interested in understanding the focusing mechanisms for the so-called Salerno model
$$
i \frac{d u_n}{d t}=-h^{-2}\left(u_{n+1}+u_{n-1}\right)-\left|u_n\right|^2\left(2\mu u_n+(1-\mu)\left(u_{n+1}+u_{n-1}\right)\right)
$$
where $h$ is the lattice spacing the Salerno parameter $\mu$ interpolates from the so-called Ablowitz-Ladik (AL) model $(\mu=0)$ to the so-called discrete nonlinear Schrodinger (DNLS) equation $(\mu=1)$. The AL model is integrable which makes this a testbed for our favorite mathematical methods.

We found that in the continuous case, we could engineer Thomas-Fermi-like initial conditions and use these as proxies for the Talanov semi-circular initial condition; this initial condition was shown to generate infinite order rogue waves. The procedure for generating the TF state is described in
https://arxiv.org/abs/2406.06869
The parameters (nonlinear strength, trap strength, chemical potential) I used to generate the TF state for the Salerno model is not that critical. The main point is to fix a family of initial conditions and to track the dynamics of the model as the lattice stepsize $h$ is varied.

## Observed Dynamics
All plots are rescaled by the value of the ground state at the origin so that what we observe the magnification. For the Ablowitz-Ladik model I observe
For DNLS
For Salerno

We see there is an $h$ value where the dynamics of the Ablowitz-Ladik model abruptly change into a strongly focusing wave. I want to understand this and I propose to do so by computing the nonlinear spectrum of the Ablowitz-Ladik model over $h$. Then, maybe I can use this as a proxy to understand Salerno and somehow continue through the parameter $\mu$ from zero to one.

## Reproducing Chen and Pelinovsky, 2023
Luckily, rogue waves on the Ablowitz-Ladik lattice where studied extensively by Chen and Pelinovsky
https://arxiv.org/abs/2303.17052
Indeed, I am quite fortunate that they did work here before I did. Like every good Pelinovsky paper on an integrable model, he states *very* cleanly what the Lax pair is:
\begin{equation*}
U\left(u_n, \lambda\right)=\frac{1}{\sqrt{1+\left|u_n\right|^2}}\left(\begin{array}{cc}
\lambda & u_n \\
-\bar{u}_n & \lambda^{-1}
\end{array}\right)
\end{equation*}

\begin{equation*}
\begin{aligned}
V\left(u_n, \lambda\right):
V_{11}&=\frac{i}{2}\left(\lambda^2+\lambda^{-2}+u_n \bar{u}_{n-1}+\bar{u}_n u_{n-1}\right) \\
V_{12}&=i\lambda u_n-i\lambda^{-1} u_{n-1} \\
V_{21}&=-i\lambda \bar{u}_{n-1}+i\lambda^{-1} \bar{u}_n \\
V_{22}&=-\frac{i}{2}\left(\lambda^2+\lambda^{-2}+u_n \bar{u}_{n-1}+\bar{u}_n u_{n-1}\right)
\end{aligned}
\end{equation*}

Compatibility of the operators $V$ and $U$ can be shown to reproduce the $h=1$ Ablowitz-Ladik model written above. As written, computing the spectrum of $U\left(u_n, \lambda\right)$ amounts to solving an operator pencil problem. Originally, we implemented a Newton trace method with deflation (to remove multiple roots) that numerically solves this nonlinear root-finding problem. Although I had a lot of fun in operator pencil land and even had a working implementation, this method is too computationally inefficient for this nonlinear waves study. Chen and Pelinovsky of course found a much easier approach. Just bit of algebra and reshifting indices leads to
$$
\left\{\begin{array}{l}
\sqrt{1+|u_n|^2} P_{n+1}+\sqrt{1+|u_{n-1}|^2} P_{n-1}-\left(U_n-U_{n-1}\right) Q_n=z P_n \\
\sqrt{1+|u_n|^2} Q_{n+1}+\sqrt{1+|u_{n-1}|^2} Q_{n-1}+\left(\bar{u}_n-\bar{u}_{n-1}\right) P_n=z Q_n
\end{array}\right.
$$
where $z:=\lambda+\lambda^{-1}$. This is now an eigenvalue problem in $z$, and much easier to solve than an operator pencil problem!

Since I'm using periodic boundary conditions, the problem I'm left with is
\begin{equation*}
M\xi=z \xi
\end{equation*}
where
\begin{equation*}
M=\left[\begin{array}{l|l}
A & B \\
\hline C & A
\end{array}\right]
\end{equation*}
\begin{equation*}
a_n=\left(1+\left|u_{n-1}\right|^2\right),\qquad b_n=\left(u_{n-1}-u_n\right),\qquad c_n=\left(\bar{u}_n-\bar{u}_{n-1}\right)
\end{equation*}
\begin{equation*}
B=\operatorname{diag}(\vec{b}), \qquad C=\operatorname{diag}(\vec{c})
\end{equation*}
$$
A = \begin{pmatrix}
0 & 1 & 0 & 0 & \dots & a_{1} \\
a_2 & 0 & 1 & 0 & \dots & 0 \\
0 & a_3 & 0 &1 & \dots & 0 \\
\vdots & \ddots & \ddots & \ddots & & \vdots \\
1 & 0 & \dots & & a_N & 0
\end{pmatrix}
$$
Using the Floquet-Hill method to fill in the continuous spectrum turns out to be really easy; this 20 year old article shows how easy it is
https://www.sciencedirect.com/science/article/pii/S0021999106001665

We now perform a test to compute the nonlinear AL spectra of the dnoidal

$$u_{\rm dn}=\frac{\operatorname{sn}(\alpha ; k)}{\operatorname{cn}(\alpha ; k)} \operatorname{dn}(\alpha n ; k),$$
and cnoidal
$$u_{\rm cn}=\frac{k \operatorname{sn}(\alpha ; k)}{\operatorname{dn}(\alpha ; k)} \operatorname{cn}(\alpha n ; k)$$
waves, where $\alpha \in(0, K(k))$ and $k \in(0,1)$, where $k$ is the elliptic modulus.

Compare this with Chen and Pelinovsky, Figure 3

Compare this with Chen and Pelinovsky, Figure 5 top

Compare this with Chen and Pelinovsky, Figure 5 bottom


Looks like my implementation is aligned with theirs. Here's the whole MATLAB code so that you can run it yourself and see how satisfying it is to use Floquet-Hill to fill out the continuous spectrum

%% A Test to Compute the nonlinear Ablowitz-Ladik Spectrum

% spatial domain is one period of a discrete dnoidal wave
ss=0;
N=40+ss*40;
x=-(N-1)/2:(N+1)/2;
x(end)=[];

% wave parameters consistent with Chen and Pelinovsky Section 4.2
M=20;
k=.8;
K= ellipke(k);
alpha=K/M;

% here is the dnoidal wave (Equation 2.8)
[S,C,~] = ellipj(alpha,k);
[~,~,D] = ellipj(alpha*x,k);
u0=S.*D./C;

% here is the cnoidal wave (Equation 2.11)
% [S,~,D] = ellipj(alpha,k);
% [~,C,~] = ellipj(alpha*x,k);
% u0=S/D.*C;

% use the Floquet-Hill Method and loop over Floquet parameters
for theta=0:pi/M/8:pi

% Compute the Ablowitz Ladik Spectrum for one theta
[lambda,z]=getLaxSpectrum(u0,theta);

% visualize the spectrum (compare with Figure 3)
subplot(3,1,1)
plot(real(lambda),imag(lambda),'.k'), axis tight, hold on
xlabel('$\Re\{\lambda\}$'),ylabel('$\Im\{\lambda\}$')
title('Nonlinear Ablowitz-Ladik Spectrum')
xlim([.9 1.1]), ylim([-.08 .08])
subplot(3,1,2)
plot(real(lambda),imag(lambda),'.k'), axis tight, hold on
xlabel('$\Re\{\lambda\}$'),ylabel('$\Im\{\lambda\}$')
xlim([-1.1 -.9]), ylim([-.08 .08])
subplot(3,1,3)
plot(real(lambda),imag(lambda),'.k'), axis tight, hold on
xlabel('$\Re\{\lambda\}$'),ylabel('$\Im\{\lambda\}$')
axis equal
drawnow
end

function [lambda,z]=getLaxSpectrum(u0,theta)

N=length(u0);

% make the u0 dependent functions appearing in the block matrix 
for k=2:N
    a(k)=sqrt(1+abs(u0(k-1))^2);
end

for k=1:N
    d(k)=sqrt(1+abs(u0(k))^2);
end

% enforces periodicity
a(1)=a(end);

% note that these coefficients equal zero when k=1 (as they should)
for k=2:N
    b(k)=(u0(k-1)-u0(k));
    c(k)=(conj(u0(k))-conj(u0(k-1)));
end

% make the tridiagonal matrix A
A=spdiags(a(2:end).',-1,N,N)+spdiags(d.',-1,N,N).';

% enforce periodicity in A (this is the only Floquet-Hill modification)
A(end,1)=d(end).*exp(-1i*theta);
A(1,end)=a(1).*exp(1i*theta);

% construct the inhomogeneous diagonal matrices B and C
B=spdiags(b.',0,N,N);
C=spdiags(c.',0,N,N);

% the full matrix is
M=[A B;C A];

% compute the spectrum
[~,DD]=eigs(M,2*N);
z=diag(DD);

% infer lambda via the quadratic formula
lambdap=(z+sqrt(z.^2-4))/2;
lambdam=(z-sqrt(z.^2-4))/2;
lambda=([lambdam; lambdap]);

end

## Computing the Nonlinear Spectrum for the TF States
First of all, I need to verify that the nonlinear spectrum of the dynamics is invariant over time. Indeed, if I compute the nonlinear spectrum of $h\vec{u}(t)$ using the above MATLAB code from some $\vec{u}_0$ then I observe isospectrality of the spectral operator. I was lucky that the rescaling was simple and I found it. Great.

Here is what the spectra look like for AL ground states

DNLS ground states

and Salerno ($\mu=1/2$) ground states

I think there's some interesting mathematics that could probably be done here. I could probably resolve space so that $h$ values change more continuously. By the way, here is the entire MATLAB code that gets the ground states and computes the nonlinear spectrum of each one.

% spatial domain 
L=3*pi;
mu=0;

for N=3:2:51
x=linspace(-L,L,N)';

%lattice spacing
h=x(2)-x(1);

% build the matrix for the derivative
A=diag(ones(1,N-1),-1) + diag(ones(1,N-1),1);

% use periodic boundary conditions
A(1,N)=1;
A(N,1)=1;

% time domain to get the Ground State
T=8;dt=2^-14;
t=0:dt:T;

% T-F like Initial Condition Computed as the GS with potential V
omega=3;
V=(omega^2/2*x.^2);
tol=1e-12;
u=4.5*exp(-x.^2);
u0=getGroundState(x,u,A,V,mu,h,t,tol);
% compute the nonlinear spectrum using Floquet-Hill
    for theta=0:pi/80:pi
        % Compute the Spectrum for one Theta
        [lambda,z]=getLaxSpectrum(h*u0,theta);

        % visualize
        plot(real(lambda),imag(lambda),'.k'), axis tight, xlim([-2 2]), hold on
        xlabel('$\Re\{\lambda\}$'),ylabel('$\Im\{\lambda\}$')
        axis equal
        title(['Nonlinear AL Spectrum at h=' num2str(h)])
    end
    % draw it now that the continuous spectrum has been filled in
    drawnow
    
    % if you want some control, uncomment to pause
    % pause
    clf
end

function u=getGroundState(x,u0,A,V,mu,h,t,tol)

nt=length(t);
dt=t(2)-t(1);

norm=sqrt(trapz(x,abs(u0).^2));
MatExp=expm(dt*A*h^(-2));

% set the initial condition
u=u0;

%Time-Stepping Loop
for ktime=2:nt

    % Full-Step of DNLS 
    u=u-2*dt*mu*u.^3-dt/2*V.*u;

    % Full-Step of AL NL
    u=u-dt*(1-mu)*abs(u).^2.*A*u-dt/2*V.*u;

    % Full Step of Linear
    u=MatExp*u;

    % Renormalize
    u=norm*u/sqrt(trapz(x,abs(u).^2));

    if trapz(x,abs(u-u0).^2)<tol
        % we've converged to the steady state
        break
    end
    u0=u;
end

end

For completeness, here's the MATLAB code that computes the focusing Salerno dynamics

function [u,U]=SolveFocusing(u,mu,A,MatExp,x,t,nplt,plots)
nt=length(t);
dt=t(2)-t(1);
counter=1;
U(:,counter)=u;

%Time-Stepping Loop
for ktime=2:nt

    % Full Step of Linear
    u=MatExp*u;

    % Full-Step of DNLS NL
    theta=angle(u)+2*mu*dt*abs(u).^2;
    u=abs(u).*exp(1i*theta);

    % Full-Step of AL NL
    u=u+1i*(1-mu)*dt*abs(u).^2.*A*u;

    % visualize
    if mod(ktime,nplt)==0
        counter=counter+1;
        U(:,counter)=u;
        if plots
        plot(x,abs(u))
        xlim([-10 10]), ylim([0 5]), pause(1), hold off
        end
    end

end
end

## Further Studies
I noticed that can interpolate from a "Christmas tree pattern" to a KM breather. Interesting. The way I do this is by rescaling the dynamics through a single parameter $a$, I did away with the $h$ dependence, such a way the the continuum dynamics are acheived via the ansatz
$$
u_n(t)=a \mathfrak{u}\left(a n, a^2 t\right) e^{2 i t}
$$
where $\mathfrak{u}$ asymptotically approaches the solution the continuum NLS
$$
i \mathfrak{u}_T+\mathfrak{u}_{X X}+2|\mathfrak{u}|^2 \mathfrak{u}=0
$$
as $a\to0$.


I would like to make quantitative statements about what I am observing. It would be great if I could compute the spectrum of exact solutions and make fits of my computed dynamics to them, etc...

By the way, I found this interesting picture of the spectrum over specific values of $a$. I don't know how to interpret it, but it looks like a cnoidal spectrum colliding and merging into the main branch of the continuous spectrum.