Anomalous fracture in two-dimensional rhenium disulfide

See allHide authors and affiliations

Science Advances  18 Nov 2020:
Vol. 6, no. 47, eabc2282
DOI: 10.1126/sciadv.abc2282


Low-dimensional materials usually exhibit mechanical properties from those of their bulk counterparts. Here, we show in two-dimensional (2D) rhenium disulfide (ReS2) that the fracture processes are dominated by a variety of previously unidentified phenomena, which are not present in bulk materials. Through direct transmission electron microscopy observations at the atomic scale, the structures close to the brittle crack tip zones are clearly revealed. Notably, the lattice reconstructions initiated at the postcrack edges can impose additional strain on the crack tips, modifying the fracture toughness of this material. Moreover, the monatomic thickness allows the restacking of postcrack edges in the shear strain–dominated cracks, which is potentially useful for the rational design of 2D stacking contacts in atomic width. Our studies provide critical insights into the atomistic processes of fracture and unveil the origin of the brittleness in the 2D materials.


Although two-dimensional (2D) materials are known to have extraordinary intrinsic mechanical properties (1, 2), such as ultraflexibility and ultrahigh strength in graphene and 2D transition metal dichalcogenides (TMDs) (3, 4), in general they are brittle and prone to fracture. Their high strengths are based on the fact that most 2D materials have much less defects than bulks, with some being entirely defectless (5). In the presence of defects or initiated cracks, their strengths can be markedly lowered, and the 2D materials usually exhibit sudden and unstable failure (610). Therefore, to realize the claimed applications and intriguing usage of 2D materials in the future, it is important to understand their failure mechanism thoroughly, in particular the development of cracks in the 2D materials.

Fracture of materials can be classified into two main groups: brittle and plastic (11). Plasticity associated with fracture, including dislocation, twining, amorphorization, etc., is initiated in the stress-focused crack tip zone when the local stress level exceeds the respective thresholds (12). The plastic deformation zone in crack tips can relieve the stress and strain, enlarging the mechanical energy absorption area during fracture, like in 2D MoS2 as we have observed (13). However, as we will present here, in the anisotropic 2D rhenium disulfide (ReS2), the plastic lattice reconstructions initiated from the postcrack edges (not at the maximum strain position) may conversely increase the stress in crack tip, causing abnormal “self-persistent” fracture in this material.

The three fracture modes (I, II, and III) (14), corresponding to opening, shearing, and tearing loads, are unambiguously defined for the fracture in 3D bulk material. In terms of 2D atomic layers, here we show that the shear strain along the crack line (hence mode II) can trigger crystal rotation and atomically parallel restacking of the postcrack edges. This particular phenomenon is only allowed in such ultrathin flexible 2D materials, where the restacking of these atomic layers does not cost a high level of bending energy.

Despite the importance of fracture mechanics in 2D, currently, most of the studies used atomistic simulations [molecular dynamics or density functional theory (DFT)] to understand the fracture processes (1520). Previous microscopic or mesoscopic fracture or tensile experiments on 2D materials, either by atomic force microscopy indentation or microelectromechanical system stages (68), have not directly addressed the atomic-scale processes. Therefore, the real stress/strain in the area of interest and the stress/strain calculated by the boundary conditions subjected to external loading in experiments might be different, limiting the accuracy of results from mechanical tests. In the present work, the postcrack lattice reconstructions notable modify the strain inside or close to the crack tips, so that the experimental fracture toughness measured by mechanical tests cannot be directly correlated with the real atomic-scale processes inside the crack tip, which is entirely brittle.


