Temporal Evolution of Bradford Curves in Academic Library Contexts

Abstract

Bradford’s law of bibliographic scattering is a fundamental principle in bibliometrics, offering valuable guidance for academic libraries in literature search and procurement. However, Bradford curves can exhibit various shapes over time, and predicting these shapes remains a challenge due to a lack of causal explanation. This paper attributes the deviations from the theoretical J-shape to integer constraints on the number of journals and articles, extending Leimkuhler’s function to encompass highly productive core journals, where the theoretical journal number falls below one. Using the Simon–Yule model, key parameters of the extended formulas are identified and analyzed. The paper explains the reasons for the Groos droop and examines the critical points for shape changes. The proposed formulas are validated with empirical data from the literature, demonstrating that this method can effectively predict the evolution of Bradford curves, providing academic libraries with a valuable tool for evaluating journal coverage, optimizing resource allocation, and refining Collection Development Policies (CDP).

1. Introduction

1.1. Background and Significance

Bradford’s law, a foundational principle in bibliometrics, holds significant theoretical value in modeling and simulating the dynamics of scientific knowledge production [1,2,3,4,5,6,7,8]. By integrating it into the framework of science dynamics, this research reveals how micro-level behaviors influence macro-level trends in knowledge growth and shifts in core literature. However, generating Bradford curves can be labor-intensive, especially for journals with few relevant papers. Additionally, as the scientific literature in a given discipline often grows exponentially or passes through various developmental stages [9], a Bradford curve created at one point in time may become outdated without adjustments. Some mathematical models predict a J-shaped curve; but, in practice, Bradford curves can take on at least six different shapes, including the S-shaped curve with the so-called Groos droop [10]. There is a lack of causal explanations for these variations, and comprehensive empirical studies are limited [11]. As such, predicting the evolution of Bradford curves remains an open question requiring further investigation.

This paper attributes the different shapes of Bradford curves to the integer constraints of journal and paper numbers. If journal productivity $n$ is high enough that the corresponding theoretical journal number $f_t\left(n\right)=C/n^{\alpha }$ falls below one, the actual journal number $f_e\left(n\right)$ can only be zero or one. This discrete nature causes deviations in the core zone from the theoretical predictions of models like Lotka or Simon–Yule. To address this, the paper proposes two distinct formulas for the core and the normal zones, which are analyzed through theoretical methods and Monte Carlo simulations of the Simon–Yule model. The causes of the Groos droop are explained, and the critical points for shape changes are identified. Finally, the proposed formulas are validated using empirical data, demonstrating that this method can predict the evolution of Bradford curves.

While Bradford’s law may not directly guide acquisition strategies in modern academic libraries, it still provides a valuable framework for evaluating journal coverage and understanding the distribution of core literature within specific fields [12,13,14,15,16]. This suggests its potential for use in ensuring comprehensive and balanced collections, contributing to a more strategic approach to library resource management.

1.2. Literature Review

Bradford’s law was first proposed by Bradford in 1934 [17] but did not gain wide recognition until Vickery further developed the theory in 1948 [18]. According to Bradford’s law, if journals are arranged in descending order of productivity and divided into $p$ groups with the same number of papers, the number of journals in each group $n_i$ follows the ratio $n_1:n_2:\cdots :n_p=1:k:\cdots :k^{p-1}$, where $k$ is the Bradford multiplier. Besides this verbal form, Bradford’s law can also be depicted as a J-shaped curve by plotting the accumulated productivity $R\left(r\right)$ of the first $r$ journals against the natural logarithm of the journal rank $r$. Leimkuhler proposed the mathematical formula for Bradford curve in 1967 [19], and Egghe developed a method for determining the parameters of this formula in 1990 [20]. In Leimkuhler’s function $R\left(r\right)=a\mathrm{l}\mathrm{o}\mathrm{g}\left(1+br\right)$, where the key parameters $a$ and $b$ can be calculated from the article number $A$, journal number $T$, and the productivity $y_m$ of the most productive journal. Although Leimkuhler’s function matches well with many bibliographies, it corresponds to a J-shaped curve, which deviates from those with a Groos droop [20].

Incomplete bibliographies were initially believed to cause the Groos droop, but further research refuted this hypothesis [21]. Egghe demonstrated that if the ranking of each journal $r$ is transformed into $r^{\prime }=r+r_0$ by adding a large constant $r_0>1/b$, then the new curve will concave downwards, showing a Groos droop [22]. The merging of different bibliographies, each with a different maximum journal productivity $y_m^{\left(i\right)}$, could explain the large constant $r_0$ [22]. However, it is also likely that the large core regions (regions with the most productive journals where $f_t\left(n_i\right)<1$) of some bibliographies contribute to the large $r_0$ [23]. Essentially, $y_m$ in Leimkuhler’s function denotes the journal productivity, where $f_t\left(y_m\right)=C/y_m^{\alpha }\approx 1$ [24], rather than the maximum yield $X_1$ of a journal as claimed by Egghe himself. Thus, if the total number of these journals $T_0$ exceeds the critical value $r_0=1/b$, a Groos droop will emerge. This paper adopts this explanation and extends Leimkuhler’s function to predict the evolution of Bradford curves.

