## Introduction

Co-citation, “the frequency with which two documents from the earlier literature are cited together in the later literature”, was first described in 1973 [marshakova-shaikevich_system_1973, small_co_citation_1973]. As noted by [small_co_citation_1973], co-citation patterns differ from bibliographic coupling patterns [kessler_bibliographic_1963] but align with patterns of direct citation and frequently co-cited publications must have high individual citations.

Co-citation has been the subject of further study and characterization, for example, comparisons to bibliographic coupling and direct citation [boyack_cocitation_2010], the study of invisible colleges [gmur_co-citation_2003, noma_co-citation_1984], construction of networks by co-citation [small_clustering_1985, small_clustering_1985-1], evaluation of clusters in combination with textual analysis [braam_mapping_1991], textual similarity at the article and other levels [colavizza_closer_2018], and the fractal nature of publications aggregated by co-citations [vanraan_fractal_1990].

Co-citations provide details of the relationship between key (highly cited) ideas, and changes in co-citation patterns over time may provide insight into the mechanism with which new schools of thought develop. Implicit in the definition of co-citation is novel combinations of existing ideas, but only some frequently co-cited article pairs reflect surprising combinations. For example, two publications presenting the leading methods for the same computational problem may be highly co-cited, but this does not reflect a novel combination of ideas. Similarly, two publications describing methods that often constitute part of the same workflow may be highly co-cited, but these co-citations are also not surprising. On the other hand, for two articles in different fields, frequent co-citation is generally unexpected.

Novel, atypical, or otherwise unusual combinations of co-cited articles have been explored at the journal-level [wang_2017, bradley_co-citations_2019, boyack_vs_uzzi_2014, uzzi_atypical_2013]. However, journal-level classifications have limited resolution relative to article-level studies, which may better represent the actual structure and aggregations of the scientific literature [shu_comparing_2019, article_boyack_topic, waltman_new_2012, article_stasa, article_gomez_journal]. Accordingly, we sought to discover measurable characteristics of frequently co-cited publications from an article-level perspective.

To study frequently co-cited articles, we have developed a novel graph-theoretic approach that reflects the citation neighborhood of a given pair of articles. In seeking to determine the degree to which a co-cited pair of papers represented a surprising combination, we wished to avoid journal-based field classifications, which present challenges. Instead, we attempted to use citation history to produce an estimate of the probability that a given pair of publications would be co-cited. Since we focus on the activity before they are first co-cited, the “probability" of co-citation is zero, by definition, since there are no co-citations yet. Hence, we approximated co-citation probabilities: we treat an article that cites one member of a co-cited pair and also cites at least one article that cites the other member as a proxy for co-citation. Specifically, given a pair of publications , we construct a directed bipartite graph whose vertex set contains all publications that cite either or previous to their first co-citation. We then compute , a normalized count of such proxies, and use it to predict the probability of co-citation between and . This approach enables an evaluation that is specific to the given pair of articles, and does so without substantial computational cost, while avoiding definitions of disciplines derived from journals or having to measure disciplinary distances.

To support our analysis, we constructed a dataset of articles from Scopus [scopus_ref] that were published in the eleven year period, 1985-1995, and extracted the cited references in these articles. Recognizing that frequently co-cited publications must derive from highly-cited publications [small_co_citation_1973], we identified those reference pairs (33.6 million pairs) for each article in the dataset that are drawn from the top 1% most cited articles in Scopus and measured their frequency of co-citation.

To investigate which statistical distributions might best describe the co-citation frequencies in these 33.6 million co-cited pairs, we reviewed prior work on distributions of citation frequency [radicchi_statistical_2008, eom_2011, price_1965, price_general_1976, newman_structure_2003, wang_quantifying_2013, stringer_statistical_2008, stringer_statistical_2010, redner_statistical_2005]. This research has fit the frequency distribution of citation strength sometimes to a power law distribution and other times to a lognormal distribution. A graph of the analogous co-citation data suggests that power law or lognormal distributions are candidates for describing co-citation strength as well and so we, accordingly, investigated that conjecture. Interestingly, [Mitzenmacher_2003] notes the debate between the appropriateness of power law versus lognormal distributions is not confined to bibliometrics, but has been at issue in many disciplines and contexts.

