text:
\section{Introduction}
The $S=1$ Heisenberg antiferromagnetic chain,
\begin{equation} H=J\sum_j\vec S_j\cdot \vec S_{j+1},\end{equation}
has been a subject of intense theoretical and experimental study since
Haldane's observation that it has an excitation gap above a singlet
ground state to a triplet excited state, quite unlike the S=1/2 case. Much has been learned
about the time independent properties using a combination of analytic and numerical methods,
in particular the density matrix renormalization group (DMRG).\cite{dmrg}  Our knowledge of dynamical
properties, while substantial, leaves more room for improvement.\cite{kuhner}
Recently, DMRG has been
significantly extended to allow direct calculation of time dependence (tDMRG)\cite{vidal-time,whitefeiguin,daley,feiguinwhite} and
time-dependent correlation
functions in particular. In this paper we will utilize tDMRG to calculate high-resolution spectral
functions for the $S=1$ chain for a broad range of momenta.  Some of the numerical techniques developed here were
briefly described in Ref. \onlinecite{pereira}.

The spectral function we shall focus on, in the isotropic case, is
\begin{equation}  S(k,\omega )=\sum_{j=-\infty}^\infty e^{-ikj}\int_{-\infty}^\infty dt \ e^{i\omega t}
\langle 0|S^a_j(t)S^a_{0}(0)|0\rangle \label{S}\end{equation}
where
$S(k,\omega) \equiv S^{aa}(k,\omega )$, with $a = x$, $y$, or $z$, and the subscripts
on $S$ indicate sites.  The tDMRG method
calculates the space-time dependent expectation values appearing in Eq. (\ref{S}) directly,
and then one performs  the Fourier transforms (FTs) in Eq. (\ref{S})
to obtain $S(k,\omega )$.  The crucial practical
issue for this approach, which we discuss in detail in this paper, is dealing with the finite range
of times available from a tDMRG simulation.  We find that two different methods for extrapolation in time to increase the
range of times used in the FTs  are both very useful for increasing the frequency resolution.

Much insight into the model can be obtained from an approximate mapping,
based on the large S limit,
onto the non-linear $\sigma$-model (NL$\sigma$M) with Lagrangian density:
\begin{equation} \mathcal{L}={1\over 2gv}[(\partial_t\vec \phi )^2-v^2\partial_x\vec \phi )^2],\end{equation}
where the lattice spin operators are represented as:
\begin{equation} \vec S_j\approx {1\over vg}\vec \phi \times \partial_t\vec \phi +(-1)^js\vec \phi ,\end{equation}
$\phi^2=1$, $v$ is the spin velocity and $g\approx 2/s$ is the coupling constant. $\vec \phi $
is the antiferromagnetic order parameter and $\vec l\equiv (1/vg)\vec \phi \times \partial_0\vec \phi$
is the uniform magnetization density. This low energy representation is only
valid for wave-vectors near $0$ and $\pi$. This field theory is known
to have a singlet ground state in one space dimension, with a gap, $\Delta$, to
a massive spin-triplet of excitations.  There are no bound states
or any other single particle excitations besides this triplet.  This
implies that a single particle ($\delta$-function) peak should appear in $S(k,\omega)$
near $k=\pi$ at energy $\Delta$ but near $k\approx 0$ the lowest excitations are a 2-particle continuum
starting at $2\Delta$.
Interaction effects mix $\vec \phi$ with $(\vec \phi )^3$ [but not $(\vec \phi )^2$
by symmetry] so that the theory
also predicts a 3-particle continuum at $k\approx \pi$ starting at $3\Delta$, a 4-particle
continuum at $k\approx 0$, et cetera.
Continuity between $k\approx 0$ and $k\approx \pi$ then implies a
qualitative sketch like Fig. \ref{fig:spec_wt} for the region of non-zero spectral weight.  In particular,
the single particle peak must emerge out of the 2-particle continuum at some critical
wave-vector $k_c$.  Integrability of the NL$\sigma$M allows for the
calculation of exact form factors and hence predictions for the
detailed shape of the 2-particle contribution to $S(k,\omega )$ near $k\approx 0$
and the 3-particle contribution near $k\approx \pi$.\cite{weston,horton}

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.6\hsize,scale=1.0]{dispschematic}
\caption{General features of the spectral function. The four and higher magnon continua are not shown.
}
\label{fig:spec_wt}
\end{center}
\end{figure}

Since the field theory
is based on $S\to \infty$ and is only valid for $k$ very close to $0$ and $\pi$
it is unclear how well any of these predictions should describe the S=1 case.
Various comparisons of NL$\sigma$M predictions with numerical results
have been performed before, but these have necessarily focused primarily on equal time
correlators and thermodynamic quantities.  Experiments
on quasi-1-dimensional antiferromagnets have clearly confirmed
the Haldane gap but the 2-particle nature of the small $k$ excitations
and the existence of a gap at $k=\pi$ from the single particle
excitation at $\Delta$ to the bottom of the 3-particle continuum at $3\Delta$
have not been confirmed and have led to some questioning of the validity
of these field theory predictions.

One of the purposes of this
paper is to compare our tDMRG calculations of $S(k,\omega )$  with the NL$\sigma$M predictions.
In addition, we will
determine $k_c$ and examine the behavior near this wave-vector.

In Sec. II we review the time dependent DMRG algorithm and discuss methods for
extrapolating the time-dependent correlation function to longer times before Fourier transforming.
In Sec. III we discuss properties of the single magnon
excitation.  In Sec. IV we discuss $S$ near $k\approx 0$ and compare with
the  2-magnon NL$\sigma$M predictions.  In Sec. V
we study $S$ near $k\approx \pi$ and compare with
the 1-magnon plus 3-magnons predictions.  In Sec. VI
we study intermediate $k$ in the vicinity of $k_c$
where the single magnon peak first emerges.
Sec. VII contains our conclusions.

\section{Time dependent DMRG and extrapolation in time}
The time dependent DMRG method\cite{whitefeiguin,daley,vidal-time} has been described elsewhere; here
we describe some of the practical issues related to obtaining dynamical spectra.  First,
we describe the specific tDMRG algorithms used; second, time step and truncation errors;
and third,  extrapolation in
time and windowing for time Fourier transforms.

The first step of a calculation is to use ordinary DMRG to find the ground state of a finite open system to high
accuracy, with  typical lengths being $L=200-400$.
To avoid $S=1/2$ end states, we put real $S=1/2$ spins on the first and last
sites. Let $i$ be one of the two center sites. After obtaining the ground state $\phi$,
we apply  the operator $S^+_i$ or $S^z_i$ to create a state which
is a mixture of excited states, $\psi(t=0)$.
We target and time evolve  $\psi$, and also, in order to minimize time step errors in
correlation functions, $\phi$.  As the time evolution occurs, we measure the correlation function
$S(j-i,t)$ when the tDMRG step is at site $j$, accumulating the function for a wide range of separations and times.\cite{whitefeiguin}

For the time evolution we use only the  Suzuki-Trotter  decomposition methods,
which are very efficient for
chain systems with only nearest neighbor interactions.  We use two different variations.
First, we use the original method of White
and Feiguin,\cite{whitefeiguin} which has second order Trotter errors.   An advantage of
this method is that every DMRG step has a time evolution bond operator applied, whereas Trotter
methods splitting the links into even and odd groups\cite{daley}
 apply a time evolution operator on
only half the steps of a sweep, resulting in roughly twice as much total truncation error per
unit time evolved.  To obtain a measurement for a specific time, we perform the measurements
in a half sweep without any time evolution. Thus,
including measurements, we evolve in the repeated 6 half-sweep pattern:
evolve left-to-right, evolve right-to-left,
left-to-right measurement half-sweep (without time evolution), evolve right-to-left, evolve left-to-right, right-to-left measurement half-sweep.  During each half-sweep time evolution step, we
evolve a time step $\tau$, so that measurements are available with a time step of $2 \tau$.
Typically we use $\tau=0.1$. (We set $J=1$.)

The second method we use is a fourth order Trotter method,\cite{forestruth,feiguinwhite}
which with $\tau=0.1$ virtually eliminates time-step errors, at the expense of more sweeps for a given time, and
consequently larger accumulated truncation error.
We show below that the {\it second order} Trotter errors from the first method with a time step of  $\tau=0.1$
primarily lead to modest frequency shifts, shifting the Haldane gap by about 1\%, for example.
(The rate of decay in time of the correlation functions seems
not to be strongly effected.)  For high accuracy studies it is more efficient to use the
fourth order method rather than simply reducing $\tau$.
The decomposition we use is\cite{forestruth}
\begin{eqnarray}
&e^{(A+B)\tau + {\cal O}(\tau^5)} =
e^{A \theta \frac{\tau}{2}}
e^{B \theta \tau}
e^{A (1-\theta) \frac{\tau}{2}}
e^{B (1-2\theta) \tau}\nonumber \\
&\times e^{A (1-\theta) \frac{\tau}{2}}
e^{B \theta \tau}
e^{A \theta \frac{\tau}{2}}
\end{eqnarray}
with $\theta=1/(2-2^{1/3}) \approx 1.35$. Here $A$ would represent, say, the odd bonds,
and $B$ the even, and each of the seven terms is applied in a half-sweep.  Adding to this
a measurement half-sweep, a total of
eight half-sweeps are needed to evolve by $\tau$, with
measurements available with a time step of $\tau$.

\newcommand{t_{\rm max}}{t_{\rm max}}

Generally we specify a desired truncation error for each step, and
vary the number of states kept $m$ to achieve this truncation error. However, we also
constrain $m$ to be no larger than a specified $m_{\rm limit}$ (typically 1000-2000),
and no smaller than
a minimum $m_{\rm min}$ (typically 100-150). The purpose of $m_{\rm min}$ is to
reduce the truncation error to near zero, at little cost, in steps where the state has very small
entanglement. (These small entanglement steps occur outside the ``light cone'' of the initial disturbance
at site $i$, and time 0.) We specify a $m_{\rm limit}$ to avoid memory limitations and to avoid a
handful of steps (near site $i$) taking a large fraction of the computer time.

The total accumulated truncation error, $\varepsilon_{\rm tot}$,
summed over all DMRG steps since the start
of the time evolution, is readily available and useful.  At each step, a small part of each wavefunction
is discarded; the truncation error is the magnitude squared of these small parts.   In
Fig. \ref{truncerrors} we show that
 $\varepsilon_{\rm tot}$ gives a rough estimate of the
