text:
\section{Introduction}

A key insight of statistical mechanics is that {\em
equilibrium} states can be accurately described in terms of only a
small number of thermodynamic variables, such as temperature and
pressure. For non-equilibrium systems such as glasses no similar
simplification exists {\em a
priori}; the whole past history of a sample is in principle required
to specify its state at a given time. This complexity makes
theoretical analysis awkward, and one is led
instead to look for a description of non-equilibrium states in
terms of a few effective thermodynamic parameters. Much work in recent
years has focussed on one such parameter, the {\em effective
temperature}. This can be defined on the basis of fluctuation-dissipation
(FD) relations between correlation and response functions and has
proved to be very fruitful in mean field
systems~\cite{CugKurPel97,CriRit03}.

The use of FD relations to quantify the out-of-equilibrium dynamics in
glassy systems is motivated by the occurrence of {\em
aging}~\cite{BouCugKurMez98}: the time scale of response to an
external perturbation increases with the age (time since preparation)
$t_{\rm w}$ of the system. As a consequence, time translational invariance
and the equilibrium fluctuation-dissipation
theorem~\cite{Reichl80} (FDT) relating correlation and response
functions break down. To quantify this, one considers the correlation
function of two generic observable $A$ and $B$ of a system, defined as
\begin{equation}
C(t,t_{\rm w})=\left\langle A(t)B(t_{\rm w})\right\rangle -\left\langle A(t) \right\rangle \left\langle B(t_{\rm w}) \right\rangle
\label{corr}
\end{equation}
The associated (impulse) response function can be defined as
\[
R(t,t_{\rm w})=\left.\frac{\delta \left\langle A(t)\right\rangle}{\delta h_B(t_{\rm w})}\right|_{h_B=0}
\]
and gives the linear response of $A$ at time $t$ to a small impulse in
the field $h_B$ conjugate to $B$ at the earlier time $t_{\rm w}$. (The
latter is normally thought of as a ``waiting time'' since preparation
of the system at time 0.) Equivalently one can work with the
susceptibility
\begin{equation}
\chi(t,t_{\rm w})=\int_{t_{\rm w}}^t\! dt'\, R(t,t')
\label{switch_on}
\end{equation}
which encodes the response of $A(t)$ to a small step
$h_B(t)=h_B\Theta(t-t_{\rm w})$ in the field starting at $t_{\rm w}$. In
equilibrium, FDT implies that
$-\partial_{t_{\rm w}} \chi(t,t_{\rm w}) = R(t,t_{\rm w}) =
{T}^{-1}\partial_{t_{\rm w}}C(t,t_{\rm w})$. Out of equilibrium, the
violation of FDT can be measured by an FD ratio (FDR), $X$, defined
through~\cite{CugKur93,CugKur94}
\begin{equation}
\label{eqn:non_eq_fdt}
-\partial_{t_{\rm w}} \chi(t,t_{\rm w}) = R(t,t_{\rm w}) =
\frac{X(t,t_{\rm w})}{T}\partial_{t_{\rm w}}C(t,t_{\rm w})
\end{equation}
This implies that $X$ can be read off from the slope $-X/T$ of a
parametric FD plot showing $\chi$ vs $C$, at fixed $t$ and with $t_{\rm w}$
as the curve parameter. This remains the case if both axes are
normalized by the equal-time variance of $A$, $C(t,t)$, a procedure
which is helpful
in fixing the scale of the plot in situations where $C(t,t)$ varies
significantly with time~\cite{FieSol02,SolFieMay02}. In equilibrium, the FD
plot is a straight line with slope $-1/T$.

In mean-field spin glasses~\cite{CugKurPel97,CugKur93,CugKur94} one
finds that FD plots of autocorrelation and response of local spins and
similar observables approach a limiting shape for large $t$. This is
typically composed of
two straight line segments. In the first of these one finds $X=1$,
corresponding to
quasi-equilibrium dynamics for time differences $t-t_{\rm w}$ that do not
grow with the age of the system. The second line segment has $X<1$ and
reflects the dynamics on aging timescales, i.e.\ time differences
growing (in the simplest case linearly) with $t_{\rm w}$. One can use this
to define a non-equilibrium {\em effective temperature} $T_{\rm
eff}=T/X$, which has been shown to have many of the properties of a
thermodynamic temperature~\cite{CugKurPel97,CugKur93,CugKur94}.

How well this physically very attractive mean-field scenario transfers
to more realistic non-equilibrium systems with short-range
interactions has been a matter of intense research
recently~\cite{CriRit03}. A useful class of systems for studying this
question in detail is provided by ferromagnets quenched from high
temperature to the critical temperature (see e.g.\
Refs.~\cite{GodLuc00b,MayBerGarSol03} and the recent
review~\cite{CalGam05}) or below. The system then coarsens -- by the
growth of domains with the equilibrium magnetization, for $T<T_c$ --
and exhibits aging; in an infinite system equilibrium is never
reached. The aging is clearly related to the growth of a
lengthscale~\cite{Bray94} (domain size for $T<T_{\rm c}$, or correlation
length for $T=T_{\rm c}$), and this makes ferromagnets attractive
``laboratory'' systems for understanding general properties of
non-equilibrium dynamics. They are of course not completely generic;
compared to e.g.\ glasses they lack features such as thermal
activation over energetic or entropic barriers.

We focus in this paper mostly on ferromagnets quenched to $T_{\rm c}$, i.e.\
on critical coarsening dynamics. Some care is needed in this case with
the interpretation of limiting FD plots: while in the mean-field
situation $X$ becomes at long times a function of $C$ only, as implied
by the existence of a nontrivial limit plot, in critical coarsening
$X$ approaches a function of $t/t_{\rm w}$~\cite{GodLuc00b}. In the
interesting regime where $t/t_{\rm w}$ is finite but $>1$, time differences
$t-t_{\rm w}={{\mathcal O}}(t_{\rm w})$ are then large, and e.g.\ spin autocorrelation
functions have decayed to a small part of their initial value. In the
limit $t,t_{\rm w}\to\infty$ the FD plot then assumes a pseudo-equilibrium
shape, with all nontrivial detail compressed into a vanishing region
around $C=0$.

The fact that the FDR is a smooth function of $t/t_{\rm w}$ makes the
interpretation of $T/X$ as an effective temperature less obvious than
in mean-field spin glasses, where $T/X$ is constant within each time
sector ($t-t_{\rm w}={{\mathcal O}}(1)$ vs $t-t_{\rm w}$ growing with $t_{\rm w}$). To eliminate
the time-dependence one can consider the limit of times that are both
large and well-separated. This defines an {\em asymptotic FDR}
\begin{equation}
X^\infty=\lim_{t_{\rm w}\to\infty}\lim_{t\to\infty}X(t,t_{\rm w})
\end{equation}
An important property of this quantity is that it should be {\em
universal}~\cite{GodLuc00b,CalGam05} in the sense that its value is
the same for different systems falling into the same universality
class of critical non-equilibrium dynamics. This makes a study of
$X^\infty$ interesting in its own right, even without an
interpretation in terms of effective temperatures.

If one nevertheless wants to pursue such an interpretation, the
resulting value of
the effective temperature $T/X^\infty$ should be the same for all (or
at least a large class of) observables $A$. The observable-dependence
of $X^\infty$ therefore becomes a key
question~\cite{FieSol02,SolFieMay02,MayBerGarSol03,CalGam04}.
Conventionally, much work on non-equilibrium ferromagnets has focussed
on the local spin autocorrelation function and associated response. An
obvious alternative is the long-wavelength analogue, i.e.\ the
correlation function of the fluctuating magnetization. Exact
calculations for the Ising chain~\cite{MayBerGarSol03,MaySol04} as
well as numerical simulations~\cite{MayBerGarSol03,MayBerGarSol04} in
dimension $d=2$ show that the resulting global $X^\infty$ is always
identical to the local version. This local--global correspondence,
which can also be obtained by field-theoretic
arguments~\cite{CalGam05,CalGam02c,CalGam02}, arises physically
because the long wavelength Fourier components of the spins are
slowest to relax and dominate the long-time behaviour of both local
and global quantities.

The local--global correspondence does of course not address the full
range of observable-dependence of the asymptotic FDR; one might ask
about other observables which are linear combinations not of spins
but for example products of interacting spins. In the critical Ising
model in $d=2$, numerical
simulations~\cite{MayBerGarSol03,MayBerGarSol04} suggest that even
these give the same $X^\infty$, so that an interpretation of
$T/X^\infty$ in terms of an effective temperature appears
plausible. One of the motivations for the current study was to verify
whether this {\em observable-independence} of $X^\infty$ across
different types of observables holds in an exactly solvable model, the
spherical ferromagnet~\cite{BerKac52,Joyce72}. In addition, we will study what
effect different {\em initial conditions} have on $X^\infty$. This is
motivated by our recent study of Ising models in the classical regime
of large $d$,
where critical fluctuations are irrelevant~\cite{GarSolPagRit05}. It
turned out that magnetized initial do states produce a different value of
$X^\infty$, so that critical coarsening in the presence of a nonzero
magnetization is in a new dynamical universality class even though the
magnetization does decay to zero at long times.

We begin in Sec.~\ref{sec:Gaussian} with a brief review of the
standard setup for the dynamics of the spherical model, as used in
e.g.~\cite{GodLuc00b}. Fluctuations in an effective Lagrange
multiplier enforcing the spherical constraint are neglected, leading
to a theory where all spins are Gaussian random variables. In
Sec.~\ref{sec:finite-range} this is applied to various observables of
finite range, by which we mean correlations and responses probing
lengthscales that can be large but remain small compared to the
overall system size. For spin observables we show that the expected
equality of $X^\infty$ between local and long-range quantities holds
(Sec.~\ref{sec:spin_observables}). We check observable-independence of
$X^\infty$ further by considering bond and spin product observables,
in Sec.~\ref{sec:bond_observables} and~\ref{sec:product_observables}
respectively.

The major part of the paper is then devoted to a study of FDRs for
{\em global} observables, with a focus on the energy, i.e.\ the global
bond observable. Because of the weak infinite-range interaction
generated by the spherical constraint, such observables behave
differently from their long-range analogues in the spherical
model. Calculations of correlation and response functions are
technically substantially more difficult because Lagrange multiplier
fluctuations can no longer be neglected. To account for them we construct
in Sec.~\ref{sec:setup} a systematic expansion of the dynamically
evolving spins in $N^{-1/2}$. This allows us to calculate the leading
non-Gaussian corrections that we need for global correlations, as
shown for the case of the energy in
Sec.~\ref{sec:energy_general}. After a brief digression to equilibrium
dynamics, we evaluate the resulting expressions in
Sec.~\ref{sec:noneq_large_d} for $d$ above the critical dimension
$d_{\rm c}=4$, and in Sec.~\ref{sec:noneq_small_d} for
$d<4$. Importantly, we will find that in the latter case the
asymptotic FDR is different from those for finite-range
observable. This means that an effective temperature interpretation of
$X^\infty$ is possible at best in a very restricted sense. However, we
will find that our results are in agreement with recent
renormalization group (RG) calculations near $d=4$~\cite{CalGam04} in
the $O(n\to\infty)$-model. This suggests that the non-Gaussian
effects captured in global observables are important for linking the
spherical model to more realistic systems with only short-range
interactions. Finally, we turn in Sec.~\ref{sec:magnetized} to
critical coarsening starting from magnetized initial conditions. Here
already the global {\em spin} observable is affected by non-Gaussian
corrections. Once these are accounted for, we find $X^\infty=4/5$ for
$d>4$ as in the Ising case~\cite{GarSolPagRit05}. For $d<4$ we provide
the first exact values of the asymptotic FDR in the presence of a
nonzero magnetization; these turn out to be highly nontrivial even to
first order in $4-d$ and $d-2$. Our results are summarized and
discussed in Sec.~\ref{sec:conclusion}. Technical details are relegated to two
appendices.

\section{Langevin dynamics and Gaussian theory}
\label{sec:Gaussian}

We consider the standard spherical model Hamiltonian
\begin{equation}
H = \frac{1}{2} \sum_{(ij)} (S_i-S_j)^2
\label{H_spherical}
\end{equation}
The sum runs over all nearest neighbour (n.n.) pairs on a
$d$-dimensional (hyper-)cubic lattice; the lattice constant is taken
as the unit length. At each of the $N$ lattice sites $\mathbf{r}_i$ there is a
real-valued spin $S_i$. The spherical constraint $\sum_i S_i^2=N$ is
imposed, which can be motivated by analogy with Ising spins $S_i=\pm
1$~\cite{BerKac52}.

The Langevin dynamics for this model can be written as
\begin{equation}
\partial_t S_i = - \frac{\partial H}{\partial S_i} + \xi_i - \frac{1}{N}
\sum_k S_k\left(-\frac{\partial H}{\partial S_k} + \xi_k\right)S_i
\label{eqn_motion}
\end{equation}
with $\xi_i$ Gaussian white noise with zero mean and covariance $\left\langle
\xi_i(t)\xi_j(t')\right\rangle = 2T\delta_{ij} \delta(t-t')$. The last term
in\eq{eqn_motion}, i.e.\ the sum over $k$, enforces the
spherical constraint at all times by removing the component of the
velocity vector
$(\partial_t S_1,\ldots,\partial_t S_N)$ along $(S_1,\ldots,S_N)$. We
use here the Stratonovic
convention for products like $S_k\xi_k$. This allows the ordinary
rules of calculus to be used when evaluating derivatives such as
$\partial_t S_i^2$. Physically it corresponds to the intuitively reasonable
scenario where the noise $\xi_i$ is regarded as a smooth random
process but with a correlation time much shorter than any other
dynamical timescale.

The prefactor of $S_i$ in the last term of\eq{eqn_motion}, being an
average of $N$ contributions, will have fluctuations of ${{\mathcal O}}(N^{-1/2})$.
Conventionally one ignores these and approximates the equation of
motion as
\begin{equation}
\partial_t S_i = - \frac{\partial H}{\partial S_i} + \xi_i - z(t) S_i
\label{eqn_motion_Gaussian}
\end{equation}
where $z(t)$ can be viewed as an effective time-dependent Lagrange multiplier
implementing the spherical constraint. This approximation works for
local quantities, but as we will see can give incorrect results when
one considers {\it e.g.}\ fluctuations of the magnetization or the energy,
which involve correlations across the entire system.  One can see
directly that\eq{eqn_motion_Gaussian} is an approximation from the
fact that it corresponds to Langevin dynamics with the effective
Hamiltonian $H+\frac{1}{2} z(t)\sum_i S_i^2$. Since the latter is
time-dependent, this dynamics does not satisfy detailed balance. It is
simple to check, on the other hand, that the original equation of
motion\eq{eqn_motion} does satisfy detailed balance and leads to the
correct equilibrium distribution $P\eql(\{S_i\})\propto \exp(-\beta
H)\delta(\sum S_i^2-N)$ where $\beta=1/T$ is the inverse temperature
as usual.