To study how the best-fit distributional function and parameters for co-citation might vary with , we stratified co-citation frequency data. We also measured whether a direct link exists between two members of a co-cited pair (i.e., whether one member of a pair cites the other) and how this property is related to co-citation frequencies. We find that the distribution of co-citation frequencies varies with and that a power law distribution fits co-citation frequencies more often when is small, whereas a lognormal distribution fits more often for large .

A pertinent aspect of co-citation is the rate at which frequencies accumulate. While citation dynamics of individual publications have been fairly well studied by others, for example, [wallace_2009, eom_2011], the dynamics of co-cited articles are less well studied. Our interest was the special case analogous to the Sleeping Beauty phenomenon [vanraansleeping2004, ke_defining_2015], which may reflect delayed recognition of scientific discovery and the causes attributed to it [merton_1963, garfield_1970, garfield_1980, cole_1970, barber_1961, glanzel_myth_2004]. Thus, we also identified co-cited pairs that featured a period of dormancy before accumulating co-citations.

## Materials & Methods

*Data* Citation counts were computed for all Scopus articles (88,639,980 records) updated through December 2019, as implemented in the
ERNIE project [GithubERNIE2019]

. Records with corrupted or missing publication years or classified as ‘dummy’ by the vendor were then removed, resulting in a dataset of 76,572,284 publications. Hazen percentiles of citation counts, grouped by year of publication, were calculated for the these data

[bornmann_use_2013]. The top 1% of highly cited publications from each year were combined into a set of highly cited publications consisting of 768,993 publications.Publications of type ‘article’, each containing at least five cited references and published in the 11 year period from 1985-1995, were subset from Scopus to form a dataset of 3,394,799 publications and 51,801,106 references (8,397,935 unique). For each of these publications, all possible reference pairs were generated and then restricted to those pairs where both members were in the set of highly cited publications (above).

For example, the data for 1985 consisted of 223,485 articles after processing as described above. Computing all reference pairs (that were also members of the highly cited publication set of 768,993) from these 223,485 articles gave rise to 2,600,101 reference pairs (Table 1) that ranged in co-citation frequency from 1 to 874 within the 1985 dataset; from 1 to 11,949 across the 11 year period 1985-1995; and from 1 to 35,755 across all of Scopus. Collectively, the publications in our 1985-1995 dataset generated 33,641,395 unique co-citation pairs, for which we computed co-citation frequencies across all of Scopus.

Year | Articles | References | Co-cited Pairs |
---|---|---|---|

1985 | 223,485 | 1,796,502 | 2,600,101 |

1986 | 238,096 | 1,920,225 | 2,840,557 |

1987 | 250,575 | 2,037,654 | 3,180,261 |

1988 | 269,219 | 2,182,571 | 3,406,902 |

1989 | 285,873 | 2,303,481 | 3,793,986 |

1990 | 305,010 | 2,490,909 | 4,546,915 |

1991 | 325,782. | 2,662,005 | 5,039,334 |

1992 | 343,239. | 2,846,607 | 5,622,164 |

1993 | 360,916 | 3,006,374 | 6,121,147 |

1994 | 387,062. | 3,228,240 | 7,022,499 |

1995 | 405,503. | 3,432,228 | 7,626,684 |

*Derivation of *
We now show how we define our prior on the probability of and being co-cited, based on the citation graph restricted to publications that cite either or (but not both) up to the year of their
first co-citation.
Recall that we defined a proxy co-citation of and to be
an article that cites one member of the co-cited pair and also cites at least one article that cites the other member.
The idea behind this definition is that
we consider papers that cite as proxies for , and papers that cite as proxies for . Thus, if a paper cites both and (where is a proxy for ), then it is a proxy for a
co-citation of and . Similarly, if a paper cites both and (where is a proxy for ), it is also a proxy for a co-citation of and .
This motivates the graph-theoretic formulation, which we now formally present.