typical errors to be expected in measurements at that time.
The plot compares $\varepsilon_{\rm tot}$ with the errors in
$S(x,t)$ for $x=j-i=0$ and $x=1$.  The errors in $S(x,t)$ were estimated
by running with two different truncation error parameters $\varepsilon$,
$\varepsilon=2\times10^{-10}$ and $\varepsilon=4\times10^{-10}$. We also
see that errors
of order $10^{-5}$ are computationally feasible for times $t<10$, but that the errors steadily
grow with time. In this case the errors grow roughly linearly with time because
a target truncation error per step was specified, with the upper limit $m_{\rm limit}$ playing
a small role.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.7\hsize,scale=1.0]{truncerrors}
\caption{Comparison of the errors in measurements of $S(x,t)$ and the total
accumulated truncation error $\varepsilon_{\rm tot}$. Two runs were made using the
4th order Trotter decomposition, both with $\tau=0.1$, with $\varepsilon=2\times10^{-10}$ and $\varepsilon=4\times10^{-10}$.  The values of $\Delta S$ were the difference in the results
between these
two runs. The value for $\varepsilon_{\rm tot}$ was for the larger $\varepsilon$ run.
}
\label{truncerrors}
\end{center}
\end{figure}

To clarify further how this works in Fig. \ref{maxmplot} we compare
two runs, one a second order Trotter method, the other
the more accurate run of Fig.  \ref{truncerrors}. Both runs took comparable amounts of
computer time. We define $m_{\rm max}$  to be the
maximum number of states kept over all the steps in a half-sweep. Clearly
$m_{\rm max} \le m_{\rm limit}$.  In each half-sweep, the largest values of $m$ were
for steps near the center of the
 system,  where the spin operator was applied. After a moderate time, for these steps
 near the center $m$ was limited by $m_{\rm limit}$.  The relatively small number of
 these steps made the effect of $m_{\rm limit}$  on $\varepsilon_{\rm tot}$ small.
 The calculation time for a step is proportional to $m^3$, so limiting $m$  is important
 for efficiency.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.8\hsize,scale=1.0]{fig3}
\caption{Maximum number of states kept $m_{\rm max}$
versus time (upper-left curves), and total discarded weight
$\varepsilon_{\rm tot}$ versus time (lower right curves)
 for the systems of Fig. \ref{Nwtau}.
The simulation adjusted $m$ at each step  to try to achieve a total discarded weight of $\varepsilon$ for that step, subject to a maximum $m$ of either 1300 or 1800 (flat portions).
The figure
shows that in this case $\varepsilon_{\rm tot}$ depends more strongly on $\varepsilon$ than on the maximum
allowed $m$.  }
\label{maxmplot}
\end{center}
\end{figure}

The tDMRG yields directly the space-time dependent correlation function directly. The signal initiated by
the application of the spin operator in the center of the system spreads out with time.  Because the
system has a finite correlation length, the correlation function is only non negligible within a range
of about $|v t|$ of the center, where $v$ is the maximum spin velocity.
We always keep
the maximum time $t_{\rm max}$ of the simulation small enough so that the signal has not reached the edges at
the end of the simulation.  Hence there are essentially no spatial finite size effects. One can spatially Fourier
transform (FT) $x\to k$  for any $k$; the available values of $k$ are continuous, not discrete.
On the other hand, the feasible $t_{\rm max}$ is strictly limited by the available computer time,
and the correlation functions decay
slowly or not at all in time.
Truncation of the signal following by an FT $t \to \omega$ would result in severe ``ringing''. The
standard approach is to multiply the signal by a windowing function,
which typically resembles a Gaussian centered at $t=0$ but which
vanishes exactly at $\pm t_{\rm max}$.  A drawback of this approach is that most of the
data gets ``thrown away'', and the
frequency broadening of the spectrum is large.

To avoid this over-broadening, we have developed an alternative approach based
on linear prediction.\cite{linearprediction}
Before the time-frequency FT, we extrapolate the time signal to long times using
linear prediction. We then
apply a broad window which does not throw away a significant amount of the original data.
Linear prediction
extrapolates a discrete equally spaced time series $\{y_i\}$ as
\begin{equation} y_i = \sum_{j=1}^n d_j y_{i-j}\end{equation}
The coefficients $d_j$ are determined by the known data points $\{y_i\}$ by requiring
that their prediction for
each point $y_i$, based on $y_{i-n} \ldots y_{i-1}$,  vary as little as possible
from the actual value $y_i$, using
a least-squares criterion. One finds that the $d_j$ are determined from
correlation functions $\langle y_i y_{i+j}\rangle$,
where the average is over $i$, and the principle computational work in determining $d_j$
is the inverse of a $n\times n$ matrix. In our work we have used $n=20$,
so that the numerical work involved in the extrapolation is
negligible.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.8\hsize,scale=1.0]{fig1}
\caption{DMRG results for the real part of the onsite spin-spin correlation function $S^{zz}(x=0,t)$ for
the more accurate run of Fig. \ref{truncerrors}.
The red curve represents DMRG data out to $t=22$; the black curve is the DMRG data for $t \le 10$ and a linear
extrapolation for $t>10$.}
\label{reSpm}
\end{center}
\end{figure}

In Fig. \ref{reSpm} we show the effectiveness of this extrapolation. For this run which had $t_{\rm max}=22$,
the data for $t \le 10$ was extrapolated to longer times and compared with the DMRG results. The
error in the extrapolation starts out small and grows reasonable slowly. Provided one does not rely
on too long an extrapolation to try to achieve higher frequency resolution than the data supports,
this method performs much better than ordinary windowing methods.

In Fig. \ref{Nwtau} we compare the density of states, defined as $N(\omega) = S^{zz}(x=0,\omega)$,
for the two systems of Fig. \ref{maxmplot}. After the linear extrapolation, the results were
multiplied by a Gaussian window $\exp[-t^2/(20 t_{\rm max}^2)]$ before Fourier transforming. On this scale
the effects of the Trotter errors in the 2nd order data are almost not visible.  The results are considerably
sharper than with the alternative method which did not use extrapolation. The sharp peaks are
broadened square root singularities from the top and the bottom of the single-magnon dispersion.
Above the top of the band, two- and three-magnon contributions to the spectrum are visible in a small tail.
Using a larger window width
would increase the resolution at the cost of increasing the likelihood of artifacts from the extrapolation.  We
will consistently use this window width. Assuming the extrapolation is accurate, this means that our spectra
should look like the exact spectra, but broadened by convolving with a gaussian
\begin{equation}
\exp(-\frac{1}{2} \omega^2 / W^2)
\end{equation}
where the frequency resolution is $W = (t_{\rm max} \sqrt{10})^{-1}$.

Another approach is to fit the moderate time data to
the correct asymptotic long-time form, in this case stemming from the square root singularities,
to extrapolate to long times very accurately.  The linear extrapolation asymptotically gives exponential
decays in time, an incorrect assumption in this case,
so for very long times the fitting method can be more accurate.  A disadvantage of the fitting
method is that it assumes one has some understanding of the results analytically.  Another
disadvantage is that the fitting process can take much more computer time than the linear prediction method,
although still much less than tDMRG simulation itself.
In cases where
one does not know what sort of spectra to expect, one can first fit with the linear prediction method,
and then guess an asymptotic form for fitting.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.7\hsize,scale=1.0]{Nwtau}
\caption{Results for the density of states  $N(\omega)$ for the two systems of
Fig. \ref{maxmplot}, with linear extrapolation utilizing only the data out to $t_{\rm max}=20$. The two curves
are almost identical.
The third curve used the 2nd order data but did not
utilize extrapolation; instead, the data out to $t_{\rm max}$ was multiplied by a simple finite window based
on the cosine function. }
\label{Nwtau}
\end{center}
\end{figure}

For the case of $N(\omega)$, we have leading singularities of the form
\begin{equation} a \theta(\omega-\Delta) (\omega-\Delta)^{-1/2} + b \theta(\Omega-\omega) (\Omega-\omega)^{-1/2}
\label{sing}\end{equation}
where $\Omega$ is the maximum in the single magnon dispersion relation near $\pi/2$ of about 2.725.
This Fourier transforms to long time tails with leading terms of the form
\begin{equation} A \exp(-i \Delta t) t^{-1/2} + B \exp(-i \Omega t) t^{-1/2}.
\label{singt}\end{equation}
A more convenient form for fitting comes from the Fourier transform
integral identities
\begin{eqnarray} \int_{-\infty}^\infty d\omega e^{-i \omega t} &\theta(\omega-b) e^{-a (\omega-b)}
 (\omega-b)^g \nonumber\\=
&\Gamma(1+g)e^{-i b t} (a + i t)^{-1-g} \label{integral}\end{eqnarray}
and
\begin{eqnarray} \int_{-\infty}^\infty d\omega e^{-i \omega t} &\theta(b-\omega) e^{-a (b-\omega)}
 (b-\omega)^g\nonumber \\=
&\Gamma(1+g)e^{-i b t} (a - i t)^{-1-g}. \label{integraltwo}\end{eqnarray}
These identities are useful because the frequency expressions have only a single one-sided
singularity,  and convenient well-behaved Fourier
transforms.
Eq. (\ref{integral}) is useful for describing a singularity on a lower edge, while
Eq. (\ref{integraltwo})
describes an upper edge. Thus to fit to the tDMRG data for the Fourier transform of
$N(\omega)$, we use the asymptotic form
\begin{equation} A \exp(-i \Delta t) (a+it)^{-1/2} + B \exp(-i \Omega t) (b-it)^{-1/2}.
\label{singtt}\end{equation}
Fitting the accurate 4th order data over the time range $10-20$ with this form, we
find that the
fit within this range matches the data very accurately, with a typical absolute
deviation of about $2\times 10^{-4}$, or a relative error of about $10^{-3}$.

Using the fitting parameters and asymptotic form, one extends the data to large times.
We did not use the identities Eqs. (\ref{integral}) and (\ref{integraltwo}) to help perform
the Fourier transform (although this might be convenient);
we simply used a fast Fourier transform over a very large range of times (e.g. $-10000 \le t \le 10000$).
Over the fitting region in time a smooth transition is made from the data for
small times to the fit for large times. The results
of this procedure for $N(\omega)$ are shown in Fig. \ref{Nwfit}. The fits allowed
both $\Delta$ and $\Omega$ to vary; the result for the gap for the time range $10-20$ was
$\Delta = 0.4104327$, accurate to 4 digits (see next section).  The resulting spectra, based on fitting
over different time ranges, appears to be accurate to the line width in the figure.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.9\hsize,scale=1.0]{Nwfit}
\caption{Results for the density of states  $N(\omega)$ from fitting to the asymptotic
form Eq. (\ref{singtt}) over the time range indicated in the caption, and then extending
to very large times with the fit before Fourier transforming.  The two curves overlap,
indicating that the error in the determination of $N(\omega)$ is less than the line width. }
\label{Nwfit}
\end{center}
\end{figure}