In the 1990s, research interest in Bradford’s law shifted from the static presentation of data at a particular time to its dynamic and evolutionary aspects [25]. Oluić-Vuković studied how the increase in productivity of core journals affected the shape of the distribution curve over time [26]. By analyzing the research output of Croatian scholars in different subjects, she concluded that the Groos droop or S-shaped curve is caused by an increase in the concentration/dispersal disparity, reflected by the rise in the core/periphery ratio [27]. The dynamic evolution of Bradford curves and the emergence of the Groos droop were presented in her 1992 study [28], and other similar empirical studies partitioning bibliographies over time were conducted by Garg [29], Wagner-Döbler [11] and Sen [30].

Meanwhile, stochastic models like the Simon–Yule model have increasingly been used to study the dynamic characteristics of bibliometric laws [25,31]. Initially introduced by Yule in 1924 for studying the distribution of biological genera by species number, the Simon–Yule model gained recognition when Simon expanded it in 1955 to analyze the frequency distributions of words in writing samples [32]. Besides employing theoretical methods for precisely solving the constant entry rate $\alpha$ of new sources [32], Monte Carlo simulations have been used to explore more complex scenarios, such as declining entry rates ${\alpha }_t$ [33] and autocorrelated growth rates $\gamma$ of established journals (also referred to as aging or obsolescence rate) [34]. Chen et al. [35,36,37] first used the Simon–Yule model to study the evolution of Lotka’s and Bradford’s laws over time. They found the entry rate ${\alpha }_t$ and the autocorrelated growth rate $\gamma$ have significant yet opposite effects on the Bradford curves, offering an explanation for the various types of Bradford curves [36].

Later, Oluić-Vuković also explored the dynamics of Bradford distributions using the Simon–Yule model but found that its steady-state solution was too restrictive to handle time variations, limiting its applicability [25,31]. This paper also utilizes the Simon–Yule model to examine different scenarios’ effects on key parameters (e.g., journal number $T_0$, article number $A_0$, and maximum productivity $X_1$ of the core region) of the extended Leimkuhler’s function. However, it is not used directly to forecast the evolution of Bradford curves or to compare them with empirical data. Instead, key parameters are estimated from past empirical data to improve predictions of Bradford curve evolution in the future.

2. Theoretical Study

2.1. Simon-Yule Model

Simon’s generating mechanism for the Bradford distribution is based on the following two assumptions, where $f_t\left(n,t\right)$ denotes the number of journals that have published exactly $n$ papers in the first $t$ published papers:

• There is a constant probability α that the $\left(t+1\right)$-th paper is published in a new journal—a journal that has not published in the first t papers.
• The probability that the $\left(t+1\right)$-th paper is published in a journal that has published $n$ papers is proportional to $nf\left(n,t\right)$—that is, to the total number of papers of all journals that have published exactly $n$ papers.

Therefore, if there are $A$ papers at a given time, then the corresponding journal number $T$ is approximately $T=A\alpha$. Based on Simon’s two assumptions, the steady-state solution of the Bradford distribution can be written as [37]:

$$f_t\left(n\right)=\rho B\left(n,\rho +1\right)\approx \rho \mathsf{\Gamma }\left(\rho +1\right)n^{-\left(\rho +1\right)}$$

(1)

where $B$ is the beta function, $\mathsf{\Gamma }$ is the gamma function, and $\rho$ is a function of the entry rate of new journal $\alpha$, defined as $\rho =1/\left(1-\alpha \right)$. Equation (1) suggests that the analytical outcome of the Simon–Yule model aligns with Lotka’s law when $\rho \approx 1$.

In addition to the analytical solutions, Monte Carlo simulations are conducted for $\alpha =0.15$, and the results are compared with the theoretical results of Equation (1), as shown in Figure 1. Detailed procedures for these simulations can be found in [33] and the data and MATLAB code used are available in the Supplementary Materials. To reduce inherent randomness, each case is simulated $N={10}^4$ times, utilizing only the medians of these simulations as the final outputs.

Figure 1 illustrates distinct zones in the simulation results: a normal zone (blue circles) and a core zone (red squares). This distinction arises from the necessity for actual journal numbers $f_e\left(n\right)$ to be integers, unable to fall below one. Thus, when the journal productivity $n$ is high enough for the theoretical journal number $f_t\left(n\right)$ to drop below one, actual journal numbers $f_e\left(n\right)$ are constrained to zero or one, deviating from theoretical predictions as shown by the red squares and black dots in Figure 1. Moreover, despite the smaller number of journals in the core region, their contribution to the paper count is significant, as shown in Figure 1b. Hence, accurate prediction of the paper count for each journal $X_r$ in the core zone is crucial for depicting the Bradford curve faithfully.

