Universal patterns in passenger flight departure delays

Departure delays are a major cause of economic loss and inefficiency in the growing industry of passenger flights. A departure delay of a current flight is inevitably affected by the late arrival of the flight immediately preceding it with the same aircraft. We seek to understand the mechanisms of such propagated delays, and to obtain universal metrics by which to evaluate an airline’s operational effectiveness in delay alleviation. Here we use big data collected by the American Bureau of Transportation Statistics to design models of flight delays. Offering two dynamic models of delay propagation, we divided all carriers into two groups exhibiting a shifted power law or an exponentially truncated shifted power law delay distribution, revealing two universal delay propagation classes. Three model parameters, extracted directly from dual data mining, help characterize each airline’s operational efficiency in delay mitigation. Therefore, our modeling framework provides the crucially lacking evaluation indicators for airlines, potentially contributing to the mitigation of future departure delays.

direct insight into the operational efficiency of the airlines in mitigating and absorbing the potential propagation of departure delays.

Results
We collected data on flight delays from the BTS 20 , as detailed in Table S1 of the Supporting Material (SM), together, covering the activity of 14 US carriers over a period of 20 years. To quantify departure delays we calculate the difference l (minutes) between the actual versus the scheduled departure time of each flight, hence the greater is l the more severe is the delay. We characterize the delay patterns through the probability density p(l) to observe a delay of magnitude l, as shown in Fig. 1. We find that the 14 airlines condense around two distinct Groups, 1 and 2. To observe this we present in Figs. 3A and 4A the complementary cumulative distribution functions (CCDF) [27][28][29] .
l 0 ∫ > = ∞ which capture the probability to observe a delay in excess of duration l. Group 1 (Fig. 3A). In 6 of the airlines (AA, MQ, F9, DL, HA, AS) the CCDF is best captured by a shift power law (SPL) [30][31][32] of the form P(l 0 > l):(l + β) −α . Group 2 (Fig. 4A). For the remaining 8 airlines we observe an exponentially truncated shift power law (ETSPL) 33 , in which the tail exhibits an exponential cutoff. Below we show that these two classes represent two distinct models of propagated delays, providing a window into the airlines' operational efficiency in delay mitigation.
Decreased type of propagated delay. The effects of propagation are not readily discernible from the primary data 34 . On the one hand, given two successive delayed flights carried by the same aircraft, delay time of a current departure may be not merely attributed to the late arrival of the flight immediately preceding it, but also be attributed to one or more other factors, e.g., extreme weather on arrival airport. On the other hand, the role of non-propagation factors is likely diminished in the cases where the current departure delay is shorter than the one immediately preceding it with the same aircraft. This decrease of delay time expresses the efforts of the airline's operational staff in countering PD, hence characterizing the carrier's ability to absorb delays. Therefore, we concentrate below on decreasing type of propagated delays (DTPD), in which a propagated delay of an aircraft at j was shorter than its immediately preceding one at i.

