# Understanding the bi-modal fatigue behavior of the case-hardened M50NiL steel

Anicleson Raony Alves de Moura, Franck Morel, Etienne Pessard, Daniel

Bellett, Louis Augustins, Damien Herisson

► To cite this version:

Anicleson Raony Alves de Moura, Franck Morel, Etienne Pessard, Daniel Bellett, Louis Augustins, et al.. Understanding the bi-modal fatigue behavior of the case-hardened M50NiL steel. International Journal of Fatigue, 2026, 207, pp.109495. 10.1016/j.ijfatigue.2026.109495. hal-05465795

HAL Id: hal-05465795

https://hal.science/hal-05465795v1

Submitted on 19 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Understanding the bi-modal fatigue behavior of the case-hardened M50NiL steel

Aniclelsion Raony Alves de Moura<sup>a,b,1,2,*</sup>, Franck Morel<sup>a</sup>, Etienne Pessard<sup>a,1,2</sup>, Daniel Bellett<sup>a,1,2</sup>, Louis Augustins<sup>b</sup>, Damien Herisson<sup>b</sup>

## ARTICLE INFO

**Keywords:**

Case-hardened steel

Material property gradients

Bi-modal fatigue behavior

Crack initiation

Probabilistic model

HCF and VHCF

Gears

Notched specimens

## ABSTRACT

The gears used in aircraft engines are typically made from high-strength steels reinforced by thermochemical treatments (TCT). These treatments increase surface fatigue strength through microstructural modifications, enhancing hardness and adding compressive residual stresses. In some cases, the combination of material, TCT, and applied stress can lead to a bi-modal fatigue behavior, notably in failures at the gear tooth root. This work investigates the bi-modal fatigue response of M50NiL case-hardened steel by characterizing and analyzing crack initiation mechanisms to propose a relevant fatigue modeling approach. A comprehensive experimental fatigue test campaign was carried out on notched specimens under plane bending and on gear specimens using a Single Tooth Bending Fatigue (STBF) method. The resulting Wöhler diagram shows significant scatter in fatigue life for several stress levels, suggesting a bi-modal behavior with two distinct populations. Fractographic analyses confirmed the competition between two different crack initiation mechanisms depending on stress level and number of cycles to failure. A statistical analysis using a mixture model also indicates that a bi-modal distribution best represents the results. Accordingly, a probabilistic model is proposed to describe the bi-modal fatigue behavior from a global perspective, based on the maximum applied hot-spot surface stress for a fixed stress ratio. Finally, a complementary local stress analysis shows that the combined effect of stress and material property distributions significantly influences local maximum stress variation. Correcting for these factors reduces scatter in the bi-modal stress levels.

## 1. Introduction

Gears are critical mechanical components in power transmission systems, directly influencing the capacity, reliability, and service life of the system [1,2]. The gears used in aircraft engines are typically made from high-performance steels reinforced by thermochemical treatments (TCT), such as carburizing and nitriding [3]. These treatments increase the local fatigue strength in the surface layer through microstructural modifications, enhancing the surface hardness and adding a compressive residual stress gradient [4].

The most common fatigue failure modes in gears are due to repeated contact between flanks (such as wear, scuffing, pitting and contact fatigue), and tooth root bending fatigue [5]. This latter failure mode is one of the most frequent, accounting for over 30% of observed failures, and often leads to complete system breakdown [6-8]. In this case, the fatigue life of the teeth is greatly influenced by the properties of the surface layers, such as hardness, residual stress, microstructure, and roughness.

Many testing methodologies are available to evaluate the tooth root fatigue strength, including running gear tests, back-to-back tests, single tooth bending fatigue (STBF) tests, and notched specimen tests [9]. STBF tests are among the most widely used due to their efficiency: multiple tests can be performed on the same gear, costs are lower compared to full-gear testing, and they can be carried out using standard universal testing machines [5,10]. When carefully designed to replicate the stress state at the tooth root, notched specimens offer additional advantages, such as higher test frequencies, lower cost, and improved compatibility with instrumentation [11].

In some cases, the relationship between the applied stress distribution, the material property distributions close to the surface, in particular the hardness and residual stress profiles, the inclusionary properties and the microstructure can lead to a bi-modal behavior in terms of fatigue life [9]. It is believed that this behavior occurs due to competition between different crack initiation mechanisms,

<footnote><sup>a</sup>Arts et Métiers Institute of Technology, LAMPA, 2 Bd du Ronceray, 49035 Angers, France</footnote>

<footnote><sup>b</sup>SafranTech, Materials and Processes Department, Rue des Jeunes Bois, Chateaufort, 78772 Magny les Hameaux, France</footnote>

<footnote>* Corresponding author at: Arts et Métiers Institute of Technology, LAMPA, 2 Bd du Ronceray, 49035 Angers, France. E-mail address: aniclelson.alves_de_moura@ensam.eu (A.R. Alves de Moura).</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Micrography of a tooth of a carburized gear (etched with a 4% nitral solution).</figcaption></figure>

characterized by very different crack initiation kinetics. In such scenarios, classical fatigue life prediction models are often inadequate to accurately describe the material behavior.

Several studies in the literature have investigated the fatigue behavior of materials and loading conditions that lead to significant variability in fatigue life. In steels subjected to rotating bending, such as the SUJ2 steel, Sakai et al. [12] and Goto et al. [13] observed the formation of duplex S-N curves, associated with the transition from surface crack initiation (at high stress levels) to internal initiation (at lower stress levels). They noted that the phenomenon is strongly influenced by surface property gradients. In tension-compression fatigue tests, authors such as Bayraktar et al. [14], Akiniwa et al. [15], Furuya [16] Sakai [17], and Li et al. [18] also reported this type of competition, where the presence or absence of residual stresses significantly affects the dominant initiation mechanism. The role of inclusions was emphasized by Yang et al. [19], who showed that internal fatigue crack nucleation is more frequent in steels with high inclusion content, while in steels of greater purity, crack initiation is more likely to occur at microstructural inhomogeneities. Similar findings were reported by Li et al. [20] and Zhang et al. [21] in high strength steels, where inclusions or local microstructural anomalies proved to be critical in determining the transition between mechanisms.

This bi-modal behavior in the fatigue life has also been observed in non-ferrous materials. In aluminum alloys [22], nickel-based superalloys [23–25], and titanium alloys [26–28], the coexistence of surface and internal initiation mechanisms was reported across different stress ranges, depending on microstructure, loading type, and the presence of inclusions, highlighting the multivariate nature of the phenomenon. For gears, similar behavior has been identified in carburized steels subjected to tooth root bending fatigue, such as the 18CrNiMo7-6 steel [2] and the 16NiCrMo13 steel [11]. These findings reinforce the need for fatigue models that account for multiple crack initiation mechanisms, in order to accurately predict bi-modal fatigue life behavior, particularly in the very high cycle fatigue (VHCF) regime.

This work aims to investigate the presence of bi-modal fatigue behavior in gears made of case-hardened M50NiL steel, focusing especially on the origins of the high scatter in the fatigue life and the impact of the residual stress gradient. To achieve this, single tooth bending fatigue (STBF) tests and plane bending fatigue tests on notched specimens were conducted. In previous work [11], the authors demonstrated that the notched specimens are representative of the gear tooth fatigue behavior. The crack initiation mechanisms were characterized and analyzed in order to propose a relevant fatigue modeling approach. A statistical analysis of the experimental data is also performed to verify the presence of bi-modal behavior from a probabilistic standpoint. A novel probabilistic approach is introduced to interpret the macroscopic fatigue response of the M50NiL case-hardened steel, taking into account its bi-modal fatigue behavior, based on the dominant mechanisms involved in crack initiation. For the sake of clarity, the probabilistic analysis is initially presented using the applied hot-spot

<figure><img src="image_4.png" /><figcaption>Fig. 2. Hardness profile in the critical fatigue zone at the tooth root.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Residual stress profile in the critical fatigue zone at the tooth root.</figcaption></figure>

<table><caption>Table 1
Chemical composition of M50NiL steel (wt%) [30,33,34].</caption><thead><tr><th></th><th>C</th><th>Cr</th><th>Mo</th><th>Ni</th><th>V</th><th>Mn</th><th>Si</th><th>Fe</th></tr></thead><tbody><tr><td>Minimum</td><td>0.11</td><td>4.00</td><td>4.00</td><td>3.20</td><td>1.13</td><td>0.15</td><td>0.10</td><td>Base</td></tr><tr><td>Maximum</td><td>0.15</td><td>4.25</td><td>4.50</td><td>3.60</td><td>1.33</td><td>0.35</td><td>0.25</td><td>Base</td></tr></tbody></table>

stress at the notch tip. In the final discussion, a local perspective is adopted, incorporating the influence of stress distributions on both lifetime and stress level scatter, as well as on crack initiation sites. This approach enables a local assessment of how these distributions affect the bi-modal behavior.

## 2. Material and experimental procedures

### 2.1. Material and thermochemical treatment

The material investigated in this study is M50NiL steel, a high-strength, martensitic alloy characterized by a high molybdenum content. It is a modified version of M50 steel, distinguished by a higher nickel content and lower carbon content [29,30]. This steel is commonly used in the aerospace industry for the manufacturing of engine components such as bearings and gears [31]. In the quenched and tempered state, it exhibits a yield strength of approximately 1200 MPa, an ultimate tensile strength (UTS) of about 1400 MPa, and a hardness of approximately 450 HV [32]. The chemical composition of the material is provided in Table 1.