\section{Single magnon results}

Fourier transforming the DMRG $S(x,t)$ data $x \to k$, we obtain $S(k,t)$ for any $k$.
The time Fourier transform
gives spectra with (potentially) single magnon and multimagnon contributions. As we discuss
below, for $k_c  < k < \pi$, with $k_c \approx 0.23\pi-0.24 \pi$, a single magnon delta function peak
is present, plus multimagnon continua at higher frequencies. In this section we focus on
this well-defined single-magnon mode.

First consider the band minimum, at $k=\pi$, where the excitation energy is the Haldane gap.
The most accurate method to determine the Haldane gap is still ground state DMRG, where
we have used the same method as in Ref. \onlinecite{huse}, but with up to $m=500$ and $L=400$, to
determine the Haldane gap to very high accuracy, $\Delta = 0.41047925(4)$.  (The
end coupling used to make the lowest excitation have $k=\pm\pi$ was $J_{\rm end}=0.50865$,
compared to $0.5088$ in Ref. \onlinecite{huse}.)

In Fig. \ref{Spiclose}, we show results for $S(\pi,\omega)$ near $\omega=\Delta$.
Linear extrapolation plus Fourier transforming as described above
give rather narrow Gaussian-shaped peaks.
The 2nd order peak is narrower because of a larger $t_{\rm max}$, and it is shifted from
the exact result because of Trotter error.  Because two and three magnon contributions are very
weak and separated in frequency from $\Delta$, a least squares fit to a pure exponential
is almost identical to the maximum of the broadened peaks. In fact, neglecting any Trotter
shifts, either the maximum or the fit
frequency provide  much more accurate determinations of the exact magnon energy than  the
peak widths would indicate.  For a set of $k$'s spaced $0.01 \pi$ apart we have fit the time data either
to a  pure complex exponential (for $k>1$) and for
$k_c<k<1$ to a complex exponential
 plus an asymptotic form describing the near-by two-magnon edge. These latter more
 complicated fits are discussed in Section VI.
The  frequency
of the exponential term
determines the dispersion
$\varepsilon(k)$ for $k>k_c$. With this fitting approach, a larger $t_{\rm max}$ is not very important compared to
the Trotter error, so we utilize the 4th order data. The corresponding result for $\varepsilon(\pi)$
is 0.41050; the error is only a few times $10^{-5}$. For smaller $k$ the error is expected to
be larger because the multimagnon continuum is  larger  relative to the single magnon peak,
and the continuum is closer in frequency,  but the errors for $k>k_c$ are probably no bigger
than $10^{-3}$.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.8\hsize,scale=1.0]{Spiclose}
\caption{Results for   $S(\pi,\omega)$ for the two systems of
Fig. \ref{maxmplot}. The Gaussian shaped curves come from linear extrapolation and Fourier transforming.  The vertical lines represent delta functions coming from a
least squares fit of $A e^{i \omega t}$ to $S(\pi,t)$. The 4th order fit frequency and the ``exact''
ground state DMRG result are indistinguishable in this plot.
}
\label{Spiclose}
\end{center}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.8\hsize,scale=1.0]{dispersion}
\caption{Single magnon dispersion. The curve labeled ``DMRG''  comes from the fit of a pure
exponential to the DMRG $S(k,t)$ data. The curve labeled ``Fit'' is the analytic expression
Eq. (\ref{fitdispersion}).
The final curve is the two magnon band minimum at $2 \varepsilon(\pi-k/2)$.
The inset shows the region near $k_c$, where the magnon line enters the two magnon continuum.
The first two curves are not very meaningful well below $k_c$, since the there is no single
magnon delta function. }
\label{dispersion}
\end{center}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.6\hsize,scale=1.0]{disp_near_pi}
\caption{Comparison of the single magnon dispersion from DMRG and the relativistic approximation,
Eq. \ref{Ek}. }
\label{disp_near_pi}
\end{center}
\end{figure}

In Fig. \ref{dispersion} we show $\varepsilon(k)$ from the 4th order run, along with an
analytic fit motivated by the NL$\sigma$ model.
The massive triplet excitations of the NL$\sigma$M have the
relativistic dispersion relation:
\begin{equation} \epsilon_0(\tilde k)=\sqrt{\Delta^2+v^2\tilde k^2}.\label{Ek}\end{equation}
Here the momentum $\tilde k$ of the NL$\sigma$M is $k-\pi$ for the spin chain.
In the NL$\sigma$M, $\tilde k$ can take any real value.  Of course,
in the spin chain, crystal momenta lie in the Brillouin zone, $|k|<\pi$.
This discrepancy limits the validity of the NL$\sigma$M,
especially when we consider multi-particle excitations.
$\epsilon_0(\tilde k)$ is only an approximation to the exact single magnon
dispersion relation, $\varepsilon(\tilde k)$.
 We expect a perfectly stable single
magnon excitation to exist for $\tilde k<\pi -k_c$ with this dispersion relation.
Strictly speaking $\varepsilon(\tilde k)$ is not defined for $\tilde k>\pi -k_c$.
The ``Fit'' curve shown in Fig. \ref{dispersion} is based on the
expression
\begin{equation} \epsilon (k) \approx
 \Delta\sqrt{1+\sum_{n=1}^5a_n\{ 1-\cos [n(\pi -k)]\}}.\label{fitdispersion}\end{equation}
with
the parameters $a_n$ given in Table 1.  The gap, at $k=\pi$,  from the data used in the fit
is $\Delta =0.410504$; one could also use the more accurate value $\Delta =0.41047925$.
$\epsilon (k)$ goes through a maximum of $2.72551$ at $k\approx .476\pi$
and has the value $1.96$ at $k\approx 0.23 \pi$ near where the single magnon excitation becomes
unstable. There is an inflection point ($d^2\epsilon /dk^2=0$) at $k_{\rm in}\approx 0.868 \pi$.
As shown in Fig. \ref{disp_near_pi}
$\epsilon (\tilde k)$, agrees quite well with the Lorentz invariant
approximation, $\epsilon_0(\tilde k )$, for $\tilde k <0.1\pi$ ,
and reasonably well for $\tilde k <0.2\pi$,
with $v\approx 2.472$.
Also shown in Fig. \ref{dispersion} is the two magnon band minimum, which for
$|k| < 2(\pi-k_{\rm in}) \approx 0.265 \pi$ is given
by $2 \varepsilon(\pi-k/2)$. Near
$\tilde k\approx 0$, this is approximately
\begin{equation} \epsilon (\tilde k)\to \Delta + {(v\tilde k)^2\over 2\Delta }.\end{equation}
The intersection of the two magnon minimum and the single magnon dispersion line
determines $k_c$, as show in the inset of Fig. \ref{dispersion}.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.6\hsize,scale=1.0]{magnon_minmax}
\caption{Single magnon line (black solid line), two magnon band (horizontal line fill, red),
and three magnon (dotted fill, green) bands. }
\label{magnon_minmax}
\end{center}
\end{figure}

The single magnon dispersion can be used to construct the multimagnon band minima and maxima,
assuming that the magnon-magnon interactions are negligible, and adding the single magnon energies.
The resulting bands are shown in
Fig. \ref{magnon_minmax}. In this construction, it was assumed that the single magnon line stops abruptly
at $0.23\pi \approx k_c$. This resulted in slope discontinuities visible in the two magnon band. However,
at $k_c$, as we discuss below, the sharp "single magnon peak"
can be described as either a single magnon or two magnon feature.  This ambiguity effectively blurs the
distinction between the two and three magnon bands near the slope discontinuities. If one changed the construction
to include the broadened ``single magnon peak" to values of $k$ below $k_c$, the two magnon bandwidth would
be broadened near the slope discontinuities. This suggests that, for example, at $k=\pi$ one might expect to
see the two magnon band maximum as a visible feature of the spectrum, but one might not see anything
for the two magnon minimum. (The spectrum at $k=\pi$ is shown in Section V.)

\begin{table}[b]
\caption{Coefficients in the fit of the single magnon dispersion relation, Eq. (\ref{fitdispersion}). \label{tableI}}
\begin{tabular}{cccccccc}
\hline\hline
$n$&
$a_n$\\
1 & 1.96615 \\
2 & 21.20162 \\
3 & -1.61279 \\
4 & -0.04766\\
5 & 0.02407\\
\hline
\hline\hline
\end{tabular}
\end{table}

The amplitude of the single magnon peak is approximated well by\cite{weston}
\begin{equation} S(q,t=0) \approx \frac{v Z }{\sqrt{v^2 k^2 + \Delta^2}}
\end{equation}
with $Z=1.26$. A related  quantity is the fraction of the spectral weight in the multimagnon
continua, given by
\begin{equation} f = \frac{S(q,t=0)-A}{S(q,t=0)} \end{equation}
where $A$ is the amplitude of the single magnon peak. This quantity is shown in Fig. \ref{fratio}.
It is interesting that $f$ shows non monotonic behavior with $k$, with  minima near
$k=\pi/2$ and $k=\pi$. This non monotonic behavior is roughly correlated with the gap between
the single magnon line and the lowest multi-magnon band minimum.  Naively, one might
expect that a nearby continuum has an easier time than a faraway one does in
taking spectral weight from the single magnon line.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.6\hsize,scale=1.0]{ff}
\caption{Fraction of the spectral weight in the multimagnon continua for $k>k_c$,
shown by the solid black line.  The gap between the single magnon line and the
lowest multimagnon band is shown by the dashed red line.}
\label{fratio}
\end{center}
\end{figure}

\section{$k\approx 0$}
We now discuss the properties of the spectral function for small $k$. That
$S\propto k^2$ as $k \to 0$
can easily be proven to be exactly true using the fact that
$\sum_jS^z_j|0\rangle =0$ (singlet ground state) and Taylor expanding
$S(k,\omega )$ in Eq. (\ref{S}) to second order in $k$.
On general symmetry grounds, based on the NL$\sigma$M, we expect that
near $k\approx 0$, S will contain only multi-particle continua
corresponding to  even numbers of bosons.  At a finite small $k$,
the lowest energy 2-magnon state
with total momentum $k$ is one in which each magnon has momentum $\pi+k/2$.
Therefore, since the dispersion is even about $\pi$,
 the bottom of the 2-magnon continuum should be