Model 1 -Flight delay distributions with shift power law. The conception of a PD is rooted in an
immediately preceding departure delay, caused by any of the aforementioned reasons, from category (i) to (v), which then propagates to cause a current delay l with the same aircraft. In the case of DTDP, the initial delay can be of any length longer than l, here we seek the probability per delay interval near l, q(l), that a delay in the range (l, ∞) propagates to cause a delay of length l. Denoting by n B (l) the total number of DTPD instances per delay interval near l, we write l represents the total number of departure delays with a lag greater than l; the delay number density n(l) is the number of delays per delay interval near l. The transfer function q(l) captures the probability density for longer delays with number N(l) to generate a DTPD -hence providing the contribution of N(l) that transfers into n B (l). The reasoning behind Eq. (2) is analogous to the classic mean-field approach of statistical physics, where current delays are assumed to be affected by all previous delays, via the l-dependent transfer density function q(l). We obtain the density function q(l) directly from data mining, by comparing the ratio of DTPD number density n B (l) and the total number of delays N(l) with the delayed intervals longer than l. We find that the probability of a flight www.nature.com/scientificreports www.nature.com/scientificreports/ having a delay l by DTPD decreases with l ( Fig. 2), and therefore the observed transfer density functions q(l) can be effectively approximated by α + β = an inverse linear function, in which the parameters α and β represent airline -specific phenomenological parameters.
Next, we assume that DTPD represents a constant fraction k of all delays, allowing us to write  www.nature.com/scientificreports www.nature.com/scientificreports/ where we used (4) to express q(l). We, therefore, obtain where α 0 = α/k, and c is an integration constant. As a result, the analytical CCDF of departure delays is found to follow where N 0 is the yearly total number of flights of an airline, and the constant c 0 is defined as c 0 = c/N 0 (see Table S1 in the SM). This, indeed, represents the SPL 30-32 form of P(l) observed in Fig. 1a for the six airlines in group one.
To further examine the relevance of Model 1 above, in Fig. 3A we fit the CCDFs(red lines) from Eq. (9) to the real data from BTS (colored filled circles) for the six airlines (Group 1) in Fig. 3A (for Year 2014). At least in three decades, in the form of shift power-law, red lines fit well to empirical data which are independent of Model 1. As the result, a unique pair of parameters β and α 0 are obtained for each airline. Hence, while the form of the CCDF is universal, across all six airlines, the fitted parameters β and α 0 represent airline specific phenomenological values, whose meaning we discuss below. Results of the remaining 19 years of data are shown in Figs. S1-S10 of SM.
The relevance of our assumptions, as expressed in Eqs. (2)(3)(4)(5), can be reexamined by refitting the parameter β 2 , this time from the transfer density function q(l). This will provide us with two fitted evaluations of β: β 1 extracted above from the CCDF, i.e., integrated p(l) of Fig. 1, and β 2 extracted from Eq. (2). Their consistency, as well as any  2) and (10). To obtain the compensation interval m we fit the observed F(l) with the the integrated q(l) with compensation m (solid red lines). By appropriately tuning m we obtain the best fit in which the estimator β 2 approaches its CCDF-based counterpart β 2 . In all panels we set Δl = 5 min. (2020) 10:6890 | https://doi.org/10.1038/s41598-020-62871-6 www.nature.com/scientificreports www.nature.com/scientificreports/ discrepancies between these two fits, can gain us an additional validation and insight pertaining to our model analysis. Hence, in Fig. 3B we present the integrated function .. versus the shifted delays l + β 2 for all six airlines classified in model 1. Using Eq. (4) we write 2 2 allowing us to extract β 2 directly from the data, and thus yielding the desired independent evaluation of β in Eqs. (4) and (9). The obtained values for β 2 appear in Fig. 3B, showing similar trends to the β 1 values obtained from fitting Eq. (9) (red lines) to real data (colored filled circles) in Fig. 3A. Results of the remaining 19 years are shown in Figs. S30-S39 in the SM. These further support the relevance of our model and their analyses to the PD statistics. When we fit Eq. (9) to practical distributions CCDFs, i.e., P(l 0 > l) integrated from p(l) of Fig. 1a, we are calibrating the contributions from all five categories of delaying factors into DTPD assumed by Eq. (2), which yields parameters (α 1 , β 1 ) for each airline (shown in Table S2 of SM for year 2014). On the other hand, β 2 , is extracted from the transfer density function q(l) capturing only the propagation instances that had a decreasing effect with l, i.e. DTPD. Such decreases in delay time are likely the results of all efforts from airlines, airports, air controllers, and so on. to overcome the PD. Hence the ratio of these two estimators, β 2 /β 1 , can help us characterize an airline's operational effectiveness, quantifying how effectively that airline absorbs delays and attenuates their propagation.
Taken together, Model 1 is validated in Fig. 3 upon approximately three orders of magnitude, covering delay scales over several hundreds of minutes. Beyond that, all empirical curves (colored symbols) are truncated, deviating from the CCDF of Eq. (9) (red lines), indicating that the effects of non-propagation factors begin to dominate the delay dynamics. The empirically obtained F(l) curves saturate in the limit of large l, as q(l) approaches zero. In such cases, F(l) begins to deviate from the Eq. (10), as our approximate assumptions in Eq. (2) no longer hold.