The specimens were machined from forged bars, austenitized, and then subjected to low-pressure carburizing in a carbon-rich atmosphere, which promoted surface carbon absorption and diffusion into the subsurface. The carburizing process was carried out at 950 °C for approximately 12 h to achieve a Case Hardening Depth (CHD) of 1 mm. This treatment produces a carbon concentration gradient from the surface toward the core. The dark layer that can be seen in Fig. 1 near the surface represents the treated zone, which fades toward the core, indicating the treatment depth. Subsequently, the specimens underwent gas
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Sketch of a gear tooth and the corresponding notched specimen. The red arrow shows the position and direction where residual stress and hardness profiles were measured. 0 refers to the surface/origin, and x represents the depth variation.</figcaption></figure>

quenching at 1095 °C to form a martensitic microstructure, followed by a cryogenic treatment at -75°C for 2 h to reduce retained austenite. Finally, tempering was performed at 535 °C for three cycles of 2 h each to relieve internal stresses. This TCT sequence produced a hardness gradient, with surface hardness reaching approximately 730 HV while the core remained at 450 HV. A compressive residual stress gradient is also observed, with a maximum value at the surface. The residual stress and hardness profiles are shown in Figs. 2 and 3, respectively. The origin (0) and direction along which these profiles were measured are indicated by the red arrow in Fig. 4. The residual stresses are approximately biaxial, but only the component perpendicular to the measurement direction is shown in Fig. 3, as it is the direction most relevant to the applied loading.

## 2.2. Specimen geometries

Two specimen types have been tested in fatigue: gear specimens and notched specimens. The notched specimens were designed to replicate the stress conditions at the tooth root of the gear specimens, allowing testing at higher frequencies. This geometry (Fig. 4) has already been experimentally and numerically validated in previous studies, demonstrating equivalent mechanical properties, fatigue behavior, and crack initiation mechanisms when compared to the gear specimens [9, 11]. The geometric characteristics of the gear specimens are listed in Table 2.

To ensure identical metallurgical properties, the notched specimens were machined from the same M50NiL steel bar as the gear specimens, and identical thermochemical treatment parameters were applied. After the TCT process, the critical fatigue zones of the specimens were ground to a depth of 300 µm in order to obtain at the surface an arithmetic mean height roughness (Sa) of approximately 0.4 µm.

The gears specimens received an additional surface treatment consisting of shot-blasting (ASR130, 125% media coverage) followed by tribofinishing (i.e. a combination of chemical and mechanical polishing). This combined treatment generated a compressive residual stress peak at the surface of approximately -800 MPa and a Ra of 0.4 µm. Finally, to accommodate the load limitations of the fatigue test machines, the notched specimen width was set to 6 mm, whereas the gear tooth width was 10 mm. Given the results presented in the following section, this difference in thickness is considered to have a negligible size effect on the fatigue behavior.

<table><caption>**Table 2**
Main characteristics of gear specimens.</caption><thead><tr><td></td><td>Symbol</td><td>Unit</td><td>Value</td></tr></thead><tbody><tr><td>Module</td><td>m</td><td>mm</td><td>4.5</td></tr><tr><td>Number of teeth</td><td>z</td><td>–</td><td>24</td></tr><tr><td>Face width</td><td>b</td><td>mm</td><td>10</td></tr><tr><td>Pressure angle</td><td>α</td><td>deg</td><td>20</td></tr><tr><td>Helix angle</td><td>β</td><td>deg</td><td>0</td></tr><tr><td>Root fillet</td><td>r</td><td>mm</td><td>1.71</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Gear specimen mounted in the STBF bench [11].</figcaption></figure>

## 2.3. Fatigue testing

The gear tests were carried out on a Single Tooth Bending Fatigue (STBF) bench integrated within a standard MTS 100 kN servo-hydraulic fatigue testing machine, as shown in Fig. 5. A sinusoidal load was applied at 40 Hz with a constant stress ratio $R_{ref} = 0.05$. The tests were stopped if the tooth did not fail before $2 \times 10^7$ cycles or after a stiffness drop corresponding to a drop of 0.02 mm in the minimum actuator position. Six evenly spaced teeth were removed from each gear to accommodate the testing bench, allowing up to five fatigue tests per gear. This modification has no significant effect on the stress distribution at the tooth root adjacent to the removed tooth.

The notched specimens were tested in plane bending, using RUMUL Cracktronic resonant testing machines (Fig. 6) operating at approximately 120 Hz, also with a stress ratio $R_{ref} = 0.05$. Crack initiation was detected by monitoring the drop in the natural frequency ($\Delta f$) resulting from the decrease in specimen stiffness due the presence of a crack. The stopping criterion was defined as a frequency drop of 0.5 Hz or after $10^8$ cycles without specimen failure. For both specimen types, the stopping criteria correspond to the presence of a crack which is approximately 1 mm deep. These two test benches are the same as those used in previous work [11]. Compared to the servo-hydraulic machine, the resonant machine is simpler and more versatile, allowing fatigue testing at negative stress ratios.

# 3. Results

## 3.1. STBF and plane bending fatigue tests

A total of 74 fatigue tests were conducted, comprising 32 plane bending tests and 42 STBF tests. Among the plane bending specimens, 4 were subjected to the same additional surface treatments as the gear specimens (i.e. shot-blasting and tribofinishing). The Table 3 summarizes the surface conditions applied to all tested specimens, listed in the sequence of their application, and highlights the differences between the two configurations.
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Notched specimen mounted in the Plane bending bench.</figcaption></figure>

<table><caption>**Table 3**

Surface treatments applied to the tested specimens.</caption><thead><tr><th>Specimen</th><th>Qty</th><th>TCT</th><th>Grinding</th><th>Shot-blasting</th><th>tribofinishing</th></tr></thead><tbody><tr><td>Gear</td><td>42</td><td>Carburizing</td><td>300 µm</td><td>Yes</td><td>Yes</td></tr><tr><td rowspan="2">Notched</td><td>24</td><td>Carburizing</td><td>300 µm</td><td>—</td><td>—</td></tr><tr><td>4</td><td>Carburizing</td><td>300 µm</td><td>Yes</td><td>Yes</td></tr></tbody></table>

<table><caption>**Table 4**

Number of tests performed at each stress level for the two specimen geometries.</caption><thead><tr><th rowspan="2">σ<sub>max</sub> [MPa]</th><th colspan="2">Gear specimens</th><th colspan="2">Notched specimens</th><th rowspan="2">Total</th></tr><tr><th>Failed</th><th>Run-out</th><th>Failed</th><th>Run-out</th></tr></thead><tbody><tr><td>2400</td><td>4</td><td>—</td><td>—</td><td>—</td><td>4</td></tr><tr><td>2300</td><td>—</td><td>—</td><td>2</td><td>—</td><td>2</td></tr><tr><td>2200</td><td>10</td><td>—</td><td>1</td><td>—</td><td>11</td></tr><tr><td>2100</td><td>—</td><td>—</td><td>4</td><td>—</td><td>4</td></tr><tr><td>2000</td><td>10</td><td>—</td><td>—</td><td>—</td><td>10</td></tr><tr><td>1900</td><td>—</td><td>—</td><td>13</td><td>—</td><td>13</td></tr><tr><td>1800</td><td>5</td><td>5</td><td>5</td><td>—</td><td>15</td></tr><tr><td>1700</td><td>—</td><td>6</td><td>1</td><td>3</td><td>10</td></tr><tr><td>1600</td><td>—</td><td>2</td><td>—</td><td>2</td><td>4</td></tr><tr><td>1500</td><td>—</td><td>—</td><td>—</td><td>1</td><td>1</td></tr></tbody></table>

The first specimen was tested using the step-loading method, starting at a maximum applied hot-spot stress of $\sigma_{max} = 1500$ MPa. The value of $\sigma_{max}$ at the notch surface (position $x = 0$ in the direction indicated in Fig. 4) was determined following the method described in [9]. This procedure is based on the distribution of the first principal stress ($\sigma_1$) along a path normal to the hot-spot surface for a given applied force, which is usually estimated at the HPSTC<fnref n="1" /> according to the standard ISO 6336-3 [35]. If no failure occurred after $10^7$ cycles, the maximum stress was incrementally increased by 100 MPa. The stress levels and the number of specimens tested at each level were progressively adjusted based on the results of previous tests. This adaptive approach was adopted to better capture the potential presence of bi-modal fatigue behavior and to construct a representative S-N diagram (Wöhler curve) for the material. Testing order was randomized to eliminate fabrication bias. The number of specimens tested at each stress level is presented in Table 4.

The results of the fatigue tests are presented in Fig. 7 in the form of an S-N curve, showing the maximum applied hot-spot stress at the

surface of the tooth or notch root, as a function of the number of cycles to failure.

The fatigue resistance exhibited by the specimens is remarkably high, with no failures occurring even at stress levels below 1700 MPa, a value significantly higher than the UTS of the core material - 1400 MPa. This underscores the improvement in the mechanical performance provided by the thermochemical treatment.

It is important to note that very large scatter in the fatigue life is observed for different stress levels. For example, at 1900 MPa, the fatigue life ranges from approximately $9 \times 10^3$ to $3 \times 10^7$ cycles. The origin of this dispersion will be discussed later, after the fractographic analyses of the fatigue failure surfaces have been presented.

Two distinct S-N data-sets can be seen in Fig. 7 and they can be identified based on the variability in the fatigue life. The first group corresponds to specimens with a low number of cycles to failure, typically less than $2 \times 10^5$ cycles, while the second group includes specimens with a high number of cycles to failure, generally greater than $5 \times 10^5$ cycles, along with the run-out specimens.

Additionally, it can be observed that, in the short-life regime, the STBF data-points (blue points) tend to exhibit slightly longer fatigue lives compared to the notched specimens (black points). In contrast, in the long-life regime, both the STBF and specimen data-points appear to display similar fatigue lives. This difference for short lives may be attributed to the additional surface treatments applied to the gear specimens (i.e. shotblasting and tribofinishing), which are known to enhance the surface residual compressive stresses. Thus, improving the fatigue resistance as least when crack initiation occurs at the surface.