2D ReS2 in ambient condition has a monoclinic crystal structure, and the 2D [monatomic layer (1L)] form is distorted from the normal tetragonal (T) phase of TMDs (DFT relaxed structure is shown in Fig. 1A; also, see the definition of a, b, and i for basis vectors) (21). The 1L ReS2 in our experiments was synthesized via the chemical vapor deposition (CVD) method (see Materials and Methods) and transferred onto quantifoil transmission electron microscopy (TEM) grids by the polymethyl methacrylate (PMMA) method (see Materials and Methods) (22). During the cooldown process of CVD, a high level of internal strain in the ReS2 flakes caused by thermal expansion mismatch can trigger cracking in the samples (Fig. 1B). The cracks preferentially follow the low-index crystal planes of ReS2 and show atomically smooth postcrack edges on both sides. Unexpectedly, we have frequently found the two opposite fresh crack edges in different crystal orientations (Fig. 1, C and D). Our 1L ReS2 samples are grown as single-crystal flakes on ca. 20-μm scale (fig. S1), and the selected area diffraction pattern (SAED) shows the single-crystalline structures in pristine samples (fig. S2). Moreover, we have used ex situ TEM and frequently observed the two sides of the postcracked surfaces having different orientations, with only short time (<10 s) low-dose electron beam exposure under low-magnification observations (<5M×), suggesting that the crystalline reorientation observed is not associated with electron beam effect. Referring to previous publications (3, 9, 23) and our own observations (fig. S3), cracks do not necessarily follow the defective boundaries in 2D materials. In addition, the directly joint fresh crack edges with orientations shown in Fig. 1 (C and D) can have large lattice mismatches on the interface, which is less likely. Therefore, the different crystal orientations in opposite fresh crack edges should be attributed to the straining process during fracture.

Fig. 1 The cracks in ReS2.

(A) Stick-bond model of 1L ReS2 and the DFT-calculated structure (unit cell) with basis vector definitions. (B) Low-magnification STEM image showing cracks in 1L ReS2 in the mesoscopic scale. Scale bar, 200 nm. (C) High-magnification STEM image for postcrack edges in 1L ReS2. The lattice orientations for different cracked parts have been modified from the original single crystal (false colored to highlight different domains). Scale bar, 2 nm. (D) Another crack image showing the lattice reconstructions in right side crack edge. Scale bar, 2 nm.

By using the in situ scanning TEM (STEM) technique (see Materials and Methods), we have intentionally created cracks in 1L ReS2 samples and tracked the crack propagation inside the samples by in situ observations. The original ReS2 samples are single crystalline, and the lattice reconstructions (24, 25) can be clearly observed around the crack tip zones (within 20-nm ranges) (Fig. 2A). To minimize the energy of newly formed boundaries between the original and reconstructed domains, the boundaries are predominantly twin-like boundaries. However, the twinning here is different from the ordinary mechanically induced twining. It actually involves the lattice reconstruction and swop between the a, b, and i directions. Figure S4 presents the magnified image for the lattice reconstructed zone. The relaxed twin boundary structures show a relatively small shear angle (1° to 2.5°) differences between the respective lattice directions of the original and reconstructed domains, implying that a small shear stress can enable these lattice reconstructions (Fig. 2B). In addition, the fresh crack edges provide ideal free surfaces for the nucleation and growth of these reconstructed domains.

Fig. 2 In situ STEM observation of mode I cracking in 1L ReS2.

(A) Snapshots of a crack propagation in 1L ReS2 and accompanying lattice reconstructions. Different domains are highlighted in red (reconstructed) and yellow (original). Scale bar, 1 nm. (B) DFT-obtained angle differences after lattice reconstructions along the twin boundary along the i direction. Only Re atoms are shown. “Diamond chain” directions (a direction) are highlighted by blue stripes. Separated domain structures are highlighted in yellow and purple. Open and closing states corresponding to mode I and mode III cracking. Other types of boundaries with variable angle changes of lattices are possible. (C) STEM image for the area near the crack tip showing lattice reconstructions stem from postcrack surfaces and modify the strain at the crack tip. Scale bar, 2 nm. (D) Typical geometric phase analysis (GPA) strain analysis (see Materials and Methods) result of normal strain perpendicular to the crack line for the mode I cracks in 1L ReS2. The local enhancement of open stress induced by the lattice reconstructed zone is highlighted. Scale bar, 2 nm. (E) Serial STEM images during the 30-s in situ observation (image positions aligned vertically); the crack extension is highlighted by white dashed lines. Lattice reconstructed zones are false colored. Scale bar, 1 nm.