exactly at $2\epsilon (\pi-k/2)$. This remains true up to
$k=2(\pi -k_{\rm in})\approx .265\pi$ (recall $k_{\rm in}\approx 0.868\pi$ is the inflection point).
For a range of energies, at small enough $k$, the {\it only}
possible excitations have 2 magnons.  This is true
up to $\omega = 4\epsilon (k/4)$, for $k<(\pi -k_{\rm in})$.

A simplified, ``mean field'' version of the NL$\sigma$M is
a free massive boson model with Lagrangian:
\begin{equation} \mathcal{L}={1\over 2v}[(\partial_t\vec \phi )^2-v^2(\partial_x\vec \phi )^2-\Delta^2(\vec \phi )^2],
\label{Lfree}\end{equation}
and no constraint on $\vec \phi$.
Expanding $\vec \phi$ in boson creation and annihilation operators one finds the
free boson result
\begin{equation} S_0(k,\omega )\approx
{k^2\sqrt{\omega^2-(kv)^2-4\Delta^2}\over v[\omega^2-(kv)^2]^{3/2}}\theta (\omega^2-(vk)^2-4\Delta^2).
\label{S0}\end{equation}
As expected from general principles, this vanishes quadratically as $k\to 0$,
and also vanishes below the 2-magnon threshold, $\omega_{\rm th} = 2\sqrt{\Delta^2+(vk/2)^2}$.
The exact 2-magnon expression for $S(k,\omega )$, in the NL$\sigma$M is
known exactly and can be written:
\begin{equation} S_{0\sigma}(k,\omega )=S_0(k,\omega ){\pi^4\over 64}{1+(\theta /\pi )^2\over 1+(\theta /2\pi )^2}
\left({\tanh \theta /2\over \theta /2}\right)^2,\label{Ssig}\end{equation}
where the rapidity, $\theta$, is defined by:
\begin{equation} \theta = 2\cosh^{-1}[(\omega^2-v^2k^2)/(4\Delta^2)].\end{equation}

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.7\hsize,scale=1.0]{Stentht}
\caption{Comparison of  $|S(k,t)|$ for $k=\pi/10$ for two different runs, with 2nd order and
4th order Trotter decomposition runs of Fig. 2, and the asymptotic form expected analytically.}
\label{Stentht}
\end{center}
\end{figure}

Just above $\omega_{\rm th}$, both
Eqs. (\ref{S0}) and (\ref{Ssig}) rise as $(\omega-\omega_{\rm th})^{1/2}$, which would lead to
the asymptotic time behavior
\begin{equation} S(k,t) \sim e^{-i \omega_{\rm th} t} t^{-3/2}.\end{equation}
In Fig. \ref{Stentht} we compare $|S(k,t)|$ for a typical small value of $k$
with  $A t^{-3/2}$, finding good agreement for large $t$
with the empirical parameter $A=0.6$.
Using the linear prediction method, we Fourier transformed the $S(k,t)$ DMRG results for both
runs, and in  Fig. \ref{ktenthcomp}, they are compared to Eqs. (\ref{S0}) and (\ref{Ssig}).
They look qualitatively similar. In particular,
the threshold singularity appears the same and
they have peaks at similar frequencies.
However, the peak is about a factor of 2 larger in the DMRG
data than in the NL$\sigma$M.
 Furthermore, the field
theory results drop off much more slowly at large $\omega$. This
latter feature is to be expected since magnons of arbitrarily
high momentum are included in the field theory while
there is a cut off at the Brillouin zone boundary in reality.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.7\hsize,scale=1.0]{ktenthcomp}
\caption{$S(k,\omega)$ for $k=\pi/10$ for two different runs, with different accuracies and different
total times, both Fourier transformed using linear extrapolation.  For comparison, two analytic
results, based on Eqs. (\ref{S0}) and (\ref{Ssig}), are shown.}
\label{ktenthcomp}
\end{center}
\end{figure}

A simple way to improve the high frequency behavior
in the free boson approximation is to replace the relativistic
model of Eq. (\ref{Lfree}), by a Hamiltonian:
\begin{equation} H={1\over 2}\sum_k[\vec \Pi_k\cdot \vec \Pi_{-k}+
\epsilon^2(k)\vec \phi_k\cdot \vec \phi_{-k}],\end{equation}
where $\Pi^a_k$ is canonically conjugate to $\phi^a_k$
and $\epsilon (k)$ is the numerically determined single
magnon dispersion relation. The appropriate form
of the small $k$ spin operators:
\begin{equation} \vec S_k\approx \sum_{k'}\vec \phi_{k'}\times \vec \Pi_{k-k'},\end{equation}
is determined by the requirement that the spin commutation relations
are obeyed and that $\vec S_0$ commute with the Hamiltonian. This
changes $S_0(k,\omega )$ to:
\begin{equation} S_0\to {[\epsilon (k')-\epsilon (k-k')]^2\over
2\epsilon (k')\epsilon (k-k')|\epsilon '(k')-\epsilon '(k-k')|},\label{S0re}\end{equation}
where $\epsilon (k)$ is the exact (numerically determined)
dispersion relation and $\epsilon '(k)$ denotes its derivative. $k'$,
the momentum of one of the 2 magnons,
is determined from $\omega$ and $k$ by energy-momentum conservation:
\begin{equation} \epsilon (k')+\epsilon (k-k')=\omega .\label{k'}\end{equation}
For low enough energy, there is only one pair of solutions to Eq. (\ref{k'}),
and one element of the pair should be chosen in evaluating Eq. (\ref{S0re}).
The resulting improvement of the free boson result is shown in Fig. \ref{ktenthcompw}
As expected, there is little change near the peak and threshold,
but the high energy tail is cut off. (Here we restricted
each boson to have $|\tilde k|<\pi /2$.)

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.7\hsize,scale=1.0]{ktenthcompw}
\caption{$S(k,\omega)$ for $k=\pi/10$, comparing the near-exact DMRG results with the modified
free boson result.}
\label{ktenthcompw}
\end{center}
\end{figure}

\section{$k\approx \pi$}
In addition to the single magnon mode, which has most of the spectral
weight near $k=\pi$ there is also a contribution from
$3$, $5$, $\ldots$ magnons. Of these, the largest is expected
to be the $3$-magnon contribution. The lower threshold
for this, corresponding to each magnon having momentum $\tilde k/3$
is at $3\epsilon (\tilde k/3)\approx 3\sqrt{\Delta^2+(v\tilde k/3)^2}$.
At small $\tilde k$ the exact lower threshold approaches $3\Delta + (v\tilde k)^2/(6\Delta )$.
Note that the presence of a multi-magnon contribution is a consequence
of inter-magnon interactions; it vanishes for the non-interacting model
of Eq. (\ref{Lfree}).
The exact 3-magnon form factor is known for the NL$\sigma$M and the resulting
3-magnon contribution to $S$ can be expressed
in terms of an elementary integral. The result is compared to our
DMRG results in Fig. \ref{kpi3mag}.  Again there is a qualitatively
similarity, with a peak at a similar energy, but now the NL$\sigma$M
peak is about $3$ times too low and there is far too much
spectral weight at high energies. The total  spectral
weight in the 3-magnon peak compared to that
in the single magnon is found from DMRG to be 2.7\%.
The lower edge of the multimagnon band is given by the three magnon edge, c.f. Fig. \ref{magnon_minmax}.
The band has a sharp dropoff at the two magnon band maximum.  The three magnon band above that
is rather small. As discussed earlier, one does not expect a sharp feature for the two magnon minimum
(nominally near $\omega/\Delta \approx 9$), and none is visible.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.7\hsize,scale=1.0]{kpi3mag}
\caption{$S(k,\omega)$ for $k=\pi$ in the multimagnon frequency regime
for two different runs. Also shown are the results from the  NL$\sigma$M model.
The vertical line near $\omega/\Delta \approx 13$ is two times the maximum of the single magnon
dispersion, roughly locating the top of the two magnon band.
The tiny bump near $\omega/\Delta=2$ is an artifact of the FT of the much larger
single magnon peak at
$\omega/\Delta=1$, while we believe the small tail above $\omega/\Delta \approx 13$ is real.}
\label{kpi3mag}
\end{center}
\end{figure}

\section{$k\approx k_c$}
A remarkable feature of $S(k,\omega )$ which is completely
missed by the NL$\sigma$M approach is the merging of the
single particle peak into the 2-particle continuum at $k=k_c\approx 0.23\pi-0.24\pi$.
Results for $S(k,\omega)$ from
the linear prediction method near
$k_c$ are shown in Fig. \ref{Snearkc}. Above $k_c$, one sees the separate single magnon peak, broadened by
the finite run time and Fourier transform. By $k=0.2\pi$ the peak has disappeared, and one sees a characteristic
small $k$ line shape. Close to $k_c$ the broadening from the finite maximum time obscures the details of
the spectrum.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.7\hsize,scale=1.0]{Snearkc}
\caption{$S(k,\omega)$ near $k_c \approx 0.23\pi-0.24 \pi$ for the 4th order run using the linear prediction method.
The vertical line associated with each $k$ is  the two magnon lower band edge. }
\label{Snearkc}
\end{center}
\end{figure}

The fitting method is very useful to capture the behavior of $S(k,\omega)$ near $k_c$ more accurately.
Just above $k_c$ we expect a combination of a single magnon $\delta$ function peak and a two-magnon
continuum similar to that at small $k$. To fit the time data as accurately as possible, we assume
that the two magnon band starts exactly at the expected threshold
$\omega_{\rm th}= \min_{k'} \epsilon(\pi-k/2+k') + \epsilon(\pi-k/2-k')$. The relevant values of $\epsilon(k)$
used to determine $\omega_{\rm th}$ are far from $k_c$ and  a simple exponential fit determines the
magnon peak location very accurately. We also assume that the small $k$ threshold behavior
$(\omega-\omega_{\rm th})^{1/2}$ applies.  Including, in addition, the next expansion term
$(\omega-\omega_{\rm th})^{3/2}$, we utilize the fitting form (cf Eq. (\ref{singtt}))
\begin{eqnarray}
&A \exp(-i \bar \omega t) +
B \exp(-i \omega_{\rm th} t) (b+it)^{-3/2} \nonumber\\&+ C \exp(-i \omega_{\rm th} t) (c+it)^{-5/2}.
\label{fitkc}\end{eqnarray}
Results from this fitting followed by Fourier transforming are shown in  Fig. \ref{Snearkcfit}. The
fits deviated from the data over the range $t=10-20$ typically by a few times $10^{-5}$, and the
magnitude of the data points fitted to was typically near $0.1$---an excellent fit, making a convincing case
that the assumed asymptotic form is correct and that the results for $S(k,\omega)$ are very accurate.
The peak locations from these fittings were used in the determination of the dispersion relation of
Section III near $k_c$.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.7\hsize,scale=1.0]{Snearkcfit}
\caption{$S(k,\omega)$ near $k_c \approx 0.23\pi-0.24 \pi$ for the 4th order run using fitting to the asymptotic
time decay. }
\label{Snearkcfit}
\end{center}
\end{figure}