To confirm the influence of the additional treatments, four notched specimens were subjected to the same additional treatments (i.e. shotblasting and tribofinishing) and subsequently tested in fatigue (red points in Fig. 7). As can be seen in the S-N diagram, the behavior of these four notched specimens is closer to that of the gear specimens than to that of the untreated notched specimens. This further demonstrates that the notched specimens are representative of the STBF tests in terms of both fatigue strength and lifetime variability.

## 3.2. Initiation mechanisms

After failure, in order to observe the fatigue failure surfaces, the specimens were completely broken in half by immersing it in liquid nitrogen and then hitting it with a hammer to induce brittle fracture. The fractographic analyses of the failure surfaces was conducted using a Scanning Electron Microscope (SEM). All specimens were observed, the carburized layer and the core material are easily distinguishable, as well as the crack initiation sites, as shown in Fig. 8.

Both the gear specimens and notched specimens exhibit similar crack initiation mechanisms, which have been classified into two main categories: surface initiation or internal initiation.

### 3.2.1. Surface initiation

This failure mode is characterized by crack initiation at surface defects, such as scratches (Fig. 9(a)) and corrosion pits (Fig. 10c), followed by transgranular propagation. In these cases, transgranular propagation is identified by fatigue striations converging toward the surface along the semi-elliptical crack front. All striations converge to a single point, as illustrated in Fig. 10(a).

A total of 11 notched specimens and 24 gear teeth exhibited failures initiated at the surface, across various stress levels. The surface defects exhibited $\sqrt{\text{area}}$ values between 3 and 8 $\mu$m, while the depth of the observed scratches ranged from 0.5 to 3 $\mu$m. Despite the improved fatigue life in the short-life regime, the shotblasting and tribofinishing treatments did not significantly affect the nature or size of the defects, which were similar in all specimen types.

<footnote n="1">¹ The Highest Point of Single-Tooth Contact is the highest point on the line of action where only one pair of teeth are engaged.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** STBF and plane bending fatigue specimen results shown on a Wöhler diagram: maximum applied hot-spot stress at the surface vs. number of cycles for a stress ratio of $R_{ref} = 0.05$. (SB + TF) = shotblasted and tribofinished.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Fatigue failure surfaces for (a) a gear tooth and (b) a notched specimen.</figcaption></figure>

### 3.2.2. Internal initiation

In many cases crack initiation occurred below the surface, associated with the presence of non-metallic inclusions, resulting in the formation of characteristic fish-eye morphologies. Chemical analyses of the inclusions indicated they consisted of aluminum, sulfur, and calcium compounds ($Al_2O_3$ + CaS). In all observed fish-eye features, no Fine Granular Area (FGA) was detected, or it was extremely small, as illustrated in Figs. 10 and 9(b).b.

This mechanism was identified in 5 gear specimens and 12 notched specimens. The depth of the inclusions ranged from 55 to 150 µm in the gears, and from 20 to 65 µm in the notched specimens. In both cases, the minimum initiation depth was located in regions where the compressive residual stress gradient began to stabilize, at a value of approximately -300 MPa, as shown in Fig. 11. At these depths, the hardness remains close to surface values.

The depth at which this stabilization occurs varies between the two specimen types due to differences in surface treatments, and consequently, residual stress profiles. The observed difference of approximately 35 µm in initiation depth is directly related to the intensity of the compressive residual stress peak.

In both specimen types, the inclusion sizes, expressed in terms of the $\sqrt{area}$ parameter, ranged from 2 to 10 µm. By comparing the

inclusion depth, inclusion size, and the corresponding fish-eye size, as shown in Fig. 12, a linear correlation between the three parameters was identified. The fish-eye size is defined by the $\sqrt{area}$ of the nearly circular region surrounding the inclusion, as illustrated in Fig. 9(b) (the fish-eye is enclosed within the yellow circle). The fish-eye always touches the specimen surface, indicating that the inclusion depth corresponds to approximately the fish-eye radius. Deeper initiation sites and larger inclusions result in larger fish-eyes. Although the same trend was observed in both specimen types, a clear shift was noted between the two data sets. Assuming that both follow linear relations with similar slopes, the shift between them can be estimated at approximately 35 µm, which corresponds to the observed difference in the initiation depth.

These two crack initiation mechanisms (surface or internal) are strongly correlated with the observed lifetime variability, as shown in Fig. 10. This is further evidence of a bi-modal fatigue behavior. The only exception was one notched specimen, which exhibited surface initiation at a maximum stress of $\sigma_{max} = 1800$ MPa and failed after 14,258,394 cycles.

Despite the different lifetimes and crack initiation mechanisms, the range of defect sizes that lead to failure in terms of $\sqrt{area}$ are approximately the same for both specimen types.
<figure><img src="image_11.png" /><figcaption>Fig. 9. Fatigue failure surfaces of notched specimens.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Wöhler diagram of M50NiL showing the identified initiation mechanisms. (a) Surface initiation site, in transgranular mode (specimen M15-146 175 cycles at 1700 MPa), (b) Surface initiation site, in transgranular mode on a pit (specimen G115-08-44 390 cycles at 2000 MPa), (c) Internal initiation site, on an inclusion (specimen G116-18-4 263 394 cycles at 2000 MPa) and (d) Surface initiation site, in transgranular mode on a pit (specimen G127-14-5 926 704 cycles at 1800 MPa).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 11. Depth of fish eyes as a function of the evolution of the residual stress profiles. To facilitate data visualization, the gear points were shifted by -50 MPa.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Relationship between inclusion depths, size and fish eye area.</figcaption></figure>

This competition between crack initiation mechanisms in case hardened M50NiL gears in terms of bending fatigue can also be observed in the work of Lisle [36]. In this study, the lifetime prediction was presented with the endurance limit defined by a horizontal line. A representation that the authors themselves describe as valid "if one believes the endless mortality of a steel component under infinite cycles".

When a bi-modal behavior is observed, classical fatigue models and traditional S-N curves, such as those proposed by Basquin and Stromeyer, become inadequate to accurately capture the material's fatigue response. A dedicated modeling approach, specifically adapted to account for these effects, is therefore required.

## 4. Statistical analysis of fatigue data

When the results of a fatigue study lead to the formation of a duplex S-N curve, traditional modeling methods become inadequate. Typically, the lifetime (or the stress) populations are modeled using a single statistical distribution. However, in the presence of multiple subpopulations, modeling based on a single distribution tends to overestimate the scatter.

Similarly, propagation-based models are inadequate for simultaneously capturing the kinetics associated with different crack-initiation mechanisms. In addition, they do not account for the evolution and

variability of local material properties, such as hardness and residual stress, during crack growth.

In such cases, probabilistic models can be employed to better describe the observed behavior. In statistics, the term “bi-modality” refers to populations composed of two subpopulations, each following a distinct probability distribution. The most frequent value of a random variable is referred to as the mode, thus, a bi-modal population exhibits a primary mode and a secondary mode, each corresponding to one of the two subgroups.

In the literature, two main approaches are identified for modeling datasets that exhibit multimodality: methods that combine several probability laws (such as the mixture model, competing risks model, multiplicative model, segmented linear model, and nested model) or bi-modal distributions.

In a previous study, Argoud et al. [9] conducted a comprehensive analysis of different methods available in the literature to determine the most appropriate approach for modeling bi-modal fatigue behavior. The study concluded that a mixture model provides one of the most suitable frameworks for capturing bi-modality. Several examples of the application of mixture models to fatigue data are available for various metallic materials [37–41]. While the log-normal distribution is most commonly used, Sakai et al. [12], Harlow [42] employed three-parameter and two-parameter Weibull distributions, respectively, to model each sub-population. Therefore, in this work, mixture models based on both the log-normal and Weibull distributions, will be investigated.

The S-N data and the fractographic analyses clearly indicate the formation of two sub-populations. However, statistical bi-modality has not yet been demonstrated, which is a prerequisite for the application of a mixture model.

In the following sections, the modeling procedure and results for a stress level of $\sigma_{max} = 1900$ MPa for the notched specimens and 2000 MPa for the gear specimens are presented. These stress levels were selected as examples due to the larger number of available results. In both cases, the analyzed stress corresponds to the maximum applied stress at the surface of the notched root or at the gear tooth root. The data is analyzed using the cumulative distribution function (CDF) of the fatigue life $N_f$. This function, $F(N)$, associates each number of cycles $N$ with the probability that $N_f \le N$. Thus, $F(N)$ is defined as $F(N) = P(N_f \le N)$.

### 4.1. Empirical cumulative distribution function

Initially, the empirical cumulative distribution function (CDF) $\hat{F}(N)$ was calculated using Eq. (1). This formulation, also known as Chegodyev's formula, is considered less biased than the classical Weibull plotting position ($i/(n+1)$), according to Rosbjerg et al. [43]. Furthermore, as reported by Veysseyre [44], it provides an estimation error below 1%, whatever the sample size.

$$ \hat{F}(x_i) = \frac{i - 0.3}{n + 0.4} \qquad (1) $$

where $\hat{F}$ is a cumulative distribution function that assigns a probability to each of the $n$ results of an experiment, denoted as $x_1, \dots, x_i, \dots, x_n$, and ordered in ascending order.

The results for the stress levels of 1900 MPa and 2000 MPa are presented in Fig. 13. A gap between the data points can be observed around $N \approx 2 \times 10^5$ cycles, accompanied by a change in the slope of the curve. This change in slope likely indicates that the points follow a different probability distribution. The region where this transition occurs coincides with the fatigue life where a change in the crack initiation mechanisms is observed.
<figure><img src="image_15.png" /><figcaption>Fig. 13. Cumulative failure probability $\hat{F}(N)$ as a function of the number of cycles to failure. (a) for notched specimens tested at 1900 MPa and (b) for gear specimens tested at 2000 MPa.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Weibull distribution $\Psi(\hat{F})$ as a function of N, calculated with the fatigue results and associated linear regression.</figcaption></figure>

