\section{Introduction}
\label{sec:introduction}

Understanding the factors and reasons determining the effectiveness and acceptance of given pieces of artistic or scientific work represents a continuing challenge in artificial intelligence (e.g.,~\cite{wang2013quantifying,barabasi2018formula,lee2018predicting,tohalino2022predicting,cetinic2019deep}). As is often the case with complex systems, not only a large number of possible factors is potentially involved, but their individual and combined effects also tend to be highly non-linear. In this manner, small effects can lead to considerable impacts, being also likely to vary along time and space in modes that are hard to predict.

Among the several aspects that are more likely to influence the visibility and accomplishment of an artistic piece, we have its intrinsic \emph{quality}, \emph{innovation}, and \emph{affinity} with the main trends, interests, and expectations predominating in a given period and place. All these three main aspects are not only challenging to \emph{define}, but even more so to \emph{predict}, which has motivated growing interest from the scientific community (e.g.,~\cite{wang2019success}).

A better understanding of the motivations why an artistic piece becomes successful constitutes a particularly interesting objective for a handful of reasons: (i) this type of study can motivate the development of new concepts and methods capable of quantifying the three main aspects identified above, namely quality, innovation, and affinity of an artistic piece; (ii) that kind of research has great potential for revealing important aspects of the mechanisms underlying human preferences for specific subjects and styles along time and space; (iii) such developments can lead to strategies for predicting the acceptance of certain types of works, which may provide subsidies and motivation for developing new and more effective artistic pieces.

The present work aims at studying whether it is feasible to characterize and identify stories and narratives listed as best sellers by combining full-text content information and machine learning models. In this regard, the textual content of a set of books was modeled, and a series of experiments assessed the possibility of automatically differentiating a best seller from an ordinary book. In particular, we employed a dataset encompassing the full-text content of literary works collected from the Project Gutenberg platform. The dataset was split into two categories: \emph{success} (books that appear at least once in the \emph{Publishers Weekly Bestseller Lists}) and \emph{others}. After applying a preprocessing step (removal of stopwords, lemmatization, and tokenization), the content of each book was embodied in terms of a word embedding representation by using the bag-of-words~\cite{manning1999foundations} and doc2vec~\cite{le2014distributed} approaches. Finally, we employed different strategies to assess the prediction of the success of books in terms of their embedding representations, including: (i) visualization approaches, namely the linear discriminant analysis (LDA)~\cite{ghojogh2019linear} and SemAxis~\cite{an2018semaxis} techniques; (ii) classification approaches, encompassing different models and cross-validation strategies.

{In contrast to previous studies, here we rely on one of the prime published sources of best sellers book lists, namely the \emph{Publishers Weekly Bestsellers Lists}, which comprises the best-selling books every year since 1885. Although its criteria to define a book as an absolute success is not entirely specified, it is established that every considered paperbound book sold at least 2,000,000 copies, and every selected hardbound book sold 750,000 copies or more. It is also settled that Publishers Weekly only regards books distributed through the trade -- that is, bookstores and libraries --, not including those sold by mail or book clubs~\cite{hackett1977}. Besides that, our work compounds the list of few studies which analyzed the success factor by analyzing the full-text content of the texts, posthumously modeling it through embeddings, and analyzing it both qualitatively (applying visualization and seeking for words that lead to discrimination) and quantitatively (involving supervised classifiers).}

{The obtained results suggest that it is infeasible to predict the success of a literary work with high accuracy by using only its full-text content. The best classification accuracy acquired throughput the value of 0.75, combining a bag-of-words representation with a logistic regression model, which is a fair-to-middling outcome. Nonetheless, our experiments evince that the subject of the books does not seem to be a core factor for a title becoming a best seller and that there are words more typically found in this category of books.}

This work is organized as follows. Section~\ref{sec:related_works} presents and discusses the related works. In Section~\ref{sec:res_questions}, we present the research questions. Section~\ref{sec:dataset} describes the used datasets. Section~\ref{sec:methodology} describes the methodology adopted to analyze the books, including text preprocessing, representation, visualization and classification. The results and discussions are reported in Section~\ref{sec:results}. Finally, in Section~\ref{sec:conclusions}, we present the conclusions and future works.

\section{Related Works}
\label{sec:related_works}

The study conducted in~\cite{yucesoy2018success} analyzed the success of books using as reference the \emph{The New York Times Best Sellers}, which includes a list of best-selling books in the United States. The authors considered the books appearing on the list between August 2008 and March 2016. As additional information, the sales patterns of books were also considered by using data from \emph{NPD BookScan}~\cite{yucesoy2018success}.
Several interesting results were reported. Fiction books were found to be more likely to become best sellers, while nonfiction books tended to be sold with lower intensity. The authors also proposed a model that can accurately measure long-term impact since it can predict the number of copies sold by best sellers short after their release. The proposed description was found to be consistent with a previous model devised to describe the attention received by scientific papers~\cite{wang2013quantifying}. The authors argue, therefore, that the underlying processes of attention are similar -- despite the differences in time scale.

A model to predict book sales was proposed in~\cite{wang2019success}. The authors used as a dataset the \emph{NPD Bookscan}, focusing on a list of the 10 thousand top-selling books in a given period. A machine learning approach was proposed using different book features. Authors' visibility was taken into account by measuring the public interest in authors via Wikipedia page views. Previous sales were also considered as a feature to measure the previous success of authors. Book features included genre (e.g., horror and science fiction) and topic information (as provided by readers). In addition, publishers' information was used. All features were combined in the so-called Learning to Place (L2P) machine learning algorithm~\cite{wang2019l2p}, which aims at classifying a new instance (i.e., predicting book sales) within a sequence of previously published books. This study found that in fiction and nonfiction books, the publisher quality tends to play an important role in the prediction. The visibility of authors was also found to be an important feature, as more visible authors potentially are more likely to sell more copies. Finally, the other factors related to the text content itself (e.g., genre and topic information) were found to play relatively a minor role in the prediction model.

Differently from previous works that did not take into account the textual content~\cite{yucesoy2018success,wang2019success}, the relevance of writing style was analyzed in~\cite{ashok2013success}. The authors analyzed full books from different genres (e.g. adventure, mystery, fiction). The dataset was collected from the \emph{Project Gutenberg} repository. Several linguist marks of writing style were used to characterize the texts. Examples include lexical features, distribution of grammar rules, and sentiment analysis. The authors used SVM as classifier~\cite{fan2008liblinear}, and download counts were used as a surrogate for the visibility of books. Additional information such as award recipients and the number of copies sold was also used to quantify success. The authors concluded that the used stylistic metrics are effective to quantify the success of novels.

{Because only a few works have analyzed the content of books to predict if they will become best sellers, in the current study we focus our analysis on full-textual features to discriminate between best sellers and ordinary literary works.}

\section{Research questions}
\label{sec:res_questions}

This study aims to probe whether the full-text content of the book alone can indicate if it will become a best seller. While there are several ways to represent a text, we focused on the most common approaches devoted to representing \emph{long texts}. For this reason, here we also investigate which text representation better grasps the information about a book becoming a success. Finally, to recognize patterns in the textual data, we also examined which classifier is the most appropriate for discriminating between successes and ordinary books.

Briefly, the main research questions here are:

\begin{enumerate}
    \item Is it possible to predict the inclusion of books into best sellers lists by analyzing only their full-text content?
    \item Can one use bag-of-words and neural network embeddings to detect informative attributes for identifying best sellers?
    \item Can the abovementioned embeddings be influenced by the subject headings available in the dataset (such as genre or literary class)?
    \item How different is the performance of supervised classifiers in discriminating between the two categories of books analyzed?
\end{enumerate}

\section{Dataset}
\label{sec:dataset}

As the main objective of this work is to understand whether it is possible to identify and characterize styles and stories classified as best sellers, our dataset was composed of two categories: \emph{success} and \emph{others}. In the first, we included all books considered best sellers; in the latter, literary works not listed as such (at least not in the analyzed period and the consulted list).