We fix the pair and we define to be the set of all publications that cite (but do not also cite ), and are published no later than the year of the first co-citation of and . We similarly define . We define a directed bipartite graph with vertex set . Note that if cites then , and similarly for the case where cites . Note also that since we have restricted and that . We now describe how the directed edge set is constructed. For any pair of articles where and , if cites then we include the directed edge in . Similarly, we include edge if cites . Finally, if a pair of articles both cite each other, then the graph has parallel edges. By construction, this graph is bipartite, which means that all the edges go between the two sets and (i.e., no edges exist between two vertices in , nor between two vertices in ).

Note that by the definition, every edge in arises because of a proxy co-citation, so that the number of proxy co-citations is the number of directed edges in . Consider the situation where a publication cites (so that ) and also cites in : this defines three directed edges from to nodes of . We count this as three proxy co-citations, not as one proxy co-citation. Similarly, if we have a publication that cites and also cites in , then there are four directed edges that go from to nodes in and we will count each of those directed edges as a different proxy co-citation.

Accordingly, letting denote the cardinality of a set , we note , i.e., the number of directed edges that go between and , is the number of proxy co-citations between and . If no parallel edges are permitted, the maximum number of possible proxy co-citations is . Under the assumption that both and each have at least one article, we define , our prior on the probability of and being co-cited, as follows:

Note that if parallel edges do not occur in the graph, then , but that otherwise the value can be greater than . Note also that if (i.e., if there are no proxy co-citations) and that if every possible proxy co-citation occurs.

To efficiently calculate , we used the following pipeline. We copied Scopus data from a relational schema in PostgreSQL into a citation graph from Scopus into the Neo4j 3.5 graph database using an automated Extract Transform Load (ETL) pipeline that combined Postgres CSV export and the Neo4j Bulk Import tool. The graph vertex set is all publications, each with a publication year attribute, and the edge set is all citations between the publications. A Cypher index was created on the publication year. We developed Cypher queries to calculate and tuned performance by splitting input publication pairs into small batches and processing them in parallel, using parallelization in Bash and GNU Parallel. Batch size, the number of parallel job slots, and other parameters were tuned for performance, with best results achieved on batch sizes varying from 20 to 100 pairs. The results of calculations were cross-checked using SQL calculations. In the small number of cases where computed to (above) it was set to 1 for the purpose of this study.

*Statistical Calculations* We denote the observed co-citation frequency data by the multi-set

where is the total number of pairs of articles and is the observed frequency of the pair of papers being co-cited. Note that this is in general a multi-set, as different pairs of articles can have the same co-citation frequency. Let be the number of times that appears in (equivalent, is the number of pairs of articles that are co-cited times), and let denote the total number of pairs of articles that are co-cited at least times. Then

(1) |

where is a parameter we use to analyze the distribution’s right tail starting at varying frequencies. We describe in this subsection (i) the statistical computations for fitting lognormal and power law distributions to right tails of the observed co-citation frequency distributions as defined by (1) for various and (ii) how we assessed the quality of those fits. Further, we performed such analyses for various slices of the data, stratifying by and other parameters, as is described in the Results section.

We used a discrete version of a lognormal distribution to represent integer co-citation frequencies, , following [stringer_statistical_2008] and [stringer_statistical_2010], while appropriately normalizing for our conditional assessment of the right tail commencing at :

(2) | |||||

where and

are the mean and standard deviation, respectively, of the underlying normal distribution. These probabilities can be computed with the cumulative normal distribution,

using the well-known error function.

We fit distributions to the co-citation frequency data for various extremities of the right tail, as parameterized by , using a maximum (log) likelihood estimator (MLE). We solved for the best-fit distributional parameters for the lognormal distribution, and , by modifying a multi-dimensional interval search algorithm from [press_statistical_2007] and following [stringer_statistical_2010]. A compiled version of this code using the C++ header file, amoeba.h, is available on our Github site [GithubERNIE2019].