<table><caption>**Table 5**

Linear form $Y = AX + B$ of probability distributions used in the construction of probability plots.</caption><thead><tr><th>Distribution</th><th>Linear form</th><th>X</th><th>Y</th><th>A</th><th>B</th></tr></thead><tbody><tr><td>Log-normal</td><td>$\Phi^{-1}(F(x)) = \frac{1}{\sigma}\ln(x) - \frac{\mu}{\sigma}$</td><td>$\ln(x)$</td><td>$\Phi^{-1}$</td><td>$\frac{1}{\sigma}$</td><td>$-\frac{\mu}{\sigma}$</td></tr><tr><td>Weibull</td><td>$\Psi(F(x)) = k \ln(x) - k \ln(\lambda)$</td><td>$\ln(x)$</td><td>$\Psi$</td><td>$k$</td><td>$-k \ln(\lambda)$</td></tr></tbody></table>

Note : $\Psi(F) = \ln(-\ln(1 - F))$.

## 4.2. Analytical cumulative distribution function

In order to determine which probability distribution best describes the empirical CDF data (Log-normal or Weibull), and to determine their respective parameters, the linearized forms of these distributions will be used, as presented in Table 5.

The empirical CDF is used to calculate the analytical functions $\Psi(\hat{F})$ and $\Phi^{-1}(\hat{F})$ for the log-normal and Weibull distributions, respectively. The results for the Weibull distribution are presented in Fig. 14.

The change in slope observed in the empirical CDF remains evident for both applied distributions. A linear regression was performed on the data to assess the fit accuracy of the models across all points, resulting

in a determination coefficient of $R^2 \approx 0.92$ for the point at 1900 MPa and $R^2 \approx 0.72$ for 2000 MPa.

Based on the conclusions of the experimental study and the variations in slope observed, the test results were divided into two sub-populations. The division was made based on the fatigue life and the two populations identified through the fractographic analysis (i.e. short and long fatigue lives). Each sub-population was then analyzed separately by recalculating the functions $\Psi(\hat{F})$ and $\Phi^{-1}(\hat{F})$ individually.

Fig. 15 confirms that, once the data are separated into sub-populations, each group exhibits much better linearity, supporting the hypothesis of bi-modal fatigue life, related to distinct fatigue behaviors associated with different crack initiation mechanisms.

The linear regression lines, applied to each subgroup show a significantly improved fit to the experimental data. The regression coefficients are systematically greater than $R^2 > 0.97$, indicating a substantial improvement compared to the analysis without data separation.

This analysis confirms the bi-modal nature of the fatigue life distributions and the link between the competing crack initiation mechanisms and the scatter in the fatigue life. The fatigue design of a structure using only a uni-modal distribution function overestimates the scatter and leads to an overly conservative design. It is therefore necessary to
<figure><img src="image_17.png" /><figcaption>Fig. 15. Weibull distribution $\Psi(\hat{F})$ as a function of N, calculated with the fatigue results and associated linear regression. The data are separated into sub-populations.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Cumulative distribution functions of the lifetime of carburized M50NiL, identified from a mixture model of Weibull laws (denoted as W) and a lognormal laws (denoted as Log-N).</figcaption></figure>

divide the data into two groups based on fatigue life and to model them separately in order to more accurately represent the fatigue behavior of the material.

With these results, all the necessary elements for the identification of a mixture law are available.

### 4.3. Bi-modal mixture model

Mixture models are statistical tools used to represent data from multiple sub-populations within a larger population, without requiring explicit identification of the sub-population to which each observation belongs. These models assume that the observed data are generated by a combination of different probability distributions, with each component representing a distinct sub-population [45]. The general formula of a mixture model is a weighted sum of the probability density functions of each sub-population, and can be written as a linear combination, as shown in Eq. (2).

$$F(x; p, \Theta) = \sum_{i=1}^{n} p_i F_i(x; \Theta_i), \quad p_i \ge 0, \quad \sum_{i=1}^{n} p_i = 1 \qquad (2)$$

where $p = [p_1, p_2, ..., p_n]$ and $\Theta = [\Theta_1, \Theta_2, ..., \Theta_n]$, with $p_i$ representing the weight coefficient of the distribution function $F_i$ with parameters $\Theta_i$.

The weight coefficients are calculated based on the proportion of data points in the two sub-populations. In general, the classification can be guided by fractographic observations, although explicit information concerning the different crack initiation types is not required to construct the model.

This type of modeling is performed at a single stress level at a time. However, applying this modeling across multiple stress levels enables the reconstruction of Wöhler curves associated with different failure probabilities. The data generated by the model also allows the construction of cumulative distribution functions at various stress levels [12].

#### 4.3.1. Identification of a mixture model

The bi-modal mixture model is defined by the parameter $p$, which represents the weight of one of the two distributions (and $1 - p$ the weight of the other distribution), as well as by the type of distribution employed. Both the Weibull and Log-normal distributions are investigated to construct the model, which is identified based on the fatigue
<figure><img src="image_19.png" /><figcaption>Fig. 17. Cumulative distribution functions of fatigue life for all stress levels. Gray points indicate the loading levels where only notched specimens were tested, while the three blue colors are the levels where only gears were analyzed.</figcaption></figure>

test results. Thus, the cumulative distribution function of a bi-modal mixture model can be expressed by Eqs. (3) and (4), for the Weibull and log-normal distributions, respectively.

$$F_W(N, \sigma_{\max}) = pF_{W1}(k_1, \lambda_1) + (1-p)F_{W2}(k_2, \lambda_2) \quad (3)$$

$$F_{LN}(N, \sigma_{\max}) = pF_{LN1}(\mu_1, \sigma_1^2) + (1-p)F_{LN2}(\mu_2, \sigma_2^2) \quad (4)$$

Each bi-modal cumulative distribution function requires the identification of five parameters: p is determined from the fatigue test results, while the parameters of the individual distributions are identified through linear regressions applied to the functions $\Psi(\hat{F})$ and $\Phi^{-1}(\hat{F})$.

The identified bi-modal cumulative distribution functions for the 1900 MPa and 2000 MPa stress levels are shown in Fig. 16.

The use of either the Log-normal or the Weibull distributions leads to very similar results, which is consistent with the fact that the $R^2$ coefficients were comparable for both cases. Fig. 16 shows that good agreement between the experimental results and the bi-modal mixture model is observed, demonstrating that a bi-modal distribution provides the most suitable description and confirming the existence of bi-modal fatigue lives for the carburized M50NiL steel.

By applying this modeling approach to all stress levels, without mixing the data from gears and notched specimens, the results shown in Fig. 17 are obtained.

The behavior observed at different stress levels is similar to that reported by Sakai [12] for the SUJ2 bearing steel tested in rotating bending. In that study, a mixed model based on a Weibull distribution was shown to effectively describe the fatigue behavior at stress levels exhibiting bi-modal life distributions.

It is important to recall that the weights used in the mixture model depend on the analyzed stress level, as shown in Fig. 18. These weights are influenced by the number of experimental points available at each stress level and by how these points are distributed between the sub-populations. As a result, the relative weight of each sub-population may vary from one stress level to another. For certain levels, the first sub-population dominates, while for others, the second is more significant. At very high or very low stress levels, the bi-modal model tends to converge toward an uni-modal behavior. Thus, it is not possible to use this type of model to simultaneously describe all stress levels.

One limitation of this model is that the identification of the parameters governing each sub-population is performed through linear

<figure><img src="image_20.png" /><figcaption>Fig. 18. Dependence of the weight of each sub-population in the bi-modal distribution.</figcaption></figure>

regression. As a result, at least two data points are required in each sub-population to model a given stress level. Therefore, bi-modal stress levels such as 2100 MPa and 1800 MPa, for which only one point is available in one of the sub-populations, cannot be modeled using this approach. This modeling approach is also limited by its dependence on the weighting of the distributions, which varies with the stress level.

In addition, the estimation of the mixture weights carries non-negligible uncertainty. At very high or low stress levels, the analytical results may exhibit significant dispersion due to the limited number of tested specimens. Consequently, additional experimental data are required to improve the robustness of the estimated occurrence probabilities

Finally, it only considers the statistical distribution of the data, without accounting for the influence of mechanical properties on the formation of the different populations. Therefore, it is necessary to introduce a more robust model that incorporates the influence of the material property distributions on the fatigue behavior at different stress levels.
# 5. Fatigue behavior model

## 5.1. Probabilistic model for bi-modal behavior in fatigue

To model Wöhler curves for a material exhibiting bi-modal fatigue behavior at one or more stress levels, it is necessary to explicitly represent the physical mechanisms involved in crack initiation. These mechanisms must be considered in terms of their distinct natures and activation conditions. The present model incorporates assumptions supported by both experimental results and prior studies in the literature.

The following hypotheses are considered in the formulation of the model:

1. The fatigue life distributions may be bi-modal at certain stress levels and uni-modal at others;

2. The observed bi-modality is due either to a competition between different crack initiation mechanisms or to the same mechanism(s) occurring over very different numbers of cycles;

3. The activation of a given mechanism depends not only on its dominance but also on the presence of favorable local conditions, such as inclusions or surface topography features;

4. In the short-life regime, the fatigue life is often dominated by the crack propagation phase;

5. In the long-life regime, the propagation time is negligible when compared to the total fatigue life, and failure is governed primarily by crack initiation;

