\section{Introduction}
\label{sec:intro}
\begin{figure*}
    \centering
    \begin{minipage}[t]{0.49\linewidth}
        \centering
        \includegraphics[height=5cm]{figs/dspn_fig_v21.pdf}
        \vspace{-4.5em}
        \caption{
            \glsreset{DSPN}\Gls{DSPN} starts with a deterministic set that is gradually changed into the desired prediction with gradient descent on a learned loss function.
        }
        \label{fig:dspn}
    \end{minipage}
    \hfill
    \begin{minipage}[t]{0.49\linewidth}
        \centering
        \includegraphics[height=5cm]{figs/transformer_dspn_fig.pdf}
        \vspace{-5.5em}
        \caption{
            Our \glsreset{TSPN}\gls{TSPN} extends \gls{DSPN} to use more expressive permutation-equivariant set transformations based on the Transformer \cite{Vaswani2017transformer}.
            Additionally, \gls{TSPN} explicitly learns the set cardinality, which allows it to generalize to much larger sets.
        }
        \label{fig:transformer}
    \end{minipage}
\end{figure*}

It is natural to reason about a group of objects as a set.
Therefore many machine learning tasks involving predicting objects or their properties can be cast as a set prediction problem.
These predictions are usually conditioned on some input feature that can take the form of a vector, a matrix or a set.
Some examples include predicting future states for a group of molecules in a simulation \cite{No2020molecular}, object detection from images \cite{Carion2020detr} and generating correlated samples for sequential Monte Carlo in object tracking \cite{Zhu2020diffresampling,Neiswanger2014dependent}.
Elements of a set are unordered, which brings about two challenges that set prediction faces.
First, the model must be permutation-equivariant; that is, the generation of a particular permutation of the set elements must be equally probable to any other permutation.
Second, training a generative model for sets typically involves comparing a predicted set against a ground-truth set.
Since the result of this comparison should not depend on the permutation of the elements of either set, the loss function used for training must be permutation-invariant.
While it is possible to create a set model that violates either or both of these requirements, such a model has to \emph{learn} to meet them, which is likely to result in lower performance.

Permutation equivariance imposes a constraint on the structure of the model \cite{Bloemreddy2018exchangable,Bloemreddy2019probabilistic}, and therefore sets are often treated as an ordered collection of items, which allows using standard machine learning models.
For example assuming that a set has a fixed size, we can treat it as a tensor and turn set-prediction into multivariate regression \cite{Achlioptas2018aae}.
If the ordering is fixed but the size is not, we can treat set prediction as a sequence prediction problem \cite{Vinyals2016seq2set}.
Both approaches require using permutation-invariant loss functions to allow the model to learn a deterministic ordering policy \cite{Eslami2016air}.
However imposing such an ordering can lead to a pathology that is commonly referred to as the \textit{responsibility problem} \citep{Zhang2019dspn, Zhang2020fspool}; there exist points in the output set space where a small change in set space (as measured by a set loss) requires a large change in the generative model's output. This can lead to sub-optimal performance, as shown in \citet{Zhang2020fspool}.
Some approaches choose to learn the ordering of set elements \cite{Rezatofighi2018dpspn}, but this also suffers the same problem as well as adding further complexity to the set prediction problem.

Recently, \citet{Zhang2019dspn} introduced the \glsreset{DSPN}\gls{DSPN}---a model that generates sets in a permutation-equivariant manner using permutation-invariant loss functions.
\Gls{DSPN} relies on the observation that the gradient of a permutation-invariant function is equivariant with respect to the permutation of its inputs, also noticed by \citet{Papamakarios2019flowreview}.
\Gls{DSPN} uses this insight to generate a set by gradient-descent on a learned loss function with respect to an initially-guessed set.
\Gls{DSPN} has several limitations, however.
The functional form of the update step is limited, as the gradient information is only used to tanslate the set elements.
This, in turn, means that the method can be computationally costly: not only is the backward pass expensive, but many such passes might be needed to arrive at an accurate prediction.

In this paper we develop the \glsreset{TSPN}\gls{TSPN}, where we replace the gradient-based updates of \gls{DSPN} with a Transformer \cite{Vaswani2017transformer}, that is also permutation-equivariant and learns to jointly update the elements of the initial set.
We make the following contributions:
\begin{itemize}[leftmargin=*]
    \setlength\itemsep{-.25em}
    \item We show that the set-cardinality learning method of \citet{Zhang2019dspn} is prone to falling into local minima.
    \item We thus introduce an alternative, principled method for learning set cardinality.
    \item We learn a distribution over the elements of the initial set (as opposed to \gls{DSPN} that learns a fixed initial set).
    This allows one to directly generate sets with cardinality determined by the model (by sampling the correct number of points), and to dynamically change the size of the generated sets as test time.
    \item We demonstrate that \gls{TSPN} outperforms \gls{DSPN} on conditional point-cloud generation and object-detection tasks.
\end{itemize}
We show that our model is not only more expressive than the \gls{DSPN}, but can also generalize at test-time to sets of vastly different cardinality than the sets encountered during training.
We evaluate our model on auto-encoding \gls{SETMNIST} \cite{Lecun2010mnist} and on object detection on \gls{CLEVR} \cite{Johnson2017clevr}.
We now proceed to describe \gls{DSPN} and our method \gls{TSPN} in detail, followed by experiments.

\section{Permutation-Equivariant Set Generation}
\label{sec:background}

\subsection{Permutation-Equivariant Generator}

The \glsreset{DSPN}\gls{DSPN} iteratively transforms an initial set into the final prediction, and the transformation is conditioned on some input.
That is, given a conditioning $\by \in \RR^{d_y}$ and an initial set $\bx \coloneqq \set{\bx_i^0}_{i=1}^N$ of $N$ points in $\RR^d$, \gls{DSPN} iteratively applies a permutation-equivariant $f$ to transform this initial set into the final prediction $\bx^K \coloneqq \set{\bx_i^K}_{i=1}^N$ over $K$ iterations:
\begin{equation}
    \set{\bx_i^K}_{i=1}^N = f^K\left(\set{\bx_i^0}_{i=1}^N, \by\right)\,.%
\end{equation}
In \gls{DSPN}, the points $\bx_i^0$ are initialised randomly as model parameters and learned, hence the model assumes fixed set cardinality.
For handling variable set sizes, each element of the predicted set is augmented with a \textit{presence} variable $p_i \in \interval{0, 1}$, which are transformed by $f$ along with the $x_i$ and then thresholded to give the final prediction.
The ground-truth sets are padded with all-zero vectors to the same size.
Note, however, that this mechanism does not allow to extrapolate to set sizes beyond the maximum size encountered in training.

\Gls{DSPN} employs a permutation-invariant set encoder (e.g.\ deep-sets \cite{Zaheer2017deepset}, relation networks \cite{Santoro2017relnet}) to produce a set embedding, and updates the initial prediction using the gradients of an embedding loss, arriving at a permutation-equivariant final prediction.

This leads to the following set-generation procedure.
Given the input embedding $\bm{h} = \operatorname{\mathrm{input\_encoder}}(\by)$, an initial set\footnote{We drop the explicit indication of cardinality to avoid clutter.} $\bx^0$, and a permutation-invariant $\mathrm{set\_encoder}$, we can arrive at the final prediction by performing gradient descent,
\begin{align}
    \widehat{\bm{h}}^k &= \operatorname{\mathrm{set\_encoder}} \left( \bx^k \right)\,,\\
    \bx^{k} = \bx^{k-1} - \lambda &\nabla_\set{\bx} \operatorname{d}(\bm{h}, \widehat{\bm{h}}^{k-1}) \text{ for } k=1,\ldots,K\,, %
\end{align}
with step size $\lambda~\in~\RR_+$, distance function $d$ and number of iterations $K$.
The final prediction is $\bx^K$.
This set-generator can be used as a decoder of an autoencoding framework (with a permutation-invariant $\mathrm{input\_encoder}$) as well as for any conditional set-prediction task (\eg in object-detection where $\by$ is an image and $\mathrm{input\_encoder}$ is a \textsc{CNN}).

\subsection{Permutation-Invariant Loss} \label{sec:loss}

The model is trained using a set loss, i.e.\ a permutation-invariant loss. Common choices are the Chamfer loss or the Hungarian loss:
\begin{align}
    \loss[\mathrm{cham}](A, B) &= \sum_{a \in A} \min_{b \in B} \operatorname{d}(a, b) + \sum_{b \in B} \min_{a \in A} \operatorname{d}(a, b) \label{eq:chamfer}\\
    \loss[\mathrm{hung}](A, B) &= \min_{\pi \in P} \sum_{a_i \in A} \operatorname{d}(a_i, b_{\operatorname{\pi}(i)}) \label{eq:hungarian}\,,
\end{align}
where $\pi$ is a permutation in the space of all possible permutations $P$ and $d$ can be any distance or loss function defined on pairs of set points.
Note that the computational complexity of the Chamfer loss is $O(N^2)$ for sets of size $N$, whereas the Hungarian loss is $O(N^3)$---it uses the Hungarian algorithm to compute the optimal permutation, whose complexity is $\mathcal{O}(N^3)$ \cite{Bayati2008mwn}.
Hence the Chamfer loss is suitable for larger sets, and the Hungarian loss for smaller sets.
For $d$, \citet{Zhang2019dspn} use the Huber loss defined as $\operatorname{d}(\bm{a}, \bm{b}) = \sum_i \min\left(\nicefrac{1}{2}(a_i-b_i)^2, |a_i-b_i| - \nicefrac{1}{2}\right)$\,.

Recall that in the implementation of \gls{DSPN}, the ground truth set is padded with zero vectors so that all sets have the same size. Padding a set $A$ to a fixed size with constant elements $\widehat{a}$ turns it into a \textit{multiset} $\widehat{A}$. A Multiset is a set that contains repeated elements, that can be represented by an augmented set where each unique element is paired with its multiplicity.
If we use a multiset in its default form (\ie with repeated elements) as the ground-truth in the Chamfer loss, then it is enough for the model to predict a set $B$ containing exactly one element $b_i=\widehat{a}$ equal to the repeated element of the set $\widehat{A}$ in order to account for all its repetitions in the first term of \Cref{eq:chamfer}.
The remaining superfluous elements $b_j \in B$ predicted by the model can match \textit{any} other element in $A$ without increasing the second term of \Cref{eq:chamfer}.
This implies that padding a ground-truth set of size $N$ with $M$ constant elements creates $\binom{N+1}{M-1}$ predictions that are all optimal and hence indistinguishable under the Chamfer loss.
These predictions have a set size that varies from $N$ to $N+M-1$, and hence the model is likely to fail to learn the correct set cardinality---an effect clearly visible in our experiments, \cf \Cref{sec:set_mnist} and \Cref{tab:set_mnist}.

\Gls{DSPN} uses an additional regularization term {\small $\loss[\mathrm{repr}]{\bm{h}, \widehat{\bm{h}}} = \operatorname{d}(\bm{h}, \widehat{\bm{h}})$} defined as the distance between the model's input features and the encoding of the generated set, which improves the performance of the internal optimization loop.