We fit a discrete power law distribution to the data for various values of , which was normalized for our conditional observations of the right tail:

(3) |

where the Hurwitz zeta function,

is a generalization of the Riemann zeta function, , as is needed for analysis of the right tail.

We solved first-order conditions for the (log) MLE to find the best-fit distributional exponent ,

(4) |

as described in [clauset_power-law_2009] and [goldstein_problems_2004], where , are the observed co-citations with frequencies at least as great as and is the number of such co-citations. We solved (4) to find using a bisection algorithm.

We used the goodness of fit (

) and the Kolmogorov-Smirnov (K-S) tests to assess the null hypothesis that the distribution of the observed co-citation frequencies and the best-fit lognormal distribution are the same, and similarly for the best-fit power law distribution. We also computed the Kullback-Leibler Divergence (K-L) between the observed data and the best-fit distributions.

Both the and K-S tests employed the null hypothesis that the observed co-citation frequencies, for , were sampled from the best-fit lognormal or power law distributions, which we denote by for , while suppressing the parameters specific to each of the distributions.

The usual statistic was computed by, first, grouping each of the observed co-citation frequencies into bins, denoted by for , and then computing

where is the observed number of co-citations having frequencies associated with the -th bin,

and is the expected number of observations for frequencies in bin , if the null hypothesis was true, in a sample with size equal to the number of observed data points, :

If the null hypothesis was true, then we would expect and to be approximately equal, with deviations owing to variability due to sampling.

Constructing the bins requires only that for every . Test outcomes are sometimes sensitive to the minimum permitted, which we will denote by

, and so we tested with multiple thresholds, including 10, 20, 50, and 70. Furthermore, statistical tests are stochastic: these multiple tests permitted a reduction in the probability of erroneously rejecting or accepting the null hypothesis based on a single test. The distribution of observed co-citation frequencies was skewed right with a long tail, so that aggregating bins to satisfy

was most critical in the right tail. This motivated a bin construction algorithm that aggregated frequencies in reverse order, starting with the extreme right tail. Algorithm 1 requires a set of the unique observed co-citation frequencies, , which includes the elements of the multiset without repetition. While Algorithm 1 does not guarantee in general that all bins satisfy , that criterion was satisfied for the observed data.We implemented a K-S test using simulation to generate a sampling distribution to account for the discrete frequency observations [internetks_bootstrap]. We denote the cumulative distribution of observed co-citation frequencies by , and the best-fit cumulative distribution by . The K-S test involves testing the maximum absolute difference between the observed and theorized cumulative distributions,

where is the number of observations giving rise to , against the distribution of such differences between samples from the theorized distribution with the same number of observations, ,

where is the empirical distribution of sample of size (notation suppressed) drawn from

. We generated 100 such random variables

for each test. We reject the null hypothesis if is larger than substantially all of the , say all but 5%, for equivalence with a -value of 0.05. The number of samples drawn yields a -value with a resolution of 1%.We computed the K-L Divergence two ways due to its asymmetry:

Separate from the tests above, we tested whether the distribution of co-citation frequencies was independent of using a test, using the null hypothesis that the co-citation frequency distribution was independent of

. We initially created a contingency table on

and co-citation frequency using these bins for , , and logarithmic bins for frequency to accommodate the skewed distributions:We, subsequently, aggregated these bins to have an expected number of co-citations in each bin equal to or greater than 5 to account for a decreasing number of observations as and frequency increased by having just two intervals for frequency: .

*Kinetics of Co-citation*
We extended prior work on delayed recognition and the Sleeping Beauty phenomemon [ke_defining_2015, vanraansleeping2004, li_distinguishing_2016, glanzel_myth_2004] towards co-citation.
We have modified the beauty coefficient (B) of [ke_defining_2015] to address co-citations by:
(i) counting citations to a pair of publications (co-citations) rather than citations to individual papers,
(ii) setting (age zero) to the first year in which a pair of publications could be co-cited (i.e., the publication year of the more recently published member of a co-cited pair), and
(iii) setting to the number of co-citations occurring in year . Rather than calculate awakening time as in [ke_defining_2015], we opted to measure the simpler length of time between and the first year in which a co-citation was recorded; we label this measurement as the timelag , so that if a co-citation was recorded in .