Although the atomic structure is usually kept original inside mode I opening crack tips (within 5 nm to the crack tip), the lattice reconstruction around the crack tip zones will impactfully alter the stress/strain level in the crack tips. Particularly, although the local strains near the coherent twin boundaries are negligible, the lattice axis switching from a to b along the crack line has yielded ca. −1.6% normal strain perpendicular to the crack line for cracks along the a direction (Fig. 2C). That is, perpendicular to the crack line, the reconstructed domain virtually shrinks compared to the original and adds on extra tensile strains in the adjacent crack tip zone (Fig. 2D). This phenomenon contradicts the previous concept of plastic fracture, in which the plasticity mainly relieves the strain in crack tips. The origin should be attributed to the remarkable lattice reconstructions in such 2D materials due to the absence of out-of-basal-plane constraint.

As a result, we have observed continuous crack extension by in situ observations (Fig. 2E), which is primarily driven by such lattice reconstructions. It should be noted that these cracks are not subject to external loadings during our TEM experiments, and are not affected by the electron beam damage; the S atoms on the fresh crack edges, which are easiest to be knocked out by electron beam, still can be seen; and no loss of S atom is found (fig. S5); hence, the lattice reconstructions play a major role for the crack extension. Because lattice reconstructions themselves are generated and developed on fresh crack edges, the positive feedback on crack extension eventually leads to a kind of self-persistent fracture. Basically, the lattice direction switch from a to b during cracking causes the lattice spacing contraction perpendicular to the crack lines; hence, the mode I loading is enhanced (fig. S6). Crack extension further pushes the reconstructed domain boundaries moving forward and eventually converts the entire crack edge into new lattice orientations.

Mode I cracks along different crystal directions in 1L ReS2 are directly compared (Fig. 3A). The roughness of the crack edges is variable with different directions: a < b < i < other directions. a direction crack edges have the lowest roughness (also the lowest fracture toughness) because of the local open stress enhancement in mode I crack tip by lattice reconstruction as we have discussed above. More lattice reconstructions can be seen in cracks along directions other than the a direction (figs. S7 and S8). The introduced domain boundary and defects by lattice reconstructions can instantly modify the strain field around the crack tip zones and alter the crack modes and crack directions. Moreover, the emission of two partial dislocations near the crack tip has been captured (Fig. 3B). For the original mode I crack along the a direction, if the open stress is turned oblique midway, similar to our previous report (13), then full dislocations or partial dislocations could be emitted from the crack front and blunt the crack tips. The dislocations resided in the “plastic deformation zones” with less than 2-nm sizes.

Fig. 3 Crystal direction dependence and dislocations.

(A) Mode I cracks along different crystal directions in 1L ReS2. Central panel: Schematic crack paths along the different crystal directions of ReS2 (crystal direction indicated by the arrows and atomic models in the bottom left). Peripheral panels: Representative crack images corresponding to directions 1 to 6 in the central panel. All the images rotated in the same crystal direction according to the central. (B) Serial STEM images for the emission of two partial dislocations near the crack tip. Scale bar, 1 nm.

In the ex situ 1L ReS2 samples, we were also able to observe the anomalous straight cracks steadily propagating in parallel with atomically stacked postcrack edges in ReS2 samples (Fig. 4, A to C). After cracking, the fresh edges behind the crack front were immediately stacked by the van der Waals (vdW) forces so that persistent and steady cracking was achieved. The stacking order still follows the regular crystal planes, except for some point-like defective positions. High-resolution images were captured and analyzed for strain distributions at the crack tips (Fig. 4, D and E). Here, the major driving force was the in-plane shear stress (Fig. 4E)—originating from faulted stacking by half unit cell, suggesting that these cracks are mainly loaded by mode II type. The shear stress is proportional to the magnitude of mismatch in stacking, and the driving force for this crack was gradually reduced by leaving point-like “dislocations” in the crack path (Fig. 4C).

Fig. 4 Mode II crack in 1L ReS2.