\section{Size-Conditioned Set Generation with Transformers}
\label{sec:method}
We follow \gls{DSPN} in that we generate sets by transforming an initial set.
However, we employ Transformers \cite{Vaswani2017transformer,Lee2019settransformer} as a learnable set-transformation, which can readily account for any interactions between the set elements.
This change implies that the elements of the initial set can have a different (higher) dimensionality to the elements of the final predicted set, increasing the flexibility of the model.
Moreover, instead of concatenating a presence variable to each set point, \gls{TSPN} uses a \gls{MLP} to directly predict the number of points from the input embedding $\bm{h}$, which then chooses the size $N$ of the initial set $\set{\bx}_{i=1}^N$.
Given this size, we must generate the elements of the initial set such that the final predicted set is permutation-equivariant. One solution is to sample the intial set from a permutation-invariant distribution \eg \textit{iid} samples from a fixed distribution
, for which we choose $\gauss{\alpha\bm{1}, \operatorname{diag}(\beta\bm{1})}$ with learnable $\alpha \in \RR$, $\beta \in \RR_+$.
This makes the model stochastic, but also allows one to choose an arbitrary size for the generated set, even to sizes that lie outside the range of set sizes encountered in training.
Formally, let $\bm{h}$ be the input embedding; the generating process reads as follows:
\begin{align}
    N &= \operatorname{\textsc{mlp}}(\bm{h})\,,\\
    \bx_i \sim&~\gauss{\alpha\bm{1}, \operatorname{diag}(\beta\bm{1})}, \quad i=1,\ldots,N\,,\\
    \widehat{\bx}_i &= \operatorname{\mathrm{concatenate}}(\bx_i, \by)\,,\\
    \set{\bx} &= \operatorname{\mathrm{transformer}}\left(\set{\widehat{\bx}}\right)\,.%
\end{align}
While training \gls{TSPN}, we use the ground-truth set-cardinality to instantiate the initial set, and we separately train the \gls{MLP} by minimizing categorical cross-entropy with the ground-truth set sizes.
The cardinality-\gls{MLP} is used only at test-time.
Note that, in contrast to \gls{DSPN}, \gls{TSPN} does not require any additional regularization terms applied to its representations.

We describe our work in relation to \gls{DSPN} \cite{Zhang2019dspn}, but we note that there are concurrent works that share the ideas presented here.
Both \textsc{detr}~\cite{Carion2020detr} and Slot Attention~\cite{Locatello2020slot} use a variant of the transformer \cite{Vaswani2017transformer} for predicting a set of object properties.
\textsc{Detr} uses an object-specific query initialization and is, therefore, not equivariant to permutations, similarly to \cite{Zhang2019dspn}.
Slot Attention is perhaps the most similar to our work---transforms a randomly sampled point-cloud, same as \gls{TSPN}, but it uses attention normalized along the query axis instead the key axis.
\section{Experiments}
\label{sec:experiments}
We apply \gls{TSPN} and our implementations of \gls{DSPN} and a \gls{CDSPN} to two tasks: point-cloud prediction on \gls{SETMNIST} and object detection on \gls{CLEVR}.
Point-cloud prediction is an autoencoding task, where we use a set encoder to produce vector embeddings of point clouds of varying sizes, and explore the use of the above set prediction networks for reconstructing the point clouds conditioned on these embeddings.
Object detection, instead, requires predicting sets of bounding boxes conditioned on image features.
While generating point-clouds requires predicting large numbers of points, which are often assumed to be conditionally independent of each other, detecting objects typically requires generating much smaller sets.
Due to possible overlaps between neighbouring bounding boxes, however, it is essential to take relations between different objects into account.
We implemented all models in \textsc{jax} \cite{Jax2018github} and \textsc{haiku} \cite{Haiku2020github} and run experiments on a single Tesla V100 GPU.
Models are optimized with \textsc{adam} \cite{Kingma2015adam} with default $\beta_1, \beta_2$ parameters.
We used batch size $=32$ and trained all models for $100$ epochs.

\subsection{Point-Cloud Generation on \textsc{Set-Mnist}}
\label{sec:set_mnist}

    \begin{figure}
        \centering
        \includegraphics[width=\linewidth]{figs/set_mnist.png}
        \caption{Set reconstructions (rows 1-5) and inputs (row 6) from different models, annotated with the number of points (above).
        \Gls{DSPN} fails to learn set cardinality, see \Cref{sec:loss} for an explanation of why this happens.
        \Gls{CDSPN} learns the cardinality well, but produces reconstructions of lower fidelity than \gls{TSPN}.
        }
        \label{fig:mnist}
    \end{figure}

    \begin{table}
        \centering

        \begin{tabular}{c|c|c}
            \toprule
            Model & Chamfer $\left[10^{-5}\right]$ $\downarrow$ & Set Size \textsc{Rmse} $\downarrow$\\
            \midrule
            \textsc{dspn} & $8.2 \pm 0.66$ & $165 \pm 54.7$\\
            \textsc{dspn-big} & $8.5 \pm 1.83$ & $182 \pm 29.2$\\
            \textsc{c-dspn} & $8.92 \pm 0.78 $ & \bm{$0.3 \pm 0.13$}\\
            \textsc{c-dspn-big} & $9.63 \pm 1.17$ & \bm{$0.3 \pm 0.09$}\\
            \textsc{tspn}& \bm{$5.42 \pm 0.15$} & $0.8 \pm 0.08$\\
            \bottomrule
        \end{tabular}

        \caption{
            \Gls{SETMNIST} results, averaged over 5 runs.
            \textsc{Rmse} is root-mean-square error.
            Note that our \gls{DSPN} results match the Chamfer loss reported in \cite{Zhang2019dspn}.
            \textsc{C-dspn} stands for the \gls{DSPN} conditioned on the number of points.
            The \textsc{-big} variants use wider layers with the same number of parameters as \gls{TSPN}.
            For conditional models, sets were generated with the number of points inferred by the model.
        }
        \label{tab:set_mnist}
    \end{table}

We convert \gls{MNIST} into point clouds by thresholding pixel values with the mean value in the training set; then normalize coordinates of remaining points to be in $\interval{0, 1}$.

We use the same hyperparameters for \gls{DSPN} as reported in \cite{Zhang2019dspn}.
\gls{CDSPN} uses the same settings, except instead of presence variables $p_i$ it uses an \gls{MLP} with one hidden layer of size $128$ to predict the set size, as used in \gls{TSPN}.
\Gls{TSPN} uses a three-layer Transformer, with parameters shared between layers\footnote{This gives 383k parameters compared to 190k for \gls{DSPN}, which shares parameters between its encoder and gradient-based decoder. Not sharing parameters between layers does not improve results, but significantly increases parameter count. Sharing parameters and increasing layer width to match the number of parameters does not increase performance, either.}.
Each layer has $256$ hidden units and $4$ attention heads.
The initial set is two-dimensional (same as the output), but it is linearly-projected into $256$-dimensional space.
The outputs of the transformer layers are kept at $256$ dimensions, too, and we project them back into $2$ dimensions after the last layer.
\Gls{TSPN} uses the same input encoder as (\textsc{c}-)\gls{DSPN}, that is a two-layer \textsc{deepset} with \textsc{fspool} pooling layer \cite{Zhang2020fspool}.
All models are trained using the Chamfer loss with learning rate~$=10^{-3}$.

\Cref{tab:set_mnist} shows quantitative results.
Conditioning on the set size does not improve the Chamfer loss for \gls{CDSPN} but it does significantly improve the accuracy of set-size prediction.
Further, replacing the decoder with the Transformer (\Gls{TSPN}) leads to a significant loss reduction with respect to \gls{DSPN}.
\Cref{fig:mnist} shows inputs set and model reconstructions.
Note that in our experiments, \gls{DSPN} almost always predicts the same number of points ($323$), which is close to the maximum number of points we used (here, $342$).
Interestingly, \gls{TSPN} performs very well while extrapolating to much bigger sets than the ones encountered in training: \Cref{fig:tspn_mnist_extrapolation} in \Cref{sec:mnist_app} shows reconstructions, where we manually change the desired set size up to $1000$ points.
This is in contrast to \gls{CDSPN}, whose performance decreases significantly when we require it to generate a set whose size differs only slightly from the input, \cf \Cref{fig:dspn_mnist_extrapolation} in \Cref{sec:mnist_app}.
We conjecture that this is caused by how \textsc{fspool} handles sets of different sizes, which leads to incompatibility between embeddings of sets of different cardinality.
This, in turn, causes the distances in the latent space to be ill-defined, and breaks the internal gradient-based optimization.

\subsection{Object Detection on \textsc{Clevr}}

      \begin{figure}
        \centering
        \includegraphics[width=\linewidth]{figs/clevr_obj_det.png}
        \caption{
            Object detection on \gls{CLEVR}, with the number of detected objects in the lower-left corner of each image.
            \Gls{DSPN} fails to detect the correct number of objects.
            The ground-truth bounding boxes are approximate.
        }
        \label{fig:clevr}
        \vspace{-1em}
    \end{figure}

       \begin{table*}
        \centering

        \begin{tabular}{c|c|c|c|c|c|c|c}
            \toprule
            Model & $\textsc{ap}_{50}$ & $\textsc{ap}_{95}$ & $\textsc{ap}_{95}$ & $\textsc{ap}_{98}$ & $\textsc{ap}_{99}$ & Set Size \textsc{Rmse} $\downarrow$ &  \#Params\\
            \midrule
            \textsc{dspn} & $67.7 \pm 5.49$ & $7.4 \pm 0.91$ & $0.6 \pm 0.10$ & $0.0 \pm 0.01$ & $0.0 \pm 0.00$ & $2.53 \pm 0.221$ & $0.3$\,M\\
            \textsc{c-dspn} & $71.6 \pm 3.40$ & $10.8 \pm 1.50$ & $0.9 \pm 0.21$ & $0.0 \pm 0.01$ & $0.0 \pm 0.00$ & $1.74 \pm 0.301$ & $0.3$\,M\\
            \textsc{tspn}& \bm{$81.2 \pm 1.03$} & \bm{$20.7 \pm 0.16$} & \bm{$3.0 \pm 0.20$} & \bm{$0.1 \pm 0.02$} & $0.0 \pm 0.00$ & \bm{$0.58 \pm 0.046$} & $1.9$\,M\\

            \bottomrule
        \end{tabular}
        \caption{
            \Gls{CLEVR} object detection results averaged over 5 runs.
            These results are not directly comparable to the ones reported in \citet{Zhang2019dspn}, since all our models were trained using Chamfer loss. We subtract the resnet parameters ($21.8$\,M) when reporting the number of model parameters.
        }
        \label{tab:clevr}
    \end{table*}

\Gls{CLEVR} images consist of up to 10 rendered objects on plain backgrounds.
Following \citet{Zhang2019dspn}, we use the \gls{CLEVR} dataset to test the efficacy of our models for object detection in a simple setting, which might pave the way for more advanced object-based inference in future works.
We use the same hyperparameters for \gls{CDSPN} as \citet{Zhang2019dspn}.
\Gls{TSPN} uses four layers with four attention heads and $256$ neurons each, without parameter sharing between layers.
We apply layer normalization \cite{Xiong2020Oonln} before attention as in \cite{Ba2016ln}.
The last transformer layer is followed by an \gls{MLP} with a single hidden layer of %
All models use a \textsc{resnet34} \cite{He2016resnet} as the input encoder, and are trained with the Chamfer loss; \gls{DSPN}-based models are trained for $200$ epochs with learning rate~$=3\times10^{-5}$; \gls{TSPN} uses learning rate~$=10^{-4}$ and is trained for $1200$ epochs.
Longer training for \gls{DSPN}-based models lead to overfititng and decreased validation performance.
Note that \gls{DSPN} in \cite{Zhang2019dspn} use the Hungarian loss (\Cref{eq:hungarian}), which leads to better results than using Chamfer.
We report the average precision scores at different thresholds and the set size root-mean-square error in \Cref{tab:clevr}.
Qualitative results are available in \Cref{fig:clevr}.
We see that \gls{TSPN} outperforms \gls{CDSPN} and produces bounding boxes that are better aligned with ground-truth, although these results are worse than the ones reported in \citet{Zhang2019dspn}---we expect improvements when using the Hungarian loss as well.
\section{Related Work}
\label{sec:related_work}