## Results and Discussion

Our base dataset, described in Table 1, consists of the 33,641,395 co-cited reference pairs (33.6 million pairs) and their co-citation frequencies, gathered from Scopus during the 11-year period from 1985-1995 (Materials and Methods). A striking distribution of co-citation frequencies with a long right tail is observed with a minimum co-citation of 1, a median of 2, and a maximum co-citation frequency of 51,567 (Figure 2). Approximately 33.3 of 33.6 million pairs (99% of observations) have co-citation frequencies ranging from 1–67 and the remaining 1% have co-citation frequencies ranging from 68–51,567. Since the focus of our study was co-citations of frequently cited publications, we further restricted this dataset to those pairs with a co-citation frequency of at least 10, which resulted in a smaller dataset of 4,119,324 co-cited pairs (4.1 million pairs) with minimum co-citation frequency of 10, median of 18, and a maximum co-citation frequency of 51,567. In order to focus on co-citations derived from highly cited publications, was calculated for all pairs with a co-citation frequency of at least 10. We also note whether one article in a co-citation pair cites the other (connectedness).

Influenced by the use of linked co-citations for clustering [small_clustering_1985], we also examined the extent to which members of a co-cited pair were also found in other co-cited pairs. We found that 205,543 articles contributed to 4.12 million co-cited pairs. The highest frequency observed in our dataset, 51,567 co-citations, was for a pair of articles from the field of physical chemistry: Becke (1993) [becke_densityfunctional_1993] and Lee, Yang, and Parr (1988) [lee_development_1988]. The members of this pair are not connected and are found in a total of 1,504 co-cited pairs with frequencies ranging from 10 to 51,567. The second highest frequency, 28,407 co-citations, was for another pair of articles from the field of biochemistry: [laemmli_cleavage_1970, bradford_rapid_1976]. Members of this pair are not connected and are found in 41,909 co-cited pairs, 24,558 for the Laemmli gel electrophoresis article and 17,352 for the Bradford protein estimation article. In terms of this second pair, both articles describe methods heavily used in biochemistry and molecular biology, an area with strong referencing activity, so this result is not entirely surprising.

Having developed as a prediction of the probability that articles and would be co-cited, we first tested whether the distribution of co-citation frequencies was independent of (Materials and Methods). The null hypothesis that the co-citation frequency distribution was independent of was rejected with a very small -value: the statistical software indicated a -value with no significant non-zero digits. We next investigated what distribution functions might fit the frequencies of co-citation as varied.

Based on the long tails of citation frequencies, prior research has assessed the fit of lognormal and power law distributions [stringer_statistical_2008, radicchi_statistical_2008, stringer_statistical_2010]. We noted long right tails in co-citation frequencies, which, similarly, motivated us to assess the fit of lognormal and power law distributions to co-citation data. Further, we stratified the data according to (i) the minimum frequency for the right tail , (ii) , and (iii) whether the two members of each co-citation pair were connected. Figure 3 shows which distribution, if either, fits the data in each slice, based on tests of statistical significance. Note that there were no circumstances where both distributions fit: if one fit, then the other did not.

Statistical tests were not possible for some slices due to an insufficient number of data points. This was the case for certain combinations of large , large , and co-citations that were not connected. The number of data points obviously decreases as increases, and we found the decrease in the number of data points to be more precipitous when was large and co-citations were unconnected due to the lighter right tails for these parameter combinations. The graph in the right panel of Figure 4, which has a logarithmic -axis, shows that the number of data points per interval analyzed decreases most often by more than an order of magnitude from one interval to the next as increases. Most pairs of publications that are co-cited at least ten times, therefore, have small values of .