6. The relative proportion of short-life and long-life events depends on the applied maximum stress level.

It is important to note that the distinction between short and long fatigue lives is material-dependent and not always clearly observable.

If these conditions are present at least at one stress level, the fatigue behavior can be described as a combination of two failure probabilities: $P_A = P(A)$, corresponding to short-life events, and $P_B = P(B)$, corresponding to long-life events. Thus, the probability of failure before $N$ cycles is given by $P_f = P(N_f \le N)$, which can be expressed as a function of $P_A$ and $P_B$.

### 5.1.1. Short-life population ($P_A$)

For $P_A$, the fatigue life is essentially governed by the number of cycles required for a crack to propagate from its initial size to a critical length at which failure occurs according to the chosen failure criterion. For this propagation stage, it is assumed that the hypotheses of Linear Elastic Fracture Mechanics (LEFM) are valid.

The modeling of event A considers two successive subprocesses:

* Crack initiation ($A_i$): formation of an initial crack of size $c_i$ (demi-width) during the early loading cycles;

* Crack propagation ($A_p$): growth of the crack from $c_i$ to a critical size $c_f$ over a number of cycles $N$.

The occurrence of event A requires that both subprocesses take place, such that:

$$P_A = P(A_p \cap A_i) = P(A_i)P(A_p | A_i) \quad (5)$$

Here, $P(A_p | A_i)$ represents the conditional probability of propagation given that a crack has been initiated. These events are dependent, as propagation cannot occur without a pre-existing crack.

The initiation probability $P(A_i)$ depends on the maximum applied hot-spot stress $\sigma_{max}$, and is denoted by $P_A(\sigma_{max}) = P(A_i)$. In turn, the conditional propagation probability $P(A_p | A_i)$, denoted by $P_{A_p|A_i}(N, \Delta\sigma)$, depends on the stress range $\Delta\sigma$, that can be expressed as a function of the maximum applied hot-spot stress and the stress ratio: $\Delta\sigma = \sigma_{max}(1 - R)$. Inspired by LEFM, propagation is only possible if $\Delta\sigma$ is sufficiently high for the stress intensity factor to exceed the threshold for crack growth.

### 5.1.2. Long-life population ($P_B$)

Long fatigue lives are associated with situations in which the initiation of a crack capable of propagating and satisfying the failure criterion requires a large number of loading cycles. In such cases, the number of cycles to crack initiation becomes the dominant factor in the total fatigue life of the specimen. The event $B_i$ is thus defined as the crack initiation occurring after a significant number of cycles $N$. The corresponding probability is denoted by $P_{B_i}(\sigma_{max}) = P(B_i)$, reflecting its dependence on the applied maximum stress $\sigma_{max}$, then:

$$P_B = P_{B_i} \quad (6)$$

### 5.1.3. Probability of failure $P_f$

The total probability of failure $P_f$ can be expressed by combining the competing events A and B using a competing risks model (also known as the weakest-link model), resulting in:

$$P_f = 1 - (1 - P_A)(1 - P_B) = P_A + P_B - P_A P_B \quad (7)$$

Substituting the expressions for $P_A$ and $P_B$, we obtain:

$$P_f = P_{B_i} + P_{A_i} P_{A_i|A_i} [1 - P_{B_i}] \quad (8)$$

And in expanded form, as a function of the input parameters:

$$P_f(N, \Delta\sigma, \sigma_{max}) = P_{B_i}(N, \Delta\sigma) + P_{A_i}(\sigma_{max}) P_{A_i|A_i}(N, \Delta\sigma) [1 - P_{B_i}(N, \Delta\sigma)] \quad (9)$$

This expression is conceptually similar to the one proposed by Paolino et al. [46], although it is based on different initial assumptions.

This modeling approach provides a global representation of the fatigue behavior based on the maximum applied hot-spot stress at the specimen surface. In the following sections, for the development of this model, the stress ratio will be considered constant and equal to 0.05.

## 5.2. Probabilistic S-N curves

### 5.2.1. Identification of probability $P_{A_i}$

For the modeling of $P_{A_i}$, it is assumed that the number of active crack initiation sites in a specimen is proportional to the applied maximum principal stress. In some cases, these initiation sites may be activated as early as the first loading cycle [47,48].

Assuming that the number of initiation sites depends on the maximum applied stress $\sigma_{max}$, and that their spatial density, denoted by $\lambda_{A_i}$, follows a power-law distribution of the form:

$$\lambda_{A_i}(\sigma_{max}) = \lambda_0 \left( \frac{\sigma_{max}}{\sigma_0} \right)^{k_{A_i}} \quad (10)$$

where $k_{A_i}$ is the shape parameter and $\lambda_0$ is the scale parameter. It is further assumed that the probability of having $n_{is}$ initiation sites within a domain $\Omega$ of size $V_\Omega$ follows a Poisson distribution:

$$P_{n_{is}}(\Omega, \sigma_{max}) = \frac{[n(\Omega)]^{n_{is}}}{n_{is}!} \exp[-n(\Omega)] \quad (11)$$

where $n(\Omega) = V_\Omega \lambda_{A_i}$ is the expected number of initiation sites in the domain $\Omega$.

The probability that at least one initiation site occurs within the volume, leading to crack initiation, is then given by:

$$P_{A_i}(\sigma_{max}) = \mathbb{P}(n_{is} \ge 1) = 1 - \mathbb{P}(n_{is} = 0) = 1 - \exp[-n(\Omega)] \quad (12)$$

Substituting $n(\Omega)$ and adopting $\lambda_0 = 1/V_0$, the expression becomes:

$$P_{A_i}(\sigma_{max}) = 1 - \exp \left[ - \frac{V_\Omega}{V_0} \left( \frac{\sigma_{max}}{\sigma_0} \right)^{k_{A_i}} \right] \quad (13)$$

This expression corresponds to a two-parameter Weibull distribution, where $k_{A_i}$ is the shape parameter, and the scale parameter is given by $\sigma_{A_i} = V_\Omega/V_0\sigma_0^{k_{A_i}}$.
<figure><img src="image_21.png" /><figcaption>Fig. 19. $P_{Ai}$ distribution function for short-life failures as a function of $\sigma_{max}$, based on the observed proportion of short-life points at each stress level.</figcaption></figure>

This approach for modeling the formation of crack initiation sites has been adopted in several studies [49–51]. However, the available data do not allow for the direct estimation of the spatial density $\lambda_{Ai}$, or the scale parameter $\sigma_{Ai}$ as defined in Eqs. (10) and (13).

Nevertheless, considering that crack initiation during the first load cycles necessarily leads to a short fatigue life, it is proposed to approximate $P_{Ai}$ by the observed proportion of short-life points at a given stress level.

Based on this assumption, a two-parameter Weibull distribution can be fitted to the experimental data, enabling the identification of the corresponding shape and scale parameters.

The gear specimens exhibit higher residual compressive stresses at the surface. These stresses can influence the crack propagation kinetics and, consequently, affect the modeling of fatigue life in the short-life regime. For this reason, the two types of specimens investigated in this study were modeled separately. The results obtained for the short-life fatigue regime in the notched and the gear specimens are shown in Figs. 19(a) and 19(b), respectively.

### 5.2.2. Identification of probability $P_{A_p|A_i}$

The conditional probability $P_{A_p|A_i}$ represents the likelihood that a crack with initial size $c_i$ propagates to a critical size $c_f$ after N cycles, given that the crack has already been initiated. To describe this probability, a probabilistic crack growth model is required, where the random variable may be the initial crack size $c_i$, one or more parameters of the propagation law, or both.

In this context, Ciavarella and Papangelo [52] proposed models to address such cases. They showed that if the initial crack size $c_i$ is the only random variable and follows a Fréchet distribution:

$$F(c_i) = \exp \left[ - \left( \frac{\beta_{c_i}}{c_i} \right)^{\alpha_{c_i}} \right] \quad (14)$$

then the propagation life follows a two-parameter Weibull distribution. By integrating the Paris law and neglecting the effect of the crack size $c_f$, as in the approach proposed in [52], the failure probability due to propagation is given by:

$$P_{f_p}(N, \Delta\sigma) = W \left( \frac{2\alpha_{c_i}}{m-2}, \frac{\beta_{c_i}^{1-m/2}}{Y^m \pi^{m/2}(m/2-1)C \Delta\sigma^m} \right) \quad (15)$$

Where $\alpha_{c_i}$ and $\beta_{c_i}$ are the parameters of the initial crack size distribution, m and C are the parameters of the Paris law. Y is a geometric factor used in the calculation of the stress intensity factor ($\Delta K = Y \Delta \sigma \sqrt{\pi a}$).

Thus, the expression for $P_{f_p}$ can be simplified to:

$$P_{f_p}(N, \Delta\sigma) = W \left( \frac{2\alpha_{c_i}}{m-2}, \frac{\xi}{(\Delta\sigma)^m} \right) \quad (16)$$

<figure><img src="image_22.png" /></figure>

Hence, the conditional probability $P_{A_p|A_i}$ can be expressed as a Weibull distribution with two parameters:

$$P_{A_p|A_i}(N, \Delta\sigma) = 1 - \exp \left[ - \left( \frac{N}{N_{A_p}(\Delta\sigma)} \right)^{k_{A_p}} \right] \quad (17)$$

where:

$$k_{A_p} = \frac{2\alpha_{c_i}}{m-2}, \quad N_{A_p}(\Delta\sigma) = \frac{\xi_A}{(\Delta\sigma)^{m_A}} \quad (18)$$

Here, $N_{A_p}$ represents the propagation life and has a form analogous to the Basquin law.

If the variability in the initial crack size $c_i$ is not sufficient to describe the scatter in the propagation life, the Paris law parameter C may also be treated as a random variable [52].