The key advantage of the approximation\eq{eqn_motion_Gaussian} is, of
course, that the spins are Gaussian random variables at all times as
long as the initial condition is of this form. Explicitly, if we
define a matrix $\Omega$ with $\Omega_{ij}=-1$ for n.n.\ sites $i,j$ and
$\Omega_{ii}=2d$, the Gaussian equation of motion is
\begin{equation}
\partial_t S_i = -\sum_j \Omega_{ij}S_j - z(t) S_i + \xi_i
\label{real_space_dotSi}
\end{equation}
We review briefly how this is solved (see e.g.~\cite{GodLuc00b} and
references therein), since these results form the
basis for all later developments. In terms of the Fourier components
$S_\mathbf{q} = \sum_i S_i \exp(-i\mathbf{q}\cdot\mathbf{r}_i)$ of the spins,
equation\eq{real_space_dotSi} reads
\begin{equation}
\partial_t S_\mathbf{q} = -(\omega_\mathbf{q}+z(t))S_\mathbf{q} + \xi_\mathbf{q}
\label{dotSq}
\end{equation}
where $\omega_\mathbf{q} = 2\sum_{a=1}^d (1-\cos q_a)$; we mostly write just
$\omega$. The Fourier mode response function can be read off as
\begin{equation}
R_\mathbf{q}(t,t_{\rm w}) = \exp\left(-\omega(t-t_{\rm w})-\int_{t_{\rm w}}^t dt'\,
z(t')\right) \equiv \sqrt{\frac{g(t_{\rm w})}{g(t)}}e^{-\omega(t-t_{\rm w})}
\label{Rq}
\end{equation}
where
\begin{equation}
g(t) = \exp\left(2\int_{0}^t dt'\, z(t')\right)
\label{g_def}
\end{equation}
In terms of this, the time-dependence of the $S_\mathbf{q}$ becomes
\begin{equation}
S_\mathbf{q}(t) = R_\mathbf{q}(t,0) S_\mathbf{q}(0) + \int_0^t dt'\,R_\mathbf{q}(t,t')\xi_\mathbf{q}(t')
\label{Sqt}
\end{equation}
The equal-time correlator $C_\mathbf{q}(t,t) = (1/N) \left\langle
S_\mathbf{q}(t)S_\mathbf{q}^*(t)\right\rangle$ follows as
\begin{eqnarray}
C_\mathbf{q}(t,t) &=& C_\mathbf{q}(0,0) R_\mathbf{q}^2(t,0) + 2T \int_0^t dt'\,R_\mathbf{q}^2(t,t') \\
&=& \frac{C_\mathbf{q}(0,0)}{g(t)}e^{-2\omega t} + 2T \int_0^t dt'\,
\frac{g(t')}{g(t)} e^{-2\omega(t-t')}
\label{Cqtt}
\end{eqnarray}
and we note for later the identity
\begin{equation}
\partial_t C_\mathbf{q}(t,t) = 2T -
\left(2\omega+\frac{g'(t)}{g(t)}\right) C_\mathbf{q}(t,t)
\label{Cqtt_deriv}
\end{equation}
The two-time correlator $C_\mathbf{q}(t,t_{\rm w}) = (1/N) \left\langle
S_\mathbf{q}(t)S_\mathbf{q}^*(t_{\rm w})\right\rangle$ can be deduced from the analogue of\eq{Sqt}
for initial time $t_{\rm w}$
\begin{equation}
S_\mathbf{q}(t) = R_\mathbf{q}(t,t_{\rm w}) S_\mathbf{q}(t_{\rm w}) + \int_{t_{\rm w}}^t dt'\,R_\mathbf{q}(t,t')\xi_\mathbf{q}(t')
\end{equation}
as
\begin{equation}
C_\mathbf{q}(t,t_{\rm w}) = R_\mathbf{q}(t,t_{\rm w}) C_\mathbf{q}(t_{\rm w},t_{\rm w})
\label{C_twotime}
\end{equation}
The position-dependent correlation and response functions
$C_{ij}(t,t_{\rm w})$ and $R_{ij}(t,t_{\rm w})$ are then just the inverse Fourier
transforms of $C_\mathbf{q}(t,t_{\rm w})$ and $R_\mathbf{q}(t,t_{\rm w})$, respectively, with
$\mathbf{q}$ conjugate to $\mathbf{r}_j-\mathbf{r}_i$.

\subsection{The function $g(t)$}

The calculations outlined above show that the Gaussian dynamics is
fully specified once the function $g(t)$ is known. The latter can be
found from the spherical constraint, which imposes $\int(dq)\,
C_\mathbf{q}(t,t)=1$. Here and below we abbreviate $(dq) \equiv
d\mathbf{q}/(2\pi)^d$, where the integrals runs over the first Brillouin zone
of the hypercubic lattice, i.e.\ $\mathbf{q}\in[-\pi,\pi]^d$. Using\eq{Cqtt}
this constraint gives an integral equation for $g(t)$:
\begin{equation}
g(t) =
\int(dq)\, C_\mathbf{q}(0,0) e^{-2\omega t} + 2T \int_0^t dt'\,  f(t-t')g(t')
\label{g_eq}
\end{equation}
where
\begin{equation}
f(t) = \int(dq)\,  e^{-2\omega t} = [e^{-4t}I_0(4t)]^d \approx
(8\pi t)^{-d/2}
\label{f_def}
\end{equation}
Here $I_0$ denotes a modified Bessel function and the final expression
gives the asymptotic behaviour for large $t$. In terms of Laplace transforms
$\hat{g}(s)=\int_0^\infty dt\, \exp(-st)g(t)$, eq.\eq{g_eq} then has the
solution
\begin{equation}
\hat{g}(s) = \frac{1}{1-2T\hat{f}(s)} \int(dq)\, \frac{C_\mathbf{q}(0,0)}{s+2\omega}
\label{hat_g_general}
\end{equation}

With the exception of Sec.~\ref{sec:magnetized}, we focus in this
paper on random initial conditions, $C_\mathbf{q}(0,0)=1$, corresponding to a
quench at time $t=0$ from equilibrium at infinite temperature. In this
case the $\mathbf{q}$-integral in the last equation is just $\hat{f}(s)$, so
that
\begin{equation}
\hat{g}(s) = \frac{\hat{f}(s)}{1-2T\hat{f}(s)}
\end{equation}
The asymptotics of the corresponding $g(t)$ are well known; see
{\it e.g.}~\cite{GodLuc00b,CorLipZan03}. For $T$ above the critical
temperature $T_{\rm c}$, which is given by
\begin{equation}
T_{\rm c}^{-1}=2 \hat{f}(0) = \int(dq)\, \frac{1}{\omega}
\label{Tc}
\end{equation}
there is a pole in $\hat{g}(s)$ at $s=2z\eql$. Here $z\eql$ is found from the
condition $2T\hat{f}(2z\eql)=1$ or
\begin{equation}
\int(dq)\, \frac{T}{z\eql+\omega} = 1
\label{zeql_cond}
\end{equation}
The presence of this pole tells us that $g(t)\sim \exp(2z\eql t)$ for
long times,
implying that the Lagrange multiplier $z(t)$ approaches $z\eql$ for
$t\to\infty$. Correspondingly, the condition\eq{zeql_cond} is just the
spherical constraint at equilibrium, bearing in mind that
$C^{\rm eq}_\mathbf{q}(t,t)=T/(z\eql+\omega)$ from\eq{dotSq}. Because
$\omega\approx q^2$
for small $q=|\mathbf{q}|$, the phase space factor in the $\mathbf{q}$-integrals is
$(dq) \sim d\omega\, \omega^{d/2-1}$ for small $q$ or $\omega$. This
shows that $T_{\rm c}$ as given by\eq{Tc} vanishes as $d\to 2$ from above;
consequently we will always restrict ourselves to dimensions $d$ above this
lower critical dimension.

At criticality ($T=T_{\rm c}$), $z\eql$ vanishes, and $g(t)$ therefore no
longer grows exponentially; instead one
finds~\cite{GodLuc00b,CorLipZan03}
\begin{equation}
g(t) \sim t^{-\kappa}, \quad \kappa =
\left\{\begin{array}{lll}
(4-d)/2 & \mbox{for} & 2<d<4\\
0 & \mbox{for} & d>4
\end{array}\right.
\label{g_asympt}
\end{equation}
It is this case, of a quench to the critical temperature, that we will
concentrate on throughout most of this paper. This is because here the
FDR has the most interesting behaviour.

We note briefly that, in principle, $\int(dq)\,$ should be written as
$(1/N)\sum_\mathbf{q}$, with the sum running over all $\mathbf{q}$ whose components
are integers in the range $-L/2\ldots -1,0,1,\ldots L/2-1$ (assuming
$L$ is even) multiplied by an overall factor $2\pi/L$; there are $N$
such $\mathbf{q}$. When considering continuous functions of $\mathbf{q}$ this sum
can be replaced by the integral $\int(dq)\,$, and this will almost always be
the case in our analysis. Exceptions are situations with a nonzero
magnetization, where the wavevector $\mathbf{q}=\mathbf{0}$ is special and has to be
treated separately. This is relevant in equilibrium below $T_{\rm c}$, which
we discuss briefly in Sec.~\ref{sec:energy_FDT_eq}, and for
non-equilibrium dynamics starting from magnetized initial states
(Sec.~\ref{sec:magnetized}).

\subsection{Long-time scaling of $C_\mathbf{q}$}

It will be useful later to have a simplified long-time expression for
$C_\mathbf{q}(t_{\rm w},t_{\rm w})$ for the case of a critical quench. At zero wavevector
one has
\begin{equation}
C_\mathbf{0}(t_{\rm w},t_{\rm w}) = \frac{1}{g(t_{\rm w})}\left(1 + 2T \int_0^{t_{\rm w}} dt'\,
g(t')\right) \approx \frac{2Tt_{\rm w}}{1-\kappa}
\label{C0}
\end{equation}
where the last approximation is based on\eq{g_asympt} and is valid for large
$t_{\rm w}$. For nonzero $\mathbf{q}$, on the other hand,
\begin{equation}
C_\mathbf{q}(t_{\rm w},t_{\rm w}) = \frac{1}{g(t_{\rm w})}\left(e^{-2\omega t_{\rm w}}
+ 2T \int_0^{t_{\rm w}} dt'\, g(t') e^{-2\omega (t_{\rm w}-t')}\right)
\approx \frac{T}{\omega}
\label{Cq}
\end{equation}
which is as expected since all nonzero Fourier modes eventually
equilibrate. The crossover between the two limits takes place when
$\omegat_{\rm w} \sim 1$, or $q\sim t_{\rm w}^{-1/2}$; physically this represents
the growth of the time-dependent correlation length as $\sim t_{\rm w}^{1/2}$.
We therefore introduce the scaling variable $w=\omega t_{\rm w}$:
\begin{equation}
\frac{C_\mathbf{q}(t_{\rm w},t_{\rm w})}{C_\mathbf{0}(t_{\rm w},t_{\rm w})}
= \frac{e^{-2w} + 2T\omega^{-1} \int_0^{w} dy\,
g(y/\omega) e^{-2(w-y)}}
{1+ 2T\omega^{-1} \int_0^{w} dy\,
g(y/\omega)}
\end{equation}
Now keep $w$ constant and let $t_{\rm w}\to\infty$, {\it i.e.}\ $\omega\to 0$.
Then $g(y/\omega) \sim (y/\omega)^{-\kappa}$ and the second terms dominate
in denumerator and nominator to give
\begin{equation}
\frac{C_\mathbf{q}(t_{\rm w},t_{\rm w})}{C_\mathbf{0}(t_{\rm w},t_{\rm w})} = (1-\kappa)\int_0^1 dy\,y^{-\kappa}
e^{-2w(1-y)}
\label{Cratio_scaling}
\end{equation}
Combining\eq{Cratio_scaling} with\eq{C0} then gives the desired
long-time scaling form
\begin{equation}
C_\mathbf{q}(t_{\rm w},t_{\rm w}) = \frac{T}{\omega}\sc{C}(\omega t_{\rm w}), \quad
\sc{C}(w) = 2w\int_0^1 dy\,y^{-\kappa} e^{-2w(1-y)}
\label{C_scaling}
\end{equation}
For $d>4$ ($\kappa=0$) this simplifies to $\sc{C}(w)=1-e^{-2w}$. As
the derivation shows, eqs.\eq{Cratio_scaling} and\eq{C_scaling} are
valid whenever $t_{\rm w}\gg 1$, even for $\omega={{\mathcal O}}(1)$. The latter
case corresponds to $w\to\infty$ and gives $\sc{C}(w)=1$, which
is indeed consistent with\eq{Cq}.

For quantities such as $C_\mathbf{q}(t_{\rm w},t_{\rm w})$ that depend only on a single
time variable, what is meant by the long-time limit is
unambiguous. For two-time quantities like $C_\mathbf{q}(t,t_{\rm w})$ we use the
following terminology: the {\em long-time} limit refers to the regime
$t\gg 1$ and $t_{\rm w}\gg 1$ but without any restriction on $t-t_{\rm w}$, which
in particular is allowed to be short, i.e.\ of ${{\mathcal O}}(1)$. The {\em aging}
regime indicates more specifically the limit $t_{\rm w}\to\infty$ at fixed
ratio $x=t/t_{\rm w}>1$, which implies that also $t-t_{\rm w}$ is large, of
${{\mathcal O}}(t_{\rm w})$. Occasionally we specialize further to the regime of
{\em well-separated} times, which corresponds to $t\ggt_{\rm w}\gg1$, i.e.\
the asymptotic behaviour of the aging limit for $x\gg 1$.

To illustrate the difference, consider which wavevectors dominate the
integral $\int(dq)\, C_\mathbf{q}(t,t_{\rm w})$. In the long-time limit at equal times
$t=t_{\rm w}$, the scaling $C_\mathbf{q}(t_{\rm w},t_{\rm w}) \sim 1/\omega$ for $\omega\gg
1/t_{\rm w}$ combined with $(dq) \sim d\omega\, \omega^{d/2-1}$ for small
$\omega$ shows that the integral is divergent at the upper end of the
frequency regime $\omega={{\mathcal O}}(t_{\rm w}^{-1})$ for all $d>2$; in other
words, it is always dominated by values of $\omega$ (and therefore
$q$) of ${{\mathcal O}}(1)$. This remains true for two-time correlations, as
long as $t-t_{\rm w}={{\mathcal O}}(1)$. In the aging limit, however, we have
$t-t_{\rm w}={{\mathcal O}}(t_{\rm w})\gg 1$ and the exponential factor from $R_\mathbf{q}$
in\eq{C_twotime} then ensures that only values of
$\omega<(t-t_{\rm w})^{-1}={{\mathcal O}}(t_{\rm w}^{-1})$ have to be considered in the
integral.

\section{Fluctuation-dissipation relations for finite-range
observables}
\label{sec:finite-range}

In this section we consider FD relations for observables that probe
correlations over a lengthscale that can be much larger than the
lattice spacing, but remains much smaller than the system size. The
latter can then be taken to infinity independently, so that the
${{\mathcal O}}(N^{-1/2})$-fluctuations of the Lagrange multiplier $z$ become irrelevant.
We begin by considering briefly spin observables, and then discuss
bond observables in some more detail.

\subsection{Spin observables}
\label{sec:spin_observables}

Since all observables that are linear in the spins can be written as
superpositions of the Fourier modes $S_\mathbf{q}$, the basic ingredient for
understanding the FD behaviour is the FDR for the
latter. Using\eq{Cqtt_deriv}, this follows after a couple of lines as
($C'\equiv \partial_{t_{\rm w}} C$)
\begin{equation}
X_\mathbf{q}(t,t_{\rm w}) = \frac{TR_\mathbf{q}(t,t_{\rm w})}{C'_\mathbf{q}(t,t_{\rm w})} =
T\left[2T - \left(\omega + \frac{g'(t_{\rm w})}{2g(t_{\rm w})}\right)
C_\mathbf{q}(t_{\rm w},t_{\rm w}) \right]^{-1}
\label{Xq}
\end{equation}
This is {\em independent} of the later time $t$, a feature that is
commonly observed in simple non-equilibrium models~\cite{CriRit03}.

The fluctuating magnetization is simply $S_\mathbf{0}/N$, so setting
$\mathbf{q}=\mathbf{0}$ in\eq{Xq} gives directly the FDR for the magnetization
\begin{equation}
X_\mathbf{0}(t,t_{\rm w}) = T\left[2T - \frac{g'(t_{\rm w})}{2g^2(t_{\rm w})}
\left(1+2T \int_0^{t_{\rm w}} dt'\, g(t')\right) \right]^{-1}
\end{equation}
As $t_{\rm w}$ increases this converges on an ${{\mathcal O}}(1)$ timescale to the
limit-FDR
\begin{equation}
X^\infty = \frac{T}{2T+(\kappa/2)[2T/(1-\kappa)]} = \frac{1-\kappa}{2-\kappa} =
\left\{\begin{array}{lll}
1/2 & (d>4)\\
1-2/d & (d<4)
\end{array}\right.
\label{X_baseline}
\end{equation}
which is identical to the value obtained from the local
magnetization~\cite{GodLuc00b} as one would expect on general
grounds. Without working out the susceptibility explicitly, it is
clear from the $t$-independence of $X_\mathbf{0}(t,t_{\rm w})$ and its fast
convergence to $X^\infty$ that the limiting FD plot is a straight
line. Both of these observations are exactly as in the Ising model in
$d=1$~\cite{MayBerGarSol03}. Simulations have shown that also in the
$d=2$ Ising case the local--global correspondence holds for spin
observables; the limiting FD plot is numerically indistinguishable
from a straight line, though renormalization group arguments suggest
that it should deviate
slightly~\cite{MayBerGarSol04,CalGam02,CalGam05}.

We should clarify that the Gaussian theory above applies directly not
to the FDR for $S_\mathbf{0}$ but to the one for $S_\mathbf{q}$ with $q\ll
t_{\rm w}^{-1/2}$ but $q\gg L^{-1}$, where $L=N^{1/d}$ is the linear system
size. The corresponding physical observable is a ``block''
magnetization, i.e.\ the average of the spins within a block of size
$\ell\sim 1/q$ much larger than the time-dependent correlation length
$\simt_{\rm w}^{1/2}$ but still small compared to the overall system
size. For $\mathbf{q}=\mathbf{0}$, i.e.\ $\ell=L$, one would in principle need to
account for the non-Gaussian fluctuations. However, it turns out that
these are negligible as long as the system is not magnetized on
average (see Sec.~\ref{sec:magnetized}), so that the above results
remain correct even for the global magnetization itself.

More generally, the FDR for any finite-range spin observable can be
expressed as a superposition of those for the Fourier modes; this can
be seen by arguments parallelling those in the $d=1$ Ising
case~\cite{MayBerGarSol03}. As there, one can then show that the
asymptotic FDR that is approached for well-separated times $t\ggt_{\rm w}\gg
1$ is dominated by the
contribution from $\mathbf{q}=\mathbf{0}$, and hence identical to $X^\infty$
calculated above~\cite{CalGam05}. At equal times, on the other hand,
equilibrated
modes with $q={{\mathcal O}}(1)$ dominate and give $X=1$. The crossover
between these two regimes takes place when $t-t_{\rm w}={{\mathcal O}}(t_{\rm w})$ and
follows (by superposition) from the corresponding crossover at
$q={{\mathcal O}}(t_{\rm w}^{-1})$ in the Fourier mode FDRs. From\eq{C_scaling}
and\eq{Xq} the latter can be expressed as
\begin{equation}
\fl X_\mathbf{q}(t,t_{\rm w}) = \sc{X}(\omegat_{\rm w}), \quad \sc{X}^{-1}(w) = 2-(2w-\kappa)
\int_0^1 dy\,y^{-\kappa} e^{-2w(1-y)}
\label{Xq_scaling}
\end{equation}
in the long-time limit, providing the expected interpolation between
$X=X^\infty=1/[2+\kappa/(1-\kappa)]$ for $w\to 0$ and $X=1$ for
$w\to\infty$.

\subsection{Bond observables}
\label{sec:bond_observables}

Next we consider bond energy observables, $\frac{1}{2}(S_i-S_j)^2$, where
$i$ and $j$ are n.n.\ sites. Since all variables are Gaussian, the
connected correlations follow by Wick's theorem. For the correlation
of bond energies one gets
\begin{eqnarray}
\fl C_{ij,kl}(t,t_{\rm w}) &=& 2\frac{1}{4}\left\langle
[S_i(t)-S_j(t)][S_k(t_{\rm w})-S_l(t_{\rm w})]\right\rangle^2
= \frac{1}{2}\left[C_{ik}-C_{il}-C_{jk}+C_{jl}\right]^2
\label{C_bond}
\end{eqnarray}
where time arguments have been
left implicit. For the local case $(i,j)=(k,l)$, this simplifies to
\begin{equation}
C_{ij,ij} = 2[1-C_{ij}]^2
\label{C_bond_local}
\end{equation}
which tends to a nonzero constant for $t=t_{\rm w}\to\infty$ since $C_{ij}$
then approaches its equilibrium value, which is $<1$.

Next we turn to the response function. In general, if one perturbs the
Hamiltonian $H$ by $-hB\delta(t-t_{\rm w})$, then the equation of motion for
$S_i$ acquires an extra term $h(\partial B/\partial S_i)\delta(t-t_{\rm w})$. So
the perturbation in $S_i$ is
\begin{equation}
\delta S_i(t) = h\sum_j R_{ij}(t,t_{\rm w})\frac{\partial B}{\partial S_j}(t_{\rm w})
\end{equation}
Thus the perturbation of an observable $A$ is
\begin{equation} \delta A(t) = h\sum_{ij} \frac{\partial A}{\partial S_i}(t)
R_{ij}(t,t_{\rm w})\frac{\partial B}{\partial S_j}(t_{\rm w}) \end{equation}
giving the response function~\cite{CalGam04}
\begin{equation}
R_{AB}(t,t_{\rm w}) = \sum_{ij} R_{ij}(t,t_{\rm w}) \left\langle \frac{\partial
A}{\partial S_i}(t)\frac{\partial B}{\partial S_j}(t_{\rm w})\right\rangle
\end{equation}
For $A=\frac{1}{2}(S_i-S_j)^2$, $B=\frac{1}{2}(S_k-S_l)^2$ this yields
\begin{equation}
R_{ij,kl}(t,t_{\rm w}) =
[R_{ik}-R_{il}-R_{jk}+R_{jl}][C_{ik}-C_{il}-C_{jk}+C_{jl}]
\end{equation}

We now analyse the scaling of correlation, response and the resulting
FDR. In terms of $C_\mathbf{q}(t,t_{\rm w})$, the bond correlation\eq{C_bond} is
\begin{equation}
\fl C_{ij,kl}(t,t_{\rm w}) = \frac{1}{2}\left[\int(dq)\, C_\mathbf{q}
\left(
 e^{i\mathbf{q}\cdot(\mathbf{r}_i-\mathbf{r}_k)}
-e^{i\mathbf{q}\cdot(\mathbf{r}_i-\mathbf{r}_l)}
-e^{i\mathbf{q}\cdot(\mathbf{r}_j-\mathbf{r}_k)}
+e^{i\mathbf{q}\cdot(\mathbf{r}_j-\mathbf{r}_l)}
\right)
\right]^2
\label{C_bond_ijkl}
\end{equation}
We can take out a factor $\exp(i\mathbf{q}\cdot\Delta\mathbf{r})$ from all the
exponentials, where $\Delta\mathbf{r} = \frac{1}{2} (\mathbf{r}_i+\mathbf{r}_j)+\frac{1}{2}(\mathbf{r}_k-\mathbf{r}_l)$ is
the distance vector between the bond midpoints. In the remaining
exponentials, $\mathbf{q}$ is multiplied by vectors with lengths of order
unity.

Now assume $t-t_{\rm w}\gg 1$. As explained above, integrals of two-time
quantities over $\mathbf{q}$ are then dominated by the small-$q$ regime,
$q^2\approx \omega<(t-t_{\rm w})^{-1}$. We can therefore Taylor expand in
$\mathbf{q}$ and get, using the equivalence of the $d$ lattice directions,
\begin{equation}
C_{ij,kl}
= \frac{1}{2}\left[d^{-1}(\mathbf{r}_i-\mathbf{r}_j)\cdot
(\mathbf{r}_k-\mathbf{r}_l) \int(dq)\, C_\mathbf{q} q^2 e^{i\mathbf{q}\cdot \Delta\mathbf{r}}\right]^2
\label{C_bond_Taylor}
\end{equation}
Similarly one finds for the response
\begin{equation}
\fl R_{ij,kl}
 = \left[d^{-1}(\mathbf{r}_i-\mathbf{r}_j)\cdot
(\mathbf{r}_k-\mathbf{r}_l)\right]^2 \int(dq)\, R_\mathbf{q} q^2 e^{i\mathbf{q}\cdot \Delta\mathbf{r}}
\int(dq)\, C_\mathbf{q} q^2 e^{i\mathbf{q}\cdot \Delta\mathbf{r}}
\end{equation}
For the {\em local} bond-bond correlation and response one sets
$\Delta\mathbf{r}=\mathbf{0}$ and has $(\mathbf{r}_i-\mathbf{r}_j)\cdot (\mathbf{r}_i-\mathbf{r}_j)=1$, which gives
for the FDR
\begin{equation}
X_{\rm bond}^{\rm loc}(t,t_{\rm w}) = \frac{\int(dq)\, TR_\mathbf{q} q^2} {\int(dq)\, C'_\mathbf{q} q^2}
= \frac{\int(dq)\, R_\mathbf{q} q^2} {\int(dq)\, X^{-1}_\mathbf{q} R_\mathbf{q} q^2 }
\label{X_local_bond}
\end{equation}
So $1/X_{\rm bond}^{\rm loc}(t,t_{\rm w})$ can be thought of as an average of
$X_\mathbf{q}^{-1}(t,t_{\rm w})$ over $\mathbf{q}$, with the weight $R_\mathbf{q}(t,t_{\rm w})
q^{-2}$. The factor $R_\mathbf{q}$ ensures that significant contributions
come only from wavevectors $\mathbf{q}$ up to length $q\sim
(t-t_{\rm w})^{-1/2}$, i.e.\ up to $\omega(t-t_{\rm w})\approx 1$. Thus, when
$t-t_{\rm w}\ll t_{\rm w}$, the result is dominated by the regime $\omega t_{\rm w}\gg
1$, where $X_\mathbf{q} = 1$. For $t-t_{\rm w}\gg t_{\rm w}$, meanwhile, one only gets
contributions from $\omegat_{\rm w}\ll 1$, where $X_\mathbf{q} =X^\infty$. So the
FDR\eq{X_local_bond} for local bond observables is a scaling function
interpolating between $1$ and $X^\infty$, with the same $X^\infty$ as
for the magnetization. Explicitly one has, using\eq{Rq} and changing
integration variable from $\mathbf{q}$ to $w=\omegat_{\rm w}$,
\begin{equation}
X_{\rm bond}^{\rm loc}(t,t_{\rm w}) = \frac{\int_0^{\infty} dw\, w^{d/2} e^{-(x-1)w}}
{\int_0^{\infty} dw\, w^{d/2} e^{-(x-1)w} \sc{X}^{-1}(w)}
\label{X_bond_local}
\end{equation}
with $x=t/t_{\rm w}$ and $\sc{X}$ the scaling form\eq{Xq_scaling} of
$X_\mathbf{q}$. To find the shape of the FD plot, recall that the equal-time
value of the local-bond correlation\eq{C_bond_local} is a constant in
the long-time limit. For $t-t_{\rm w}\gg 1$, on the other hand,
eq.\eq{C_bond_Taylor} shows that $C_{\rm bond}^{\rm loc}$ scales as
\begin{eqnarray}
C_{\rm bond}^{\rm loc}(t,t_{\rm w})
&\sim& \frac{g(t_{\rm w})}{g(t)}\left[\int d\omega\,
\omega^{d/2-1}\frac{T}{\omega}\sc{C}(\omegat_{\rm w}) e^{-\omega(t-t_{\rm w})}
\omega\right]^2 \\
&\sim& \frac{g(t_{\rm w})}{g(t)} t_{\rm w}^{-d}\left[\int dw\,
w^{d/2-1}\sc{C}(w) e^{-w(t-t_{\rm w})/t_{\rm w}}\right]^2
\end{eqnarray}
Since $\sc{C}(w)\to 1$ for $w\to\infty$, the $w$-integral would be
divergent without the exponential cutoff and scales as
$[(t-t_{\rm w})/t_{\rm w}]^{-d/2}$ for $t-t_{\rm w}\ll t_{\rm w}$, so that $C_{\rm bond}^{\rm loc}(t,t_{\rm w}) \sim
(t-t_{\rm w})^{-d}$ in this regime. The regime $t-t_{\rm w}>t_{\rm w}$ where
$X_{\rm bond}^{\rm loc}(t,t_{\rm w})\neq 1$ is therefore compressed into the region where
$C_{\rm bond}^{\rm loc}$ is of order $t_{\rm w}^{-d}$, so that the long-time limit of the
FD plot is a straight line with equilibrium slope. Qualitatively one
thus has the same behaviour as for local bond observables in the Ising
model~\cite{MayBerGarSol03}.

Next consider {\em long-range} bond observables, where we sum $(ij)$ and $(kl)$
over all bonds. The same proviso as above for the magnetization
applies here, i.e.\ by applying the Gaussian theory we are effectively
considering the bond energies averaged over a block that is large but has to
remain nonetheless small compared to the system size. One can show that the
resulting equal-time correlation again approaches a constant value for
$t_{\rm w}\to\infty$. (This follows because for large $\Delta\mathbf{r}$, one can use the
small $q$-expansion\eq{C_bond_Taylor} even for equal times. From
$C_\mathbf{q}(t_{\rm w},t_{\rm w}) \approx T/q^2$ one gets $C(\Delta\mathbf{r};t_{\rm w},t_{\rm w})\sim
|\Delta\mathbf{r}|^{2-d}$ for large $\Delta\mathbf{r}$ and so $\int(dq)\, C_\mathbf{q} q^2 e^{i\mathbf{q}\cdot\Delta\mathbf{r}}
\sim \nabla^2 |\Delta\mathbf{r}|^{2-d} \sim |\Delta\mathbf{r}|^{-d}$. The square $|\Delta\mathbf{r}|^{-2d}$
then yields a convergent sum over $\Delta\mathbf{r}$.) So we focus directly on the
regime $t-t_{\rm w}\gg 1$, where the expansion\eq{C_bond_Taylor} is again
valid. Keeping the bond $(ij)$ fixed, the scalar product
$(\mathbf{r}_i-\mathbf{r}_j)\cdot (\mathbf{r}_k-\mathbf{r}_l)$ means that only bonds $(kl)$
parallel to $(ij)$ contribute, so that the sum over $(kl)$ becomes a
sum over $\Delta\mathbf{r}$, running over all lattice vectors. (For non-parallel
bonds, $\Delta\mathbf{r}$ could also assume other values not corresponding to
lattice vectors.) The sum over $(ij)$ then just gives an overall
factor of $Nd$. Normalizing by $N$, the block bond correlation function is
\begin{equation}
C_{\rm bond}^{\rm block}(t,t_{\rm w}) = \frac{1}{2d}\frac{1}{N}\sum_{\Delta\mathbf{r}} \left[
\int(dq)\, C_\mathbf{q} q^2 e^{i\mathbf{q}\cdot \Delta\mathbf{r}}\right]^2
= \frac{1}{2d} \int(dq)\, C^2_{\mathbf{q}} q^4
\label{C_block_bond}
\end{equation}
and similar arguments give for the (normalized) response
\begin{equation}
R_{\rm bond}^{\rm block}(t,t_{\rm w}) = \frac{1}{d} \int(dq)\, R_\mathbf{q} C_\mathbf{q} q^4
\end{equation}
so that the FDR becomes
\begin{equation}
X_{\rm bond}^{\rm block}(t,t_{\rm w}) = \frac{\int(dq)\, TR_\mathbf{q} C_\mathbf{q} q^4} {\int(dq)\,
C'_\mathbf{q} C_\mathbf{q} q^4}
= \frac{\int(dq)\, R_\mathbf{q} C_\mathbf{q} q^4} {\int(dq)\, X^{-1}_\mathbf{q} R_\mathbf{q} C_\mathbf{q} q^4 }
\label{X_bond_long-range}
\end{equation}
Again, this is the inverse of a weighted average of $X^{-1}_\mathbf{q}$, now
with weight $R_\mathbf{q} C_\mathbf{q} q^4$. The same arguments as
for\eq{X_local_bond} then show that $X^{\rm block}(t,t_{\rm w})$ scales with
$x=t/t_{\rm w}$ and interpolates between $X=1$ for $x=1$ and
$X=X^\infty$ for $x\to\infty$. The value of the
correlator\eq{C_block_bond} decays from ${{\mathcal O}}(1)$ at $t=t_{\rm w}$ to
${{\mathcal O}}(t_{\rm w}^{-d/2})$ at the point $x-1\approx 1$ where aging effects
appear. While this is larger than for the local bond observables, it
still decreases to zero for $t_{\rm w}\to\infty$, so that the limiting FD
plot is again of pseudo-equilibrium form. This is different to the
case of the Ising model, where the global bond observables give
nontrivial limiting FD plots~\cite{MayBerGarSol03}.

In more detail, the scaling of the block bond correlator\eq{C_block_bond}
in the aging regime $t-t_{\rm w}\gg 1$ is
\begin{eqnarray}
C_{\rm bond}^{\rm block}(t,t_{\rm w})
&\sim& \frac{g(t_{\rm w})}{g(t)}
\int\!d\omega\, \omega^{d/2-1} \frac{T_{\rm c}^2}{\omega^2}
\sc{C}^2(\omegat_{\rm w}) e^{-2\omega(t-t_{\rm w})}\omega^2
\\
&\sim& \left(\frac{t}{t_{\rm w}}\right)^{\kappa} t_{\rm w}^{-d/2}
\int\!dw \,w^{d/2-1} \sc{C}^2(w) e^{-2(x-1)w}
\end{eqnarray}
The integral scales as $(x-1)^{-d/2}$ for $x\approx 1$, so
$C_{\rm bond}^{\rm block}(t,t_{\rm w}) \sim (t-t_{\rm w})^{-d/2}$ there. For $x\gg 1$, on the other
hand, the integral becomes $\sim x^{-(d+4)/2}$, so $C_{\rm bond}^{\rm block}(t,t_{\rm w})\sim
t_{\rm w}^{-\kappa+2} t^{\kappa-(d+4)/2}$. Explicitly, in this $t\ggt_{\rm w}$ regime,
$C_{\rm bond}^{\rm block} \sim t_{\rm w}^2 t^{-(d+4)/2}$ for $d>4$, and $C_{\rm bond}^{\rm block} \sim
t_{\rm w}^{d/2} t^{-d}$ for $d<4$. The response function scales in the same
way as $\partial_{t_{\rm w}}C_{\rm bond}^{\rm block}$, because $X$ is everywhere of order unity.

\subsection{Product observables}
\label{sec:product_observables}

Instead of the bond observables $\frac{1}{2}(S_i-S_j)^2$ we could consider
the spin products $A=S_i S_j$, $B=S_k S_l$. The correlations are then
\begin{eqnarray}
\fl C_{ij,kl}(t,t_{\rm w})
&=& \left\langle S_i(t)S_j(t)S_k(t_{\rm w})S_l(t_{\rm w})\right\rangle - \left\langle
S_i(t)S_j(t)\right\rangle \left\langle S_k(t_{\rm w})S_l(t_{\rm w})\right\rangle\\
\fl &=& C_{ik}(t,t_{\rm w})C_{jl}(t,t_{\rm w}) +
C_{il}(t,t_{\rm w})C_{jk}(t,t_{\rm w})
\label{C_prod}
\end{eqnarray}
The local equal-time correlation function $C_{ij,ij}(t_{\rm w},t_{\rm w})$ thus
approaches $2$ for $t_{\rm w}\to\infty$. The corresponding response function
is
\begin{equation}
R_{ij,kl}(t,t_{\rm w}) = R_{ik}C_{jl} + R_{il}C_{jk} + R_{jk}C_{il} + R_{jl}C_{ik}
\end{equation}
In the {\em local} case, one can replace all
functions by local ones in the aging regime -- there are no
cancellations leading to extra
factors of $q^2$ as was the case for bond observables, compare
e.g.\eq{C_bond_ijkl} and\eq{C_bond_Taylor} -- so that the FDR
\begin{equation}
X_{\rm prod}^{\rm loc}(t,t_{\rm w}) = \frac{4T_{\rm c} R_{ii}C_{ii}}{4C'_{ii}C_{ii}}
= \frac{T_{\rm c} R_{ii}}{C'_{ii}}
\label{X_prod_local}
\end{equation}
becomes identical to the one for the spin autocorrelation and
response. In particular, one again gets the same $X^\infty$.

For the global (block) case, we can write
\begin{eqnarray}
\fl C_{\rm prod}^{\rm block}(t,t_{\rm w}) &=& \sum_{(ij),(kl)}C_{ij,kl}(t,t_{\rm w})
\\
\fl &=& \frac{1}{4}\sum_{ijkl} n_{ij} n_{kl} \left[
C_{ik}(t,t_{\rm w})C_{jl}(t,t_{\rm w}) +
C_{il}(t,t_{\rm w})C_{jk}(t,t_{\rm w})\right]
\\
\fl &=& \frac{1}{2}\sum_{ijkl} n_{ij} n_{kl}
C_{ik}(t,t_{\rm w})C_{jl}(t,t_{\rm w})
= \frac{1}{2} \int(dq)\, n^2_\mathbf{q} C_\mathbf{q}^2(t,t_{\rm w})
\label{C_global_product}
\end{eqnarray}
where $n_{ij}=1$ if $i$ and $j$ are nearest neighbours and 0
otherwise, and $n_\mathbf{q}$ is its Fourier transform. For the response
one has similarly
\begin{equation}
R_{\rm prod}^{\rm block}(t,t_{\rm w}) = \int(dq)\, n^2_\mathbf{q} R_\mathbf{q}(t,t_{\rm w}) C_\mathbf{q}(t,t_{\rm w})
\end{equation}
In the aging regime, where $t-t_{\rm w}\gg 1$, the integrals are dominated
by small $q$, where $n_\mathbf{q}$ can be approximated by the constant
$n_\mathbf{0}=2d$. This cancels from the FDR, which becomes
\begin{equation}
X_{\rm prod}^{\rm block}(t,t_{\rm w}) = \frac{\int(dq)\, T_{\rm c} R_\mathbf{q} C_\mathbf{q}}{\int(dq)\, C'_\mathbf{q} C_\mathbf{q}}
\label{Xbl_large_d}
\end{equation}
This is the inverse of the average of $X_\mathbf{q}^{-1}$ with weight $R_\mathbf{q}
C_\mathbf{q}$. Again, this is a scaling function of $t/t_{\rm w}$ interpolating
between $1$ and the same $X^\infty$ as for spin observables.

The scaling of the block product correlation function\eq{C_global_product}
itself is a little more complicated than for the bond observables and
depends on dimensionality. Focussing again on $t-t_{\rm w}\gg 1$ one has
$C_{\rm prod}^{\rm block}(t,t_{\rm w}) \approx 2d^2 \int(dq)\, C_\mathbf{q}^2(t,t_{\rm w})$. The integral defines
the function $CC(t,t_{\rm w})$ discussed in Sec.~\ref{sec:noneq_large_d}
for $d>4$ and Sec.~\ref{sec:noneq_small_d} for $d<4$. In the former
case, one has from (\ref{CC_decomp}--\ref{CC_scaling}) that
$CC(t,t_{\rm w}) = CC\eql(t-t_{\rm w})\sc{CC}(t/t_{\rm w})$ where $CC\eql(t-t_{\rm w})\sim
(t-t_{\rm w})^{(4-d)/2}$ asymptotically; this equilibrium contribution
governs the behaviour of $C_{\rm prod}^{\rm block}(t,t_{\rm w})$ for $t-t_{\rm w}\ll t_{\rm w}$.  Where
aging effects appear ($t-t_{\rm w} \sim t_{\rm w}$), $C_{\rm prod}^{\rm block}\sim t_{\rm w}^{(4-d)/2}$
and so one gets a limiting pseudo-equilibrium FD plot. In the regime
of well-separated times $x\gg 1$, the scaling function $\sc{CC}(x)$ decays as
$x^{-2}$ so that $C_{\rm prod}^{\rm block}(t,t_{\rm w}) \sim t_{\rm w}^{2} t^{-d/2}$. These
scalings, though not the overall magnitude of $C_{\rm prod}^{\rm block}$, are the same as
for the energy correlation function $C_E$ in\eq{C4_final} below: both
functions are proportional to $CC(t,t_{\rm w})$ in the aging regime.

In the opposite case $d<4$, the equal-time value $CC(t,t)$ (and
therefore $C_{\rm prod}^{\rm block}(t,t)$) diverges as $t^{(4-d)/2}$,
see\eq{CC_equal_time}. The normalized correlator $CC(t,t_{\rm w})/CC(t,t)$
is a scaling function ${\mathcal{G}}(x)$ of $x=t/t_{\rm w}$, implying that the
normalized FD plot will approach a nontrivial limit form, with
asymptotic slope $X^\infty$ as shown above. Quantitatively, because
${\mathcal{G}}(x)\sim x^{-d/2}$ for $x\gg 1$, one has $C_{\rm prod}^{\rm block}(t,t_{\rm w}) \sim
t^{(4-d)/2}(t_{\rm w}/t)^{d/2} \sim t_{\rm w}^{d/2} t^{2-d}$ for $t\ggt_{\rm w}$.

\section{Correlation and response for global observables}
\label{sec:setup}

We now ask what happens if we go from block observables to
truly global ones, which reflect properties averaged over the entire
system; the total energy is an important example. We anticipate that
here non-Gaussian fluctuations are important. Indeed, the results
above show that this must be case. Otherwise we could directly extend
the Gaussian theory results from block to global observables, with no
change to correlation and response functions. The global {\em bond}
observable is just the energy. Using the spherical constraint, this
can be written as
\begin{equation}
E=\sum_{(ij)}\frac{1}{2}(S_i-S_j)^2=N-\sum_{(ij)}S_i S_j
\end{equation}
and so is identical to the global {\em spin product} observable, up to
a trivial additive constant and sign. So the global bond and product
observables must have identical correlation and response functions;
but we saw above that this requirement is not satisfied by the
Gaussian theory. Thus, non-Gaussian fluctuations are essential to get
correct results for global observables.

Physically, the origin of the distinction between block observables
and global ones is the effective infinite-range interaction induced by the
spherical constraint. In a model with short-range interactions, block
observables will show identical behaviour to global ones whenever the
block size $\ell$ is larger than any correlation length in the system,
whether or not $\ell\ll L$: the behaviour of any large subsystem is
equivalent to that of the system as a whole. In the spherical model,
the infinite-range interaction breaks this connection, and global
correlation and response functions cannot be deduced from those for
block observables.

\subsection{Non-Gaussian fluctuations}

To make progress, we need to return to the original equation of
motion\eq{eqn_motion}. This can be written as
\begin{equation}
\partial_t S_i = -\sum_{j} \Omega_{ij}S_j + \xi_i - (z(t)+N^{-1/2}\Delta z)S_i
\label{eqn_motion2}
\end{equation}
where the notation emphasizes that the fluctuating contribution to the
Lagrange multiplier is of ${{\mathcal O}}(N^{-1/2})$. The latter induces
non-Gaussian fluctuations in the $S_i$ of the same order. This shows
quantitatively why the Gaussian theory works for block observables: as
long as one considers correlations of a number of spins that is $\ll
N$, fluctuations of ${{\mathcal O}}(N^{-1/2})$ can be neglected. For global
observables, on the other hand, we require the correlations of all $N$
spins and the Gaussian approximation then becomes invalid.

To account systematically for non-Gaussian effects we represent the
spins as $S_i = s_i + N^{-1/2} r_i$, where $s_i$ gives the limiting
result for $N\to\infty$, which has purely Gaussian statistics, and
$N^{-1/2}r_i$ is a leading-order fluctuation correction which will be
non-Gaussian. Inserting this decomposition into\eq{eqn_motion2} and
collecting terms of ${{\mathcal O}}(1)$ and ${{\mathcal O}}(N^{-1/2})$ gives $\partial_t s_i =
-\Omega_{ij} s_j -z(t) s_i + \xi_i$ as expected; to lighten the
notation we use the summation convention for repeated indices from now
on. For the non-Gaussian corrections one gets the equation of motion
\begin{equation}
\partial_t r_i = -\Omega_{ij} r_j -z(t)r_i - \Delta z\, s_i
\end{equation}
with solution
\begin{equation}
r_i(t) = R_{ij}(t,0)r_j(0) - \int_0^t dt'\, R_{ij}(t,t')s_j(t')\Delta z(t')
\label{yt}
\end{equation}
The properties of $\Delta z(t')$ can now be determined from the requirement
that, due to the spherical constraint, $N^{-1}\sum_i S^2_i(t)=1$ at
all times. Inserting $S_i = s_i + N^{-1/2} r_i$ and expanding to the
leading order in $N^{-1/2}$ gives the condition
\begin{equation}
\frac{1}{N} \sum_i s_i(t)r_i(t) = -\frac{1}{2} N^{-1/2} \sum_i (s_i^2(t)-1)
\equiv -\frac{1}{2} \Delta(t)
\label{norm_condition}
\end{equation}
where the last equality defines $\Delta(t)$, a fluctuating quantity of
${{\mathcal O}}(1)$ that describes the (normalized) fluctuations of the
squared length of the Gaussian spin vector $s_i$. At $t=0$, the
condition\eq{norm_condition} is solved
to leading order by setting $r_i(0) = -\frac{1}{2} \Delta(0)s_i(0)$, since
$(1/N) \sum_i s_i^2(0)=1+{{\mathcal O}}(N^{-1/2})$. With this assignment, and
setting $a(t)=2\Delta z(t)+\Delta(0)\delta(t)$, eq.\eq{yt} reads
\begin{equation}
r_i(t) = - \frac{1}{2} \int dt'\, R_{ij}(t,t')s_j(t')a(t')
\label{yt2}
\end{equation}
We have left the integral limits unspecified here: the factor $R_{ij}$
automatically enforces $t'<t$, and we use the convention $a(t')=0$ for
$t'<0$. The spherical constraint condition\eq{norm_condition} then becomes
\begin{equation}
 \int dt'\, \frac{1}{N} s_i(t)R_{ij}(t,t')s_j(t')a(t')
= \Delta(t)
\label{norm_condition2}
\end{equation}
Now, up to fluctuations of ${{\mathcal O}}(N^{-1/2})$ which are negligible to
leading order (even if they are correlated with $a(t')$), we can
replace $(1/N) s_i(t)R_{ij}(t,t')s_j(t')$ by its average
\begin{equation}
\fl K(t,t') \equiv \frac{1}{N} \left\langles_i(t)R_{ij}(t,t')s_j(t')\right\rangle =
\frac{1}{N} R_{ij}(t,t')C_{ij}(t,t')
= \int(dq)\, R_\mathbf{q}(t,t')C_\mathbf{q}(t,t')
\label{p_def}
\end{equation}
If we then define the inverse operator, $L$, of $K$ via
\begin{equation}
\int dt'\, K(t,t')L(t',t_{\rm w}) = \delta(t-t_{\rm w})
\label{pinv_def}
\end{equation}
for $t_{\rm w}\geq 0$, then the solution to\eq{norm_condition2} is
\begin{equation}
a(t) = \int dt'\, L(t,t')\Delta(t')
\end{equation}
where for consistency we adopt the convention $\Delta(t')=0$ for
$t'<0$. With\eq{yt2} we then get an explicit expression for the
non-Gaussian ${{\mathcal O}}(N^{-1/2})$-corrections to the spins,
\begin{equation}
r_i(t) = - \frac{1}{2} \int dt' dt''\, R_{ij}(t,t')s_j(t')L(t',t'')\Delta(t'')
\label{yt3}
\end{equation}
in terms of the properties of the uncorrected Gaussian spins $s_i$.

\subsection{The functions $K$ and $L$}
\label{sec:KandL}

Before proceeding, we analyse the properties of $K$ and
$L$. From\eq{p_def}, $K(t,t')$ vanishes for $t<t'$ while its
limit value for $t\to t'^{+}$ is
$(1/N)\delta_{ij}C_{ij}(t,t)=(1/N)C_{ii}(t,t)=1$. Inserting\eq{Rq},
\eq{Cqtt_deriv} and\eq{C_twotime} into\eq{p_def}, one also finds that
the equal-time slope has the simple value $\partial_{t'} K
\left.(t,t')\right|_{t=t'^{+}}=2T$. From these properties and the
definition\eq{pinv_def} it follows that
\begin{equation}
\fl L(t,t') = \delta'(t-t') + L^{(1)}(t,t'), \qquad
L^{(1)}(t,t') = 2T\delta(t-t') - L^{(2)}(t,t')
\label{pinv_structure}
\end{equation}
where $L^{(2)}(t,t')$ vanishes for $t<t'$ and jumps to a finite value at
$t=t'^{+}$; otherwise it is smooth and, as we will later see,
positive. The structure of\eq{pinv_structure} can be easily verified
e.g.\ for the limit of equilibrium at high temperature $T$, where
$z\eql=T$ and all $\omega$ can be neglected compared to $z\eql$. One then
has $K(t,t')=\exp(-2T(t-t'))$ and the inverse\eq{pinv_def} can be
calculated by Laplace transform. Since
the Laplace transform of $K(t,t')$ is
$\hatK(s)=1/(s+2T)$ this gives $\hat{L}(s)=s+2T$, which
corresponds to\eq{pinv_structure} with $L^{(2)}\equiv 0$.

We next determine the long-time forms of $K$ and $L^{(2)}$ for
quenches to criticality.  In both cases it is useful to factor out the
equilibrium contribution. For $K$ this is, from\eq{p_def} and
using\eq{Rq} and\eq{C_twotime},
\begin{equation}
K\eql(t-t_{\rm w}) = \int(dq)\, e^{-2\omega(t-t_{\rm w})} \frac{T_{\rm c}}{\omega}
\label{p_eql}
\end{equation}
Apart from the factor of 2 in the time argument, this is just the
(critical) equilibrium spin-spin autocorrelation function. One can
also write $K(t)=2T_{\rm c}\int_t^\infty dt'\, f(t')$ from\eq{f_def} and this
shows that $K\eql(t)\sim t^{(2-d)/2}$ for large time differences. The
ratio $K(t,t_{\rm w})/K\eql(t-t_{\rm w})$ will show deviations from 1 when aging
effects appear, i.e.\ when $t-t_{\rm w}\sim t_{\rm w}$. The form of these
deviations can be worked out by using the scaling form\eq{C_scaling} of
$C_\mathbf{q}(t_{\rm w},t_{\rm w})$ and recalling that only the small $q$-regime
contributes, where $(dq)\sim d\omega\, \omega^{d/2-1}$. Changing
integration variable to $w=\omega t_{\rm w}$ gives
\begin{eqnarray}
K(t,t_{\rm w})&=&K\eql(t-t_{\rm w})\sc{K}(t/t_{\rm w})
\label{p_scaling_general}
\\
\sc{K}(x) &=& x^\kappa \frac{\int dw \, w ^{(d-4)/2} e^{-2(x-1)w
} \sc{C}(w)}
{\int dw \, w ^{(d-4)/2} e^{-2(x-1)w }}
\label{p_scaling}
\end{eqnarray}
where the first factor in\eq{p_scaling} arises from the two factors of
$[g(t_{\rm w})/g(t)]^{1/2}$ contributed by $R_\mathbf{q}(t,t_{\rm w})$ and
$C_\mathbf{q}(t,t_{\rm w})$, respectively. By construction, $\sc{K}(x)$ should
approach $1$ for $x\to 1$; indeed, in this
limit the $w$-integrals are dominated by large values of $w\sim 1/(x-1)$,
for which $\sc{C}=1$. The decay for large $x$ follows from
$\sc{C}(w)\sim w$ for small $w$ as $\sc{K}(x)\sim
x^{\kappa-1}$. Explicitly, one finds by using\eq{C_scaling} and carrying
out the $w$-integrals that
\begin{equation}
\sc{K}(x) = \frac{d-2}{2}(x-1)^{(d-2)/2}x^{\kappa} \int_0^1 dy\,
y^{-\kappa}(x-y)^{-d/2}
\label{p_scaling_fn_y_integral}
\end{equation}
For $d>4$, where $\kappa=0$, this gives
\begin{equation}
\sc{K}(x)=1-\left(\frac{x-1}{x}\right)^{(d-2)/2} \qquad (d>4)
\label{FKd_gt_4}
\end{equation}
while for $d<4$ the required indefinite integral is $[(d-2)/2]\int dy\,
y^{(d-4)/2}(x-y)^{-d/2} = -x^{-1}(x/y-1)^{(2-d)/2}$ and one gets simply
\begin{equation}
\sc{K}(x)
=x^{(2-d)/2} \qquad (d<4)
\label{FKd_lt_4}
\end{equation}

Next we determine $L^{(2)}$. Combining\eq{pinv_def}
and\eq{pinv_structure}, the defining equation is
\begin{equation}
\int dt'\, K(t,t')L^{(2)}(t',t_{\rm w}) = 2T K(t,t_{\rm w}) -
\partial_{t_{\rm w}}K(t,t_{\rm w})
\label{ptwo_def}
\end{equation}
Again it makes sense to extract the equilibrium part of $L^{(2)}$. This
is defined by
\begin{equation}
\int dt'\, K\eql(t-t')L^{(2)}\eql(t'-t_{\rm w}) = 2T K\eql(\Delta t) + K\eql'(\Delta t)
\label{ptwo_eql_condition}
\end{equation}
where $\Delta t=t-t_{\rm w}$. Solving by Laplace transform gives
\begin{equation}
\hatL^{(2)}\eql(s) =
\frac{2T\hatK\eql(s)+(s\hatK\eql(s)-1)}{\hatK\eql(s)} =
s + 2T - \frac{1}{\hatK\eql(s)}
\label{ptwo_eql}
\end{equation}
where from\eq{p_eql}, at criticality,
\begin{equation}
\hatK\eql(s) = T_{\rm c} \int(dq)\, \frac{1}{\omega(s+2\omega)}
\label{p_eql_LT}
\end{equation}
The leading small-$s$ behaviour of this is
$\hatK\eql(0)-\hatK\eql(s) \sim s^{(d-4)/2}$ for $d>4$ (plus, for
$d>6$, additional analytic terms of integer order in $s$ which are
irrelevant for us). For $d<4$, on the other hand, $\hatK\eql(s) \sim
s^{(d-4)/2}$ is divergent for $s\to 0$. Inserting these scalings
into\eq{ptwo_eql} and inverting the Laplace transform gives for the
asymptotic behaviour of $L^{(2)}\eql$
\begin{equation}
L^{(2)}\eql(t) \sim \left\{ \begin{array}{ll}
t^{(2-d)/2} & \mbox{($d>4$)} \\
t^{(d-6)/2} & \mbox{($d<4$)}
\end{array}\right.
\label{ptwo_eql_asymptotics}
\end{equation}
It will be important below that, for $d>4$, $K\eql(t)$ and
$L^{(2)}\eql(t)$ both decay asymptotically as $t^{(2-d)/2}$. The ratio
between them can be worked out from\eq{ptwo_eql}, by expanding
for small $s$ as $1/\hatK\eql(s) \approx 1/[\hatK\eql(0)-cs^{(d-4)/2}] =
1/\hatK\eql(0) + cs^{(d-4)/2}/\hatK^2\eql(0)$ where $c$ is some
constant; comparing with $\hatK\eql(s) \approx
\hatK\eql(0)-cs^{(d-4)/2}$ gives
\begin{equation}
L^{(2)}\eql(t)=K\eql(t)/\hatK\eql^2(0)
\label{ptwo_p_link}
\end{equation}
for large time differences $t$.

The integral of $L^{(2)}\eql(t)$ over all times follows
from\eq{ptwo_eql} as
\begin{equation}
\fl \hatL^{(2)}\eql(0) = \int_0^\infty dt\,L^{(2)}\eql(t) = \left\{
\begin{array}{ll}
2T_{\rm c} - \hatK\eql^{-1}(0) = 2T_{\rm c}[1-1/\int(dq)\, (T_{\rm c}/\omega)^{2}]
& \mbox{($d>4$)}\\
2T_{\rm c} & \mbox{($d<4$)}
\end{array}
\right.
\label{L2hat0}
\end{equation}
Using the fact that $\int(dq)\, (T_{\rm c}/\omega)=1$, one has $\int(dq)\,
(T_{\rm c}/\omega)^2>1$ so that $\hatL^{(2)}\eql(0)$ is positive independently
of $d$. This is consistent with the intuition that, with the sign as
chosen in\eq{pinv_structure}, the function $L^{(2)}$ is positive.

With the equilibrium part of $L^{(2)}$ determined we make a long-time
scaling ansatz for $L^{(2)}$,
\begin{equation}
L^{(2)}(t,t_{\rm w}) = L^{(2)}\eql(t-t_{\rm w})\sc{L}(t/t_{\rm w})
\label{L_scaling}
\end{equation}
so that\eq{ptwo_def} becomes
\begin{eqnarray}
\fl\lefteqn{\int dt'\,
K\eql(t-t')L^{(2)}\eql(t'-t_{\rm w}) \sc{K}(t/t')
\sc{L}(t'/t_{\rm w}) = }
\nonumber\\
&=&2T_{\rm c} K\eql(\Delta t)\sc{K}(x)+
K\eql'(\Delta t)\sc{K}(x)+
\frac{t}{t_{\rm w}^2}K\eql(\Delta t)\sc{K}'(x)
\label{ptwo_scaling_cond}
\end{eqnarray}
where $\Delta t=t-t_{\rm w}$ and $x=t/t_{\rm w}$ as before. We now take the aging limit
of large
$t_{\rm w}$ with $\Delta t={{\mathcal O}}(t_{\rm w})$ to determine $\sc{L}$. The second and
third terms on the r.h.s.\ are then smaller by factors of order
$1/t_{\rm w}$ than the first, and can be neglected to leading order. The
second term on the r.h.s.\ of\eq{ptwo_eql_condition} is likewise
subdominant, and this can be used to rewrite the dominant first term
on the r.h.s.\ of\eq{ptwo_scaling_cond}, giving
\begin{eqnarray}
\fl \sc{K}\left(x\right)
&=& \frac{\int dt'\,K\eql(t-t')L^{(2)}\eql(t'-t_{\rm w})
          \sc{K}(t/t')\sc{L}(t'/t_{\rm w})}
         {\int dt'\, K\eql(t-t')L^{(2)}\eql(t'-t_{\rm w})}
\label{unscaled_ptwo}
\end{eqnarray}
We consider first $d>4$. Then both the functions $K\eql(\Delta t)$ and
$L^{(2)}\eql(\Delta t)$ have finite integrals $\hatK\eql(0)$ and
$\hatL^{(2)}\eql(0)$, respectively, over $\Delta t=0\ldots\infty$. In the
aging limit, the factors
$K\eql(t-t')$ and $L^{(2)}\eql(t'-t_{\rm w})$ therefore act to concentrate
the mass of the integrals appearing in\eq{unscaled_ptwo} around
$t'=t_{\rm w}$ and $t'=t$. This can be seen more formally by changing to
$y=t'/t_{\rm w}$ as the integration variable and taking $t_{\rm w}$
large. Then the factors $K\eql(t_{\rm w}(x-y))$ and $L^{(2)}\eql(t_{\rm w}(y-1))$
produce singularities $\sim(x-y)^{(2-d)/2}$ for $y\to x$ and $\sim
(y-1)^{(2-d)/2}$ for $y\to 1$, respectively, and because these are
non-integrable they dominate the integral for $t_{\rm w}\to\infty$. All
other factors in the integrals are slowly varying near the relevant
endpoints and can be replaced by their value there. In the aging limit we
can therefore write\eq{unscaled_ptwo} as
\begin{equation}
\fl \sc{K}(x)
= \frac{\hatK\eql(0)L^{(2)}\eql(\Delta t) \sc{K}(1)\sc{L}(x)+
K\eql(\Delta t)\hatL^{(2)}\eql(0) \sc{K}(x)\sc{L}(1)}
{\hatK\eql(0)L^{(2)}\eql(\Delta t)+K\eql(\Delta t)\hatL^{(2)}\eql(0)}
\end{equation}
Eq.\eq{ptwo_p_link} tells us that the $\Delta t$-dependent factors cancel,
giving together with\eq{L2hat0} and $\sc{K}(1)=\sc{L}(1)=1$
\begin{equation}
\fl \sc{K}(x)
= \frac{\hatK\eql^{-1}(0) \sc{L}(x)+\hatL^{(2)}\eql(0) \sc{K}(x)}
{\hatK\eql^{-1}(0)+\hatL^{(2)}\eql(0)}
= \frac{\sc{L}(x)+[2T_{\rm c}\hatK\eql(0)-1]\sc{K}(x)}
{2T_{\rm c}\hatK\eql(0)}
\label{ptwo_simple_scaling}
\end{equation}
In $d>4$, where $\hatK\eql(0)$ is finite, we therefore have the
simple result that the scaling functions for $K$ and $L^{(2)}$ are
identical,
\begin{equation}
\sc{L}(x) = \sc{K}(x)
\label{same_sc_fn}
\end{equation}
But in the limit $d\to 4$ from above, $\hatK\eql(0)$ diverges
and\eq{ptwo_simple_scaling} gives no information about
$\sc{L}$. For $d<4$ a different approach is therefore needed
to determine $\sc{L}$. One subtracts from\eq{ptwo_scaling_cond}
the first and second terms on its r.h.s., using\eq{ptwo_eql_condition} to
rewrite them as an integral and changing integration variable from
$t'$ to $y=t'/t_{\rm w}$. This gives
\begin{eqnarray}
\fl\lefteqn{
t_{\rm w}\int_1^x dy\,
K\eql(t_{\rm w}(x-y))L^{(2)}\eql(t_{\rm w}(y-1))\left[\sc{K}\left({x}/{y}\right)
\sc{L}\left(y\right)-\sc{K}(x)\right]
= }\nonumber\\
&=&\frac{t}{t_{\rm w}^2}K\eql(\Delta t)\sc{K}'(x)
\label{L_condition_almost}
\end{eqnarray}
For $y\to x$, $K\eql(t_{\rm w}(x-y))$ contributes a singularity $\sim
(x-y)^{(2-d)/2}$ which is integrable in $d<4$. For $y\to 1$, the terms
in square brackets vanish as $\sim y-1$ since $\sc{L}(y)$ is
smooth at $y=1$ as we will see below, in the sense that
$\sc{L}'(1)$ is finite. These terms combine with the $\sim
(y-1)^{(d-6)/2}$ from $L^{(2)}\eql$ to give an integrable $\sim
(y-1)^{(d-4)/2}$. The contributions from the short time behaviour of
$K\eql$ and $L^{(2)}\eql$ are therefore unimportant in the aging
limit and we can replace these functions by their power-law
asymptotes. Up to overall $d$-dependent numerical factors the
condition\eq{L_condition_almost} then becomes
\begin{eqnarray}
\fl\lefteqn{
t_{\rm w}^{-1}\int_1^x dy\,
(x-y)^{(2-d)/2}(y-1)^{(d-6)/2}\left[\sc{K}\left({x}/{y}\right)
\sc{L}\left(y\right)-\sc{K}(x)\right]
= }\nonumber\\
&=&\frac{t}{t_{\rm w}^2}\Delta t^{(2-d)/2}\sc{K}'(x)
\end{eqnarray}
In the aging limit $\Delta t$ scales as $t_{\rm w}$, and so the l.h.s.\ of this
equation ($\simt_{\rm w}^{-1}$) becomes large compared to the r.h.s.\ ($\sim
t_{\rm w}^{-d/2})$ unless the $y$-integral vanishes. The required
condition for $\sc{L}$ is therefore
\begin{equation}
\fl
\int_1^x dy\,
(x-y)^{(2-d)/2}(y-1)^{(d-6)/2}\left[\sc{K}\left({x}/{y}\right)
\sc{L}\left(y\right)-\sc{K}(x)\right] = 0
\label{pinv_cond}
\end{equation}
This is in principle an integral equation for
$\sc{L}$. Fortunately, however, the solution is the naive
extension of\eq{same_sc_fn} to $d<4$: with $\sc{K}(x)=x^{(2-d)/2}$
from\eq{FKd_lt_4} one sees that for
$\sc{L}(x)=\sc{K}(x)=x^{(2-d)/2}$ the square bracket
in\eq{pinv_cond} vanishes identically. The identity\eq{same_sc_fn}
therefore holds both for $d<4$ and for $d>4$.

In summary, we have determined long-time scaling forms for $K$ and $L$ for
quenches to criticality. For $K$, the result is\eq{p_scaling_general}
with\eq{p_eql} and (\ref{FKd_gt_4},\ref{FKd_lt_4}); for $L$, we
have\eq{L_scaling} with\eqq{ptwo_eql}{ptwo_eql_asymptotics}
and\eq{same_sc_fn}. Combined with\eq{yt3}, this fully determines the
leading non-Gaussian corrections to the spherical model dynamics (at
long times, and after a quench to criticality from a random initial state).

\section{General expressions for energy correlation and response}
\label{sec:energy_general}

In this section we derive general expressions for the two-time correlation and
response functions of the energy, taking into account non-Gaussian
fluctuations. The results will be valid for arbitrary quenches since
we will leave $K$ and $L$ unspecified.

\subsection{Energy correlation function}
\label{sec:energy_corr}

We can write the energy\eq{H_spherical} as $H=\frac{1}{2} S_i\Omega_{ij}
S_j$. Inserting $S_i = s_i + N^{-1/2} r_i$, the energy correlation
function (normalized by $N$) is to leading order
\begin{equation}
\fl C_E(t,t_{\rm w})=\frac{1}{4N}\left\langle \left.\left(
s_i \Omega_{ij} s_j + 2N^{-1/2} s_i \Omega_{ij} r_j
\right)\right|_{t}
\left.\left(
s_k \Omega_{kl} s_l + 2N^{-1/2} s_k \Omega_{kl} r_l
\right)\right|_{t_{\rm w}} \right\rangle'
\label{CE_basic}
\end{equation}
Using\eq{yt3}, all quantities involved can be expressed in terms of
the Gaussian variables $s_i$ so that the average can be performed using
Wick's theorem, i.e.\ by taking products of all possible pairings. We use
the prime on the average to indicate the connected correlation
function. This just means that in the Wick expansion all terms not
containing any pairings of a variable at $t$ with one at $t_{\rm w}$ have to
be discarded, since these terms give the disconnected contribution $\left\langle
\left. (\ldots)\right|_{t}\right\rangle \left\langle\left. (\ldots)\right|_{t_{\rm w}}
\right\rangle$. Multiplying out\eq{CE_basic} one obtains four
contributions. The first one is
\begin{equation}
\fl 4C_E^{(1)} = \frac{1}{N}\Omega_{ij}\Omega_{kl}\left\langle
s_i(t)s_j(t)s_k(t_{\rm w})s_l(t_{\rm w})\right\rangle' =
\frac{2}{N}\Omega_{ij}\Omega_{kl}C_{jk}(t,t_{\rm w})C_{il}(t,t_{\rm w})
\label{C1_first}
\end{equation}
To eliminate one of the factors of $\Omega$, note from\eq{Rq} that
$(\partial_t+z(t)) R_\mathbf{q}(t,t_{\rm w}) = -\omega R_\mathbf{q}(t,t_{\rm w})$ for
$t>t_{\rm w}$. In real space, this reads
\begin{equation}
\left(\partial_t+z(t)\right)R_{ik}(t,t_{\rm w}) =
-\Omega_{ij}R_{jk}(t,t_{\rm w}) =
-R_{ij}(t,t_{\rm w})\Omega_{jk}
\label{Om_elimination}
\end{equation}
and from\eq{C_twotime} an exactly analogous relation holds for
$C_{ij}(t,t_{\rm w})$. Thus
\begin{eqnarray}
\fl 4C_E^{(1)}
&=&
- \frac{2}{N}[(\partial_t+z(t))C_{ik}(t,t_{\rm w})]\Omega_{kl}C_{il}(t,t_{\rm w})
\\
\fl &=&
- \frac{1}{N}(\partial_t+2z(t))C_{ik}(t,t_{\rm w})\Omega_{kl}C_{il}(t,t_{\rm w}) \equiv
- (\partial_t+2z(t))C\Om C(t,t_{\rm w})
\label{C1}
\end{eqnarray}
The last equality defines $C\Om C$, which is just the normalized trace
of the product of the matrices $C_{ik}(t,t_{\rm w})$, $\Omega_{kl}$ and
$C_{il}=C_{li}$; in Fourier space, $C\Om C(t,t_{\rm w}) = \int(dq)\, \omega
C_\mathbf{q}^2(t,t_{\rm w})$.

\begin{figure}
\centerline{\includegraphics[width=7.0cm,clip=true]{C2.eps}}
\caption{Illustration of Wick pairings for $C_E^{(2)}$. (a) The six
Gaussian spins that need to be paired in\eq{C2_raw} are indicated by
circles with site labels. Spins arising from the expansion of $H(t)$
and $H(t_{\rm w})$ are to the left and right of the vertical dashed line,
respectively; any pairing contributing to the connected correlator
must have links across this line. Dotted lines connect spins that
are already coupled in space: $s_i$ and $s_m$ via the factor
$\Omega_{ij}R_{jm}$, and $s_k$ and $s_l$ via $\Omega_{kl}$. The vertical
lines between $s_n$ and $s_n$ indicate that pairings which couple these
spins are not allowed. (b) The solid lines show the only Wick pairing
that contributes to
leading order in $1/N$: it gives two independent groups of spins. (c)
only gives one group and so is subleading. (d) has two groups but is
excluded from the connected correlator because there are no pairs
across the dashed line.
\label{fig:C2}
}
\end{figure}
The second contribution to $C_E$ is
\begin{eqnarray}
4C_E^{(2)} &=& \frac{2}{N^{3/2}}\Omega_{ij}\Omega_{kl}\left\langle
s_i(t)r_j(t)s_k(t_{\rm w})s_l(t_{\rm w})\right\rangle' \\
&=&
-\frac{1}{N} \int dt' dt''\, \Omega_{ij}\Omega_{kl} R_{jm}(t,t')L(t',t'')
\nonumber\\
& &\times \left\langle s_i(t)s_m(t')N^{-1/2}\Delta(t'') s_k(t_{\rm w})s_l(t_{\rm w}) \right\rangle'
\end{eqnarray}
where we have inserted\eq{yt3}. Writing $N^{-1/2}\Delta(t'')$ explicitly this
becomes
\begin{eqnarray}
4C_E^{(2)} &=&
-\frac{1}{N^2} \int dt' dt''\,
\Omega_{ij}\Omega_{kl}R_{jm}(t,t')L(t',t'') \nonumber\\
& &\times
\left\langle s_i(t)s_m(t')\left[s_n^2(t'')-\left\langle
s_n^2(t'')\right\rangle\right] s_k(t_{\rm w})s_l(t_{\rm w}) \right\rangle'
\label{C2_raw}
\end{eqnarray}
We now need to perform the Wick expansion of the average. The
subtraction $s_n^2-\left\langle s_n^2\right\rangle$ means that all terms which would
pair $s_n$ with $s_n$ are excluded; $s_k$ and $s_l$ also cannot be
paired because we are considering the connected correlation. We can
reduce the number of pairings further by considering that we need to
get an overall result of ${{\mathcal O}}(1)$. The index $j$ does not need to
be considered further: after summing over $j$, $\Omega_{ij}R_{jm}$ is
some translationally invariant function of the distance vector
between spins $s_i$ and $s_m$. If the remaining indices $i,k,l,m,n$
where unrestricted, then together with the $1/N^2$ prefactor we would
maximally get an ${{\mathcal O}}(N^3)$ result. Each of the factors $\Omega_{ij}R_{jm}$
and $\Omega_{kl}$ couples two indices and so reduces the order of
the result by $1/N$. Having already got two such couplings outside
the average, we can only ``afford'' one extra coupling from the Wick
pairings to get a contribution of ${{\mathcal O}}(1)$. After some reflection
one sees that this only leaves the pairing $[im][kn][ln]$: $[im]$
introduces no further coupling beyond $\Omega_{ij}R_{jm}$, and
$[kn][ln]$ gives only one further coupling beyond
$\Omega_{kl}$. Alternatively, we can think of this pairing as having the
indices $i,m$ and $k,l,n,n$ in two independent groups; each group
gives a factor of $N$ and this just cancels the $1/N^2$ prefactor. All
other allowed Wick pairings give smaller terms, as illustrated
graphically in Fig.~\ref{fig:C2}a,b. For example,
$[ik][mn][ln]$ together with $\Omega_{ij}R_{jm}$ and $\Omega_{kl}$
couples {\em all} indices into a single group and thus gives a term of
only ${{\mathcal O}}(1/N)$ (Fig.~\ref{fig:C2}c). The pairing
$[in][mn][kl]$ would give two
independent groups and thus an ${{\mathcal O}}(1)$ term, but is excluded
because $k$ and $l$ cannot be paired in the connected correlator
(Fig.~\ref{fig:C2}d). Bearing in mind that our
dominant pairing has a symmetry factor of 2 because the $s_n$'s in
$[kn][ln]$ can be swapped we have thus finally
\begin{eqnarray}
\fl 4C_E^{(2)} &=&
-2 \int dt' dt''\,
\frac{1}{N} \Omega_{ij} R_{jm}(t,t') C_{im}(t,t') L(t',t'')
\nonumber\\
\fl & &\times
\frac{1}{N}\Omega_{kl}C_{nk}(t'',t_{\rm w})C_{nl}(t'',t_{\rm w}) \\
\fl &=&
\int dt' dt''\,
[(\partial_t+2z(t))K(t,t')-\delta(t-t')]
 L(t',t'') C\Om C(t'',t_{\rm w})
\end{eqnarray}
In going to the last line we have exploited\eq{Om_elimination} to
eliminate $\Omega_{ij}$. Since $t'$ is an integration variable, we have
also been careful here to subtract off with the $\delta(t-t')$ the
spurious contribution which the $\partial_t$ applied to the step
discontinuity in $K(t,t')$ would otherwise give. The $t'$-integration
can now be carried out using\eq{pinv_def} and we get
\begin{eqnarray}
4C_E^{(2)}
&=&
(\partial_t+2z(t)) C\Om C(t,t_{\rm w})-\int dt' L(t,t') C\Om C(t',t_{\rm w})
\\
&=&
2z(t) C\Om C(t,t_{\rm w})-\int dt' L^{(1)}(t,t') C\Om C(t',t_{\rm w})
\label{C2}
\end{eqnarray}
Note that we could use the same trick here as
in\eqq{Om_elimination}{C1} to write
$C\Om C(t',t_{\rm w})=-(\frac{1}{2}\partial_{t'}+z(t'))CC(t',t_{\rm w})$, with
$CC(t,t_{\rm w})$ defined in the obvious manner as $CC = N^{-1}
C_{ik}C_{ki} = \int(dq)\, \omega C_\mathbf{q}^2$.  However, this reduction from
$C\Om C$ to $CC$ applies only for $t'>t_{\rm w}$, while for $t'<t_{\rm w}$ one
would need to take a time derivative w.r.t.\ $t_{\rm w}$ instead of $t'$.
This case distinction would make evaluation of\eq{C2} awkward, so we
retain $C\Om C$ here and below.

The third contribution to $C_E$ is obtained by simply swapping the
roles of $t$ and $t_{\rm w}$ in\eq{C2} and remembering that $C\Om C$ is
symmetric in its time arguments,
\begin{eqnarray}
4C_E^{(3)} &=&
\frac{2}{N^{3/2}}\Omega_{ij}\Omega_{kl}\left\langle s_i(t)s_j(t)s_k(t_{\rm w})r_l(t_{\rm w})\right\rangle' \\
&=&
2z(t_{\rm w}) C\Om C(t,t_{\rm w})-\int dt' L^{(1)}(t_{\rm w},t') C\Om C(t,t')
\label{C3}
\end{eqnarray}

\begin{figure}
\centerline{\includegraphics[width=9.0cm,clip=true]{C4.eps}}
\caption{Wick pairings for $C_E^{(4)}$. (a) represents the constraints
for the possible pairings in\eq{C4_raw}. (b) is the only pairing that
contributes to leading order, forming three independent groups of
spins.
\label{fig:C4}
}
\end{figure}
The fourth and last contribution to $C_E$ is, using again\eq{yt3} and
writing out $\Delta(t'')$ and $\Delta(t_{\rm w}'')$
\begin{eqnarray}
\fl 4C_E^{(4)} &=& \frac{4}{N^2}\Omega_{ij}\Omega_{kl}\left\langle
s_i(t)r_j(t)s_k(t_{\rm w})r_l(t_{\rm w})\right\rangle' \\
\fl &=&
\frac{1}{N^3} \int dt' dt'' dt_{\rm w}' dt_{\rm w}''\,
\Omega_{ij}\Omega_{kl}
R_{jm}(t,t')L(t',t'')R_{ln}(t_{\rm w},t_{\rm w}')L(t_{\rm w}',t_{\rm w}'')
\nonumber\\
\fl & & \times \left\langle
s_i(t)  s_m(t')  s_o^2(t'')
s_k(t_{\rm w})s_n(t_{\rm w}')s_p^2(t_{\rm w}'')
\right\rangle'
\label{C4_raw}
\end{eqnarray}
where it is understood that, because of subtractions which we have not
written explicitly, pairings of $s_o$ with itself and of $s_p$ with
itself are to be excluded. The only pairing that gives an overall
${{\mathcal O}}(1)$ contribution turns out to be $[im][kn][op][op]$ and gives
the required 3 independent index groups ($i,m$; $k,n$; $o,o,p,p$) to
cancel the $1/N^3$ prefactor; see Fig.~\ref{fig:C4}. With the symmetry
factor 2 for the possible swap of the $[op][op]$ pairings one gets
\begin{eqnarray}
\fl 4C_E^{(4)}
&=&
2 \int dt' dt'' dt_{\rm w}' dt_{\rm w}''\,
\frac{1}{N}\Omega_{ij}R_{jm}(t,t')C_{im}(t,t')
\nonumber\\
\fl & & \times
\frac{1}{N}\Omega_{kl}R_{ln}(t_{\rm w},t_{\rm w}')C_{kn}(t_{\rm w},t_{\rm w}')
L(t',t'')L(t_{\rm w}',t_{\rm w}'')\frac{1}{N} C_{op}^2(t'',t_{\rm w}'')\\
\fl &=&
\frac{1}{2} \int dt' dt'' dt_{\rm w}' dt_{\rm w}''\,
[(\partial_t+2z(t))K(t,t')-\delta(t-t')]L(t',t'')
\nonumber\\
\fl & & \times
[(\partial_{t_{\rm w}}+2z(t_{\rm w}))K(t_{\rm w},t_{\rm w}')-\delta(t_{\rm w}-t_{\rm w}')]
L(t_{\rm w}',t_{\rm w}'') CC(t'',t_{\rm w}'')
\end{eqnarray}
Using\eq{pinv_def} one can carry out
two of the time integrations to get
\begin{eqnarray}
4C_E^{(4)}
&=&
\frac{1}{2} \int dt' dt_{\rm w}' \,
[2z(t)\delta(t-t')-L^{(1)}(t,t')]
\nonumber\\
& & \times
[2z(t_{\rm w})\delta(t_{\rm w}-t_{\rm w}')-L^{(1)}(t_{\rm w},t_{\rm w}')]
CC(t',t_{\rm w}')
\label{C4}
\end{eqnarray}

\subsection{Energy response function}

To find the energy response, consider perturbing the Hamiltonian by a
term $-h\delta(t-t_{\rm w})H=-(h/2)\delta(t-t_{\rm w})S_i\Omega_{ij}S_j$, where $h$
is the field conjugate to the energy. The equation of motion in the
presence of the perturbation is therefore
\begin{equation}
\partial_t S_i = - \Omega_{ij}S_j - (z(t)+h\Delta z(t))S_i +h\delta(t-t_{\rm w})\Omega_{ij}S_j
\end{equation}
where $h \Delta z$ is now the change in the Lagrange multiplier induced by
the perturbation. The fluctuating component of $z$ of ${{\mathcal O}}(N^{-1/2})$
is in principle still present, but negligible compared to $h\Delta z$ for
field strengths that are ${{\mathcal O}}(N^0)$. Inserting a corresponding
expansion of the spins, $S_i=s_i+hr_i$, gives for $s_i$ the unperturbed
equation of motion and for the perturbed component $r_i$
\begin{equation}
\partial_tr_i = -\Omega_{ij}r_j - z(t)r_i -\Delta z(t)s_i +\delta(t-t_{\rm w})\Omega_{ij}s_j
\label{dot_ri}
\end{equation}
The solution of this is
\begin{equation}
r_i(t) = R_{ik}(t,t_{\rm w})\Omega_{kl}s_l(t_{\rm w}) - \int_{t_{\rm w}}^t dt'\,
R_{ik}(t,t') \Delta z(t') s_k(t')
\label{yt_resp}
\end{equation}
One now needs to determine $\Delta z$. This can be done by noting that the
normalized length of $S_i$ is
\begin{equation}
\frac{1}{N}\sum_i S_i^2 = \frac{1}{N}\sum_i s_i^2 + 2h \frac{1}{N}
\sum_i s_i r_i + {{\mathcal O}}(h^2)
\end{equation}
The change to first order in $h$ must vanish to preserve the spherical
constraint, giving the condition $(1/N)\sum_i \left\langle s_i r_i\right\rangle=0$
or, using\eq{yt_resp},
\begin{equation}
\frac{1}{N} R_{ik}(t,t_{\rm w})\Omega_{kl}C_{il}(t,t_{\rm w}) =
\int_{t_{\rm w}}^t dt'\, \frac{1}{N} R_{ik}(t,t')C_{ik}(t,t')\Delta z(t')
\end{equation}
In the integrand one recognizes the definition\eq{p_def} of $K$, so
that one can write the solution of this as
\begin{equation}
\Delta z(t) = \int dt'\, L(t,t') R\Om C(t',t_{\rm w})
\end{equation}
with obvious notation for $R\Om C$.

Now we can find the change in the energy, $1/(2N) \left\langle S_i\Omega_{ij}
S_j\right\rangle$, which is given by $(h/N) \left\langle r_i \Omega_{ij}s_j\right\rangle$ to linear order
in $h$. Dividing by $h$ and using\eq{yt_resp} then gives the energy
response function
\begin{eqnarray}
R_E(t,t_{\rm w}) &=& \frac{1}{N}\left\langle  r_i \Omega_{ij}s_j\right\rangle\\
&=&
\frac{1}{N} R_{ik}(t,t_{\rm w})\Omega_{kl}\Omega_{ij}C_{jl}(t,t_{\rm w})
\nonumber\\
& &
{} - {}\int dt' dt''\,
R\Om C(t,t')L(t',t'') R\Om C(t'',t_{\rm w})
\label{RE_general}
\end{eqnarray}
One can eliminate one of the time integrals by
using\eq{Om_elimination}, being careful to remove the unwanted
contribution from differentiating the step discontinuity in
$R\Om C(t,t_{\rm w})$. Using
also\eq{pinv_def} and\eq{pinv_structure} then gives
\begin{eqnarray}
\fl 2R_E
&=&
(-\partial_t-2z(t))R\Om C(t,t_{\rm w}) + \delta(t-t_{\rm w})R\Om C(t_{\rm w}^{+},t_{\rm w})
\nonumber\\
\fl & &{} - {}\int dt' dt''\,
[(-\partial_t-2z(t))K(t,t')+\delta(t-t')]L(t',t'') R\Om C(t'',t_{\rm w})
\\
\fl &=&
- \int dt' L(t,t') R\Om C(t',t_{\rm w})
+ \delta(t-t_{\rm w})R\Om C(t_{\rm w}^{+},t_{\rm w})
\label{RE}
\end{eqnarray}

\subsection{Equilibrium}
\label{sec:energy_FDT_eq}

Above, we derived general expressions for the energy two-time
correlation and response, in terms of the known correlation and
response functions for the Gaussian spins which in turn determine $K$ and
$L$. Before looking at non-equilibrium, we consider briefly the
equilibrium situation; even here the results for the dynamics are new
as far as we are aware.

For the response function one uses that at equilibrium $R\Omega^aC(t) =
\int(dq)\, e^{-2(z\eql+\omega)t}\omega^a [T/(z\eql+\omega)]$ for
$a=1,2$. Here we have retained a possible nonzero equilibrium value $z\eql$ of
the Lagrange multiplier, to include the
case of equilibrium at $T\neq T_{\rm c}$. Inserting into\eq{RE_general} and
taking LTs gives
\begin{eqnarray}
\hat R_{E}^{\rm eq}(s) &=& \int(dq)\,
\frac{T\omega^2}{(z\eql+\omega)[s+2(z\eql+\omega)]} \nonumber\\
& & {} - {} \frac{1}{\hatK\eql(s)} \left(\int(dq)\,
\frac{T\omega}{(z\eql+\omega)[s+2(z\eql+\omega)]}\right)^2
\end{eqnarray}
where $\hatK\eql$ is generalized from\eq{p_eql_LT} to
\begin{equation}
\hatK\eql(s) = T \int(dq)\, \frac{1}{(z\eql+\omega)[s+2(z\eql+\omega)]}
\label{phat_eql}
\end{equation}
Using the spherical constraint condition
$T\int(dq)\, (z\eql+\omega)^{-1} = 1$ one then shows, after a few lines of
algebra, that
\begin{equation}
\hat R_{E}^{\rm eq}(s)=\frac{1}{4}\left(s + 2T -
\frac{1}{\hatK\eql(s)}\right) = \frac{1}{4}\hatL^{(2)}\eql(s)
\label{R_Eql}
\end{equation}
Remarkably, therefore, the equilibrium energy response function
$R_{E}^{\rm eq}(t)$ is directly proportional to the inverse kernel
$\hatL^{(2)}\eql(t)$. Its asymptotics are then given
by\eq{ptwo_eql_asymptotics}. The long-time equilibrium susceptibility
which encodes the response to a
step change in the field is, using $\hatK\eql(0)=\int(dq)\, T/[2(z\eql+\omega)^2]$
and the generalization of\eq{L2hat0} to $T\neq T_{\rm c}$,
\begin{equation}
\fl \chi_{E}^{\rm eq}=\int_0^\infty dt\, R_{E}^{\rm eq}(t) = \frac{1}{4}\hat
L^{(2)}\eql(0) = \frac{T}{2}\left(1-\frac{1}{\int(dq)\, [T/(z\eql+\omega)]^2}\right)
\label{chi_Eql}
\end{equation}
It is easily shown that this is consistent with the known result for
the temperature dependence of the equilibrium energy, $E=\langle
H\rangle = \int(dq)\, \frac{1}{2}\omega[T/(z\eql+\omega)]$: one confirms $\chi_{E}^{\rm eq}
= T\,dE/dT$ as it should be. The factor $T$ arises because our field
$h$ is introduced via $H\to H-hH=(1-h)H$ and so corresponds to a
temperature change of $T/(1-h)-T=hT$ to linear order in $h$. The
inclusion of subleading non-Gaussian fluctuations is crucial for achieving this
consistency; as discussed at the beginning of Sec.~\ref{sec:setup},
the Gaussian theory does not even give the same answers for the
fluctuations of $H$ in its two representations as a bond or spin
product observable. The same phenomenon occurs in a purely static
calculation of the energy fluctuations and response~\cite{Joyce72}.

The temperature-dependence of the susceptibility\eq{chi_Eql} deserves
some comment. As $T$ approaches $T_{\rm c}$ from above, one has $z\eql\to
0$. For $d<4$, the denominator of the second term in\eq{chi_Eql} then
diverges, and $\chi_{E}^{\rm eq}/T$ smoothly approaches the value $1/2$ and
remains constant for $T<T_{\rm c}$. For $d>4$, the denominator has a finite
limit for $z\eql\to 0$. This produces the well-known discontinuity in
$\chi_{E}^{\rm eq}/T$ at $T=T_{\rm c}$, since for $T<T_{\rm c}$ the second term
in\eq{chi_Eql} again does not contribute. To see this explicitly, one
notes that for $T<T_{\rm c}$ the $\omega=0$ term in the spherical constraint
condition, with its weight $1/N$, has to be treated separately:
\begin{equation}
T\left(\int(dq)\, \frac{1}{z\eql+\omega} + \frac{1}{Nz\eql} \right)=1
\end{equation}
For $z\eql\to 0$ (on the scale ${{\mathcal O}}(N^0)$) the first integral is
$\int(dq)\, 1/\omega = 1/T_{\rm c}$ and so $z\eql = (1/N)TT_{\rm c}/(T_{\rm c}-T)$. This then
gives for the denominator integral in\eq{chi_Eql}
\begin{equation}
\int(dq)\, \frac{T^2}{(z\eql+\omega)^2} + \frac{T^2}{Nz\eql^2} \approx
N\frac{(T_{\rm c}-T)^2}{T_{\rm c}^2}
\end{equation}
which diverges for $N\to\infty$ at $T<T_{\rm c}$ as claimed.

An interesting question is how the discontinuity at $T=T_{\rm c}$ in $d>4$
of the susceptibility $\chi_{E}^{\rm eq}$, i.e.\ the {\em integrated} response,
is related to the temperature variation of the {\em time-dependent}
$R_{E}^{\rm eq}(t)$. In\eq{phat_eql} one notes that, for $T<T_{\rm c}$,
$\hatK\eql(s)$ acquires a distinct contribution from the $\omega=0$
mode
\begin{eqnarray}
\hatK\eql(s) &=&
\int(dq)\, \frac{T}{\omega(s+2\omega)} +
\frac{1}{N}\frac{T}{z\eql(s+2z\eql)}\\
& = & \int(dq)\, \frac{T}{\omega(s+2\omega)} + \frac{T_{\rm c}-T}{T_{\rm c} s}
\end{eqnarray}
where in the second line we have neglected $z\eql={{\mathcal O}}(1/N)$
against $s$. From\eq{R_Eql} one sees that the additional contribution
to $\hatK\eql(s)$ produces an extra pole in $\hat R_{E}^{\rm eq}(s)$, which
for small $(T_{\rm c}-T)/T_{\rm c}$ is located at $s=1/t_0$ with
$t_0=[TT_{\rm c}/4(T_{\rm c}-T)]\int(dq)\, \omega^{-2}$. Transforming to the time domain,
one finds that this gives an extra contribution of $\sim 1/t_0
\exp(-t/t_0)$ to $R_{E}^{\rm eq}(t)$. Crucially, this has a finite integral
even in the limit $T\to T_{\rm c}$, where $t_0$ diverges, and the appearance
of this term causes the discontinuity in $\chi_{E}^{\rm eq}$ at
$T=T_{\rm c}$. Translating these results to the time-dependent susceptibility
\begin{equation}
\chi_{E}^{\rm eq}(t) = \int_0^t dt' R_{E}^{\rm eq}(t')
\label{chi_Eql_t}
\end{equation}
one has that, for fixed finite $t$, $\chi_{E}^{\rm eq}(t)$ depends smoothly on
temperature around $T_{\rm c}$. For large $t$, $\chi_{E}^{\rm eq}(t)$ approaches a
plateau value equal to the equilibrium susceptibility at $T=T_{\rm c}^{+}$;
the approach to this plateau is as $\sim t^{(4-d)/2}$. For $T<T_{\rm c}$,
however, $\chi_{E}^{\rm eq}(t)$ eventually increases further on a diverging
timescale $t\sim t_0 \sim 1/(T_{\rm c}-T)$ to approach a higher plateau
value given by the susceptibility at $T=T_{\rm c}^{-}$. By FDT (see below)
the energy correlation function correspondingly shows a power-law
decay to a plateau for $t\ll t_0$, from which it falls to zero only
for $t>t_0$.

A final check on our results is that in equilibrium the energy
response and correlation function should be related by the usual
FDT. This is indeed the case. One combines the results\eq{C1},
\eq{C2}, \eq{C3}, \eq{C4} for the constituent parts of the correlation
function and decouples all the convolution integrals using temporal
Fourier transforms. It is easiest to do this starting from expressions
where all factors of $\Omega$ have been preserved, e.g.\ for $C^{(1)}_E$
one uses\eq{C1_first} rather than\eq{C1}. After some straightforward
but lengthy algebra one then indeed finds the Fourier domain version of
equilibrium FDT, $\int_{-\infty}^{\infty} dt\,C_{E}^{\rm eq}(t)e^{-i\nu t}=T[\hat
R_{E}^{\rm eq}(-i\nu)-\hat R_{E}^{\rm eq}(i\nu)]/(i\nu)$.

\section{Energy correlation and response: Non-equilibrium, $d>4$}
\label{sec:noneq_large_d}

We now evaluate the behaviour of the energy correlation and response
for the out-of-equilibrium dynamics after a quench to criticality. For
the correlation function, we need the long-time behaviour of
\begin{equation}
C\Om C(t,t_{\rm w}) = \int(dq)\, \frac{T_{\rm c}^2}{\omega^2}\sc{C}^2(\omega t_{\rm w})
\frac{g(t_{\rm w})}{g(t)} e^{-2\omega(t-t_{\rm w})}\omega
\end{equation}
At equilibrium, where $\sc{C}=1$ and $g(t)=$ const., this gives
$T_{\rm c}K\eql(t-t_{\rm w})$. In the non-equilibrium case one gets a correction
factor which becomes relevant in the aging regime, where  $t-t_{\rm w}\sim
t_{\rm w}\gg 1$:
\begin{eqnarray}
C\Om C(t,t_{\rm w}) &=& T_{\rm c}K\eql(t-t_{\rm w})\sc{C\Om C}\left(t/t_{\rm w}\right)
\label{COC}
\\
\sc{COC}(x) &=& \left(\frac{t}{t_{\rm w}}\right)^{\kappa}
\frac{\int dw \, w ^{(d-4)/2} e^{-2(x-1)w } \sc{C}^2(w)}
{\int dw \, w ^{(d-4)/2} e^{-2(x-1)w }}
\end{eqnarray}
This is valid for all dimensions $d>2$, but in the rest of this
section we focus on the case $d>4$, where $\kappa=0$. The regime $2<d<4$
is more complicated and treated separately in the next section.

Before assembling the four parts of our expression for $C_E(t,t_{\rm w})$,
it is useful to have a guide on what to expect overall. As discussed
above, the equilibrium energy fluctuations remain finite at
criticality. One therefore expects that, in the short-time regime
$t-t_{\rm w}\sim{{\mathcal O}}(1)$, $C_E(t,t_{\rm w})$ will decay as in
equilibrium. From\eq{ptwo_eql_asymptotics}, \eq{R_Eql} and FDT it
follows that this decay becomes a power law, $\sim(t-t_{\rm w})^{(4-d)/2}$
as $t-t_{\rm w}$ becomes large. Aging effects should then appear when
$t-t_{\rm w}\sim t_{\rm w}$ and manifest themselves through a scaling function of
$t/t_{\rm w}$. Overall, one expects in the aging regime
\begin{equation}
C_E(t,t_{\rm w})=(t-t_{\rm w})^{(4-d)/2}\times[\mbox{scaling function of
}t/t_{\rm w}]
\label{CE_expect}
\end{equation}
and this is indeed what we will find. Consider now $C_E^{(1)}$, for
which we obtained in\eq{C1} that $4C_E^{(1)} = -(\partial_t+2z(t))C\Om C(t,t_{\rm w})$.
Now from\eq{g_def} and the fact that $g(t)\to$ const.\ for
$t\to\infty$ and $d>4$ it follows that $z(t)=g'(t)/[2g(t)]$
decreases more quickly than $1/t$. Using also\eq{COC}, we see that the
$z(t)$ term is negligible for large times and that $4C_E^{(1)}$
decays as $(t-t_{\rm w})^{-d/2}$ times an aging correction. This is
negligible compared to\eq{CE_expect}, so we do not need to analyse
this contribution further.

For $C_E^{(2)}$, we have from\eq{C2}
\begin{equation}
\fl 4C_E^{(2)} = 2z(t) C\Om C(t,t_{\rm w})-2T_{\rm c} C\Om C(t,t_{\rm w}) + \int dt' L^{(2)}(t,t')
C\Om C(t',t_{\rm w})
\label{C2_large_d}
\end{equation}
Because of the factor $C\Om C(t',t_{\rm w})$ it is useful to split the
$t'$-integral into the regimes $t'>t_{\rm w}$ and $t'<t_{\rm w}$. The first regime
gives
\begin{equation}
T_{\rm c}\int_{t_{\rm w}} dt' L^{(2)}\eql(t-t')
K\eql(t'-t_{\rm w})
\sc{L}(t/t')
\sc{C\Om C}(t'/t_{\rm w})
\label{C2_part}
\end{equation}
As in the analysis of\eq{ptwo_scaling_cond} one can now argue that,
for large $t-t_{\rm w}$, the integral will be dominated by the regions
$t'\approx t$ and $t'\approx t_{\rm w}$. The weights contributed by these
regions are $K\eql(t-t_{\rm w})\int_0^{\infty} dt'L^{(2)}\eql(t')$ and
$L^{(2)}\eql(t-t_{\rm w})\int_0^{\infty} dt'K\eql(t')$, respectively; the
integrals are both finite, so that both terms scale as
$(t-t_{\rm w})^{(2-d)/2}$. These weights are then multiplied by the relevant
values of scaling functions in the integrand of\eq{C2_part}. The
overall result is smaller by $\sim 1/(t-t_{\rm w})$ than the leading
term\eq{CE_expect} in $C_E$, and can be neglected. The $t'<t_{\rm w}$ part
of the integral from\eq{C2_large_d} is
\begin{equation}
T_{\rm c}\int^{t_{\rm w}} dt' L^{(2)}\eql(t-t')
K\eql(t_{\rm w}-t') \sc{L}(t/t') \sc{C\Om C}(t_{\rm w}/t')
\end{equation}
The factor $K\eql(t_{\rm w}-t')$ concentrates the weight of the integrand
near $t'=t_{\rm w}$, so that for $t-t_{\rm w}$ and $t_{\rm w}\gg 1$ all other factors
can be replaced by their values at $t'=t_{\rm w}$, giving the result
$T_{\rm c}\hatK\eql(0) L^{(2)}\eql(t-t_{\rm w})\sc{L}(x)$. This scales as
$(t-t_{\rm w})^{(2-d)/2}$ times a scaling function of $x$ and so is also
negligible compared to the leading contribution\eq{CE_expect}.
From\eq{COC} the second term on the r.h.s.\ of\eq{C2_large_d} has the
same subleading scaling, and the first term is even smaller. Overall,
$C_E^{(2)}$ therefore gives a subleading contribution to $C_E$ in the
aging regime. By very similar arguments one shows that $C_E^{(3)}$
from\eq{C3} is also subleading.

In the aging regime we therefore only need to consider $C_E^{(4)}$
from\eq{C4}. The long-time behaviour of the function $CC(t,t_{\rm w})$ is
easily worked out as
\begin{eqnarray}
\flCC(t,t_{\rm w}) &=& CC\eql(t-t_{\rm w})\sc{CC}(t/t_{\rm w})
\label{CC_decomp}
\\
\flCC\eql(t-t_{\rm w}) &=& \int(dq)\, \frac{T_{\rm c}^2}{\omega^2}e^{-2\omega(t-t_{\rm w})} = 2T_{\rm c}
\int_{t-t_{\rm w}}^\infty dt'\, K\eql(t') \sim (t-t_{\rm w})^{(4-d)/2}
\label{CC_eql}
\\
\fl\sc{CC}(x) &=& \frac{\int dw \, w ^{(d-6)/2} e^{-2(x-1)w }\sc{C}^2(w )}
{\int dw \, w^{(d-6)/2} e^{-2(x-1)w}}
\label{CC_scaling}
\end{eqnarray}
where the scaling function $\sc{CC}(x)\sim x^{-2}$ for $x\gg 1$. Now
consider the $t'$-integral from\eq{C4},
\begin{equation}
I(t,t_{\rm w}') = \int dt' \,
[(2z(t)-2T_{\rm c})\delta(t-t')+L^{(2)}(t,t')]CC(t',t_{\rm w}')
\label{Idef}
\end{equation}
The equilibrium part $L^{(2)}\eql(t-t')$ of $L^{(2)}(t,t')$ ensures that
this has its mass concentrated near $t'=t$. Because we are looking at
the aging regime, we have $t-t_{\rm w}'(>t-t_{\rm w})\gg 1$, and\eq{CC_eql} then
shows that the function $CC(t',t_{\rm w}')$ is slowly varying near
$t'=t$. It can thus be replaced by its value there, $CC(t,t_{\rm w}')$, to
give to leading order
\begin{equation}
I(t,t_{\rm w}') = (2z(t)-2T_{\rm c}+\hatL^{(2)}\eql(0))CC(t,t_{\rm w}') =
\frac{1}{\hatK\eql(0)}CC(t,t_{\rm w}')
\label{I_approx}
\end{equation}
where in the last step we have used $z(t)\ll 1$. One might suspect
that in the $t'$-integral of\eq{Idef}, the region around $t'\approx
t_{\rm w}'$ also contributes, since this is where $CC(t',t_{\rm w}')$ is
largest. However, it can be shown that this region only gives a
subleading contribution compared to\eq{I_approx}.

In terms of $I(t,t_{\rm w}')$ we can write $C_E^{(4)}$ using\eq{C4} as
\begin{equation}
4C_E^{(4)} =
\frac{1}{2} \int dt_{\rm w}' [(2z(t_{\rm w})-2T_{\rm c})\delta(t_{\rm w}-t_{\rm w}')+L^{(2)}(t_{\rm w},t_{\rm w}')]I(t,t_{\rm w}')
\label{C4_I}
\end{equation}
Again, the integral is dominated by the region $t_{\rm w}'\approx t_{\rm w}$
because of the factor $L^{(2)}\eql(t_{\rm w}-t_{\rm w}')$ and we thus get our final
expression
\begin{equation}
4C_E^{(4)} =
\frac{1}{2\hatK\eql(0)}I(t,t_{\rm w}) =
\frac{1}{2\hatK\eql^2(0)}CC(t,t_{\rm w})
\label{C4_final}
\end{equation}
As argued above, in the aging regime the full energy correlator will
be identical to this.

We can now turn to the evaluation of the out-of-equilibrium response
function. We will not write explicitly the second term in\eq{RE},
which just removes the unwanted $\delta(t-t_{\rm w})$ contribution arising
from the derivative applied to the first term's step
discontinuity. Thus,
\begin{equation}
2R_E = (-\partial_t-2T_{\rm c})R\Om C(t,t_{\rm w}) + \int dt' L^{(2)}(t,t') R\Om C(t',t_{\rm w})
\label{RE_again}
\end{equation}
The long-time scaling of the function $R\Om C$ is
\begin{eqnarray}
\fl R\Om C(t,t_{\rm w}) &=& R\Om C\eql(t-t_{\rm w})\sc{R\Om C}(t/t_{\rm w})
\\
\fl R\Om C\eql(t-t_{\rm w}) &=& T_{\rm c} \int(dq)\, e^{-2\omega(t-t_{\rm w})} = T_{\rm c} f(t-t_{\rm w}) \sim
(t-t_{\rm w})^{-d/2}
\label{ROC_eql}
\\
\fl \sc{R\Om C}(x) &=& \frac{\int dw\, w^{(d-2)/2} e^{-2(x-1)w}\sc{C}(w)}
{\int dw\, w^{(d-2)/2} e^{-2(x-1)w}}
\label{ROC_scaling}
\end{eqnarray}
The integrand in\eq{RE_again} has its mass concentrated near $t'=t$,
from the factor $L^{(2)}\eql(t-t')$, and near $t'=t_{\rm w}$, from the factor
$R\Om C\eql(t'-t_{\rm w})$. This gives
\begin{equation}
\fl 2R_E = (-\partial_t-2T_{\rm c})R\Om C(t,t_{\rm w}) + \hatL^{(2)}\eql(0) R\Om C(t,t_{\rm w}) +
\frac{1}{2}L^{(2)}\eql(t-t_{\rm w})\sc{L}(t/t_{\rm w})
\label{RE_aux}
\end{equation}
using that $\int_0^\infty dt'\, R\Om C\eql(t') = \frac{1}{2} 2T_{\rm c} \int_0^\infty
dt'\, f(t') = \frac{1}{2} K\eql(0) = \frac{1}{2}$ (see after\eq{p_eql}). The last
term in\eq{RE_aux}
scales as $(t-t_{\rm w})^{(2-d)/2}$ times a function of $x=t/t_{\rm w}$, while the
other terms scale at most as $(t-t_{\rm w})^{-d/2}$ and so are
subdominant. Using also that $K\eql$ and $L^{(2)}\eql$ are
asymptotically proportional for $d>4$ from\eq{ptwo_p_link} and that
$\sc{L}=\sc{K}$ then gives the final aging regime expression
\begin{equation}
R_E =
\frac{1}{4}L^{(2)}\eql(t-t_{\rm w})\sc{L}(t/t_{\rm w}) =
\frac{1}{4}L^{(2)}(t,t_{\rm w})
\label{RE_final}
\end{equation}
Interestingly, the simple relation\eq{R_Eql} between the energy
response and $L^{(2)}$ therefore holds not only at equilibrium but also in
the aging regime of the non-equilibrium dynamics, and thus overall
across the entire long-time regime.

We can now, finally, find the FDR for the energy correlation and
response. In the aging regime of interest here, it is useful to
rewrite the response function\eq{RE_final} as
\begin{equation}
R_E(t,t_{\rm w}) = \frac{1}{4\hatK\eql^2(0)}K\eql(t-t_{\rm w})\sc{K}(x) =
\frac{1}{4\hatK\eql^2(0)}K(t,t_{\rm w})
\label{RE_final_K}
\end{equation}
using\eq{ptwo_p_link} and $\sc{L}=\sc{K}$. Recalling the
definition\eq{p_def}, and the fact that $C_E$ is given by\eq{C4_final}
to leading order in the aging limit, we then obtain
\begin{equation}
X_E(t,t_{\rm w}) = \frac{T_{\rm c} R_E(t,t_{\rm w})}{C'_E(t,t_{\rm w})} = \frac{\int(dq)\, T_{\rm c} R_\mathbf{q}
C_\mathbf{q}}{\int(dq)\, C'_\mathbf{q} C_\mathbf{q}}
\label{XE}
\end{equation}
Remarkably, this is in fact {\em identical} to the FDR\eq{Xbl_large_d}
for large blocks of spin product observables. In particular, $X_E$
tends to the limit value $X^\infty=1/2$ for $t\ggt_{\rm w}$, and this is
identical to the one we found for spin and bond observables: for
$d>4$, all observables we have considered lead to a unique value of
$X^\infty=1/2$.

The FD plot for the energy has a limiting pseudo-equilibrium
form; this follows from\eq{CC_eql} and\eq{C4_final}, which show that
$C_E(t,t_{\rm w})$ has decayed to a small value $\simt_{\rm w}^{(4-d)/2}$ at the
point where aging effects begin to appear. More generally, the
correlation function scales as $C_E \sim (t-t_{\rm w})^{(4-d)/2}$ for
$x\approx 1$, and for $x\gg1$, where $\sc{CC}(x) \sim x^{-2}$
from\eq{CC_scaling}, as $C_E \sim t^{(4-d)/2}x^{-2} =
t_{\rm w}^{2}t^{-d/2}$.

It should be pointed out that while the FDR for the energy matches
that for blocks of product observables at long times, the correlation
and response functions themselves do not. This follows from the
nontrivial proportionality factors $1/{\hatK\eql^2(0)}$
in\eq{C4_final} and\eq{RE_final_K}. The latter are required to match
the limiting behaviour for $t-t_{\rm w}\llt_{\rm w}$ of the aging regime results
to the asymptotics of the equilibrium results for $t-t_{\rm w}\gg
1$. Indeed, by combining the effective equilibrium behaviour for
$t-t_{\rm w}={{\mathcal O}}(1)$ with the above aging expressions we can write
\begin{eqnarray}
C_E(t,t_{\rm w}) &=& C_{E}^{\rm eq}(t-t_{\rm w})\sc{CC}(x)
\label{CE_d_gt_4_longtime}
\\
R_E(t,t_{\rm w}) &=& R_{E}^{\rm eq}(t-t_{\rm w})\sc{K}(x)
\label{RE_d_gt_4_longtime}
\end{eqnarray}
and these are now valid throughout the long-time regime, i.e.\ for
large $t_{\rm w}$ but independently of whether $t-t_{\rm w}$ is large or not. As
promised they match with the aging expressions for $1\ll \Delta t
\llt_{\rm w}$. For the response this is obvious from\eq{R_Eql}. For $C_E$,
\eq{ptwo_p_link} and\eq{R_Eql} show that $-\partial_{\Delta t} C_{E}^{\rm eq}(\Delta t) =
(T_{\rm c}/4)L^{(2)}\eql(\Delta t)\approx [T_{\rm c}/4\hatK\eql^2(0)]K\eql(\Delta t)$ for
large $\Delta t$; from\eq{CC_eql} this agrees with the corresponding
derivative of the result in\eq{C4_final},
$-\partial_{\Delta t}CC\eql(\Delta t)/[8\hatK\eql^2(0)] =
[2T_{\rm c}/8\hatK\eql^2(0)]K\eql(\Delta t)$. Note from the discussion
after\eq{chi_Eql_t} that the function $C_{E}^{\rm eq}(\Delta t)$ discontinuously
acquires an additive constant (non-decaying) part as $T$ crosses $T_{\rm c}$ from
above in $d>4$. What we mean in\eq{CE_d_gt_4_longtime} is the limiting
form for $T\searrow T_{\rm c}$ {\em from above}, which does not contain this
plateau. That this is the correct choice can be seen as follows: the
non-decaying part of $C_{E}^{\rm eq}(\Delta t)$ at $T=T_{\rm c}^{-}$ arises from the
fluctuations of the $\mathbf{q}=\mathbf{0}$ Fourier mode, i.e.\ of the magnetization,
which are larger by ${{\mathcal O}}(N)$ than those of the other Fourier
modes. In the context of our non-equilibrium calculation, all Fourier
modes have fluctuations of comparable order (in system size), so that
this contribution is absent; it would appear only for times $t_{\rm w}$ that
scale with system size $N$.

Finally, the long-time
expressions\eqq{CE_d_gt_4_longtime}{RE_d_gt_4_longtime} also show that
an energy FD plot would have a pseudo-equilibrium form at long times,
because e.g.\ in $C_E(t,t_{\rm w})$ the equilibrium factor has already
decayed to $\sim t_{\rm w}^{(4-d)/2}$ of its initial value when aging
effects appear around $\Delta t\sim t_{\rm w}$.

\section{Energy correlation and response: Non-equilibrium, $d<4$}
\label{sec:noneq_small_d}

In dimension $d<4$, the evaluation of the energy correlation function in the
aging regime is somewhat more complicated. One can nevertheless show
that, as before, the dominant contribution to $C_E$ is $C_E^{(4)}$,
with the other three terms being subleading. We thus again need to
consider the function $I(t',t_{\rm w})$ defined in\eq{Idef}, and then work
out $C_E$ using\eq{C4_I}. One can see clearly that the approach
leading above to\eq{I_approx} no longer works: in $d<4$,
$\hatL^{(2)}\eql(0)=2T_{\rm c}$, so the leading terms in\eq{I_approx} actually
cancel. To treat this cancellation accurately, we
use\eq{ptwo_eql_condition} to write the coefficient $2T_{\rm c}$ in the
following way
\begin{equation}
2T_{\rm c} = - \frac{K\eql'(t)}{K\eql(t)} + \int dt'\,
\frac{K\eql(t')}{K\eql(t)} L^{(2)}\eql(t-t')
\end{equation}
This allows us to express\eq{Idef} as
\begin{eqnarray}
\fl I(t,t_{\rm w}') &=& \left[2z(t)+\frac{K\eql'(t)}{K\eql(t)}\right]
CC(t,t_{\rm w}')
\label{I_subtract}
\\
\fl & &  {}+{} \int_0^t
dt'\,L^{(2)}\eql(t-t') \left[\sc{L}(t/t')CC(t',t_{\rm w}')
-\frac{K\eql(t')}{K\eql(t)}CC(t,t_{\rm w}')\right]
\nonumber
\end{eqnarray}
To make progress, we need the behaviour of $CC$ for $d<4$. For long
times we have
\begin{equation}
\fl CC(t,t_{\rm w}) = \int(dq)\, C^2_\mathbf{q}(t,t_{\rm w}) = T_{\rm c}^2 \frac{g(t_{\rm w})}{g(t)}\int(dq)\,
\omega^{-2} \sc{C}^2(\omegat_{\rm w}) e^{-2\omega(t-t_{\rm w})}
\label{CC_gen_dlt4}
\end{equation}
The equal-time value thus scales as
\begin{equation}
CC(t,t) \sim \int d\omega\, \omega^{d/2-3} \sc{C}^2(\omega t) =
t^{(4-d)/2}\int dw \, w ^{(d-6)/2} \sc{C}^2(w )
\label{CC_equal_time}
\end{equation}
hence $CC(t,t)=\gamma_d t^{(4-d)/2}$ with some constant $\gamma_d$. A scaling
function is then obtained if we normalize $CC(t,t_{\rm w})$ by this
equal-time value:
\begin{equation}
\frac{CC(t,t_{\rm w})}{CC(t,t)} =
\frac{g(t_{\rm w})}{g(t)}
\frac{t_{\rm w}^{(4-d)/2}}{t^{(4-d)/2}}
\frac{\int dw \, w ^{(d-6)/2} \sc{C}^2(w ) e^{-2w (t-t_{\rm w})/t_{\rm w}}}
     {\int dw \, w ^{(d-6)/2} \sc{C}^2(w )}
\label{sc_CC_ord}
\end{equation}
This is a function ${\mathcal{G}}(x)$ of $x=t/t_{\rm w}$; note that the first two
fractions cancel since $g(t)\sim t^{-\kappa}=t^{(d-4)/2}$. So far this
scaling expression holds for $t>t_{\rm w}$. To be able to use it also for
non-ordered times, note that for $t<t_{\rm w}$,
\begin{equation}
\frac{CC(t,t_{\rm w})}{CC(t,t)}=\left(\frac{t_{\rm w}}{t}\right)^{(4-d)/2}
\frac{CC(t_{\rm w},t)}{CC(t_{\rm w},t_{\rm w})}=x^{(d-4)/2}
{\mathcal{G}}(t_{\rm w}/t)
\label{sc_CC_dis}
\end{equation}
So we have overall
\begin{equation}
\fl\frac{CC(t,t_{\rm w})}{CC(t,t)} ={\mathcal{G}}(t/t_{\rm w}),
\quad
{\mathcal{G}}(x)=\left\{
\begin{array}{ll}
{\displaystyle\frac{\int dw\, w^{(d-6)/2} \sc{C}^2(w) e^{-2(x-1)w}}
{\int dw\, w^{(d-6)/2} \sc{C}^2(w)}} & \mbox{for $x\geq 1$} \\
x^{(d-4)/2}{\mathcal{G}}(1/x) & \mbox{for $x\leq 1$}
\end{array}
\right.
\label{CC_d_lt_4}
\end{equation}
One easily works out the asymptotics of ${\mathcal{G}}$: ${\mathcal{G}}(x)\sim x^{-d/2}$ for
$x\to\infty$, ${\mathcal{G}}(x)\sim x^{d-2}$ for $x\to 0$. For $|x-1|\ll1$, on
the other hand, $1-{\mathcal{G}}(x)\sim |x-1|^{(4-d)/2}$. This corresponds to
$CC(t,t)-CC(t,t_{\rm w})\sim |t-t_{\rm w}|^{(4-d)/2}$ for $1\ll |t-t_{\rm w}|\ll
t_{\rm w}$. (The behaviour of this difference for smaller time intervals
$t-t_{\rm w}={{\mathcal O}}(1)$ is not captured accurately by the scaling
form\eq{CC_d_lt_4}, but integrals over this regime contribute
only subleading corrections to the results below.)

We can now insert the scaling form\eq{CC_d_lt_4} of $CC$
into\eq{I_subtract}. The non-integral terms turn out to be subleading
(see below), so
\begin{eqnarray}
\fl I(t,t_{\rm w}') &=& \int dt' \, L^{(2)}\eql(t-t')
\nonumber\\
\fl & &\times\left[\sc{L}(t/t')CC(t',t')
{\mathcal{G}}(t'/t_{\rm w}')
-\frac{K\eql(t')}{K\eql(t)}CC(t,t){\mathcal{G}}(t/t_{\rm w}')\right]
\\
\fl &=& CC(t,t) \int dt' \,
L^{(2)}\eql(t-t')\left[\frac{t'}{t}{\mathcal{G}}(t'/t_{\rm w}') - \frac{K\eql(t')}{K\eql(t)}
{\mathcal{G}}(t/t_{\rm w}')\right]
\label{I_sub2}
\end{eqnarray}
where we have used that, from\eq{FKd_lt_4} and\eq{same_sc_fn},
$\sc{L}(t/t')CC(t',t')/CC(t,t) = (t'/t)^{(d-2)/2}
(t'/t)^{(4-d)/2} = t'/t$. The remaining ``subtracted'' integral now no
longer has its mass concentrated near $t'=t$ because the terms in
square brackets give a factor $t-t'$ there. The whole integration
range contributes, so that we can replace $L^{(2)}\eql(t-t')$ by its
asymptotic form\eq{ptwo_eql_asymptotics}, $L^{(2)}\eql(t-t')\approx\lambda_d
(t-t')^{(d-6)/2}$ with a $d$-dependent constant $\lambda_d$. Similarly, the
ratio ${K\eql(t')}/{K\eql(t)}$ can be replaced by $(t'/t)^{(2-d)/2}$
to leading order. Scaling the integration variable as $y=t'/t_{\rm w}'$ then
gives
\begin{equation}
\fl I(t,t_{\rm w}') = I(x')
= \gamma_d\lambda_d x'^{-2}\int_0^{x'} dy\,(1-y/x')^{(d-6)/2}y\left[{\mathcal{G}}(y) -
y^{-d/2}{\mathcal{G}}(x')x'^{d/2}\right]
\label{I_sub3}
\end{equation}
This shows that in the aging regime $I$ depends on $x'=t/t_{\rm w}'$
only. One can now also see that the neglected terms
from\eq{I_subtract} are indeed subleading: they scale as
$t^{-1}CC(t,t_{\rm w}')\sim t^{(2-d)/2}{\mathcal{G}}(t/t_{\rm w}')$.

We can now proceed to simplifying\eq{C4_I} in the aging regime. Using
the same subtraction method as above, and remembering that
$C_E=C_E^{(4)}$ to leading order, one finds by analogy with\eq{I_subtract}
\begin{equation}
\fl 8C_E(t,t_{\rm w}) = \int dt_{\rm w}' \,
L^{(2)}\eql(t_{\rm w}-t_{\rm w}')\left[\sc{L}(t_{\rm w}/t_{\rm w}')
I(t,t_{\rm w}') - \frac{K\eql(t_{\rm w}')}{K\eql(t_{\rm w})}I(t,t_{\rm w})\right]
\label{C4_sub}
\end{equation}
Here subleading terms similar to those in the first line
of\eq{I_subtract} have already been neglected. With the scaled
integration variable $y=t_{\rm w}'/t_{\rm w}$, and using the asymptotic forms of
$L^{(2)}\eql(t_{\rm w}-t_{\rm w}')$ and $K\eql(t_{\rm w}')$ as in\eq{I_sub3} one gets
\begin{equation}
\fl 8C_E(t,t_{\rm w}') = \lambda_d t_{\rm w}^{(d-4)/2}
\int_0^1 dy \, (1-y)^{(d-6)/2}\left[y^{(d-2)/2} I(x/y) - y^{(2-d)/2}I(x)\right]
\label{CE_sub}
\end{equation}
This shows that $C_E$ scales as $t_{\rm w}^{(d-4)/2}$ times a function of
$x=t/t_{\rm w}$. It is difficult to work out the whole functional dependence
on $x$. We therefore focus below on the asymptotic behaviour for
$x\to\infty$, which gives the asymptotic FDR. First, though, we check that in
the limit $x\to 1$, where aging effects are unimportant, our
expression matches with the equilibrium result
$C_{E}^{\rm eq}(t-t_{\rm w})$. From\eq{R_Eql} the latter behaves as $C_{E}^{\rm eq}(\Delta t) =
(T_{\rm c}/4)\int_{\Delta t}^\infty dt'\,L^{(2)}\eql(t') \approx
[\LdT_{\rm c}/2(4-d)]\Delta t^{(d-4)/2}$ for large $\Delta t$. To compare with the
prediction from\eq{CE_sub}, one uses that $I(x') \sim \ln(x'-1)$ for
$x'\approx1$; this follows from the behaviour of ${\mathcal{G}}(x)$ for $x\approx
1$. (Note that the proportionality constant in $I\sim \ln(x'-1)$ is
{\em positive}, so that $I$ itself is -- in contrast to the case $d>4$
-- {\em negative}.) Inserting into\eq{CE_sub} one then finds that for
$x\approx 1$ the integral scales as $(x-1)^{(d-4)/2}$. Overall, one
gets $C_E(t,t_{\rm w}) \sim t_{\rm w}^{(d-4)/2}(x-1)^{(d-4)/2} =
(t-t_{\rm w})^{(d-4)/2}$, consistent with the expectation from the
equilibrium result. One can work out the prefactor of the power law
and finds that this, too, agrees as it should.

Turning now to the behaviour of\eq{CE_sub} for large $x$, we need the
asymptotics of $I(x')$.
The leading tail $\sim y^{-d/2}$ of ${\mathcal{G}}(y)$ is subtracted off in the
expression in square brackets in\eq{I_sub3}, leaving as the next term
$y^{-(d+2)/2}$. Even together with the additional factor $y$ this
makes the integral convergent at the upper end for $x'\to\infty$. In
the limit we thus get $I(x') \approx -\alpha_d\gamma_d\lambda_d x'^{-2}$, with
\begin{equation}
\alpha_d = \int_0^\infty dy\,y\left[g_d y^{-d/2}-{\mathcal{G}}(y)\right]
\label{Ad}
\end{equation}
and $g_d=\lim_{x\to\infty} {\mathcal{G}}(x)x^{d/2}$.

Inserting this inverse-square asymptotic behaviour of $I(x')$ gives
for the integral in\eq{CE_sub} the scaling $\alpha_d\beta_d\gamma_d\lambda_d x^{-2}$ with
\begin{equation}
\fl \beta_d = \int_0^1 dy \,
(1-y)^{(d-6)/2}\left(y^{(2-d)/2}-y^{(d+2)/2}\right)
= -
\frac{\Gamma\left(\frac{d-4}{2}\right)\Gamma\left(\frac{d+4}{2}\right)}
{\Gamma(d)}
\label{CE_factor}
\end{equation}
It then follows, finally, that $C_E(t,t_{\rm w})=(\alpha_d\beta_d\gamma_d\lambda_d^2/8)
t_{\rm w}^{d/2} t^{-2}$ for $t\ggt_{\rm w}$, and
\begin{equation}
C'_E(t,t_{\rm w}) = (\alpha_d\beta_d\gamma_d\lambda_d^2 d/16) t_{\rm w}^{(d-2)/2} t^{-2}
\label{CE_prime}
\end{equation}

To get the FDR, we now need the response function $R_E$. This can be
evaluated very similarly to the case $d>4$ and one finds that the last
term in\eq{RE_aux} is again the dominant one, giving
\begin{equation}
R_E(t,t_{\rm w}) = \frac{1}{4}L^{(2)}\eql(t-t_{\rm w})\sc{L}(t/t_{\rm w}) =
\frac{1}{4}L^{(2)}(t,t_{\rm w})
\label{RE_d_lt_4}
\end{equation}
(The $R\Om C(t,t_{\rm w})$ terms in\eq{RE_aux} look dangerous: they scale as
$(t-t_{\rm w})^{-d/2}$ and are thus {\em larger} than $L^{(2)}\eql(t-t_{\rm w})\sim
(t-t_{\rm w})^{(d-6)/2}$ in $d<3$. However, their prefactors
$-2T_{\rm c}+\hatL^{(2)}\eql(0)$ cancel. Treating this cancellation more
carefully then shows that these terms do remain subleading compared
to\eq{RE_d_lt_4}.) As before, $4R_E$ equals $L^{(2)}$ in the aging
regime, and this then holds across the whole long-time
regime since for $t-t_{\rm w}={{\mathcal O}}(1)$ it matches the equilibrium
behaviour $4R_E(t,t_{\rm w})=L^{(2)}\eql(t-t_{\rm w})$. For $t\ggt_{\rm w}$, on the other
hand, the response\eq{RE_d_lt_4} becomes
\begin{equation}
R_E(t,t_{\rm w}) = \frac{1}{4}\lambda_d t^{(d-6)/2}
\left(\frac{t}{t_{\rm w}}\right)^{(2-d)/2} = \frac{\lambda_d}{4} t_{\rm w}^{(2-d)/2} t^{-2}
\end{equation}
Comparing with\eq{CE_prime} then finally gives for the asymptotic FDR
for the energy in $d<4$,
\begin{equation}
X_E^\infty = \frac{4T_{\rm c}}{\alpha_d\beta_d\gamma_d\lambda_d d}
\label{Xinf_dlt4}
\end{equation}
After evaluating the various numerical factors (see
\ref{sec:X_E_infty}) this can be written as
\begin{equation}
X_E^\infty = \frac{2}{d\tilde{\alpha}_d}
\frac{\Gamma(d)\Gamma{\textstyle\left(\frac{4-d}{2}\right)}}
{\Gamma\left(\frac{d+4}{2}\right)}
\label{Xinf_dlt4_explicit}
\end{equation}
where
\begin{equation}
\fl\tilde{\alpha}_d =  \frac{\pi}{\sin[\pi(4-d)/2]}
+\int_1^\infty dx
\int_0^1 dy\, y^{(d-4)/2} \frac{(x-y)^{(2-d)/2}}{1+x-y}
\left(1-y-x^{-(d+2)/2}\right)
\label{Atd}
\end{equation}
Near $d=4$, one can expand in $\epsilon=4-d$.  It can be shown by
explicit calculation that the integral in\eq{Atd} is exactly zero for
$d=4$, giving $\tilde{\alpha}_d = 2/\epsilon+{{\mathcal O}}(\epsilon)$ and so
\begin{equation}
X_E^\infty = \frac{\epsilon
\Gamma(4-\epsilon)\Gamma(\epsilon/2)}{(4-\epsilon)\Gamma(4-\epsilon/2)}
+ {{\mathcal O}}(\epsilon^2) = \frac{1}{2}-\frac{\epsilon}{3} +
{{\mathcal O}}(\epsilon^2)
\end{equation}
Notably, this is {\em different} from the FDR
$X^\infty=1-2/d=1/2-\epsilon/8+{{\mathcal O}}(\epsilon^2)$ for all the other,
finite-range, observables that we considered previously in $d<4$. It is,
however, consistent with RG calculations to
${{\mathcal O}}(\epsilon)$ for the $O(n)$ model in the limit $n\to\infty$, for
an analogous choice of observable~\cite{CalGam04}. Remarkably,
therefore, the
non-Gaussian fluctuations induced by the weak infinite-range interaction
in the spherical model seem to mimick precisely the effects that are
seen in more realistic models such as the $O(n)$, even though in the
latter all interactions are short-ranged and there is no difference
between the behaviour of block observables and global ones.

\begin{figure}
\centerline{\includegraphics[width=11.0cm,clip=true]{XE.eps}}
\caption{Asymptotic FDR $X^\infty$ vs dimension $d$ for finite-range
observables (equation~(\protect\ref{X_baseline}), dashed) and energy
(equation~(\protect\ref{Atd}), solid). Dotted lines indicate the
corresponding linear expansions $1/2-\epsilon/8$ and $1/2-\epsilon/3$ in
$\epsilon=4-d$.
\label{fig:X_normal_and_E}
}
\end{figure}
We show in Fig.~\ref{fig:X_normal_and_E} the dependence of the
asymptotic energy-FDR $X_E^\infty$ on dimension $d$ and compare with
the result for finite-range observables. They agree in $d\geq 4$, but the
difference between them grows as $d$ decreases towards $d=2$, with the
energy FDR always having the lower value. In the limit $d\to 2$, both
FDRs converge to zero, but while the finite-range FDR
$X^\infty=\epsilon'/2+{{\mathcal O}}(\epsilon'^2)$ does so linearly in
$\epsilon'=d-2$, the energy FDR vanishes quadratically as
$X_E^\infty=\epsilon'^2/8+{{\mathcal O}}(\epsilon'^3)$, due to the divergence
of $\tilde{\alpha}_d=4/\epsilon'^2+{{\mathcal O}}(1/\epsilon')$.

As in the case $d>4$, an energy FD plot would have a
pseudo-equilibrium form which hides all non-equilibrium effects at
long times. Indeed, one could write\eq{CE_sub} in the form
$C_E(t,t_{\rm w})=C_{E}^{\rm eq}(t-t_{\rm w})\sc{{C_E}}(x)$, and the decay of the
equilibrium factor $C_{E}^{\rm eq}$ squeezes all details about the aging factor
$\sc{{C_E}}(x)$ into a vanishingly small region of the FD plot for
long times. While we have not calculated $\sc{{C_E}}$ explicitly, the
discussion after\eq{CE_sub} shows that $\sc{{C_E}}(1)=1$ as it should
be. For large $x$, on the other hand, we saw in \eq{CE_prime} that
$C_E(t,t_{\rm w})\sim t_{\rm w}^{d/2} t^{-2}$, implying $\sc{{C_E}}(x)\sim
t_{\rm w}^{d/2}t^{-2}(t-t_{\rm w})^{(4-d)/2} \sim x^{-d/2}$. This matches
continuously at $d=4$ with the prediction\eq{CE_d_gt_4_longtime} for
$d>4$, where the aging correction decays as $\sc{CC}(x) \sim x^{-2}$.

\section{Magnetized initial states}
\label{sec:magnetized}

In this final section we consider the dynamics for initial
configurations with nonzero
magnetization, focussing as before on the non-equilibrium
dynamics that results when the system is subsequently brought to the
critical temperature $T_{\rm c}$. Physically, such a situation could arise
in an ``up-quench'', where the system is equilibrated at $T<T_{\rm c}$ and
temperature is then increased to $T=T_{\rm c}$. As explained in the
introduction, our interest in this
scenario arises from recent results~\cite{GarSolPagRit05} which show
that such initial conditions produce FDRs that differ nontrivially
from those for zero magnetization. The analysis
of~\cite{GarSolPagRit05} was limited to high dimensions $d$ or
infinite-range interactions, however; the calculation below will allow us
to see explicitly how the results change in finite dimension. In
particular, we will obtain exact FDRs for magnetized coarsening
below the upper critical dimension, i.e.\ $d<4$ in the spherical
model.

We will continue to use the notation $C_\mathbf{q}(t,t_{\rm w})=(1/N) \left\langle
S_\mathbf{q}(t)S_\mathbf{q}^*(t)\right\rangle$ for Fourier mode correlations. For $\mathbf{q}=\mathbf{0}$
this is now a full, unsubtracted correlator, with
$C_\mathbf{0}(t,t)=Nm^2(t)+{{\mathcal O}}(1)$ and $m(t)=(1/N)\left\langle S_\mathbf{0}(t)\right\rangle$ the
time-dependent magnetization. The difference $\Cc_\zv(t,t_{\rm w}) =
C_\mathbf{0}(t,t_{\rm w})-Nm(t)m(t_{\rm w})$ is the connected correlator, which has
values of ${{\mathcal O}}(1)$ and is the relevant quantity for analysing the
FD behaviour. For $\mathbf{q}\neq\mathbf{0}$, connected and full
correlators coincide:
\begin{equation}
\tilde{C}_\mathbf{q}(t,t_{\rm w}) = C_\mathbf{q}(t,t_{\rm w})-(1/N)\left\langle S_\mathbf{q}(t)\right\rangle \langle
S^*_\mathbf{q}(t_{\rm w})\rangle = C_\mathbf{q}(t,t_{\rm w})
\end{equation}

We now need to check how the analysis in the previous sections is
modified by the
presence of a nonzero magnetization. The Fourier space equation of
motion\eq{dotSq} remains valid, and so do the resulting expressions
for the response function $R_\mathbf{q}(t,t_{\rm w})$\eq{Rq} and the full
correlator $C_\mathbf{q}(t,t_{\rm w})$\eqq{Cqtt}{C_twotime}.
The expression\eq{hat_g_general} for the Laplace transform of $g(t)$ that
results from the spherical constraint also still holds, but the
solution is now different. In the $\mathbf{q}$-integral,
the $\mathbf{q}=\mathbf{0}$ contribution $(1/N)C_\mathbf{0}(0,0)/s =
m^2(0)/s$ has to be treated separately. In fact, one sees that for
$s\to 0$ this term always dominates the rest of the integral, which
diverges less strongly. At criticality, where $T=T_{\rm c}=[\int(dq)\,
1/\omega]^{-1}=[2\hat{f}(0)]^{-1}$, one thus has for $s\to 0$
\begin{equation}
2T_{\rm c} \hat{g}(s) = \frac{m^2(0)}{s} \left[\int(dq)\,
\left(\frac{1}{2\omega}-\frac{1}{s+2\omega}\right)\right]^{-1}
\label{gs}
\end{equation}
which using\eq{p_eql_LT} can be rearranged into
\begin{equation}
\frac{\hat{g}(s)}{m^2(0)}=s^{-2}\left[\int (dq)
\frac{T_{\rm c}}{\omega(2\omega+s)}\right]^{-1} = [s^{2}\hat{K}\eql(s)]^{-1}
\label{gs2}
\end{equation}
For $d>4$, $\hat{K}\eql(0)$ is finite so this scales as $s^{-2}$; for
$d<4$, on the other hand, $\hat{K}\eql(s)$ diverges as $s^{(d-4)/2}$
so that $\hat{g}(s) \sim s^{-d/2}$. Translating back to the time
domain, $g(t)$ behaves for large $t$ as
\begin{equation}
g(t)\sim m^2(0)\,t^\alpha, \qquad \alpha = \left\{\begin{array}{lll}
1 & \mbox{ } & (d>4)\\ (d-2)/2 & \mbox{ } & (d<4) \end{array}\right.
\label{g_scaling_magn}
\end{equation}
Note that this asymptotic behaviour is independent of any details of
the initial condition except for the presence of a nonzero $m(0)$; it
depends on the actual value of $m(0)$ only through the prefactor
$m^2(0)$. For the time-dependence of $m(t)$, one gets by taking an
average of\eq{Sqt}
\begin{equation}
m(t) = R_\mathbf{0}(t,0) m(0) = \frac{m(0)}{\sqrt{g(t)}}
\label{m_in_terms_of_g}
\end{equation}
Because of the proportionality of $g(t)$ to $m^2(0)$ for large $t$,
the asymptotic decay of $m(t)$ is independent of the initial
conditions, in terms of {\em both} the decay exponent {\em and} the prefactor.

\subsection{Finite-range spin observables}

We first analyse the correlation and response functions for
observables that relate to a number of spins that is much smaller
than $N$. As for the unmagnetized case, the fluctuations of the
Lagrange multiplier $z$ can then be neglected. To understand the
magnetized case, it is useful to shift the spin variables by $m(t)$.
We will see that the equations of motion then acquire the same form as
before, so that we can directly transfer the main results from the
unmagnetized
case. Explicitly, we consider the following decomposition of the spin
variables
\begin{equation}
S_i=m+U_i
\end{equation}
where $U_i$ is a zero-mean variable. The equation of
motion\eq{real_space_dotSi} for $S_i$ then gives
\begin{equation}
\partial_t m+\partial_t{U_i}=-\Omega_{ij}(m+U_j)+\xi_i-z(m+U_i)
\end{equation}
From\eq{m_in_terms_of_g} and the definition\eq{g_def}, $\partial_t\ln m(t) =
-(1/2)\partial_t\ln g(t) = -z(t)$, so $\partial_t m = -zm$. Also $\sum_j
\Omega_{ij}=0$, giving
\begin{equation}
\partial_t{U_i} = -\Omega_{ij}U_j + \xi_i - zU_i
\label{dot_Ui}
\end{equation}
This is the same as the equation for $S_i$ in the unmagnetized case,
and so one can deduce directly the solutions for the correlation
functions of the $U_i$; these are the connected correlations
$\tilde{C}$. The initial values again become unimportant for long times,
allowing us to work out the scaling of the $\tilde{C}$, and then together
with the response $R$ also the FDR $X$. It is clear from the
description in terms of the subtracted spins $U_i$ that there is
nothing special about the case $\mathbf{q}=\mathbf{0}$, and all results will have a
smooth limit as $\mathbf{q}\to \mathbf{0}$. Because we are neglecting the
fluctuations of the Lagrange multiplier $z$, this limit again has to be
understood as that of a block magnetization calculated over a
region much larger than the correlation length (in the time regime
being considered) but much smaller than the linear system size, so
that $q\gg 1/L$.

Applying\eq{Cqtt}, we can now write down directly the connected
correlation function at equal times as
\begin{equation}
\tilde{C}_\mathbf{q}(t_{\rm w},t_{\rm w})=\tilde{C}_\mathbf{q}(0,0)\frac{e^{-2\omegat_{\rm w}}}{g(t_{\rm w})}
+2T_{\rm c}\frac{e^{-2\omegat_{\rm w}}}{g(t_{\rm w})}\int_0^{t_{\rm w}}dt'e^{2\omega t'}g(t')
\label{Cqtt_magn}
\end{equation}
At long times, the first term is subleading due to the
scaling\eq{g_scaling_magn}, and one has the behaviour
\begin{equation}
\tilde{C}_\mathbf{q}(t_{\rm w},t_{\rm w}) = \frac{T_{\rm c}}{\omega}\sc{C}(\omega t_{\rm w}), \quad
\sc{C}(w) = 2w\int_0^1 dy\,y^{\alpha} e^{-2w(1-y)}
\label{C_scaling_magn}
\end{equation}
This result is of course the same as\eq{Xq_scaling}, except for the
replacement of $-\kappa$ by $\alpha$ which reflects the difference in the
asymptotic behaviour of $g(t)$.

The two-time connected correlations are
$\tilde{C}_\mathbf{q}(t,t_{\rm w})=R_\mathbf{q}(t,t_{\rm w})\tilde{C}_\mathbf{q}(t_{\rm w},t_{\rm w})$ with $R_\mathbf{q}(t,t_{\rm w})$
given by\eq{Rq} as before. As a consequence, the expression\eq{Xq} for
the FDR $X_\mathbf{q}$ also remains valid, and one finds the scaling form
$X_\mathbf{q}(t,t_{\rm w}) = \sc{X}(\omegat_{\rm w})$ with
\begin{equation}
\sc{X}^{-1}(w) = 2-(2w+\alpha)\int_0^1 dy\
y^{\alpha}e^{-2w(1-y)}
\end{equation}
which is directly analogous to\eq{Xq_scaling}. In the limit
$w=\omegat_{\rm w}\rightarrow\infty$, $X_\mathbf{q}=1$, which means that all
modes with $\mathbf{q}\neq 0$ equilibrate once $t_{\rm w}\gg 1/\omega$. In the
opposite limit $w\rightarrow 0$, corresponding to $t_{\rm w}\ll
1/\omega$,
\begin{equation}
X_\mathbf{q}=\frac{\alpha+1}{\alpha+2} = \left\{\begin{array}{lll}
2/3 & \mbox{ } & (d>4) \\
d/(d+2) & \mbox{ } & (d<4)
\end{array}
\right.
\label{Xqneq}
\end{equation}
For $\mathbf{q}\to\mathbf{0}$ this result applies independently of the value of
$t_{\rm w}$ as long as $t_{\rm w}\gg 1$, so that the FDR for the block
magnetization will be a straight line with slope\eq{Xqneq}. This is as
for the unmagnetized case, but the actual value of the FDR is now {\em
different}. It is also different from the value $X^\infty=4/5$
predicted for Ising models in the limit of large
$d$~\cite{GarSolPagRit05}; we will see below that the latter value is
obtained for the {\em global} magnetization, which is affected by
local spin fluctuations of ${{\mathcal O}}(N^{-1/2})$.

For later reference we write down the long-time forms of the
correlation and response functions for $\mathbf{q}\to\mathbf{0}$. By setting
$\omega=0$ in\eq{Cqtt_magn} and taking the long-time limit where the
first term becomes negligible, we find for the connected equal-time
correlator
\begin{equation}
\tilde{C}_\mathbf{0}(t_{\rm w},t_{\rm w}) = \frac{2\Tct_{\rm w}}{1+\alpha}
\end{equation}
The response function is, from\eq{Rq} and\eq{m_in_terms_of_g},
\begin{equation}
R_\mathbf{0}(t,t_{\rm w}) = \sqrt{\frac{g(t_{\rm w})}{g(t)}} = \frac{m(t)}{m(t_{\rm w})}
= \left(\frac{t_{\rm w}}{t}\right)^{\alpha/2}
\label{R_simplification_magn}
\end{equation}
where the last equality holds for long times. The two-time correlator
is therefore
\begin{equation}
\tilde{C}_\mathbf{0}(t,t_{\rm w})=\frac{2\Tct_{\rm w}}{1+\alpha}\left(\frac{t_{\rm w}}{t}\right)^{\alpha/2}
\label{Cc_0}
\end{equation}
From these results one of course retrieves the long-time FDR
$X_\mathbf{0}(t,t_{\rm w})={T_{\rm c} R_\mathbf{0}(t,t_{\rm w})}/{\tilde{C}_\mathbf{0}'(t,t_{\rm w})}=({\alpha+1})/({\alpha+2})$,
obtained in\eq{Xqneq} via the limit $\omegat_{\rm w}\to 0$. As explained
above, these results apply in the regime $1/L\ll q\ll
1$. For $\mathbf{q}=\mathbf{0}$ itself, they capture only the Gaussian part of the
spin-fluctuations, and non-Gaussian corrections become relevant as
discussed in the next section.

\subsection{General expressions for magnetization correlation and response}

We now turn to the FD behaviour of the global magnetization,
corresponding to $\mathbf{q}=\mathbf{0}$ rather than $q\gg {1}/{L}$. All $N$ spins
are now involved and one needs to account for the fluctuating
contribution of the Lagrange multiplier, which we write as
$z+N^{-1/2}\Delta z$ as before. To understand why this is necessary in
the magnetized case, but was not in the unmagnetized scenario,
consider the equation of motion\eq{dotSq} for the zero-wavevector
Fourier component of the spins,
\begin{equation}
\partial_t S_\mathbf{0} = -(z+N^{-1/2}\Delta z)S_\mathbf{0} + \xi_\mathbf{0}
\end{equation}
In the unmagnetized case, $S_\mathbf{0}$ is a zero-mean quantity of
${{\mathcal O}}(N^{1/2})$. The $\Delta z$-term then contributes only subleading
${{\mathcal O}}(1)$ fluctuations. For nonzero magnetization, on the other hand, the
mean of $S_\mathbf{0}$ is $Nm$, with fluctuations around this of
${{\mathcal O}}(N^{1/2})$. The coupling of $\Delta z$ with $m$ then gives an
${{\mathcal O}}(N^{1/2})$ contribution to $\partial_t S_\mathbf{0}$, which is no
longer negligible.

To find the resulting non-Gaussian fluctuations in $S_\mathbf{0}$ explicitly,
we make the decomposition $S_i=s_i+N^{-1/2}r_i$ as before.  The
discussion in Sec.~\ref{sec:setup} then goes through, and we
retrieve\eq{yt3} for the ${{\mathcal O}}(N^{-1/2})$-corrections $r_i$ to the
spins. For the zero Fourier mode, in particular, we have
\begin{equation}
r_\mathbf{0}(t)=-\frac{1}{2} \int dt' dt''R_\mathbf{0}(t,t') s_\mathbf{0}(t')
L(t',t'')\Delta(t'')
\label{r0}
\end{equation}
To simplify the calculation of connected correlations, we now
decompose the Gaussian part of the spins into $s_i = m + u_i$, so
that the $u_i$ are zero-mean Gaussian variables. This corresponds to a
decomposition $U_i=u_i+\sqrt{N}r_i$ of the fluctuating parts of the
spins into leading Gaussian terms and small non-Gaussian corrections,
in analogy to the representation $S_i=s_i+\sqrt{N}r_i$ in the
unmagnetized case. The $u_i$ obey the equation of motion\eq{dot_Ui},
and their correlation and response functions are the $\tilde{C}_\mathbf{q}$ and
$R_\mathbf{q}$ calculated previously.

We will write the connected correlation function for the global
magnetization which {\em includes} non-Gaussian corrections as
$C_{\rm m}(t,t_{\rm w})=(1/N)\left\langle U_\mathbf{0}(t)U_\mathbf{0}(t_{\rm w})\right\rangle$. Making the
decomposition into Gaussian and non-Gaussian parts, this reads
\begin{equation}
C_{\rm m}(t,t_{\rm w})=\frac{1}{N}
\langle[u_\mathbf{0}(t)+N^{-1/2}r_\mathbf{0}(t)][u_\mathbf{0}(t_{\rm w})+N^{-1/2}r_\mathbf{0}(t_{\rm w})]\rangle
\label{concor}
\end{equation}
with, from\eq{r0},
\begin{eqnarray}
r_\mathbf{0}(t) &=& -\frac{1}{2}\int dt'\,dt''\,
R_\mathbf{0}(t,t')[Nm(t')+u_\mathbf{0}(t')]L(t',t'')\Delta(t'') \\
         &=& -\frac{N}{2}\int dt'\,M(t,t')\Delta(t')
\label{r}
\end{eqnarray}
Here we have defined
\begin{equation}
M(t,t_{\rm w})=\int dt'\,R_\mathbf{0}(t,t') m(t') L(t',t_{\rm w}) =
m(t) \int^t dt'\,L(t',t_{\rm w})
\label{Mdef}
\end{equation}
where the second form follows from\eq{R_simplification_magn}. $M$ is,
like $R$ and $L$, causal and so vanishes for $t<t_{\rm w}$. In\eq{r}
we have also discarded the Gaussian fluctuation term $u_\mathbf{0}$, which is
of ${{\mathcal O}}(N^{1/2})$ and so negligible against the term $Nm$. This is
in line with the intuition discussed earlier that non-Gaussian
fluctuations arise only from the coupling of $\Delta z$ to $m$. Note
also that $r_\mathbf{0}$ is ${{\mathcal O}}(N)$, so that in\eq{concor} the
non-Gaussian correction $N^{-1/2}r_\mathbf{0}$ is of the same order as the
Gaussian fluctuation $u_\mathbf{0}$, again as expected.

Substituting\eq{r} into\eq{concor}, we see that we need the two-time
correlations of $u_\mathbf{0}$ and $\Delta$. In the presence of a nonzero $m$
the latter becomes
\begin{equation}
\Delta=\frac{1}{\sqrt{N}}\sum_i(s_i^2-1) = \frac{1}{\sqrt{N}}\sum_i
(m^2-1+2u_i m+u_i^2)
\label{Delta_m_nonzero}
\end{equation}
The required correlations are therefore $\left\langle u_\mathbf{0}(t)u_\mathbf{0}(t_{\rm w})\right\rangle =
N\tilde{C}_\mathbf{0}(t,t_{\rm w})$ and
\begin{eqnarray}
\fl \langle u_\mathbf{0}(t)\Delta(t')\rangle &=& \frac{1}{\sqrt{N}}\sum_{ij}
\langle
u_i(t)[m^2(t')-1+2u_j(t')m(t')+u_j^2(t')]\rangle
\\
\fl
&=& \frac{2}{\sqrt{N}} m(t')
\sum_{ij}\langle u_i(t)u_j(t')\rangle= 2m(t')\sqrt{N}\tilde{C}_\mathbf{0}(t,t')
\end{eqnarray}
For the autocorrelation of $\Delta$, we can exploit the fact
that $\left\langle \Delta\right\rangle = 0$ to write\eq{Delta_m_nonzero} as $\Delta = N^{-1/2}
\sum_i (2u_i m+u_i^2 - \left\langle u_i^2\right\rangle)$. This gives
\begin{eqnarray}
\fl \langle \Delta(t')\Delta (t_{\rm w}')\rangle
&=&
\frac{1}{N}\sum_{ik}\langle[2u_i(t')m(t')+u_i^2(t')-\langle
u_i^2(t')\rangle][\cdots
t'\rightarrow t_{\rm w}' \cdots]\rangle
\\
\fl &=& \frac{4}{N}m(t')m(t_{\rm w}')\sum_{ij}\langle u_i(t')u_j(t_{\rm w}')\rangle
+\frac{2}{N}\sum_{ij}\langle u_i(t')u_j(t_{\rm w}')\rangle^2
\\
\fl &=& 4m(t')m(t_{\rm w}')\tilde{C}_\mathbf{0}(t',t_{\rm w}')+2\int(dq)\,\tilde{C}_\mathbf{q}^2(t',t_{\rm w}')
\end{eqnarray}
where we have used Wick's theorem to simplify the fourth-order average
$\langle (u_i^2-\langle u_i^2\rangle)(u_j^2-\langle
u_j^2\rangle)\rangle = \langle u_i^2 u_j^2\rangle - \langle u_i^2
\rangle \langle u_j^2\rangle = 2\langle u_i
u_j\rangle^2$. Abbreviating the $q$-integral as $\tilde{C}\Cc(t',t_{\rm w}')$, the
full connected correlation function\eq{concor} can thus be written as
\begin{eqnarray}
\fl C_{\rm m}(t,t_{\rm w})&=&
\tilde{C}_\mathbf{0}(t,t_{\rm w})
-\frac{1}{2\sqrt{N}} \int dt'\,M(t,t')\langle u_\mathbf{0}(t_{\rm w})\Delta(t')\rangle
\nonumber
\\
\fl &&{}-{}
\frac{1}{2\sqrt{N}} \int dt'\, M(t_{\rm w},t')
\langle u_\mathbf{0}(t)\Delta(t')\rangle
\nonumber
\\
\fl &&{}+{}\frac{1}{4}\int\,dt'\,dt_{\rm w}'\,
M(t,t') M(t_{\rm w},t_{\rm w}') \langle \Delta(t_{\rm w}')\Delta(t')\rangle
\\
\fl &=&
C_{\rm m}^{(1)}(t,t_{\rm w})+C_{\rm m}^{(2)}(t,t_{\rm w})
\label{DM}
\end{eqnarray}
where
\begin{eqnarray}
\fl C_{\rm m}^{(1)}(t,t_{\rm w})&=&\tilde{C}_\mathbf{0}(t,t_{\rm w})-\int dt'\,
    [M(t,t')\tilde{C}_\mathbf{0}(t_{\rm w},t')+M(t_{\rm w},t')\tilde{C}_\mathbf{0}(t,t')]m(t')
\nonumber\\
\fl & &{}+{}\int\, dt'\,dt_{\rm w}'\, M(t,t')M(t_{\rm w},t_{\rm w}')m(t')m(t_{\rm w}')\tilde{C}_\mathbf{0}(t',t_{\rm w}')
\label{c_m1_unscaled}
\\
\fl &=& \int\, dt'\,dt_{\rm w}'\, [\delta(t-t')-M(t,t')m(t')]
\nonumber\\
\fl & & \times [\delta(t_{\rm w}-t_{\rm w}')-M(t_{\rm w},t_{\rm w}')m(t_{\rm w}')]\tilde{C}_\mathbf{0}(t',t_{\rm w}')
\label{DM_part}
\end{eqnarray}
and
\begin{equation}
\fl C_{\rm m}^{(2)}(t,t_{\rm w})=\frac{1}{2}\int\, dt'\,dt_{\rm w}'\, M(t,t')M(t_{\rm w},t_{\rm w}')\tilde{C}\Cc(t',t_{\rm w}')
\label{c_m2}
\end{equation}

Next we derive an expression for the corresponding magnetization response
function. To this purpose we expand the spins for small fields as
\begin{equation}
S_i=s_i+hr_i
\end{equation}
where $s_i$ are the unperturbed spins and we neglect the
${{\mathcal O}}(N^{-1/2})$ non-Gaussian corrections as irrelevant, as in the
unmagnetized case. The Lagrange multiplier is similarly written as
$z+h\Delta z$. By collecting the ${{\mathcal O}}(h)$ terms from the equation
of motion for the $S_i$, we then find by analogy with\eq{dot_ri} that
the $r_i$ obey
\begin{equation}
\partial_t{r_i}=-\Omega_{ij}r_j-zr_i-\Delta z\, s_i+\delta(t-t_{\rm w})
\label{dot_ri_magn}
\end{equation}
Here the last term represents a field impulse at time $t_{\rm w}$, uniform
over all sites as is appropriate for the field conjugate to the
global magnetization. Since $r_i(t)=0$ before the field is applied,
i.e.\ for $t<t_{\rm w}$, this impulse perturbation gives
$r_i(t_{\rm w}^+)=1$. Starting from this value we can then
integrate\eq{dot_ri_magn} forward in time to get
\begin{equation}
r_i(t)=\sum_j R_{ij}(t,t_{\rm w})-\int_{t_{\rm w}}^t dt'\, R_{ij}(t,t')\Delta z(t')s_j(t')
\label{ri_sol_magn}
\end{equation}
The condition we need to impose in order to get $\Delta z$ is
that the spherical constraint $(1/N)\sum_i(s_i+hr_i)^2=1$ needs to be
satisfied to linear order in $h$, giving the condition
$({1}/{N})\sum_i\langle r_is_i\rangle=0$. Inserting\eq{ri_sol_magn} into
this yields
\begin{equation}
R_\mathbf{0}(t,t_{\rm w})m(t)=\int_{t_{\rm w}}^{t}dt'\,K(t,t')\Delta z(t')
\end{equation}
where we have used the definition\eq{p_def} of $K(t,t')$. Applying the
inverse kernel $L$ gives
\begin{equation}
\Delta z(t)=\int_{t_{\rm w}}^t dt'\,L(t,t')m(t')R_\mathbf{0}(t',t_{\rm w})
\end{equation}
Note that this result vanishes when $m=0$, consistent with the fact
that we did not need to consider perturbations of $z$ in our
calculation of the magnetization response in the unmagnetized case. We
can now write down the magnetization response function, which we
denote by $R_{\rm m}(t,t_{\rm w})$. It is given by $R_{\rm m}=({1}/{N})\sum_i\langle r_i
\rangle$; inserting the result for $\Delta z$ into\eq{ri_sol_magn}, we
get explicitly
\begin{eqnarray}
\fl R_{\rm m}(t,t_{\rm w}) &=& R_\mathbf{0}(t,t_{\rm w})-\int
dt''\,dt'\,R_\mathbf{0}(t,t'')m(t'')L(t'',t')m(t')R_\mathbf{0}(t',t_{\rm w})
\\
\fl &=& \int dt'\, [\delta(t-t')-M(t,t')m(t')]R_\mathbf{0}(t',t_{\rm w})
\label{response_magn}
\end{eqnarray}

This completes the derivation of the general expressions for the
magnetization correlation and response. To make progress, we need to
find the kernel $M(t,t')$. This requires $L(t,t')$, which is the
inverse of $K(t,t')=\int(dq)\, R_\mathbf{q}(t,t')C_\mathbf{q}(t,t')$. As is clear from the
discussion in Sec.~\ref{sec:setup}, the correlator occurring here is
the {\em unsubtracted} one. Because $C_\mathbf{0}(t,t')=Nm(t)m(t')$ is
${{\mathcal O}}(N)$, the $\mathbf{q}=\mathbf{0}$ term needs to be treated separately in
spite of its vanishing weight $1/N$. It makes a contribution
$(1/N)R_\mathbf{0}(t,t')Nm(t)m(t')= m^2(t)\theta(t-t')$, where we have
simplified using\eq{R_simplification_magn}. We can thus write
\begin{equation}
K(t,t') = \tilde{K}(t,t') + m^2(t)\theta(t-t')
\label{K_decomp}
\end{equation}
with the $\mathbf{q}\neq \mathbf{0}$ contribution
\begin{equation}
\tilde{K}(t,t') = \int(dq)\, R_\mathbf{q}(t,t')\tilde{C}_\mathbf{q}(t,t')
\end{equation}
We have switched to the connected correlator here; this makes no
difference for $\mathbf{q}\neq 0$, but allows us to include $\mathbf{q}=\mathbf{0}$ in the
integral because $\tilde{C}_\mathbf{0}={{\mathcal O}}(1)$. To say more, we will need to
distinguish between dimensions $d>4$ and $d<4$.

\subsection{Magnetization correlation and response: Non-equilibrium, $d>4$}

The scaling of the connected part $\tilde{K}$ of $K$ can be analysed exactly
as in the case of zero magnetization: it consists of the same equilibrium
time dependence modulated by an aging function,
$\tilde{K}(t,t')=K\eql(t-t')\sc{\tilde{K}}(t/t')$. The aging part can be worked out
exactly as in\eq{p_scaling} with the only difference arising from the
changed asymptotic behaviour of $g(t)\sim t^{\alpha}$ rather than
$t^{-\kappa}$. For $\sc{\tilde{K}}$ we can therefore use
directly\eq{p_scaling_fn_y_integral}, with $\kappa$ replaced by $-\alpha$:
\begin{eqnarray}
\sc{\tilde{K}}(x)&=&
\frac{d-2}{2}(x-1)^{(d-2)/2}x^{-\alpha} \int_0^1 dy\,
y^{\alpha}(x-y)^{-d/2}
\label{sc_Kc}
\end{eqnarray}
In $d>4$, where $\alpha=1$, the integral can be computed explicitly
to give
\begin{equation}
\sc{\tilde{K}}(x)=1-\frac{d-2}{d-4}\left(\frac{x-1}{x}\right)+
\frac{2}{d-4}\left(\frac{x-1}{x}\right)^{\frac{d-2}{2}}
\end{equation}
We will see below that the precise behaviour of this function does not
affect the results. Briefly though, for $x-1\ll 1$ the second term on the
r.h.s.\ is leading so that $\sc{\tilde{K}}$ decreases linearly with $x-1$,
while for large $x$ one finds by expanding in $1/x$ that $\sc{\tilde{K}} \approx
[(d-2)/4]x^{-2}$.

To find the inverse kernel $L$, consider how $K(t,t')$ varies with
$t$. The first part in\eq{K_decomp} starts off close to unity and
decays on ${{\mathcal O}}(1)$ timescales $t-t'$ as $K\eql(t-t')$, with a
modulation by the aging factor $\sc{\tilde{K}}(t/t')$ once $t-t'$ becomes
comparable to $t'$. The second part, on the other hand, is small
initially but only decays on aging timescales. Comparing
$K\eql(t-t')\sim (t-t')^{(2-d)/2}$ to $m^2(t)\sim 1/t$, this second
term therefore eventually becomes dominant, for $t-t'\sim
t^{2/(d-2)}$.
This discussion
suggests that also the inverse kernel $L$ should be composed of two
parts with distinct long-time behaviour. We therefore write
\begin{equation}
L = \tilde{L}+L_\zv
\end{equation}
where $\tilde{L}$ is the inverse of $\tilde{K}$ and $L_\zv$ arises from the
zero-wavevector part of $K$. The continuous part of $L$ is then
similarly decomposed as $L^{(2)} = \tilde{L}^{(2)} + \Ltwo_\zv$.

We proceed by writing the defining equations for $L^{(2)}$ and
$\tilde{L}^{(2)}$. The full inverse $L$ is defined by\eq{pinv_def} and as before has
singular parts which are related to the behaviour of $K(t,t')$ for
$t\to t'$. One can show directly from the definition of $K$, and
exactly as in the unmagnetized case, that
\begin{equation}
K(t'^+,t') = 1, \qquad
\left.\partial_{t'} K(t,t')\right|_{t\to t'^+} = 2T_{\rm c}
\label{K_initial}
\end{equation}
The decomposition\eq{pinv_structure} of the inverse kernel $L$
therefore also remains valid, and from\eq{pinv_def} and\eq{K_decomp}
we get the following equation for its continuous part $L^{(2)}$
\begin{equation}
\fl \int dt'\,[\tilde{K}(t,t')+m^2(t)]L^{(2)}(t',t_{\rm w})
=2T_{\rm c}\tilde{K}(t,t_{\rm w})+2T_{\rm c} m^2(t) -\partial_{t_{\rm w}} \tilde{K}(t,t_{\rm w})
\label{N.6}
\end{equation}
This is the analogue of the relation\eq{ptwo_def} for the case $m=0$. We can
argue similarly for $\tilde{L}$, which is defined by
\begin{equation}
\int dt'\,\tilde{K}(t,t')\tilde{L}(t',t_{\rm w})=\delta(t-t_{\rm w})
\label{defi}
\end{equation}
From\eq{K_decomp} and\eq{K_initial},
\begin{equation}
\tilde{K}(t'^+,t') = 1-m^2(t), \qquad
\left.\partial_{t'}\tilde{K}(t,t')\right|_{t\to t'^+} = 2T_{\rm c}
\end{equation}
and this initial behaviour implies that $\tilde{L}$ can be decomposed as
\begin{equation}
\tilde{L}(t',t_{\rm w})=\frac{\delta'(t'-t_{\rm w})}{1-m^2(t_{\rm w})}
+\frac{2T_{\rm c}}{[1-m^2(t_{\rm w})]^2}\delta(t'-t_{\rm w})-\tilde{L}^{(2)}(t',t_{\rm w})
\label{lc}
\end{equation}
Inserting into the definition\eq{defi} gives for the continuous part
$\tilde{L}^{(2)}$
\begin{equation}
\fl \int dt'\,\tilde{K}(t,t')\tilde{L}^{(2)}(t',t_{\rm w})=
\frac{2T_{\rm c}}{[1-m^2(t_{\rm w})]^2} \tilde{K}(t,t_{\rm w})
-\frac{1}{1-m^2(t_{\rm w})}\partial_{t_{\rm w}}\tilde{K}(t,t_{\rm w})
\label{lc2}
\end{equation}
Now for long times, we can approximate $1-m^2(t_{\rm w})\approx
1$. Then\eq{lc2} becomes identical to the relation\eq{ptwo_def} which
determined $L^{(2)}$ in the unmagnetized case. Since $\tilde{K}$ has the same
scaling form as $K$ in\eq{ptwo_def}, except for the replacement of
$\sc{K}$ by $\sc{\tilde{K}}$, the solution for $\tilde{L}^{(2)}$ can be found in
exactly the same way. In particular, the scaling functions describing
the aging corrections in $\tilde{K}$ and $\tilde{L}^{(2)}$ are again identical, and
we can write directly
\begin{equation}
\tilde{L}^{(2)}(t,t')=L^{(2)}\eql(t-t')\sc{\tilde{K}}(t/t')
\end{equation}
as the long-time form of $\tilde{L}^{(2)}$. Here $L^{(2)}\eql$ is the same function
as in the unmagnetized case, with Laplace transform\eq{ptwo_eql}.

It now remains to find $\Ltwo_\zv$. Subtracting\eq{lc2} from\eq{N.6} gives
\begin{eqnarray}
\fl\lefteqn{
\int_{t_{\rm w}}^t dt'\,\tilde{K}(t,t')\Ltwo_\zv(t',t_{\rm w})+m^2(t)\int_{t_{\rm w}}^t
dt'\,\Ltwo_\zv(t',t_{\rm w}) =
2T_{\rm c} \frac{-2m^2(t_{\rm w})+m^4(t_{\rm w})}{[1-m^2(t_{\rm w})]^2} \tilde{K}(t,t_{\rm w}) }
\nonumber
\\
& &
{}+{}m^2(t)\left[2T_{\rm c}-\int_{t_{\rm w}}^t dt'\,\tilde{L}^{(2)}(t',t_{\rm w})\right]
+\frac{m^2(t_{\rm w})}{1-m^2(t_{\rm w})}\partial_{t_{\rm w}}\tilde{K}(t,t_{\rm w})
\label{deltal2}
\end{eqnarray}
To make progress we assume that, by analogy with the
structure\eq{K_decomp} of $K$, $\Ltwo_\zv$ varies only on aging
timescales; we will find this confirmed {\em a posteriori}. We can
then concentrate on the aging regime $t-t_{\rm w}\simt_{\rm w}\gg 1$. In this
regime, the integral $\int_{t_{\rm w}}^t dt'\,\tilde{L}^{(2)}(t',t_{\rm w}) = \int_{t_{\rm w}}^t dt'\,
L^{(2)}\eql(t'-t_{\rm w}) \sc{\tilde{K}}(t'/t_{\rm w})$ on the r.h.s.\ of\eq{deltal2}
becomes to leading order $\int_{t_{\rm w}}^\infty
dt'\,L^{(2)}\eql(t'-t_{\rm w})=\hatL^{(2)}\eql(0)$; the aging correction
$\sc{\tilde{K}}$ is unimportant because the integral converges for
$t'-t_{\rm w}={{\mathcal O}}(1)$, and the upper integration limit can be set to
infinity for the same reason. From\eq{ptwo_eql},
$\hatL^{(2)}\eql(0)=2T_{\rm c}-1/\hatK\eql(0)$, and so the square bracket on
the r.h.s.\ of\eq{deltal2} becomes simply the constant
$1/\hatK\eql(0)$. In the first and third term on the r.h.s., on the
other hand, $\tilde{K}$ and $\partial_{t_{\rm w}}\tilde{K}$ scale as
$(t-t_{\rm w})^{-(d-2)/2}$ and $(t-t_{\rm w})^{-d/2}$, respectively, and are
negligible compared to the second term in the aging regime.
Disregarding these subleading terms, equation\eq{deltal2}
is transformed to
\begin{equation}
\int_{t_{\rm w}}^t
dt'\,\tilde{K}(t,t')\Ltwo_\zv(t',t_{\rm w})+m^2(t)\int_{t_{\rm w}}^t dt'\,\Ltwo_\zv(t',t_{\rm w})
=\frac{m^2(t)}{\hat{K}\eql(0)}
\end{equation}
In the first integral, if as assumed $\Ltwo_\zv$ varies only on aging
timescales, the integral is dominated by the region $t'\approx t$
because of the factor $K\eql(t-t')$ in $\tilde{K}(t,t')$. This factor again
makes the integral convergent within a region $t-t'={{\mathcal O}}(1)$, and so
we can approximate it by $\hat{K}\eql(0)\Ltwo_\zv(t,t_{\rm w})$. This gives
\begin{equation}
\int_{t_{\rm w}}^{t} dt'\,\Ltwo_\zv(t',t_{\rm w}) =
\frac{1}{\hat{K}\eql(0)}-
\frac{\hat{K}\eql(0)}{m^2(t)}\Ltwo_\zv(t,t_{\rm w})
\label{dl2}
\end{equation}
and specifically in the limit $t/t_{\rm w}\to 1$
\begin{equation}
\Ltwo_\zv(t_{\rm w},t_{\rm w})=\hat{K}\eql^{-2}(0)m^2(t_{\rm w})
\label{dLtwo_initial}
\end{equation}
To find $\Ltwo_\zv$ for $t/t_{\rm w}>1$, we use that $m^2(t)=\mu_d/t$ for large
$t$ and in $d>4$, see\eqq{g_scaling_magn}{m_in_terms_of_g},
with some dimension-dependent coefficient $\mu_d$.
Taking a derivative of\eq{dl2} with respect to $t$ then gives
\begin{equation}
\Ltwo_\zv(t,t_{\rm w})\left[1+\mu_d^{-1}\hat{K}\eql(0)\right]=
-t {\mu_d}^{-1}\hat{K}\eql(0)\partial_t\Ltwo_\zv(t,t_{\rm w})
\end{equation}
This implies $\partial(\ln\Ltwo_\zv)/\partial(\ln t) =
-[1+\mu_d/\hat{K}\eql(0)]$ and so together with\eq{dLtwo_initial}
\begin{equation}
\Ltwo_\zv(t,t_{\rm w}) = \hat{K}\eql^{-2}(0)
\frac{\mu_d}{t_{\rm w}}\left(\frac{t}{t_{\rm w}}\right)^{-[1+\mu_d/\hat{K}\eql(0)]}
\label{incremento}
\end{equation}
This can be simplified because in fact $\mu_d=\hat{K}\eql(0)$. To see
this, note from\eq{gs2} that $\hat{g}(s)/m^2(0)=1/[\hat{K}\eql(0)s^2]$
for small $s$; here we have used that $\hat{K}\eql(0)$ is finite for
$d>4$. Transforming back to the time domain gives
$g(t)/m^2(0)=m^{-2}(t)=t/\hat{K}\eql(0)$ for large $t$, i.e.\
$m^2(t)=\hat{K}\eql(0)/t$.
Thus\eq{incremento} simplifies to
\begin{equation}
\Ltwo_\zv(t,t_{\rm w}) = \frac{t_{\rm w}}{\mu_d t^2}
\label{dLtwo}
\end{equation}
This result is consistent with our assumption that $\Ltwo_\zv$
only varies on aging timescales. Overall, we have thus found for
$L^{(2)}(t,t_{\rm w})$ the following long-time form
\begin{equation}
L^{(2)}(t,t_{\rm w})=L^{(2)}\eql(t-t_{\rm w})\sc{\tilde{K}}\left({t}/{t_{\rm w}}\right)
+\frac{t_{\rm w}}{\mu_d t^2}
\label{Ltwo_high_d}
\end{equation}
Of course, for time differences $t-t_{\rm w}={{\mathcal O}}(1)$, $\Ltwo_\zv(t,t_{\rm w})$
will deviate from the form\eq{dLtwo} derived for the aging
regime. However, one can verify by expanding both sides of\eq{deltal2} to
${{\mathcal O}}(t-t_{\rm w})$ that, even for $t=t_{\rm w}$, $\Ltwo_\zv$ remains of order $1/t_{\rm w}$,
so that these small deviations are always negligible compared to the
first term in\eq{Ltwo_high_d}.

We can now proceed to find $M$ as defined in\eq{Mdef}, and from there
finally the explicit forms of the magnetization correlation and
response functions. Using the general structure\eq{pinv_structure} of
$L$ we have
\begin{equation}
m^{-1}(t)M(t,t_{\rm w})=\delta(t-t_{\rm w})+2T_{\rm c}-\int_{t_{\rm w}}^{t}dt'\, L^{(2)}(t',t_{\rm w})
\label{M_general}
\end{equation}
The integral can be separated into the contributions from the two
parts of $L^{(2)}$ as given in\eq{Ltwo_high_d}. The first part yields an
integral that converges for $t'-t_{\rm w}={{\mathcal O}}(1)$; for $t-t_{\rm w}\gg 1$, it
therefore gives $\hat{L}\eql(0)=2T_{\rm c}-1/\hat{K}\eql(0)=2T_{\rm c}-1/\mu_d$ to leading
order. The second part, on the other hand, yields explicitly
$\int_{t_{\rm w}}^{t}dt' (t_{\rm w}/\mu_d t'^2) = \mu_d^{-1}(1-t_{\rm w}/t)$, so that
\begin{equation}
\fl m^{-1}(t)M(t,t_{\rm w})=\delta(t-t_{\rm w})+\frac{1}{\mu_d} -
\frac{1}{\mu_d}\left(1-\frac{t_{\rm w}}{t}\right) = \delta(t-t_{\rm w}) +
\frac{1}{\mu_d}\frac{t_{\rm w}}{t}
\label{M}
\end{equation}
This result applies for $t-t_{\rm w}\gg 1$. For $t-t_{\rm w}={{\mathcal O}}(1)$ it is not
accurate; e.g.\ at $t=t_{\rm w}$ the continuous part of $m^{-1}(t)M(t,t_{\rm w})$
is, from\eq{M_general}, $2T_{\rm c}$ rather than $1/\mu_d$. However, this
deviation over an ${{\mathcal O}}(1)$ time-range only gives subleading corrections in
the integrals over $M$ that we need below, as indeed does the
$\delta(t-t_{\rm w})$-term.
This can be seen in\eq{DM_part} and\eq{response_magn}, where only the
combination $\delta(t-t_{\rm w})-M(t,t_{\rm w})m(t_{\rm w})$ occurs; the latter can be written as
\begin{equation}
\fl\delta(t-t_{\rm w})-M(t,t_{\rm w})m(t_{\rm w}) = \delta(t-t_{\rm w})[1-m(t)m(t_{\rm w})] -
\frac{t_{\rm w}^{1/2}}{t^{3/2}} = \delta(t-t_{\rm w}) - \frac{t_{\rm w}^{1/2}}{t^{3/2}}
\label{1_minus_Mm}
\end{equation}
The second form holds for long times, where the $m(t)m(t_{\rm w})$ term that
originated from the $\delta$-term in\eq{M} is negligible.

We can now work out the expression\eq{DM} for the full connected
correlation function. One can show that the contribution $C_{\rm m}^{(2)}$
involving $\tilde{C}\Cc$ is negligible in the long-time limit.
In the expression\eq{c_m1_unscaled} for the remainder $C_{\rm m}^{(1)}$,
let us call the second and third term $I_2$ and $I_3$.
We need $\tilde{C}_\mathbf{0}$, which
from\eq{Cc_0} reads $\tilde{C}_\mathbf{0}(t',t_{\rm w}')=T_{\rm c} t_{\rm w}'^{3/2}t'^{-1/2}$ for $t'>t_{\rm w}'$;
because $\tilde{C}_\mathbf{0}$ is symmetric in time this then implies
$\tilde{C}_\mathbf{0}(t',t_{\rm w}')=T_{\rm c} t'^{3/2}t_{\rm w}'^{-1/2}$ for $t'<t_{\rm w}'$. Paying
due attention to this temporal ordering of the arguments of $\tilde{C}_\mathbf{0}$ and
using\eq{1_minus_Mm}, one finds
\begin{eqnarray}
\fl I_2(t,t_{\rm w}) &=& T_{\rm c}\left[\int_0^{t_{\rm w}}dt'\, \frac{t'^{1/2}}{t^{3/2}}
\frac{t'^{3/2}}{t_{\rm w}^{1/2}}+\int_{t_{\rm w}}^{t}dt'\,
\frac{t'^{1/2}}{t^{3/2}}\frac{t_{\rm w}^{3/2}}{t'^{1/2}}
\right] +  T_{\rm c} \int_0^{t_{\rm w}}dt'\, \frac{t'^{1/2}}{t_{\rm w}^{3/2}}
\frac{t'^{3/2}}{t^{1/2}}
\\
\fl &=& T_{\rm c}\left(\frac{t_{\rm w}}{t}\right)^{3/2}\left(\frac{4}{3}t-\frac{2}{3}t_{\rm w}\right)
\end{eqnarray}
Similarly, the double integral in\eq{c_m1_unscaled} can be evaluated as
\begin{eqnarray}
\fl I_3(t,t_{\rm w})&=&\int_0^t dt'\int_0^{t_{\rm w}}\!\!dt_{\rm w}'\,
\frac{t'^{1/2}}{t^{3/2}}\frac{t_{\rm w}'^{1/2}}{t_{\rm w}^{3/2}}\,\tilde{C}_\mathbf{0}(t_{\rm w}',t')
\\
\fl &=&T_{\rm c}\int_{t_{\rm w}}^t dt'\int_0^{t_{\rm w}}\!\!dt_{\rm w}'\, \frac{t'^{1/2}}{t^{3/2}}
\frac{t_{\rm w}'^{1/2}}{t_{\rm w}^{3/2}}\frac{t_{\rm w}'^{3/2}}{t'^{1/2}}
+2T_{\rm c}\int_0^{t_{\rm w}}\!\! dt'\int_0^{t'}\!\!dt_{\rm w}'\, \frac{t'^{1/2}}{t^{3/2}}
\frac{t_{\rm w}'^{1/2}}{t_{\rm w}^{3/2}}\frac{t_{\rm w}'^{3/2}}{t'^{1/2}}
\\
\fl &=&\frac{T_{\rm c}}{3}\left(\frac{t_{\rm w}}{t}\right)^{3/2}(t-t_{\rm w})
+ \frac{2T_{\rm c}}{3}\int_0^{t_{\rm w}} dt'\frac{t'^3}{t^{3/2}t_{\rm w}^{3/2}}
\\
\fl &=&\frac{T_{\rm c}}{6}\left(\frac{t_{\rm w}}{t}\right)^{3/2}(2t-t_{\rm w})
\end{eqnarray}
Our final long-time result for the connected magnetization correlator
including non-Gaussian corrections is then
\begin{eqnarray}
\fl C_{\rm m}(t,t_{\rm w})&=&
T_{\rm c}\frac{t_{\rm w}^{3/2}}{t^{1/2}}+I_2(t,t_{\rm w})+I_3(t,t_{\rm w})\\
\fl &=&
T_{\rm c}\left(\frac{t_{\rm w}}{t}\right)^{3/2}\left(t-\frac{4}{3}t+\frac{2}{3}t_{\rm w}
+\frac{1}{3}t-\frac{1}{6}t_{\rm w}\right)
\ = \ \frac{\Tct_{\rm w}}{2}
\left(\frac{t_{\rm w}}{t}\right)^{3/2}
\label{Cm_dgt4}
\end{eqnarray}
For the conjugate magnetization response function we get
from\eq{response_magn} and\eq{M}, together with
$R_\mathbf{0}(t,t_{\rm w})=(t_{\rm w}/t)^{1/2}$,
\begin{eqnarray}
\fl R_{\rm m}(t,t_{\rm w}) &=&
R_\mathbf{0}(t,t_{\rm w})-\int_{t_{\rm w}}^t dt'\,
\frac{t'^{1/2}}{t^{3/2}}R_\mathbf{0}(t',t_{\rm w})
\\
\fl &=&\left(\frac{t_{\rm w}}{t}\right)^{1/2}-\frac{t_{\rm w}^{1/2}}{t^{3/2}}(t-t_{\rm w})
=\left(\frac{t_{\rm w}}{t}\right)^{3/2}
\label{Rm_dgt4}
\end{eqnarray}
The FDR follows finally as
\begin{equation}
X_{\rm m}(t,t_{\rm w})=\frac{T_{\rm c}R_{\rm m}(t,t_{\rm w})}{C_{\rm m}'(t,t_{\rm w})}=\frac{4}{5}
\label{X_m_d_gt_4}
\end{equation}
Interestingly, this agrees exactly with the result for the Ising
ferromagnet in the limit of large dimensionality
$d$~\cite{GarSolPagRit05}. As in the unmagnetized case, we see
therefore that it is the {\em global} observables in the spherical
model, which are strongly affected by non-Gaussian fluctuations, that
behave like their analogues in short-range models. In fact, even the
expressions for the correlation and response functions we find here
are identical to those in the large-$d$ Ising case, implying
``universality'' at a more detailed level than one might have expected.

It is worth noting
that the effect of the non-Gaussian corrections is very large:
compared to the Gaussian result
$\tilde{C}_\mathbf{0}(t,t_{\rm w})\sim t_{\rm w}(t_{\rm w}/t)^{1/2}$, the corrections increase the
decay exponent
to $C_{\rm m}(t,t_{\rm w})\simt_{\rm w}(t_{\rm w}/t)^{3/2}$, so that $C_{\rm m}/\tilde{C}_\mathbf{0} \sim t_{\rm w}/t
\ll 1 $ for $t\gg t_{\rm w}$: there is an almost perfect cancellation of
Gaussian terms and non-Gaussian corrections for well-separated
times. Similar comments apply to the response. The
overall effect of the non-Gaussian corrections on the FD relation is
to leave this as a straight line (since $X$ is constant), but to
increase the slope from 2/3 to 4/5.

\subsection{Magnetization correlation and response: Non-equilibrium, $d<4$}

We now consider systems below the upper critical dimension, $d<4$;
here there are no predictions yet from other models for the
non-equilibrium FD behaviour following a quench of a magnetized
initial state to $T_{\rm c}$ (but see Sec.~\ref{sec:conclusion}). As in the
case of the energy correlations for
unmagnetized initial states, leading order cancellation effects have
to be taken care of for these low values of $d$.

We again need to know the scaling of $K$ to determine $L^{(2)}$; from
this we then get $M$ and finally the correlation and response
functions. The long-time scaling of the connected part of $K$ is
$\tilde{K}(t,t')=K\eql(t-t')\sc{\tilde{K}}({t}/{t'})$ as before, with $\sc{\tilde{K}}$
given by \eq{sc_Kc} but now $\alpha=(d-2)/2$.  The contribution to $K$
from $\mathbf{q}=\mathbf{0}$, given by second term in\eq{K_decomp}, is negligible
relative to $\tilde{K}(t,t')$ for $t-t'={{\mathcal O}}(1)$. However, for $t-t'\gg
1$, it becomes comparable and has the same overall time scaling as
$\tilde{K}(t,t')$ in the aging regime. To see this explicitly, recall
from\eqq{g_scaling_magn}{m_in_terms_of_g} that
the square of the magnetization decays asymptotically as $m^2(t)=\mu_d
t^{-\alpha}=\mu_d t^{(2-d)/2}$ with some constant $\mu_d$. Similarly,
the equilibrium part of $\tilde{K}$ behaves as $K\eql(t-t') =
k_d(t-t')^{(2-d)/2}$ for $t-t'\gg 1$ (see after\eq{p_eql}). This gives
\begin{equation}
\frac{m^2(t)}{K\eql(t-t')}
=
\frac{\mu_d(t-t')^{(d-2)/2}}{k_dt^{(d-2)/2}}=
\frac{\mu_d}{k_d}\left(\frac{t/t'-1}{t/t'}\right)^{(d-2)/2}
\label{scaling}
\end{equation}
in the aging regime where both $t'$ and $t-t'$ are large.
(For $t-t'={{\mathcal O}}(1)$ this
expression is not accurate but this is irrelevant because there the
term $m^2(t)$ is subleading compared to $\tilde{K}(t,t')$ anyway.)
We thus have the overall scaling of $K$
\begin{equation}
\fl K(t,t')=K\eql(t-t')\sc{K}(t/t'),
\quad
\sc{K}(x) = \sc{\tilde{K}}(x)+\frac{\mu_d}{k_d}
\left(\frac{x-1}{x}\right)^{(d-2)/2}
\label{Ktot}
\end{equation}
To simplify the scaling function, we integrate the
expression\eq{sc_Kc} by parts and rescale $y\to xy$, bearing in mind
that $\alpha=(d-2)/2$:
\begin{eqnarray}
\sc{\tilde{K}}(x)&=&\frac{d-2}{2} \left(\frac{x-1}{x}\right)^{(d-2)/2}
\left[\frac{2}{d-2}(x-1)^{(2-d)/2}\right.
\nonumber\\
& & {}-{}\left.\int_0^{1/x}dy\,y^{(d-4)/2}(1-y)^{(2-d)/2}\right]
\end{eqnarray}
For $x\to 1$ the integral becomes a Beta-function which evaluates to
$\Gamma((d-2)/2) \Gamma((4-d)/2)$, and extracting this term gives
\begin{eqnarray}
\fl \sc{\tilde{K}}(x)&=& x^{(2-d)/2}-\Gamma\left({\textstyle\frac{d}{2}}\right)
\Gamma\left({\textstyle\frac{4-d}{2}}\right) \left(\frac{x-1}{x}\right)^{(d-2)/2}
\nonumber\\
\fl & & {}+{} \frac{d-2}{2}
\left(\frac{x-1}{x}\right)^{(d-2)/2}\int_{1/x}^1 dy\,
y^{(d-4)/2}(1-y)^{(2-d)/2}
\end{eqnarray}
The second term has the same dependence on $x$ as the additional
contribution from zero wavevector in\eq{Ktot}. In fact, it turns out
that these two terms cancel exactly: from our
definitions of $k_d$ and $\mu_d$ we have $K\eql(t)=k_dt^{(2-d)/2}$ for
large $t$, while $m^{-2}(t)=g(t)/m^2(0)=\mu_d^{-1}t^{(d-2)/2}$. Laplace
transforming gives $\hat{K}\eql(s)=k_d \Gamma((4-d)/2) s^{(d-4)/2}$
and $\hat{g}(s)/m^2(0) = \mu_d^{-1}\Gamma(d/2)s^{-d/2}$ to leading
order for small $s$. But then\eq{gs2} shows that
$\mu_d^{-1}\Gamma(d/2)=[k_d\Gamma((4-d)/2)]^{-1}$, or $\mu_d/k_d =
\Gamma(d/2)\Gamma((4-d)/2)$, proving the cancellation anticipated
above. Overall, we thus have for the scaling function of $K$
\begin{equation}
\fl\sc{K}(x)=x^{(2-d)/2}+\frac{d-2}{2}\left(\frac{x-1}{x}\right)^{(d-2)/2}
\int_{1/x}^1 dy\,y^{(d-4)/2}(1-y)^{(2-d)/2}
\label{Ktot_final}
\end{equation}
Expanding for $x\approx 1$, one sees that the leading order variation
is linear in $x-1$:
\begin{eqnarray}
\fl \sc{K}(x)&\approx&1+\frac{2-d}{2}(x-1)+
\frac{d-2}{2}\left(\frac{x-1}{x}\right)^{(d-2)/2}
\int_{1/x}^1dy\, (1-y)^{(2-d)/2}
\nonumber\\
\fl &\approx&1+\frac{1}{2}\frac{(d-2)^2}{4-d}(x-1)
\end{eqnarray}
Note that the prefactor is positive, so that $\sc{K}(x)$ {\em
increases} with $x$ in the current scenario. This trend persists for
all $x$, not just $x\approx 1$, and the scaling function monotonically
approaches a limit value for $x\to\infty$. The latter
follows from\eq{Ktot} as $\mu_d/k_d=\Gamma(d/2)\Gamma((4-d)/2)$,
since the connected contribution $\sc{\tilde{K}}(x)$ decays to zero for
$x\to\infty$.

With the single overall scaling\eq{Ktot} of $K$ we no longer need to
decompose $L^{(2)}$ into $\tilde{L}^{(2)}$ and $\Ltwo_\zv$ as we did in $d>4$;
instead the long-time behaviour of $L^{(2)}$ will have the same
structure as in the unmagnetized case,
\begin{equation}
L^{(2)}(t,t_{\rm w})=L^{(2)}\eql(t-t_{\rm w})\sc{L}(t/t_{\rm w})
\label{Ltot}
\end{equation}
One can then follow exactly the discussion in Sec.~\ref{sec:KandL} to
arrive at the integral equation\eq{pinv_cond}
for $\sc{L}(x)$. Solving the latter
looks a rather formidable task, given that $\sc{K}$ itself has the
relatively complicated form\eq{Ktot_final}. Remarkably, however,
the solution can be found in closed form and is given simply by
\begin{equation}
\sc{L}(x)=\frac{2}{4-d}\,x^{(2-d)/2}+\frac{2-d}{4-d}\,x^{-d/2}
\label{solution_L}
\end{equation}
We were led to this result initially by a systematic series expansion
of both $\sc{K}(x)$ and $\sc{L}(x)$ in terms of $(x-1)/x$. We do not
detail this here, but verify in~\ref{sec:L_solution} by
direct calculation that\eq{solution_L} does indeed
solve\eq{pinv_cond}.

We next calculate the kernel $M(t,t_{\rm w})$. One inserts the
scaling\eq{Ltot} into\eq{M_general}, subtracting off and adding back
on the contribution from the equilibrium part of $L^{(2)}$:
\begin{eqnarray}
m^{-1}(t)M(t,t_{\rm w})&=&\delta(t-t_{\rm w})+2T_{\rm c}-\int_{t_{\rm w}}^t
dt'\,L^{(2)}\eql(t'-t_{\rm w})
\nonumber\\
& &{}+{}\int_{t_{\rm w}}^{t} dt'\,L^{(2)}\eql(t'-t_{\rm w})
\left[1-\sc{L}\left(\frac{t'}{t_{\rm w}}\right)\right]
\label{Mt}
\end{eqnarray}
This is done to account for the leading order cancellation of the
second and third terms:
\begin{equation}
\fl 2T_{\rm c}-\int_{t_{\rm w}}^t\!\!
dt'\,L^{(2)}\eql(t'-t_{\rm w}) =
2T_{\rm c}-\hatL^{(2)}\eql(0)+\int_t^\infty \!\!dt'\,L^{(2)}\eql(t'-t_{\rm w}) =
\frac{2\lambda_d}{4-d}
(t-t_{\rm w})^{({d-4})/{2}}
\end{equation}
where in the last step we have restricted ourselves to the aging
regime $t-t_{\rm w}\gg 1$ and used the asymptotic
behaviour\eq{ptwo_eql_asymptotics}, $L^{(2)}\eql(t'-t_{\rm w})=\lambda_d
(t'-t_{\rm w})^{(d-6)/2}$. In the remaining integral in\eq{Mt}, the factor
$1-\sc{L}(t'/t_{\rm w})\sim (t'-t_{\rm w})/t_{\rm w}$ ensures that the integration no
longer has its weight concentrated near $t'=t_{\rm w}$; put differently,
after scaling the integration variable by $t_{\rm w}$ to $y=t'/t_{\rm w}$ the
integral is convergent at the lower end. We thus obtain in the aging
regime
\begin{eqnarray}
\fl m^{-1}(t)M(t,t_{\rm w})&=&\delta(t-t_{\rm w}) + \lambda_dt_{\rm w}^{(d-4)/2}
\left\{\frac{2}{4-d}(x-1)^{(d-4)/2} \right.
\nonumber\\
& & {}+{} \left. \int_1^x dy\,(y-1)^{(d-6)/2}[1-\sc{L}(y)]\right\}
\label{initial_M}
\end{eqnarray}
Inserting\eq{solution_L} for $\sc{L}$, the integral in\eq{initial_M}
can be done explicitly to give
$[2/(4-d)]x^{(2-d)/2}(x-1)^{(d-4)/2}$ for the sum of the terms in
curly brackets. A little care is needed
here because the separate integrals over $(y-1)^{(d-6)/2}$ and
$(y-1)^{(d-6)/2}\sc{L}(y)$ are divergent at the lower end.
One can avoid this by analytical continuation from $d>4$, where these
divergences are absent, or by integrating from $1+\epsilon$ to $x$ and
taking $\epsilon\to 0$ at the
end.
The $\delta$-term is again subleading for long times in the relevant
combination\eq{1_minus_Mm}, which we can write as
\begin{eqnarray}
\fl \delta(t-t_{\rm w})-M(t,t_{\rm w})m(t_{\rm w}) &=& \delta(t-t_{\rm w}) -
 \frac{2}{4-d}\frac{\lambda_d\mu_dt_{\rm w}^{(d-4)/2}}{t_{\rm w}^{(d-2)/4}t^{(d-2)/4}}x^{(2-d)/2}(x-1)^{(d-4)/2}
\nonumber\\
\fl &=& \delta(t-t_{\rm w})-\frac{1}{t}\sc{M}(x)
\label{combination}
\end{eqnarray}
with
\begin{equation}
\sc{M}(x) = \frac{d-2}{2}\, x^{(2-d)/4}\left(\frac{x-1}{x}\right)^{(d-4)/2}
\label{M_scaling_final}
\end{equation}
Here we have eliminated the constants $\lambda_d$ and
$\mu_d$ using the following argument. From\eq{gs2},
$\hat{g}(s)/m^2(0)=\hat{L}\eql(s)/s^2$ for small
$s$, i.e.\ $\hat{L}\eql(s)=s^2\hat{g}(s)/m^2(0)$. In the time domain
this gives at long times
$L^{(2)}\eql(t)=-L\eql(t)=-(\partial_t)^2m^{-2}(t)=-\mu_d^{-1}(\partial_t)^2
t^{(d-2)/2} = -\mu_d^{-1}[(d-2)/2][(d-4)/2]t^{(d-6)/2}$, so that
$\lambda_d\mu_d=(d-2)(4-d)/4$.

Reassuringly, for $d\to 4$ the result\eq{M_scaling_final} tends to
$\sc{M}(x)=x^{-1/2}$,
matching smoothly onto the result\eq{1_minus_Mm} we found earlier in $d>4$.
For $d<4$, the scaling function diverges as $x\to 1$. Since
from\eq{M_general} the continuous part of $M(t,t_{\rm w})$ is exactly given
by $2T_{\rm c} m(t)$ for equal times, this indicates that the above aging
regime expression must break down eventually when $t-t_{\rm w}$ becomes
small, as expected. In the integrals where $M$ appears below such
effects can be neglected, however, because they only give subleading
corrections.

With the scaling of $\sc{M}$ in hand we can now compute the connected
correlation function $C_{\rm m}(t,t_{\rm w})$ including non-Gaussian corrections,
given by\eq{DM}. After rescaling the integration variables
the first part\eq{DM_part} can be written as
\begin{eqnarray}
\fl C_{\rm m}^{(1)}(t,t_{\rm w})&=& \int_0^x dy\,
[\delta(y-x)-x^{-1}\sc{M}(x/y)]
\nonumber\\
\fl & & \times \int_0^1 dy_{\rm w}\,
[\delta(y_{\rm w}-1)-\sc{M}(1/y_{\rm w})] \tilde{C}_\mathbf{0}(t_{\rm w} y,t_{\rm w} y')
\label{Ione}
\end{eqnarray}
where the Gaussian magnetization correlator is
\begin{equation}
\fl \tilde{C}_\mathbf{0}(t_{\rm w} y,t_{\rm w} y')
= ({4\Tct_{\rm w}}/{d})\min\left\{y(y/y_{\rm w})^{(d-2)/4},y_{\rm w}(y_{\rm w}/y)^{(d-2)/4}\right\}
\end{equation}
from\eq{Cc_0}.
At first sight\eq{Ione} suggests, e.g.\ from the $\delta(y-x)$ term,
an asymptotic decay of $C_{\rm m}^{(1)}\sim t_{\rm w} x^{(2-d)/4}$
for large $x$. But this would not match continuously with the
result\eq{Cm_dgt4} we found for $d>4$. A cancellation of such leading
order terms must therefore occur for $x\to\infty$. To show this
explicitly, we verify from\eq{M_scaling_final} the identity
\begin{eqnarray}
\fl\lefteqn{\int_0^x dy\,
  \left[\delta(y-x)-x^{-1}\sc{M}\left({x}/{y}\right)\right]y^{-(d-2)/4} = }
\nonumber\\
&=& x^{(2-d)/4}-\frac{d-2}{2x}\int_0^x
 dy\,(x/y)^{(2-d)/4}(1-y/x)^{(d-4)/2}y^{(2-d)/4}\
\\
&=& x^{(2-d)/4}-\frac{d-2}{2}x^{(2-d)/4}\int_0^1 dz\,(1-z)^{(d-4)/2} =  0
\label{identity}
\end{eqnarray}
Multiplying this by $({4\Tct_{\rm w}}/d) \int_0^1 dy_{\rm w}\, [\delta(y_{\rm w}-1)-
\sc{M}(1/y_{\rm w})]y_{\rm w}^{(d+2)/4}$ and subtracting from\eq{Ione} exactly
cancels all contributions in the range $y>y_{\rm w}$, giving
\begin{eqnarray}
\fl C_{\rm m}^{(1)}(t,t_{\rm w})&=&\frac{4T_{\rm c} t_{\rm w}}{dx} \int_0^1 dy_{\rm w}\,
\int_0^{y_{\rm w}} dy\, \sc{M}\left({x}/{y}\right)
\left[\delta(y_{\rm w}-1)-\sc{M}\left({1}/{y_{\rm w}}\right)\right]
\nonumber\\
\fl & & \times \left[y_{\rm w}\left(\frac{y_{\rm w}}{y}\right)^{(d-2)/4}
- y\left(\frac{y}{y_{\rm w}}\right)^{(d-2)/4}\right]
\label{C_m1}
\end{eqnarray}
This shows that $C_{\rm m}^{(1)}/t_{\rm w}$ is a scaling function of
$x=t/t_{\rm w}$. Its full $x$-dependence has to be found numerically
from\eq{C_m1} or via series
expansion~\cite{Annibale_thesis},
but we can obtain the large-$x$ behaviour that is required for the
asymptotic FDR $X^\infty$ in closed form. For $x\rightarrow\infty$
one can replace the function
$\sc{M}(x/y)$ with its asymptotic form
$[(d-2)/2](x/y)^{(2-d)/4}$ from\eq{M_scaling_final} to get
\begin{eqnarray}
\fl \frac{C_{\rm m}^{(1)}(t,t_{\rm w})}{\Tct_{\rm w}}
&=&\frac{2(d-2)}{d} x^{-(d+2)/4}
\int_0^1 dy_{\rm w}\, \int_0^{y_{\rm w}} dy\,y^{(d-2)/4}
\left[\delta(y_{\rm w}-1)-\sc{M}\left(\frac{1}{y_{\rm w}}\right)\right]
\nonumber\\
\fl& & \times \left[y_{\rm w}\left(\frac{y_{\rm w}}{y}\right)^{(d-2)/4}
-y\left(\frac{y}{y_{\rm w}}\right)^{(d-2)/4}\right]
\\
\fl&=&\frac{2(d-2)}{d+2} x^{-(d+2)/4}
\left[1-\frac{d-2}{2} \frac{\Gamma\left(\frac{d+4}{2}\right)
\Gamma\left(\frac{d-2}{2}\right)}{\Gamma(d+1)}\right]
\label{C_m1inf}
\end{eqnarray}
This exhibits the expected leading order cancellation for large
$x$, which gives an additional factor of $1/x$ compared to the naive
result $x^{(2-d)/4}$.

To complete the calculation of the correlation function we need to
evaluate $C_{\rm m}^{(2)}$ from\eq{c_m2}, which cannot be neglected for
$d<4$.  This requires the long-time behaviour of $\tilde{C}\Cc(t,t_{\rm w})$, which is
given by\eq{sc_CC_ord} and\eq{sc_CC_dis} for $t>t_{\rm w}$ and
$t<t_{\rm w}$, respectively, as for the unmagnetized case. The only
modification arises from the different behaviour of $g(t)$. One thus finds
\begin{equation}
\fl \frac{\tilde{C}\Cc(t,t_{\rm w})}{\tilde{C}\Cc(t,t)}={\mathcal{G}}(t/t_{\rm w}),
\quad
{\mathcal{G}}(x)=\left\{
\begin{array}{ll}
{\displaystyle \frac{\int dw\, w^{(d-6)/2} \sc{C}^2(w) e^{-2(x-1)w}}
{x\int dw\, w^{(d-6)/2} \sc{C}^2(w)}} & \mbox{for $x\geq 1$} \\
x^{(d-6)/2}{\mathcal{G}}(1/x) & \mbox{for $x\leq 1$}
\end{array}
\right.
\label{CCt_scaling}
\end{equation}
The scaling of the equal-time value of $\tilde{C}\Cc$
is, from\eq{sc_CC_ord}, $\tilde{C}\Cc(t,t)=\tilde\gamma_d t^{(4-d)/2}$ with
$\tilde\gamma_d = \sigma_d T_{\rm c}^2\int dw\, w^{(d-6)/2} \sc{C}^2(w)$;
compare\eq{CCd}. Inserting this into\eq{c_m2} gives
\begin{eqnarray}
\fl\lefteqn {C_{\rm m}^{(2)}(t,t_{\rm w})=
\frac{1}{2} \int\, dt'\,dt_{\rm w}'\, M(t,t')M(t_{\rm w},t_{\rm w}') \tilde{C}\Cc(t',t_{\rm w}')}
\\
\fl &=&\frac{1}{2} \int\, dt'\,dt_{\rm w}'\, \frac{1}{m(t')m(t_{\rm w}')}
\frac{1}{t}\sc{M}(t/t')\frac{1}{t_{\rm w}}\sc{M}(t_{\rm w}/t_{\rm w}')
\tilde{C}\Cc(t',t') \mathcal{G}(t',t_{\rm w}')
\\
\fl &=&\frac{t_{\rm w} \tilde\gamma_d}{2\mu_d x}
\int_0^1 dy_{\rm w}\,y_{\rm w}^{(d-2)/4}\sc{M}(1/y_{\rm w})
\int_0^{x}dy\, \sc{M}(x/y) y^{(6-d)/4} {\mathcal{G}}(y/y_{\rm w})
\\
\fl &=&\frac{t_{\rm w}}{x} \int_0^1 dy_{\rm w}\, y_{\rm w}^2 \sc{M}(1/y_{\rm w})
\int_0^{x/y_{\rm w}}du\, \sc{M}(x/uy_{\rm w})u^{(6-d)/4}
\frac{\tilde\gamma_d}{2\mu_d}{\mathcal{G}}(u)
\label{general_Cm2}
\end{eqnarray}
In the second line we have used\eq{combination}
to write $M(t,t')=m(t)\delta(t-t')-t^{-1}m^{-1}(t')\sc{M}(t/t')$ up to
negligible corrections, and then
immediately discarded the $\delta$-function contributions, which give
subleading corrections.

Let us denote by $U$ the value of the $u$-integral
in\eq{general_Cm2}. Since $\mathcal{G}(u)$ is defined separately for
$u>1$ and $u<1$ in\eq{CCt_scaling}, one splits the integral accordingly:
\begin{eqnarray}
\fl \frac{2\mu_d U}{\tilde\gamma_d} &=&
\int_0^1\!\! du\, \sc{M}(x/uy_{\rm w})u^{(d-6)/4} {\mathcal{G}}(1/u)+
\int_1^{x/y_{\rm w}}\!\! du\, \sc{M}(x/uy_{\rm w})u^{(6-d)/4} {\mathcal{G}}(u)
\\
\fl &=& \int_1^\infty\!\! du\, \sc{M}(xu/y_{\rm w})u^{-(d+2)/4} {\mathcal{G}}(u)+
\int_1^{x/y_{\rm w}}\!\! du\, \sc{M}(x/uy_{\rm w})u^{(6-d)/4} {\mathcal{G}}(u)
\label{decomposition}
\end{eqnarray}
We now need ${\mathcal{G}}(u)$ for $u>1$. The denominator
in\eq{CCt_scaling} is $\tilde\gamma_d/(\sigma_d T_{\rm c}^2)$,
and bearing in mind the
definition\eq{C_scaling_magn} of $\sc{C}$ with $\alpha=(d-2)/2$ gives
\begin{eqnarray}
\fl \frac{\tilde\gamma_d}{\sigma_d T_{\rm c}^2} u{\mathcal{G}}(u) &=&
\int dw\, w^{(d-6)/2}\sc{C}^2(w)e^{-2(u-1)w}
\\
\fl &=& 4\int dw\, w^{(d-2)/2}\int_0^1 dy\,\int_0^1 dy'\,
(yy')^{(d-2)/2}e^{-2w(1-y-y'+u)}
\\
\fl &=&2^{(4-d)/2}\Gamma\left({\textstyle\frac{d}{2}}\right)
\int_0^1 dy\,\int_0^1 dy'\, (yy')^{(d-2)/2}(1-y-y'+u)^{-d/2}
\label{3_dim_int}
\end{eqnarray}
We now insert this into\eq{decomposition} and simplify the
numerical prefactors by using $\mu_d=\lambda_d^{-1}(d-2)(4-d)/4$ and
the explicit expression\eq{Ld} for $\lambda_d$ to get
\begin{eqnarray}
\fl
U&=& \frac{T_{\rm c}}{\Gamma(\frac{d-2}{2}) \Gamma(\frac{4-d}{2})}
\left[\int_1^{\infty}du\, \sc{M}(xu/y_{\rm w})u^{-(d+6)/4} \right. \nonumber\\
\fl & &\times \int_0^1 dy \int_0^1 dy'\,(yy')^{(d-2)/2}(1-y-y'+u)^{-d/2}
\nonumber\\
\fl &&{}+{}\left. \int_1^{x/y_{\rm w}}du\, \sc{M}(x/uy_{\rm w})u^{(2-d)/4}
\int_0^1 dy \int_0^1 dy'\, \ldots\right]
\label{exact_I}
\end{eqnarray}
This is the $u$-integral from\eq{general_Cm2} and so overall we have
a four-dimensional integral over $y_{\rm w}, u, y, y'$ for $C_{\rm m}^{(2)}$. In
general this cannot be evaluated in closed form; a series expansion is
given
in~\cite{Annibale_thesis}.
The large-$x$ behaviour, which will give
us the asymptotic FDR, is easier to extract. In the first $u$-integral
of\eq{exact_I} one can directly use the asymptotic form of
$\sc{M}(xu/y_{\rm w})$. One can show that for large $x$ the same replacement
can be made in the second integral, and the upper integration limit sent to
infinity thereafter.
This gives for the large-$x$ behaviour of $U$
\begin{equation}
U= \frac{d-2}{2}\frac{T_{\rm c} V_d}{\Gamma(\frac{d-2}{2}) \Gamma(\frac{4-d}{2})}
x^{(2-d)/4} y_{\rm w}^{(d-2)/4}
\label{U_large_x}
\end{equation}
where $V_d$ is a $d$-dependent numerical constant given by
\begin{equation}
\fl V_d=\int_1^{\infty}\!\!du\, (u^{-(d+2)/2}+1) \int_0^1\! dy \int_0^1\! dy'\,
(yy')^{(d-2)/2}(1-y-y'+u)^{-d/2}
\label{integral_sums}
\end{equation}
Inserting\eq{U_large_x} into\eq{general_Cm2}, the remaining
$y_{\rm w}$-integral can be done explicitly to give
\begin{eqnarray}
\fl C_{\rm m}^{(2)}(t,t_{\rm w})=\left(\frac{d-2}{2}\right)^2
\frac{\Gamma(\frac{d+4}{2}) V_d}{\Gamma(d+1)\Gamma(\frac{4-d}{2})}
\Tct_{\rm w} x^{-(d+2)/4}
\label{C_m2inf}
\end{eqnarray}
As anticipated this has the same scaling as the first
contribution\eq{C_m1inf} to the correlation function, so that overall
for large $x$
\begin{eqnarray}
\fl C_{\rm m}(t,t_{\rm w}) = C_{\rm m}^{(1)}+C_{\rm m}^{(2)}&=&\frac{d-2}{2}\Tct_{\rm w}
x^{-(d+2)/4}\left[\frac{4}{d+2}
\left(1-\frac{d-2}{2}\frac{\Gamma(\frac{d+4}{2})\Gamma(\frac{d-2}{2})}
{\Gamma(d+1)}\right)\right.
\nonumber\\
& &{}+{}\left.\frac{d-2}{2}\frac{\Gamma(\frac{d+4}{2})V_d}
{\Gamma(d+1)\Gamma(\frac{4-d}{2})}\right]
\label{Cm_tot_infty}
\end{eqnarray}

The magnetization {\em response} function is rather easier to
find, by using\eq{R_simplification_magn} and\eq{combination}
in\eq{response_magn}
and rescaling the integration variable to $y=t'/t_{\rm w}$ as usual:
\begin{equation}
R_{\rm m}(t,t_{\rm w}) = \int_1^x dy\,\left[\delta(y-x)-x^{-1}\sc{M}(x/y)\right]y^{-(d-2)/4}
\end{equation}
The structure of this is rather similar to $C_{\rm m}^{(1)}$, and
by subtracting the vanishing term\eq{identity} one again gets a
significant cancellation,
\begin{equation}
R_{\rm m}(t,t_{\rm w}) = \frac{1}{x}\int_0^1 dy\,\sc{M}(x/y)y^{(2-d)/4}
\end{equation}
Inserting the explicit form\eq{M_scaling_final} of $\sc{M}$ one then
gets simply
\begin{equation}
R_{\rm m}(t,t_{\rm w}) = x^{(2-d)/4}\left[1-\left(1-\frac{1}{x}\right)^{(d-2)/2}\right]
\label{R_finite}
\end{equation}
which for $x\rightarrow\infty$ behaves as
\begin{equation}
R_{\rm m}(t,t_{\rm w}) = \frac{d-2}{2}x^{-(2+d)/4}
\label{R_infty}
\end{equation}

With the expressions $\eq{Cm_tot_infty}$ and $\eq{R_infty}$ for
correlation and response in the limit of long, well-separated times we
can now finally compute the asymptotic FDR defined by
\begin{eqnarray}
\fl X_{\rm m}^\infty &=&
\lim_{t\ggt_{\rm w}\gg 1}\frac{T_{\rm c}R_{\rm m}(t,t_{\rm w})}{C_{\rm m}'(t,t_{\rm w})}
=\frac{4}{d+6}\left[\frac{4}{d+2}
\left(1-\frac{d-2}{2}\frac{\Gamma(\frac{d+4}{2})\Gamma(\frac{d-2}{2})}
{\Gamma(d+1)}\right)\right.
\nonumber\\
\fl & &{}+{}\left.\frac{d-2}{2}\frac{\Gamma(\frac{d+4}{2}) V_d}
{\Gamma(d+1)\Gamma(\frac{4-d}{2})}\right]^{-1}
\label{Xz}
\end{eqnarray}
where $V_d$ is given by\eq{integral_sums} and the prefactor $4/(d+6)$ accounts
for the fact that $C_{\rm m}\sim t_{\rm w} x^{-(2+d)/4} \simt_{\rm w}^{(d+6)/4}$ and hence
$C_{\rm m}'(t,t_{\rm w})=[(d+6)/4t_{\rm w}]C_{\rm m}(t,t_{\rm w})$.

Before exploring this result, let us comment briefly on the limit as
$d\to 4$, which should make contact with our results in the
previous subsection. The contribution from $C_{\rm m}^{(2)}$, which appears
in the second line of\eq{Cm_tot_infty} and\eq{Xz}, vanishes linearly
in $4-d$ in this limit because of the factor $\Gamma^{-1}((4-d)/2)$;
$V_d$ stays finite as we show below. This is consistent with the fact
that in dimension $d>4$ this term does not contribute.  In the limit
$d\to 4$ one has, from\eq{Cm_tot_infty}, $C_{\rm m}=\Tct_{\rm w} x^{-3/2}/2$ for
large $x$ which matches precisely\eq{Cm_dgt4} for
$d>4$. Similarly, the large-$x$
magnetization response\eq{R_finite} for $d\to 4$ is $R_{\rm m}=x^{-3/2}$ in
agreement with\eq{Rm_dgt4}.

We now look in more detail at the $d$-dependence of the asymptotic FDR\eq{Xz}
for the magnetization. Expanding in $\epsilon=4-d$ one has
\begin{eqnarray}
\fl X_{\rm m}^\infty
&=&\left(\frac{2}{5}+\frac{\epsilon}{25}+{{\mathcal O}}(\epsilon^2)\right)
\left[\left(\frac{2}{3}+\frac{\epsilon}{9}\right)
\left(\frac{3}{4}-\frac{\epsilon}{6}\right)+
\frac{\epsilon}{8} V_4+{{\mathcal O}}(\epsilon^2)\right]^{-1}
\nonumber\\
&=&\frac{4}{5}+\left(\frac{28}{225}-\frac{V_4}{5}\right)\epsilon
+{{\mathcal O}}(\epsilon^2)
\end{eqnarray}
where $V_4$ is the limiting value of $V_d$ for $d\to 4$, which can be
worked out explicitly as
\begin{equation}
V_4=-\frac{2}{3}\ln 2+\frac{11}{12}+\frac{\pi^2}{24}
\end{equation}
so that
\begin{equation}
\fl X_{\rm m}^\infty=\frac{4}{5}+\left(\frac{2\ln 2}{15}
-\frac{53}{900}-\frac{\pi^2}{120}\right)\epsilon+{{\mathcal O}}(\epsilon^2)
\label{X_m_epsilon_expansion}
\end{equation}
It is remarkable that in a system as simple as the spherical model,
where standard critical exponents are rational functions of the
dimension $d$, the magnetization FDR for magnetized initial states is
very much more complicated, and irrational already {\em to first order} in
$\epsilon$.

\begin{figure}
\begin{center}
\centerline{\includegraphics[width=10.0cm,clip=true]{Xinf_new.eps}}
\caption{Asymptotic FDR $X^\infty_{\rm m}$ for the magnetization, for
critical coarsening with nonzero initial magnetization. Solid line:
Full theory~(\protect\ref{Xz}) including non-Gaussian corrections;
dotted lines indicate the first-order expansions near $d=2$ and
4. Dashed line: Gaussian theory~(\protect\ref{Xqneq}).
\label{fig:X_magnetized}
}
\end{center}
\end{figure}
In the opposite limit $d\to 2$, the $u$-integral in the
definition\eq{integral_sums} of $V_d$ diverges at the upper end;
dropping all non-divergent corrections gives the leading divergence as
\begin{equation}
\fl V_d\approx\int_1^{\infty}du\,\int_0^1 dy \int_0^1 dy'\,(1-y-y'+u)^{-d/2}
\approx \int_1^{\infty}du\,u^{-d/2} = \frac{2}{d-2}
\label{limit_d2}
\end{equation}
This divergence balances the vanishing prefactor
$(d-2)/2$ in\eq{Cm_tot_infty} and\eq{Xz} whereas the contribution in
the round brackets coming from $C_{\rm m}^{(1)}$ vanishes linearly with
$d-2$. Consequently the asymptotic behaviour of the correlation
function is, for $d$ close to $2$, $C_{\rm m}=[(d-2)/2]\Tct_{\rm w} x^{-1}$ and
the FDR becomes $X_{\rm m}^\infty=4/(d+6)=1/2$. One can also obtain
the leading order correction in $\epsilon'=d-2$, which is given by
\begin{equation}
X_{\rm m}^\infty=
\frac{1}{2}+\left(\frac{1}{16}+\frac{\pi^2}{48}\right)\epsilon'
+{{\mathcal O}}(\epsilon'^2)
\label{X_m_epsilon_prime_expansion}
\end{equation}
The only subtlety here is working out the subleading term of
$V_d$. Setting $V_d=2/\epsilon'+a_0+\ldots$, $a_0$ can be obtained as
the limit for $d\to 2$ of\eq{integral_sums} with $u^{-d/2}$ subtracted
from the integrand. The limit can be taken in the integrand itself for
finite $u$, giving a convergent integral with value
$3/2-\pi^2/12$. But one has to account separately for the large-$u$
tail $u^{-d/2}[-1+\int dy\,dy' (yy')^{(d-2)/2}]$ which integrates to
$(-1+4/d^2)[2/(d-2)]\to -2$ for $d\to 2$, giving $a_0=-1/2-\pi^2/12$
overall.

In summary, the asymptotic magnetization FDR $X_{\rm m}^\infty$ decays
from $4/5$ in $d=4$ to $1/2$ in $d=2$. Fig.~\ref{fig:X_magnetized} also
shows numerical values for intermediate $d$. $X_{\rm m}^\infty$ is
larger than the
FDR\eq{Xqneq} from the Gaussian theory except in the limit $d\to2$; it also
remains larger throughout than the FDR\eq{X_baseline} for
unmagnetized initial states, shown in Fig.~\ref{fig:X_normal_and_E}.

\begin{figure}
\begin{center}
\centerline{\includegraphics[width=9.0cm,clip=true]{d234_new.eps}}
\caption{Normalized magnetization FD plot for $d=2, 3, 4$, showing
normalized susceptibility $\tilde{\chi}_{\rm m}$ versus normalized
correlation $\tildeC_{\rm m}$ in the limit of long times. For $d=4$ the plot
is a straight line with (negative) slope $4/5$ as expected. Increasing
deviations from this appear as $d$ decreases towards $2$.
\label{fig:mFD_plot}
}
\end{center}
\end{figure}
We next turn to the shape of the FD plot for the magnetization. As
explained in the introduction, this
is obtained by plotting the normalized susceptibility
$\tilde{\chi}_{\rm m}(t,t_{\rm w}) = T_{\rm c}\chi_{\rm m}(t,t_{\rm w})/C_{\rm m}(t,t)$ versus
$\tildeC_{\rm m}(t,t_{\rm w})=C_{\rm m}(t,t_{\rm w})/C_{\rm m}(t,t)$.
In the limit $d\to 4$ the FD
plot must be a straight line with slope $X_{\rm m} = 4/5$, by
continuity with the results for $d>4$. Numerical evaluation
(see Fig.~\ref{fig:mFD_plot}%
) shows that, as the dimensionality decreases,
the FD plots deviate progressively from this straight lines. The most
extreme case is the limit $d\to 2$, where analytical forms can be
found.

To find the correlation function for $d\to 2$, it is useful to note that
the scaling function $\sc{M}(x)$ becomes equal to $\delta(x-1)$ in the
limit. Formally, one sees easily from\eq{M_scaling_final} that for
any smooth bounded function $f(x)$ and fixed $c>1$,
$\int_1^c dx\,
\sc{M}(x) f(x) \to f(1)$ because the divergence of $\sc{M}(x)$ at
$x=1$ becomes non-integrable in $d=2$. We now exploit this to
simplify $U$ from\eq{exact_I}.
In the first term in square brackets, the argument $xu/y_{\rm w}$ of
$\sc{M}$ is always larger than $x/y_{\rm w}$ and hence than $1$ (except at
the irrelevant boundary $y_{\rm w}=x$ of the $y_{\rm w}$-integral); in the limit
$d\to 2$, this contribution therefore vanishes. In the second term,
replacing $\sc{M}(x/uy_{\rm w})$ by $\delta(x/uy_{\rm w} - 1)$ and expanding the
prefactor to leading order in $(d-2)/2$ gives
\begin{equation}
\frac{U}{T_{\rm c}} = \frac{d-2}{2} \frac{x}{y_{\rm w}}
\int_0^1 dy \int_0^1 dy'\,(1-y-y'+x/y_{\rm w})^{-1}
\label{d2_expansion}
\end{equation}
Performing the integrals over $y$ and $y'$, one has
\begin{equation}
\frac{U}{T_{\rm c}}=\frac{d-2}{2} \frac{x}{y_{\rm w}}
\left[\frac{x}{y_{\rm w}}\ln\left(1-\frac{y_{\rm w}^2}{x^2}\right)+
\ln\left(\frac{x+y_{\rm w}}{x-y_{\rm w}}\right)\right]
\end{equation}
We can now use this limit form of $U$ to get the contribution
$C_{\rm m}^{(2)}$ to the correlation function; recall that $U$ was defined
as the $u$-integral in\eq{general_Cm2}. In the remaining
$y_{\rm w}$-integral in this equation one can again replace $\sc{M}(1/y_{\rm w})$ by
$\delta(1/y_{\rm w} -1)$ so that only $y_{\rm w}=1$ contributes, giving
\begin{equation}
C_{\rm m} = C_{\rm m}^{(2)}=\frac{d-2}{2}T_{\rm c} t_{\rm w} \left[x\ln\left(1-\frac{1}{x^2}\right)+
\ln\left(\frac{x+1}{x-1}\right)\right]
\label{Cm_d2}
\end{equation}
One can show that the other contribution $C_{\rm m}^{(1)}$ to the
magnetization correlator, given by\eq{C_m1}, vanishes as $\sim(d-2)^2$ for
$d\to 2$ so that to leading order $C_{\rm m}=C_{\rm m}^{(2)}$ as anticipated
in writing\eq{Cm_d2}. The corresponding response function
is found by expanding\eq{R_finite} to leading order in $d-2$:
\begin{equation}
R_{\rm m} = \frac{d-2}{2}\ln\left(\frac{x}{x-1}\right)
\end{equation}
To get the FDR it only remains to work out the $t_{\rm w}$-derivative of
$C_{\rm m}$:
\begin{eqnarray}
\fl C_{\rm m}' \ \ \equiv \ \ \partial_{t_{\rm w}}C_{\rm m}&=&
\frac{d-2}{2}\,T_{\rm c}\ln\left(\frac{x+1}{x-1}\right)
\end{eqnarray}
We therefore get, for the full dependence of the limiting FDR for
$d\to 2$ on scaled time $x=t/t_{\rm w}$
\begin{equation}
X_{\rm m}(x)=\frac{T_{\rm c}R_{\rm m}}{C_{\rm m}'} = \ln\left(\frac{x}{x-1}\right)
\left[\ln\!\left(\frac{x+1}{x-1}\right)\right]^{-1}
\end{equation}
For $x\to\infty$ this gives $X_{\rm m}^\infty=1/2$ consistent with the
discussion above. In the limit $x\to 1$ of comparable times, on the
other hand, $X_{\rm m}(x)$ approaches 1, logarithmically slowly; the FD
plot for $d\to 2$ therefore starts off with a pseudo-equilibrium
slope. Interestingly, this implies that the trends of the slope with
$d$ are different at the two ends of the plot: for well-separated
times ($x\to\infty$), the slope {\em decreases} from $4/5$ to $1/2$ as
$d$ decreases from 4 to 2; for comparable times ($x\to 1$) it {\em
increases} from 4/5 to 1.

To get the FD plot itself we need the susceptibility, which is found
by integration of the response\eq{R_finite} as
\begin{eqnarray}
\chi_{\rm m}(t,t_{\rm w})&=&\int_{t_{\rm w}}^t dt'\, R_{\rm m}(t/t')=t\int_{1/x}^1dz\, R_{\rm m}(1/z)
\\
&=& t\int_{1/x}^1 dz\,z^{(d-2)/4}
\left[1-\left(1-z\right)^{(d-2)/2}\right]
\label{chi_fin}
\end{eqnarray}
Expanding to linear order in $d-2$ and integrating gives
\begin{equation}
\chi_{\rm m}(t,t_{\rm w})
=\frac{d-2}{2}\,t\,\frac{x-1}{x}\left[1-\ln\left(\frac{x-1}{x}\right)\right]
\end{equation}
The normalized susceptibility $\tilde\chi_{\rm m}$ is obtained by
dividing by $C_{\rm m}(t,t)/T_{\rm c}$, which from\eq{Cm_d2} for $x\to 1$ equals
$2\ln 2[(d-2)/2]\,t$.
The resulting FD plot for $d\to 2$ is shown in
Fig.~\ref{fig:mFD_plot}, together with the ones for $d=3$
(determined numerically%
) and $d=4$.  For $d\to
2$ the approach of the slope to the equilibrium value $X=1$ for $x\to
1$ is difficult to see because it is logarithmically slow. As
expected from the trends with $d$ in the initial and final slopes of
the FD plot, the curve for $d\to 2$ is the most strongly curved, while
in $d=4$ we have the anticipated straight line with slope $4/5$
that is required by continuity with the results for $d>4$.

To quantify the shape of the FD plot further one can also consider the
axis ratio $Y$, defined as the limiting value of $\tilde\chi_{\rm
m}(t,t_{\rm w})=T_{\rm c}\chi_{\rm m}(t,t_{\rm w})/C_{\rm m}(t,t)$ for well-separated times
$t\ggt_{\rm w}\gg 1$. In equilibrium this would correspond to the FDT for
the static quantities, i.e. equal-time fluctuations and static
susceptibilities. Out of equilibrium, if the FD plot is straight then
$Y$ coincides with $X$; if it is not, then is has been argued that in
some circumstances $Y$ can be more relevant for characterizing
effective temperatures than $X$~\cite{OHeLiuNag04}. For the
magnetization FD plot in the current scenario of magnetized initial
states, we see from Fig.~\ref{fig:mFD_plot} that $Y$ decreases along
with $X_{\rm m}^\infty$ from $4/5$ as the dimension is lowered below
$d=4$. The two quantities begin to
differ more noticeably as $d$ decreases further, with $X_{\rm
m}^\infty=1/2$ and $Y=1/(2\ln 2)=0.7213\ldots$ in the limit $d\to 2$.

\section{Summary and discussion}
\label{sec:conclusion}

In this paper we have considered the non-equilibrium dynamics of the
spherical ferromagnet after a quench to the critical temperature
$T_{\rm c}$. Our focus has been the calculation of correlation and response
functions and the associated fluctuation-dissipation ratios (FDRs)
$X(t,t_{\rm w})$. The key quantity that can be extracted from the latter is
the asymptotic FDR $X^\infty$ for large and well-separated times
$t\ggt_{\rm w}\gg 1$; it is independent of model-specific details within a
given dynamical universality
class. We were motivated by two questions: how does $X^\infty$
depend on the observable considered, both with regard to the
lengthscale and the type of observable (spin, bond, spin product)? And
what is the effect of initial conditions, in particular the presence
of a nonzero magnetization in the initial state? The first question
has implications for the interpretation of $T/X^\infty$ as an
effective temperature, which is plausible only if this quantity is
observable-independent. The second one allowed us to uncover whether
different initial conditions can lead to different universality
classes of critical coarsening.

A peculiarity of the spherical model is the weak infinite-range interaction
produced by the spherical constraint. This requires that one
distinguishes between long-range or ``block'' observables, which probe
lengthscales large compared to the (time-dependent) correlation length
but small compared to the system size, and global observables whose
behaviour depends on correlations across the entire
system. Technically, the first case is much easier to treat because
the standard theory where the spins have Gaussian statistics can be
used. Global correlation and response functions, on the other hand,
require non-Gaussian corrections arising from the fluctuations of the
effective Lagrange multiplier.

We dealt with the case of finite-range (i.e.\ either local or
long-range) observables in Sec.~\ref{sec:finite-range}. For spin
observables, we found in the long-range case\eq{X_baseline} the same
$X^\infty$ as for local spin correlations and
response~\cite{GodLuc00b}. This was as expected from the general
correspondence between local and long-range observables discussed in
the introduction. The FD plot for the long-range spin observable,
i.e.\ the magnetization, is a straight line in the long-time
limit. This is as in the Ising case in $d=1$~\cite{MayBerGarSol03},
but it is interesting to note that here it holds for all dimensions
$d>2$. On the other hand, in the Ising model with $d\geq 2$, RG
arguments have been adduced~\cite{CalGam05,MayBerGarSol04,CalGam02} to
suggest that the magnetization FD plot should not be straight, though
with deviations that are likely too small to be detectable
numerically~\cite{MayBerGarSol04}. It is likely that the Gaussian
statistics of the spherical model are responsible for producing a
simpler, straight-line magnetization FD plot in all dimensions,
although it would be interesting to know whether any other models have
this property.

We then looked at the effect of the type of observable on $X^\infty$,
considering both bond and spin product observables in either the local
or long-range versions. The results in equations~(\ref{X_bond_local}),
(\ref{X_bond_long-range}), (\ref{X_prod_local}), (\ref{Xbl_large_d})
show that, although the precise time-dependence of $X(t,t_{\rm w})$ varies,
the asymptotic FDR $X^\infty$ is the same in all cases. This is
consistent with general arguments~\cite{CalGam04} suggesting that for
a Gaussian theory all observables should yield the same $X^\infty$. In
contrast to the Ising case~\cite{MayBerGarSol03}, not all long-range
observables give nontrivial FD plots; in fact, only the block product
observable does so, and only for $d<4$, while all others produce
pseudo-equilibrium FD plots for long times.

The bulk of the paper was concerned with the more challenging analysis
of global correlation and response functions, focussing mostly on the
energy as a key observable. In Sec.~\ref{sec:setup} we constructed a
framework for calculating non-Gaussian corrections to the
spins, which are ${{\mathcal O}}(N^{-1/2})$ to leading order. This lead to the
general expression\eq{yt3} for these leading-order corrections. It
involves a two-time kernel $L(t,t_{\rm w})$ which from\eq{pinv_def} is the
functional inverse of $K(t,t_{\rm w})$ defined in\eq{p_def}. The basis of
all subsequent calculations is the determination of the long-time scaling of
these two functions, as summarized at the end of Sec.~\ref{sec:KandL}.

In Sec.~\ref{sec:energy_general} we obtained general expressions for
energy correlation and response functions, in terms of the kernel $L$
and other quantities known from the Gaussian theory; the results can
be found in eqs.~(\ref{C1},\ref{C2},\ref{C3},\ref{C4}) and\eq{RE}.
Evaluating these first for the equilibrium case, we found that the
energy correlation and susceptibility display a plateau for $T$ just
below $T_{\rm c}$ and $d>4$; this is caused by the $\mathbf{q}=\mathbf{0}$ wavevector,
i.e.\ by the slow relaxation of the global magnetization.  In
Sec.~\ref{sec:noneq_large_d}, we proceeded to the long-time analysis
of energy FD behaviour in the non-equilibrium case for $d>4$; the key
results are\eq{CE_d_gt_4_longtime} and\eq{RE_d_gt_4_longtime}.  The
associated FDR is given explicitly in\eq{XE} and has the {\em same}
asymptotic value $X^\infty=1/2$ as for all other (finite-range)
observables in $d>4$. The analysis of the case $d<4$ is more
difficult, and we were able to find closed-form
results\eqq{CE_prime}{RE_d_lt_4} only in the limit of well-separated
times $t/t_{\rm w}\gg 1$. This is, however, enough to determine $X^\infty$,
with the result\eq{Xinf_dlt4_explicit}. Evaluating this, both
numerically and by expansion in $4-d$ and $d-2$, the crucial
conclusion is that it does not coincide with the asymptotic FDR for
finite-range observables; see Fig.~\ref{fig:X_normal_and_E}. A naive
interpretation of $T/X^\infty$ as an
effective temperature for critical coarsening dynamics is therefore
ruled out, since such a temperaure ought to be observable-independent.
On the other hand, to first order in $4-d$ the result agrees with an
RG calculation~\cite{CalGam04} for the $O(n)$-model. We conclude that
non-Gaussian corrections to the FD behaviour of global observables in
the spherical model capture genuine physical effects that have close
counterparts in more realistic systems with only short-range
interactions.

Finally, in Sec.~\ref{sec:magnetized} we turned our attention to
critical coarsening starting from magnetized initial states;
physically this situation could be produced by an up-quench from an
equilibrated state at a starting temperature $T<T_{\rm c}$. We concentrated
on the simpler spin observables and found that already for them, the
presence of a nonzero magnetization makes global properties sensitive
to non-Gaussian corrections. As with the energy fluctuations, it is the
{\em global} correlation and response functions that make contact with
the results for short-range models, as obtained recently for the Ising
case~\cite{GarSolPagRit05}: we find $X^\infty_{\rm m}=4/5$ for $d>4$,
eq.\eq{X_m_d_gt_4}. This is distinct from the value $X^\infty=1/2$ for
the unmagnetized case, indicating that magnetized critical coarsening
is in a separate dynamical universality class. Surprisingly, even the
expressions
for correlation and response functions themselves, which are not
expected to be universal, coincide with those for the Ising case.
It remains to be understood whether this is accidental or has more
profound origins. For the case $d<4$, we obtained new {\em exact} values for
the asymptotic FDR of magnetized critical coarsening. The
magnetization response\eq{R_finite} can be found
explicitly for long times, while for the magnetization
correlator only the asymptotics for well-separated
times\eq{Cm_tot_infty} can be written in closed form. The resulting
$X^\infty_{\rm m}$, eq.\eq{Xz}, is surprisingly nontrivial: while it
matches continuously with $X^\infty_{\rm m}=4/5$ in $d>4$ and
approaches the simple value $X^\infty_{\rm m}=1/2$ for $d\to 2$ as
shown in Fig.~\ref{fig:X_magnetized}, it is irrational already to
first order in an expansion in $4-d$, eq.\eq{X_m_epsilon_expansion}, or
$d-2$, eq.\eq{X_m_epsilon_prime_expansion}.

While the conclusion of our calculation as regards the existence of a
well-defined effective temperature for critical coarsening is
negative, the issue of dynamic universality classes and new asymptotic
FDRs due to magnetized (and possibly other, different) initial
conditions clearly deserves further study. Results for systems with
short-range interactions, such as the $O(n)$ and $n$-vector models,
would be particularly welcome. After the present work was completed we
became aware that a first step in this direction has recently
been taken by the authors of Ref.~\cite{FedTri06}, who calculated the
FDR for the $n$-vector model with a magnetized initial state within an
$\epsilon$-expansion around $d=2$. Intriguingly, their result
$X^\infty=1/2$ for $d=2$ itself agrees with ours, but the first-order
correction in $d-2$ remains {\em rational} even for $n\to\infty$. It
therefore disagrees with our spherical model
result\eq{X_m_epsilon_prime_expansion}. This appears to be the first
example of genuine differences between the spherical and $n$-vector
(with $n\to\infty$) models, which are known to have identical
properties in equilibrium~\cite{Stanley68} and within a Gaussian
theory of the dynamics.

As regards future work, we note first that a complete classification
of dynamical universality classes within critical coarsening remains
to be achieved. An earlier study of the spherical model considered
initial conditions with long-range correlations but no overall
magnetization; this yields no new (non-zero) values of the asymptotic
FDR $X^\infty$~\cite{PicHen02}. The presence of a non-zero
magnetization thus appears to be important for observing new
phenomena, and is reflected in our calculation by the fact that
non-Gaussian fluctuations become important. Whether there are yet
other initial conditions that could give rise to distinct values of
$X^\infty$ is an open problem.

Our general framework for treating non-Gaussian
corrections to the dynamics can also be applied in other contexts.
For example, it can be used to analyse the {\em
fluctuations} across thermal histories of correlation and response
functions that have been coarse-grained across a finite-sized
system. The properties of these fluctuations should be useful for
understanding dynamical heterogeneities in coarsening
dynamics~\cite{CasChaCugIguKen03}, and we will report on the
results of such a study shortly. We have also extended our approach to
non-Gaussian corrections for the dynamics of e.g.\ the $O(n)$-model
with large but finite $n$, opening up the attractive prospect of
obtaining exact results analogous to the ones in this paper for models
with exclusively short-range interactions.

{\bf Acknowledgements}: We thank P Calabrese and A Gambassi for
sharing their field theoretic results with us before publication, and
for helpful discussions. PS acknowledges the hospitality of the KITP,
Santa Barbara, where this research was supported in part by the NSF
under grant no.\ PHY99-07949.

redpajama_set_name: RedPajamaArXiv