One can view
the split-off of the single magnon peak from the two particle continuum at $k_c$ in two different ways.
First, in an approach motivated by the small $k$ two-magnon nature of the spectrum, one can
regard the splitoff  as due to the formation
of a sharp 2-magnon bound state for $k>k_c$.\cite{huse}
Second, if we imagine that $k_c$ is ``large'' so that the two magnon description
is inappropriate, then we can regard the single magnon peak as surviving for $k$ near but less than $k_c$, but with
a broadening caused by decays into two magnons.  Here we will consider in more detail these two pictures.

Consider first the two-magnon bound state picture.
Since the quantum numbers of the single magnon state with $S_z=1$
and a two magnon bound
state composed of a $S_z=0$ and a $S_z=1$ magnon are identical, we are free to regard the excitation which splits
off in either way.  Well above $k_c$ the two-magnon bound state picture is clearly
not very useful, since its formation would imply large
magnon-magnon interactions which are not otherwise observed.
The NL$\sigma$M
does not have any bound states.  Furthermore,
in {\it any} Lorentz invariant theory, increasing the centre of mass
momentum can never lead to bound state formation. As was
discussed in Ref. \onlinecite{huse}, since the bosonic magnon excitations
form a triplet, and any excitation produced from the singlet ground state by
the spin operators must also be a triplet, it follows that
the bound state wave-function must be antisymmetric in its
spatial coordinates. In  Ref. \onlinecite{huse}, it was found that the
effective magnon-magnon interaction is {\it attractive}
in the antisymmetric channel. On the other hand, it is
apparently repulsive in the symmetric channel. This is
related to the BEC picture of the transition at a
critical magnetic field\cite{sorensen} where the Zeeman energy equals $\Delta$.

Even in one dimension, an arbitrarily weak attraction does
not produce a bound state in the antisymmetric channel (although
it does produce one in the symmetric channel.) It was
observed in Ref. \onlinecite{huse} that as $k$ is increased from zero,
the momenta of the two magnons forming the bound state near
the threshold $\omega_{min}(k)$, which are near
$\pi + k/2$, approach the inflection points $k_{in}\approx 1.131\pi$.
To study the bound state near the threshold, we
can expand the dispersion relation, $\epsilon $, near
$\pi +k/2$.  The momenta of the two bosons are:
$\pi +k/2\pm q$.  Expanding the total kinetic
energy in powers of $q$ gives
 the effective
kinetic energy :
\begin{equation} T\approx 2\epsilon +
\epsilon^{(2)}q^2
+{\epsilon^{(4)}\over 12}q^4+\ldots ,\end{equation}
where $\epsilon$ and its derivatives are evaluated at $\pi +k/2$.
We see that the effective mass for the centre of mass motion
is given by:
\begin{equation} {1\over 2m}=\epsilon^{(2)}(\pi +k/2).\end{equation}
$(1/m)$ vanishes at the inflection point, $k\to 2 k_{in}-2\pi\approx
0.262\pi$
and the effective kinetic energy becomes quartic. At
$k_{in}$, the coefficient of the quartic term is
$\epsilon^{(4)}/12\approx 3.341>0$. For a quartic
kinetic energy, an arbitrarily weak attraction leads to
a bound state, in both symmetric {\it and} antisymmetric channel.
(For the antisymmetric case this can be seen by considering
a trial wave-function of $\propto xe^{-x^2/(2w^2)}$.
The potential energy can be assumed to be everywhere
less than a square well of depth $v_0$ and width $a$.
For $w\gg a$, the potential energy is less than
a quantity $\propto -v_0(a/w)^3$,
while the kinetic energy is $\propto (1/w)^4$. For
small $v_0$, this has a negative minimum at $w\propto 1/v_0\gg a$,
proving the existence of an antisymmetric bound state with
binding energy $\propto v_0^4$.)
This argument implies that $k_c<2|\pi -k_{in}|\approx .262\pi$,
since for any attractive potential, as the effective mass
diverges the bound state will eventually form. It is
interesting to note that our DMRG estimate of $k_c\approx 0.23\pi-0.24\pi$
is only very slightly less than $2|\pi -k_{in}|$, suggesting
that the attractive interaction between magnons is weak.
In addition, the splitting of the peak from the two magnon continuum varies as $(m - m_0)^2$ (where $m_0$ is the
mass at $k=k_c$) and thus
as $(k-k_c)^2$ within this bound state picture.

Now consider the second picture of the split-off, that of a magnon entering the continuum, but surviving
in a broadened form near $k_c$. In this case one would expect for $k>k_c$ the single magnon peak and the two
magnon continuum would vary independently with $k$ and that the splitting of the peak from
the two magnon continuum would be linear in $k-k_c$. In Fig. \ref{splitnearkc}, we show this splitting near $k_c$.
Indeed, the splitting appears to be linear in $k-k_c$.  Below $k_c$, the broadening grows very rapidly. If the
splitting is quadratic in $k-k_c$, it must be so only very close to $k_c$.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.7\hsize,scale=1.0]{Delta12}
\caption{Splitting between the single magnon peak and the bottom of the two magnon band
near $k_c$ determined from a fit of the time dependent DMRG data. For $\omega \ge 0.24$, the data was fit well
assuming the single magnon peak was a delta function; below a  better fit was obtained assuming  the peaks was
broadened as a Gaussian. The ``error bars'' indicate the width of the Gaussian. }
\label{splitnearkc}
\end{center}
\end{figure}

The behavior of the splitting versus $k-k_c$ favors the second picture of the split-off.  However, we cannot
rule out that the two-magnon bound state picture applies very close to $k_c$.
Note that
within the
two magnon bound state picture, there does not seem to be a compelling reason for more than
a small fraction of the spectral weight to appear in the bound state. In fact, as shown in Fig. \ref{fratio},
about 50\% of the
spectral weight appears in the single magnon peak near $k\approx 0.25$, but the weight in the peak
is rapidly falling as $k$ is decreased.

It is also interesting to examine the line shape for $k$
slightly less than $k_c$.  In order to compare line shapes for different $k$'s, in Fig. \ref{Sbelowkc} we have
shifted the curves to make the two-magnon thresholds identical, and have scaled them to make them identical at
the  arbitrary point
$\omega_{\rm th}+0.25$.  These curves were made using the linear prediction method. The curve for $k=0.2\pi$
shows a sharp resonance persists below but near $k_c$.
This resonance disappears by the time
$k$ is reduced to $0.14 \pi$.

\begin{figure}[htbp]
\begin{center}
\includegraphics*[width=0.7\hsize,scale=1.0]{Sbelowkc}
\caption{$S(k,\omega)$ near and below $k_c \approx 0.23\pi-0.24 \pi$ for the 2nd order run. Each curve is shifted
by the two-magnon threshold energy $2 \varepsilon(k/2)$, and scaled by an arbitrary factor to make
$S(k,\omega-\omega_{\rm th}=0.25)$
identical in each of the three curves.}
\label{Sbelowkc}
\end{center}
\end{figure}

\section{Conclusions}

The combination of time dependent DMRG and extrapolation of the time dependent correlation functions has
proved to be an extremely effective method for calculating spectral functions for the $S=1$ chain. We have been
able to study fine details of the spectra with much  greater resolution and accuracy than with any
previous method.  In comparing with free boson  and nonlinear sigma model predictions for features of
the spectra near $k=0$ and $k=\pi$, we find good qualitative agreement, but quantitative disagreements
in the overall magnitude of the spectrum and in the high frequency tails.  Our results near $k_c$ where the single
magnon peak enters the two magnon continuum are better described in terms of a single magnon exhibiting
decay and  scattering
below $k_c$ rather than viewing the single magnon peak as the formation of a two magnon bound state above
$k_c$.

We acknowledge very helpful discussions with David Huse.
We acknowledge  support from the NSF under grant DMR-0605444 (SRW), from
NSERC (IA),  and from CIfAR (IA).

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

A widely studied occurrence of the classical
divisor function $d(n)=\sum_{d|n}1$ is in the additive divisor problem, in which one
investigates the asymptotic behavior of the sum
\begin{equation}\label{divisori}
D(x;m)=\sum_{n\leq x}d(n)d(n+m)
\end{equation}
as $x$ tends to infinity and $m$ is a given positive integer.
Ingham was the first to find an asymptotic formula for the above sum in 1927, and in 1931 Estermann brought the Kloosterman sums
into discussion, improving the earlier result. The importance of uniformity in the shift $m$ was noted ten years later by Atkinson \cite{Atkinson},
who further used his findings in studying the fourth power mean of the Riemann
zeta-function on the critical line. In 1982 Deshouillers and Iwaniec applied Kuznetsov's trace formulas into the problem,
obtaining a notable improvement to the upper bound of the error term. Finally, especially crucial for
this paper, Motohashi \cite{Motohashi4} worked out an explicit spectral decomposition for a weighted convolution sum, separating a main term and obtaining an
upper bound for the error term uniformly in the shift $m$. We mention shortly papers by Heath-Brown, Iwaniec \cite{Iwaniec1} and Ivi\'{c}-Motohashi \cite{Ivic2}, \cite{Ivic},
who obtained mean value results for the error term, and \cite{Kuznetsov}, \cite{Meurman} and \cite{Tahtadjan} for other related results.
For a more thorough discussion of the history and the references of the above mentioned classical papers we refer to Motohashi's comprehensive paper \cite{Motohashi4}.

