Nearly every random vibration fatigue course teaches that the response of a single-degree-of-freedom system to broadband base excitation has Rayleigh-distributed peaks. The claim is repeated so often that it has hardened into doctrine. It is a good approximation for the largest peaks, which is why it has survived. But it is quantitatively wrong for the small peaks, and the error is not small. This post measures the discrepancy, shows why Dirlik’s one-exponential-plus-two-Rayleigh form captures it, and then asks the question that actually matters: does it change the answer? For classical S-N fatigue on a single mode the honest answer is barely — but that reprieve disappears the moment a second mode participates. For fracture mechanics crack growth the honest answer is more complicated, and more interesting, than the usual telling.
Where the Rayleigh Assumption Comes From
Consider a stationary Gaussian random process with zero mean and standard deviation $\sigma$. If the process is ideally narrowband — energy concentrated in a vanishingly thin band around a single frequency — then the time history looks like a sine wave with a slowly varying envelope. Every zero upcrossing is followed by exactly one maximum, and the amplitude of that maximum equals the envelope value. The envelope of a narrowband Gaussian process is Rayleigh distributed, and therefore so are the peaks:
$$p(\eta) \; = \; \eta \, e^{-\eta^2/2}, \qquad \eta = \frac{A}{\sigma}, \qquad \eta \ge 0$$
This gives the familiar reference values that show up in every peak-statistics report: the mean absolute peak is $\sqrt{\pi/2}\,\sigma = 1.2533\,\sigma$, and the RMS peak is $\sqrt{2}\,\sigma = 1.4142\,\sigma$. When a measured or synthesized time history returns a mean-to-sigma ratio noticeably below 1.2533, that is the first clue that the narrowband idealization has broken down.
The logic is sound. The problem is the premise. A lightly damped SDOF driven by broadband base input is not ideally narrowband, and the departure is larger than most engineers expect.
The Exact Result: Rice, Not Rayleigh
For any stationary Gaussian process, the distribution of maxima is known exactly. It was derived by S. O. Rice in 1945 and it depends on a single bandwidth parameter. Define the spectral moments of the one-sided response power spectral density $G(f)$:
$$m_n \; = \; \int_0^\infty f^{\,n} \, G(f) \, df$$
The rate of mean-level upcrossings, the rate of maxima, and the irregularity factor are then
$$n_0 = \sqrt{\frac{m_2}{m_0}}, \qquad E[P] = \sqrt{\frac{m_4}{m_2}}, \qquad \gamma = \frac{n_0}{E[P]} = \frac{m_2}{\sqrt{m_0 \, m_4}}$$
The irregularity factor $\gamma$ is the fraction of maxima that are also upcrossings. It runs from 0 for broadband white noise to 1 for an ideal narrowband process. Rice’s distribution of maxima is a weighted sum of a Gaussian term and a Rayleigh term:
$$p(\eta) \; = \; \frac{\sqrt{1-\gamma^2}}{\sqrt{2\pi}} \, \exp\!\left[ \frac{-\eta^2}{2(1-\gamma^2)} \right] \; + \; \gamma \, \eta \, e^{-\eta^2/2} \, \Phi\!\left( \frac{\gamma \eta}{\sqrt{1-\gamma^2}} \right)$$
where $\Phi$ is the standard normal cumulative distribution function. Setting $\gamma = 1$ kills the Gaussian term and recovers Rayleigh exactly. Setting $\gamma = 0$ leaves a pure Gaussian, which is correct for white noise where maxima occur everywhere including well below the mean. So the Rayleigh model is not merely an approximation — it is the $\gamma = 1$ endpoint of a one-parameter family, and the whole question is how far the real system sits from that endpoint.
The critical distinction. The Gaussian term in the Rice distribution has support on the negative axis. A finite fraction of the maxima of a real SDOF response lie below zero — they are local maxima on the downslope of a larger oscillation. The Rayleigh model assigns those events zero probability. They are not rare. They are not a rounding error. In the example below they account for roughly one maximum in ten.
Why a Real SDOF Is Never Narrowband
The relative displacement of a base-excited SDOF, which is the quantity proportional to spring stress, has the frequency response function
$$\frac{Z(\omega)}{\ddot{Y}(\omega)} \; = \; \frac{-1}{\left( \omega_n^2 – \omega^2 \right) + j \, 2 \zeta \omega_n \omega}$$
At frequencies well above resonance this rolls off as $1/\omega^2$, so the response PSD rolls off as $f^{-4}$. Now look at the fourth spectral moment. Its integrand is $f^4 G(f)$, which above resonance behaves as $f^4 \cdot f^{-4} = \text{constant}$. The integral does not converge. As the upper band limit of the base input is pushed higher, $m_4$ grows without bound, $E[P]$ grows without bound, and $\gamma$ marches toward zero.
This is not an academic curiosity. It means the number of peaks in an SDOF response is not a property of the SDOF at all — it is set by the bandwidth of the input. Two analysts using the same natural frequency and the same Q, but specifying the base input to 2000 Hz versus 10000 Hz, will count very different numbers of maxima and will get very different irregularity factors from the same physical structure.
The left panel makes the point. For Q = 10 and an input band limit of 500 Hz, the irregularity factor is 0.93 and the Rayleigh assumption is defensible. Extend the same input to 20000 Hz and $\gamma$ falls to 0.37. The right panel shows the Q dependence at a fixed 20 to 2000 Hz input band: even at Q = 60, $\gamma$ only reaches about 0.9. The narrowband ideal is approached slowly and never arrived at.
For absolute acceleration response the situation is worse still. That transfer function rolls off only as $1/\omega$, so the response PSD rolls off as $f^{-2}$ and even $m_2$ fails to converge. The upcrossing rate itself becomes unbounded. Whenever peak statistics are computed on an absolute acceleration signal rather than on a stress-proportional relative displacement, the numbers carry an additional bandwidth artifact on top of everything discussed here.
A Numerical Experiment
The remainder of this post uses a single controlled case. The base input is band-limited Gaussian white noise, flat at 0.10 G²/Hz from 20 to 2000 Hz, for an overall level of 14.07 GRMS. The SDOF has a 200 Hz natural frequency. The primary case is Q = 10, with Q = 25 and Q = 3 carried along for contrast. The record is 900 seconds at 8192 samples per second, and the response quantity is relative displacement. Peaks were extracted with no hysteresis gate, and cycles were counted with a self-contained ASTM E1049 three-point rainflow algorithm with the residual stack counted as half cycles.
One hundred milliseconds of the Q = 10 response is shown above, with maxima marked. The large oscillation at 200 Hz is the resonant response everyone expects. Riding on it are the small maxima marked in red — the direct, non-resonant transmission of the broadband input through the same structure. These are the low peaks. They are perfectly real, they are counted by any peak-picking or rainflow routine, and the Rayleigh model says they do not exist.
The measured distribution of maxima is compared with both models above. Rayleigh fails badly below about $1.5\sigma$ and is respectable above $2.5\sigma$. The Rice distribution, using the irregularity factor computed from the spectral moments with no fitting whatsoever, tracks the data across the full range including the negative-peak population highlighted in the log-scale panel. Note that a negative peak value in the Rice curve means: a local high point that happened while the signal was still below zero.
| Peak Statistic | Rayleigh | Q = 25 | Q = 10 | Q = 3 |
|---|---|---|---|---|
| Irregularity factor $\gamma$ | 1.000 | 0.894 | 0.785 | 0.574 |
| Maxima per second | 200.0 | 220.8 | 248.8 | 334.3 |
| Mean absolute peak / $\sigma$ | 1.2533 | 1.1563 | 1.0731 | 0.9446 |
| RMS peak / $\sigma$ | 1.4142 | 1.3477 | 1.2816 | 1.1641 |
| Fraction of maxima below zero | 0.0 | 0.047 | 0.098 | 0.202 |
| Fraction of maxima below $0.5\sigma$ | 0.118 | 0.195 | 0.268 | 0.394 |
| Largest peak / $\sigma$ | — | 4.80 | 4.83 | 5.00 |
Note the mean absolute peak column. Even at Q = 25, an unusually lightly damped structure, the ratio is 1.156 rather than 1.2533. At Q = 10 it is 1.073, a 14 percent shortfall. Field data is typically further off still, because real structures are multi-modal and real inputs are not flat. Measured aerospace buffet response records routinely return mean-to-sigma ratios in the 0.81 to 0.86 range, corresponding to irregularity factors far below anything an SDOF idealization would suggest.
From Peaks to Cycles: Dirlik’s Formula
A point worth making explicitly, because it is a frequent source of confusion: the Rice distribution describes peaks, while Dirlik’s formula describes rainflow ranges. They are different quantities answering different questions, but both exist for the same reason — the low-amplitude population that the narrowband ideal ignores.
Rainflow counting pairs maxima with minima according to the hysteresis-loop closure rules, so the resulting cycle set is not simply the peak set. Dirlik obtained an empirical closed-form expression for the rainflow range distribution by regression against a large set of Monte Carlo simulations. In terms of the normalized range $Z = S / (2\sqrt{m_0})$, which is simply amplitude divided by $\sigma$:
$$p(Z) \; = \; \frac{D_1}{\bar{Q}} e^{-Z/\bar{Q}} \; + \; \frac{D_2 \, Z}{R^2} e^{-Z^2/(2R^2)} \; + \; D_3 \, Z \, e^{-Z^2/2}$$
with coefficients determined entirely by the first four spectral moments:
$$x_m = \frac{m_1}{m_0}\sqrt{\frac{m_2}{m_4}}, \quad D_1 = \frac{2\left(x_m – \gamma^2\right)}{1+\gamma^2}, \quad R = \frac{\gamma – x_m – D_1^2}{1 – \gamma – D_1 + D_1^2}$$
$$D_2 = \frac{1 – \gamma – D_1 + D_1^2}{1 – R}, \quad D_3 = 1 – D_1 – D_2, \quad \bar{Q} = \frac{1.25\left(\gamma – D_3 – D_2 R\right)}{D_1}$$
The structure is exactly the one exponential plus two Rayleigh terms of the usual description. The third term, with coefficient $D_3$, is the classical narrowband Rayleigh contribution. The exponential term supplies the low-amplitude population. The second term interpolates between them. For the Q = 10 case the coefficients come out as $\gamma = 0.7846$, $x_m = 0.7650$, $D_1 = 0.1849$, $D_2 = 0.0528$, $D_3 = 0.7623$, $R = -0.2253$, and $\bar{Q} = 0.2311$.
The individual Dirlik terms are drawn as dotted curves in the left panel. The Rayleigh model overshoots the counted data through the entire $0.5\sigma$ to $2\sigma$ region and undershoots below $0.4\sigma$. Dirlik follows the data closely across four decades of probability density.
A caution on the Dirlik coefficients. Note that $R$ came out negative for this case. The formula still evaluates, since only $R^2$ appears in the second term and $D_2$ is small here, but a negative $R$ signals that the input spectrum sits near the edge of the parameter space over which Dirlik performed his regression. Dirlik’s expression is an empirical fit, not a derivation, and it should be treated as such whenever the response spectrum is strongly bimodal or has an unusual bandwidth. This limitation is exactly what motivated the Meta-Dirlik correction method.
Effect on Classical Fatigue Damage
Here is where the story requires honesty rather than advocacy. Take the standard Palmgren-Miner accumulation with an S-N curve of slope $b$, so that relative damage is proportional to $\sum n_i S_i^{\,b}$. Compare three cycle sets over the same 900 second record: the counted rainflow, the Dirlik prediction, and the narrowband Rayleigh assumption of $n_0 T$ cycles with Rayleigh-distributed amplitudes.
| S-N Slope $b$ | Dirlik / Counted | Narrowband Rayleigh / Counted |
|---|---|---|
| 3.0 | 1.003 | 1.026 |
| 4.0 | 1.012 | 1.038 |
| 6.4 | 1.041 | 1.068 |
| 9.0 | 1.082 | 1.110 |
The narrowband Rayleigh assumption is conservative by 3 to 11 percent depending on slope. That is all. Against the factor-of-two and factor-of-three scatter inherent in S-N data, a 7 percent bias at $b = 6.4$ is not the thing that will sink a fatigue assessment.
The reason is the $S^{\,b}$ weighting. With $b$ in the range 6 to 9, damage is overwhelmingly controlled by the high-amplitude tail, and in that tail the two models agree closely. They must agree: the rate of crossings of a high level is governed by $n_0$, which both models share. The low peaks that Rayleigh gets wrong are precisely the cycles that a steep S-N exponent renders irrelevant. Two errors of opposite sign — Rayleigh predicts fewer total cycles but assigns them larger amplitudes — also partially cancel.
So if the only question is classical high-cycle fatigue with a steep S-N slope, the Rayleigh assumption survives. That is why it has survived.
Extension to Multi-Degree-of-Freedom Systems
Everything above generalizes to real hardware without modification, because the linear response of an MDOF system to Gaussian input is still Gaussian. The recipe is unchanged: recover the stress PSD at the critical location, integrate its first four spectral moments, and evaluate $\gamma$, the Rice distribution, and the Dirlik coefficients from those moments. No new theory is required. What changes is the answer, and it changes far more than the SDOF results would lead you to expect.
Consider two modes at 60 Hz and 400 Hz, each with Q = 15, combined at equal RMS contribution and driven by the same 20 to 2000 Hz base input. The table below adds a third mode at 1200 Hz, with a single 60 Hz mode carried as the reference case.
| Case, Q = 15 per mode | $\gamma$ | Mean abs peak / $\sigma$ | Std dev of cycle means | Rayleigh / counted, $b = 6.4$ | Dirlik / counted, $b = 6.4$ |
|---|---|---|---|---|---|
| Single mode, 60 Hz | 0.648 | 0.983 | $0.343\,\sigma$ | 1.002 | 0.967 |
| Two modes, 60 + 400 Hz | 0.651 | 0.976 | $0.599\,\sigma$ | 1.515 | 0.906 |
| Three modes, 60 + 400 + 1200 Hz | 0.609 | 0.937 | $0.647\,\sigma$ | 1.727 | 0.891 |
The first row is the important control. A single mode at 60 Hz has an irregularity factor of only 0.648, well below the narrowband ideal, yet the Rayleigh damage estimate is correct to within 0.2 percent. The cancellation described in the previous section is doing all the work: Rayleigh predicts too few cycles and assigns them amplitudes that are too large, and for a unimodal spectrum those errors offset almost exactly. Add a second mode and the cancellation collapses. The Rayleigh estimate becomes non-conservative by 52 percent, and by 73 percent with three modes, while Dirlik stays within roughly 10 percent. The mechanism is straightforward: the large stress cycles are produced at the rate of the lowest participating mode, near 60 Hz, but the narrowband model manufactures full-amplitude cycles at the mean upcrossing rate of the combined response, which is 285 Hz for the two-mode case. It counts nearly five times too many large cycles. Note also that $\gamma$ alone does not warn you — the single-mode and two-mode cases have essentially the same irregularity factor and completely different answers. Bandwidth is not the same as bimodality.
Two consequences follow for the crack growth discussion below. First, the small-riding-on-large geometry that produces the low peaks is dramatically amplified when a high-frequency mode rides on a low-frequency one, and with it the cycle-mean spread: the standard deviation of rainflow cycle means grows from $0.241\sigma$ for the earlier SDOF case to $0.599\sigma$ here, and the fraction of cycles with $R > 0$ rises from 18.5 to 33.7 percent. Second, a bimodal or widely separated modal response is precisely where Dirlik’s empirical regression is weakest, which is why the bimodal-specific spectral methods of Jiao and Moan, Fu and Cebon, Sakai and Okamura, Benasciutti and Tovo, and Zhao and Baker exist at all, and what the Meta-Dirlik correction is intended to address. For an MDOF assessment the practical instruction is short: never assume the peak distribution, recover the stress PSD at the critical element from the random response solution, and compute $\gamma$ and the modal participation before choosing a method.
Why Crack Growth Is a Different Problem
Fracture mechanics changes the weighting. The Paris law is
$$\frac{da}{dN} \; = \; C \left( \Delta K \right)^{n}, \qquad \Delta K = \Delta S \, Y \sqrt{\pi a}$$
with $n$ typically between 3 and 4 for aluminum alloys. In terms of stress range, growth per cycle scales as $\Delta S^{\,n}$ — equivalent to an S-N slope of about 3.2 rather than 6.4. Halving the exponent moves the center of gravity of the calculation dramatically toward the mid- and low-amplitude cycles.
| Cycles with amplitude below | Fraction of cycle count | Fraction of $\Delta S^{3.2}$ (crack growth) | Fraction of $S^{6.4}$ (S-N) |
|---|---|---|---|
| $1.0\,\sigma$ | 0.513 | 0.032 | 0.001 |
| $1.5\,\sigma$ | 0.743 | 0.174 | 0.022 |
| $2.0\,\sigma$ | 0.896 | 0.438 | 0.127 |
| $2.5\,\sigma$ | 0.967 | 0.710 | 0.358 |
| $3.0\,\sigma$ | 0.992 | 0.888 | 0.638 |
Cycles below $2\sigma$ drive 44 percent of crack extension but only 13 percent of S-N damage. Below $1.5\sigma$ the split is 17 percent versus 2 percent. Crack growth draws roughly eight times as much of its driving from the amplitude band where the Rayleigh assumption is demonstrably wrong.
That said, a straight Paris integration with no threshold, run on all three cycle sets from the same 900 second record with identical initial and critical crack sizes, gives predicted lives within 3 percent of one another. The Rayleigh model returns 0.973 times the counted-rainflow life; Dirlik returns 0.997. The composition of the calculation is completely different, but the total comes out nearly the same, for the same cancellation reason as before.
The sensitivity appears once the model stops being a plain power law. Real crack growth analysis never is.
Three Mechanisms That Expose the Difference
1. The threshold is a moving amplitude filter. Below $\Delta K_{th}$ a cycle produces no growth at all. Since $\Delta K = \Delta S \, Y \sqrt{\pi a}$, the stress range required to reach threshold falls as $1/\sqrt{a}$. A cycle that is inert when the flaw is 0.010 inch becomes damaging once the crack reaches 0.100 inch. The threshold therefore sweeps down through the amplitude distribution over the life of the part, and it sweeps through exactly the region where the two models disagree. For the Q = 10 spectrum with a 2.5 ksi-in½ threshold and a 0.010 inch initial flaw, the sub-threshold fraction at $a_0$ is 46 percent of counted cycles but only 33 percent of Rayleigh cycles. By the time the crack reaches 0.400 inch, the counted spectrum still has 21 percent of its cycles below threshold while the Rayleigh spectrum has only 2 percent. The two models describe entirely different populations of inert cycles.
2. Mean stress and R-ratio. This is the mechanism most likely to bite in practice, and it is invisible in any amplitude-only comparison. Both the Rayleigh and Dirlik descriptions are amplitude distributions. They carry no information about cycle mean. If a spectrum is generated from them, every cycle is implicitly fully reversed at $R = -1$. Counted rainflow from an actual time history is nothing like that, even for a zero-mean process:
| Rainflow cycle-mean statistic, Q = 10, zero-mean process | Value |
|---|---|
| Standard deviation of cycle means | $0.241\,\sigma$ |
| Cycles with $\lvert \text{mean} \rvert > 0.25\,\sigma$ | 18.0 % |
| Cycles with $\lvert \text{mean} \rvert > 0.50\,\sigma$ | 5.1 % |
| Median R-ratio | −0.937 |
| Cycles with $R > 0$ (both ends same sign) | 18.5 % |
| Mean $\lvert \text{mean} \rvert$, small cycles (range $< 1\sigma$) | $0.284\,\sigma$ |
| Mean $\lvert \text{mean} \rvert$, large cycles (range $> 3\sigma$) | $0.108\,\sigma$ |
Nearly one cycle in five has both endpoints on the same side of zero, a condition that cannot occur under the fully-reversed idealization. And the correlation between range and absolute mean is negative: the small cycles are the ones carrying the large mean offsets, because they ride on the flanks of the large oscillations. This is not a coincidence — it is the same physics that produces the low peaks in the first place. NASGRO and every other modern crack growth code applies an R-ratio correction to $\Delta K_{th}$ and to the growth rate. Feed it a Rayleigh-derived spectrum and it sees $R = -1$ throughout. Feed it a counted rainflow spectrum from real response data and it sees a broad R distribution with a substantial positive-R population. The two runs will not agree, and the disagreement will be blamed on mean stress when its actual origin is the peak distribution.
3. Sequence effects. Retardation models such as Willenborg and Wheeler, and any crack-closure formulation, depend on the order in which cycles are applied. An overload creates a plastic zone that suppresses growth for many subsequent cycles. A rainflow histogram has thrown that information away. A spectrum synthesized from a Rayleigh or Dirlik amplitude distribution has never had it. Only an ordered turning-point history preserves it, which is a strong argument for handing the crack growth code the turning points and letting it do its own counting rather than delivering a pre-counted cycle table.
One more warning from the threshold study. When the stress level is low enough that growth is threshold-dominated, predicted life becomes extraordinarily sensitive to everything. Dropping the RMS stress from 2.0 to 1.5 ksi in the example case extended predicted life from 15 hours to 421 hours — a factor of 28 for a 25 percent stress reduction. In that regime the answer is controlled by a handful of extreme cycles, and neither Rayleigh, nor Dirlik, nor a 900-second sample of a random process estimates that tail with any reliability. The right response is not a better amplitude distribution. It is extreme value theory, a longer record, and a candid statement of uncertainty.
Practical Recommendations
Report the irregularity factor. It costs nothing, it is computed from spectral moments that are already available, and it immediately tells the reader whether the narrowband idealization is defensible. A value above roughly 0.95 justifies Rayleigh. A value of 0.8 does not, and a value of 0.2 means the analysis is in different territory altogether.
Report the mean absolute peak divided by sigma alongside it. Comparing that number to 1.2533 is the fastest sanity check available on any peak-statistics summary.
State the input bandwidth whenever a peak count or an irregularity factor is quoted. Because $m_4$ does not converge for this transfer function, those quantities are meaningless without it. Two analyses of the same hardware that differ only in the upper band limit will not be comparable.
For S-N fatigue on a single-mode-dominated response, use Dirlik or use Rayleigh; the difference is a few percent and Rayleigh errs conservative. For anything multi-modal, do not use Rayleigh. The 52 percent non-conservative error in the two-mode case above is not a refinement, and the irregularity factor by itself will not alert you to it — check the modal participation in the stress PSD, not just the bandwidth.
For crack growth, do not synthesize the spectrum from an amplitude distribution. Count the actual response time history, preserve the cycle means, and where retardation matters, deliver ordered turning points rather than a counted table. If a truncation threshold is applied to keep the spectrum file manageable, justify it by the fraction of damage discarded at the relevant exponent, not by the fraction of cycles — those two numbers differ by orders of magnitude, and the exponent that matters for crack growth is around 3, not 6.
Finally, when a crack growth prediction disagrees with expectation, check the R-ratio distribution before blaming the material data. The mean-stress content of a counted spectrum is a direct consequence of the same low-peak population discussed throughout this post, and it is the mechanism most often misattributed.
Summary
The Rayleigh peak distribution is the $\gamma = 1$ limiting case of the Rice distribution, and a real SDOF driven by broadband base input does not reach that limit. At Q = 10 with a 20 to 2000 Hz input, the irregularity factor is 0.78, one maximum in ten falls below zero, and the mean absolute peak is 14 percent below the Rayleigh value. Dirlik’s one exponential plus two Rayleigh terms reproduces the rainflow amplitude distribution well because the exponential term supplies exactly the low-amplitude population that the narrowband ideal omits.
For classical S-N fatigue on a single-mode-dominated response the practical consequence is modest, a few percent in the conservative direction, because a steep exponent discards the cycles in question. That reprieve does not extend to multi-modal systems: adding a second participating mode breaks the error cancellation that protects the unimodal case and leaves the Rayleigh estimate non-conservative by roughly 50 percent, with no change in the irregularity factor to warn of it. For fracture mechanics the consequence is larger in principle — cycles below $2\sigma$ carry 44 percent of crack extension versus 13 percent of S-N damage — but the plain Paris integration still lands within 3 percent. The real exposure comes from the features that make crack growth analysis what it is: a threshold that sweeps down through the amplitude distribution as the crack grows, an R-ratio dependence that amplitude-only models cannot represent at all, and sequence effects that a counted histogram has already destroyed. On each of those, the low peaks are not a detail. They are the whole question.
Free ebooks covering shock and vibration response spectra, stress-velocity, fatigue, and related topics are available at https://blog.vibrationdata.com/2025/11/27/toms-ebooks/