Model 2 -Flight delay distributions with exponentially truncated shift power law.
We now turn to analyze the remaining 8 airlines (Group 2), whose CCDF features an exponential truncation in the limit l → ∞. Such behavior can arise from Eq. (7) if we assume the appropriate dependence of k on l, i.e. substitute the constant k by k(l). Such generalization of Model 1 relaxes the assumption that decreasing type of propagation delays (DTPD) occupy a constant fraction k of whole numbers of delays, and rather introduces an l-dependence, where, e.g., longer delays are less likely to be of the DTPD. Hence, we now replace k in Eq. (7) by which approaches a constant (1/h) in the limit of small l, but tends to zero as l → ∞. Model 2, therefore, converges to Model 1 in the limit where g is vanishingly small. As above, the parameters g and h represent phenomenological constants, extracted for each airline from the data. Using (11) in (7), we obtain  Performance metrics for delay mitigation. Our analysis allows us to extract three empirically accessible parameters with direct relation to an airline's treatment of PD.
Compensation parameter m. Consider the parameter β in (4). We can evaluate it through two independent empirical functions: the CCDF in (13), as done in Fig. 4A, yielding the estimator β 1 ; or, alternatively, the cumulative function F(l) in (10), providing the estimator β 2 (Fig. 4B). While both are aimed at evaluating the same parameter β, these two estimators are naturally inconsistent. In the CCDF, since we cannot distinguish the effects of these two kinds of delay factors 34 , all types of delays, including even cancelations 35,36 , are intermixed. In contrast, F(l), based on q(l) from Eq. (2), is constructed exclusively from DTPD. Therefore, we expect an inevitable discrepancy between the two estimates β 1 and β 2 . To correct for this discrepancy we re-evaluate F(l), to include the impact of increasing type of propagation delays (ITPD), thus bringing it closer to the observed CCDF. We achieve this by redefining n B (l) in Eq. (2) to include also contributions of ITPD, namely it now counts the number density of flights with delay around l, whose immediately preceding delay l′ ≥ l − m. The compensation parameter m allows us to include ITPD and impacts of non -propagation factors in our PD analysis. Therefore, in the limit m → 0, n B (l) is restricted to just DTPD, i.e. l′ ≥ l, reverting back to its original definition in Eq. (2). However, for a larger m, n B (l), and hence also q(l) and F(l), account for a growing contribution of ITPD. By tuning the parameter m we can now force the estimated β 2 to become consistent with β 1 , as we demonstrate in Fig. 4B. For example, for (2020) 10:6890 | https://doi.org/10.1038/s41598-020-62871-6 www.nature.com/scientificreports www.nature.com/scientificreports/ airline B6 we set m = 15 to reach β 2 = 8.70, a value that is within 0.3% of the estimated β 2 , which is evaluated at 8.73; Fig. 4A(a),B(a).
The parameter m characterizes the impact of ITPD together with that of DTPD in shaping the delay distribution P(l). Indeed, a small m indicates a dominant role of DTPD, while for a larger m, a greater factor of ITPD is introduced in order to reconcile F(l) with P(l 0 > l). Therefore, large m characterizes an airline in which ITPD plays a significant role, that is to say, a poor performance in mitigating PD. Consequently, m offers an inverse metric of the performance to evaluate an airline's success in absorbing delays and attenuating their propagating impact.
shift parameter β. A crucial quantity in the context of DTPD is the average flight delay absorption time L = l p − l c (minutes), capturing the level of decrease in delay time between the preceding departure delay (l p ) and the current propagated delay (l c ) with the same aircraft 37,38 . A successful delay mitigation is captured by a large L, i.e. a significant decrease in delay from preceding to current flight. We find, from the data, that L can be directly related to the shift parameter β 1 through an approximate linear relationship (Fig. 5A)  www.nature.com/scientificreports www.nature.com/scientificreports/ magnitude difference. Therefore, on average, the smaller β 1 is, the greater L is, and hence the more effective the airline in absorbing delays is. In Fig. 5A we exclude two airlines from the linear fit: in panel (a) the carrier HA is a distinctive outlier. This is rooted in its unique point-to-point operation style which is different from the common hub-and-spoke style 39 . Therefore, HA has intrinsically different propagation patterns. In panel (b) VX has a smaller flight volume compared to all other carriers (see Table S1 in SM) and it, therefore, falls below the mainstream statistical line. Equation (14) based on Fig. 5A show that airlines with smaller values β 1 are more effective in absorbing delays from the immediately preceding ones, which can help us to understand that the majority of flights with shorter delays (l < λ in Fig. 5B), where n B (l) in formula (1) are cumulated into histograms of N B1 (l) with each width set as λ 1 4 , where l is the mean value of l 1 and l 2 . Average absorption time L in Group 2 behaves more sensitive than those in Group 1, since L changes over about 10 minutes against the β 1 variation over the range of 50 minutes in Fig. 5A(a) for Group 1, while L changes over about 15 minutes against the β 1 variation over only about 11 minutes in Fig. 5A(b) for Group 2. The comparison is effective within each panel (a) and (b) of Fig. 5A, respectively. For instance, airline AA behaves more effective in absorbing the preceding delays than airline AS in Group 1; while EV behaves more effective in delay -absorption than WN in group 2, during 2014. Plots of L(β 1 ) for airlines operated in other 19 years are shown in Figs. S59-S69 of SM.
Critical delay λ. Now we turn to the role played by parameter λ. It acts as the critical delay of airlines in Group 2. CCDFs by Eq. (13) fit empirical data well because they are obtained with the aid of one more phenomenological function k(l). However, practical q(l) mined from primary data goes to zero not in the manner assumed by Eq. (4). Actually, empirical q(l) often touches zero beyond 240 minutes and fluctuates more or less above it (see Fig. 2). Therefore, cumulated q(l), i.e., colored filled circles obtained from data mining shown in Fig. 4B drastically deviate from the integrated transfer density function F(l) (red lines) from Eq. (4) starting at the critical delays (l = λ) (see panels in Fig. 4B). The exact fitting between analytical and empirical F(l) in the range l < λ indicates that the assumption of DTPD is successful; while deviations shown in Fig. 4B mean that F(l) integrated from q(l) becomes less effective as l exceeds λ, indicating a gradually increasing impacts of ITPD and non -propagation factors, since DTPD becomes less dominant as l increases. Therefore, λ serves as the critical delay separating DTPD dominated range from the range dominated by above mentioned two others. Moreover, the growing error effect starting at the deviations (l = λ) can be estimated with Fig. 5B. The cumulated flight number N B1 (l) often occupies the fraction larger than 90 percent as seen in it, manifesting successful degree of the assumption DTPD. Furthermore, we have intuitive results from the comparisons based on simple observations. Comparing panel(a) and (e) in Fig. 4B, also, taking the reference of values λ in counterparts of Fig. 4A, we see that red lines from Model 2 deviate more from real F(l) (filled colored circle) with larger m but similar values of λ; The same result was concluded for panel (d) and (f); While comparing panel (c) with (g),red line deviates more from real data for smaller λ but with almost equal values m; Comparing four central panels, ((b) and (f) with (c) and (g)), we see that red lines with both larger λ and smaller m from Model 2 fit better with real data of airlines in Group 2. Actually, the effect of the exponential factor in Eq. (13) becomes stronger as l grows beyond λ. The larger λ an airline has, the larger it has the dominant range of DTPD, which is the positive effect in delay mitigation. In this way, combining the observation of behaviours of λ and m in corresponding panels of both Fig. 4A,B, and the observation of the behaviors of them in other 19 years shown in Figs. S11-S29, and Figs. S40-S58, we capture an obvious tendency: the smaller m and the larger λ an airline has, the better for Model 2 to fit to the results mined from big data, and the better it operates, since DTPD governs larger range of delay interval (0, λ) supported by all efforts to counter the propagated delays.