Estimating journal productivity $X_r$ of the core region involves determining the journal count $T_0$ and the paper count $A_0$ first. Since Figure 1 shows a close match between the journal number $f\left(n\right)$ and the paper number $nf\left(n\right)$ in the normal region, the total journal count $T_1$ and the total paper count $A_1$ of the normal region can be directly obtained by summing all journal and paper numbers. These are calculated as $T_1=\sum _{n=1}^{y_m}f\left(n\right)$ and $A_1=\sum _{n=1}^{y_m}nf\left(n\right)$, where $y_m$ is the journal productivity when $f_t\left(y_m\right)\approx 1$. In this context, $y_m$ can be seen as the productivity of the most productive journal in the normal region, as indicated by the black diamond in Figure 1. According to Equation (1), the analytical expression for $y_m$ can be derived as:

$$y_m={\left[A\left(\rho -1\right)\mathsf{\Gamma }\left(\rho +1\right)\right]}^{{\displaystyle \frac{1}{\rho +1}}}$$

(2)

Once $y_m$ is calculated, the total journal count $T_0$ and the paper count $A_0$ of the core region can be calculated as $T_0=T-T_1$ and $A_0=A-A_1$. Alternatively, they can be directly calculated as:

$$T_0\approx {\int }_{y_m}^{+\infty }Tf\left(n\right)\mathrm{d}n={\displaystyle \frac{y_m}{\rho }}$$

(3)

$$A_0\approx {\int }_{y_m}^{+\infty }Tnf\left(n\right)\mathrm{d}n={\displaystyle \frac{y_m^2}{\rho -1}}$$

(4)

where $y_m$ is calculated using Equation (2).

In the Simon–Yule model with constant entry rate $\alpha$, the maximum number of papers one journal can have, $X_1$, can be estimated using Gumbel’s $r$-th characteristic extreme theory [38,39,40]:

$$G\left(X_r\right)\approx {\int }_{X_r}^{+\infty }f\left(i\right)\mathrm{d}i={\displaystyle \frac{r}{A}}$$

(5)

By solving this equation, it can be derived that the productivity of the most productive journal, $X_1$, can be expressed as:

$$X_1={\left[A\mathsf{\Gamma }\left(\rho +1\right)\right]}^{{\displaystyle \frac{1}{\rho }}}={\left(\rho -1\right)}^{-{\displaystyle \frac{1}{\rho }}}{y}_m^{{\displaystyle \frac{\rho +1}{\rho }}}$$

(6)

The productivity of the $r$-th most productive journal, $X_r$, is related to $X_1$ by $X_r=X_1{r}^{-1/\rho }$. Figure 2 compares Gumbel’s $r$-th characteristic extreme values with the medians of the simulation results. While Gumbel’s theory effectively predicts the largest paper number $X_1$, it falls short in estimating other paper numbers in the core zone $X_r$, for $r=2,3,\cdots ,T_0$. Thus, an alternative method is assumed in this paper, where all other $X_r$, for $r=2,3,\cdots ,T_0$, are related to the largest paper number $X_1$ through the equation:

$${\displaystyle \frac{X_1}{X_r}}=k\left(r-1\right)+1$$

(7)

where $k$ is the only parameter waiting to be determined. The validity of this equation is supported by Figure 2b, where the blue circles represent the simulation results, and the blue dashed lines show the linear fitting results. Hence, the productivity of the $r$-th most productive journal can be derived from Equation (7), and the cumulative productivity of the first $r$ most productive journals can be written as:

$$R_c\left(r\right)=\sum _{i=1}^r{\displaystyle \frac{X_1}{k\left(i-1\right)+1}}$$

(8)

If there are $T_0$ journals with $A_0$ papers in the core region and the numbers of $T_0$ and $A_0$ are known, then the parameter $k$ can be calculated from the equation $R_c\left(T_0\right)=A_0$. Then, Equation (8) can be used to predict the evolution of the core regions ($r\le T_0$) of Bradford curves.

2.2. Leimkuhler’s Function

After removing the $T_0$ journals and the $A_0$ papers of the core region, the remaining $T_1$ journals and $A_1$ papers align well with the theoretical results predicted by Equation (1). Consequently, they follow Lotka’s law, and their Bradford curve can be predicted using the revised Leimkuhler’s function [20]:

$$R\left(r_1\right)=a\mathrm{l}\mathrm{o}\mathrm{g}\left(1+b{r}_1\right)$$

(9)

where the key parameters $a$ and $b$ are defined as:

$$a={\displaystyle \frac{A_1}{\mathrm{l}\mathrm{o}\mathrm{g}\left(e^{\gamma }{y}_m\right)}}$$

(10)

$$b={\displaystyle \frac{e^{\gamma }{y}_m-1}{T_1}}$$

(11)

where $\gamma$ is the Euler–Mascheroni constant, $\gamma \approx 0.5772$, and $y_m$ is the journal productivity when the corresponding theoretical journal number $f_t\left(y_m\right)\approx 1$. While $y_m$ can be directly calculated from Equation (2), it can also be estimated using the following equation if the values of $X_1$, $T_0$ and $A_0$ are known:

$$y_m\approx {\displaystyle \frac{X_1}{k\left(T_0-1\right)+1}}$$

(12)

Since the core region’s journal productivity is higher than the normal region’s, these significant journals rank lower. Consequently, the Bradford curve for the normal region starts at the point $\left(T_0,A_0\right)$. Each rank $r_1$ in the normal region should be transformed to $r=r_1+T_0$, and the cumulative productivity of the first $r$ journals $R\left(r\right)$ should be transformed into $R_n\left(r\right)=R\left(r_1\right)+A_0$. The revised Leimkuhler’s function for the normal region is then written as:

$$R_n\left(r\right)=R\left(r-T_0\right)+A_0=a\mathrm{l}\mathrm{o}\mathrm{g}\left[1+b\left(r-T_0\right)\right]+A_0$$

(13)

Equation (13) can be used to predict the dynamic evolution of the normal regions ($T_0<r<T$) of the Bradford curve. Therefore, Equations (8) and (13) together can be used to predict the dynamic evolution of the Bradford curves.

Figure 3a displays the Bradford curve. The blue circles represent the normal zone, while the red squares indicate the core zone. The blue dashed lines show the prediction results of Equation (13), and the red dotted lines show the prediction results of Equation (8). The black upper triangle, the diamond, and the lower triangle represent the points $\left(1,X_1\right)$, $\left(T_0,A_0\right)$, and the $\left(T,A\right)$, respectively. From the discussion above, it is evident that these three points and two lines are crucial for predicting the evolution of the Bradford curve.

2.3. Groos Droop

Groos was the first to observe that in some datasets, when the journal productivity is low, the Bradford curve tends to bend downwards [10]. Egghe [22] explained the cause of the Groos droop as merging datasets. However, this section demonstrates that the core region’s existence causes the Groos droop in the normal region.

The first and second derivatives of $R_c\left(r\right)$ for the core region can be derived from Equation (8):

$${\displaystyle \frac{\partial R_c\left(r\right)}{\partial \left(\mathrm{l}\mathrm{o}\mathrm{g}r\right)}}={\displaystyle \frac{X_{1}r}{k\left(r-1\right)+1}}$$

(14)

$${\displaystyle \frac{{\partial }^2{R}_c\left(r\right)}{\partial {\left(\mathrm{l}\mathrm{o}\mathrm{g}r\right)}^2}}={\displaystyle \frac{X_1\left(1-k\right)r}{{\left[k\left(r-1\right)+1\right]}^2}}$$

(15)

From Equation (15), it can be noted that when $k>1$, ${\displaystyle \frac{{\partial }^2{R}_c\left(r\right)}{\partial {\left(\mathrm{l}\mathrm{o}\mathrm{g}r\right)}^2}}<0$, the Bradford curve for the core region concaves downwards. Conversely, when $k<1$, ${\displaystyle \frac{{\partial }^2{R}_c\left(r\right)}{\partial {\left(\mathrm{l}\mathrm{o}\mathrm{g}r\right)}^2}}>0$, the curve concaves upwards. As the entry rate of new journals $\alpha$ increases, the number of journals $T$ rises, and the distribution of articles become more dispersed, leading to a decrease in the largest journal productivity $X_1$. From Equation (8) and $R_c\left(T_0\right)=A_0$, we know that a lower $X_1$ results in a lower $k$ if $A_0$ is relatively constant. Therefore, as $\alpha$ increases, the Bradford curve for the core region will gradually concave upwards, as shown in Figure 3.

Similarly, the first and second derivatives of $R_n\left(r\right)$ for the normal region can be derived from Equation (13):

$${\displaystyle \frac{\partial R_n\left(r\right)}{\partial \left(\mathrm{l}\mathrm{o}\mathrm{g}r\right)}}={\displaystyle \frac{abr}{b\left(r-T_0\right)+1}}$$

(16)

$${\displaystyle \frac{{\partial }^2{R}_n\left(r\right)}{\partial {\left(\mathrm{l}\mathrm{o}\mathrm{g}r\right)}^2}}={\displaystyle \frac{ab\left(1-b{T}_0\right)r}{{\left[b\left(r-T_0\right)+1\right]}^2}}$$

(17)

From Equation (17), it can be noted that when $T_0>1/b$, ${\displaystyle \frac{{\partial }^2{R}_c\left(r\right)}{\partial {\left(\mathrm{l}\mathrm{o}\mathrm{g}r\right)}^2}}<0$, the Bradford curve for the normal region will concave downwards, showing a Groos droop. When $T_0<1/b$, ${\displaystyle \frac{{\partial }^2{R}_c\left(r\right)}{\partial {\left(\mathrm{l}\mathrm{o}\mathrm{g}r\right)}^2}}>0$, the curve will concave upwards, forming a J-shaped curve. As the entry rate of new journals $\alpha$ increases, Figure 3a shows that the $T_0$ will eventually fall below $1/b$, causing the Bradford curve for the normal region to concave upwards, similar to the core region.

