Abstract In this study we aim for a deeper understanding of the power-law slope, α, of waiting time distributions. Statistically independent events with linear behavior can be characterized by binomial, Gaussian, exponential, or Poissonian size distribution functions. In contrast, physical processes with nonlinear behavior exhibit spatiotemporal coherence (or memory) and “fat tails” in their size distributions that fit power-law-like functions, as a consequence of the time variability of the mean event rate, as demonstrated by means of Bayesian block decomposition in the work of Wheatland et al. In this study we conduct numerical simulations of waiting time distributions N(τ) in a large parameter space for various (polynomial, sinusoidal, Gaussian) event rate functions λ(t), parameterized with an exponent p that expresses the degree of the polynomial function λ(t) ∝ t p . We derive an analytical exact solution of the waiting time distribution function in terms of the incomplete gamma function, which is similar to a Pareto type II function and has a power-law slope of α = 2 + 1/p, in the asymptotic limit of large waiting times. Numerically simulated random distributions reproduce this theoretical prediction accurately. Numerical simulations in the nonlinear regime (p ≥ 2) predict power-law slopes in the range of 2.0 ≤ α ≤ 2.5. The self-organized criticality model yields a prediction of α = 2. Observations of solar flares and coronal mass ejections (over at least a half solar cycle) are found in the range of α obs ≈ 2.1–2.4. Deviations from strict power-law functions are expected due to the variability of the flare event rate λ(t), and deviations from theoretically predicted slope values α occur due to the Poissonian weighting bias of power-law fits.