Figure 3 indicates when the null hypothesis of a best-fit lognormal or power law fitting the observed data can not be rejected. We computed two types of statistics for evaluating the null hypothesis ( and K-S) and, moreover, we computed the statistic for four binning strategies. Figure 3 indicates a distributional fit, specifically, if either the K-S -value is greater than 0.05 or if two or more of the statistics are greater than 0.05. While we computed the K-L Divergence (see supplementary material), we did not use these computations for formal statements of distributional fit because they are neither a norm nor do they determine statistical significance. These K-L computations did, however, support the findings based on formal tests of statistical significance.

Power law distributions fit most often when co-citations are connected (Fig. 3), when more extreme right tails are considered, and when co-citations have small values of . Lognormal distributions fit, conversely, in some circumstances, when a greater portion of the right tail is considered. These observations support the existence of heavy tails for small, even if a lognormal distribution fits the observed data more broadly. This observation is consistent with our observations of the most frequent co-citations having small values, as shown in the scatter plot in the left panel of Figure 4.

Mitzenmacher [Mitzenmacher_2003] shows a close relationship between the power law and lognormal distributions vis-à-vis subtle variations in generative mechanisms that determine whether the resulting distribution is power law or lognormal. The stratified layers in Figure 3 where a lognormal distribution fits for some portion of the right tail and, in the same instance, a power law describes the more extreme tail, may, therefore, be due to a generative mechanism whose parameters are close to those for a power law distribution as well as those for a lognormal distribution.

Right-tail cutoff () | Power law exponent () | |
---|---|---|

200 | 3.26 | |

200 | 3.37 | |

250 | 3.27 | |

250 | 3.37 | |

300 | 3.22 | |

300 | 3.35 |

Table 2 shows the exponents of the best-fit power law distributions when statistical tests indicated that a power law was a good fit and where comparisons were possible among the intervals of : these were possible for intervals of and , for connected co-citations, and right tails commencing at . The power law exponent in these comparisons was less for than for , indicating heavier tails for small and, therefore, a greater chance of extreme co-citation frequency. Figure 5 shows a log-log plot of the number of co-citations (-axis) exhibiting the counts on the -axis, for in the interval (note that both axes employ log scaling). The pattern for points below the 99th percentile clearly indicate that the number of co-citations referenced at a given frequency decreases greatly as the frequency increases. Also, the broadening of the scatter where fewer co-citations are cited more frequently is indicative of a long right tail, as has been observed in other research where lognormal or power law distributions have been fit to data, as in [montebruno_tale_2019].

Perline [perline_2005] warns against fitting a power law function to truncated data. Informally, a portion of the entire data set can appear linear on a log-log plot, while the entire data set would not. He cites instances where researchers have mistakenly characterized an entire data set as following a power law due to an analysis of only a portion of the data, when a lognormal distribution might provide a better fit to the entire data set. Indeed, the scatter plot in Figure 5 is not linear and so, as Figure 3 shows, a power law does not fit the entire data set. This is what Perline calls a weak power law where a power law distribution function fits the tail, but not the entire distribution. Our concern, however, is not with characterizing the distributional function for the entire data set, but with characterizing the features of high frequency co-citations, which by definition means we were concerned with the right tail of the distribution. Moreover, the results avoid confusion between lognormal and power law distribution functions because we have shown not only that a power law provides a statistically significant fit, but also that a lognormal distribution function does not fit.

Our analysis found particularly heavy tails that were well fit by power law distributions for small , in the intervals and , and for co-citations whose constituents are connected, as shown in Fig. 3. The closely related Matthew Effect [merton_1968], cumulative advantage [price_general_1976], and the preferential attachment class of models [barabasi_albert_2002] provide a possible explanation for citation frequencies following a power law distribution for some sufficiently extreme portion of the right tail. For greater values of , insufficient data in the right tails precludes a definitive assessment in this regard, although one might argue that the lack of observations in the tails is counter to the existence of a power law relationship. It is also noteworthy that the exponents we found for co-citations (Table 2) are close in value to those reported for citations by [price_general_1976] and [radicchi_statistical_2008].