Figure 3b illustrates the variation of key parameters $T_0$, $1/b$ and $k$ with the entry rate $\alpha$. When $A={10}^4$, the normal region will start to concave upwards at critical point ${\alpha }_n\approx 0.2$, while the core region will do so at ${\alpha }_c\approx 0.3$. Therefore, when $\alpha <0.2$, the entire Bradford curve will concave downwards; when $\alpha >0.3$, it will concave upwards. For $0.2<\alpha <0.3$, the Bradford curve will exhibit a reversed S-shape, with the core region concaving downwards and the normal region concaving upwards. Figure 3a shows the three shapes of Bradford curves. In this specific case, since ${\alpha }_n<{\alpha }_c$, there is no S-shaped Bradford curve. This is because the aging of journals is not considered here, making the largest journal productivity $X_1$ relatively large. When aging effects are considered, which will be discussed in detail in Section 3.2, $X_1$ decreases significantly, which results in a lower $k$ and thus a much lower ${\alpha }_c$. If ${\alpha }_c<{\alpha }_n$, an S-shaped Bradford curve will appear for ${\alpha }_c<\alpha <{\alpha }_n$, with the core region concaving upwards and the normal region concaving downwards.

2.4. Bradford Dynamics

Given the analytical expression of $T_0$, $A_0$ and $X_1$ (Equations (3), (4) and (6)), the core region of the Bradford curve at any time can be predicted using Equation (8). The parameters for the normal regions can be derived from these factors through $T_1=T-T_0$, $A_1=A-A_0$ and Equation (12), allowing the normal region of the Bradford curve to be predicted using Equation (13).

The Bradford curves for the constant entry rate scenario are shown in Figure 4a. Here, the red squares and blue circles represent the simulation results for the core and the normal regions, respectively, while the red dashed lines and the blue dotted lines represent the theoretical results of Equations (8) and (13). The key points $\left(1,X_1\right)$, $\left(T_0,A_0\right)$, and $\left(T,A\right)$ are also shown as the black upper triangle, the diamond and the lower triangle in Figure 4a. It is notable that although the core region contains far fewer journals than the normal region, its representation is significant due to the x-axis’s log scale. Consequently, journals with lower ranks are better represented in Figure 4a.

Figure 3 shows that when $\alpha =0.15$, the entire Bradford curve concaves downwards for ${10}^3<A<{10}^4$. Figure 4a depicts the evolution of Bradford curves as the paper number $A$ increases from ${10}^3$ to ${10}^4$, aligning well with theoretical predictions.

Figure 4b presents the simulation and the analytical results for $T_0$, $A_0$ and $X_1$. Hollow symbols represent the simulation results, while solid symbols denote the analytical ones. It can be observed that these three key factors are linear functions of the paper number $A$, as indicated by Equation (18):

$$\mathrm{l}\mathrm{o}\mathrm{g}\left(Y\right)=a_{\rho }+b_{\rho }\mathrm{l}\mathrm{o}\mathrm{g}\left(A\right)$$

(18)

where the constants $a_{\rho }$ and $b_{\rho }$ are functions of $\rho$. Since $\rho$ is approximately one, $a_{\rho }$ and $b_{\rho }$ can be considered constants. The close match between the analytical and the numerical results confirms the validity of Equations (3), (4) and (6). The analytical result for $T_0$ is slightly lower than the numerical ones because, in the numerical results, all journals with productivity $n$ that satisfy $f\left(n\right)=1$ are considered part of the core region; whereas, in theory, some of them belong to the normal region.

3. Numerical Study

3.1. Decreasing Entry Rate

Assume the probability of adding a new journal decreases linearly over the time

$$\alpha \left(t\right)={\alpha }_s-kt$$

(19)

where $k$ is a constant, $k=\left({\alpha }_s-{\alpha }_f\right)/A_f$, where ${\alpha }_f$ and $A_f$ are the entry rate of new journals and the total number of articles in the final state, respectively. The accumulated number of journals can then be expressed as:

$$T=\sum _{t=1}^A\alpha \left(t\right)={\alpha }_{s}A-{\displaystyle \frac{1}{2}}k{A}^2$$

(20)

Using Equation (20), a quadratic fitting of $T$ and $A$ allows us to determine the values of ${\alpha }_s$ and ${\alpha }_f$. Once these are known, the average entry rate $\overline{\alpha }=\left({\alpha }_s+{\alpha }_f\right)/2$ can be used to calculate the analytical results.

Figure 5 shows the dynamic evolution of Bradford curves and the variations of key parameters when the entry rate decreases linearly from 0.2 to 0.1. The proposed method effectively predicts these variations, with analytical results using $\overline{\alpha }$ matching well with the simulation results. Although the numerical results for $A_0$ and $X_1$ are slightly lower than the analytical ones, this suggests that a decreasing entry rate has a slight negative impact on their increase. However, the effect on the overall shape of the Bradford curve and key parameters is relatively insignificant, indicating that the analytical results for a constant entry rate $\overline{\alpha }$ can still be used to predict key parameters without significant errors.