Discussion and Conclusions
Discussion of techniques used. The mean-field approach in temporal regime has not been fully explored in previous works. Typically, the mean-field approach is used in the following setting: given a particle, the effect of interactions from all surrounding particles is represented instead by an equivalent external spatial field. In this paper, we use an analogous approach given in Eq. (9). Specifically, the delay (in l minutes) of a current flight is the result of delays transferred from all previous flights with delays longer than l. The phenomenological function q(l) acts to transfer intensity. In this way, we model the main stream tendency of air transportation as a DTPD. Note that both Model 1 and Model 2 have this main assumption. Moreover, an airline with a better model fit indicates a stronger tendency to be DTPD in terms of its operations, and thus its ability to mitigate departure delays. We observed this behavior through empirical flight delay data from the US.
We note that assumption 1 (Eqs. (2) and (3)) is only partially correct since we incorporated the impacts from four other categories (i-iv) as well as ITPD into that of the DTPD. However, we keep with this assumption because admits analytical derivations with good fits to real data. Moreover, we can distinguish airline-specific differences between β 1 and β 2 , as well as the emergence of metrics m and λ. Specifically, this allows us to check the precision of our models and validation ranges. Note that the CCDFs are model-free; they are obtained purely from data.

Conclusions.
Using the large scale data from BTS of US, we exposed two universal patterns of propagated departure delays of passenger flights, separating airlines into two distinctive groups -one with an SPL and the other with ETSPL CCDF. These groups exhibit two different mechanisms of delay propagation, as captured by Model 1 and Model 2. In the majority of cases the airlines remained consistent along this divide, however, few exceptions, in which airlines changed classification over the years have been observed (EV, AS, DL)(see SM) -a phenomenon that deserves further investigation.
Our analysis, investigating the relative roles of DTPD vs. ITPD, identifies three parameters that can help characterize an airline's operational performance in avoiding PD. The shift parameter β, which negatively correlates (2020) 10:6890 | https://doi.org/10.1038/s41598-020-62871-6 www.nature.com/scientificreports www.nature.com/scientificreports/ with the average delay absorption time L; the compensation interval m that quantifies the participation of ITPD as compared to DTPD; and the critical delay λ that separates the delay interval that is dominated by DTPD from that of ITPD and non-propagation factors. Together, they offer quantitative empirically accessible metrics for performance evaluation over time and across airlines.
Tracking the temporal behavior of these three metrics over time we can observe the impact of the US instituted Airline Passenger Protection rules (also known as Passenger Bill of Rights) in 2009, as well as their update in 2012 40 . These rules reduced the number of commercial aircraft delays in excess of 240 minutes, prompting the carriers to cancel such flights and reroute passengers. This is observed through the enhancement of operational performance exhibited by most airlines directly after 2008, and, again, after 2012 (Fig. S69 in SM). This not only indicates the positive impact of these US regulations instituted in 2009 and 2012, but also provides independent validation for the relevance of our proposed metrics, that were, indeed, able to detect this regulatory shift. Hence, we offer these three metrics as a basis for air transportation assessment in both the US and other countries.
In a broader perspective, we believe that our statistical physics inspired approach can advance our understanding of systems that are well beyond the realm of standard physical systems 41 . Indeed, focusing on statistical properties, and constructing simplified models, we were able to characterize a seemingly unpredictable and highly complex phenomenon, such as passenger flight delays, and expose its universally recurring and consistently predictable patterns.