*Delayed Co-citations* The delayed onset of citations to a well cited publication, also referred to as ‘Delayed Recognition’ and ’Sleeping Beauty’, has been studied by Garfield, van Raan, and others [garfield_1970, vanraan_2019, ke_defining_2015, vanraansleeping2004, li_distinguishing_2016, glanzel_myth_2004, bornmann_identifying_2018].
We sought to extend this concept to frequently co-cited articles. As an initial step, we calculated two parameters (Materials and Methods): (1) the beauty coefficient [ke_defining_2015] modified for
co-cited articles and (2) timelag , the length of time between first possible year of co-citation and the first year in which a co-citation was recorded. We further focused our consideration of delayed co-citations
to the 95th percentile or greater of co-citation frequencies in our dataset of 4.1 million co-cited pairs. Within the bounds of this restriction, 24 co-cited pairs have a beauty coefficient of 1,000 or greater and all 24 are
in the 99th percentile of co-citation frequencies. Thus, very high beauty coefficients are associated with high co-citation frequencies.

We also examined the relationship of with co-citation frequencies (Fig. 6) and observed that high values were associated with lower co-citation frequencies. These data in appear to be consistent with a report from van Raan and Winnink [vanraan_2019], who conclude that ‘probability of awakening after a period of deep sleep is becoming rapidly smaller for longer sleeping periods’. Further, when two articles are connected, they tend to have smaller values compared to pairs that are not connected in the same frequency range.

## Figures

(a) Co-citation Scopus frequency versus | (b) Number of co-cited pairs per interval |

## Conclusions

In this article, we report on our exploration of features that impact the frequency of co-citations. In particular, we wished to examine article pairs with high co-citation frequencies with respect to whether they originated from the same school(s) of thought or represented novel combinations of existing ideas. However, defining a discipline is challenging, and determining the discipline(s) relevant to specific publications remains a challenging problem. Journal-level classifications of disciplines have known limitations and while article-level approaches offer some advantages, they are not free of their own limitations [article_stasa].

Consequently, we designed , a statistic that examines the citation neighborhood of a pair of articles and to estimate the probability that they would be co-cited. Our approach has advantages compared to alternate approaches: it avoids the challenges of journal-level analyses, it does not require a definition of “discipline" (or “disciplinary distance"), it does not require assignment of disciplines to articles, it is computationally feasible, and, most importantly, it enables an evaluation that is specific to a given pair of articles.

We note that when and are from the same sub-field, then may be very large, and conversely, when and are from very different fields, it might be reasonable to expect that will be small. Thus, in a sense, may correlate with disciplinary similarity, with large values for reflecting conditions where the two publications are in the same (or very close) sub-disciplines, and small values for reflecting that the disciplines for the two publications are very distantly related. We also comment that in this initial study, we have not considered second-degree information, that is publications that cite publications that cite an article of interest.

Our data indicate that the most frequent co-citations occur when co-citations have small values of , as shown in Figure 4. Our study considered the hypothesis that the frequency distribution is independent of , but our statistical tests rejected this hypothesis, and showed instead that the frequency distribution is best characterized by a power law for small values of and connected publications, and in many other regions is best characterized by a lognormal distribution.

The observation that power laws are consistent with small values of and connected co-citations is consistent with the theory of preferential attachment for these parameter settings. To the extent that preferential attachment is the mechanism giving rise to a power law, this suggests that preferential attachment is, at least, stronger for small values and connected co-citations than for other parameter combinations, or that preferential attachment is not applicable to other parameter values.

Observing power laws, heavy tails, and pairs with extreme co-citation strength for small values of (i.e., pairs that have small a priori probabilities of being co-cited) may seem, on its face, paradoxical. One possible explanation of the pairs in the extreme right tail with both small and large co-citation strength is that those pairs represent novel combinations of ideas that, when recognized within the research community, catalyze an increased citation rate, consistent with preferential attachment coupled to time-dependent initial attractiveness [eom_2011] as an underlying generative mechanism. However, small values of do not guarantee a high co-citation count: indeed, even for small values of , co-citations with a power law predominantly have relatively low co-citation strength.