(A) Scheme of the mode II crack in 1L ReS2. (B) Scheme for restacking of crack edges; gray area shows the bilayer stacked part, thicker lines in the middle show the crack edges, and dislocations are lattice faulted–stacked zones. To lower the energy cost and minimize the shear strain, lattice rotation is allowed, which is more energy favorable than lattice distortion, so the number of dislocations in the crack tip is proportional to the width of the overlapping (gray zone) during mode II crack propagation. (C) Typical mode II crack in ReS2; crack edges are stacked and thus have brighter contrast in HAADF. Driving forces (shear) are gradually consumed by emission of partial dislocations (marked). The bright spot (on the left) is caused by PMMA residue during specimen transfer. Scale bar, 2 nm. (D) Magnified STEM image corresponding to the dashed rectangle in (C). Scale bar, 1 nm. (E) In-plane shear strain analyzed by the GPA method on (D); the dashed circle highlights the crack front. Scale bar, 1 nm. (F and G) STEM images for the crack paths and crack edges with different overlapping widths. Scale bars, 1 nm (F) and 2 nm (G).

The reason for restacking of postcrack surfaces (edges) is mechanical. Shear strain under mode II loads can be in part released by lattice rotation if restacked (Fig. 4B), and the stack widths of crack edges are proportional to the mismatch of lattices along crack lines. Therefore, the driving shear stress level can be straightforwardly determined by the restacked widths. The real shear stress during these mode II cracks on growth substrates is higher than currently measured. It should be noted that the shear strain is applied by the growth substrate (sapphire) during the CVD cooling down process, and the crack tip strain measured by STEM images has been partially released after the samples are exfoliated from substrates and transferred to the TEM grids.

Normally, the fresh crack surfaces are not allowed to overlap because of the space exclusion concept for opposite crack surfaces in bulk materials; however, in 2D materials, the ultrathin thickness allows the vertical overlapping, and the flexibility of 2D materials significantly lowers the bending energy cost in the crack tip zone. The continuous cracks in mode I, mode II, and even mixed mode I and II can be connected and form turning and loops in continuous fracture (Fig. 1B).

These mode II cracks can make parallel atomic stacking with variable widths, from 0.5 to 15 unit cells (Fig. 4, F and G, and fig. S9). As mentioned, the stack width depends on the shear stress/strain magnitude. Analogous to the mode I cracks above, the lattice reconstructions could sometimes change the lattice orientations in opposite crack edges (Fig. 4, F and G). Although the stacking is mediated by the vdW forces over small areas, it is quite robust against thermal vibration (in micrometer-scale suspended sample; see Fig. 1B) and wet transfer process by the PMMA method, and even survives intensive electron beam irradiation during our in situ TEM observations. Electron tunneling governs the electrical transport through such vertically stacked atomic layers and has direct dependence on the contact atomic configuration and area (width) (26); in this instance, the atomic stacking by controllable mode II cracking can be a general preparation method for new tunneling devices in 2D materials.

Furthermore, the mode III cracks in 2D materials (Fig. 5A) have been identified in our experiments as well. By in situ TEM, occasionally the out-of-plane shear could be supplied by the bended crack edges (Fig. 5B). The out-of-plane tilting in the samples can be corroborated with in-plane shear angle/strain measured near the crack tip (Fig. 5, C and D). Because of basal plane tilting (Fig. 5A), the open stress as in mode I loading is changed to closing stress in mode III loading (Fig. 5E), and the lattice switching direction of the a/b axis is correspondingly changed (Fig. 5F); hence, the mode I and mode III cracks can be clearly distinguished.

Fig. 5 Mode III crack in 1L ReS2.

(A) Scheme of mode III crack; the original basal (pink) plane is tilted (blue) after the crack. Right panel shows the original and tilted basal plane near mode III crack tip. (B) Low-magnification image showing the origin of out-of-plane shear. Scale bar, 100 nm. (C) Example of mode III crack. Scale bar, 1 nm. (D) The shear strain analysis by GPA method on (C), showing the tilting in the right side of crack. Scale bar, 1 nm. (E) Scheme of mode III crack viewing in tilted basal plane [blue plane in (A)] near crack edges, showing closing stress near edge. (F) One mode III crack in 1L ReS2 viewing in original basal plane [pink plane in (A)], lattice reconstruction, and domain boundary induced by closing stress. Scale bar, 1 nm.

The anomalous fracture behavior in 2D ReS2 can be seen in the above easy twinning process during fracture compared with 3D bulk materials, and also in the differences of mode II and mode III fractures between 2D and 3D bulk materials. For the qualitative evaluation on the trend of twining in 1L ReS2, we have applied the DFT calculations (see Materials and Methods) to estimate the onset twinning stress in 2D ReS2 as low as 4 MPa (for the twinning structures as shown in Fig. 2B), which is much lower than the twinning stresses of conventional bulk hexagonal and cubic metals that are usually between 30 and 200 MPa (27).