\section{Conclusions}
\label{sec:conclusion}
We introduced the \glsreset{TSPN}\gls{TSPN}---a transformer-based model for conditional set prediction.
\Gls{TSPN} infers the cardinality of the set, randomly samples an initial set of the desired size and applies a transformer to generate the final prediction.
Set prediction in \gls{TSPN} is permutation-equivariant, and the model can be applied to any set-prediction tasks.
Interesting directions include scaling the model to large-scale point-clouds and object detection (\eg similar to \citet{Carion2020detr}), as well as turning this model into a generative model in either the \textsc{vae} or \textsc{gan} framework.

\section*{Acknowledgements}
We would like to thank George Papamakarios, Karl Stelzner, Thomas Kipf, Teophane Weber and Yee Whye Teh for helpful discussions.

\section{Set Size Extrapolation on \textsc{Set-Mnist}}
\label{sec:mnist_app}

\begin{figure*}[b]
    \centering
    \begin{minipage}[t]{.49\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figs/tspn_mnist_extrapolation.pdf}
        \caption{
            \Gls{TSPN} extrapolates to far bigger set sizes than encountered in training.
            Here, the model was trained with up to 342 points, and yet can generate sets of up to 1000 points.
            The bottom row contains the ground-truth annotated with the number of points at the bottom.
            This figure uses a smaller marker size than \Cref{fig:mnist}, hence ground-truth appears less dense.
        }
        \label{fig:tspn_mnist_extrapolation}
    \end{minipage}
    \hfill
    \begin{minipage}[t]{.49\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figs/dspn_mnist_extrapolation.pdf}
        \caption{
            The size-conditional \gls{DSPN} fails to reconstruct sets when we slightly change the output set size.
            Here, we take the size of the reconstructed set to be within $\interval{-4, 5}$ points of the input size.
            As the size of the reconstruction strays from the original size, the reconstruction quality quickly deteriorates.
            The bottom row contains the input sets with their respective sizes.
        }
        \label{fig:dspn_mnist_extrapolation}
    \end{minipage}
\end{figure*}


\section{Introduction}

Macroscopic systems, in absence of an external drive, equilibrate with
the environment.
However, relaxation may be slow,
i.e. with a relaxation time which exceeds any attainable observation time~\cite{Palmer}.
In that case,
only dynamical properties are accessible to observation
and the question naturally arises of what can
be learnt about equilibrium from dynamics. Paradigmatical examples of slow relaxation
are glassy systems~\cite{BCKM} or systems undergoing
phase ordering after a sudden temperature quench from above
to below the critical point~\cite{BCKM,Bray}.
Here we shall look at the problem in the latter
context, whose prototypical instance is the quench of a
ferromagnetic system. In order to make the
presentation as simple as possible, we shall mostly concentrate on
the Ising model. The extension to other phase-ordering systems will be
discussed at the end of the paper, with the example of the spherical model.

Phase ordering in the Ising model by now is a mature subject, generally considered
to be well understood. For reviews see Refs.~\cite{Bray,Puri,Jo,Henkel}.
Among the many interesting features of the process, in this paper we
shall be primarily concerned with the lack of equilibration in any finite time,
if the system is infinite. This is frequently referred
to with the catchy expression that the system remains permanently out of equilibrium,
whose meaning, however, has never been fully clarified.
For instance, a similar circumstance arises also when the quench is made to the
critical temperature $T_c$, because, due to critical slowing down,
again equilibrium is not reached in any
finite time. Nonetheless, in that case, the process cannot be regarded as
substantially different from one of equilibration,
because as time grows the system gets
closer and closer to the equilibrium critical state, which is unique in the sense that
in the thermodynamic limit it is independent of the boundary conditions (BC).
Instead, in the quench to below $T_c$ the picture is qualitatively different,
because, although the state extrapolated from dynamics is unique, the same cannot
be said of the equilibrium state, which
depends on BC even in the thermodynamic limit. This we have shown
in Ref.~\cite{FCZ} (to be referred
to as I in the following), where we have investigated
the nature of the equilibrium state in the Ising model below
$T_c$, under different symmetry-preserving BC.
We have found that while periodic boundary conditions (PBC)
lead to the usual ferromagnetic ordering, due to
the breaking of ergodicity with the consequential spontaneous breaking
of the $\mathbb{Z}_2$ up-down symmetry,
the scenario changes dramatically with antiperiodic boundary conditions
(APBC), because
ergodicity breaking is precluded. Then,
the system cannot order and complies with the requirement of the transition
by remaining critical also below $T_c$,  all the way down to $T=0$.
We have argued that this new transition, without spontaneous symmetry breaking
and without ordering,
consists in the condensation of fluctuations. In the $1d$ case, since $T_c=0$, the low temperature
phase is shrunk to just $T=0$.

Motivated by the existence of such diversity in the equilibrium properties,
in this paper we address
the next natural question, formulated in the title of the paper, of matching
statics and dynamics.
Using the equal-time correlation function as the probing observable, we shall see
that the asymptotic state, extrapolated
from dynamics, that is by taking the
$t \to \infty$ limit after the thermodynamic limit, is unique and {\it critical}.
Now, the point is that this, which we may call the time-asymptotic state
and which, we emphasize, is the same for both choices of BC,
is found to coincide with the bona fide equilibrium state, i.e. the one
computed from equilibrium statistical mechanics, in the APBC case but
to be remote from it in the PBC case.
Thus, we have the one and the same
dynamical evolution which, although not reaching equilibrium in any finite time,
turns out to be informative of the true equilibrium state in one case
(APBC), but not in the other (PBC).
It is, then, appropriate to regard the APBC case as one in which
equilibrium is approached,
just as in a quench to $T_c$, while the PBC case offers an instance of a system remaining
permanently out of equilibrium. The poor performance in approaching equilibrium
with PBC is traceable to the presence
of ergodicity breaking at the working temperature, which, instead, is preserved when
APBC are applied.
At the end of the paper we shall argue that the connection
between the presence/absence of ergodicity breaking and the absence/presence of equilibration goes beyond the Ising example, by showing that it takes place
with the same features in the rather different context of
the spherical and mean-spherical model.

The paper is organized as follows: in section~\ref{II} we formulate the problem.
In section~\ref{III}  the relation between equilibrium and relaxation in the quench
to above $T_c$ is analyzed by using scaling
arguments. The cases of the quench to $T_c$ with $d=2$, to $T=0$ with $d=1$
and to below $T_c$ with $d=2$ are analyzed in sections~\ref{IV}, \ref{V} and
\ref{VI}, respectively. The spherical and mean spherical model are introduced and
investigated in section~\ref{SM}. Concluding remarks are made in section~\ref{CR}.

\section{The problem}
\label{II}

We are concerned with the relaxation dynamics of a system initially prepared in
an equilibrium state at the temperature $T_I$ and suddenly quenched to the
lower temperature $T_F$. We consider
the Ising model on a lattice of size $V=L^d$, with the usual nearest neighbours
interaction
\begin{equation}
{\cal H}(\boldsymbol{s}) = -J\sum_{<ij>} s_is_j,
\label{Ham.1}
\end{equation}
where $J > 0$ is the ferromagnetic coupling,
$\boldsymbol{s} = [s_i]$ is a configuration of spin variables $s_i=\pm1$
and $<ij>$ is a pair of nearest neighbours.
We shall study the $d=1$ and $d=2$ cases, where
in the thermodynamic limit there is a critical point at $T_c=0$ and
$T_c=2.269J$, respectively.
Since the system's size is finite, BC must be specified and, because of the major role
that these will play in the following developments, it is necessary to enter in some
detail from the outset.
As anticipated in the
Introduction, we shall consider PBC and APBC (precisely cylindrical
antiperiodic BC) implemented by adding to the interaction an extra term
${\cal B}(\boldsymbol{s})$ with couplings among spins on the boundary~\cite{Gallavotti,Antal,FCZ}.
In the $d=2$ case spins on opposite edges are coupled ferromagnetically, just
like spins in the bulk, if PBC are applied.
Instead, in the APBC case, spins on one pair of opposite edges
are coupled ferromagnetically, while those on the other pair antiferromagnetically. Hence,
the boundary term reads
\begin{equation}
{\cal H}_{b}(\boldsymbol{s}) =
-J \sum_{y=1}^L s_{1,y} s_{L,y} - b J \sum_{x=1}^L s_{x,1}s_{x,L},
\label{bdr.1}
\end{equation}
where we have denoted by
$b=\pm$ the sign of $J$, which identifies
PBC $(+)$ or APBC $(-)$.
In the $d=1$ case this term simplifies to
 \begin{equation}
{\cal H}_b (\boldsymbol{s}) =
- bJ s_{1}s_{L},
\label{bdr.2}
\end{equation}
where $L$ is the length of the chain. It is important to note that both these BC preserve
the up-down symmetry of the Ising interaction.

Taking as it is customary $T_I = \infty$ in order to have an uncorrelated initial state,
the system is put in contact with a thermal reservoir at the lower and
finite temperature $T_F$
and let to evolve according to a dynamical rule which does not conserve the
order parameter, like Glauber or Metropolis. This simply corresponds to
running a Markov chain at the fixed temperature $T_F$, with the so-called hot start,
that is with a uniformly random initial condition.
The relaxation process is monitored through the equal-time spin-spin correlation function
\begin{equation}
\mathcal{C}(r, \epsilon,t^{-1},L^{-1}; b)  =   \left [ \langle s_{i}(t)s_{j}(t)\rangle -
\langle s_{i}(t) \rangle \langle s_{j}(t) \rangle \right ],
\label{corr.0}
\end{equation}
where the angular brackets denote averages taken over the noisy dynamics
and the initial conditions, while the square brackets stand for the average over all
pairs of sites $(i,j)$ keeping fixed the distance $r$ between $i$ and $j$.
In the set of control parameters, $\epsilon=T_F-T_c$ is the temperature difference from
criticality, $t^{-1}$ is the inverse time and $L^{-1}$ is the inverse linear size.

We are interested in taking both
the large-time and the thermodynamic limit
of the above quantity, and then to compare the outcomes,
depending on the order in which these two limits have been taken. Letting
$t^{-1} \to 0$ first, while keeping $L$ fixed,
the equilibrium correlation function is obtained
\begin{equation}
\lim_{t^{-1} \to 0}\mathcal{C}(r,\epsilon,t^{-1},L^{-1};b)  =
\mathcal{C}_{\rm eq}(r,\epsilon,L^{-1};b)  =
\left [ \langle s_{i}s_{j}\rangle_{\rm eq} -
\langle s_{i} \rangle_{\rm eq} \langle s_{j} \rangle_{\rm eq} \right ],
\label{D3.7}
\end{equation}
where now the angular brackets stand for the Gibbs ensemble average
and the square brackets have the same meaning as in
Eq.~(\ref{corr.0}). Then, the subsequent
thermodynamic limit implements the prescription~\cite{Gallavotti} for the construction of
the equilibrium correlation function in the infinite system
\begin{equation}
\lim_{L^{-1} \to 0} \lim_{t^{-1} \to 0}\mathcal{C}(r,\epsilon,t^{-1},L^{-1};b) =
C_{\rm eq}(r,\epsilon;b).
\label{D3.6}
\end{equation}
The crux of the matter is that, after reversing the order of these
limits, the end result might not be the same as the one above,
because the large-time limit of the
time-dependent correlation function for the infinite system
\begin{equation}
\lim_{t^{-1} \to 0} \lim_{L^{-1} \to 0} \mathcal{C}(r,\epsilon,t^{-1},L^{-1};b) =
C^*(r,\epsilon;b),
\label{D3.5}
\end{equation}
exists but does not necessarily coincide with
$C_{\rm eq}(r,\epsilon;b)$.
Referring to $C^*(r,\epsilon;b)$ as the time-asymptotic correlation
function, if it matches $C_{\rm eq}(r,\epsilon;b)$ then
the infinite system equilibrates. If not, it remains permanently out
of equilibrium.
Which is the case depends on $T_F$ and on the choice of BC.
In the
quench to $T_F \geq T_c$ both $C_{\rm eq}(r,\epsilon)$ and
$C^*(r,\epsilon)$ are independent of the BC choice and do coincide,
signaling equilibration. Instead, in the
quench to below $T_c$, as we shall see,
$C^*(r,\epsilon)$ does not depend on
$b$, while $C_{\rm eq}(r,\epsilon;b)$ retains this dependence,
implying that equilibration can be achieved at most with one of the two BC, but
certainly not with both. As anticipated
in the Introduction, the equilibration condition is fulfilled with APBC, but
not with PBC.