The resolution of the propagation probability described in Eqs. (17) and (18) can be approached in two distinct ways.

The first is an analytical method, which involves the experimental identification of the parameters ($k_{A_p}$, $m_A$ and $\xi_A$). These parameters characterize the crack propagation kinetics. To determine them, dedicated crack growth tests must be carried out, and a Paris law can be identified that adequately represents crack growth in the material. However, this approach can become complex due to the influence of the material property gradients generated by the thermo-chemical treatment, as well as stress distributions generated by the specimen geometry and bending loading condition.

A second and more empirical approach consists of using statistical methods such as Maximum Likelihood Estimation (MLE) to fit the distribution directly to the experimental fatigue data ($\Delta\sigma, N$) obtained in the short-life regime.

In the following, the model parameters ($m_A$, $\xi_A$, $k_{A_p}$) were estimated by maximizing the likelihood function, which was minimized numerically using the Nelder-Mead algorithm. The identified values are presented in Table 6, and the corresponding life estimations are illustrated in Fig. 20.

This statistical approach shows good agreement with the experimental data and offers a practical and robust alternative for fatigue life estimation, especially in situations where detailed crack growth data is not available.

### 5.2.3. Identification of probability $P_{B_i}$

In this case, the number of cycles to crack initiation becomes the dominant factor in the total fatigue life. Due to the limited number of results in the long-life regime ($N > 2 \times 10^5$) and the fact that these data reflect a combination of two distinct mechanisms, the probability
<figure><img src="image_23.png" /><figcaption>Fig. 20. Simulated fatigue lives for case hardened M50NiL steel based on the identified $P_{A|A_i}$ and $P_{B_i}$ functions.</figcaption></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 21. Failure probability $P_f$ for case hardened M50NiL steel as a function of the number of cycles and the applied stress, according to the proposed model. The model is valid only within the experimental data domain.</figcaption></figure>

<figure><img src="image_26.png" /></figure>

associated with event $P_{B_i}$ was modeled using a log-normal distribution, expressed as:

$$P_{B_i}(N, \Delta\sigma) = \frac{1}{2} \left[ 1 + \operatorname{erf} \left( \frac{N - \mu_B}{s_B \sqrt{2}} \right) \right], \qquad (19)$$

where $s_B$ is the standard deviation, $\mu_B$ is the mean, and erf is the error function. To account for the variation of the mean with the stress amplitude, $\mu_B$ is assumed to follow a Basquin-type relation:

$$\mu_B = \frac{\xi_B}{(\Delta\sigma)^{m_B}}, \qquad (20)$$

where $\xi_B$ and $m_B$ are two parameters to be identified.

The identification of these parameters was also performed using the Maximum Likelihood Estimation (MLE) method, which was applied to fit the probabilistic fatigue model defined by Eqs. (19) and (20) to the experimental data. The estimated parameters are presented in Table 6, and the corresponding fatigue life predictions are shown in Fig. 20.

<figure><img src="image_27.png" /><figcaption>Fig. 22. Stress profiles as a function of depth for a maximum applied hot-spot stress of 1900 MPa: notch stress distribution, residual stress distribution (notched and gear specimens), local maximum stress, and corresponding local stress ratio.</figcaption></figure>

### 5.2.4. Identification of probability $P_f$

The total probability of failure $P_f(N, \sigma_{max})$ is formulated as a weighted combination of two failure modes: mode A, which dominates in short-life regimes, and mode B, associated with long-life failures. The estimated values for the parameters of the different probability functions are presented in Table 6.

Based on these results, the total failure probability $P_f$ can be calculated for different failure probabilities. Accordingly, Wöhler curves representing the fatigue behavior of the carburized M50NiL steel were plotted for failure probabilities of 10%, 50%, and 90%. The resulting S-N diagram is shown in Fig. 21.

The proposed model with the identified parameters results in good reproduction of the experimental fatigue behavior. Most of the experimental results fall within the range $P_f \in [0.1, 0.9]$, indicating that the model is capable of capturing the variability of fatigue life across different stress levels while accounting for bi-modal behavior.
<figure><img src="image_28.png" /><figcaption>**Fig. 23.** STBF and notched specimen results in the form of a Wöhler diagram based on the local maximum stress and local stress ratio. The red vertical line indicates the division between the two populations (short and long lives).</figcaption></figure>

<table><caption>**Table 6**
The model parameter values estimated by the maximum likelihood method.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="2">P<sub>Ai</sub></th><th colspan="3">P<sub>Ap|Ai</sub></th><th colspan="3">P<sub>Bi</sub></th></tr><tr><th>k<sub>Ai</sub></th><th>σ<sub>Ai</sub></th><th>k<sub>Ap</sub></th><th>ξ<sub>A</sub></th><th>m<sub>A</sub></th><th>s<sub>B</sub></th><th>ξ<sub>B</sub></th><th>m<sub>B</sub></th></tr></thead><tbody><tr><td>Notched</td><td>10,27</td><td>2032</td><td>0,99</td><td>1,77 × 10<sup>17</sup></td><td>3,97</td><td>1,34</td><td>1465</td><td>0,60</td></tr><tr><td>Gear</td><td>26,72</td><td>1965</td><td>2,49</td><td>3,57 × 10<sup>25</sup></td><td>6,35</td><td>0,64</td><td>134</td><td>0,29</td></tr></tbody></table>

The model's response varies according to the values of the different estimated parameters. The probability $P_{B_i}$ follows a three-parameter lognormal distribution. The standard deviation $s_B$ governs the scatter in fatigue life, particularly at long lifetimes, i.e., the distance between the 10% and 90% Wöhler curves. The mean of the distribution follows a Basquin-type relation, dependent on the applied stress amplitude and controlled by two parameters: $\xi_B$, shifts the Wöhler curves (higher values correspond to longer predicted life); $m_B$, defines the slope of the Basquin curve (higher values lead to steeper curves).

The conditional probability $P_{A_i|A_i}$ also follows a Weibull distribution with two parameters, analogous in behavior to those of $P_{B_i}$. The scale parameter follows a Basquin-type law: $k_{A_i}$ controls the statistical scatter of the propagation life, $\xi_A$ shifts the life prediction curves, $m_A$ governs the life reduction between high and low stress levels.

The probability $P_{A_i}$ is described by a two-parameter Weibull distribution defined by $k_{A_i}$ and $\sigma_{A_i}$, which captures the bi-modal nature of the Wöhler curves: $\sigma_{A_i}$ influences the stress level at which the transition between short-life and long-life regimes occurs, $k_{A_i}$ defines the slope of the distribution. Higher values result in a steeper transition and a narrower stress range exhibiting bi-modality.

The proposed model makes it possible to construct a complete Wöhler diagram for materials exhibiting competing crack initiation mechanisms, which result in bi-modal fatigue lives. It provides robust predictions within the stress range covered by the experimental data, where the parameters have been identified (see the validity domain in Fig. 21). However, outside this domain — particularly in the very high cycle fatigue (VHCF) regime — extrapolation of the model is not recommended. In this regime, it appears necessary to introduce a fatigue threshold, highlighting the need for additional experimental data to validate and extend the predictive capacity of the model. Furthermore, a local analysis of how different property gradients influence the local fatigue behavior should be carried out.

## 6. Local analysis of fatigue data

The proposed probabilistic model investigates the bi-modal fatigue behavior from a global perspective, considering the maximum applied hot-spot stress at the notch root surface ($x = 0$) under a fixed reference

<figure><img src="image_29.png" /><figcaption>**Fig. 24.** S-N diagram based on local maximum stress, after correction to a common stress ratio. The red vertical line indicates the division between the two data populations (short and long lives).</figcaption></figure>

stress ratio ($R_{ref}$). However, due to the applied stress distribution and the compressive residual stress field, both the maximum stress and the stress ratio vary locally with the crack initiation depth ($x$).

To capture this variability, the local maximum stress $\sigma_{max}^*(x)$ and the corresponding local stress ratio $R^*(x)$ were determined through a two-step procedure.

**Step 1.** Based on the crack initiation depth ($x$) and $R_{ref}$, the maximum stress $\sigma_{max}(x)$ was obtained from the notch stress distribution along a path normal to the hot-spot surface. The minimum stress followed as $\sigma_{min}(x) = R_{ref}\sigma_{max}(x)$. The corresponding mean and alternating stresses were then derived:

$$ \sigma_m(x) = \frac{\sigma_{max}(x) + \sigma_{min}(x)}{2}, \quad \sigma_a(x) = \frac{\sigma_{max}(x) - \sigma_{min}(x)}{2} \quad (21) $$

**Step 2.** The effect of the residual stress field was then included. At the crack initiation depth, the local residual stress $\sigma_{res}(x)$ was added to the applied mean stress, giving:

$$ \sigma_m^*(x) = \sigma_m(x) + \sigma_{res}(x). \quad (22) $$

Using this corrected mean stress, the local maximum and minimum stresses were reconstructed as:

$$ \sigma_{max}^*(x) = \sigma_m^*(x) + \sigma_a^*(x), \quad \sigma_{min}^*(x) = \sigma_m^*(x) - \sigma_a^*(x) \quad (23) $$

Where $\sigma_a^*(x) = \sigma_a(x)$. Finally, the local stress ratio was obtained as:

$$ R^*(x) = \frac{\sigma_{min}^*(x)}{\sigma_{max}^*(x)} \quad (24) $$

**Fig. 22** illustrates the evolution of $\sigma_{max}^*(x)$ and $R^*(x)$ with depth from the notch-tip hot-spot for both specimen types (i.e. notched and gear specimens), assuming an applied maximum hot-spot stress of 1900 MPa at the surface. The figure also shows the residual stress profile $\sigma_{res}(x)$ and the notch stress distribution $\sigma_{max}(x)$ used in the calculations.