According to our TEM observations presented above and in previous works (13, 28), the atomic defects such as vacancies in 2D TMD materials (such as MoS2 and ReS2) will facilitate the dislocation dynamics in crack tips; thus, it can enhance the fracture toughness of these 2D materials. However, it should also be noted that atomic vacancies in 2D materials can expediate the fracture themselves. Second, while the low-angle grain boundaries can be a possible way to enhance the strength, they are usually not directly enrolled in the fracture path. Third, for the anisotropic 2D materials such as ReS2, the formation of coherent twinning domains might reduce the fracture toughness and make atomically smooth crack edges, but it can be avoided by introducing atomic vacancies that might trigger dislocation emissions in crack tips. Overall, properly controlled addition of atomic vacancies in 2D TMD materials can be one feasible choice for toughening the 2D TMD materials.


The fracture in 2D ReS2, including mode I, II, and III cracks, have been thoroughly investigated down to the atomic scale by using direct TEM observations. As compared with 3D bulks, we found a much easier twinning in anisotropic 2D materials, and it has led to some interesting new insights near the crack tip zones. Our study revealed the notable lattice reconstructions around the crack tips, which significantly change the loading stress in crack tips. In contrast to the ordinary plastic deformation in crack tip zones, the lattice reconstructions in 2D ReS2 may conversely enhance the stress intensity for the corresponding fracture modes, increasing the brittleness and rendering self-persistence in cracking. In addition, the anomalous mode II cracks in 2D materials that make atomic parallel restacked crack edges have been discovered. Owing to the ultimate thin thickness and entirely free surfaces of 2D materials, new peculiar fracture behaviors are made possible. The present work shed light on the future mechanical engineering on 2D materials, facilitating the high strength and high toughness as well as the novel functionalities in 2D materials.


Synthesis of ReS2

Monolayer (1L) ReS2 was synthesized on polished c-face sapphire substrate in atmospheric pressure using two-zone tube furnace. The precursors with weight ratio 1:50 of ammonium perrhenate (NH4ReO4) (99.999%; Sigma-Aldrich) and sulfur (99.998%; Sigma-Aldrich) were used as received. Sulfur source was put in the upstream zone, while Re source was in the downstream zone. The sapphire substrate was placed on top of the Re source with the polished face upside down. Upstream and downstream zones were both heated up to 200° and 850°C at the same time, respectively, and held for 10 min. Before heating, the tube was flushed by 300 standard cubic centimeter per minute (sccm) of Ar for 10 min, and afterward, Ar flow was adjusted to 80 sccm and kept constantly in the whole synthesis process. The triclinic and anisotropic ambient stable Td phase of ReS2 stems from the typical tetragonal (T) phase in TMDs with “diamond”-type superlattices in two directions (labeled as a and b) and slightly distorted. Figure S1 shows the major characterization results of our 1L ReS2 membrane.

(S)TEM specimen preparation

By the PMMA-assisted technique, the CVD-grown ReS2 was transferred onto a TEM grid. Initially, thin-layer PMMA was spin coated on the as-grown sapphire substrate at 3000 revolutions per minute for 50 s. With the aid of PMMA, the ReS2 detached from the substrate by immersing in 75°C deionized water for 2 hours. Next, the PMMA/ReS2 layer was transferred onto a TEM grid and dried at ambient temperature. Acetone vapor was introduced for gentle removal of PMMA film.

(S)TEM characterizations

Aberration-corrected STEM under 60-kV accelerating voltage was applied to improve the image resolution and prevent beam damage. To prevent knock-on damage on ReS2, 60 kV energy lower than threshold for ReS2 was used. A JEM-ARM200F transmission electron microscope with a Corrected Electron Optical Systems (CEOS) probe spherical (Cs) aberration corrector was used. The vacuum value during measurement was around 1.3 × 10−7 mbar, together with electron beam current of 13.1 μA. The beam size for STEM was approximately 1.5 Å. The camera length of STEM was 120 mm during imaging. The defocus was −4 nm, and acquisition time of the high-angle annular dark field (HAADF) image was 19 μs per pixel to minimize damage and obtain the figures with lower drift. The 256 × 256 or 512 × 512 pixel images were acquired with condenser lens aperture of 40 μm, and the range of collection angle was 45 to 180 mrad. In postprocessing, Wiener filtering was applied on HAADF images for reduction of noises.