In the next section we shall substantiate the above statements with results for the
$1d$  and $2d$ Ising model.
We shall take the aforementioned limits, after setting up the
general scaling scheme which unifies static and dynamic phenomena
into one single framework encompassing both.
In order to do this, it is convenient to treat separately the three cases:
$\epsilon > 0$, $\epsilon = 0$ and $\epsilon < 0$.

\section{Statics and Dynamics: $\boldsymbol{\epsilon > 0}$}
\label{III}

\begin{figure}[ht]
\centering
\includegraphics[width=8cm]{FigPPT1SXa.pdf} \hspace{-2cm}
\includegraphics[width=8cm]{FigPPT1DXb.pdf}
        \caption{Parameter space of the $1d$ model (a)
        and of the $2d$ model (b).}
\label{fig:pspace}
\end{figure}

Let us assume that at a generic point in the $\epsilon > 0$ sector
of the three-dimensional space of the parameters
$(\epsilon,t^{-1},L^{-1})$, depicted in Fig.\ref{fig:pspace},
the correlation function obeys scaling in the form~\footnote{This is a finite-size
extension of the scaling form derived in Ref.~\cite{Janssen}.}
\begin{equation}
\mathcal{C}(r,\epsilon,t^{-1},L^{-1};b)
= \frac{1}{r^a} \mathcal{F} \left ( \frac{r}{\ell}, \frac{\ell}{R},\frac{R}{L};b \right ),
\label{anml.1}
\end{equation}
where the exponent $a$ is related to the anomalous dimension exponent $\eta$
by $a = d-2+\eta$ and to the fractal dimensionality $D$ of the
Coniglio-Klein (CK)~\cite{CK,FK} correlated clusters~\cite{CF} by
\begin{equation}
a = 2(d-D).
\label{anml.2}
\end{equation}
From the exact results~\cite{Stanley,Goldenfeld}
\begin{equation}
\eta = \left \{ \begin{array}{ll}
         1, \;\; $for$ \;\; d=1,\\
         1/4,\;\; $for$ \;\; d=2,
        \end{array}
        \right .
        \label{z.0}
        \end{equation}
follows
\begin{equation}
a = \left \{ \begin{array}{ll}
         0, \;\; $for$ \;\; d=1,\\
         1/4,\;\; $for$ \;\; d=2,
        \end{array}
        \right .
        \label{z.1}
        \end{equation}
and
 \begin{equation}
D = \left \{ \begin{array}{ll}
         1, \;\; $for$ \;\; d=1,\\
         15/8,\;\; $for$ \;\; d=2,
        \end{array}
        \right .
        \label{z.2}
        \end{equation}
which shows that the CK clusters are compact in $1d$ and fractals in $2d$.
Up to a proportionality constant, the scaling variable $\ell$
is the equilibrium correlation length of the infinite system, given by~\cite{Stanley}
\begin{equation}
\ell(\epsilon)
= \left \{ \begin{array}{ll}
         -[\ln \tanh (J/\epsilon)]^{-1}, \;\; $for$ \;\; d=1,\\
         \epsilon^{-\nu},\;\; $with$ \;\; \nu=1 \;\; $for$ \;\; d=2.
        \end{array}
        \right .
        \label{anml.3}
        \end{equation}
The other characteristic length $R(t)$ obeys the power law~\cite{Janssen}
\begin{equation}
R(t) = t^{1/z},
\label{anml.3bis}
\end{equation}
with the dynamical exponent~\cite{1d,Nightingale}
\begin{equation}
z
= \left \{ \begin{array}{ll}
         2, \;\; $for$ \;\; d=1,\\
         2.16,\;\; $for$ \;\; d=2.
        \end{array}
        \right .
        \label{z.1}
        \end{equation}
The connection between $R(t)$ and the time dependent correlation length will be
clarified shortly and is summarized in Fig.\ref{fig:R}.
Both lengths diverge as the critical point, which is at the origin of the
reference frame in Fig.\ref{fig:pspace}, is
approached along the $\epsilon$ axis and the $t^{-1}$ axis, respectively.

The scaling ansatz~(\ref{anml.1}) is dense of information and allows to predict what
should be expected in different regions of the parameter space. The foremost relevant
features are the power-law decay $1/r^{a}$ of correlations at
short distance and the large-distance cutoff enforced by the scaling
function. The separation between short and large distances is fixed by the
correlation length
\begin{equation}
\xi(\epsilon,t^{-1},L^{-1};b) =
\left [ \frac{\int_0^L dr \, r^2 \, \mathcal{C}(r,\epsilon,t^{-1},L^{-1};b) }
{\int_0^L dr \,  \mathcal{C}(r,\epsilon,t^{-1},L^{-1};b) } \right ]^{1/2},
\label{corrl.0}
\end{equation}
which scales as
\begin{equation}
\xi(\epsilon,t^{-1},L^{-1};b) = R f \left ( \frac{R}{\ell}, \frac{\ell}{L};b \right ).
\label{crrl.1}
\end{equation}
The behaviour of $\xi$, as parameters are changed,
can be unraveled by the following argument. Suppose that $\ell$ and $L$ are fixed
in a region where $\ell \ll L$ and let us survey what happens
as the quench unfolds and $R$ grows.
\begin{figure}[ht]
\centering
\includegraphics[width=8cm]{FigPPT2SXa.pdf} \hspace{-2cm}
\includegraphics[width=8cm]{FigPPT2DXb.pdf}
        \caption{Schematic representation of the
        saturation of $\xi$ vs $R$ for $\ell \ll L$ (a)
        and for $\ell \gg L$ (b).}
\label{fig:R}
\end{figure}
Approximating the above equation by
\begin{equation}
\xi(\epsilon,t^{-1}) \simeq R f \left ( \frac{R}{\ell},0 \right ),
\label{crrl.101}
\end{equation}
in the early stage of the quench, when $R \ll \ell$, it can be further reduced to
\begin{equation}
\xi \sim R,
\label{xiR.1}
\end{equation}
because the system behaves as if it was approaching the critical
point along the $t^{-1}$ axis. As $R$ is let to grow further, equilibrium is
eventually reached
when $R \sim \ell$ and the correlation length saturates to the limiting
value
\begin{equation}
\xi \sim \ell,
\label{elle.1}
\end{equation}
as illustrated in the left panel of Fig.\ref{fig:R}, with the equilibration time
given by $t_{\rm eq} = \ell^z$.
The BC are immaterial throughout, because $\xi$ remains always much smaller than $L$
so that the system as a whole behaves as a collection of independent
finite systems, on which the far away BC have no effect.
By the same reasoning, in the region where $\ell \gg L$ we still have $\xi \sim R$
in the early stage, when $R \ll L$, with independence from BC. But then BC come into
play when the system equilibrates and
$\xi$ saturates to the limiting value
$L$, as illustrated in the right panel of Fig.\ref{fig:R},
since correlations extend up to distances where the BC are effective.
In this connection see Ref.~\cite{Das}.
Summarizing, $\xi$ is given by the shortest of the three lengths
$(\ell, R,L)$, that is
\begin{equation}
\xi(\epsilon,t^{-1},L^{-1}) \sim \min(\ell,R,L),
\label{crrl.4}
\end{equation}
in the regions of the parameter space where one of the three is considerably
shorter than the other two,
with crossovers connecting these regions.
It is clear from Fig.\ref{fig:R} that $\xi$ and $R$ coincide at
all times if both $\ell$ and $L$ are infinite.
Next to the correlation length, it is useful to keep track also of the
susceptibility
\begin{equation}
\chi(\epsilon,t^{-1},L^{-1};b) =
\int d^d r \,  \mathcal{C}(r,\epsilon,t^{-1},L^{-1};b),
\label{susc.0}
\end{equation}
which is related to the correlation length by
\begin{equation}
\chi \sim \xi^{2D-d}.
\label{susc.1}
\end{equation}
This is an important relation, because it is independent of the direction of approach
to the critical point and depends only on the geometrical nature of the correlated clusters
through $D$.

According to the above reasoning, when the limits $t^{-1} \to 0$ and $L^{-1} \to 0$
are taken in the $\epsilon > 0$ sector, we necessarily have $\xi \sim \ell$,
independently of the order in which these limits are taken, because $\ell$
is finite. Moreover, the finite correlation
length guarantees that the system equilibrates
with independence from BC
\begin{equation}
\lim_{L^{-1} \to 0} \lim_{t^{-1} \to 0}
\mathcal{F} \left ( \frac{r}{\ell}, \frac{\ell}{R},\frac{R}{L};b \right ) =
\lim_{t^{-1} \to 0} \lim_{L^{-1} \to 0}
\mathcal{F} \left ( \frac{r}{\ell}, \frac{\ell}{R},\frac{R}{L};b \right ) =
F_{\rm eq} \left (\frac{r}{\ell} \right ).
\label{D3.6bis}
\end{equation}

\subsection{Example: $1d$ system}

As an example, let us check the above statements against exact results
in the particular case of the
$t^{-1} \to 0$ limit of the $1d$ model with finite $L$.
The equilibrium correlation function is given by
\begin{equation}
\mathcal{C}_{\rm eq}(r,\epsilon,L^{-1};b)
= \frac{1}{r^a} \mathcal{F}_{\rm eq} \left ( \frac{r}{L}, \frac{L}{\ell};b \right ),
\label{anml.4}
\end{equation}
where $a=0$ according to Eq.~(\ref{z.1}), while
the two explicit forms of the scaling function (see I and Ref.~\cite{Antal}) read
\begin{equation}
\mathcal{F}^{(p)}_{\rm eq}(z,\zeta) = \frac{\cosh[\zeta(1-z)]}{\cosh(\zeta)},
\label{anml.6}
\end{equation}
\begin{equation}
\mathcal{F}^{(a)}_{\rm eq}(z,\zeta) = \frac{\sinh[\zeta(1-z)]}{\sinh(\zeta)},
\label{anml.7}
\end{equation}
where we have set
\begin{equation}
z=r/L, \quad \zeta = L/\ell.
\label{anml.8}
\end{equation}
We have considered a chain of length $2L$ in order to simplify notation.
The superscripts $(p)$ and $(a)$ have been used for PBC and for APBC, respectively.
The equilibrium correlation length, defined through the second moment as
in Eq.~(\ref{corrl.0}),
scales as
\begin{equation}
\xi_{\rm eq}(\epsilon,L^{-1};b) = \ell f_{\rm eq}(\zeta;b),
\label{crrl.1bis}
\end{equation}
with the scaling functions
\begin{equation}
 f_{\rm eq}^{(p)} (\zeta) =  \left [2-\frac{2\zeta}{\sinh (\zeta)} \right ]^{1/2},
 \quad f_{\rm eq}^{(a)} (\zeta) =  \left [2-\frac{\zeta^2}{\cosh(\zeta) - 1} \right ]^{1/2},
\label{crrl.5}
\end{equation}
from which follows
\begin{equation}
\xi_{\rm eq}^{(p)}(\epsilon,L^{-1}) =  \left \{ \begin{array}{ll}
         \sqrt{2} \, \ell, \;\; $for$, \;\; \ell \ll L,\\
         \frac{1}{\sqrt{3}} \, L,\;\; $for$, \;\; L \ll \ell,
        \end{array}
        \right .
        \quad
\xi_{\rm eq}^{(a)}(\epsilon,L^{-1}) =  \left \{ \begin{array}{ll}
         \sqrt{2} \, \ell, \;\; $for$, \;\; \ell \ll L,\\
         \frac{1}{\sqrt{6}} \, L,\;\; $for$, \;\; L \ll \ell,
        \end{array}
        \right .
        \label{crrl.10}
        \end{equation}
showing that in the regimes $\ell \ll \L$ and $\ell \gg L$, indeed
one has $\xi \sim \min(\ell,L)$.
Completing, next, the sequence of limits by letting
$L^{-1} \to 0$, it is straightforward to check that the dependence on BC disappears,
yielding
\begin{equation}
\lim_{L^{-1} \to 0} \mathcal{F}^{(p)}_{\rm eq}(z,\zeta)
= \lim_{L^{-1} \to 0} \mathcal{F}^{(a)}_{\rm eq}(z,\zeta) = e^{-r/\ell}.
\label{anml.9}
\end{equation}
Using the definition, it is immediate to verify that also $\chi \sim \ell$ and,
therefore, that moving toward the critical point along the $\epsilon$ axis one has
\begin{equation}
\chi \sim \xi,
\label{susc.2}
\end{equation}
in agreement with Eq.~(\ref{susc.1}), because $D=1$ when $d=1$.

\section{Statics and Dynamics: $\boldsymbol{\epsilon = 0, d=2}$}
\label{IV}

When $\epsilon = 0$, the $1d$ and $2d$ cases are quite different and need
to be treated separately. In the latter one, which we shall now consider,
$T_c > 0$ and ergodicity does not break. In the former, instead, $T_c = 0$
and ergodicity may break, depending on BC. This makes it more akin to
the $2d$ quench to below $T_c$. So, it will be dealt with in the next section.

\begin{figure}[ht]
\centering
\includegraphics[width=10cm]{fig3.pdf}
        \caption{Scaling function of the time-dependent correlation
        function in the quench to $T_c$ of the $2d$ model,
        demonstrating independence from BC, in the system with $L=256$.
        PBC (black symbols) and APBC (empty symbols).}
\label{fig:collapseTc}
\end{figure}

The specificity of the quench to $\epsilon = 0$ is that $\ell$ diverges and, consequently,
that $\xi$ can be limited only by $R$ or $L$. Thus, when the $t^{-1} \to 0$ limit is taken first
and $L$ is kept fixed,
$\xi$ crosses over from $R$ to $L$ in the finite time $t_{\rm eq} \sim L^z$,
as in the right panel of Fig.\ref{fig:R}, and the system equilibrates to
\begin{equation}
\mathcal{C}_{\rm eq}(r,L^{-1};b) = \frac{1}{r^{1/4}}
\mathcal{F}_{\rm eq} \left ( \frac{r}{L};b \right ),
\label{sat.1}
\end{equation}
which depends on BC because correlations extend up to the boundary.
Letting next $L^{-1} \to 0$, the BC dependence disappears from the
critical correlation function of the infinite system
\begin{equation}
C_{\rm eq}(r) \sim \frac{1}{r^{1/4}}.
\label{sat.2}
\end{equation}

Instead, if the $L^{-1} \to 0$ limit is taken first, $R$ is the only length
left in the problem. This implies $\xi \sim R$ at all times,  so that there is no
finite equilibration time. However, the time-dependent
correlation function
\begin{equation}
C(r,t^{-1}) = \frac{1}{r^{1/4}}  F_c (r/R),
\label{sat.3}
\end{equation}
which is BC independent (see Fig.\ref{fig:collapseTc}),
gets arbitrarily close to the equilibrium counterpart~(\ref{sat.2}) as time
grows, because the limit
\begin{equation}
\lim_{t^{-1} \to 0} C(r,t^{-1}) = C^*(r) \sim \frac{1}{r^{1/4}},
\label{sat.4}
\end{equation}
coincides with it.
In summary, in the quench to $\epsilon = 0$ of the $2d$ system, like in
the $\epsilon > 0$ case previously considered, the equilibrium correlation
function of the infinite system $C_{\rm eq}(r)$ does not depend on
BC and coincides with the time-asymptotic one $C^*(r)$,
warranting the conclusion that the system can get arbitrarily close to
equilibrium by waiting long enough.

Comparing Eqs.~(\ref{sat.1}) and~(\ref{sat.3}), it is evident that the scaling
structure is the same, the only difference being in the specific forms of the
scaling functions, which is inessential for the present considerations. This shows
that the time direction along the $t^{-1}$ axis, as far as scaling is concerned,
is just another direction
of approach to the critical point, on the same footing with the other two.
In addition, from the formal similarity of the two scaling expressions follows
straightforwardly that the susceptibility satisfies Eq.~(\ref{susc.1}) in the form
\begin{equation}
\chi \sim \xi^{7/4},
\label{suscett.1}
\end{equation}
irrespective of the direction of approach, with $\xi \sim L$ along the $L^{-1}$
axis and $\xi \sim R$ along the $t^{-1}$ axis.

\section{Statics and Dynamics: $\boldsymbol{\epsilon = 0, d=1}$}
\label{V}

As mentioned above and explained at length in I,
in the $1d$ system at $\epsilon = 0$ we are confronted with a radically different situation,
because ergodicity, which holds for both BC above $T_c$, is now broken with PBC,
but not with APBC. In order to ease the comparison, and to highlight
the contrast, with the less familiar case of a transition without ergodicity breaking,
let us first briefly summarize the well-established concept of ergodicity breaking~\cite{Palmer}. In the
PBC case there are two degenerate ground states: the two ordered configurations
with all spins either up $\boldsymbol{s}_+ = [s_i=+1]$ or
down $\boldsymbol{s}_- = [s_i=-1]$. These, by themselves, form
two absolutely-confining ergodic components, which are
dynamically disconnected because the activated moves  needed to go from one to the
other are forbidden at zero temperature. Consequently, time averages
coincide with ensemble averages taken with either one of the two
broken-symmetry ferromagnetic pure states
$P_-(\boldsymbol{s})=\delta_{\boldsymbol{s},\boldsymbol{s}_-}, P_+(\boldsymbol{s})=\delta_{\boldsymbol{s},\boldsymbol{s}_+}$
and {\it do not} coincide with the symmetric
ensemble averages taken in the Gibbs state, which is the
even mixture of the pure states
\begin{equation}
P^{(p)}(\boldsymbol{s}) = \frac{1}{2}[P_-(\boldsymbol{s}) + P_+(\boldsymbol{s})].
\label{mixt.0}
\end{equation}
In such a situation, only time averages are physically meaningful.
Conversely, in the APBC case all the $4L$ degenerate ground-state
configurations with one defect (or domain wall) belong to the same
ergodic component, because the defect can freely sweep
the whole system at no energy cost. Then, in this case
time and ensemble averages coincide.
The qualitative difference between the two zero-temperature states is
well illustrated (see Fig.\ref{fig:1}) by the
probability distribution $P_b(m)$ of the magnetization
density $m=\frac{1}{2L}\sum_i s_i$, which is demonstrated~\cite{Antal} to be
double peaked in the PBC case
\begin{equation}
P^{(p)}(m) = \frac{1}{2} [\delta(m+1) + \delta(m-1)],
\label{mxt.1}
\end{equation}
\begin{figure}[ht]
\centering
\includegraphics[width=5cm]{fig1a-sublabel.pdf}
\includegraphics[width=5cm]{fig1b-sublabel.pdf}
\caption{Magnetization density distributions at $\epsilon=0$ in the $1d$ Ising model,
	with $m_\pm = \pm 1$. The spikes in the panel (a) stand for $\delta$ functions.}
\label{fig:1}
\end{figure}
and uniform over the $[-1,1]$ interval in the APBC case
\begin{equation}
P^{(a)}(m) \to  \left \{ \begin{array}{ll}
         1/2, \;\; $for$, \;\; m \in [-1,1],\\
         0,\;\; $for$, \;\; m \notin [-1,1].
        \end{array}
        \right .
        \label{mxt.2}
        \end{equation}

So, if we now take the $t^{-1} \to 0$ limit while keeping $L^{-1}$ fixed,
we find  BC-related differences in the results.
With PBC, as explained above, the meaningful averages are those
in the broken symmetry pure states, yielding
\begin{equation}
\mathcal{C}^{(p)}_{\rm eq}(r,L^{-1}) = [\langle s_is_j \rangle_\pm - \langle s_i \rangle_\pm \langle s_j \rangle_\pm] = 0,
\label{1dI.0}
\end{equation}
where the angular brackets stand for
the average with respect  to $P_-(\boldsymbol{s})$ or $P_+(\boldsymbol{s})$.
The vanishing of correlations
for any $r$ holds independently of $L$ and clearly implies
that also the correlation length vanishes.
Notice that
if the correlation function had been extracted by taking the $\epsilon \to 0$ limit
of the Gibbs average taken with the distribution~(\ref{mixt.0}), as in
 Eq.~(\ref{anml.6}), the result would have been
\begin{equation}
\mathcal{C}^{(p)}_{\rm eq}(r,L^{-1}) = 1,
\label{1dI.16}
\end{equation}
which is independent of $r$ and does not decay. However, this would have been
just an artefact of the mixture.

Conversely, in the APBC case the ensemble Gibbs average, as calculated
in Eq.~(\ref{anml.7}), gives the correct time-average result
because ergodicity is not broken. Hence, in the $\epsilon \to 0$ limit,
from Eq.~(\ref{anml.7}) one has
\begin{equation}
\mathcal{C}^{(a)}_{\rm eq}(r,L^{-1}) = (1-r/L).
\label{1dI.17}
\end{equation}
The dependence on $r/L$ in the above expression reveals
that correlations extend over a distance of order L in agreement with the
general argument expounded in section~\ref{III}.
Hence, by letting $L^{-1} \to 0$  the correlation length diverges, leading
to the conclusion that
the state at the origin of the parameter space
is a {\it critical point} for the APBC system, where the correlation function
displays the constant behaviour
\begin{equation}
C^{(a)}_{\rm eq}(r) = 1.
\label{corrl.01}
\end{equation}
Contrary to Eq.~(\ref{1dI.16}), now the lack of decay is a real physical effect, which corresponds to the critical power law decay $1/r^{a}$ with a vanishing
exponent $a$, due to the compactness of the CK
correlated clusters.

When the sequence of limits is reversed, after taking the thermodynamic
limit we are again in the situation in which $R$ is
the only length in the problem. Therefore $\xi \sim R$, as in the previous section,
and we get the BC independent result
\begin{equation}
C(r,t^{-1}) = \frac{1}{r^a} F (r/R ),
\label{q.1}
\end{equation}
with $a=0$. The function $F(x)$
is known from exact analytical computation with
PBC~\cite{1d,Bray} and is given by
\begin{equation}
F(x) = \mathrm{erfc}\,(x), \quad \text{with} \quad x=r/2R.
\label{quench.1}
\end{equation}
That the same scaling function applies also to the case of APBC is demonstrated
by the numerical data displayed in Fig.~\ref{fig:scaledC}, which have been
obtained by simulating the quench dynamics with the
Metropolis algorithm on a system with $L=10^5$, after imposing PBC and APBC.
The plot shows that the above result indeed holds irrespective of the
BC choice, because the PBC and APBC data superimpose to the theoretical
curve of Eq.~(\ref{quench.1}) with great accuracy, as long as
$R(t) \ll L$.
The existence
of an endlessly growing correlation length $R(t)$ means that the relaxation dynamics
along the $t^{-1}$axis drives both systems, with PBC and with APBC,
toward the same asymptotic critical state at the origin
\begin{figure}[ht]
\centering
\includegraphics[width=10cm]{fig5.pdf}
\caption{Collapse on the master curve of Eq.~(\ref{quench.1}) of the data for
$C(r,t^{-1},L^{-1})$ in the time regime $R \ll L$. PBC (black symbols)
and APBC (empty symbols). The data for $R/L=10^{-4},10^{-3}$ have been obtained with
$L=10^5$, those with $R/L=10^{-2}$ with $L=10^4$.}
\label{fig:scaledC}
\end{figure}
with the unique time-asymptotic correlation function given by
\begin{equation}
\lim_{t^{-1} \to 0} C(r,t^{-1}) = C^*(r)= 1,
\label{elle.2}
\end{equation}
which coincides with the APBC equilibrium result in Eq.(\ref{corrl.01}).
So, if we compare the asymptotic result of Eq.~(\ref{elle.2})
with the APBC static one of Eq.~(\ref{corrl.01}),
and with the PBC equilibrium result of Eq.~(\ref{1dI.0}),
we see, as stated in the Introduction,
that the APBC system tends toward equilibrium, although with an infinite
relaxation time, while the PBC system remains permanently out of equilibrium.
It is evident that the origin of the diversity of behaviours is in the presence or
absence of ergodicity breaking. In fact, we shall see in the next section that
the same behaviour occurs in the quench of the $2d$ system to below the
critical point.

In order to complete the picture of critical behaviour, let us
check on the validity of Eq.~(\ref{susc.1}).  From
Eqs.~(\ref{anml.9},\ref{1dI.17},\ref{quench.1}) follows that along the three
directions one has $\xi \sim \ell, \xi \sim R, \xi \sim L$, as well as
$\chi \sim \ell, \chi \sim R, \chi \sim L$, yielding
\begin{equation}
\chi \sim \xi,
\label{susc.10}
\end{equation}
independently of the direction of approach to the critical point, as it should be
since $D=1$.

\section{Statics and Dynamics: $\boldsymbol{\epsilon < 0, 2d}$}
\label{VI}

As in the previous
case, the nature of the equilibrium state
of the $2d$ model below $T_c$ depends strongly on BC, even in the thermodynamic limit.
In I we have shown
that the segment with $\epsilon < 0$ in the parameter space (see right panel of
Fig.\ref{fig:pspace})
is the coexistence line of states spontaneously-magnetized in opposite
directions, when PBC are imposed,
while it is a line of critical points with APBC. Since this is a crucial point, let
us overview the equilibrium picture before turning to the discussion of the
quench dynamics.

\begin{figure}[ht]
\centering
	\includegraphics[width=10cm]{pannello2-sublabel.pdf}
        \caption{Typical equilibrium configurations below $T_c$.
        In the PBC case (a) one black domain of up spins fills
        the entire systems. Thermal fluctuations produce the small white domains of
        down spins. In the APBC case (b) there are two large domains separated
        by one interface cutting across the system. Within each domain there are
        the small patches of reversed spins due to thermal fluctuations.}
\label{fig:conf}
\end{figure}

\subsection{{\bf Equilibrium with PBC}}
\label{PBCsotto}

When PBC are imposed, two confining components of spins aligned
either prevalently up or prevalently down, are formed in phase space.
A configuration typical of the up component
is shown in the left panel of Fig.\ref{fig:conf}.
In the thermodynamic limit these components become absolutely
confining, ergodicity breaks down and, therefore, we are confronted with the same
situation discussed in the $1d$ case at $T=0$. Namely,
the Gibbs state becomes
the even mixture of the two broken-symmetry pure states like in Eq.~(\ref{mixt.0}),
that is
\begin{equation}
P_{\rm eq}^{(p)}(\boldsymbol{s}) =
\sum_\alpha p(\alpha) P_\alpha (\boldsymbol{s}).
\label{mix.1}
\end{equation}
Here, $\alpha = \pm$ is the component label, the mixing
probability is uniform $p(\alpha)=1/2$ and $P_\alpha (\boldsymbol{s})$ is the
ferromagnetic pure state. The nonvanishing spontaneous
magnetization density
$m_\alpha$ is given by
\begin{equation}
m_-=-m_+,  \quad |m_\alpha| = |\epsilon|^\beta, \quad \beta=1/8.
\label{mix.2}
\end{equation}
Using the above definitions and rewriting the Gibbs average in terms of the component
averages, i.e.  $\langle \cdot \rangle_{\rm eq} = \sum_\alpha p(\alpha)
\langle \cdot \rangle_\alpha$, the equilibrium correlation
function can be rearranged in the form
\begin{equation}
C_{\rm eq}^{(p)}(r,\epsilon) = \overline{ \langle (s_i -m_\alpha)
(s_{i+r} -m_\alpha) \rangle_\alpha} \; + \;
\overline{[\langle s_i \rangle_\alpha - \overline{m_\alpha}]
[\langle s_{i+r} \rangle_\alpha - \overline{m_\alpha}]},
\label{mix.4}
\end{equation}
where the overline denotes averaging with respect to $p(\alpha)$.
The first contribution is the average over components of the {\it intra}-component
correlation function
$\langle \psi_i \psi_{i+r} \rangle_\alpha$, where the variables $\psi_i = s_i - m_\alpha$
represent the thermal fluctuations in the pure state
$P_\alpha (\boldsymbol{s})$. As it is intuitively clear, deviations from the
average by symmetry do not depend on $\alpha$, so we shall use the notation
$G_{\rm eq}(r,\epsilon)$ for $\langle \psi_i \psi_{i+r} \rangle_\alpha$.
At low $T_F$ this quantity is short ranged,
since in the broken-symmetry state the correlation length $\xi_\psi$
of the $\psi$ variables vanishes as $T_F \to 0$.
The second term, instead, represents the {\it inter}-components
contribution, which reduces to $m^2_\alpha$, since
$\overline{m_\alpha}=0$ and $m^2_\alpha$ is independent of $\alpha$. Thus,
in the end, from the Gibbs average we have
\begin{equation}
C_{\rm eq}^{(p)}(r,\epsilon) = G_{\rm eq}(r,\epsilon) +m^2_\alpha.
\label{mix.5}
\end{equation}
It is important, for what follows, to keep in mind that the constant term
$m^2_\alpha$, which is the variance of the variable $m_\alpha$ distributed according
to $p(\alpha)$, arises exclusively from the mixing as
the constant term in Eq.~(\ref{1dI.16}).
Therefore, in the PBC case the only dynamical variables
are the $\psi_i$, which means that the dynamical rule updates $\psi_i$, but not $m_\alpha$.
The magnetization distribution exhibits
the double peak structure~\cite{Binder,Bruce}
which, in the thermodynamic limit, becomes the sum of the two $\delta$ functions
\begin{equation}
P^{(p)}(m|\epsilon) = \frac{1}{2} [ \delta(m-m_-) + \delta(m-m_+)].
\label{mix.60}
\end{equation}
Hence, as explained in the $1d$ case, the meaningful averages are those taken
with the broken symmetry ensembles $P_\alpha (\boldsymbol{s})$, which coincide with
time averages and give
\begin{equation}
C_{\alpha,\rm eq}^{(p)}(r,\epsilon) = G_{\rm eq}(r,\epsilon).
\label{mix.5bis}
\end{equation}

\subsection{{\bf Equilibrium with APBC}}
\label{APBCeq}

When APBC are imposed, like in the $1d$ case ergodicity does
not break. As explained in I, there is only one ergodic component,
whose typical configurations at sufficiently low $T_F$ are composed of two large
ordered domains, separated by one interface cutting across the system and
sweeping through it, as illustrated in the
right panel of Fig.\ref{fig:conf}. This suggests to split the
spin variable into the sum of two independent
components
\begin{equation}
s_i = m_{\alpha(i)}+ \psi_i,
\label{split.1}
\end{equation}
where $\alpha(i)=\pm$ is the label of the
domain to which the site $i$ belongs and $\psi_i=s_i - m_{\alpha(i)}$ is, as before,
the thermal
fluctuation variable. The significant difference with respect to the previous
case is that now $\psi_i$ and $m_{\alpha(i)}$ are both dynamical variables,
since the fluctuations of  the latter one are not due to the mixing of pure states,
but to the transit
of the interface through the site $i$, which means that the dynamical rule updates
both $\psi_i$ and $m_{\alpha(i)}$.
Using the independence of these variables and
the vanishing of averages $\langle s_i \rangle_{\rm eq} =
\langle m_{\alpha(i)} \rangle_{\rm eq} =\langle \psi_i \rangle_{\rm eq} =0$,
the correlation function can be written as the sum of two contributions
\begin{equation}
C_{\rm eq}^{(a)}(r,\epsilon,L^{-1}) =  G_{\rm eq}(r,\epsilon)
+ D_{\rm eq}(r,\epsilon,L^{-1}),
\label{mix.6}
\end{equation}
which have quite different properties.
The first one, which is the same as in Eq.~(\ref{mix.5}), is short ranged.
The $L$ dependence has been neglected, because we may always assume
that the conditions for $\xi_\psi \ll L$ are realized.
The second one, which contains the correlations of the
background variables $m_{\alpha(i)}$, i.e.
\begin{equation}
D_{\rm eq}(r,\epsilon,L^{-1}) =
\frac{1}{V} \sum_i  \langle m_{\alpha(i)} m_{\alpha(i+r)} \rangle_{\rm eq},
\label{mix.7}
\end{equation}
has been studied numerically in I and scales as
\begin{equation}
D_{\rm eq}(r,\epsilon,L^{-1}) = \frac{1}{r^a}Y(\epsilon,r/L), \quad \text{where}
\quad Y(\epsilon,x) = m^2_\alpha (1-x).
\label{mix.8}
\end{equation}
We have retained the power-law prefactor $r^{-a}$ in front, even though now
$a=0$ because the correlated clusters of the background variables are compact,
in order to emphasize the
similarity with Eq.~(\ref{sat.1}) and to render it evident by inspection
that the correlation length $\xi_m$ of these variables
coincides with $L$. Notice that $\epsilon$ does not enter the scaling function but
only its amplitude through $m^2_\alpha$.

From the divergence of $\xi_m$ in the thermodynamic limit, there
follows that the whole segment on the $\epsilon$ axis with $\epsilon < 0$
is a locus of critical points, as anticipated above. The corresponding critical
properties can be extracted by
using $L^{-1}$ as the parameter of approach to criticality.
It should be clear that this is bulk criticality, in no way related to
the properties of the interface, to which the attention of previous studies of the
APBC model was primarily directed.
In I we have shown that
the exponents satisfy the relations $\dot{\beta}/\dot{\nu} = 0$ and
$\dot{\gamma}/\dot{\nu}=d$, where the dots identify the exponents with respect to $L^{-1}$,
e.g. from $\xi_m \sim  L$,
follows $\dot{\nu} = 1$. This implies $\dot{\beta} = 0$ and
$\dot{\gamma} = d$. Hence, the
hyperscaling relation $2\dot{\beta} + \dot{\gamma} = \dot{\nu}d$ is satisfied,
suggesting that the upper critical dimensionality might diverge.
So, if now we take the thermodynamic limit,
from Eqs.~(\ref{mix.6}) and~(\ref{mix.8}) we get
\begin{equation}
C_{\rm eq}^{(a)}(r,\epsilon) =  G_{\rm eq}(r,\epsilon)
+ \frac{m^2_\alpha}{r^a},
\label{mix.6bis}
\end{equation}
and, consequently, the susceptibility of the background variables
$\chi^{(a)}_m$ diverges like
\begin{equation}
\chi^{(a)}_m(\epsilon,L^{-1}) \sim  L^d,
\label{mix.9bis}
\end{equation}
in agreement with Eq.~(\ref{susc.1}), the $m$-CK
clusters being compact.
The strong magnetization fluctuations,
implied by the divergence of the susceptibility, are indeed exhibited by the
distribution $P^{(a)}(m)$ which, instead of being double peaked like
in Eq.~(\ref{mix.60}), has been shown in I to be uniform over the interval
$[m_-,m_+]$. The qualitative difference between
$P^{(p)}(m)$ and $P^{(a)}(m)$ is the same previously analyzed
in the $1d$ case and schematically represented in Fig.\ref{fig:1}.
The uniformity of $P^{(a)}(m)$ is the distinctive feature
which highlights the difference between condensation of fluctuations and the
usual ordering transition associated to the double-peak structure of Eq.~(\ref{mix.60}).

\subsection{Relaxation Dynamics}

When the relaxation of the infinite system is studied, by taking first
the thermodynamic limit, the dependence on BC
is expected to disappear,
because at any finite time the correlation length is limited by $R$. This is
confirmed by the snapshots of the typical configurations (see Fig.\ref{fig:conf1})
taken after the quench to $T_F/T_c = 0.79$.
The top panel depicts the PBC case and the
bottom panel the APBC one. In each panel time increases from left to right.
The first three snapshots, taken at $t=1,10,100$, display the
self-similar morphology characteristic of coarsening domains, which does not show
to be affected by the
type of the imposed BC, because $R \ll L$. The BC influence is evident, instead, in the fourth snapshot taken
at $t=10^5$, when $R \gtrsim L$ and the system has equilibrated.
\begin{figure}[ht]
\centering
\includegraphics[width=10cm]{pannello-sublabel.pdf}

	\caption{The two panels show the sequence of snapshots taken at
	$t=1, 10, 100, 10^5$ with PBC (a) and APBC (b), after a quench
	at $T_F=1.8$, with $L=256$. Time increases from left to right.
	The first three configurations belong to the coarsening regime
	and show independence from BC. The last pair of configurations is
	morphologically similar to those in Fig.~\ref{fig:conf}  and shows that at
	$t=10^5$ the system has equilibrated.}
\label{fig:conf1}
\end{figure}

The configurations morphology, with large compact growing domains containing
in their interior small patches of thermal fluctuations,
suggests to generalize
to the off-equilibrium regime the split of variables~(\ref{split.1}) by
$s_i(t) = m_{\alpha(i,t)}+ \psi_i(t)$, where
$\alpha(i,t)$ is the label of the domain to which the site $i$ belongs at
the time $t$. Then, as in Eq.~(\ref{mix.6}), the correlation
function separates into the sum of two contributions
\begin{equation}
C(r,\epsilon,t^{-1}) =  G_{\rm eq}(r,\epsilon)
+ D(r;\epsilon,t^{-1}),
\label{gen.19}
\end{equation}
where the first one is BC-independent, time-independent
and identical to the analogous term appearing in Eqs.~(\ref{mix.5}) and~(\ref{mix.6}),
because thermal fluctuations equilibrate quickly.
The second contribution contains the correlations of the background
variables and obeys scaling in the form
\begin{equation}
D(r,\epsilon,t^{-1}) = \frac{m^2_\alpha}{r^a} F(r/R),
\label{gen.20}
\end{equation}
where $a=0$, due to the compactness of domains, the growth law $R(t) = t^{1/z}$
is the same of Eq.~(\ref{anml.3bis})
with $z=2$ and the $\epsilon$ dependence has been factorized in the amplitude $m^2_\alpha$.
Comparing with Eq.~(\ref{sat.3}), we see that the same behavior as in the quench to $T_c$
is obtained, apart for the change of the exponents $z$ and $a$, and for
the specific forms of the functions $F_c(x)$ and $F(x)$.

\begin{figure}[ht]
\centering
\includegraphics[width=10cm]{fig8.pdf}
\caption{Collapse on the master curve of Eq.~(\ref{quench.1}) of the data for
$C(r,t^{-1},L^{-1})$ in the time regime $R \ll L$, system size $L=256$.
PBC (black symbols) and  APBC (empty symbols).
The continuous line is the plot of the Ohta-Jasnow-Kawasaki function
defined in Eq.~(\ref{OJK}).}
\label{fig:scaledCbis}
\end{figure}

The above statements are substantiated by the plot in Fig.\ref{fig:scaledCbis}
of the numerical data
for the equal-time correlation function, generated for a quench to
$T_F/T_c=0.666$, which corresponds to $\epsilon = -0.9$, with $L=256$ and with
both PBC and APBC. The APBC data have been
circularly averaged to smooth out the anisotropy induced by the cylindrical BC.
The good collapse of the data, in the time
regime such that $R \ll L$, shows that for the chosen value of $T_F$
the thermal fluctuations contribution is negligible.
Moreover, the master curve $F(x)$ compares well with the
Otha-Jasnow-Kawasaki~\cite{OJK} approximate result
\begin{equation}
F(x) = \left (\frac{2}{\pi} \right ) \arcsin (\gamma), \quad \gamma = \exp(-x^2/b),
\label{OJK}
\end{equation}
where $b$ is a constant, as it is demonstrated by Fig.\ref{fig:scaledCbis}.
Therefore, the relaxation to below $T_c$ is not qualitatively different from the one
to $T_c$. Both are coarsening processes and both do not depend on the
imposed BC. Differences between the two are in the quantitative details,
like the values of the $z$ exponent, the dimensionality of correlated clusters
and the shape of the scaling functions.
The implication is
that also in the quench to below $T_c$ the system tends toward a
critical state, because the time-dependent correlation length $R$ diverges,
eventually yielding the time-asymptotic critical correlation function
\begin{equation}
C^*(r,\epsilon)= G_{\rm eq}(r,\epsilon) + \frac{m^2_\alpha}{r^a}.
\label{dyna.5}
\end{equation}
It is then evident, according to the discussion made at the end
of subsection~\ref{APBCeq}, that this asymptotic form, which
we emphasize once more is the same for both choices of BC,
matches $C^{(a)}_{\rm eq}(r;\epsilon)$ but not $C^{(p)}_{\rm eq}(r;\epsilon)$.
Finally, recalling that $a=0$, it is straightforward to see from Eq.~(\ref{gen.20})
that the background susceptibility scales like
\begin{equation}
\chi (\epsilon,t^{-1}) \sim  R^d,
\label{dyna.6}
\end{equation}
in agreement with the result~(\ref{mix.9bis}) for $\chi^{(a)}_m(\epsilon,L^{-1})$.

In conclusion, in the PBC case, as anticipated in the Introduction,
the asymptotic state and the equilibrium one are remote one from the other, and the system
may be
regarded as remaining strongly out of equilibrium, because in the former one there are
long-range correlations, which are absent in the second one.
In the APBC case, instead, both the asymptotic and the equilibrium state are
critical and with the same universal properties, hence the system equilibrates
although with an infinite equilibration time, just as in the quench to $T_c$.

\subsection{Summary}

So far we have shown that when the Ising model is quenched in the
two-phase region, i.e. to below $T_c$ for $d=2$ and to $T_F=0$ for $d=1$, the APBC system equilibrates and the PBC one remains off equilibrium.
The basic elements of the mechanism underlying this phenomenology
are as follows:

\begin{enumerate}

\item To different BC, in principle, there correspond different statistical ensembles.

\item These ensembles become {\it equivalent} when the limits are taken according
to the sequence: $L^{-1} \to 0$ first and then $t^{-1} \to 0$, for all temperatures $T_F$.

\item Instead, the ensembles {\it may} become {\it non equivalent}, depending on BC,
when the limits are taken in the reverse sequence: $t^{-1} \to 0$ first
and then $L^{-1} \to 0$ with $T_F < T_c$.

\item
Equivalence fails with PBC because of ergodicity breaking, and holds
with APBC since ergodicity is preserved.

\item Ergodicity breaking induces spontaneous symmetry breaking, which
makes correlations short-ranged.

\item Instead, when ergodicity holds an unusual type of
criticality sets in, with long-range correlations and  {\it compact} correlated domains.

\end{enumerate}

In the next section we shall show that this is not just a peculiarity of
the Ising model, but that it is a more general phenomenon, since it
takes place with the same characteristics also in the quench to below the critical point
of the spherical model, without invoking the imposition of different types of BC.
In fact, two different ensembles arise not from the choice of BC,
which is taken to be the standard PBC
one, but from enforcing the spherical constraint
either sharply or smoothly. These ensembles
turn out to be equivalent or non equivalent, just as in the Ising case, depending on
the order of the $L^{-1} \to 0$ and $t^{-1} \to 0$ limits.

\section{Spherical models}
\label{SM}

\subsection{Equilibrium}

Let us briefly recall what the spherical model is about starting from
equilibrium, which means that the $t^{-1} \to 0$ limit has been taken beforehand.
Consider a classical paramagnet in the volume $V=L^d$ and
with the energy function~\cite{Ma}
\begin{equation}
{\cal H}(\boldsymbol{\varphi}) = \int_V d\vec r \, \varphi(\vec r) \left
(-\frac{1}{2}\nabla^2 \right )\varphi(\vec r),
\label{gauss.1}
\end{equation}
where $\boldsymbol{\varphi}$ stands for a configuration  of the local, continuous
and unbounded spin variable $\varphi(\vec r)$. PBC are understood throughout.
Due to its bilinear character, the above Hamiltonian can be diagonalized by Fourier transform
\begin{equation}
{\cal H} = \frac{1}{2V}\sum_{\vec k} k^2|\varphi_{\vec k}|^2.
 \label{gauss.1bis}
\end{equation}
In the spherical model (SM) of Berlin and Kac~\cite{BK} a coupling among the modes is
induced by the imposition of an overall sharp constraint on
the square magnetization density
\begin{equation}
\mathit{s}(\boldsymbol{\varphi}) = \frac{1}{V} \int_V d \vec r \,
\varphi^2(\vec r) = \frac{1}{V^2}\sum_{\vec k} |\varphi_{\vec k}|^2 = 1.
\label{gauss.2}
\end{equation}
Then, in thermal equilibrium the statistical ensemble is given by
\begin{equation}
P_\textrm{SM}(\boldsymbol{\varphi}) = \frac{1}{Z_\textrm{SM}}
e^{-\beta {\cal H}(\boldsymbol{\varphi})}
\, \delta \left (\mathit{s}(\boldsymbol{\varphi})-1 \right ),
\label{Gauss.4}
\end{equation}
where $Z_\textrm{SM}$ is the partition function.
A variant of the model, called
the mean-spherical model (MSM)~\cite{LW,KT}, is obtained by imposing the constraint in the mean: An exponential bias is introduced in place of the $\delta$ function
\begin{equation}
P_\textrm{MSM}(\boldsymbol{\varphi}) = \frac{1}{Z_\textrm{MSM}}e^{-\beta [{\cal H}(\boldsymbol{\varphi}) +\frac{\kappa}{2} {\cal S}(\boldsymbol{\varphi})]},
\label{Gauss.3}
\end{equation}
where ${\cal S}(\boldsymbol{\varphi}) = V\mathit{s}(\boldsymbol{\varphi})$ and the parameter $\kappa$ must be so adjusted to satisfy the requirement
\begin{equation}
\langle \mathit{s}(\boldsymbol{\varphi}) \rangle_\textrm{MSM} = 1.
\label{msph.1}
\end{equation}
Although it is the common usage to refer to these  as models, it should be clear from Eqs. (\ref{Gauss.4}) and (\ref{Gauss.3}) that we are dealing with two conjugate ensembles, distinguished by conserving or letting to fluctuate the density $\mathit{s}$.

In both models there exists a phase transition at the same critical
temperature $T_c$, above which they are equivalent and below which they are not,
which means that the nature of the low temperature phase is different.
It is worth, here, to go in some detail~\cite{CSZ} because the point is quite illuminating
on the equivalence or lack-of issue.
Let us separate in $\mathit{s}$ the excitations from the ground-state
contribution
\begin{equation}
\mathit{s} = \mathit{s}_0 +  \mathit{s}^*, \quad \text{with} \quad
\mathit{s}_0 = \frac{1}{V^2}\varphi_0^2,  \quad
\quad \mathit{s}^*= \frac{1}{V^2} \sum_{\vec k \neq 0} |\varphi_{\vec k}|^2.
\label{s.1}
\end{equation}
Then, taking the average in either ensemble, from the spherical
constraint follows the sum rule
\begin{equation}
\langle  \mathit{s}_0 \rangle +  \langle \mathit{s}^* \rangle = 1,
\label{s.2}
\end{equation}
which must be satisfied at all temperatures and it is the motor of the transition.
In fact, in the thermodynamic limit the excitations contribution is superiorly
bounded~\cite{CSZ} by
\begin{equation}
\langle \mathit{s}^* \rangle \leq TB,
\label{s.3}
\end{equation}
where $B$ is a dimensionality-dependent positive constant, which is finite
for $d > 2$ and diverges at $d=2$. Therefore, by enforcing the
constraint~(\ref{s.2}) there remains defined
the critical temperature
\begin{equation}
T_c = 1/B,
\label{s.4}
\end{equation}
above which the sum rule~(\ref{s.2}) is saturated without any contribution from
$\langle  \mathit{s}_0 \rangle$, while below there must necessarily be a finite contribution
from the ground state, yielding
\begin{equation}
\langle  \mathit{s}_0 \rangle  = \left \{ \begin{array}{ll}
         0, \;\; $for$, \;\; T \geq T_c,\\
         1 - T/T_c,\;\; $for$, \;\; T < T_c.
        \end{array}
        \right .
        \label{s.5}
        \end{equation}
Rewriting $\mathit{s}_0 = \psi_0^2$, where $\psi_0$ is the density $\frac{1}{V}\varphi_0$,
the question is how can there arise a finite contribution to $\langle  \mathit{s}_0 \rangle$
from this single degree of freedom and here is precisely where
the two models differ. In the SM the sharp
version~(\ref{gauss.2}) of the constraint  introduces enough nonlinearity for the
transition to take place by {\it ordering}.
This means that ergodicity breaks down inducing the spontaneous breaking of
the $\mathbb{Z}_2$ symmetry. Then, exactly like in the Ising model with PBC,
the probability
distribution of the magnetization density, that is of $\psi_0$, results from the mixture
of the two pure ferromagnetic states
\begin{equation}
P_\textrm{SM}(\psi_0) =
\frac{1}{2}[\delta (\psi_0 - m_-) +
         \delta (\psi_0 - m_+)],
\label{ensmbl.010}
        \end{equation}
where $m_\pm = \pm \sqrt{1 - T/T_c}$ is the spontaneous magnetization.
Thus, in this case $\langle  \mathit{s}_0 \rangle_{SM}$ stands for
the square of the spontaneous magnetization $m_\pm^2$.
     Instead, in the MSM ordering cannot take place, because the soft
     version~(\ref{msph.1}) of the constraint leaves the statistics Gaussian.
     Neither ergodicity nor symmetry break down, as in the Ising APBC case. Then, below $T_c$,
     the only mean to build up the finite value of $\langle  \mathit{s}_0 \rangle_{MSM}$
     needed to saturate the sum rule is by growing the fluctuations
     of $\psi_0$ through the spread out probability distribution given
     by
\begin{equation}
P_\textrm{MSM}(\psi_0) =
         \frac{e^{-\frac{\psi_0^2}{2(1-T/T_c)}}}{\sqrt{2\pi (1 - T/T_c)}}.
\label{ensmbl.001}
\end{equation}
Therefore, now $\langle  \mathit{s}_0 \rangle_{MSM}$
     stands for the macroscopic
     variance of $\psi_0$. Elsewhere~\cite{EPL,CCZ,Zannetti,Merhav,Marsili}, this type of transition, characterized by the fluctuations of an extensive quantity condensing into one microscopic component, has been referred to as {\it condensation of fluctuations}.

\begin{figure}[!tb]
\centering
\includegraphics[width=5cm]{fig9a-sublabel.pdf} \hspace{0cm}
\includegraphics[width=5cm]{fig9b-sublabel.pdf}
	\caption{Magnetization distribution in the MSM model (a) and in the SM model (b),
  for $T < T_c$. The spikes in the right panel stand for $\delta$ functions.}
\label{fig1}
\end{figure}
Comparing Figs.~\ref{fig1} and~\ref{fig:1}, it is evident that the distributions are the same in the two cases
where ergodicity breaks down, that is in the Ising model with PBC and in the SM.
In the other two cases, Ising with APBC and MSM, the distributions are not
superimposable but show the same physical phenomenon: ergodicity is preserved
by developing macroscopic fluctuations of the magnetization, which remain finite
in the thermodynamic limit and reveal
the critical nature of the low temperature phase. In fact, in the MSM the structure
factor, i.e. the Fourier transform of the correlation function, is given by~\cite{CCZ}
\begin{equation}
C_\textrm{MSM}(\vec k) = \frac{T}{k^2} + m_\pm^2 \delta(\vec k).
\label{str.1}
\end{equation}
The two terms appearing above are the anolouges
in Fourier space of those entering $C_{\rm eq}^{(a)}(r;\epsilon)$
in Eq.~(\ref{mix.6bis}), with the correspondences
\begin{equation}
G_{\rm eq}(r;\epsilon) \longleftrightarrow \frac{T}{k^2}, \quad
\frac{m^2_\alpha}{r^a} \longleftrightarrow m_\pm^2 \delta(\vec k).
\label{str.2}
\end{equation}
Notice that, as it is well known, the thermal fluctuations contribution in the MSM
is massless, i.e. is critical, at all temperatures below $T_c$. For simplicity, let
us set $T=0$ in order to get rid on this contribution and to focus on
the interesting one, which is the $\delta$-function term (Bragg peak).
We emphasize that this
is the Fourier transform of the background {\it critical} contribution with
compact correlated clusters, just as the corresponding term in the Ising APBC case.

Finally, we point out that the $d=2$ case is analogous to Ising with $d=1$,
because $T_c$ vanishes.
However, for brevity, we shall not elaborate on this case here.

\subsection{Dynamics}

Let us next consider the relaxation dynamics in the quench to $T_F=0$. In Ref.~\cite{Fusco} it was shown that,
when the thermodynamic limit is taken first, the two models are equivalent
at all times. Then, it is an exact result that the dynamical structure factor
both for the SM and MSM is given by
\begin{equation}
C(\vec k,t) = \Delta \left (1 +\frac{2r^2_0}{R^2} \right )^{d/2} R^d \, e^{-(kR)^2},
\label{str.3}
\end{equation}
where $C(\vec k,0)=\Delta$ is the spatially uncorrelated initial condition at $T_I=\infty$,
$R = \sqrt{2t}$ is the growth law for nonconserved
dynamics and $r_0 = 1/(\sqrt{2}\Lambda)$ is the microscopic length related
to the momentum cutoff $\Lambda$, which is imposed exponentially
when integrating over $\vec k$ and is responsible for the corrections to
scaling in the early regime.
From the normalization condition at $t=0$
\begin{equation}
\int \frac{d^d k}{(2\pi)^d} \, C(\vec k,0) e^{-k^2/\Lambda^2} = 1,
\label{str.4}
\end{equation}
there follows $\Delta = (2\sqrt{\pi})^d$. Inserting this into Eq.~(\ref{str.3}), with
little algebra one can verify that indeed the normalization is satisfied at all
times. Then, since the peak grows like $C(0,t) \sim R^d$, one can conclude
that the asymptotic structure factor is the $\delta$ function
\begin{equation}
\lim_{t \to \infty} C(\vec k,t) = C^*(\vec k) = \delta(\vec k),
\label{str.5}
\end{equation}
which matches the equilibrium Bragg peak~(\ref{str.1}) in the MSM. Since
the growth of $R$ implies that the asymptotic state is critical and with compact
correlated clusters, in the quench to below $T_c$ the MSM approaches
equilibrium arbitrarily close, while the SM, which is not critical in the equilibrium
state, remains permanently out of equilibrium.

In conclusion, going through all the items listed at the end of the previous
section, one can check that perfect correspondence between
Ising-PBC and SM on one side and Ising-APBC and MSM on the other
is established.

\section{Concluding Remarks}
\label{CR}

In this paper we have addressed a problem which is of basic interest in the
physics of slowly relaxing systems. Since slow relaxation means that
equilibrium is not reached in the observable time scale, relevant
questions are whether a criterion for equilibration, or for lack of, can be established
and, if so, whether the nature of the equilibrium state
can be inferred from the available dynamical information. Although the task of giving
general answers to these questions is of formidable difficulty,
we have shown that, at least in the restricted realm of phase-ordering systems,
it is possible to arrive at some definite conclusions.

By analyzing the relaxation of the Ising model after temperature
quenches, we have found that
the system does or does not equilibrate, depending on whether the dynamics at
the final temperature of the quench is ergodic or not.
This has been established by investigating the dependence
of the spin-spin correlation function upon the order of the large-time and
thermodynamic limits, when different BC are imposed. The findings are that the
APBC system equilibrates in all conditions, because the dynamics are ergodic at
all temperatures, while the PBC system does not equilibrate for $T_F < T_c$,
because that's where ergodicity does not hold.
These statements are strengthened and corroborated by
exact analytical results from
the quench of the spherical and mean spherical model, which reproduce
very closely, although in a quite different context, the picture just outlined.
We may then answer
the first question asked at the beginning of the section by saying that
it might take an infinite time to equilibrate, but nonetheless
the system can get arbitrarily close to equilibrium if the dynamics are ergodic.
Instead, if ergodicity is broken, and if the initial state is symmetric,
the system does not get close to equilibrium, no matter how long is let to
relax.

For what concerns the second question, the answer is that, yes, once it is
established that the system approaches equilibrium, then the
nature of the equilibrium state can be inferred from the dynamical information.
Consider first the quench to $T_c$, in which case the time-dependent
correlation function obeys the scaling form~(\ref{sat.3}), while in the equilibrium
state it decays according to the pure power law~(\ref{sat.2}). It is then evident that
the latter result can be reconstructed from the short distance behaviour,
i.e. for $r \ll R$, at times finite but large enough to detect a clean scaling behaviour.
The same procedure applies also in the case of the quench to below $T_c$ with
APBC, where the background component of the correlation function is given by
Eq.~(\ref{gen.20}). Then, again the short distance approximation, which in this case
is a constant term since $a=0$, reproduces correctly the form of the equilibrium
critical correlation function.

Finally, let us comment on the nature of the line of critical points on the
$\epsilon < 0$ segment in the Ising APBC case.
According to the view put forward in this paper,
$t^{-1}$ is just another relevant parameter measuring
the distance from criticality, on the same footing with $\epsilon$ and $L^{-1}$, so that
these critical points control both statics and dynamics. In subsection \ref{APBCeq}
we have pointed out
that the static critical exponents, defined with respect to $L^{-1}$,
satisfy the hyperscaling relation
$2\dot{\beta} + \dot{\gamma} = \dot{\nu }d$ for all $d$, suggesting that
the upper critical dimension is at $d=\infty$. It is then interesting to note
the concomitance with the fact that the
Otha-Jasnow-Kawasaki approximate theory, which accounts well for the time
dependent correlation function as shown in Fig.\ref{fig:scaledCbis},
becomes exact in the $d \to \infty$ limit~\cite{Bray}.

\begin{acknowledgements}
A.F. acknowledges financial support of the MIUR PRIN 2017WZFTZP "Stochastic forecasting in complex systems".
\end{acknowledgements}