Two best seller book lists were considered in the present work: \emph{The New York Times Best Sellers}, first published in 1931, and \emph{Publishers Weekly Bestseller Lists}, published for the first time in 1895.  In the case of the first list, in the period from 1931 to the present day, only 18 titles were available on the \emph{Project Gutenberg}~\footnote{Project Gutenberg. \url{https://www.gutenberg.org}, Accessed: August 23, 2020} platform (a digital library whose collection is composed of full texts of books in the public domain). In the case of the second list, we mapped 110 titles — published from 1895 to 1924 — which therefore became part of our dataset, composing the \emph{success} category.

To select the titles of the \emph{other} category, we considered the collection of books (a) published in the same period as the selected successful ones and (b) not included in the best sellers lists of \emph{Publishers Weekly}. In this sense, if the \emph{success} class had ten titles published in 1923, the \emph{other} would have the same number of titles published in the same year — the titles randomly selected from the \emph{Gutenberg} repository. At the end of the process, this category contained 109 titles (one less than the other category, as it was infeasible to collect the same amount of titles for all the years considered).

In \cite{hackett1977} can be found the best-selling lists used in this study. It is worth mentioning that we considered only one book from each author to avoid the identification of authorship by the machine learning algorithms to be applied subsequently.

Once the dataset was ready, we cleaned up the textual content of the 219 texts to maintain only the relevant contents of the books. In this process, the header and footer included by \emph{Project Gutenberg} were removed, as well as editor/translator/author notes, captions and illustration indications, glossaries, footnotes, side-notes, annexes, and appendices.

\section{Methodology}
\label{sec:methodology}
\subsection{Text preprocessing}
Using the dataset as explained in the previous section, the preprocessing of our analysis started. First, we replaced all capital letters with their corresponding lowercase counterparts. Then, the stopwords (i.e., words that provide low or no additional meaning to the context, such as articles and connectives) were removed. Next, we performed the tokenization of the books, in which elements, like punctuations and numbers, were disregarded. Finally, the obtained words were lemmatized — being lemmatization a technique whose objective is to reduce a vocable to its canonical form and to group different forms of the same word (e.g., the term ``boys'' is reduced to ``boy'' and ``took'' becomes ``take''). Table~\ref{table:preproc} shows an example of this preprocessing.

\begin{table}
    \caption{Preprocessing of the excerpt ``\textit{It is difficult to live up to this kind of thing, and my thoughts drift to the auld schule-house and Domsie.}'', obtained from the book \emph{Beside the Bonnie Brier Bush}, by \emph{Ian Maclaren}. In the column titled \emph{Initial} is the original excerpt; in the next, the phrase without capital letters and stopwords; finally, in the last, the extract after tokenization and lemmatization processes.}
    \label{table:preproc}
    \centering
    \begin{tabular}{ P{4cm} P{5cm} P{4cm} }
    \toprule[1.0pt]
    Initial & Removed capital letters \par and stopwords & Tokenized and \par lemmatized \\
    \hline
    \textit{It is difficult} & \textit{difficult} & [`\textit{difficult}'] \\
    \textit{to live up to} & \textit{live} & [`\textit{live}'] \\
    \textit{this kind of} & \textit{kind} & [`\textit{kind}'] \\
    \textit{thing, and my} & \textit{thing,} & [`\textit{thing}'] \\
    \textit{thoughts drift} & \textit{thoughts drift} & [`\textit{thought}', `\textit{drift}'] \\
    \textit{to the auld} & \textit{auld} & [`\textit{auld}'] \\
    \textit{schule-house} & \textit{schule-house} & [`\textit{schule}', \textit{house}'] \\
    \textit{and Domsie.} & \textit{domsie.} & [`\textit{domsie}'] \\
    \bottomrule[1.0pt]
    \end{tabular}
    \end{table}

\subsection{Text embeddings}
Techniques to embed textual content have been extensively used for a variety of tasks, including grasping text similarity, sentiment analysis, and classification. Among the most widely used techniques is the bag-of-words~\cite{manning1999foundations} approach, in which the relative frequencies of words appearing in a document are organized as a vector.

Recently, other approaches, now based on neural networks, have been developed to obtain dense embedding representations of words, sentences, or entire documents, being those approaches trained to predict masked parts in texts. In this sense, among the most used techniques is word2vec, which is based on a network comprising one hidden layer and a softmax output layer. The output layer is trained for predicting the context (words appearing together) given a focus word in a sentence~\cite{mikolov2013efficient}. For a given set of sentences, such a process provides an embedding for each word.

More sophisticated techniques such as BERT~\cite{devlin2018bert} and sentence BERT~\cite{reimers2019sentence} generate embeddings that capture richer context and semantic information of words or sentences. However, these techniques are limited to a small number of tokens and can not be applied to large portions of texts, such as entire books. For this reason, we opted to use the doc2vec (D2V) method to extract a vector representation of each book~\cite{le2014distributed}.

The doc2vec approach is based on the traditional word2vec~\cite{mikolov2013efficient} pipeline with the addition of the document tags as input. More specifically, it constitutes a neural network of three layers (\emph{input}, \emph{hidden}, and \emph{softmax}), as illustrated in Figure~\ref{fig:doc2vecsemaxis}a. Just like in word2vec with a continuous bag-of-words (CBOW) architecture, the inputs are one-hot vectors representing a sequence of words from a sentence in a book. A target word is omitted from the input and used to train the neural network. In addition, the input includes an extra one-hot vector identifying the book. The model is trained to predict the target word from the context (words adjacent to the target) using a negative sampling strategy. The vectors in the hidden layer connected directly to the books encoded as one-hot are used as the book embedding. Here, we opted to use the Gensim~\cite{rehurek_lrec} software to obtain the doc2vec representations of books.

\begin{figure}[!htbp]
    \includegraphics[width=\linewidth]{scheme_doc_SemAxis.pdf}
    \caption{Representation of doc2vec and SemAxis approaches. In (a) we illustrate the neural network employed to obtain the embedding representation of books based on sequences of words (encoded as one-hot vectors) extracted from a book. The network is trained to predict a target word in the sentence based on the adjacent terms. Additionally, the original book ID is also encoded as input to the neural network, and their respective trained vectors correspond to the embedding space of books. In (b), we illustrate the SemAxis approach in which the line connecting the two categories' (\emph{success} vs. \emph{others}) centroids defines an axis to project all the books. This process results in a continuous one-dimensional (scalar) representation of books, which is employed for visualization purposes.}
    \label{fig:doc2vecsemaxis}
\end{figure}

\subsection{Visualization}

Neural network embeddings usually result in high-dimensional dense vectors that are not correlated among themselves, which limits the use of linear techniques to reduce the dimensionality of these spaces (such as PCA~\cite{gewers2021principal}). Thus, the process of visualizing such structures is usually undertaken using non-linear projections, such as t-SNE~\cite{van2008visualizing} and UMAP~\cite{mcinnes2018umap}.

However, embeddings can encode many different aspects of the data, for instance, a certain axis in a book embedding may be related to its number of pages or its adherence to the non-fiction or fantasy genres. The SemAxis approach~\cite{an2018semaxis} is a way to find an axis in a high-dimensional embedding that describes a certain aspect of the data. This is accomplished by first obtaining the centroids of two classes, e.g., small vs. larger books or non-fiction vs. fantasy books. The line connecting the two centers define an axis in which all the remaining books are projected. This process is illustrated in Figure~\ref{fig:doc2vecsemaxis}. Since in the current work we are interested in encoding the success of books, we employed the SemAxis approach to finding an axis for samples of the \emph{success} and \emph{other} classes. Similarly, in addition to SemAxis, we also employed linear discriminant analysis (LDA)~\cite{ghojogh2019linear}, which also results in an axis encoding a continuous representation of the two classes.

In contrast to neural network-based approaches, the bag-of-words embedding can result in highly correlated and sparse vectors. For instance, the frequency patterns of two close-related words can correlate strongly and rare words may only be present in a small set of documents. Nonetheless, both LDA and SemAxis are still applicable in these conditions.

\subsection{Classification}

The identification and classification of textual patterns were performed using traditional well-known machine learning classifiers~\cite{mitchell1997machine}. We considered different classifier strategies, including $k$-nearest neighbors (KNN)~\cite{fix1951knn}, logistic regression (LR)~\cite{yu2011dual}, naive Bayes (NB)~\cite{zhang2004optimality}, decision tree (DT)~\cite{breiman2017classification}, random forest (RF)~\cite{breiman2001random}, and support-vector machine (SVM)~\cite{mitchell1997machine,murthy1998automatic}. All these tests were implemented in Python language~\citep{rossum2009python} using the classifiers of Scikit-Learn~\citep{scikit-learn} library.

Following the guidelines described in related works~\cite{amancio2014systematic,rodriguez2019clustering}, we used the default parameters of the methods to classify texts. As an exception, in the case of the SVM, we changed the parameter ``max\_iter'', the maximum number of iterations, to 10,000.

\section{Results and discussions}
\label{sec:results}

This section describes the experiments performed to study the task of automatically characterizing and identifying best seller books. The proposed data analysis pipeline is illustrated in Figure~\ref{fig:schematic}. First, we obtain word embedding representations of each book by employing two distinct techniques (Figure~\ref{fig:schematic}a): bag-of-words and doc2vec (the latter with different dimensions, namely 32, 64, 128, and 256). Next, we investigate the proposed classification problem through two main approaches: visualization and classification. In the first, we employed a simple visualization pipeline to verify and illustrate the potential of using embeddings to identify best seller books (Figure~\ref{fig:schematic}b-d). The objective of this approach is to provide a preliminary and simple way to visually inspect the considered high dimensional embeddings by summarizing them into a single continuous axis.

The visualization pipeline starts with the standardization of the obtained embeddings (Figure~\ref{fig:schematic}b). To reduce the dimensionality of the embeddings, we employed SemAxis~\cite{an2018semaxis} and LDA~\cite{ghojogh2019linear}. Since these methods are supervised, the final visualizations are performed by employing the leave-one-out technique to avoid overfitting.

The second approach considered in this work is the direct application of classification methods, allowing quantitative comparison of the respective performance. For that, we employed a pipeline comprising the same embedding configurations as before, but followed by three successive stages: preprocessing, learning method, and validation, each presented as a box in Figure~\ref{fig:schematic}.  All combinations between the components of each of these boxes are considered in our evaluation.

In this sense, the following two first subsections are intended to detail the task of visualization, followed by the classification, both using the bag-of-words and then the doc2vec representation. Then, in the last subsection, we repeat these experiments to evaluate a specific variation of the constructed dataset.

\begin{figure}[!htbp]
    \includegraphics[width=\linewidth]{scheme_Results.pdf}
    \caption{Overall diagram of the main approaches, methods, representations, and validations employed in the present work. All the methods within the blue and orange boxes are applied to the two considered embeddings in a combined fashion. For example, a valid path would be: (i) embedding: doc2vec with dimension equals to 32; (ii) preprocessing: standardized; (iii) learning method: logistic regression; (iv) validation: leave-one-out.}
    \label{fig:schematic}
\end{figure}

\subsection{Bag-of-words analysis}

The first performed experiment intends to evaluate whether the frequency of words composing the books can discriminate between best sellers and ordinary literary works. For this purpose, we considered the set $S$, built based on the $3,585$ different words that appeared at least in $\frac{N}{2}$ texts of the dataset. The proportion $\frac{N}{2}$ was elected once smaller ones (such as $\frac{N}{3}$ or $\frac{N}{4}$, being $N$ the total number of books in the dataset) evoked archaic words and words not belonging to the vernacular of the English language, and higher proportions, on the contrary, led to poorer results on the experiments.

Considering each entry in $S$, we computed its frequency for all books in the dataset, resulting, in the end, in a $219 \times 3585$ matrix of frequencies, henceforward called $M$. Next, the rows of $M$ — each representing a book —   were standardized and transformed according to two approaches: LDA and SemAxis, the results being cross-validated through leave-one-out. As shown in Figure~\ref{fig:ldasemaxisbow}, such processing led to a visual separation both in~\ref{subfig:bowlda} and~\ref{subfig:bowsemaxis}, giving evidence that the bag-of-words model can provide a good — although not exact — split between the two studied categories.

\begin{figure}[!htbp]
        \begin{subfigure}{0.8\linewidth}
            \includegraphics[width=\linewidth]{ldaBoW_normalized.pdf}
            \caption{LDA}
            \label{subfig:bowlda}
        \end{subfigure}

        \vspace*{0.2cm}

        \begin{subfigure}{0.8\linewidth}
            \includegraphics[width=\linewidth]{semaxisBoW_normalized.pdf}
            \caption{SemAxis}
            \label{subfig:bowsemaxis}
        \end{subfigure}

    \caption{Kernel density estimation of the 219 investigated literary works, considering (a) LDA projection and (b) SemAxis projection of $M$. The colors on the plot correspond to the categories of the books, blue indicating \emph{success} and orange \emph{other}.}
    \label{fig:ldasemaxisbow}
\end{figure}

Moreover, to quantitatively assess the obtained separation, $M$ was used as input to supervised classification methods (videlicet: KNN, logistic regression, naive Bayes, decision tree, random forest, and SVM). We also applied leave-one-out and $k$-fold (taking k = 10) cross-validation methods and considered both the standardized and the non-standardized versions of $M$ (the standardized version denoted by $\hat{M}$). Here, we adopted the standard hyperparameters of each model — in other words, not involving tuning operations.

As shown in Table~\ref{tab:modelaccuracybow}, the linear regression model resulted as the best choice for grasping discrepancies between classes, leading to an average classification accuracy of 0.75, either for LOO or $k$-fold cross-validation. This result shows that the approach is apt — to a reasonable extent — to identify successful literary works. Furthermore, it is worth observing that the standardization positively impacts the outcomes, leading to performances as good as or better than the non-standardized case in ten out of twelve scenarios — languishing only the accuracy of the KNN model.

\begin{table}[!htbp]
\caption{Classification accuracy for different models and arrangements (i.e., whether were employed $M$ or $\hat{M}$ and leave-one-out or $k$-fold cross-validation). Highlighted in bold is the best result for each configuration.}
\renewcommand{\arraystretch}{1.3}
\centering
\begin{tabular}{@{} l
                @{\hspace*{0mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{0mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{}}
\toprule[1.0pt]
\multirow{2}{*}{}&&\multicolumn{2}{c}{$M$}&&\multicolumn{2}{c}{$\hat{M}$}\\
\cline{3-4}\cline{6-7}
&&LOO&10-fold&&LOO&10-fold\\
\hline
KNN&&0.64&$0.64\pm0.11$&&0.58&$0.56\pm0.10$\\
LR&&0.65&$0.64\pm0.13$&&\textbf{0.75}&$\boldsymbol{0.75\pm0.09}$\\
NB&&0.63&$0.62\pm0.11$&&0.63&$0.62\pm0.11$\\
DT&&0.65&$0.58\pm0.14$&&0.65&$0.58\pm0.14$\\
RF&&\textbf{0.68}&$\boldsymbol{0.68\pm0.11}$&&0.68&$0.68\pm0.11$\\
SVM&&0.66&$0.63\pm0.11$&&0.72&$0.74\pm0.09$\\
\bottomrule[1.0pt]
\end{tabular}
\label{tab:modelaccuracybow}
\end{table}%

In addition, we retrieved the 40 words of $S$ with preponderant impact onto the SemAxis projection, aiming at analyzing what sort of vocable seems to be characteristic in best sellers and in the \emph{other} books. As presented in Table~\ref{tab:semaxisrelevantwords}, the most meaningful words for successful books encompass six adjectives, nine nouns, one adverb, and four verbs; for the non-best-seller books, we have three adjectives, seven nouns, one adverb, and nine verbs. Similarly to a result formerly reported in~\cite{ashok2013success}, words referring to body parts (such as \emph{eye}, \emph{face}, and \emph{hand}) play a central role in less successful titles. Furthermore, none of the 20 most relevant terms for successes ranks among the 40 most frequent words of $S$ — however, when analyzing the non-best-seller books,  the principal words \emph{eye}, \emph{face}, \emph{hand}, and \emph{back} represent, respectively, the 5th, 6th, 10th, and 12th most common words of the dataset.

\begin{table}[!htbp]
\caption{Forty most significant words to the SemAxis projection discrimination between best sellers and others. The importance of each term for the method dictates its allocation order in the table: the element on the first row and the first column (of success/other) is the most important for the class; the one on the second row and the second column is the second most important; and so forth.}
\centering
\renewcommand{\arraystretch}{1.3}
\begin{tabular}[t]{@{} l
                  @{\hspace*{15mm}}l
                  @{\hspace*{10mm}}l
                  @{\hspace*{10mm}}l
                  @{\hspace*{10mm}}l
                  @{}}
\toprule[1.0pt]
&\multicolumn{4}{c}{Vocables}\\
\cline{2-5}
\multirow{5}{*}{Success}&\emph{ordinary}&\emph{evidence}&\emph{motive}&\emph{exhibit}\\
&\emph{grey}&\emph{substance}&\emph{improve}&\emph{copy}\\
&\emph{instruction}&\emph{contain}&\emph{examination}&\emph{practice}\\
&\emph{accordingly}&\emph{teacher}&\emph{numerous}&\emph{interesting}\\
&\emph{school}&\emph{large}&\emph{attach}&\emph{average}\\
\hline
\multirow{5}{*}{Other}&\emph{breath}&\emph{drop}&\emph{draw}&\emph{sharply}\\
&\emph{face}&\emph{eye}&\emph{turn}&\emph{stun}\\
&\emph{break}&\emph{hand}&\emph{push}&\emph{back}\\
&\emph{reckless}&\emph{caught}&\emph{arm}&\emph{shake}\\
&\emph{tone}&\emph{instant}&\emph{quick}&\emph{glance}\\
\bottomrule[1.0pt]
\end{tabular}
\label{tab:semaxisrelevantwords}
\end{table}%

\subsection{Doc2vec analysis}

The second experiment evaluates whether doc2vec's representation of literary works can grasp the dissimilarity between the two analyzed classes. With this aim, we instantiated D2V models with 32, 64, 128, and 256 dimensions (a feature commonly called vector size, hereafter referred to as $\#_2D$). We also set the minimum word count to 1 (to ignore all words with a total frequency lower than one), the window (maximum distance between the current and predicted term within a sentence) to 5, and the epochs (number of iterations over the corpus) to 40. Lastly, the model training occurred using all 219 instances of the dataset.

Next, each model vector (henceforth called $D$) — a piece representing a different book — was transformed employing LDA and SemAxis techniques along with leave-one-out cross-validation, yielding the results shown in Figure~\ref{fig:ldasemaxisd2v}. As can be observed, the method was able to characterize best seller and non-best-seller works in a contrasting fashion, both in~\ref{subfig:d2vlda} and~\ref{subfig:d2vsemaxis}. This result shows that it is possible to emphasize the differences between the two classes in two noticeably distinct approaches (either BoW or D2V).

\begin{figure}[!htbp]
        \begin{subfigure}{0.8\linewidth}
            \includegraphics[width=\linewidth]{ldaD2V_normalized.pdf}
            \caption{LDA}
            \label{subfig:d2vlda}
        \end{subfigure}

        \vspace*{0.2cm}

        \begin{subfigure}{0.8\linewidth}
            \includegraphics[width=\linewidth]{semaxisD2V_normalized.pdf}
            \caption{SemAxis}
            \label{subfig:d2vsemaxis}
        \end{subfigure}

    \caption{Kernel density estimation of the 219 investigated literary works, considering (a) LDA projection and (b) SemAxis projection of D2V representation  (adopting $\#_2D = 64$). The colors on the plot correspond to the categories of the books, blue indicating \emph{success} and orange \emph{other}.}
    \label{fig:ldasemaxisd2v}
\end{figure}

Furthermore, we used $D$ as the supervised classification methods' input to quantitatively assess the obtained separation. The models used here were the same as those applied in the BoW experiment, and we also considered LOO and 10-fold cross-validations and both the standardized and the non-standardized versions of $D$ (the standardized version denoted by $\hat{D}$). The chosen models' hyperparameters were the standard ones.

As shown in Tables~\ref{tab:modelaccuracyd2vl1o} and~\ref{tab:modelaccuracyd2vkfold}, naive Bayes was the model that best performed the task of distinguishing classes considering the D2V representation, leading to a maximum classification accuracy of $0.71$ for LOO and $0.72\pm0.12$ for 10-fold. In the LOO version, $\#_2D = 32$ raises the best results, while $\#_2D = 256$ performs better for 10-fold. Although the standardization did not affect the naive Bayes classifier, it led to the same or slightly better outcomes for the others — the exceptions being some arrangements, namely KNN (for $\#_2D = 128$) and LR (for $\#_2D = 64$ and $128$) for the LOO version and KNN ($\#_2D = 256$), LR ($\#_2D = 256$), and SVM ($\#_2D = 64$) for 10-fold.

\begin{table}[!htbp]
\caption{Classification accuracy for LOO cross-validation combined with different models and arrangements (i.e., whether were employed $D$ or $\hat{D}$ and with which D2V vector size: 32, 64, 128, or 256). Highlighted in bold is the best result for each configuration.}
\renewcommand{\arraystretch}{1.3}
\centering
\begin{tabular}{@{} l
                @{\hspace*{0mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{0mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{}}
\toprule[1.0pt]
\multirow{3}{*}{}&\multicolumn{10}{c}{LOO}\\
\cline{3-11}
&&\multicolumn{4}{c}{$D$}&&\multicolumn{4}{c}{$\hat{D}$}\\
\cline{3-6}\cline{8-11}
&&32&64&128&256&&32&64&128&256 \\
\hline
KNN&&0.65&0.65&0.65&0.67&&0.65&0.65&0.64&0.67\\
LR&&0.68&0.65&0.63&0.66&&0.70&0.64&0.61&0.67\\
NB&&\textbf{0.71}&\textbf{0.68}&\textbf{0.69}&\textbf{0.68}&&\textbf{0.71}&\textbf{0.68}&\textbf{0.69}&\textbf{0.68}\\
DT&&0.48&0.58&0.53&0.40&&0.48&0.58&0.53&0.40\\
RF&&0.66&0.63&0.67&0.64&&0.66&0.63&0.67&0.64\\
SVM&&0.68&0.63&0.62&0.63&&0.69&0.66&0.65&0.64\\
\bottomrule[1.0pt]
\end{tabular}
\label{tab:modelaccuracyd2vl1o}
\end{table}%

\begin{table}[!htbp]
\caption{Classification accuracy for 10-fold cross-validation combined with different models and arrangements (i.e., whether were employed $D$ or $\hat{D}$ and with which D2V vector size: 32, 64, 128, or 256). Highlighted in bold is the best result for each configuration.}
\renewcommand{\arraystretch}{1.3}
\centering
\begin{tabular}{@{} l
                @{\hspace*{0mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{0mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{}}
\toprule[1.0pt]
\multirow{3}{*}{}&&&\multicolumn{8}{c}{10-fold}\\
\cline{3-11}
&&&\multicolumn{8}{c}{$D$}\\
\cline{4-11}
&&&32&&64&&&128&&256 \\
\hline
KNN&&&$0.67\pm0.08$&&$0.63\pm0.06$&&&$0.65\pm0.13$&&$0.69\pm0.11$\\
LR&&&$0.66\pm0.13$&&$0.67\pm0.10$&&&$0.59\pm0.08$&&$0.71\pm0.09$\\
NB&&&$\boldsymbol{0.68\pm0.10}$&&$\boldsymbol{0.70\pm0.13}$&&&$\boldsymbol{0.70\pm0.08}$&&$\boldsymbol{0.72\pm0.12}$\\
DT&&&$0.58\pm0.12$&&$0.52\pm0.08$&&&$0.56\pm0.10$&&$0.50\pm0.06$\\
RF&&&$\boldsymbol{0.68\pm0.10}$&&$0.65\pm0.10$&&&$0.66\pm0.09$&&$0.70\pm0.10$\\
SVM&&&$0.66\pm0.13$&&$0.66\pm0.09$&&&$0.57\pm0.07$&&$0.69\pm0.10$\\
&&&\multicolumn{8}{c}{$\hat{D}$}\\
\cline{4-11}
&&&32&&64&&&128&&256 \\
\hline
KNN&&&$\boldsymbol{0.68\pm0.11}$&&$0.64\pm0.05$&&&$0.66\pm0.10$&&$0.68\pm0.08$\\
LR&&&$\boldsymbol{0.68\pm0.11}$&&$0.69\pm0.12$&&&$0.60\pm0.08$&&$0.70\pm0.10$\\
NB&&&$\boldsymbol{0.68\pm0.10}$&&$\boldsymbol{0.70\pm0.13}$&&&$\boldsymbol{0.70\pm0.08}$&&$\boldsymbol{0.72\pm0.12}$\\
DT&&&$0.58\pm0.12$&&$0.52\pm0.08$&&&$0.56\pm0.10$&&$0.50\pm0.06$\\
RF&&&$\boldsymbol{0.68\pm0.10}$&&$0.65\pm0.10$&&&$0.66\pm0.09$&&$0.70\pm0.10$\\
SVM&&&$0.67\pm0.14$&&$0.64\pm0.11$&&&$0.57\pm0.08$&&$0.70\pm0.10$\\
\bottomrule[1.0pt]
\end{tabular}
\label{tab:modelaccuracyd2vkfold}
\end{table}%

\subsection{Are the subjects being grasped by the approaches?}

What if the above mentioned approaches are only grasping the subjects of the books? That would be a valid inquiry once we did not regard this type of information during the construction of the database. To assess this possibility, we retrieved the list of subjects of each book provided by the Gutenberg platform and then analyzed the ten most common ones in the dataset. In Figure~\ref{fig:subjectVsSemaxis}, we plot those subjects against the SemAxis projection of the books' D2V vector representation (using $\#_2D = 64$), stratifying the results by category. As one can see, the only subjects with a representative number of instances are PS and PR, which also seem to explain the separation obtained through the D2V method to some degree.

\begin{figure}[!htbp]
    \includegraphics[width=\linewidth]{subjectsvssemaxis.pdf}
    \caption{On the y-axis, the ten most common subjects in the dataset. On the x-axis, the SemAxis projection of the books' D2V representation (adopting $\#_2D = 64$). The colors on the plot correspond to the categories of the books, blue indicating \emph{success} and orange \emph{other}.}
    \label{fig:subjectVsSemaxis}
\end{figure}

PR and PS are classifications used by the Library of Congress \cite{loc2013} to catalog English and British literature, respectively. In our case, the PR subject represents 102 instances of the dataset, 34 best sellers and 68 non-best-seller works. The PS one, by contrast, encompasses 98 books — 72 best sellers and 26 other types of works. In this manner, as, in principle, the success category is the only one with a limited number of instances (given its criteria), we created a new dataset (with 72 successes and 72 others, embracing the same standards stated in the creation of the former dataset) with only literary works belonging to subject PS. Then, it was selected, from the Gutenberg platform, 46 new non-best-seller titles. Using this current dataset, we repeated the previous experiments, aiming at understanding whether the fact that a book belongs to English or British literature was enough to explain the separation provided by the BoW and D2V methods. The results are presented and discussed below.

\subsubsection{Bag-of-words analysis}
Similarly to the previous experiment, we considered the set $S$, built based on the 3,257 different words that appeared in the text of at least $\frac{N}{2}$ books of the PS dataset. Then, we calculated their frequencies, resulting in a $144\times3257$ dimension matrix, (henceforth called $M_{PS}$). Next, $M_{PS}$'s rows were standardized, transformed using LDA and SemAxis, and verified via LOO cross-validation, as shown in Figure~\ref{fig:ldasemaxisbow_PS}. It is possible to observe that the separation between the classes is still perceptible — both in~\ref{subfig:bowlda_PS} and~\ref{subfig:bowsemaxis_PS} — although now solely English literature is being considered.

\begin{figure}[!htbp]
        \begin{subfigure}{0.8\linewidth}
            \includegraphics[width=\linewidth]{ldaBoW_normalized_PS.pdf}
            \caption{LDA}
            \label{subfig:bowlda_PS}
        \end{subfigure}

        \vspace*{0.2cm}

        \begin{subfigure}{0.8\linewidth}
            \includegraphics[width=\linewidth]{semaxisBoW_normalized_PS.pdf}
            \caption{SemAxis}
            \label{subfig:bowsemaxis_PS}
        \end{subfigure}

    \caption{Kernel density estimation of the 144 investigated literary works belonging to the PS subject, considering (a) LDA projection and (b) SemAxis projection of $M_{PS}$. The colors on the plot correspond to the categories of the books, blue indicating \emph{success} and orange \emph{other}.}
    \label{fig:ldasemaxisbow_PS}
\end{figure}

To quantitatively analyze the separation of the categories, we performed supervised classification methods using the standardized and non-standardized versions of $M_{PS}$ as input. The applied models, hyperparameters, and cross-validation methods were the same as the former experimentation. As shown in Table~\ref{tab:modelaccuracybow_PS}, the random forest model was the best option to distinguish the best seller and other instances, leading to an average accuracy of 0.71 — both in LOO and 10-fold cross-validations. The standardization did not affect the results. Even though this is a lower accuracy than that obtained with the previous dataset (which led to the highest accuracy of 0.75), it is worth mentioning that the PS dataset has 35\% fewer instances than the other, which makes us expect lower accuracies and higher standard deviations. Thus, it is possible to state that the BoW method is not classifying the corpus instances based predominantly or solely on their literary class.

\begin{table}[!htbp]
\caption{Classification accuracy for different models and arrangements (i.e., whether were employed $M_{PS}$ or $\hat{M}_{PS}$ and leave-one-out or $k$-fold cross-validation), considering the PS dataset. Highlighted in bold is the best result for each configuration.}
\renewcommand{\arraystretch}{1.3}
\centering
\begin{tabular}{@{} l
                @{\hspace*{0mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{0mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{}}
\toprule[1.0pt]
\multirow{2}{*}{}&&\multicolumn{2}{c}{$M_{PS}$}&&\multicolumn{2}{c}{$\hat{M}_{PS}$}\\
\cline{3-4}\cline{6-7}
&&LOO&10-fold&&LOO&10-fold\\
\hline
KNN&&0.62&$0.60\pm0.09$&&0.63&$0.62\pm0.13$\\
LR&&0.64&$0.63\pm0.15$&&0.69&$0.67\pm0.13$\\
NB&&0.65&$0.65\pm0.15$&&0.65&$0.65\pm0.15$\\
DT&&0.62&$0.57\pm0.14$&&0.62&$0.57\pm0.14$\\
RF&&\textbf{0.71}&$\boldsymbol{0.71\pm0.14}$&&\textbf{0.71}&$\boldsymbol{0.70\pm0.12}$\\
SVM&&0.61&$0.61\pm0.16$&&0.66&$0.61\pm0.16$\\
\bottomrule[1.0pt]
\end{tabular}
\label{tab:modelaccuracybow_PS}
\end{table}%

\subsubsection{Doc2vec analysis}

The D2V models were instantiated in the context of the new dataset with the same hyperparameters as the previous tests — the only exception being that we trained the model using 144 books instead of 219. We repeated the former configurations (vector sizes 32, 64, 128, and 256 and LOO and 10-fold cross-validations) and adopted the standardized and the non-standardized versions of the model vectors — called $\hat{D}_{PS}$ and $D_{PS}$, respectively. Figure~\ref{fig:ldasemaxisd2v_PS} shows the results of the transformations of the model vectors via LDA and SemAxis. The split between best seller and non-best-seller works was again observed, suggesting that the method is insensitive to the literary class.

\begin{figure}[!htbp]
        \begin{subfigure}{0.8\linewidth}
            \includegraphics[width=\linewidth]{ldaD2V_normalized_PS.pdf}
            \caption{LDA}
            \label{subfig:d2vlda_PS}
        \end{subfigure}

        \vspace*{0.2cm}

        \begin{subfigure}{0.8\linewidth}
            \includegraphics[width=\linewidth]{semaxisD2V_normalized_PS.pdf}
            \caption{SemAxis}
            \label{subfig:d2vsemaxis_PS}
        \end{subfigure}

    \caption{Kernel density estimation of the 144 investigated literary works belonging to the PS subject, considering (a) LDA projection and (b) SemAxis projection of D2V representation  (adopting $\#_2D_{PS} = 64$). The colors on the plot correspond to the categories of the books, blue indicating \emph{success} and orange \emph{other}.}
    \label{fig:ldasemaxisd2v_PS}
\end{figure}

The quantitative assessment using supervised classification led to the results shown in Tables~\ref{tab:modelaccuracyd2vl1o_PS} and~\ref{tab:modelaccuracyd2vkfold_PS}. From Table~\ref{tab:modelaccuracyd2vkfold_PS}, it is possible to conclude that the models that best performed for LOO cross-validation were logistic regression and naive Bayes — the highest accuracy (0.67) given by the latter, with $\#_2D_{PS} = 256$. In this case, the standardization did not contribute to performance improvement only in the case of the KNN model. From Table~\ref{tab:modelaccuracyd2vkfold_PS}, we conclude that no model stood out for 10-fold, with the highest accuracy of 0.67 given by the SVM model, with $\#_2D_{PS} = 32$ and non-standardized input. The standardization process induced better results on six distinct occasions, although the best-obtained accuracy counts on a non-standardized vector.

\begin{table}[!htbp]
\caption{Classification accuracy for LOO cross-validation combined with different models and arrangements (i.e., whether $D_{PS}$ or $\hat{D}_{PS}$ were employed and with which D2V vector size: 32, 64, 128, or 256), considering the PS dataset. Highlighted in bold is the best result for each configuration.}
\renewcommand{\arraystretch}{1.3}
\centering
\begin{tabular}{@{} l
                @{\hspace*{0mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{0mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{}}
\toprule[1.0pt]
\multirow{3}{*}{}&\multicolumn{10}{c}{LOO}\\
\cline{3-11}
&&\multicolumn{4}{c}{$D_{PS}$}&&\multicolumn{4}{c}{$\hat{D}_{PS}$}\\
\cline{3-6}\cline{8-11}
&&32&64&128&256&&32&64&128&256 \\
\hline
KNN&&0.56&0.59&0.60&0.56&&0.51&0.57&0.56&0.55\\
LR&&\textbf{0.62}&\textbf{0.62}&0.58&0.60&&\textbf{0.62}&\textbf{0.65}&0.60&0.60\\
NB&&0.60&\textbf{0.62}&\textbf{0.62}&\textbf{0.67}&&0.60&0.62&\textbf{0.62}&\textbf{0.67}\\
DT&&0.51&0.51&0.44&0.53&&0.51&0.51&0.44&0.53\\
RF&&0.53&0.60&0.58&0.56&&0.53&0.60&0.58&0.56\\
SVM&&\textbf{0.62}&0.61&0.60&0.58&&\textbf{0.62}&0.61&\textbf{0.62}&0.58\\
\bottomrule[1.0pt]
\end{tabular}
\label{tab:modelaccuracyd2vl1o_PS}
\end{table}%

\begin{table}[!htbp]
\caption{Classification accuracy for 10-fold cross-validation combined with different models and arrangements (i.e., whether $D_{PS}$ or $\hat{D}_{PS}$ were employed and with which D2V vector size: 32, 64, 128, or 256), considering the PS dataset. Highlighted in bold is the best result for each configuration.}
\renewcommand{\arraystretch}{1.3}
\centering
\begin{tabular}{@{} l
                @{\hspace*{0mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{0mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{\hspace*{5mm}}c
                @{}}
\toprule[1.0pt]
\multirow{3}{*}{}&&&\multicolumn{8}{c}{10-fold}\\
\cline{3-11}
&&&\multicolumn{8}{c}{$D_{PS}$}\\
\cline{4-11}
&&&32&&64&&&128&&256 \\
\hline
KNN&&&$0.59\pm0.13$&&$0.58\pm0.14$&&&$0.49\pm0.12$&&$0.60\pm0.09$\\
LR&&&$\boldsymbol{0.66\pm0.13}$&&$0.54\pm0.06$&&&$0.59\pm0.12$&&$\boldsymbol{0.66\pm0.08}$\\
NB&&&$0.61\pm0.08$&&$0.59\pm0.12$&&&$0.58\pm0.14$&&$0.61\pm0.14$\\
DT&&&$0.45\pm0.13$&&$0.49\pm0.10$&&&$\boldsymbol{0.66\pm0.13}$&&$0.48\pm0.09$\\
RF&&&$0.59\pm0.10$&&$\boldsymbol{0.60\pm0.10}$&&&$0.57\pm0.17$&&$0.59\pm0.13$\\
SVM&&&$0.65\pm0.12$&&$0.49\pm0.10$&&&$0.58\pm0.12$&&$0.63\pm0.10$\\
&&&\multicolumn{8}{c}{$\hat{D}_{PS}$}\\
\cline{4-11}
&&&32&&64&&&128&&256 \\
\hline
KNN&&&$0.60\pm0.13$&&$\boldsymbol{0.61\pm0.16}$&&&$0.53\pm0.07$&&$0.59\pm0.08$\\
LR&&&$0.66\pm0.12$&&$0.58\pm0.09$&&&$0.59\pm0.12$&&$\boldsymbol{0.65\pm0.09}$\\
NB&&&$0.61\pm0.08$&&$0.59\pm0.12$&&&$0.58\pm0.14$&&$0.61\pm0.14$\\
DT&&&$0.45\pm0.13$&&$0.49\pm0.10$&&&$\boldsymbol{0.66\pm0.13}$&&$0.48\pm0.09$\\
RF&&&$0.59\pm0.10$&&$0.60\pm0.10$&&&$0.57\pm0.17$&&$0.59\pm0.13$\\
SVM&&&$\boldsymbol{0.67\pm0.11}$&&$0.53\pm0.07$&&&$0.55\pm0.14$&&$0.63\pm0.10$\\
\bottomrule[1.0pt]
\end{tabular}
\label{tab:modelaccuracyd2vkfold_PS}
\end{table}%

For the 219-instances dataset, the best-achieved accuracy was 0.72. Again, we expected a drop in the accuracy, as the new dataset has 35\% lesser instances than the other. Thus, just as in the BoW method, it is possible to infer that the separation between classes obtained in the D2V approach does not rely on whether a book belongs to English or British literature.

\section{Conclusions}
\label{sec:conclusions}

{The study of characteristics leading to literary pieces becoming best sellers constitutes an intriguing and challenging research subject. The present work addressed this issue while considering aspects derived from the full content of a list of more and less successful books retrieved from the Gutenberg Project, based on the best seller lists of Publishers Weekly. Several alternative content representation, standardization, visualization, and classification approaches were considered, as summarized in the diagram shown in Figure~\ref{fig:schematic}.}

We started our analysis by examining the data using visualization techniques. The visualization enabled a preliminary direct inspection of the embedding by looking at a single axis that maximizes the separation between best sellers and ordinary books. Specifically, we employed SemAxis and LDA techniques -- the first providing better discrimination between classes than the latter, both for bag-of-words and doc2vec representations. Furthermore, SemAxis provided means that helped to: (i) understand the most characteristic words in best sellers and non-best-sellers; and (ii) check if the respective success was related to the subjects of the books (e.g., love stories, adventure stories, fiction, among others). {In line with earlier work~\cite{ashok2013success}, words related to body parts (like \emph{face}, \emph{eye}, and \emph{hand}) played a central role in non-best-seller books, while more varied and less common vocables (such as \emph{ordinary}, \emph{accordingly}, and \emph{examination}) were characteristic of more successful books. Moreover, we found no evidence that the subject of the books impacted the class discrimination obtained.}

For the classification tasks, we tested two strategies for preprocessing the two distinct representations: (i) standardizing and (ii) non-standardizing the embeddings. Then, we evaluated the proposed representations via different classifiers (namely: KNN, LR, NB, DT, RF, and SVM). The best-obtained result was acquired with the complete dataset (219 books), using the LR classifier with the standardized bag-of-words representation. In this case, the final classification accuracy was 0.75. Still dealing with the complete set, the best accuracy obtained for D2V embedding was 0.72, combining the standardized representation with the NB model. For the dataset considering only the PS subject (144 books), the bag-of-words approach throughput the most promising results for the standardized data inputted in the RF classifier. The D2V representation, in contrast, returned better outcomes for the standardized data combined with the NB classifier. These results agree with the tendency of the two classes' separation found in the visualization analysis. Interestingly, the standardization did not affect the results significantly in the doc2vec approach for both datasets.

The reported methodology and results pave the way for several related studies, some of which are described as follows. Firstly, it would be interesting to adapt the reported method to other types of embeddings, being one example the BERT transformer modified to work with long texts. Secondly, it would be interesting to consider the described approach for better understanding: (i) other types of documents, such as scientific books and articles, and (ii) additional types of artistic production, including music, poetry, and theater. Lastly, another point that could be explored concerns the explanation of reasons why some literary works become best sellers and others do not.

\section*{Acknowledgments}
Giovana D. da Silva acknowledges FAPESP for sponsorship (grant no. 2021/01744-0). Bárbara C. e Souza acknowledges CAPES Brazil for sponsorship -- Finance Code 001. Luciano da F. Costa thanks CNPq (grant no. 307085/2018-0) and FAPESP (grant 15/22308-2). Diego R. Amancio acknowledges financial support from CNPq (grant no. 311074/2021-9) and FAPESP (no. 2020/06271-0).

\newpage

\bibliographystyle{ieeetr}
\bibliographystyle{abbrv}


\section{Introduction}\label{s:intro}
The investigation of out-of-equilibrium quantum systems is an extremely rich and active topic in modern statistical and condensed matter physics, which aims at clarifying fundamental concepts like thermalization and the emergence of the ensemble description~\cite{srednicki1994,polkovnikov2011,rigol2008,gogolin2016,dalessio2016,deutsch2018}. Theoretical advances in recent years have highlighted the importance of quantum information scrambling~\cite{sekino2008,iyoda2018,hosur2016} as a general feature associated to the relaxation of observables: Locally encoded information spreads throughout the system during the dynamics, and, while formally it is not lost, it cannot be recovered through any local measurement at late times. A paradigmatic example of a system exhibiting scrambling behavior is the Sachdev--Ye--Kitaev (SYK) model~\cite{kitaev2015,maldacena2016,gu2020,chowdhury2022}, which received much interest in recent research due to its relevance in multiple branches of physics, ranging from the study of black holes~\cite{kitaev2015,davison2017,kitaev2018,sachdev2019,sarosi2018} to non-Fermi liquids~\cite{sachdev2015,song2017,cha2020,sachdev2010}. There are multiple versions of the SYK model, parameterized by an even number $q$ and referred to as SYK$_q$. Each one consists of fermionic modes coupled via all-to-all disordered $q/2$-body interaction amplitudes. For $q\geq 4$, the model manifests quantum chaotic behavior~\cite{kitaev2015,polchinski2016,garciagarcia2018}, diagnosed by the presence of a quantum Lyapunov exponent in out-of-time-order correlators (OTOCs)~\cite{maldacena2016bis,hashimoto2017,kobrin2021}. Specifically, this exponent is found to saturate its theoretical upper bound~\cite{maldacena2016bis}, implying that quantum scrambling occurs as rapidly as possible. While multiple properties can be derived exactly in the thermodynamic limit~\cite{kitaev2015,sachdev2015,chowdhury2022}, studying the model at finite size remains a challenging yet essential problem, especially in view of future experimental investigations that have recently been proposed~\cite{danshita2017,garciaalvarez2017,pikulin2017,chew2017,chen2018,luo2019,wei2021}.

In this work, we investigate the disorder-averaged dynamics of observables in the SYK model, highlighting the manifestation of a counter-intuitive symmetry. In contrast to the scrambling nature of the system, we prove rigorously that disorder-averaged operators do not increase in complexity throughout the time-evolution, leading to complete absence of operator growth on average. This corroborates the results of our previous work~\cite{bandyopadhyay2021}, which demonstrated numerically a non-trivial universality with respect to the choice of the initial state in the dynamics of some observables. The argument only applies to single-time observables, and, in particular, it does not hold for OTOCs: All correlators that can be used to diagnose scrambling and chaotic behavior are unaffected by our discussion, which resolves the apparent contradiction of our analysis with the well-known results from the literature. While the lack of operator growth does a priori not apply to individual realizations of the model, it is expected to manifest in self-averaging single-time physical quantities for sufficiently large system sizes. In addition, working in a framework typical of open quantum systems, we leverage these results to show how to approximate the average dynamics by means of a cumulant expansion of the effective dynamical map. Specifically, owing to the absence of operator growth, the functions approximating the dynamics of observables can be computed analytically for arbitrary system sizes. Finally, we benchmark the performance of our approximation scheme through comparison with results found using exact diagonalization.

The paper is organized as follows. Section~\ref{s:model} presents the model and the superoperator framework adopted throughout the work. In Sec.~\ref{s:operator_size}, we introduce the concept of operator size, and we prove its conservation throughout the dynamics in ensemble average. Moreover, we discuss how this result does not conflict with well-known scrambling properties of the model. Then, in Sec.~\ref{s:cumulant_expansion} we present the cumulant expansion technique that we implement to approximate the dynamics. Our theoretical findings are tested in Sec.~\ref{s:application}, in which we investigate operator growth, and we probe the performance of the cumulant expansion using exact diagonalization simulations for the SYK$_4$ model. We summarize our findings in Sec.~\ref{s:conclusions}, and we discuss their possible applications, as well as future lines of research. Appendices expand the discussion of the main work, both by providing detailed descriptions of some proofs, and by showing explicit analytic results.

\section{Model and focus of our study}\label{s:model}

We focus on the family of SYK models of $N$ complex fermion modes on a fully connected lattice, defined by the Hamiltonian
\begin{subequations}\label{hamiltonian}
\begin{equation}\label{hamiltonian-a}
    \hat{H} = \mathcal{K}_q \sum_{\substack{i_1<\dots<i_{q/2}\\j_1<\dots<j_{q/2}}} J_{i_1,\dots,i_{q/2};j_1,\dots,j_{q/2}} \hat{c}^\dag_{i_1}\dots\hat{c}^\dag_{i_{q/2}}\hat{c}_{j_1}\dots\hat{c}_{j_{q/2}},
    \end{equation}\begin{equation}
    \mathcal{K}_q = \sqrt{\frac{(q/2)!(q/2-1)!}{N^{q-1}}},
    \end{equation}
\end{subequations}
where $q$ is an even integer number, and $J_{i_1,\dots,i_{q/2};j_1,\dots,j_{q/2}}$ are complex Gaussian independent random couplings with the following statistical properties ($\eave{...}$ denotes disorder averaging):
\begin{subequations}\label{disorder_properties}
\begin{align}
\label{disorder_properties-a}
& \eave{J_{i_1,\dots,i_{q/2};j_1,\dots,j_{q/2}}}=0,\\
\label{disorder_properties-b}
& \eave{\left(J_{i_1,\dots,i_{q/2};j_1,\dots,j_{q/2}}\right)^2}=0 \\
    & \quad \text{if } \{i_1,\dots,i_{q/2}\}\neq\{j_1,\dots,j_{q/2}\}, \nonumber\\
    \label{disorder_properties-c}
 & \eave{\left|J_{i_1,\dots,i_{q/2};j_1,\dots,j_{q/2}}\right|^2} = J^2.
\end{align}
\end{subequations}
Notably, for $q>2$ the only integral of motion of any individual realization is the total charge $\hat{Q}=\sum_i \hat{n}_i$.

Our analysis relies on the study of the dynamics of observable expectation values, namely
\begin{equation}
    \braket{\hat{W}(t)}  = \Tr\left(\hat{\rho}_0e^{i\hat{H}t} \hat{W} e^{-i\hat{H}t} \right),
\end{equation}
where typically (but not necessarily) $\hat{\rho}_0 = \ket{\psi_0}\bra{\psi_0}$ is a pure-state density matrix:
This corresponds to preparing the system in a disorder-independent initial state $\ket{\psi_0}$, performing a quantum quench to the SYK$_q$ Hamiltonian, and tracking the time-evolution of physical quantities, as carried out, e.g., in~\cite{bandyopadhyay2021}.
Specifically, since the Hamiltonian involves randomness, we focus on the disorder-averaged evolution $\eave{\braket{\hat{W}(t)}}$. A well-known feature of disordered systems is that disorder-averaging introduces mixing, making the average dynamics non-unitary~\cite{kropf2016}. As a consequence, it becomes natural to adopt a superoperator-based framework, typical of dissipative open quantum systems~\cite{minganti2018,manzano2020}. We introduce the Liouvillian $\mathcal{L} \bullet = -i\comm{\hat{H}}{\bullet}$, as well as the superoperators $\ell_\alpha \bullet = -i \comm{\hat{h}_\alpha}{\bullet}$, where $\alpha = \{i_1,\dots,i_{q/2};j_1,\dots,j_{q/2}\}$ is a multi-index, and $\hat{h}_\alpha = \hat{c}^\dag_{i_1}\dots\hat{c}^\dag_{i_{q/2}}\hat{c}_{j_1}\dots\hat{c}_{j_{q/2}}$ appears in the Hamiltonian of Eq.~\eqref{hamiltonian-a}. With this notation, the ensemble-averaged dynamics reads
\begin{equation}\label{observable_dyn}
    \eave{\braket{\hat{W}(t)}} = \Tr\left(\hat{\rho}_0 \eave{e^{-\mathcal{L} t}}\hat{W} \right).
\end{equation}
In what follows, we provide a detailed analysis of the time evolution of such disorder-averaged expectation values.

\section{Operator size symmetry in the disorder-averaged ensemble}\label{s:operator_size}

In this section, we introduce the notion of operator size, and we show that it plays the role of a conserved quantum number for the average dynamics of operators. This leads us to the conclusion that no operator growth takes place for expectation values of observables. As we discuss, this result does not hold for more sophisticated correlators, e.g., OTOCs, which are able to diagnose quantum scrambling.

\subsection{Operator size symmetry}\label{ss:size_sym}

An analysis of the spectral properties of the disorder-averaged time-evolution superoperator allows us to determine some exact features of the dynamics. Specifically, we will prove further below that the average dynamics distinguishes different operator sizes. Loosely speaking, we say that an operator consisting of the product of $m$ creation and $n$ annihilation operators has operator size $(m,n)$~\cite{carrega2021,roberts2018}. This definition is appropriate if all fermionic operators refer to different lattice sites, but we must analyze more carefully the cases involving one or more number operators $\hat{n}_i$. In any charge-conserved sector of the Hilbert space, the operator $\hat{Q}=\sum_i \hat{n}_i$ is proportional to the identity, and thus its size is $(0,0)$. As a consequence, the size of the number operators $\hat{n}_i$ is not simply $(1,1)$, because each of them has a finite overlap with the identity operator. Despite this complication, it is possible to define the concept of operator size for the most general combination of creation and annihilation operators (see the discussion of App.~\ref{a:op_size}). Any linear combination of operators sharing a single common size also maintains that same size. This allows us to define an orthonormal basis of operators with well-defined sizes, over which any operator can be uniquely decomposed.

We now show that the dynamics preserves the sizes of operators on average, in the sense that it does not introduce components with operator sizes different from those present at $t=0$. We briefly present the main reasoning here, whilst a more detailed description of the following procedure can be found in App.~\ref{a:proof_op_size}.

Taking the disorder average of a simple annihilation operator $\hat{c}_i(t)$ in the Heisenberg picture gives
\begin{equation}\label{expand_exp}
    \eave{e^{-\mathcal{L} t}}\hat{c}_i = \sum_{n=0}^\infty \frac{(\mathcal{K}_q t)^{2n}}{(2n)!}\eave{J_{\alpha_1}\dots J_{\alpha_{2n}}}\ell_{\alpha_1} \dots\ell_{\alpha_{2n}}\hat{c}_i,
\end{equation}
where the superoperators $\ell_{\alpha_i}$ are defined at the end of Sec.~\ref{s:model}, and the sum over repeated multi-indices is implicit. After using Wick's theorem, the previous equation results in a sum of nested commutators involving the operators $\hat{h}_{\alpha_1},\dots, \hat{h}_{\alpha_n}$, as well as their Hermitian conjugates $\hat{h}^\dag_{\alpha_1},\dots, \hat{h}^\dag_{\alpha_n}$, appearing in a certain order. In particular, each multi-index $\alpha_k$ is shared precisely by $\hat{h}_{\alpha_k}$ and $\hat{h}^\dag_{\alpha_k}$, as follows from Eqs.~\eqref{disorder_properties-b} and~\eqref{disorder_properties-c}. If we expand each nested commutator in sequences of fermionic operators, for each $\hat{c}^\dag_k$ appearing in any string there will always be exactly one $\hat{c}_k$, with the same index, somewhere within that same string.
For each chain of fermionic operators, we may use the anticommutation relations to move $\hat{c}_i$ towards one of the edges, for example the left one. Each new sub-string originated by this process involves only the operator $\hat{c}_i$ and multiple sums over all sites of paired creation and annihilation operators. Following this argument, any initial string is equal to the original operator $\hat{c}_i$ multiplied by a complex combination of fermionic operators that, however, cannot depend on any lattice site in particular, as all indices involved are summed over: As a consequence, the latter will be a function of the number of sites $N$ and the charge operator $\hat{Q}$. In conclusion, we have proved that
\begin{equation}\label{absence_growth}
    \eave{e^{-\mathcal{L} t}}\hat{c}_i = \hat{c}_i \hat{f}(t^2;N,\hat{Q}),
\end{equation}
where the dependence on a function of $t^2$ arises from the absence of odd terms in Eq.~\eqref{expand_exp}. The above result shows that, if we limit its action to a subset of the Hilbert space with fixed $Q$, $\hat{c}_i$ is an eigen-operator of the generator of the ensemble-averaged dynamics. In other words, $\hat{c}_i(t)$ manifests no operator growth (nor hopping) on average.

Our above argument is readily extended to all fermionic strings of type $\hat{S}_{m,n}=\hat{c}^\dag_{i_1}\dots\hat{c}^\dag_{i_m}\hat{c}_{j_1}\dots\hat{c}_{j_n}$ in which all indices are distinct, i.e.,
\begin{equation}\label{absence_growth_general}
    \eave{e^{-\mathcal{L} t}}\hat{S}_{m,n} = \hat{S}_{m,n} \hat{f}_{m,n}(t^2;N,\hat{Q}).
\end{equation}
These results apply also to strings containing repeated indices, i.e., number operators, and thus they hold for any generic operator with fixed size $(m,n)$: The proof is presented in App.~\ref{a:proof_op_size}. Notice that the functions $\hat{f}_{m,n}$ are independent of the specific lattice sites appearing in $\hat{S}_{m,n}$, consistently with the $S_N$ symmetry of the SYK model under ensemble average.

In summary, we know how all operators with a well-defined size evolve on average: There is no growth, and the time-dependence of the evolution depends only on the operator size. We point out that the derivation of this conclusion is model-independent, up to the requirement that the couplings in the Hamiltonian bear no dependence on lattice indices. For this reason, not only does this proof apply to any SYK$_q$ model, but also to linear combinations of SYK Hamiltonians with different values of $q$, as well as to other fully-connected homogeneous models. In addition, we expect the result to extend also to fermionic and spin SYK-like models with non-Gaussian disorder, as well as a bosonic versions of the system with Gaussian randomness~\footnote{Gaussian disorder implies that the superoperators $\ell_\alpha$ are contracted in pairs. For fermionic systems, this is not fundamental, because the square of a creation or annihilation operator vanishes: Even if multiple operators with the same lattice index appear within a fermionic string, eventually the expression simplifies leaving at most one of them. A similar argument applies to spin models. In contrast, more-than-pairwise contractions generate new distinct operators for bosons, and the proof of operator size symmetry does not hold.}.

\subsection{Implications for correlation functions and OTOCs}\label{ss:single_realizations}

The previous results might misleadingly be interpreted as total absence of operator growth in the SYK model, which would be an incorrect conclusion. We stress that operator size symmetry holds only under disorder average, whereas individual realizations of the system do clearly manifest present operator growth (see Sec.~\ref{ss:single_realizations} for a detailed discussion). Indeed, operator growth is related to scrambling of quantum information, which is one key characteristic of the SYK model.

While it is important to distinguish physical features computed in single realizations and on average, there are cases in which the two coincide, namely for self-averaging observables. At suitably large $N$, such quantities approach their ensemble-averaged values for any single instance of the model, and thus operator size symmetry may manifest also for individual disorder realizations. Typical self-averaging quantities are expectation values of operators with extensive support~\cite{castellani2005,torresherrera2020}, which are precisely in the form of Eq.~\eqref{observable_dyn}. As a consequence, if, for instance, $\hat{W}=\hat{S}_{m,m}$ has a fixed operator size $(m,m)$, and $\hat{\rho}(0)$ has fixed number of particles $Q$, we have
\begin{equation}
    \braket{\hat{W}(t)} \approx \eave{\braket{\hat{W}(t)}} = f_{m,m}(t^2;N,Q)\Tr\left(\hat{W} \hat{\rho}_0\right),
\end{equation}
for large $N$. This implies that the functional form of the time-evolution is independent of the choice of the initial state, and the latter sets only the amplitude of the curve, as was observed in our previous work~\cite{bandyopadhyay2021}.

As anticipated previously, operator size conservation on average does not imply absence of scrambling in single physical realizations. An important aspect of the previous discussion is the limitation of our framework to the study of quantities that are linear in the time-evolution superoperator. In particular, OTOCs of the form
\begin{equation}
    \eave{\Tr\left(\comm{\hat{W}(t)}{\hat{V}(0)}^2\right)},
\end{equation}
which diagnose quantum chaotic behavior, are beyond the scope of our investigation, because they involve correlations between two superoperators $e^{-\mathcal{L} t}$. We conclude that quantum information scrambling must be encoded in these correlations, which arise because operators themselves are not self-averaging, namely $\hat{W}(t)\neq\eave{\hat{W}(t)}$ even if their quantum-mechanical expectation values become equal at large $N$. The same argument applies to the so-called Krylov complexity~\cite{parker2019,rabinovici2022,jian2021}, which is used to quantify operator growth. This is consistent with the belief that simple equal-time correlators are unable to capture the chaotic properties of the model, e.g., the quantum Lyapunov exponent. Such features are manifested only by more complicated quantities, such as, for instance, OTOCs.

To conclude this section, we point out that multiple studies in the literature focus on operator growth from the perspective of a distribution for the operator size~\cite{roberts2018,lensky2020,qi2019,lucas2020}. In such studies, operators are expanded over a basis of operators with fixed sizes, and a size distribution is obtained from the squared moduli of the expansion coefficients. Eventually, this framework allows one to compute the average size (see also~\cite{lin2019,jian2021}) and its dynamics, which can be directly related to OTOCs. Nevertheless, this bears no contradiction with our findings. The distribution has quadratic dependence on the time-dependent expansion coefficients, and thus its disorder-average probes correlations between two time-evolution superoperators; as discussed previously, our results do not apply to such a scenario.

\section{Cumulant Expansion Method}\label{s:cumulant_expansion}

In the previous section, we discussed how operators with well-defined size do not grow on ensemble average. Nevertheless, there is no straightforward way to explicitly determine the functions $\hat{f}_{m,n}(t^2;N,\hat{Q})$: The model is characterized by a single energy scale, and thus we cannot set up conventional perturbation theory starting from a known non-interacting solution. Instead, we now develop a cumulant expansion scheme to approximate the disorder-averaged dynamics of equal-time observables.

The numerical investigation presented in Ref.~\cite{bandyopadhyay2021} shows that the quench dynamics of an operator can be well approximated by a Gaussian, even though the agreement is not perfect. Motivated by this result, we look for an exponential representation of the disorder-averaged time-evolution superoperator, in such a way as to recover the Gaussian shape as the lowest order approximation, and to implement additional corrections. We formally write
\begin{equation}\label{cumulant_series}
    \eave{e^{-\mathcal{L} t}} = e^{\mathcal{C}(t)} = \exp\left(\sum_{k=1}^\infty \frac{(Jt)^k}{k!}\mathcal{C}_k\right),
\end{equation}
where $\mathcal{C}(t)$ will be referred to as a cumulant generating superoperator, in agreement with usual nomenclature in statistics. Accordingly, the superoperators $\mathcal{C}_k$ will be referred to as cumulants. Since we are considering Gaussian disorder, if $\mathcal{L}$ was just a scalar function instead of a superoperator, then all cumulants with $k>2$ would vanish. In the present case, however, the infinite series on the right-hand side of the previous equation does not terminate at finite order, because the operators that multiply different disordered couplings do not commute.

A convenient way to determine the cumulants is by expanding the right-hand side of Eq.~\eqref{cumulant_series}, and comparing it to Eq.~\eqref{expand_exp} by matching equal powers of $t$. We immediately see that all superoperators $\mathcal{C}_k$ with odd $k$ vanish. The identification of $\mathcal{C}_k$ requires the previous determination of all $\mathcal{C}_l$ with $l<k$, and thus the procedure is iterative. Here, we present the first two non-vanishing terms:
\begin{subequations}
\begin{equation}\label{c2}
    \mathcal{C}_2 = \mathcal{K}_q^2 \sum_\alpha \ell_\alpha^2,
\end{equation}
\begin{equation}\label{c4}
     \mathcal{C}_4 = \mathcal{K}_q^4 \sum_{\alpha,\beta}\left(\ell_\alpha\ell_\beta^2\ell_\alpha + \ell_\alpha\ell_\beta\ell_\alpha\ell_\beta-2\ell_\alpha^2\ell_\beta^2\right).
\end{equation}
\end{subequations}
Leading order truncation of the cumulant expansion yields precisely the time-evolution superoperator used in Ref.~\cite{bandyopadhyay2021}, which was derived alternatively by writing an effective master equation for the ensemble-averaged density matrix. Higher-order cumulants provide corrections to the results of this previous work, allowing to achieve a better approximation of the dynamics.

For the SYK model, all cumulant superoperators commute. This property is proved by observing that each superoperator of type $\ell_{\alpha_1}\dots\ell_{\alpha_{2k}}$, in which the indices are contracted in pairs and summed over, preserves operators with well-defined size, exactly as the full $\eave{e^{-\mathcal{L} t}}$ in Eq.~\eqref{absence_growth_general} (see App.~\ref{a:proof_op_size} for the details). As a consequence, we can differentiate Eq.~\eqref{cumulant_series} to obtain an effective master equation for the disorder-averaged density matrix in the Schrödinger picture~\footnote{Notice that $e^{\mathcal{C}(t)} = \eave{e^{-\mathcal{L} t}} = \eave{e^{\mathcal{L} t}}$ because all odd powers in the series expansion vanish.}, yielding
\begin{equation}\label{me}
    \dv{}{t}\eave{\hat{\rho}(t)} = \dv{}{t}e^{\mathcal{C}(t)}\hat{\rho}_0= \left(J\sum_{k=1}^\infty \frac{(Jt)^{2k-1}}{(2k-1)!}\mathcal{C}_{2k}\right)\eave{\hat{\rho}(t)}.
\end{equation}
The sum within the round brackets is the effective Liouvillian that generates the disorder-averaged dynamics. Equation~\eqref{me} is a generalization of the master equation proposed in Ref.~\cite{bandyopadhyay2021}, and it is formally exact.

For practical purposes, we can only evaluate the cumulant expansion up to some given finite order. Even if all $\mathcal{C}_k$ were known, finding an analytic expression for $\mathcal{C}(t)$ remains as hard as computing the exact disorder-averaged time-evolution superoperator; as a consequence, the series must be truncated. We point out, however, that even if only a finite number of terms is known, the resulting approximate description of the dynamics can be quite accurate over a wide time window. Indeed, using a truncated cumulant expansion scheme to characterize the average dynamics already proved successful in systems with non-Markovian noise~\cite{groszkowski2022}. Evaluating the series up to any given order is only as expensive as computing a short-time expansion up to that same order. Nevertheless, while the latter is limited to early times only, the former is potentially able to reasonably reproduce the evolution at arbitrary times. The reason lies in the exponential form of the time-evolution superoperator: Even when $\mathcal{C}(t)$ is truncated, Eq.~\eqref{cumulant_series} still involves an infinite series of powers of $t$, and thus it can represent non-polynomial time-dependence. In addition, we can argue that the impact of higher-order cumulants is relevant only on longer timescales, as they are suppressed by a factor $(2k)!$~\footnote{The timescale at which a cumulant becomes relevant is long if its magnitude in Eq.~\eqref{expand_exp} is small. Any $\mathcal{C}_{2k}$ with $k>1$ should vanish if all superoperators $\ell_\alpha$ commuted, and thus each cumulant must contain a certain number $2M_{2k}$ of superoperators of type $\ell_{\alpha_1}\dots\ell_{\alpha_{2k}}$, half of them with the positive sign, half with the negative one (for example, see Eq.~\eqref{c4}): We do not expect $\mathcal{C}_{2k}$ to grow as fast as the factor $(2k)!$, appearing at the denominator. On top of that, $M_{2k}$ is the number of irreducible Wick contractions arising from the $2k$th term of Eq.~\eqref{expand_exp}, which is smaller than $(2k-1)!!$ (it being the total number of contractions of order $2k$). By irreducible, we refer to those contractions that cannot be written as products of lower-order ones. For example, referring to Eq.~\eqref{c4}, the superoperators appearing with the plus sign are irreducible, whereas the other one is reducible. In conclusion, not only $2M_k<2(2k-1)!!<(2k)!$, but also different superoperators sum up in a destructive way.}. Appendix~\ref{a:exact_evals} presents a numerical check of this statement, which indeed supports its validity. Suppose that the disorder-averaged density matrix quickly approaches a steady state, so that $\eave{\hat{\rho}(t)}\approx\eave{\hat{\rho}(\infty)}$ for $t \gtrsim\tau$. In this case, neglecting high-order cumulants that become relevant only at times greater than $\tau$ is a good approximation, because, as seen from Eq.~\eqref{me}, their action on the steady state density matrix is practically zero (notice that all cumulants share the same steady states, because they commute). In conclusion, for a system that thermalizes quickly, a finite number of cumulants is sufficient to obtain a valid approximation of the exact dynamics. As shown in Ref.~\cite{bandyopadhyay2021}, the quench dynamics of observables in the SYK model fall within this situation, as they manifest super-exponential relaxation to stationary values.

The above cumulant expansion allows us to approximate the functions $\hat{f}_{m,n}(t^2;N,\hat{Q})$. First, as observed previously, an equation analogous to Eq.~\eqref{absence_growth_general} can be written for each individual cumulant, so that
\begin{equation}
    \mathcal{C}_{2k} \hat{S}_{m,n} = \hat{S}_{m,n}\hat{\lambda}^{(2k)}_{m,n}(N,\hat{Q}),
\end{equation}
where $\hat{S}_{m,n}$ is a generic operator in the form of Eq.~\eqref{size_operators}. The dynamical functions $\hat{f}_{m,n}(t^2;N,\hat{Q})$ are then expanded as
\begin{equation}\label{expand_fmn}
    \hat{f}_{m,n}(t^2;N,\hat{Q}) = \exp\left(\sum_{k=1}^\infty \frac{(J t)^{2k}}{(2k)!}\hat{\lambda}^{(2k)}_{m,n}(N,\hat{Q})\right).
\end{equation}
As argued before, a finite number of functions $\hat{\lambda}^{(2k)}_{m,n}$ can be sufficient to obtain a very good approximation of the exact result. The advantage of this approach is that eigen-operators of cumulants are known, and thus we can evaluate their eigenvalues analytically. Some exact expressions for the SYK$_4$ model are presented in App.~\ref{a:exact_evals}. From a practical point of view, this greatly simplifies the numerical application of our formalism: We can approximate Eq.~\eqref{expand_fmn} for arbitrarily large system sizes without resorting to exact diagonalization of cumulant superoperators.

\section{Numerical results}\label{s:application}

This Section presents numerical results that complement the previous theoretical findings. For a given operator, we first investigate the dynamics of those components  which have a well-defined size, and then we compare the exact time-evolution to the approximate one produced by the cumulant expansion. The exact dynamics is obtained by exact diagonalization of the SYK$_4$ Hamiltonian within a sector of fixed total charge $Q$. We evolve an initial disorder-independent pure state $\ket{\psi_0}$, and then we evaluate the expectation value $\braket{\hat{W}(t)}$ of an observable $\hat{W}$ that commutes with the charge $\hat{Q}$. In our numerics, we choose $\ket{\psi_0}$ to be the Néel state $\ket{1010\dots}$, where the ordering of lattice sites is arbitrary. Regarding the choice of the operator, we study the staggered magnetization $\hat{R} = \sum_{i=1}^N(-1)^{i+1}\hat{n}_i$ (the indices are assigned with the same convention used for the Néel state). Disorder-averaged quantities are computed by iterating this procedure multiple times for independent realizations of the disordered Hamiltonian given in Eq.~\eqref{hamiltonian-a}, and finally taking a statistical average over the sample.

\subsection{Operator size throughout the dynamics}\label{ss:operator_growth}

We now proceed to study the dynamics of individual components with well-defined size of a given operator. As mentioned previously, it is possible to define an orthonormal (with respect to the Hilbert--Schmidt product) basis of operators such that each basis element has a well-defined operator size. Using the notation introduced in App.~\ref{a:op_size}, we denote these basis operators as $\hat{T}^{(k)}_{m_k,n_k}$, and each of them has size $(m_k,n_k)$. We can extract their amplitudes in the expansion of $\hat{W}(t)$ by taking the Hilbert--Schmidt product
\begin{equation}\label{c_mn}
    c^{(k)}_{m_k,n_k}(t) = \Tr_Q\left(\left(\hat{T}^{(k)}_{m_k,n_k}\right)^\dag \hat{W}(t)\right),
\end{equation}
where $\Tr_Q$ denotes the trace over states with fixed charge $Q$, as defined in App.~\ref{a:op_size}. For $\hat{W}=\hat{R}$, we use exact diagonalization to study the dynamics of some coefficients $c^{(k)}_{m_k,m_k}(t)$, representative of others. We focus on $m=1,2$, and for each operator size we design two operators, one diagonal and the other off-diagonal in the basis of Fock states. Following the prescription of Eq.~\eqref{size_operators}, we introduce
\begin{subequations}\label{test_operators}
\begin{align}
\label{test_operators-a}
& \hat{T}_{1,1}^{(\text{diag})} = \mathcal{N}^{(\text{diag})}_{1,1}\left(\hat{n}_1-\frac{1}{2}\right),
\\
\label{test_operators-b}
& \hat{T}_{2,2}^{(\text{diag})} = \mathcal{N}^{(\text{diag})}_{2,2}\left(\hat{n}_1\hat{n}_2-\frac{\hat{n}_1+\hat{n}_2}{2} + \frac{N}{4(N-1)}\right),
\\
\label{test_operators-c}
& \hat{T}_{1,1}^{(\text{off-diag})} = \mathcal{N}^{(\text{off-diag})}_{1,1}\left(\hat{c}^\dag_1\hat{c}_2+\mathrm{h.c.}\right),
\\
\label{test_operators-c}
& \hat{T}_{2,2}^{(\text{off-diag})} = \mathcal{N}^{(\text{off-diag})}_{2,2}\left(\hat{c}^\dag_1\hat{c}^\dag_2\hat{c}_3\hat{c}_4+\mathrm{h.c.}\right),
\end{align}
\end{subequations}
where the prefactors $\mathcal{N}^{(k)}_{m_k,m_k}$ are such that these operators have unitary Hilbert--Schmidt norm. Figure~\ref{f:operator_growth} presents the disorder-averaged dynamics of the coefficients $c^{(k)}_{m_k,m_k}(t)$, as well as the evolution for some individual disorder realizations. As expected from the exact theoretical discussion in Sec.~\ref{ss:size_sym}, the operator $\hat{R}$ has no overlap with the size $(2,2)$. The coefficient $c^{(\text{off-diag})}_{1,1}(t)$ is also vanishing because the staggered magnetization only involves number operators $\hat{n}_i$. In contrast, $c^{(\text{diag})}_{1,1}(t)$ manifests a quick decay to zero, mirroring the super-exponential time-dependence of the expectation value $\braket{\hat{R}(t)}$ itself. This is substantiated by the inset in Fig.~\ref{f:operator_growth}a, which shows that the time-dependence of $\eave{c_{1,1}^\text{diag}(t)}$ coincides with that of $\eave{\braket{\hat{R}(t)}}$ apart from fluctuations that can be attributed to the statistical error of the averages. If we consider the dynamics of a single realization of the SYK model, in general $\hat{R}(t)$ acquires components of size different from $(1,1)$, and it does not remain diagonal in the basis of Fock states. Still, Fig.~\ref{f:operator_growth}a shows that the typical time-evolution is close to the disorder-averaged one. This effect is a consequence of self-averaging: In fact, we observe larger deviations from the average curve as we decrease the system size.
\begin{figure}[ht]
    \includegraphics[width =
    \linewidth]{operator_growth.pdf}
    \caption{Dynamics of the coefficients $c^{(k)}_{m_k,m_k}(t)$ computed in single disorder realizations (red curves) and in ensemble average (blue curves) for the operator $\hat{R}$. The inset of panel (a) portrays the difference $\Delta(t) = \left|\eave{ c_{1,1}^\text{diag}(t)/c_{1,1}^\text{diag}(0)} - \eave{\braket{\hat{R}(t)}/\braket{\hat{R}(0)}}\right|$. Data is for $N=10$ at half filling $Q=N/2$. The ensemble averages are estimated by performing statistical averages over $10000$ independent disorder realizations. The terms of size $(2,2)$ vanish on average, in agreement with the predicted absence of operator growth. In contrast, there is a finite overlap with the diagonal operator of size $(1,1)$, and the time-evolution matches the full dynamics of $\braket{\hat{R}(t)}$, as can be appreciated by the inset of panel (a). Curves for individual realizations show that components with sizes different from $(1,1)$ are in general developed throughout the evolution.}
    \label{f:operator_growth}
\end{figure}

\subsection{Application of the cumulant expansion}\label{ss:cumulant_exp}

We now test the performance of the cumulant expansion by comparing it to the dynamics obtained through numerically exact simulations. As discussed in Sec.~\ref{s:operator_size}, each constituent of an operator $\hat{W}$ with well-defined size will evolve differently. The disorder-averaged dynamics thus reads
\begin{equation}\label{decompose_dynamics}
    \eave{\braket{\hat{W}(t)}} = \sum_m c_{m,m} f_{m,m}(t^2;N,Q),
\end{equation}
where
\begin{equation}
    c_{m,m} = \Tr_Q\left(\hat{\rho}_0\mathcal{P}_{m,m}\hat{W}\right),
\end{equation}
and $\mathcal{P}_{m,m}$ is a projection superoperator over the subspace with operator size $(m,m)$. Notice that no components of size $(m,n)$ with $m\neq n$ can appear because $\hat{W}$ preserves the total charge. We explicitly see that the only role of the initial state is to set the amplitudes $c_{m,m}$. In our previous work~\cite{bandyopadhyay2021}, we provide numerical evidence that, for some observables, the time-dependence of $\eave{\braket{\hat{W}(t)}}$ manifests universality with respect to the choice of the initial state. This is readily explained by Eq.~\eqref{decompose_dynamics}: Operators that overlap with a unique operator size evolve according to a single dynamical function, and the initial state affects only the amplitude. In contrast, if two or more coefficients (excluding $c_{0,0}$) are non-zero, then changing $\hat{\rho}_0$ modifies the relative weights of different dynamical functions, and no universality is found. To summarize, the shape of the dynamics of an operator with well-defined size is independent of both the initial state and the precise definition of the operator itself.

\begin{figure*}[ht]
    \includegraphics[width =
    \linewidth]{ed_vs_ce_dynamics.pdf}
    \caption{ Comparison of disorder-averaged dynamics of $\hat{R}$ and $\hat{R}^2$ obtained using exact diagonalization (solid blue curve), and the cumulant expansion truncated at second (dashed orange curve), fourth (dashed red curve), and sixth (dashed purple curve) order. Data is for $(\mathrm{a})$ and $(\mathrm{d})$ $N=8$ with 100000 samples, $(\mathrm{b})$ and $(\mathrm{e})$ $N=12$ with 3000 samples, $(\mathrm{c})$ and $(\mathrm{f})$ $N=16$ with 500 samples, at half filling $Q=N/2$.
    For $\hat{R}^2$, in all cases the addition of an extra cumulant manifestly improves the approximation of the exact dynamics (see insets). The same is true in the case of $\hat{R}$ when adding $\mathcal{C}_4$ for $N=8,12$, whereas long-time divergences appear when adding $\mathcal{C}_6$ for all considered system sizes, and already when adding $\mathcal{C}_4$ for $N=16$. These pathologies are due to an inappropriate truncation of the cumulant at some system sizes. We also point out that generally (in the time windows preceding divergencies, if present) the variation induced by introducing the fourth cumulant is larger than that due to the sixth cumulant, which suggests quick convergence of the cumulant series.
    }
    \label{f:numerics}
\end{figure*}

We study the dynamics of the staggered magnetization $\hat{R}$ and of its square $\hat{R}^2$. The former has operator size $(1,1)$, while the latter contains both sizes $(0,0)$ and $(2,2)$. Figure~\ref{f:numerics} presents the disorder-averaged dynamics of the operators for $N=8,12,16$, comparing it to the results of the cumulant expansion method truncated at different orders. Focusing on the operator $\hat{R}^2$, for all considered system sizes, the quality of the approximation is found to improve with the number of cumulants considered. In addition, the plots suggest that the curves converge rapidly to the exact result as more cumulants are included; even though this claim cannot be rigorously confirmed with only three cumulant superoperators, it agrees with the discussion of Sec.~\ref{s:cumulant_expansion}, where we argued that high-order corrections provide negligible contributions as compared to low-order ones. These conclusions hold also for the operator $\hat{R}$, with the clear exception of some curves that manifest unphysical divergencies at late times. This happens because a cumulant eigenvalue $\lambda_{m,n}^{(2k)}$ can have a positive sign for some values of $N$, which makes the truncation of the cumulant series to that specific order insufficient. We point out that also $\lambda_{2,2}^{(4)}$ and $\lambda_{2,2}^{(6)}$, which characterize the dynamics of $\hat{R}^2$, acquire negative values for larger system sizes than those considered here. Any divergence must be compensated by the presence of a higher-order cumulant with a negative eigenvalue that restores the relaxation to the steady state.

As previously exposed in Ref.~\cite{bandyopadhyay2021}, the steady state of the ensemble-averaged dynamics corresponds to the infinite temperature state. Some deviations are observed for small system sizes, but these can be attributed to finite-size effects. This result is reproduced by the cumulant expansion approach, as the identity is an eigenoperator of all cumulants, and is unique for all cases which we studied.

\section{Discussion and conclusions}\label{s:conclusions}
In this work, we have studied the out-of-equilibrium properties of the SYK model in terms of the spectral features of the effective time-evolution superoperator for the disorder-averaged dynamics. Owing to the absence of any spatial ordering, it is possible to identify the exact eigenoperators of $\eave{e^{-\mathcal{L} t}}$, and to prove that its (time-dependent) spectrum  manifests high degeneracies, corresponding to different operator sizes. As a result, for all system sizes, the dynamics is characterized by a non-trivial symmetry, which hides the scrambling properties of the model when probing the evolution of standard expectation values of observables in the ensemble average. This conclusion is limited to correlation functions that involve a single time-evolution superoperator, and thus it does not apply to multiple-time functions nor, specifically, OTOCs. Our findings prove rigorously that the dynamics of observables with well-defined operator size manifest universal features with respect to the initial state, as highlighted by our previous numerical study~\cite{bandyopadhyay2021}.

The proof of operator size symmetry implies the existence of dynamical functions, obtained from the diagonalization of $\eave{e^{-\mathcal{L} t}}$, that completely characterize the dynamics. Nevertheless, the demonstration does not indicate a simple way to compute them explicitly, and thus we developed a cumulant expansion scheme to achieve an approximation thereof. By construction, cumulant superoperators inherit the spectral structure of the exact time-evolution generator, and thus their eigenvalues can be determined analytically, as functions of the system size $N$ and total charge $Q$. A direct comparison between this approach and exact diagonalization simulations reveals that the method is successful in reproducing the evolution to good accuracy in all time regimes. In addition, numerical results corroborate theoretical arguments that the cumulant expansion converges quickly with the number of terms included. Although this approximation scheme appears to be very effective for the system considered here and, potentially, variants thereof, we expect its usefulness to be limited for generic models: In the absence of operator size symmetry, it is impossible to determine the eigenoperators of the cumulants analytically, which implies one would need to resort to computationally expensive explicit numerical diagonalization. This difficulty is overcome for the SYK model, thus enabling the cumulant expansion method to effectively compute the dynamical functions.

We believe our investigation enriches the current understanding of the SYK model by showcasing the presence of an unexpected symmetry, and by illustrating that scrambling can be completely absent in some physical quantities even if the model is quantum chaotic. We expect our results to be relevant for future experimental implementations of the SYK model, as regular equal-time expectation values of observables are typical measurable quantities, further highlighting the need for measurement schemes that access more intricate observables~\cite{torlai2020}. In addition, the cumulant expansion approach that we developed manages to address arbitrary system sizes, thus enabling to access a regime that is out of reach of practical numerical simulations. We stress that a large enough system is fundamental to exploit the self-averaging property, which allows one to compare individual measurements to disorder-averaged quantities.

A natural topic of future investigation is whether it is possible to generalize our open system framework to study OTOCs, and thus to potentially reveal generic features also for such quantities. Moreover, as mentioned previously, the present work bears implications for other systems: We believe that similar features to the ones discussed here are also manifested by other fully-connected disordered models. In particular, the proof of operator size symmetry does not rely on the details of the SYK model, and can be extended to fermionic and spin systems with other disorder distributions, as well as to bosonic models with Gaussian disorder.

\emph{Acknowledgements.---}We thank Andrea Legramandi for reading the manuscript and giving useful comments. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 804305), the Bundesministerium für Wirtschaft und Energie through the project ``EnerQuant” (project ID 03EI1025C), Provincia Autonoma di Trento, and by Q@TN, the joint lab between University of Trento, FBK-Fondazione Bruno Kessler, INFN-National Institute for Nuclear Physics and CNR-National Research Council.
S.B. acknowledges CINECA for the use of HPC resources under ISCRA-C project ISSYK-2 (HP10CP8XXF).