3.2. Aging Rate of Journals

Simon assumes that only one paper gets published in each time period and models the probability of a journal increasing in paper number as proportional to a weighted sum of its past increments. These increments are weighted by a factor that decreases geometrically over time, with the rate of decrease denoted as $\gamma$.

Let $y_j\left(k\right)$ represent the change in the paper number of the $j$-th journal during the $k$-th time interval, where $y_j\left(k\right)$ is either 1 (indicating a unit increment) or 0 (indicating no change). The paper number of the $j$-th journal at the end of the $k$-th interval is given by $\sum _{\tau =1}^k{y}_j\left(\tau \right)$. The expected increment in the paper number during the $\left(k+1\right)$-th interval is:

$$p\left[y_j\left(k+1\right)=1\right]={\displaystyle \frac{1}{W_k}}\sum _{\tau =1}^k{y}_j\left(\tau \right){\gamma }^{k-\tau }$$

(21)

where $W_k$ is a time-dependent function consistent across all journals, defined as $W_k=\sum _{j=1}^T{w}_j\left(k\right)$ with $w_j\left(k\right)=\sum _{\tau =1}^k{y}_j\left(\tau \right){\gamma }^{k-\tau }$. The parameter $\gamma$ determines how quickly the influence of past growth is diminished and is thus referred to as the aging rate of journals in this paper.

Figure 6 illustrates the impact of the aging rate of journals $\gamma$ on the dynamics of Bradford curves and key parameters. It shows that the aging factor increases $T_0$, making the normal region concave downwards more. The aging effect also markedly reduces $X_1$ by weakening the Mathew effect, where successful journals attract more papers. As older journals lose their appeal, this “success breeds success” effect diminishes, leading to a substantial decrease in $X_1$, as shown in Figure 6a. Consequently, while the number of articles $A_0$ in the core zone remains relatively unchanged, $T_0$ must increase to offset the reduction in $X_1$. This expansion of the core region increases its share of the Bradford curve, shaping it into a J-shape due to the reduction of $k$ and $X_1$, as discussed in Section 2.3. Additionally, with the increase in $T_0$, it is more likely that $T_0$ will exceed $1/b$, further contributing to the concave downwards shape of the normal region. In summary, the aging effect of journals facilitates the Bradford curve to more easily adopt an S-shape.

3.3. Varying Entry and Aging Rates

Figure 7 shows the impact of varying entry and aging rates on Bradford curves. In real-world scenarios, both the entry rate and the aging rate often change steadily. For example, the entry rate might decrease linearly from 0.2 to 0.1, while the aging rate increases linearly from 0.95 to 1.0, and the simulation results are shown in Figure 7. Comparing Figure 5, Figure 6 and Figure 7 reveals that the effects of decreasing entry rate and increasing aging rate are similar to those observed with constant entry and aging rates (Figure 6). However, in Figure 7b, both the article number $A_0$ and journal number $T_0$ are even lower compared to Figure 6b, indicating that the decreasing entry rate further exacerbates their decrease. Consequently, the normal region of Bradford curve becomes less concave downwards, and its starting point on the y-axis is notably lower. Importantly, all three key factors continue to show linear relationships with the article number, suggesting that Equation (18) remains useful for predicting them.

4. Empirical Study

4.1. Dataset of Croatian Chemistry Research

Oluić-Vuković used the research output in chemistry by authors from Croatia to prepare full bibliographic references for a ten-year period [28]. This dataset includes only articles published in journals, comprising 2543 papers across 416 journals over a decade. The productivity of the top few (fewer than 10) most prolific journals was taken directly from Figure 1 in [28], while the productivity of other journals was taken from Tables 4 and 6 in [25].

In these tables, the journal productivity $n$ and the number of journals $f\left(n\right)$ are recorded for nine cumulative time periods. For journals with higher productivity (where $n$ correspond to $f\left(n\right)=1$), Oluić-Vuković grouped the data into ranges (e.g., $n\ge 77$), which limits the precision of the distribution for high-producing journals. To accurately represent this higher productivity range, we relied on the detailed Bradford curves provided in [28], adjusting the data to match the total number of journals $T$ and articles $A$ as indicated in Table 2 of [28].

As an example, using the data from the ninth interval (also found in Table 1 of reference [27]), we identified the highest productivity level $X_1=301$ by observing the last value in the sequence of $n$. The dividing point between the core and the normal zones was set at $n=24$ based on the fact that for $n<24$, most values of $f\left(n\right)$ were greater than 1, and $n$ values were continuous; while for $n>24$, $f\left(n\right)$ showed irregular gaps with isolated high $n$ values and lower average journal counts.