We also note the increasing proportion of connected pairs as the percentile for co-citation frequency increases (Fig. 2); this pair of parameters appears to be associated with a fertile environment where extremely high co-citation frequencies are possible. This observation raises the question of whether small values of and connected co-citations are associated with preferential attachment and, if a causal relationship exists, then how do and co-citation connection provide an environment supporting preferential attachment? A possibility is that one article in a co-cited pair citing the other makes the potential significance of the combination of their ideas apparent to researchers. The clear pattern of the highest frequency co-cited pairs typically having low values suggests that these pairs are highly cited and hence impactful because of the novelty in the ideas or fields that are combined (as reflected in low ). However, other factors should be considered, such as the prominence of authors and prestige of a journal [garfield_1980] where the first co-citation appears.

We did not apply field-normalization techniques when assembling the parent pool of 768,993 highly cited articles consisting of the top 1% of highly cited articles from each year in the Scopus bibliography. Thus, the highly co-cited pairs we observe are biased towards high-referencing areas such as biomedicine and parts of the physical sciences [small_citation_1980]. However, the dataset we analyzed has a lower bound of 10 on co-citation frequencies and includes pairs from fields other than those that are high referencing. For example, the maximum we observed in the dataset of 4.1 million pairs was 149 years, and is associated to a pair of articles independently published in 1840, establishing their eponymous Staudt-Clausen theorem [clausen_1840, vonstaudt_1840]; this pair of articles was apparently co-cited 10 times since their publication. A second pair of articles concerning electron theory of metals [drude_1900_1, drude_1900_2] was first co-cited in 1994 for a total of 109 times, with observed of 94 years. Both cases are drawn from mathematics and physics rather than the medical literature. They are also consistent with the suggestion that the probability of awakening is smaller after a period of deep sleep [vanraan_2019]. As we have defined , with its heavy penalty for early citation, we create additional sensitivity to coverage and data quality especially for pairs with low citation numbers. Indeed, for the Staudt-Clausen pair, a manual search of other sources revealed an article [carlitz_cvs_1961] in which they are co-cited. Both these articles were originally published in German and it is possible that additional co-citations were not captured. Thus, big data approaches that serve to identify trends should be accompanied by more meticulous case studies, where possible. Other approaches for examining depth of sleep and awakening time should certainly be considered [vanraansleeping2004, ke_defining_2015]. Lastly, using our approach to revisit invisible colleges [price_beaver_1966, crane_1972, small_clustering_1985] seems warranted, since it seems likely that the upper bound of a hundred members predicted by [price_beaver_1966] is likely to have increased in a global scientific enterprise with electronic publishing and social media.

Finally, we view these results as a first step towards further investigation of co-citation behavior, and we introduce a new technique based on exploring first-degree neighbors of co-cited publications; we are hopeful that this graph-theoretic study will stimulate new approaches that will provide additional insights, and prove complementary to other article level approaches.

## Acknowledgments

In addition to support through federal funding, the ERNIE project features a collaboration with Elsevier. We thank our colleagues from Elsevier for their support of the collaboration.

## Competing Interests

The authors have no competing interests. Scopus data used in this study was available to us through a collaborative agreement with Elsevier on the ERNIE project. Elsevier personnel played no role in conceptualization, experimental design, review of results, or conclusions presented. The content of this publication is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health or Elsevier. Sitaram Devarakonda’s present affiliation is Randstad USA. His contributions to this article were made while he was a full-time employee of NET ESolutions Corporation.

## Author Contributions

Conceptualization, GC, JB, SD, and TW; Methodology, AD, DK, GC, JB, SD, SL, and TW; Investigation, DL-H, GC, JB, and SD; Writing -Original Draft, GC, JB, TW; Writing- Review and Editing, AD, DK, DL-H, GC, JB, SD, SL, and TW; Funding Acquisition, GC; Resources, DK and GC; Supervision, GC. Authors are listed in alphabetic order.