In recent years the analogy between the divisor function and the Fourier coefficients of holomorphic and non-holomorphic cusp
forms has appeared repeatedly in the literature. The behaviour
of the analogous sums to \eqref{divisori},
\begin{equation}\label{holom}
A(x;m)=\sum_{n\leq x}a(n)\overline{a(n+m)}
\end{equation}
involving the Fourier coefficients of a holomorphic cusp form and
\begin{equation}\label{Maass}
T(x;m)=\sum_{n\leq x}t(n)t(n+m)
\end{equation}
involving the Hecke eigenvalues corresponding to Fourier coefficients
of a Maass form, has been studied intensively. In 1965 Selberg studied the meromorphity of a certain Dirichlet series
related to the Fourier coefficients of holomorphic cusp forms, while
Good in 1981 estimated the second moment of a modular $L$-function on the critical line using the spectral decomposition of shifted convolution sums
coming from the Fourier coefficients of a holomorphic cusp form. Jutila \cite{Jutila1}, \cite{Jutila5} derived explicit asymptotic formulae for
all three cases \eqref{divisori}, \eqref{holom} and \eqref{Maass} in a unified way, using the respective generating Dirichlet series. For a discussion of the substantial history including e.g. the works by
Blomer, Duke, Friedlander, Harcos, Iwaniec, Jutila, Michel, Motohashi and Sarnak as well as a general spectral decomposition for the weighted shifted convolution sums we refer to the paper \cite{Harcos} by Blomer and Harcos. Along with these references we mention also \cite{BlomerHarcos}, \cite{Goldfeld}, \cite{Good3}, \cite{Jutila9} and \cite{Jutila6}.
Together the three sums above go under the name of the shifted convolution problem.

As a prologue for this paper, in Lemma 3 of his paper \cite{Jutila9} Jutila studies the sum \eqref{holom}
in mean, and his argument can easily be extended to prove
\begin{equation}\label{aasumma}
\sum_{0\leq f\leq F}\sum_{1\leq n\leq N}\bigg|\sum_{1\leq l\leq L}
a(n+l)\overline{a(n+f+l)}\bigg|^{2}\ll (N+F)^{k}N^{k}L
\end{equation}
for all $N,F\geq 1$ and $1\leq L\leq N$.
Instantaneously the proof gives also an upper bound
for a sum over Hecke eigenvalues of Maass forms in mean, as stated in Lemma 6 of \cite{Jutila7} :  For $N,F\geq 1$ and $1\leq L\leq N$
\begin{equation}\label{Tsumma}
\sum_{0\leq f\leq F}\sum_{1\leq n\leq N}\bigg|\sum_{1\leq l\leq L}
t(n+l)t(n+f+l)\bigg|^ {2}\ll (N+F)^{1+\varepsilon}NL.
\end{equation}
The idea in the estimation of the triple sum lies
in its sensitiveness of the size of the innermost sum over $l$.
A similar result, sensitive for the size of the shift, is needed in the doctoral thesis of the author (\cite{Suvitie}, Lemma 3.4).
An analogous result over more general settings also appears on p. 81 in the paper \cite{Blomer} by Blomer, Harcos and Michel, leading
to the same bounds \eqref{aasumma}, \eqref{Tsumma}.
In all the three papers \cite{Jutila9}, \cite{Jutila7}, \cite{Blomer} the motivation behind these studies has been gaining information about the
upper bounds of the original shifted convolution sums \eqref{holom} and \eqref{Maass}, while the author needed her estimate in her thesis
in estimating a certain spectral sum over inner products involving holomorphic cusp forms and Maass forms.

However, to our knowledge there does not exist a similar bound for the case of the divisor function, and the extension of the proofs of the cuspidal
cases seem problematic due to the main term of the additive divisor problem. This paper is our
attempt to bring new light for this case. Hence we state the following Theorem:
\begin{thm}[Main theorem]\label{thm1}
Let $N\geq 1$, $1\leq L\leq N$ and $1\leq F\ll N^{1-\varepsilon}$. Then
\[
S=\sum_{f\sim F}\sum_{n\sim N}\left|\sum_{l\sim L}d(n+l)d(n+l+f)-\frac{6}{\pi^{2}}\int_{(L+n)/f}^{(2L+n)/f}m(x;f)dx
\right|^{2}
\]
\begin{equation}\label{paatermi}
\ll N^{2+\varepsilon}+N^{1+\varepsilon}LF.
\end{equation}
Here $m(x;f)$ is as in (1.12) in \cite{Motohashi4}.
\end{thm}
The notation $m\sim M$ stands for $M< m\leq2M$.
For the classical main term $M(N,f)$ by Atkinson (\cite{Atkinson}, p. 185) we easily obtain the following upper bound:
\begin{cor}\label{cor2}
Let $N\geq 1$, $1\leq L\leq N$ and $1\leq F\ll N^{1-\varepsilon}$. Then
\[
S=\sum_{f\sim F}\sum_{n\sim N}\left|\sum_{l\sim L}d(n+l)d(n+l+f)-M(2L+n,f)+M(L+n,f)
\right|^{2}
\]
\[
\ll N^{2+\varepsilon}+N^{1+\varepsilon}LF+N^{-1+\varepsilon}L^{2}F^{3}.
\]
\end{cor}

\begin{rem}
With small
values of the quantities $F$ and $L$, that is, if $L^{2}F=o(N)$, the trivial bound $N^{1+\varepsilon}L^{2}F$  yields
a better result than Theorem \ref{thm1} and Corollary \ref{cor2}. On the other hand, in their paper \cite{Ivic} Ivi\'{c} and Motohashi gain the upper bound
\[
\sum_{f=1}^{F}\left|\sum_{n=1}^{N}d(n)d(n+f)-\frac{6}{\pi^{2}}\int_{0}^{N/f}m(x;f)dx
\right|^{2}\ll N^{4/3+\varepsilon}F^{1/3}
\]
uniformly for $F\leq N^{1/2-\varepsilon}$. Here, according to Meurman (\cite{Meurman}, pp. 225-226)
it seems possible to extend the range to $F\leq N$ by replacing the upper bound by $N^{4/3+\varepsilon}F^{1/3}+NF$.
This result yields, up to factor $N^{\varepsilon}$, the same upper bound as in Theorem \ref{thm1} in case $F\gg N^{1/2}$ and $L=N$.
Other than this, to our knowledge, our result \eqref{paatermi} and the trivial bound together yield the genuinely best upper bound with the restriction $F\ll N^{1-\varepsilon}$.
\end{rem}

The main novelty in this paper is the methodological uniformity allowing
an analogous treatment for the cuspidal cases:
Following the analogous proof for the holomorphic cusp forms we immediately gain
\begin{thm}\label{thm3}
Let $N\geq 1$, $1\leq L\leq N$ and $1\leq F\ll N^{1-\varepsilon}$. Then
\[
\sum_{f\sim F}\sum_{n\sim N}\left|\sum_{l\sim L}a(n+l)\overline{a(n+l+f)}
\right|^{2}
\ll
N^{2k+\varepsilon}+N^{2k-1+\varepsilon}LF.
\]
\end{thm}

To compare our result to the earlier mentioned, we notice that up to a factor $N^{\varepsilon}$ and with the restriction
$F\ll N^{1-\varepsilon}$, \eqref{aasumma} yields our result only if $L\asymp 1$, and
otherwise our result is genuinely better.

In the case of non-holomorphic cusp forms we face the problem of the lack of a proper analogy for the spectral decomposition of the shifted convolution sum
in question. There exists an analogous result for our crucial Lemma \ref{Motohas}, namely Lemma 5 in \cite{Motohashi},
originally due to Jutila \cite{Jutila2}, \cite{Jutila3}. However, although the analogous part of this decomposition
results the anticipated analogous bound, there appears also an arithmetic correction term, which leads to additional problems. If we were dealing with
a genuinely oscillating weight function in our proof, then the contribution of the correction term would easily be
proven to stay under the expected bound. However, as this is not the case we shall lighten this paper and just assign the conjecture below,
leaving the proof for future study.
\begin{conj}\label{con}
Let $N\geq 1$, $1\leq L\leq N$ and $1\leq F\ll N^{1-\varepsilon}$. Then
\[
\sum_{f\sim F}\sum_{n\sim N}\left|\sum_{l\sim L}t(n+l)t(n+l+f)
\right|^{2}
\ll
N^{2+\varepsilon}+N^{1+\varepsilon}LF.
\]
\end{conj}

\begin{rem}
Note also Theorem 1 in \cite{Harcos}, which gives a more general spectral decomposition of a shifted convolution sum including
Hecke eigenvalues of the underlying cusp forms of two arbitrary cuspidal automorphic representations of $PGL_{2}(\mathbb{Z})\backslash PGL_{2}(\mathbb{R})$.
\end{rem}

\section{The needed notation and auxiliary lemmas}

\subsection{The cusp forms}

Besides their analogous behavior to the divisor function, the holomorphic and non-holomorphic cusp
forms play a role in the spectral decomposition of the shifted
convolution sum over the divisor function (Lemma \ref{Motohas}). Hence we start by introducing briefly some results concerning them.
For the proofs and for a general reference the reader is referred to Motohashi's monograph \cite{Motohashi2}.

We may confine ourselves to the cusp forms for the full modular group
$\Gamma=SL_{2}(\mathbb{Z})$ operating through M\"{o}bius
transformations on the upper half plane
$\mathbb{H}$. A
\textit{holomorphic cusp form}
$F(z):\mathbb{H}\rightarrow\mathbb{C}$ of weight $k\in\mathbb{Z}$
with respect to $\Gamma$ can be represented by its Fourier series
\[
F(z)=\sum^{\infty}_{n=1}a(n)e(nz), \;\; e(\alpha)=\exp(2\pi i
\alpha).
\]
We may assume that $k$ is even and $k\geq 12$, otherwise $F(z)$ is
trivial.
We let
\[
\{\psi_{j,k}\;|\; 1\leq j\leq\vartheta(k)\}
\]
be an orthonormal basis of the unitary space of holomorphic cusp
forms of weight $k$, and write
\[
\psi_{j,k}(z)=\sum_{n=1}^{\infty}\rho_{j,k}(n)n^{\frac{k-1}{2}}e(nz).
\]
We may suppose that the basis vectors are eigenfunctions of the
Hecke operators $T_{k}(n)$ for all positive integers $n$. Thus, in particular,
$T_{k}(n)\psi_{j,k}=t_{j,k}(n)\psi_{j,k}$ for certain real numbers
$t_{j,k}(n)$, which we call Hecke eigenvalues. Comparing Fourier
coefficients on both sides, one may verify that
$\rho_{j,k}(n)=\rho_{j,k}(1)t_{j,k}(n)$ for all $n\geq 1$, $1\leq
j\leq \vartheta(k)$.
We put $\alpha_{j,k}=16 \Gamma(2k)(4\pi)^{-2k-1}$ $\times|\rho_{j,k}(1)|^{2}$,
whence
\begin{equation}\label{thetaalpha}
\sum_{j=1}^{\vartheta(k)}\alpha_{j,k}\ll k.
\end{equation}
Furthermore, we let
\[
H_{j,k}(s)=\sum_{n=1}^{\infty}t_{j,k}(n)n^{-s}
\]
stand for the Hecke L-function attached to $\psi_{j,k}$. The
series converges absolutely for Re $s>1$ because of the bound
\begin{equation}\label{tee}
|t_{j,k}(n)|\leq d(n)\ll n^{\varepsilon}
\end{equation}
by Deligne \cite{Deligne}. Furthermore $H_{j,k}(s)$ can be continued
to an entire function, and it satisfies a functional equation, which implies that for bounded $s$
\begin{equation}\label{hoojiikoo}
H_{j,k}(s)\ll k^{c}
\end{equation}
uniformly in $j$. Here $c$ is a suitable constant, which depends only on  Re $s$.