In situ STEM on cracking

For creating initial precracks in 1L ReS2 and preventing beam damage on the ReS2 sample areas, the crack initialization zones (by intensive beam irradiation) and in situ cracking observation zones (with ultralow electron dose) are carefully separated. The electron beam sputtering used for generating the initial cracks was convergent beam with continuous electron current. High beam intensities (greater than 0.3 pA/nm2) were exposed on controlled initialization areas (~100-nm-diameter circles) of 1L ReS2 for 15 to 30 min. Then the STEM-HAADF imaging was carried out outside of the beam-exposed zones. The initialization of crack was carried out by converging the electron beam in TEM mode at 60 kV; under an intensive electron dose, it is around 100 times higher than the dose during the crack observation in the following STEM HAADF imaging. The electron beam used for the STEM observation was short time exposure. The time interval between the serial in situ STEM images was set as 1 to 20 s. The electron beam effect has been minimized during observation around crack tips, and all our observations have been selected on the crack tip zones, which are far away from the initial precrack generation zones. No damage on the S atoms on edges, which should be most likely to be sputtered, or inside the grains has been found.

Strain analysis on TEM images

The GPA strain analysis (29) (on high-resolution STEM-HAADF images) was performed with the reflexes (in reciprocal space) perpendicular to the a and b directions as the two bases, respectively. The GPA analysis resolution was set as 0.3 to 0.4 nm, and smoothing factor was set as 7.0.

DFT calculations

The simulations of 1L ReS2 atomic structures and twin domain structures were carried out by using spin-polarized DFT calculations. The Vienna ab initio simulation package (VASP) program package was used with exchange-correlation interactions described by using Perdew-Burke-Ernzerhof functional with a gradient approximation (3032). The distance of the vacuum layer was set to be more than 15 Å, which is adequately large to avoid interlayer interactions. Kinetic energy cutoff for the plane-wave basis set was 400 eV, and the electronic SCF tolerance was set to 10−4 eV. The k-point samplings are 3 × 1 × 1 in the Brillouin zone for structural optimizations. The size of the DFT supercell is 2.74 nm width (x axis) by 4.00 nm length (y axis), and the z direction is 15 Å, which can provide enough space to avoid the interactions by layers. Fully relaxed geometries were obtained by optimizing all atomic positions until the Hellmann-Feynman forces are less than 0.05 eV/Å. According to the Peierls framework for twinning (33), the critical stress (τ) for twinning can be estimated by τ = (Ψ(twinned) − Ψ(pristine))/(∆x (twinned) − ∆x (pristine)), where τ is the shear stress along the slip planes, ∆x is the share displacement of the central pair of planes, and Ψ is the energy of sheared configurations, as determined by DFT calculations.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: This work was supported by National Natural Science Foundation of China (project nos. 51872248, 21703076, and 51922113); the Hong Kong Research Grant Council under the Collaborative Research Fund (project no. C6021-14EF), General Research Fund (project no. 15302419), and Early Career Scheme (project nos. 25301018 and 21303218); the City University of Hong Kong (project no. 9610387); The Hong Kong Polytechnic University (1-ZVGH); the Natural Science Foundation of Jiangsu Province of China (project no. BK20170466); the Natural Science Foundation of the Jiangsu Higher Education Institutions of China (no. 18KJA140001); the Shenzhen Science and Technology Innovation Commission (project no. JCYJ20170818104717087); and Fund of Key Laboratory of Advanced Materials of Ministry of Education. Author contributions: J.Z. and T.H.L. conceived and led the project. L.H., F.Z., L.W.W., Q.H.T., C.-S.L., S.P.L., T.H.L., and J.Z. prepared the samples and conducted the sample characterizations and TEM experiments. Q.D. conducted the DFT simulations. Y.C. and N.W. assisted with the TEM experiments. All the authors discussed the manuscript (MS) and agreed on its contents. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article