The journal count $T_0$ and the article count $A_0$ for the core zone were calculated by summing the respective $f\left(n\right)$ and $n·f\left(n\right)$ for $n\ge 24$. Similarly, the total number of journals $T$ and articles $A$ were obtained by summing the values across all n, ensuring consistency with the overall dataset.

To predict the dynamics of the Bradford curve, the process begins by predicting the variation of the total article number $A\left(t\right)$ over time $t$. Logistic regression analysis [41] was applied to the empirical data to predict the total number of articles $A\left(t\right)$ at any given time, as shown in Figure 8a. Logistic regression was applied using MATLAB to model the growth of the article output based on empirical data, providing a reliable prediction of future article numbers.

Next, the total journal number $T$ and the entry rate of new journals $\alpha$ were estimated by plotting the total journal number $T$ against the total article number $A$ and a linear fit of $T=A\alpha$ was applied, as shown in Figure 8b. The linear fit, with an R-squared value of 0.98, confirmed a strong relationship between $T$ and $A$, providing insights into how new journals enter the system as article output increases.

Once the point $\left(T,A\right)$ is determined for any time, linear regression of Equation (18) is used to determine the three key parameters $T_0$, $A_0$, and $X_1$ on a log-log axis, with the fitting results shown in Figure 9a. These parameters reveal the structure of the core and the normal regions of the Bradford distribution.

Finally, based on the three key points $\left(T_0,A_0\right)$, $\left(1,X_1\right)$, and $\left(T,A\right)$, Equations (8) and (13) were used to model the Bradford curve for the core region and the normal region, respectively. The results, illustrated as red dashed lines and blue dotted lines in Figure 9b, show the gradual transition in the curve’s shape. This approach allows librarians and practitioners to model journal productivity trends, helping to optimize collection development strategies and predict journal coverage effectively.

The Bradford curves shown in Figure 9b align closely with the empirical data, demonstrating a strong match between the predicted and the observed outcomes. Notably, the Bradford curve transitions gradually from a J-shape to an S-shape, a transformation that is accurately captured by the analytical predictions.

4.2. Dataset of Solar Power Research

The bibliographies on solar power research for the years 1971, 1974, 1977, 1980, 1983, and 1986 were compiled by Garg et al. [29], encompassing papers published in journals from the Engineering Index. The data for this analysis were directly extracted from Tables 1–7 of the referenced study. Unlike the Croatian Chemistry dataset, Garg’s dataset is clearer, as it records every journal’s exact productivity $n$ without grouping high-productivity journals.

The dataset is divided into six cumulative periods: 1971, 1971 + 1974, 1971 + 1974 + 1977, 1971 + 1974 + 1977 + 1980, 1971 + 1974 + 1977 + 1980 + 1983, and 1971 + 1974 + 1977 + 1980 + 1983 + 1986. Although the data skips certain years (e.g., 1972 and 1973), this does not affect the Bradford curve calculation. However, when predicting article growth, care must be taken to sum the predicted outputs for each year to match the cumulative totals.

The data processing followed the same four-step method as the Croatian Chemistry dataset. First, logistic regression was applied to predict the cumulative article numbers $A\left(t\right)$ over time $t$, with data from Table 7 used to estimate the article numbers for each desired interval, as shown in Figure 10a. This step involves fitting the empirical data and adjusting for the cumulative nature of the dataset, ensuring accurate predictions for the individual years within each period.

Next, quadratic fitting (Equation (20)) was applied to the journal and article pairs to estimate the total number of journals $T$ and the entry rate of new journals $\alpha$, as shown in Figure 10b. Logistic fitting was applied to better capture the effect of a linearly decreasing journal entry rate, providing a more accurate representation of how new journals enter the system as article output increases.

The key parameters $T_0$, $A_0$, and $X_1$ were then predicted using linear fitting on the log-log axis, based on the empirical data from the earlier steps. This method, depicted in Figure 11a, was used to determine these parameters at each cumulative time point, allowing the segmentation of the core and the normal regions of the Bradford curve.

Finally, using Equations (8) and (13), Bradford curves for the core and the normal regions were plotted, as shown in Figure 11b. As with the Croatian Chemistry dataset, determining the highest productivity level $X_1$, the core zone journal count $T_0$, and the core zone article count $A_0$ was straightforward in Garg’s dataset due to the clear distribution of journals and articles. This simplicity contrasted with the Croatian dataset, where grouping high-productivity journals added complexity.

Figure 9b and Figure 11b demonstrate that while the proposed method can predict the general trend of the Bradford curves, there are inherent errors. These errors arise because Bradford’s law inherently contains uncertainties. Numerical studies reveal that the article numbers for each journal rank have a large standard deviation, making it practically impossible to predict the precise shape of the Bradford curve. Additionally, the various fitting procedures introduce errors into the process. Therefore, this method can only predict the general trend of Bradford dynamics but cannot accurately predict the article number for each journal at any given time.