A \textit{non-holomorphic cusp form} $u(z)=u(x+iy):\mathbb{H}\rightarrow\mathbb{C}$ is a
non-constant real-analytic $\Gamma$-invariant function in the
upper half-plane, square-integrable with respect to the hyperbolic
measure $d\mu(z)=\frac{dx\,dy}{y^{2}}$ over a fundamental domain
of $\Gamma$. Furthermore $u(z)$ is an eigenfunction of the non-euclidean Laplacian
 $\Delta=-y^2(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial
 y^2})$, and the corresponding eigenvalue can be written
as $1/4+\kappa^2$ with $\kappa>0$.
The Fourier series expansion for $u(z)$ is then of the form
\[
u(z)=y^{1/2}\sum_{n\neq 0}\rho(n)K_{i\kappa}(2\pi |n|y)e(nx)
\]
with $K_{\nu}$ a Bessel function of imaginary argument. We may
suppose that our cusp forms are eigenfunctions of the Hecke
operators $T(n)$ for all positive integers $n$ and that $u(x+iy)$
is even or odd as a function of $x$. Thus $T(n)u=t(n)u$ for
certain real numbers $t(n)$, which are again called Hecke eigenvalues,
and $u(-\overline{z})=\pm u(z)$. Comparing Fourier coefficients on
both sides, one may verify that $\rho(n)=\rho(1)t(n)$ and
$\rho(-n)=\pm\rho(n)$ for all $n\geq 1$.
The \textit{Maass (wave) forms} $u_{j}$
constitute an orthonormal set of non-holomorphic cusp forms
arranged so that the corresponding parameters $\kappa_{j}$
determined by the eigenvalues $1/4+\kappa_{j}^{2}$ lie in an
increasing order. We write $\rho_{j}(n)$ and $t_{j}(n)$ for the corresponding Fourier coefficients and Hecke
eigenvalues.
Now let
\[
\alpha_{j}=|\rho_{j}(1)|^{2}/\cosh(\pi\kappa_{j}),
\]
and let
\[
H_{j}(s)=\sum_{n=1}^{\infty}t_{j}(n)n^{-s}
\]
stand for the Hecke $L$-function attached to $u_{j}$. As the Ramanujan-Petersson conjecture
$t_{j}(n)\ll n^{\varepsilon}$ holds in a mean value sense, the series converges
absolutely for Re $s>1$. Again, $H_{j}(s)$ can be continued to an entire function.
We have the bound
\begin{equation}\label{hoonelja}
\sum_{\kappa_{j}\leq K}\alpha_{j}H_{j}^{4}\left(\frac{1}{2}\right)\ll K^{2}\log^{15}K.
\end{equation}
(See \cite{Motohashi2}, Theorem 3.4.)

Lastly, in this paper, the following notation will be adopted:
Vinogradov's relation $f(z)\ll g(z)$ is another notation for
$f(z)=\mathcal{O}(g(z))$. We let $\varepsilon$ stand generally for a small positive number,
not necessarily the same at each occurrence.

\subsection{The needed lemmas}

In this section we shall gather some auxiliary results which will be used
during the course of the proofs of the theorems.

For estimating the
Gamma function we shall frequently use the following
Stirling's formula without further reference:
In any fixed strip $b\leq\sigma\leq c$ we have
\[
|\Gamma(\sigma+i
t)|=\sqrt{2\pi}|t|^{\sigma-1/2}e^{-|t|\pi/2}(1+\mathcal{O}(|t|^{-1})),
\]
for $|t|\rightarrow \infty$.
For a proof, see Olver \cite{Olver} p. 294.
Further we have
\begin{equation}\label{kertoma}
\Gamma(x)= \sqrt{2\pi}x^{x-1/2}e^{-x}(1+\mathcal{O}(x^{-1})),
\end{equation}
for all positive real $x\rightarrow \infty$, see \cite{Lebedev}, Eq. (1.4.25).

The proof of the following classical estimate for the fourth moment of Riemann's zeta-function on the critical line
can be found e.g. from \cite{Titchmarsh}, p. 147:
\begin{equation}\label{zetalog}
\int_{0}^{T}\left|\zeta\left(\frac{1}{2}+it\right)\right|^{4}\,dt\ll T \log^{4} T
\end{equation}
for $T\geq 1$. Along with this estimate we have
\begin{equation}\label{zeta12}
\int_{0}^{T}\left|\zeta\left(\frac{1}{2}+it\right)\right|^{12}\,dt\ll T^{2} \log^{17} T
\end{equation}
for $T\geq 1$ by Heath-Brown \cite{Heath-Brown}.
Moreover it is known that
\[
\frac{1}{\zeta(1+it)}\ll \log |t|,
\]
as $|t|\geq 1$. For a proof, see e.g.
\cite{Titchmarsh}, p. 132.

We
have an important tool arising from spectral theory:
\begin{lem}[The spectral large sieve]
For $K\geq 1$, $1\leq \Delta\leq K$, $M\geq 1$ and any complex
numbers $a_{m}$ we have
\[
\sum_{K\leq\kappa_{j}\leq K+\Delta}\alpha_{j}\left|\sum_{m\leq
M}a_{m}t_{j}(m)\right|^{2}\ll
(K\Delta+M)(KM)^{\varepsilon}\sum_{m\leq M}|a_{m}|^{2}.
\]
\end{lem}
For a proof, see Theorem 1.1 in \cite{Jutila4} or Theorem 3.3 in
\cite{Motohashi2}.
The continuous analogy is the following:
\begin{lem}
For $K$ real, $\Delta\geq 0$, $M\geq 1$ and any complex
numbers $a_{m}$ we have
\[
\int_{K}^{K+\Delta}\left|\sum_{m\leq
M}a_{m}\sigma_{2ir}(m)m^{-ir}\right|^{2}\, dr\ll
(\Delta^{2}+M)M^{\varepsilon}\sum_{m\leq M}|a_{m}|^{2}
\]
uniformly in $K$, as $\sigma_{2ir}(f)=\sum_{d|f}d^{2ir}$.
\end{lem}
Proof can be found in \cite{Suvitie}, Lemma 1.12.

A crucial role is played by a spectral decomposition of the
shifted convolution sum over the divisor function:
\begin{lem}\label{Motohas}
Let $f$ be a positive integer and $W$ a smooth function of compact
support on $(0,\infty)$. Then
\[
\sum_{l=1}^{\infty}d(l)d(l+f)W\left(\frac{l}{f}\right)=\frac{6}{\pi^{2}}\int_{0}^{\infty}
m(x;f)W(x)\, dx
+\frac{f^{1/2}}{\pi}\int_{-\infty}^{\infty}f^{-ir}\sigma_{2ir}(f)
\]
\[
\times \frac{|\zeta(\frac{1}{2}+ir)|^{4}}{|\zeta(1+2ir)|^{2}}\Theta(r;W)\,dr
+f^{1/2}\sum_{j=1}^{\infty}\alpha_{j}t_{j}(f)H_{j}^{2}\left(\frac{1}{2}\right)\Theta(\kappa_{j};W)
\]
\begin{equation}\label{spekthajo}
+\frac{1}{4}f^{1/2}\sum_{k=6}^{\infty}\sum_{j=1}^{\vartheta(k)}(-1)^{k}\alpha_{j,k}t_{j,k}(f)H_{j,k}^{2}\left(\frac{1}{2}\right)\Xi\left(k-\frac{1}{2};W\right)
\end{equation}
Here $m(x;f)$ is defined by (1.12) in \cite{Motohashi4}, $\sigma_{2ir}(f)$ is as in the previous lemma
and
\[
\Theta(r;W) =\frac{1}{2}\textrm{Re}\left\{\left(1+\frac{i}{\sinh(\pi r)}\right)\Xi(ir;W)\right\}
\]
\begin{equation}\label{theeta}
=\frac{1}{4}\left(1+\frac{i}{\sinh(\pi r)}\right)\Xi(ir;W)+(r\mapsto -r),
\end{equation}
with
\[
\Xi(\xi;W)=\frac{\Gamma(\xi+\frac{1}{2})^{2}}{\Gamma(2\xi+1)}
\int_{0}^{\infty}x^{-1/2-\xi}
F\left(\xi+\frac{1}{2},\xi+\frac{1}{2};2\xi+1;-\frac{1}{x}\right)
\]
\begin{equation}\label{ksiii}
\times W(x)
\,dx.
\end{equation}
Here $F(. \, , .\, ; . \, ; .)$ is the hypergeometric function and $(r\mapsto -r)$ stands for an expression similar to the
preceding one, but with $r$ replaced by $-r$.
\end{lem}
For the proof, see \cite{Motohashi4}, Theorem 3.