Near the surface, $R^*(x \approx 0)$ varies more sharply due to the gradient of residual stress, especially in the gear specimens, which exhibit a greater compressive residual surfaces stresses. At greater depths, the residual stress profiles of both specimen types converge to an almost constant, leading to a similar evolution of $R^*(x)$.

For each tested specimen, the $\sigma_{max}^*(x)$ and $R^*(x)$ were calculated at the crack initiation position, determined experimentally after failure. The resulting S-N curve, constructed using the local maximum stress, is shown in Fig. 23. Note that each data point has a different local stress ratio which is presented by the color bar.

In terms of local maximum stress, in the short-life regime (i.e. the first population), the lifetime variability between the two specimen
types appears reduced compared to the initial observation. This suggests that, at the local scale, their fatigue behavior is comparable. It supports the hypothesis that the scatter initially observed in the S-N curve (Fig. 7) is largely driven by the compressive residual stress at the gear surface. The $\sigma_{\max}^*(x)$ also gives the impression that the fatigue strength of the first population is lower than that of the second. However, the local stress ratio $R^*$ varies more significantly in the first population, particularly in the gear specimens, where it ranges from -0.4 to -0.7, while in the notched specimens, the variation is approximately between -0.18 and -0.3.

For the second population (i.e. the long-life regime with internal crack initiation), the increase in crack initiation depth at inclusions, driven by surface residual stress peaks, leads to a decrease in the applied local stress and to a convergence of residual stress values for both specimen types. As a result, the variation in local stress ratio is lower than for the first population and is similar for both specimen types, ranging from -0.1 and -0.25. Only one data-point in this population has a $R^*$ close to -0.7. This point corresponds to a surface failure and exhibits behavior more consistent with the first population.

Overall, the separation of the data into two populations remains clear, highlighting the bi-modal fatigue behavior. At the local scale, each point is associated with a different stress ratio, which must be considered in the analysis. To unify the data under a common stress ratio and enable a more detailed comparison, an equivalent stress approach is proposed, incorporating the effect of the stress ratio.

Building on Murakami's formulation [53], the influence of the stress ratio variation on the stress amplitude is introduced through the factor $(\frac{1-R}{2})^\alpha$. Accordingly, the equivalent stress amplitude at $R = -1$ is given by:

$$\sigma_{a(R=-1)} = \frac{\sigma_a^*(x)}{\left(\frac{1-R^*(x)}{2}\right)^\alpha} \quad (25)$$

Where $\alpha = 0.226 + HV \times 10^{-4}$.

In the present case, to enable a unified comparison, the correction is expressed as a function of the local maximum stress. Given that $\sigma_{a(R=-1)} = \sigma_{\max(R=-1)}$ and $\sigma_a^*(x) = (\frac{1-R^*(x)}{2}) \sigma_{\max}^*(x)$, Eq. (25) can be rewritten in terms of $\sigma_{\max}^*(x)$ as follows:

$$\sigma_{\max(R=-1)} = \frac{\sigma_{\max}^*(x)}{\left(\frac{1-R^*(x)}{2}\right)^{\alpha-1}} \quad (26)$$

The value of $\alpha$ was determined by fixing the hardness at 730 HV. Given that the hardness variation across the observed crack initiation depths was less than 20 HV, its effect on $\alpha$ was considered negligible, as it would only change from 0.299 to 0.297.

Fig. 24 shows the S-N curve as a function of the equivalent maximum stress at $\sigma_{\max(R=-1)}$, where all of the data-points have been corrected to the same stress ratio.

After this correction, the dispersion in terms of local stress is significantly reduced for both populations, compared to the uncorrected case (Fig. 23). In the first population, the gear specimens still exhibit longer fatigue lives than the notched specimens, even after accounting for residual stress and local stress ratio variations. Additionally, the fatigue strength of the two populations appears much closer than in the previous case.

For the second population, the reduced scatter in terms of stress suggests the possible existence of a fatigue threshold in this regime. This stands in contrast to the impression given by the probabilistic model extrapolation, which implies a continuous decrease in life without a defined threshold. Nevertheless, the bi-modal behavior remains clearly visible.

Note that the Murakami expression for calculating $\alpha$ assumes that the hardness is homogeneous. This may underestimate its value in the present case where there are material property gradients. The experimental determination of $\alpha$ would be useful to better analyze the

effect of $R^*(x)$. This would allow a more realistic consideration of the hardness gradient.

The analysis presented above, in terms of the local stress, considers only the applied stress and residual stress distributions. However, a physically representative fatigue life model should also include a correction based on defect size, to account for the observed evolution of defect size (as shown in Fig. 12), as well as hardness and microstructural gradients. The integration of these factors will provide a more comprehensive and realistic representation of the complexity of fatigue behavior in case-hardened M50NiL steel. This extended analysis is currently in progress and will be the subject of a forthcoming publication.

## 7. Conclusions

Based on the present work, the following conclusions can be drawn regarding the fatigue behavior of gear teeth manufactured from the carburized M50NiL steel:

* The gears and notched specimens exhibit very similar fatigue behavior and crack initiation mechanisms, despite a difference in the surface residual stress of approximately -450 MPa, resulting from that fact that the gear specimens were shotblasted and tribofinished and the notched specimens were not. The peak in the compressive residual stress generated by these treatments had a beneficial effect on the fatigue life of the gears in the short-life regime. In the long-life regime, no improvement in fatigue life was observed. The effect seems to be limited to an increase in the crack initiation depth.

* The bi-modal nature of fatigue life was experimentally confirmed through fractographic analyses, indicating the competition between different crack initiation mechanisms, as well as through statistical analyses. A mixture model was identified for both datasets: gear specimens and notched specimens.

* The mixture model was able to accurately reproduce the cumulative fatigue life distributions. In this approach, modeling is performed separately at each stress level, since the weight of each population varies depending on the level analyzed. This formulation clearly reveals the presence of bi-modal behavior across different stress levels.

* A probabilistic model was proposed to describe bi-modal fatigue behavior. The choice of probability distributions associated with each term of the model is physically justified by the dominant crack initiation mechanisms in each data sub-population.

* The proposed model successfully captures the overall shape of the Wöhler diagram for a material in which multiple crack initiation mechanisms coexist and compete. Extrapolation of the model outside of the zone of validity defined by the experimental data-points, is delicate and not recommended.

* The local analysis showed that the combined effect of stress distributions (due to the notch, the bending load and the residual stress) significantly influences the local maximum stress and the local stress ratio. Applying corrections to account for these variations reduces the scatter in terms of stress. In terms of fatigue life, a slight beneficial effect of shot blasting can still be observed in the short-life regime, and there is a possible indication of a fatigue threshold in the long-life regime.

The proposed modeling framework can be applied to other materials that exhibit bimodal fatigue behavior, provided that the hypotheses considered in the formulation of the model are satisfied. While the quality of the predictions naturally depends on the amount of available experimental data, the input parameters of the initiation and conditional propagation distributions ($P_{A_i}$ and $P_{A_{p|A_i}}$) can be estimated from complementary sources, when the number of fatigue tests are limited. These include, for example, the spatial density of potential
initiation sites or independent crack-growth measurements. In contrast, the long-life distribution $P_B$ remains dependent on the available fatigue data.

The interplay between property gradients — such as hardness and residual stress — and factors like defect size and depth of inclusions strongly influences the observed bi-modal fatigue behavior. Future work will focus on systematically assessing the role of these parameters under different stress ratios and thermochemical treatments (TCTS). The resulting fatigue data will provide the basis for developing a predictive model that explicitly incorporates their combined effects.

### CRediT authorship contribution statement

**Anicleson Raony Alves de Moura:** Writing – review & editing, Writing – original draft, Validation, Project administration, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Franck Morel:** Writing – review & editing, Validation, Supervision, Resources, Project administration, Methodology, Conceptualization. **Etienne Pessard:** Writing – review & editing, Validation, Supervision, Resources, Methodology, Conceptualization. **Daniel Bellett:** Writing – review & editing, Validation, Supervision, Resources, Methodology, Conceptualization. **Louis Augustins:** Writing – review & editing, Validation, Supervision, Resources, Conceptualization. **Damien Herisson:** Writing – review & editing, Validation, Supervision, Resources, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### References

[1] He Haifeng, Zhou Ye, Liu Heli, Zhu Caichao, Xiao Bo, Zhang Taihua. Study on gear bending fatigue considering gradient characteristics: Numerical analysis and experiments. Eng Fract Mech 2023;277:108983.

[2] Brett Nick, Schurer Stefan, Tobie Thomas, Stahl Karsten, Höhn Bernd-Robert. Investigations on tooth root bending strength of case hardened gears in the range of high cycle fatigue. In: Proceedings of the American gear manufacturers association fall technical meeting, Indianapolis, in, USA. 2013, p. 15-7.

[3] Conrado Edoardo, Gorla Carlo, Davoli Piermaria, Boniardi Marco. A comparison of bending fatigue strength of carburized and nitrided gears for industrial applications. Eng Fail Anal 2017;78:41-54.

[4] Czerwinski Frank. Heat treatment: Conventional and novel applications. BoD-Books on Demand; 2012.

[5] Bonaiti Luca, Bayoumi Ahmed Bayoumi Mahmoud, Concli Franco, Rosa Francesco, Gorla Carlo. Gear root bending strength: A comparison between single tooth bending fatigue tests and meshing gears. J Mech Des 2021;143(10):103402.

[6] Farfan S, Rubio-Gonzalez C, Cervantes-Hernandez T, Mesmacque Gérard. High cycle fatigue, low cycle fatigue and failure modes of a carburized steel. Int J Fatigue 2004;26(6):673-8.