Another issue with this method is that the first derivatives of the core region (Equation (14)) and the normal region (Equation (16)) differ at the $\left(T_0,A_0\right)$ point, resulting in a non-smooth analytical curve at the intersection point. In contrast, the numerical simulation results are smooth throughout. This problem could be addressed by proposing more complex formulas for the normal region, but this would complicate the overall method. Given the difficulty in accurately predicting Bradford curve dynamics, this aspect is not explored further in this paper.

5. Discussion

5.1. Theoretical Contributions

This study makes contributions in two areas: advancing bibliometric modeling and simulation and addressing extreme value phenomena. Both are critical to understanding the dynamics of scientific knowledge production.

• Modeling and Simulation of Science or Publications

Modeling and simulations in bibliometrics help uncover micro-level behaviors that drive broader trends, such as publication growth and shifts in core literature [42]. These methods offer controlled environments to reduce data biases and explore extreme scenarios, providing valuable insights for decision-making. While fields like network science, agent-based modeling (ABM), and diffusion models have been widely used in science modeling and simulation [43], the contributions of bibliometric research in these areas have not been fully explored. This study bridges that gap by incorporating models like Bradford’s law into the broader science dynamics framework, showing their value in understanding knowledge production over time.

2.

Extreme Value and Non-Lotkaian informetrics

This study addresses extreme value phenomena—highly productive authors, highly cited papers, and influential journals—that traditional models like Lotkaian informetrics struggle to capture [44,45]. These extreme cases deviate from typical productivity patterns, necessitating new theoretical approaches. By applying the extreme value theory, order statistics, and the Simon–Yule model, this research provides a more accurate model for the core zone of Bradford curves. These rare but influential entities shape scientific output in ways that Lotkaian informetrics cannot fully explain. To address this, the study introduces a new framework—non-Lotkaian informetrics—which better accounts for extreme values and extends beyond bibliometrics to areas like patent analysis and network studies.

5.2. Practical Applications

Despite bundled licensing agreements, Bradford’s law helps academic libraries evaluate their coverage rate in specific disciplines, offering a quantitative framework to estimate the relationship between acquired and available journals for informed decision-making.

• Journal Coverage Evaluation and Resource Optimization

Bradford’s law is a key tool for evaluating journal collection coverage and optimizing resources in academic libraries. By using the Bradford curve, libraries can assess how well their subscriptions cover a specific discipline, focusing on high-impact journals like those in SCI or Scopus while deprioritizing less essential ones such as predatory journals. This method helps libraries understand their coverage rate in particular fields, ensuring comprehensive collections by focusing on core journals and guiding effective weeding strategies to retire low-usage ones [46,47,48].

Additionally, Bradford’s law supports resource optimization by helping libraries periodically evaluate journal and database subscriptions, even when subscriptions appear stable year-to-year. This assessment not only aids in making informed decisions about renewals and expansions, ensuring continued coverage of core journals as disciplines evolve, but also provides valuable insights for developing and refining Collection Development Policies (CDP). By using Bradford’s law as a reference, libraries can better balance budget efficiency with evolving research needs, ensuring that their collections remain comprehensive and aligned with institutional priorities [49,50].

2.

Broader Applications in Academic Libraries

Bradford’s law provides a versatile framework for assessing academic impact by classifying papers based on citation counts and authors by their h-index. It ranks papers into zones of influence, similar to how top-tier journals are grouped in JCR rankings, helping researchers identify the most influential papers in their field [51]. The same principle applies to authors, where those with higher h-indices are placed in the core zone, offering insights into both productivity and influence [52]. This broader application not only enhances support for impactful research but also strengthens the academic ecosystem by providing valuable insights into scholarly influence.

6. Conclusions

This paper examines how integer constraints on the number of journals $T$ and articles $A$ affect the shape of Bradford curve, dividing it into two distinct zones: the core zone and the normal zone, based on the significance of these integer effects. Using the Simon–Yule model, we derive analytical results for key parameters and distributions under a constant entry rate. Theoretical formulas for each zone are developed, and the reasons behind the various shapes of Bradford curves are analyzed. Monte Carlo simulations are employed to study the impact of decreasing entry rates of new journals and aging rates of journals on the shape of the Bradford curve and the key parameters. Finally, we validate our proposed method using empirical data from the Croatian Chemistry and the Solar Power research datasets. The main conclusions are:

• Bradford curves should be divided into two separate zones based on the significance of integer constraints on journal and article numbers. Different formulas for each zone should be derived separately;
• Bradford curves can exhibit four different shapes, determined by the second derivatives of the core and the normal zones;
• The largest productivity $X_1$, the number of journals $T_0$, and the number of articles $A_0$ are key parameters influencing the shapes of Bradford curves. Decreasing entry rates and aging rates of journals affect these parameters.
• The proposed four-step method can predict general trends in Bradford curves despite some errors.
• Bradford’s law provides a valuable framework for academic libraries to evaluate journal coverage, optimize resource allocation, and refine Collection Development Policies (CDP), ensuring comprehensive and well-balanced collections as research needs evolve.