Finally we introduce the basic inequality in the proof of the
classical large sieve:
\begin{lem}[Sobolev]\label{klassinen}
Let $a\leq u\leq a+\Delta$ and let the function $f$ be
continuously differentiable on this interval. Then
\[
|f(u)|^{2}\leq \Delta^{-1}\int_{a}^{a+\Delta}|f(x)|^{2}\,dx
+2\left(\int_{a}^{a+\Delta}|f(x)|^{2}\,dx\right)^{1/2}
\]
\[
\times \left(\int_{a}^{a+\Delta}|f'(x)|^{2}\,dx\right)^{1/2} \ll
\Delta^{-1}\int_{a}^{a+\Delta}\left(|f(x)|^{2}+\Delta^{2}|f'(x)|^{2}
\right)\,dx
\]
uniformly.
\end{lem}
For a proof, see Montgomery \cite{Montgomery}, Lemma 1.1 applied
to $f^{2}$.

\section{Proofs of the results}

\subsection{Theorem \ref{thm1}}

We start by adding a set of $\asymp\log(1/U)$ real-valued smooth weight functions
$g_{\delta}\left(\frac{l}{L}\right)$ to the $l$-sum so that
their sum produces an approximation of the characteristic function of the interval $[1,2]$ with an error of size $\ll U$
and their supports
widen step by step by factors 2 when we move away from the end points 1 and
2. Thus if the supports of the first and last weight functions are of length $U=L^{-1+\varepsilon}$, then the next
ones are of length $\asymp 2U$ and so on, and the slopes of the weight functions cancel out each other.
Hence we let $g_{\delta}(x)$ stand for a compactly supported function on some interval of length $\asymp\delta$ contained in $[1,2]$. Moreover,  $g(x)=1$ on
an interval of length $\asymp\delta$ and
$g^{(\nu)}(x)\ll_{\nu}\delta^{-\nu}$ for each $\nu\geq 0$ and $x\in \mathbb{R}$. Here $L^{-1+\varepsilon}\leq\delta\ll 1$. Now
\[
S\ll L^{\varepsilon}\sum_{f\sim F}\sum_{n\sim N}\bigg |\sum_{l=1}^{\infty}d(n+l)d(n+l+f)W\left(\frac{n+l}{f}\right)
\]
\[
-\frac{6}{\pi^{2}}\int_{0}^{\infty}m(x;f)
 W(x)dx\bigg|^{2}+N^{1+\varepsilon}FU^{2}L^{2}
\]
\[
\ll L^{\varepsilon}\sum_{f\sim F}\sum_{n\sim N}\left| e_{1}(n,f;L)+ e_{2}(n,f;L)+ e_{3}(n,f;L)\right|^{2}
+N^{1+\varepsilon}F,
\]
with $W(x)=g_{\delta}\left(\frac{fx-n}{L}\right)$. Here $e_{\nu}(n,f;L)$ ($\nu=1,2,3$) stands for the $(\nu+1)$th term on the right hand side of \eqref{spekthajo}.

The contribution of the term  $e_{2}(n,f;L)$ turns out to be the most difficult to bound. Let us denote
\[
S_{2}=\sum_{f\sim F}\sum_{n\sim N}\left| f^{1/2}\sum_{j=1}^{\infty}\alpha_{j}t_{j}(f)H_{j}^{2}\left(\frac{1}{2}\right)\Theta(\kappa_{j};W)\right|^{2}
\]
with $\Theta(r;W)$ as in \eqref{theeta}. As the support of $W$ tends to
infinity, we have in the integral in \eqref{ksiii}
\[
F\left(\frac{1}{2}+ir,\frac{1}{2}+ir;1+2ir;-\frac{1}{x}\right)
=
\left(\frac{1+\sqrt{1+x^{-1}}}{2}\right)^{-1-2ir}
\]
\begin{equation}\label{hyperg}
\times F\left(\frac{1}{2}+ir,\frac{1}{2};1+ir;\left(\frac{1-\sqrt{1+x^{-1}}}{1+\sqrt{1+x^{-1}}}\right)^{2}\right)
\end{equation}
by the quadratic transformation formula for the hypergeometric
function from Lebedev \cite{Lebedev}, Eq. (9.6.12). We write the RHS of \eqref{hyperg} by the definition of the hypergeometric series.
Now clearly the $k$th term in the series is non-oscillating and of order $x^{-2k}$. Therefore it suffices to study
only the leading term $1$, others behaving similarly until some sufficient constant, and the contribution of the rest being negligibly small.
(See also Convention 2 in \cite{JutilaMotohashi}.)

By integrating repeatedly by parts we notice that $\kappa_{j}$ can be
truncated to be $\ll N^{1+\varepsilon}\delta^{-1}L^{-1}$ with negligibly small error. We further decompose the range
$[1,N^{1+\varepsilon}\delta^{-1}L^{-1}]$ into dyadic intervals $\kappa_{j}\sim K$ with $1\leq K\ll N^{1+\varepsilon}\delta^{-1}L^{-1}$, the number
of these intervals being $\ll N^{\varepsilon}$.
Next we prepare ourselves to use the duality principle (see for example pp. 169-170 in \cite{Kowalski}).
Let $b_{f,n}$ stand for any finite sequence of complex numbers such that $\sum_{f\sim F}\sum_{n\sim N}|b_{f,n}|^{2}=1$.
Then
\[
S_{2}\ll N^{\varepsilon}L^{2}F^{-1}\sum_{\kappa_{j}\sim K}\alpha_{j}\left| \sum_{f\sim F}\sum_{n\sim N}b_{f,n}t_{j}(f)I(\kappa_{j};W)\right|^{2}
\sum_{\kappa_{j}\sim K}\alpha_{j}
H_{j}^{4}\left(\frac{1}{2}\right)\kappa_{j}^{-1}
\]
with
\[
I(\kappa_{j};W)=\int_{1}^{2}\left(\frac{zL+n}{f}\right)^{-1/2-i\kappa_{j}}
\left(\frac{1+\sqrt{1+\frac{f}{zL+n}}}{2}\right)^{-1-2i\kappa_{j}}g_{\delta}(z)\, dz.
\]
Note that  $\kappa_{j}$ replaced by $-\kappa_{j}$ can be taken care of with complex conjugation.
By \eqref{hoonelja} this is further
\begin{equation}\label{valiin}
\ll  N^{\varepsilon}KL^{2}F^{-1}\sum_{\kappa_{j}\sim K}\alpha_{j}\left| \sum_{f\sim F}\sum_{n\sim N}b_{f,n}t_{j}(f)I(\kappa_{j};W)\right|^{2}.
\end{equation}

We want to next apply the spectral large sieve, so we shall use Sobolev's lemma \ref{klassinen} as done in Jutila's paper \cite{Jutila1} on p. 454 to get rid of
the parameter
$\kappa_{j}$ in the integral $I(\kappa_{j};W)$. Therefore
the range $(K,2K]$ for $\kappa_{j}$ is
split up into segments of length $\Delta$ in such a way that the
integral $I(y;W)$ remains essentially stationary as $y$ runs
over a segment. That is
\[
\Delta \log\left(\frac{(zL+n)\left(1+\sqrt{1+\frac{f}{zL+n}}\right)^{2}}{4f}\right)\ll \log N.
\]
In this way, the second term in the upper bound in Lemma
\ref{klassinen} will be comparable to the first. An appropriate choice would be
$\Delta=1$.

Now we divide the $\kappa_{j}$-sum in \eqref{valiin} into
subsums of length $\Delta$, and apply Lemma \ref{klassinen} to each subsum
arriving at the bound
\[
N^{\varepsilon}KL^{2}F^{-1}\sum_{l=0}^{
K-1}\int_{K+l}^{K+l+1}
\sum_{K+l<\kappa_{j}\leq K+l+1}\alpha_{j}\bigg(\bigg| \sum_{f\sim F}\sum_{n\sim N}b_{f,n}t_{j}(f)I(y;W)\bigg|^{2}
\]
\[
+\bigg| \sum_{f\sim F}\sum_{n\sim N}b_{f,n}t_{j}(f)\frac{\partial}{\partial y}I(y;W)\bigg|^{2}\bigg)\, dy.
\]

We next apply the spectral large sieve to the subsum over
$\kappa_{j}$, and finally add the results
corresponding to all subsums. This leads us to the bound
\[
N^{\varepsilon}KL^{2}F^{-1}(K+F)
\sum_{f\sim F}\int_{K}^{2K}\bigg(\bigg|\sum_{n\sim N}b_{f,n}I(y;W)\bigg|^{2}
\]
\[
+\bigg| \sum_{n\sim N}b_{f,n}\frac{\partial}{\partial y}I(y;W)\bigg|^{2}\bigg)\, dy.
\]
Now we use Cauchy's inequality for the integrals $I(y;W)$ and $\frac{\partial}{\partial y}I(y;W)$, and attach
a suitable real-valued smooth weight function
$u(y)$ to the $y$-integral. We set $u$ to be compactly supported on the interval $[K/2,3K]$, $u(y)=1$ on
$[K,2K]$ and $u^{(\nu)}(y)\ll_{\nu}K^{-\nu}$ for each $\nu\geq 0$ and $y\in \mathbb{R}$. By opening the squares
we end up with a double sum $\sum_{n_{1}\sim N}\sum_{n_{2}\sim N}$. By repeated partial integration over the $y$-integral we now
see that we may truncate $n_{1}-n_{2}\ll N^{1+\varepsilon}K^{-1}$. Hence we finally have an upper bound
\[
N^{-1+\varepsilon}K^{2}L^{2}(K+F)\delta^{2}
\sum_{f\sim F}\sum_{n\sim N}|b_{f,n}|^{2}\sum_{|n-n'|\ll N^{1+\varepsilon}K^{-1}}1,
\]
and we conclude with the desired upper bound.

The contribution of $e_{1}(n,f;L)$ to $S$ can be estimated by exactly the same argumentation as above. We need the estimation
\begin{equation}\label{zeta8}
\int_{0}^{T}\left|\zeta\left(\frac{1}{2}+it\right)\right|^{8}\,dt\ll T^{3/2+\varepsilon}
\end{equation}
following readily of \eqref{zetalog} and \eqref{zeta12} by the Cauchy's inequality. Hence
\[
S_{1}\ll N^{1+\varepsilon}LF.
\]
We could use also a more straightforward argumentation, but these steps make it easier to follow the proof of Theorem \ref{thm3}.

Using \eqref{thetaalpha}, \eqref{tee}, \eqref{hoojiikoo} and \eqref{kertoma} the trivial estimates suffice for the estimation of $e_{3}(n,f;L)$.

\subsection{Corollary \ref{cor2}}

Using straightforward calculation it is easily seen that the difference of our main term and that of Atkinson's is
\[
\frac{6}{\pi^{2}}\int_{(L+n)/f}^{(2L+n)/f}m(x;f)dx-
(M(2L+n,f)-M(L+n,f))\ll FLN^{-1+\varepsilon},
\]
whence the result follows.

\subsection{Theorem \ref{thm3}}

Analogous to the proof of Theorem \ref{thm1}. For the size of the Fourier coefficients see \eqref{tee}. Instead of Lemma \ref{Motohas} we use Lemma 4 in \cite{Motohashi},
originally due to Motohashi \cite{Motohashi5}. The counterpart for \eqref{hoonelja} is Lemma 4 in \cite{Jutila1}, and instead of \eqref{zeta8} we need
Lemma 3 in \cite{Jutila1}.

\subsection*{Acknowledgment}
The author wishes to thank Professor Matti Jutila for the idea for this
subject as well as for pointing out the advantage of choosing the weight function in the proof of Theorem \ref{thm1} as it was done, leading to
considerable improvement of the results. His help in finding a proof of the duality principle easily extendable to the continuous case was
of great value. The kind support and advices of Gergely Harcos, Tom Meurman, Philippe Michel and Yoichi Motohashi are gratefully acknowledged.