[7] Alban Lester E. Number 1 gear failure—tooth bending fatigue. SAE Trans 1984;1002-8.

[8] Woods Jody Lee, Daniewicz SR, Nellums R. Increasing the bending fatigue strength of carburized spur gear teeth by presetting. Int J Fatigue 1999;21(6):549-56.

[9] Argoud Vincent. Caractérisation et modélisation du comportement en fatigue à grand nombre de cycles d'aciers cémentés : application aux systèmes de transmission de puissance pour l'aéronautique. [Ph.D. thesis], Arts et Métiers Institute of Technology; 2021.

[10] McPherson Douglas R, Rao Suren B. Mechanical testing of gears. In: Mechanical testing and evaluation. ASM International; 2000, p. 861-72.

[11] Argoud Vincent, Morel Franck, Pessard Etienne, Bellett Daniel, Thibault Simon, Gourdin Stéphane. Fatigue behaviour of gear teeth made of case hardened steel: from competing mechanisms to lifetime variability. Procedia Struct Integr 2019;19:719-28.

[12] Sakai Tatsuo, Lian Benning, Takeda Mitsuhiro, Shiozawa Kazuaki, Oguma Noriyasu, Ochi Yasuo, Nakajima Masaki, Nakamura Takashi. Statistical duplex S-N characteristics of high carbon chromium bearing steel in rotating bending in very high cycle regime. Int J Fatigue 2010;32(3):497-504.

[13] Goto M, Yamamoto T, Nisitani H, Sakai T, Kawagoishi N. Effect of grinding on the fatigue strength of a bearing steel in the super long-life field. In: ECF 13-13th European conference on fracture. 2000, p. 2000.

[14] Bayraktar Emin, Garcias Israel Marines, Bathias Claude. Failure mechanisms of automotive metallic alloys in very high cycle fatigue range. Int J Fatigue 2006;28(11):1590-602.

[15] Akiniwa Yoshiaki, Miyamoto Nobuyuki, Tsuru Hirotaka, Tanaka Keisuke. Notch effect on fatigue strength reduction of bearing steel in the very high cycle regime. Int J Fatigue 2006;28(11):1555-65.

[16] Furuya Yoshiyuki. Size effects in gigacycle fatigue of high-strength steel under ultrasonic fatigue testing. Procedia Eng 2010;2(1):485-90.

[17] Sakai Tatsuo. Review and prospects for current studies on very high cycle fatigue of metallic materials for machine structural use. J Solid Mech Mater Eng 2009;3(3):425-39.

[18] Li Wei, Yuan Huang, Sun Zhenduo, Zhang Zhenyu. Surface vs. interior failure behaviors in a structural steel under gigacycle fatigue: Failure analysis and life prediction. Int J Fatigue 2014;64:42-53.

[19] Yang ZG, Li SX, Zhang JM, Zhang JF, Li GY, Li ZB, Hui WJ, Weng YQ. The fatigue behaviors of zero-inclusion and commercial 42CrMo steels in the super-long fatigue life regime. Acta Mater 2004;52(18):5235-41.

[20] Li Wei, Sakai Tatsuo, Wakita Masami, Mimura Singo. Influence of microstructure and surface defect on very high cycle fatigue properties of clean spring steel. Int J Fatigue 2014;60:48-56.

[21] Zhang JM, Li SX, Yang ZG, Li GY, Hui WJ, Weng YQ. Influence of inclusion size on fatigue behavior of high strength steels in the gigacycle fatigue regime. Int J Fatigue 2007;29(4):765-71.

[22] Marines I, Bin X, Bathias C. An understanding of very high cycle fatigue of metals. Int J Fatigue 2003;25(9-11):1101-7.

[23] Cashman Gerald T. A mathematical model for competing failure modes in strain cycle fatigue. 2007.

[24] Caton MJ, Jha SK, Rosenberger AH, Larsen JM. Divergence of mechanisms and the effect on the fatigue life variability of rene'88dt. Superalloys 2004;2004:305-12.

[25] Alexandre Franck. Aspects probabilistes et microstructuraux de l'amorçage des fissures de fatigue dans l'alliage Inco 718 [Ph.D. thesis], École Nationale Supérieure des Mines de Paris; 2004.

[26] Jha SK, Chandran KS Ravi. An unusual fatigue phenomenon: duality of the S-N fatigue curve in the β-titanium alloy Ti-10V-2Fe-3Al. Scr Mater 2003;48(8):1207-12.

[27] Jha SK, Larsen JM, Rosenberger AH, Hartman GA. Dual fatigue failure modes in Ti-6Al-2Sn-4Zr-6Mo and consequences on probabilistic life prediction. Scr Mater 2003;48(12):1637-42.

[28] Szczepanski CJ, Jha SK, Larsen JM, Jones JW. Microstructural influences on very-high-cycle fatigue-crack initiation in Ti-6246. Met Mater Trans A 2008;39:2841-51.

[29] Qian Dongsheng, Wu Lingyan, Wang Feng, Deng Song, Yin Fei, Jiang Shaofeng. Tailored carburization gradient microstructure and enhanced wear properties of M50NiL steel via introduced prior cold rolling. Wear 2024;540:205265.

[30] Zhang CS, Yan MF, Sun Z, Wang YX, You Y, Bai B, Chen L, Long Z, Li RW. Optimizing the mechanical properties of M50NiL steel by plasma nitrocarburizing. Appl Surf Sci 2014;315:28-35.

[31] Wang Fangfang, Zhou Chungen, Zheng Lijing, Zhang Hu. Corrosion resistance of carbon ion-implanted M50NiL aerospace bearing steel. Prog Nat Sci: Mater Int 2017;27(5):615-21.

[32] Böhmer HJ. Residual stresses and material behavior of M50NiL and RBD. Astm Special Technic Public 1993;1195. 34-34.

[33] Wang Fangfang, Zhou Chungen, Zheng Lijing, Zhang Hu. Corrosion resistance of carbon ion-implanted M50NiL aerospace bearing steel. Prog Nat Sci: Mater Int 2017;27(5):615-21.

[34] Zhang Yan, Yang Ming, Long Shaolei, Li Bo, Liang Yilong, Ma Shaowei. Effect of initial state and deformation conditions on the hot deformation behavior of M50NiL steel. Materials 2020;13(23):5367.

[35] Afnor. NF ISO 6336-3 — Calculation of load capacity of spur and helical gears — Part 3: Calculation of tooth bending strength. Assoc Française Norm (AFNOR) 2020.

[36] Lisle Timothy J, Little Christian P, Aylott Chris J, Shaw Brian A. Bending fatigue strength of aerospace quality gear steels at ambient and elevated temperatures. Int J Fatigue 2022;164:107125.

[37] Fouchereau Rémy. Modélisation probabiliste des courbes SN [Ph.D. thesis], Université Paris Sud-Paris XI; 2014.

[38] Jha SK, Caton MJ, Larsen JM. Mean vs. life-limiting fatigue behavior of a nickel-based superalloy. Superalloys 2008;565-72.

[39] Jha SK, Larsen JM, Rosenberger AH. The role of competing mechanisms in the fatigue life variability of a nearly fully-lamellar γ-TiAl based alloy. Acta Mater 2005;53(5):1293-304.
[40] Jha SK, Larsen JM, Rosenberger AH. Towards a physics-based description of fatigue variability behavior in probabilistic life-prediction. Eng Fract Mech 2009;76(5):681-94.

[41] Yao Weixing, Guo Shenjie. VHCF test and life distribution of aluminum alloy LC4CS. Int J Fatigue 2008;30(1):172-7.

[42] Harlow D Gary. Statistical characterization of bimodal behavior. Acta Mater 2011;59(12):5048-53.

[43] Rosbjerc D, Corréa J, Rasmussen PF. Justification des formules de probabilité empirique basées sur la médiane de la statistique d'ordre. Rev Des Sci L'Eau 1992;5(4):529-40.

[44] Veyssyre Renée. Aide-mémoire: Statistique et probabilités pour l'ingénieur-2ème édition. Paris: Dunod / L'Usine Nouvelle; 2006.

[45] Tran Thao, Wiskow Cara, Aziz Mohammad Abdus. Skewed and flexible skewed distributions: A modern look at the distribution of BMI. Am J Undergrad Res 2017;14(2).

[46] Paolino Davide Salvatore, Chiandussi Giorgio, Rossetto Massimo. A unified statistical model for S-N fatigue curves: probabilistic definition. Fatigue Fract Eng Mater Struct 2013;36(3):187-201.

[47] Krauss George. Steels: processing, structure, and performance. ASM international; 2015.

[48] Hyde Robert S, Cohen Rachael E, Matlock David K, Krauss George. Bending fatigue crack characterization and fracture toughness of gas carburized SAE 4320 steel. Technical report, SAE Technical Paper; 1992.

[49] Doudard Cédric, Calloch Sylvain. Influence of hardening type on self-heating of metallic materials under cyclic loadings at low amplitude. Eur J Mech A Solids 2009;28(2):233-40.

[50] Ezanno Anthony, Doudard Cédric, Calloch Sylvain, Millot Thierry, Heuzé J-L. Fast characterization of high-cycle fatigue properties of a cast copper-aluminum alloy by self-heating measurements under cyclic loadings. Procedia Eng 2010;2(1):967-76.

[51] Todinov MT. Is Weibull distribution the correct model for predicting probability of failure initiated by non-interacting flaws? Int J Solids Struct 2009;46(3-4):887-901.

[52] Ciavarella Michele, Papangelo Antonio. On the distribution and scatter of fatigue lives obtained by integration of crack growth curves: Does initial crack size distribution matter? Eng Fract Mech 2018;191:111-24.

[53] Murakami Yukitaka. Metal fatigue: effects of small defects and nonmetallic inclusions. Elsevier; 2002.