# NEW FRONTIERS IN MULTISCALE MODELLING OF ADVANCED MATERIALS

EDITED BY: Simone Taioli, Maurizio Dapor and Nicola M. Pugno PUBLISHED IN: Frontiers in Materials

#### *Frontiers Copyright Statement*

*© Copyright 2007-2016 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.*

*The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.*

*Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.*

*Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.*

*As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.*

> *All copyright, and all rights therein, are protected by national and international copyright laws.*

*The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use.*

ISSN 1664-8714 ISBN 978-2-88919-755-2 DOI 10.3389/978-2-88919-755-2

# About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

# Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

# Dedication to quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

# What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# **NEW FRONTIERS IN MULTISCALE MODELLING OF ADVANCED MATERIALS**

## Topic Editors:

**Simone Taioli,** European Centre for Theoretical Studies in Nuclear Physics and Related Areas & Bruno Kessler Foundation, Italy, Charles University in Prague, Czech Republic **Maurizio Dapor,** European Centre for Theoretical Studies in Nuclear Physics and Related Areas & Bruno Kessler Foundation, Italy

**Nicola M. Pugno,** University of Trento, Italy, Queen Mary University of London, UK, Center for Materials and Microsystems, Bruno Kessler Foundation, Italy

Tight-binding model of fullerene impact with a surface using superposition of atomic wave functions and Slater-Koster interatomic matrix elements allowing the simulation of larger molecular structures.

Image taken from: Aversa L, Taioli S, Nardi MV, Tatti R, Verucchi R and Iannotta S (2015) The interaction of C60 on Si(111) 7 × 7 studied by supersonic molecular beams: interplay between precursor kinetic energy and substrate temperature in surface activated processes. Front. Mater. 2:46. doi: 10.3389/fmats.2015.00046

Atomistic simulations, based on ab-initio and semi-empirical approaches, are nowadays widespread in many areas of physics, chemistry and, more recently, biology. Improved algorithms and increased computational power widened the areas of application of these computational methods to extended materials of technological interest, in particular allowing unprecedented access to the first-principles investigation of their electronic, optical, thermodynamical and mechanical properties, even where experiments are not available.

However, for a big impact on the society, this rapidly growing field of computational approaches to materials science has to face the unfavourable scaling with the system size, and to beat the time-scale bottleneck. Indeed, many phenomena, such as crystal growth or protein folding for example, occur in a space/time scale which is normally out of reach of present simulations. Multi-scale approaches try to combine different scale algorithms along with matching procedures in order to bridge the gap between first-principles and continuum-level simulations.

This Research Topic aims at the description of recent advances and applications in these two emerging fields of ab-inito and multi-scale materials modelling for both ground and excited states. A variety of theoretical and computational techniques are included along with the application of these methods to systems at increasing level of complexity, from nano to micro. Crossing the borders between several computational, theoretical and experimental techniques, this Research Topic aims to be of interest to a broad community, including experimental and theoretical physicists, chemists and engineers interested in materials research in a broad sense.

(a) Atomic orbitals (wave functions) and DFT model of Buckminsterfullerene (C\_60). (b) ab-initio simulation of the growth of carbon-based 2D materials via supersonic impact of C\_60 molecule on Si(111) (Image taken from: Aversa L, Taioli S, Nardi MV, Tatti R, Verucchi R and Iannotta S (2015) The interaction of C60 on Si(111) 7 × 7 studied by supersonic molecular beams: interplay between precursor kinetic energy and substrate temperature in surface activated processes. Front. Mater. 2:46. doi: 10.3389/fmats.2015.00046) and corresponding molecular dynamics counterpart for the simulation of larger systems. (c) Impact simulations on macroscopic multilayer composite materials (Image Signetti S and Pugno NM (2015) Frontiers in modeling and design of bio-inspired armors. Front. Mater. 2:17. doi: 10.3389/fmats.2015.00017) whose mechanical properties can be derived from lattice-spring model mesoscale simulations, e.g., on graphene and other 2D-materials reinforced nanocomposites (Image from Brely L, Bosia F and Pugno NM (2015) A hierarchical lattice spring model to simulate the mechanics of 2-D materials-based composites. Front. Mater. 2:51. doi: 10.3389/fmats.2015.00051).

### Acknowledgments:

N.M.P. is supported by the European Research Council (ERC StG Ideas 2011 BIHSNAM n. 279985 on "BioInspired hierarchical supernanomaterials", ERC PoC 2013-1 REPLICA2 n. 619448 on "Large-area replication of biological antiadhesive nanosurfaces", ERC PoC 2013-2 KNOTOUGH n. 632277 on Supertough knotted fibers), by the European Commission under the Graphene Flag-ship (WP10 "Nanocomposites," n. 604391), and by the Provincia Autonoma di Trento ("Graphene Nanocomposites," n. S116/2012-242637 and delib. reg. n. 2266).

**Citation:** Taioli, S., Dapor, M., Pugno, N. M., eds. (2016). New Frontiers in Multiscale Modelling of Advanced Materials. Lausanne: Frontiers Media. doi: 10.3389/978-2-88919-755-2

# Table of Contents


Lucrezia Aversa, Simone Taioli, Marco Vittorio Nardi, Roberta Tatti, Roberto Verucchi and Salvatore Iannotta


Maurizio Dapor

*37 Hydrogen storage in rippled graphene: perspectives from multi-scale simulations*

Vito Dario Camiola, Riccardo Farchioni, Tommaso Cavallucci, Antonio Rossi, Vittorio Pellegrini and Valentina Tozzini

*40 Emergent models for artificial light-harvesting*

Celestino Creatore, Alex W. Chin, Michael A. Parker and Stephen Emmott

*44 Atomistic model of metal nanocrystals with line defects: contribution to diffraction line profile*

Alberto Leonardi and Paolo Scardi

*54 Molecular dynamics study on the distributed plasticity of penta-twinned silver nanowires*

Sangryun Lee and Seunghwa Ryu

*61 Frontiers in modeling and design of bio-inspired armors* Stefano Signetti and Nicola Maria Pugno

*64 A hierarchical lattice spring model to simulate the mechanics of 2-D materials-based composites*

Lucas Brely, Federico Bosia and Nicola M. Pugno


# Editorial: New Frontiers in Multiscale Modelling of advanced Materials

*Simone Taioli1,2\*, Maurizio Dapor1 \* and Nicola M. Pugno3,4,5\**

*1European Centre for Theoretical Studies in Nuclear Physics and Related Areas (ECT\*), Bruno Kessler Foundation, Trento, Italy, 2 Faculty of Mathematics and Physics, Charles University in Prague, Prague, Czech Republic, 3Department of Civil, Environmental and Mechanical Engineering, University of Trento, Trento, Italy, 4School of Engineering and Materials Science, Queen Mary University of London, London, UK, 5Center for Materials and Microsystems, Bruno Kessler Foundation, Trento, Italy*

Keywords: multiscale modelling, *ab initio* and density-functional theory quantum methods, molecular dynamics simulations, electronic and mechanical properties of solids, materials growth, graphene, materials characterization, Monte Carlo

Computer simulations, based on *ab initio*, semi-empirical, or continuum approaches, are now a widespread tool for benchmarking experimental measurements in many areas of engineering, physics, chemistry, and, more recently, biology. Indeed, improved algorithms and increased computational power have widened the areas of application of these computational methods to materials of technological interest, allowing unprecedented access to the investigation of their electronic, optical, thermodynamical, and mechanical properties.

*Edited and Reviewed by: John L. Provis, University of Sheffield, UK*

#### *\*Correspondence:*

*Simone Taioli taioli@ectstar.eu; Maurizio Dapor dapor@fbk.eu; Nicola M. Pugno Nicola.Pugno@unitn.it*

#### *Specialty section:*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials*

*Received: 17 November 2015 Accepted: 27 November 2015 Published: 14 December 2015*

#### *Citation:*

*Taioli S, Dapor M and Pugno NM (2015) Editorial: New Frontiers in Multiscale Modelling of Advanced Materials. Front. Mater. 2:71. doi: 10.3389/fmats.2015.00071*

This "third way" represented by the computational lab lies in between the use of pure theoretical models and experiments, and it can be particularly useful when experimental data are difficult to be recorded and theoretical models are too complex to be formulated. In this regard, these methods have been successfully used to investigate and predict properties of matter, ranging from the nanoscale behavior of crystalline materials to their interaction with external fields, e.g., electromagnetic and elastic, and to their macroscopic behavior.

However, they have been commonly adopted via a "single-scale" approach, in which one uses one specific method at each relevant scale. Indeed, this is the simplest approach to face the unfavorable scaling with system size and to beat the time-scale bottleneck. For example, on the one hand, *ab initio* methods, such as density functional theory, have been used to investigate accurately the nanoscale behavior of crystalline materials starting from the first-principles simulation of the electronic motion. On the other hand, coarse-grained methods and discrete mesoscopic models are adopted to assess the properties of materials at atomistic or continuum level. In the latter case, classical molecular dynamics, in connection with fiber bundle or lattice models, numerical and analytical homogenization, and upscaling techniques, as well as finite element or mesh-free approaches, have been successfully used to determine the constitutive laws of materials. Furthermore, the application of these techniques to macroscopic systems enables one to investigate mechanisms, such as fracture, crack, and grain boundary propagation.

Nevertheless, the "single-scale approach" results in a loss of information across different scales and is an important limiting factor in the science of complex materials. A comprehensive vision would allow the exploration of full ranges of materials solutions and hierarchical architectures at all relevant size scales, with the idea that the resulting product is not simply the sum of all the parts. This requires development of new science as well as novel computational methods, since multiscale tools are still in their infancy.

This collection of papers is thus devoted to this rapidly growing field of computational approaches to multiscale methods by reporting a specific set of problems related, but not limited to, materials science. Thus, the reader will find the application and point-of-views on a variety of state-of-the-art methods spanning different scales and combining different approaches to bridge the gap between first-principles and continuum-level simulations by matching procedures, both concurrent and hierarchical.

The ambition of this collection is to show how electronic, atomistic, mesoscopic, and continuum models can be combined and integrated to accelerate materials discovery, to shorten the time from "lab-to-market" from the current 10–20 years by avoiding lengthy trial-and-error loops or overlooking other materials not yet discovered for the absence of a multiscale, multiphysics vision.

Crossing the borders between several computational, theoretical, and experimental techniques, this collection should be of interest to a broad community, including experimental and theoretical physicists, chemists, and engineers interested in materials modelling in a broad sense.

# AUTHOR CONTRIBUTIONS

This contribution is an Editorial. As such the Topic Editors equally contributed to this introduction.

# FUNDING

The organizers acknowledge the European Centre for Theoretical Studies in Nuclear Physics and Related Areas (ECT\*) in Trento (Italy) for hosting the workshop "New Frontiers in Multiscale Modelling of Advanced Materials," 17–20 June, 2014. ECT\* is sponsored by the Fondazione Bruno Kessler (FBK) in collaboration with the "Assessorato alla Cultura" (Provincia Autonoma di Trento), funding agencies of EU Member and Associated States. ECT\* also receives support from various instruments of the Framework Programmes of the European Commission. NP is supported by the European Research Council (ERC StG Ideas 2011 BIHSNAM no. 279985 on "BioInspired hierarchical supernanomaterials," ERC PoC 2013-1 REPLICA2 no. 619448 on "Large-area replication of biological antiadhesive nanosurfaces," ERC PoC 2013-2 KNOTOUGH no. 632277 on Supertough knotted fibers), by the European Commission under the Graphene Flag-ship (WP10 "Nanocomposites," no. 604391), and by the Provincia Autonoma di Trento ("Graphene Nanocomposites," no. S116/2012-242637 and delib. reg. no. 2266).

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Taioli, Dapor and Pugno. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Ab initio molecular dynamics with quantum Monte Carlo

# **Ye Luo\* and Sandro Sorella**

International School for Advanced Studies (SISSA) and CRS Democritos, CNR – Istituto Nazionale per la Fisica della Materia (INFM), Trieste, Italy \*Correspondence: xw111luoye@gmail.com

#### **Edited by:**

Simone Taioli, Bruno Kessler Foundation, Italy; Charles University, Czech Republic

**Reviewed by:** Ethan Brown, ETH Zurich, Switzerland Francesco Pederiva, University of Trento, Italy

**Keywords: quantum Monte Carlo, ab initio calculations, molecular dynamics, many-body physics, high performance computing**

Quantum Monte Carlo (QMC) is one of the most accurate electronic structure methods for *ab initio* many-body calculations in a broad range of electronic systems. Thanks to fast evolving massively parallel computers, much larger simulations of realistic systems become affordable. We have developed an *ab initio* simulation scheme at finite temperature based on molecular dynamics (MD) and QMC during the past few years. In this approach, a second order Langevin MD is employed by using a statistical evaluation of the forces acting on each atom by means of QMC. Moreover, the corresponding statistical noise is also used to drive and accelerate the dynamics of ions. The accuracy and the reliability of our *ab initio* MD have been studied systematically and it has been already successfully applied to study the molecular and atomic phases of liquid hydrogen under pressure and liquid water at ambient conditions. Since this *ab initio* method provides a better resolution of the thermodynamic properties of materials, we believe that it represents a very promising tool for the most challenging applications in physics, chemistry, and biology.

Since the Car-Parrinello molecular dynamics (Car and Parrinello, 1985) was proposed three decades ago, *ab initio* molecular dynamics (MD) based on the density functional theory (DFT) has been widely accepted as an accurate and powerful tool to study the thermodynamic properties of systems at ambient and extreme conditions in the fields of physics, chemistry, and biology. Even though DFT is in principle exact, its accuracy, in practice, is affected by the quality of the approximation used in the exchange correlation functional term whose universal exact expression is not accessible, as the simplest approximate functionals, such as LDA, PBE, BLYP, and several other ones, so far cannot be systematically improved in a computationally efficient way. Quantum Monte Carlo (QMC) is the most promising successor of DFT in *ab initio* molecular dynamics. Indeed, the scaling behavior of QMC with the number of electrons is very good, from the second power (Ceperley et al., 1977) for thermodynamic properties – e.g., energy per electron, because in this case the number of samples to obtain a fixed accuracy in this quantity can be decreased when increasing the number *N*<sup>e</sup> of electrons, implying *N*<sup>3</sup> e /*N*<sup>e</sup> = *N*<sup>2</sup> e scaling with the standard algorithm – to the fourth power of the number of electrons when chemical accuracy in the total energy is required – in this case, instead the number of samples required for fixed accuracy has to increase with *N*e, yielding *N*<sup>3</sup> <sup>e</sup> × *N*<sup>e</sup> = *N*<sup>4</sup> e scaling. For this reason (a cost similar to DFT, albeit with a considerably larger pre-factor), QMC is also by far more efficient compared with other quantum chemistry methods based on multi-determinant expansions or many-body perturbation theory around the Hartree–Fock, such as, for instance, configuration interaction and coupled cluster theory, respectively. These methods provide a similar or even better accuracy than QMC for systems consisting of a few atoms, but they become computationally prohibitive when the system sizes increase. Another very important reason favoring the QMC algorithm is its well established efficient parallelization in current and upcoming supercomputers. QMC codes can be easily adapted to millions of CPU cores while DFT codes are still struggling with its parallel performance over several thousand cores. For the above

reasons, it is probably not too far from reality to claim that QMC represents the rising star in the *ab initio* MD simulations.

This opinion paper focuses on the construction of the QMC-based MD by showing the key ingredients and the main recent developments of the building blocks necessary for the efficient implementation of this method. There are four main components: the basics of QMC, the wavefunction optimization, the evaluation of forces, and the molecular dynamics scheme. The last three ingredients were efficiently developed within only 10 years, with substantial contribution given also by us (Casula et al., 2004; Sorella et al., 2007; Umrigar et al., 2007; Attaccalite and Sorella, 2008; Sorella and Capriotti, 2010; Luo et al., 2014; Mazzola et al., 2014; Zen et al., 2014). Thus, the QMC-based MD is no longer a dream.

## **BASICS OF QMC**

In QMC, the wavefunction ansatz is of fundamental importance for the efficiency, reliability, and accuracy of both variational and diffusion Monte Carlo calculations (VMC/DMC), namely the two most popular methods in QMC, the latter being more accurate and based on an approximate evaluation of the ground state properties with the restriction of having the same phases of the best reference VMC wavefunction – i.e., the so-called fixed-node approximation (Anderson, 1976; Ceperley and Alder, 1984; Umrigar et al., 1993), which is unavoidable for fermionic systems. The repeated evaluation of the wavefunction with also its gradients and Laplacian is one of the most consuming parts of the computation; the optimal wavefunction ansatz for MD should be also compact.

The most popular choice is the Jastrow factor (Jastrow, 1955) correlated Slater determinant, which is simple and generally accurate. Pseudopotentials (Trail and Needs, 2005; Burkatzki et al., 2007) are often employed like in DFT to replace chemically in-active core electrons and keep only the valence electrons in QMC calculations.

#### **WAVEFUNCTION OPTIMIZATION**

In the MD simulation, after each ionic move, the parameters in the VMC wavefunction need to be optimized in order to satisfy the Born–Oppenheimer approximation. Efficient optimization methods ensure that the energy converges to the right minimum within a sufficient accuracy in just a few steps. At present, the variance minimization method (Umrigar et al., 1988) has been recently abandoned because sometimes it failed to yield the best possible wavefunction (Snajdr and Rothstein, 2000; Bressanini et al., 2002). This method has been replaced by much more efficient and robust energy minimization methods. In our previous works, we introduced an energy minimization method, called stochastic reconfiguration (SR), which is based only on energy first derivatives, and allows a much more efficient optimization compared with other more conventional methods of this type, such as the steepest decent. In addition, the trial step length can be automatically determined by the SR with the Hessian accelerator (SRH) (Sorella, 2005) method, which makes use of the second derivatives of the energy. This method can be further improved by the recently introduced "linear method" (Umrigar and Filippi, 2005; Umrigar et al., 2007), which provides a better sampling of the Hessian matrix and represents at present the state of the art for the QMC optimization of wavefunction in realistic systems. For the application of the above techniques, to large systems, it is important to perform in a very efficient way the expensive inversion operation of a big matrix with its leading dimension (10<sup>4</sup> is possible on current supercomputers) equal to the number of variational parameters (Sorella et al., 2007; Neuscamman et al., 2012).

#### **EVALUATION OF FORCES**

In order to study thermodynamic properties, classical Monte Carlo and molecular dynamics are both widely used simulation methods. The former has a significant advantage that it does not necessarily require the computation of forces. Unfortunately, it is not straightforward to combine classical Monte Carlo and quantum Monte Carlo because all the quantities evaluated by QMC, including the energy, carry statistical errors that are very expensive to eliminate.<sup>1</sup> In order to move the ions according to the QMC energy and sample the canonical ensemble correctly, the penalty method (Ceperley and Dewing, 1999) was introduced to remove the bias due to the statistical uncertainty in the knowledge of the energy by rejecting the proposed moves more frequently than the standard Metropolis algorithm. For this reason, this method is very expensive, particularly in the low temperature regime, because of too many rejected moves, and so far applications have been limited to hydrogen with up to 54 protons in this regime (Pierleoni and Ceperley, 2006; Morales et al., 2010a,b).

Generally speaking, it is clear that, when the computational cost for the calculation of the nuclear forces is comparable to that of the energy, MD should be more efficient, because with the same cost all the atoms are moved without any rejection. For instance, in DFT, where the forces are obtained almost for free by applying the Hellmann–Feynman theorem, MD is a common practice to sample the canonical distribution, and, to our knowledge, only hybrid methods based on Monte Carlo and molecular dynamics (Rossky et al., 1978; Duane et al., 1987) can be competitive.

Unfortunately, computing forces with QMC is substantially more demanding than DFT. To obtain the forces in an accurate and efficient way within the VMC framework, several very important improvements have been done, concerning both theories and algorithms. The correlated sampling (Filippi and Umrigar, 2000) method allows the computation of the energy derivatives with errors much smaller than those obtained with a straightforward finite difference method. Together with the space warp coordinate transformation (Umrigar, 1989), forces computed by QMC were used for structural optimization (Casula et al., 2004). Another recent and important development in QMC was the solutions (Assaraf and Caffarel, 2000, 2003; Attaccalite and Sorella, 2008; Sorella and Capriotti, 2010; Zen et al., 2013) of the intrinsic infinite variance problem occurring in the calculation of nuclear forces in the simplest VMC scheme. Moreover, thanks to the algorithmic differentiation algorithm (Sorella and Capriotti, 2010), the cost of computing all the force components in a system can be afforded with a computational time at most a factor four larger than the one corresponding to the energy. This progress has led to several works, where structural optimization and highly accurate evaluations of the equilibrium configurations as well as related properties were possible even for quite large systems (Barborini and Guidoni, 2012; Coccia et al., 2013, 2014; Guareschi and Filippi, 2013; Zen et al., 2013).

#### **MOLECULAR DYNAMICS WITH QMC**

Even though the efficient wavefunction optimization and fast evaluation of forces are achieved, the application of QMC for MD simulation of materials remains challenging due to the theoretical difficulties in applying the Newton's equations of motion when the forces are given with a statistical uncertainty. For instance, the basic law of energy conservation cannot be met at all, when the forces are not exactly given at each step. Since it is prohibitive to eliminate the statistical error (for the same reason, see footnote 1), MD methods that allow the presence of random noise in the force components are certainly more suitable for QMC.

In our works, we have adopted a second order Langevin dynamics (SLD) for the sampling of the ionic configurations within the ground state Born– Oppenheimer approach. In this scheme, the statistical noise corresponding to the forces, is used to drive the dynamics at finite temperature by means of an appropriate generalized Langevin dynamics (Attaccalite and Sorella, 2008), while, during the

<sup>1</sup>The statistical error decays as 1/ √ *N*, where *N* is the number of uncorrelated samples.

dynamics, a noise correction is added to compensate the extra energy pumped into the system by the noise of the QMC forces. Later works also introduced two integration methods (Luo et al., 2014; Mazzola et al., 2014), which allowed us to work with noisy forces and large time steps, thus reducing the computational cost by a substantial fraction. The price to pay is the limitation of this approach to the evaluation of static quantities at finite temperature, because dynamical quantities, such as the diffusion constant, cannot be evaluated with the above stochastic implementation of SLD. Similar approaches (Krajewski and Parrinello, 2006; Kühne et al., 2007) have been proposed where an SLD algorithm has been devised also at the DFT level.

Last but not least, our proposed MD scheme with QMC also includes an automatic optimal tuning for the damping of the slow and the fast modes, which greatly improves the efficiency of the sampling. In the SLD equations, the damping coefficients matrix can be arbitrarily chosen to speed up the sampling, provided the fluctuation–dissipation theorem is satisfied by adding appropriate correlated noise to the force components. In our scheme, we have defined this correlation matrix in terms of the one corresponding to the noisy forces obtained with QMC. This is because we have empirically verified that the latter matrix is almost proportional to the dynamical matrix (see **Figure 1**).

Therefore, the energy modes with higher frequency vibrations can be systematically damped much more than the slow modes, and this clearly allows a faster propagation with larger time steps.

The VMC-based MD has been successfully applied to various systems from small molecules (Luo et al., 2014) to liquids, for instance, the liquid hydrogen at high pressure (Attaccalite and Sorella, 2008; Mazzola et al., 2014) and the liquid water at ambient conditions (Zen et al., 2014) for the properties of vibrational frequencies (Luo et al., 2014), atomic structures (Zen et al., 2014), and phase diagrams (Mazzola et al., 2014). It is implemented in the TurboRVB QMC package (Sorella, 2015) with a hybrid MPI and OpenMP paradigm and its performance in the production run reaches peta-scale, as it has been successfully tested up to 33K cores on BG/Q and 264K cores on the K-computer.

#### **CONCLUSION AND OUTLOOKS**

From the high performance computing perspective, the four main components that we have briefly reviewed in this work show different features in computational intensity and communications. The statistical calculation for the energy and its derivatives required by the optimization and forces puts very heavy loads on the CPU but only needs very little communication at the end of a big piece of work, which is an ideal pattern for the parallelization. In addition, each optimization step involves inversion operations of big matrices. It can be a potential bottleneck if it is solved in parallel on a network with high latency. The computational and communication cost of updating ion positions is negligible compared to the other three parts for tens and hundreds of atoms.

In conclusion, we have briefly described the main computational achievements that have allowed recently the first *ab initio* simulations of realistic materials with an efficient molecular dynamics based on accurate QMC forces. We believe that our work represents a first step toward an efficient tool for *ab initio* simulations of material properties much more reliable than DFT. Indeed, our evaluation of forces is limited to the VMC approach, and we expect some progress can be done by adopting the more reliable DMC method. Similarly, it should be possible to extend the method to more accurate wavefunctions, containing, for instance, many-body and backflow correlations, as well as implicit multi-determinant wavefunctions, such as the antisymmetrized geminal power (AGP) and the Pfaffian wavefunction. For all the above reasons, we remark once again our belief that QMC should be considered as the raising star for *ab initio* simulations of materials.

### **REFERENCES**


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 11 February 2015; paper pending published: 26 February 2015; accepted: 22 March 2015; published online: 07 April 2015.*

*Citation: Luo Y and Sorella S (2015) Ab initio molecular dynamics with quantum Monte Carlo. Front. Mater. 2:29. doi: 10.3389/fmats.2015.00029*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials.*

*Copyright © 2015 Luo and Sorella. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **The interaction of C<sup>60</sup> on Si(111) 7** *×* **7 studied by supersonic molecular beams: interplay between precursor kinetic energy and substrate temperature in surface activated processes**

#### *Lucrezia Aversa<sup>1</sup> \*, Simone Taioli 2,3,4, Marco Vittorio Nardi <sup>5</sup> , Roberta Tatti <sup>1</sup> , Roberto Verucchi <sup>1</sup> and Salvatore Iannotta<sup>6</sup>*

*1 Institute of Materials for Electronics and Magnetism (IMEM), National Research Council (CNR), Trento, Italy, <sup>2</sup> European Centre for Theoretical Studies in Nuclear Physics and Related Areas (ECT\*), Bruno Kessler Foundation (FBK), Trento, Italy, <sup>3</sup> Trento Institute for Fundamental Physics and Applications (TIFPA), National Institute for Nuclear Physics (INFN), Trento, Italy, <sup>4</sup> Faculty of Mathematics and Physics, Charles University, Prague, Czech Republic, <sup>5</sup> Institut fur Physik, Humboldt-Universitat zu Berlin, Berlin, Germany, <sup>6</sup> Institute of Materials for Electronics and Magnetism (IMEM), National Research Council (CNR), Parma, Italy*

#### *Edited by:*

*Alberto Corigliano, Politecnico di Milano, Italy*

#### *Reviewed by:*

*Massimiliano Zingales, Università degli Studi di Palermo, Italy Shangchao Lin, Florida State University, USA*

#### *\*Correspondence:*

*Lucrezia Aversa, Institute of Materials for Electronics and Magnetism (IMEM), CNR, Via Alla Cascata 56/C, Trento 38123, Italy lucrezia.aversa@cnr.it*

#### *Specialty section:*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials*

> *Received: 30 January 2015 Accepted: 27 May 2015 Published: 11 June 2015*

#### *Citation:*

*Aversa L, Taioli S, Nardi MV, Tatti R, Verucchi R and Iannotta S (2015) The interaction of C<sup>60</sup> on Si(111) 7 × 7 studied by supersonic molecular beams: interplay between precursor kinetic energy and substrate temperature in surface activated processes. Front. Mater. 2:46. doi: 10.3389/fmats.2015.00046* Buckminsterfullerene (C60) is a molecule fully formed of carbon that can be used, owing to its electronic and mechanical properties, as "clean" precursor for the growth of carbonbased materials, ranging from π-conjugated systems (graphenes) to synthesized species, e.g., carbides such as silicon carbide (*SiC*). To this goal, C<sup>60</sup> cage rupture is the main physical process that triggers material growth. Cage breaking can be obtained either thermally by heating up the substrate to high temperatures (630°C), after C<sup>60</sup> physisorption, or kinetically by using supersonic molecular beam epitaxy techniques. In this work, aiming at demonstrating the growth of *SiC* thin films by C<sup>60</sup> supersonic beams, we present the experimental investigation of C<sup>60</sup> impacts on Si(111) 7 *×* 7 kept at 500°C for translational kinetic energies (KEs) ranging from 18 to 30 eV. The attained kinetically activated synthesis of *SiC* submonolayer films is probed by *in situ* surface electron spectroscopies (X-ray photoelectron spectroscopy and ultraviolet photoelectron spectroscopy). Furthermore, in these experimental conditions, the C60-Si(111) 7 *×* 7 collision has been studied by computer simulations based on a tight-binding approximation to density-functional theory. Our theoretical and experimental findings point toward a kinetically driven growth of *SiC* on Si, where C<sup>60</sup> precursor KE plays a crucial role, while temperature is relevant only after cage rupture to enhance Si and carbon reactivity. In particular, we observe a counterintuitive effect in which for low KE (below 22 eV), C<sup>60</sup> bounces back without breaking more effectively at high temperature due to energy transfer from excited phonons. At higher KE (22 *< K <* 30 eV), for which cage rupture occurs, temperature enhances reactivity without playing a major role in the cage break. These results are in good agreement with *ab initio* molecular dynamics simulations. Supersonic molecular beam epitaxy is thus a technique able to drive material growth at low-temperature regime.

**Keywords: fullerene, silicon carbide, thin film, surface dynamics, first-principle simulations**

# **Introduction**

The synthesis of carbon-based thin films and nanostructured compounds, such as carbides and graphene, on top of semiconductor or metal surfaces represents a serious challenge to the production of electronic devices and materials coating. Several methods, mainly based on the chemical vapor deposition (CVD) of organic molecules or molecular beam epitaxy (MBE), have been successfully developed to obtain large area carbon-based materials with the common feature of working at very high temperatures. For example, the low energy deposition of fullerene on silicon has been deeply studied with experimental approaches, mainly using standard techniques (Sakamoto et al., 1999; De Seta et al., 2000; Sanvitto et al., 2000; Balooch and Hamza, 1993). MBE, in particular, has shown to be a viable approach to silicon carbide (3C–*SiC*) synthesis at about 800°C, using fullerene (C60) as carbon precursor and silicon as a growth substrate (Sanvitto et al., 2000; De Seta et al., 2000). The required high temperature, however, leads to structural defects at the nano- and micro-scale, mainly due to high lattice mismatch between *SiC* and Si and due to increasing thermal diffusion through the *SiC* film. Moreover, only amorphous Si films can be achieved at 800°C, and postdeposition thermal treatments at higher temperatures are necessary to improve crystallinity.

At room temperature, the most common process at the fullerene–surface interface, be it metallic or semiconductor, is represented by physisorption. Chemisorption is a viable way only on reconstructed silicon surfaces (both Si(111) 7 *×* 7 and Si(100) 2 *×* 1) due to the presence of several dangling bonds on the surface, the most reactive for fullerenes being the 7 *×* 7 reconstruction of Si(111). Furthermore, to induce cage breaking on the top of the surface, one needs to increase the substrate temperature. Finally, subsequent to rupture, the most common outcome on Si(111) 7 *×* 7 surfaces results in the formation of several new Si–C bonds up to the synthesis of silicon carbide (*SiC*) at 800°C.

To overcome the temperature issues, in a previous work (Verucchi et al., 2012), we reported the room temperature (about 20°C) synthesis of nanocrystalline 3C–*SiC* with fullerene on Si(111) 7 *×* 7 surface by the so-called supersonic molecular beam epitaxy (SuMBE) technique. Using this approach, the C<sup>60</sup> translational kinetic energy (KE) can reach values up to 30–35 eV by aerodynamic acceleration due to isoentropic expansion in vacuum. We demonstrated that the cage rupture after impact with the silicon surface is driven by charge excitations to the excited levels of the system and the 3C–*SiC* synthesis is activated by the extra precursor KE content for both high and low (room T) substrate temperatures (Verucchi et al., 2002; Verucchi et al., 2012). At room temperature, where thermal contribution to observed processes is negligible, out-of-thermal-equilibrium conditions achieved by SuMBE induce chemical processes at the silicon surfaces, overcoming activation barriers toward the synthesis of ordered *SiC*. Moreover, by using C<sup>60</sup> supersonic beams at 20 eV, we have shown that a polycrystalline 3C–*SiC* thick layer can be grown on Si(111) 7 *×* 7 at 750°C (Verucchi et al., 2002). This is due to the unquestionable ability of the high fullerene KE to induce cage breaking, synthesis of carbide as well as implement film structural order also in conditions where influence of temperature (and thus thermal equilibrium conditions) is significant.

In order to better identify the role played by temperature on *SiC* growth on Si by SuMBE, in this work, we present the experimental and *ab initio* investigations of supersonic C<sup>60</sup> impacts on Si(111) 7 *×* 7 at temperature of 500°C, where MBE approach is not able to induce *SiC* formation (De Seta et al., 2000). This regime is intermediate between the room temperature experiments, which gave evidence of nanocrystalline *SiC* island formation, and the full synthesis regime at 750–800°C, where all the carbon precursors were transformed into an ordered uniform layer. Synthesis of *SiC* by SuMBE is thus expected also at 500°C, but interplay between fullerene KE and substrate temperature can be better addressed.

Using photoelectron spectroscopy and low-energy electron diffraction (LEED), we explore the chemical, structural, and electronic changes induced by varying the supersonic C<sup>60</sup> beam KE at fixed substrate temperature. Experimental results point toward a *SiC* growth model driven by fullerene KE and where, surprisingly, temperature effects do not affect critically the cage rupture. Furthermore, aiming at modeling the interaction between C<sup>60</sup> and Si surface for several initial KEs and different temperatures, in this work, we also show first-principle simulations based on the density-functional tight-binding (DFTB) method (Matthew et al., 1989; Elstner et al., 1998;Frauenheim et al., 2002) within the Born-Oppenheimer approximation. Classical molecular dynamics simulations of C<sup>60</sup> impacts on Si have been already carried out by Averback (Hu et al., 2000) using a lower accuracy description of the atomic interaction through a Tersoff potential (Tersoff, 1988). However, classical molecular dynamics is not suitable for reactive problems out of thermal equilibrium, where bonds break and form as in our case. At variance, a DFTB-based approach has been recently applied to study high KE impact dynamics of fullerenes on graphite (Galli and Mauri, 1994) and room temperature growth of *SiC* (Taioli et al., 2013), displaying an accuracy comparable to density-functional theory (DFT).

Gaining insights on the chemical–physical mechanisms involved in *SiC* epitaxy by SuMBE and understanding their dependence on externally tunable parameters, such as the C<sup>60</sup> beam KE and the substrate temperature, is at the basis to control and improve *SiC* synthesis even at room temperature and to transfer this technique to produce other materials (Taioli, 2014). Indeed, the study of the C<sup>60</sup> cage breaking dynamics on metallic or semiconductor surfaces can reveal new paths toward the epitaxial growth of materials by our SuMBE approach directly on substrates already employed in the electronic device production, thus suitable for very large scale integration production.

# **Materials and Methods**

Experiments have been performed at the IMEM-CNR laboratory in Trento in a tailored deposition *in situ* characterization facility operating in ultra-high vacuum (UHV). A detailed description of the SuMBE deposition technique can be found in the study by Milani and Iannotta (1999). On the computational side, DFTB calculations have been performed using the DFTB + code (Aradi et al., 2007; Elstner et al., 1998) within the self-consistent charge framework that leads to an improved description of the Coulomb interaction between atomic partial charges.

# **Surface and Structural Characterization**

Chemical and electronic surface properties have been probed by means of X-ray photoelectron spectroscopy (XPS) and ultraviolet photoelectron spectroscopy (UPS). The UHV chamber is equipped with a CLAM2 Electron Hemisperical Analyzer, an MgKα X-ray source, and a helium discharge lamp, thus the excitation photon for XPS is at 1253.6 eV and for UPS the HeI at 21.22 eV. The total resolution allowed by the analyzer is 0.95 eV for XPS and 0.1 eV for UPS. Spectra were acquired at normal electron acceptance geometry. Core-level binding energies (BEs) have been calculated using as a reference the Au surface after sputter cleaning, i.e., the 4f 7/2 Au level at 84.0 eV, while in UPS the spectra are referred to the Fermi level of a gold foil in electric contact with the sample. Core levels have been analyzed by Voigt lineshape deconvolution after background subtraction by a Shirley function (Hüfner, 1995). The typical precision for peak energy positioning is *±*0.05 eV, uncertainty for full width at half maximum (FWHM) is less than 5% and for area evaluation is about 2.5%. Surface structural characterization has been performed by LEED operated at 50 eV.

# **Thin Film Deposition**

Substrates for C<sup>60</sup> deposition were obtained from a Si(111) wafer (resistivity 1.20 *<sup>×</sup>* <sup>10</sup>-4 <sup>Ω</sup>m) cleaned by a modified Shiraki procedure (Verucchi et al., 2002). The silicon oxide film was then removed in vacuum by several annealing cycles, until a complete Si(111) 7 *×* 7 surface reconstruction. The entire procedure was checked by analyzing LEED images, and no traces of oxygen and carbon were detected from XPS analysis. Submonolayer (sub-ML) C<sup>60</sup> films were deposited on clean and reconstructed Si(111) 7 *×* 7 surfaces kept at 500°C, using different precursor KE ranging from 18 to 30 eV. The deposition times were adjusted in order to reach a similar thickness for all the films for ease of comparison, with an average coverage of 0.5 ML (*±* 0.1 ML). The C<sup>60</sup> supersonic source is made of two coaxial quartz capillary tubes, resistively heated by a tantalum foil. The experiments were carried out using a free-jet expansions of H<sup>2</sup> in which the seeded fullerene particles are highly diluted (about 10*−*<sup>3</sup> in concentration) and can reach a KE up to 30 eV, with an average growth rate on silicon of about 0.1 Å/min having the supersonic beam perpendicular to the substrate. Fullerene KEs of about 18, 22, 25, and 30 eV have been employed, as measured from time of flight analysis (Verucchi et al., 2002).

# **Theory**

DFTB is based on a second-order expansion of the density appearing in the full DFT approach to electronic structure calculations. In this approach, the energy of a system of atoms is expressed as a sum of tight-binding-like matrix elements, a Coulomb interaction, and a repulsive pair-potential. The parameters appearing in the formulae expressing these contributions are evaluated using highlevel electronic structure methods and are highly transferable to different physical and chemical environments (Slater–Koster parameters).

The advantage of using DFTB with respect to other *ab initio* methods, such as DFT, is due to the computational cost of this approach to electronic structure that is about two orders of magnitude cheaper than the corresponding full DFT calculation (Garberoglio and Taioli, 2012). As a result of this substantial speed gain, DFTB may be used (i) to investigate much larger systems than those accessible by DFT, (ii) to follow their dynamics for much longer timescales, and (iii) makes it possible to perform several tests by tuning the main parameters affecting the impact at an affordable computational cost.

In our simulations, the Si substrate is represented by eight silicon bilayers cleaved along the (111) direction.

According to the dimer adatom stacking fault model (DAS) (Takayanagi et al., 1985), only the first silicon bilayer participates into the superficial reconstruction, resulting in a 7 *×* 7 reconstructed surface. This reconstructed surface, represented in **Figure 1**, turns out to be optimal for our applications, showing high reactivity due to the presence of 19 dangling bonds on the surface (12 adatoms, six rest atoms, and one atom at the center of the hole in the unitary cell) and to its metallic character (Cepek et al., 1999). The computational supercell used in the impact calculations is hexagonal and, after optimization of both atomic positions and lattice vectors, measures 23.21 Å along the short diagonal and 26.79 Å along the large one. Fullerene was separately optimized and initially placed at 5 Å distance on the top of the silicon substrate and the cell size along the collision direction, perpendicular to the surface plane, was set to 50 Å to avoid spurious interaction among periodic images. The Brillouin zone was sampled at the *⌈*-point only, due to the large number of atoms in the unitary cell (848). Atomic interactions between chemical species were treated by the semirelativistic, self-consistent charge Slater–Koster parameter set "matsci-0-3" (Frenzel et al., 2009). To help and enhance convergence of band structure calculations, we employed a room temperature Fermi smearing of the electronic density. Molecular dynamics simulations were performed in the micro-canonical ensemble (NVE), setting the time step to 1 fs to enforce total energy conservation and each simulation lasted 2 ps.

C<sup>60</sup> molecule.

Finally, the last silicon bilayer was kept fixed during the dynamic evolution to simulate the presence of a silicon bulk.

As a second step of our analysis, we increased the substrate temperature up to 500°C to gain some insight on the role of temperature on *SiC* growth obtained via SuMBE. In particular, our goal was to find out the temperature effect on C<sup>60</sup> cage mechanical stability.

# **Results**

### **Photoemission**

High-resolution C1s core levels for thin films deposited at increasing beam KE's are shown in **Figure 2**. They are all characterized by the presence of a main peak that is not fully symmetric but has a weak shoulder at lower BEs whose intensity increases significantly with beam KE. The lineshape analysis of core levels has been carried out by introducing the main component at 285.2 eV (FWHM of 1.1 eV) and up to two more features at low BE (with fixed FWHM of 1.0 and 0.9 eV), in particular at 284.1 and 283.2 eV. Peak positions do not change as the fullerene KE

**FIGURE 2 | C1s core level spectra (normalized in height) of submonolayer films grown at C<sup>60</sup> kinetic energies of 18, 22, 25, and 30 eV**. Green component corresponds to physisorbed fullerenes, blue component to chemisorbed cages, and red component to broken fullerenes and formation of silicon carbides.

increases. Previous experiments (Verucchi et al., 2012; Taioli et al., 2013) on the impact of supersonic fullerene on silicon at room and higher temperatures have shown that the molecule can undergo different paths, going from physisorption with intact cages to chemisorption, cage rupture and formation of new compounds, such as *SiC*, whose intensity increases with the precursor KE. The main peak at higher BE, P1, reflects the presence of unperturbed and physisorbed fullerenes, while the peaks at 284.1 and 283.2 eV can be identified as related, respectively, to chemisorbed species and to the formation of *SiC* occurring after cage breaking. The higher is the C<sup>60</sup> KE, the more intense is the contribution from "totally reacted" molecules (see **Table 1**), i.e., the fraction of fullerenes leading to new species on the surface. The contribution from physisorbed molecules is always the main one. This is quite surprising, since the 30 eV precursor KE is the same as used in SuMBE room temperature experiments (Verucchi et al., 2012), where physisorbed species are <30 and 50% of the total C1s area at coverages of 0.30 and 0.65 ML.

Chemisorbed fullerenes (peak at 284.1 eV) are present for all films and increase with precursor KE, while *SiC*-related species (peak at 283.2 eV) are present only from the C<sup>60</sup> KE of 22 eV and significantly increases only at 30 eV. The peak at 283.2 eV has a BE between that of 3C–*SiC* (Sakamoto et al., 1999; Mélinon et al., 1998) and non-stoichiometric *SiC* (Liu et al., 1995), thus it can be attributed to *SiC*s showing low structural order. The impossibility to achieve reliable results from the analysis of the Si2p core level (due to the prevailing and masking contribution coming from the bulk silicon on all surface components) makes any comment on the carbide stoichiometry purely speculative.

The valence band photoemission is reported in **Figure 3**, together with reference spectra for comparison from a clean silicon surface, a thick fullerene film and a completely synthesized *SiC* film by means of SuMBE at higher substrate temperatures. Si(111) 7 *×* 7 valence band is characterized by weak and broad structures, the one at about 0.9 eV related to the surface states, and a metallic character evidenced by the presence of signal at about 0 eV. C<sup>60</sup> thick film shows intense and different bands, representing the different π and σ molecular orbitals. In particular, we found at about 2.0 eV the highest occupied molecular orbital (HOMO). At low beam KE, the features related to the clean silicon surface completely disappear, while the HOMO and the other molecular C<sup>60</sup> bands are clearly visible. This is expected since the π and σ bands are reported to be fully developed (Ohno et al., 1991) also at this low coverage (less than 1 ML). There are, however, significant discrepancies between data found in literature, in particular with the work of Sakamoto et al. (1999). HOMO and HOMO-1 molecular levels (at 2.2 and 3.5 eV) are not fully developed and slightly shifted toward higher BEs with respect to the bulk organic film (2.0 and 3.3 eV). Moreover, the σ

**TABLE 1 | Peak area% (***±***3) on total C1s emission**.


bands in the 5–8 eV energy region loose some of their distinctive characteristics, most notably the shoulder at 5.2 eV disappears and the bands at 6.9 and 8.1 eV are no longer clearly resolved. For precursor KEs of 22 and 25 eV, a new peak appears in the UPS spectra between the HOMO and HOMO-1 structures (see arrow in **Figure 3**), probably is present also in film at 18 eV but very weak. A similar peak has been actually observed and identified as related to the strong chemical interaction with the surface (Sakamoto et al., 1999; De Seta et al., 2000). In the highest KE case, the valence band structure changes strongly and only a small residual of the original molecular π structure can be found in the 0–5 eV energy range. At BE of about 6–8 eV, the σbands disappear and a broad structure is formed, reminding of the analogous structure of not-crystalline *SiC* (Aversa et al., 2003) in which the mixed C2p-Si3s orbitals ("sp-like" band in cubic *SiC*) are predominant with respect to the p-like C and Si orbitals at about 3 eV (Mélinon et al., 1998). This suggests that *SiC* synthesis is clearly achieved at 30 eV C<sup>60</sup> KE, even if we are dealing with a partially formed carbide, as demonstrated by the differences with the *SiC* valence band in **Figure 3** (top curve).

There is a fair agreement between XPS and UPS data, suggesting that a C<sup>60</sup> KE of 22 eV highly improves chemisorption processes, while the cage rupture and *SiC* synthesis is at least partially achieved at 30 eV, where formation of defected carbides probably occurs.

### **LEED**

From a structural standpoint, the interaction of C<sup>60</sup> with silicon surface at 500°C leads to quite a few changes, but not to the appearance of new extraspots related to the formation of ordered compounds. **Figure 4** shows the LEED patterns of Si(111) clean surface and of the system after C<sup>60</sup> impact at the lowest and highest KE, respectively, taken for a primary electron energy of 50 eV. The clean surface is characterized by the well-known 7 *×* 7 reconstruction pattern of Si(111) (**Figure 4A**), with bright spots within a high contrast background, typical of a clean and ordered surface. Upon low KE impact of fullerene, the reconstruction pattern disappears and only the 1 *×* 1 hexagonal pattern is visible (**Figure 4B**), with spots having lower brightness and sharpness. At highest energy impacts, the LEED pattern undergoes further modifications, with hexagonal pattern losing spots and becoming more similar to a triangular one. The extra spots (Verucchi et al., 2002) related to the formation of ordered 3C–*SiC* islands do not appear, a result that is in complete agreement with both XPS and UPS results, where formation of a not fully developed and defected carbide is clearly detectable.

#### **Theory**

The first step of our computational analysis is the simulation of the impacts with the silicon surface kept at room temperature. Depending on the initial KE (KC60), we were able to identify three regimes: (i) for KC60 *<* 75 eV surface penetration does not occur with the cage undergoing simple distortion; (ii) between 100 and 150 eV, we found strong C<sup>60</sup> cage distortion and surface penetration; (iii) above 175 eV, the cage starts to break up and Si–C bond formation occurs with small surface damaging.

Thus, by using DFTB to describe interatomic forces in the C60 silicon surface collision, the minimum impact KE to obtain heavy damage of the silicon surface and fragmentation of the C<sup>60</sup> cage is predicted to be well above 200 eV, slightly below Averback's findings. Movies of the full trajectory for 50 and 175 eV initial KE can be found in the Videos S1 and S2 in Supplementary Material.

We point out that there is an important discrepancy between simulations and experiments in the assessment of the precursor KE inducing cage break. In this regard, it is clear that DFTB is missing a comprehensive explanation of the cage break. The reason of this discrepancy was already identified and thoroughly discussed in Taioli et al. (2013). In that analysis, non-adiabatic effects and, specifically, the failure of the Born–Oppenheimer approximation was found responsible for this important discrepancy between simulation and experiments to determine cage rupture KE threshold. In this regard, due to very fast collisional times for the KE under investigations, electronic and nuclear motion do not decouple, as assumed within the framework of the Born–Oppenheimer approximation, as electrons do not relax fast enough to the ground state of the instantaneous positions of the nuclei. Thus, electronic and nuclear dynamics should be treated on an equal quantum mechanical footing.

However, in this work, our goal is primarily to investigate the role that temperature might play in the *SiC* growth by SuMBE. While we underline that a proper non-adiabatic quantum treatment of the collision processes, intertwining electronic and nuclear motion, should be used, in this analysis of temperature effects on growth this approach can be ruled out by two factors, namely physics and computational cost. First, the effect of temperature on the electrons is provided by the Fermi–Dirac distribution. In this statistics, temperature appears in the exponential factor at denominator and plays the role of "smearing" or increasing the electronic population in the excited states of the system. However, a temperature enhancement of 800 K corresponds to about 0.06 eV. The energetic distance between electronic orbitals close to the Fermi level in fullerene is much higher (of the order of 1 eV, corresponding to about 12000 K). Thus, the effect of a small temperature increase, such as that investigated in our work, on the electronic motion can be safely neglected, while nuclear motion is quite sensitive to temperature effects. Second, the computational cost associated with non-adiabatic simulations, at least two orders of magnitude larger than using DFTB, is unaffordable for screening the C60-silicon impact for several KEs and different temperatures. Thus, the use of DFTB is totally justified to investigate temperature effects in the case of high KE fullerene impacts on a silicon surface.

MD simulations (Zhang et al., 1993; Kim and Tománek, 1994) indeed show that a break may occur if the total internal energy stored in the C<sup>60</sup> cage is between 30 and 40 eV. The internal energy is defined as the total KE of the fullerene bouncing back after the impact minus the center of mass KE. A temperature increase may help to destabilize fullerene by exciting substrate phonons, which can transfer energy to C60, eventually resulting in a lower cage breaking KE.

Thus, we performed a temperature dependent investigation of the internal KE after the impact, when fullerene bounces and leaves the silicon surface, in the same range of initial KEs previously examined at a temperature of 500°C. The system was thermostated to this temperature by a Nose–Hoover thermostat for 1 ps to reach thermal equilibrium and then NVE ensemble with a Fermi smearing corresponding to 500°C was adopted during the impact molecular dynamics run.

Surprisingly, the results point toward a not significant effect of temperature on the internal KE for all the impinging KEs.

The results are reported in **Figure 5** and demonstrate that the fullerene internal energy for the impinging KEs under investigation (from 25 to 175 eV) at room temperature conditions ranges from 13 to 28 eV (blue diamonds). For the assessment of the internal KE, we chose a molecular dynamics step corresponding to a maximum of the oscillating total KE immediately after the collision. A substrate temperature increase to 500°C (red diamonds) of course results into a higher final internal KE with respect to the room temperature case (blue diamonds) due to enhanced thermal fluctuations. However, this increase is not such to cause the rupture in the 25–175 eV energy range as much of the energy is spent by fullerene to increase its translational KE for leaving the surface, as shown by the evolution of the center of mass KE of C<sup>60</sup> after impact (see blue and red circles in **Figure 5**).

We point out that the interaction potential energy of C<sup>60</sup> as a function of the distance from the silicon surface is not a quantity at our disposal (just the total interaction potential energy is, as shown below). This analysis could be in principle performed using classical molecular mechanics, where all the interactions terms are analytically known and are singly computable within the reach of this theory. However, in the conditions at which C60–Si collision occurs, we are abundantly out of thermodynamic equilibrium with bonds breaking and forming. Thus, we are forced to use a "first-principles" technique, explicitly including electronic motion in a self-consistent way, to have a reliable description of the high KE impact. In DFTB, single terms of the potential energies are not a direct output, particularly as a function of the distance between C<sup>60</sup> and the silicon surface. In principle, one could perform a very expensive analysis to obtain this quantity, consisting in repeating all calculations at different KEs and temperatures as follows. Indeed, one could calculate separately the interaction potential energy of the isolated C<sup>60</sup> and silicon slab at each position, keeping fixed the atomic coordinates as by simulations of C<sup>60</sup> approaching

the surface. Finally, one could subtract the sum of these two separate contributions from the total interaction potential energy of the real interacting system at all distances, thus obtaining the C60–silicon slab interaction potential energy as a function of the distance. Such a calculation is of course unfeasible due to the computational scaling of our approach.

The plot in **Figure 6** represents the total interaction potential energy, a quantity computable by DFTB, for three different KEs and two different temperatures as a function of computing time, which is somehow proportional to the distance of C<sup>60</sup> from the silicon surface (see **Figure 1**). In this plot, the potential energy curves for 800 K are shifted up in energy by 600 eV for the sake of visibility. From **Figure 6** one can clearly see that, at the point where C<sup>60</sup> starts experiencing the interaction with the silicon surface (below 100 fs), the total interaction potential energy changes abruptly only in the case of high KE (350 eV, orange and cyan lines for the cases of ambient temperature and 800 K, respectively). For this KE, C<sup>60</sup> is completely fragmented. Furthermore, in this high KE regime, the difference between the total interaction potential energy at ambient temperature and 800 K is clearly visible at all times. At variance, for low KE (25 eV, green and violet curve for the cases of ambient temperature and 800 K, respectively), where C<sup>60</sup> is just physisorbed on the silicon surface, the total interaction potential energy shows little change with time and is almost independently of temperature. Finally, at intermediate regimes (175 eV, red and black curve for the cases of ambient temperature and 800 K, respectively), where C<sup>60</sup> is partially fragmented, the change of the total interaction potential energy again seems to be due more to kinetic effects than to temperature, however showing larger fluctuations with increasing temperatures with respect to the ambient temperature impact.

To further enrich our discussion about temperature effects on *SiC* growth by SuMBE, we performed the calculation of the mean squared displacement (MSD) as a function of impinging C<sup>60</sup> KE and time for different temperatures. This quantity provides of course a measure of the dispersion or diffusion of carbon atoms of the C<sup>60</sup> molecule. These results are sketched in **Figures 7A,B**. In **Figure 7A**, we report, in particular, the behavior of the MSD of C<sup>60</sup> carbon atoms at ambient (black curve) and 800 K (red curve). From this figure, one safely concludes that temperature effects are relevant only at high KE regimes, when the cage is fragmented (350 eV) or partially broken (175 and 260 eV).

In **Figure 7B**, MSD is sketched as a function of time (or, which is the same, of C<sup>60</sup> approaching the silicon surface). Again, MSD is almost constant or slightly increasing at low KE regimes (25 and 175 eV) while is divergent at 350 eV, meaning that C<sup>60</sup> is completely fragmented. Furthermore, MSD is larger at higher temperature with a sharp increase at long times. This analysis clearly points toward a cage break model that is almost independent on temperature with kinetic effects playing a major role. Temperature effects become important only after cage rupture.

# **Discussion**

Experimental results of C<sup>60</sup> supersonic impact on Si(111) 7 *×* 7 at 500°C show a picture completely different from what was observed using the same deposition technique at higher and lower

substrate temperatures, evidencing an unexpected mechanism of interaction between the precursors and the surface. From previous experiments (Verucchi et al., 2002, 2012), as well as theoretical calculations (Taioli et al., 2013), one finds a clear indication of a kinetically driven growth mechanism leading to sub- or thick *SiC* layer formation already at room temperature, with temperature factor not playing a major role. Raising the silicon substrate temperature to 500°C should generally result in a higher silicon–carbon reactivity, with increasing carbide percentage over physisorbed fullerene. Core level and valence band analysis (**Figure 2**) point toward a completely different situation with respect to our expectations. Indeed, the great majority of species present on the surface after impact are unperturbed C<sup>60</sup> up to high KE, whereby a reactivity increase is found in terms of percentage of chemisorbed fullerene and carbide formation. Valence band spectra suggest the same behavior, with molecular bands having characteristics more similar to physisorbed fullerene with C<sup>60</sup> cage breaking and carbide synthesis clearly evident only for the 30 eV KE case. Nevertheless, we achieve formation of defected *SiC* and not ordered cubic 3C–*SiC*, as observed in previous experiments at 750–800°C (Verucchi et al., 2002) and room temperature (Verucchi et al., 2012).

At 500°C, C<sup>60</sup> KE is crucial, with a threshold for *SiC* synthesis and cage breaking identified at about 30 eV. However, it is clear that chemical reactivity is different from both 750°C and room temperature cases. This may represent a surprising result, assuming that the higher the substrate temperature the easier should be the cage breaking. Indeed our calculations clearly indicate that the internal energy for a fixed value of the C<sup>60</sup> initial KE is slightly higher if temperature is raised up to 500°C (see **Figure 5**). This trend is confirmed for all KEs investigated. However, our computational analysis shows that a major part of both the increasing KE and temperature is spent to eject more effectively the fullerene at the expense of an enhanced reactivity. This can be clearly seen in the low and intermediate (below 150 eV) KE regimes, where increasing the temperature corresponds to enhance proportionally the center of mass KE of the C<sup>60</sup> bouncing back from the surface instead of the internal energy of the system, responsible for the cage break. A similar effect was observed in the collision of SF<sup>6</sup> from GaSe surface where an uptake of energy was identified (Boschetti et al., 1992). Additionally, at high KEs, when fullerene gets trapped within the silicon substrate and thus do not bounce back, temperature effects start to be more effective in transferring energy from silicon to fullerene vibrational modes, as reported in **Figure 5**, in perfect agreement with our experimental results. However, this is not enough to justify cage opening. Thus, both calculations and experimental results clearly point toward a kinetically driven fullerene cage breaking and *SiC* synthesis, with a KE threshold not critically affected by temperature below the well-known value of 750°C (Verucchi et al., 2002).

In conclusion, we have analyzed in details the impacts of energetic supersonic C<sup>60</sup> on the reconstructed Si(111) 7 *×* 7 surface at 500°C by means of experiments and *ab initio* simulations. Surface chemical and structural characterization carried out at increasing fullerene KE (18 eV to 30 eV) has shown the presence of physisorbed, chemisorbed, and reacted molecules, the latter leading by the formation of *SiC*, probably defected cubic carbide or non-stoichiometric species. Compared to previous experiments with substrate kept at room temperature, we do not find an enhancement in reactivity with increasing temperature. DFTB calculations reveal that the excess of initial KE is spent more effectively to increase the translational KE (center of mass KE) and substrate temperature helps this process via phonon transfer, at least in the low and intermediate KE regime. Only for higher KE, when cage rupture occurs, the reactivity starts to be relevant, as shown in our MSD analysis, demonstrating that this process is kinetically driven. While our discussion has been limited to the case of *SiC* growth, we believe that SuMBE may represent a more general approach to grow materials in the low-temperature regime even on metals, notably nickel and copper, in which carbide formation may compete energetically with π-conjugate structures, such as graphene (Batzill, 2012; Lu et al., 2011).

# **Acknowledgments**

The research leading to these results has received funding from the European Union Seventh Framework Programme under grant agreement no. 604391 Graphene Flagship and from project Super-Car (CMM-FBK, Fondazione Bruno Kessler). ST acknowledges support from the Istituto Nazionale di Fisica Nucleare through the "Supercalcolo" agreement with Fondazione Bruno Kessler, from the European Science Foundation under the INTELBIOMAT Exchange Grant "Interdisciplinary Approaches to Functional Electronic and Biological Materials" and the Bruno Kessler Foundation for providing economical support through the "research mobility scheme" under which this work has been accomplished. Finally, ST gratefully acknowledges the Institute of Advanced Studies in Bologna for the support given under his ISA research

# **References**


fellowship, and the high-performance computing service Archer in UK where impact calculations were performed.

# **Supplementary Material**

The Supplementary Material for this article can be found online at http://journal.frontiersin.org/article/10.3389/fmats.2015.00046/ abstract


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Aversa, Taioli, Nardi, Tatti, Verucchi and Iannotta. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Can excited electronic states of macromolecules with extended pi-systems be reliably predicted? A case study on P3HT

# **Marcus Böckmann and Nikos L. Doltsinis\***

Institute for Solid State Theory and Center for Multiscale Theory & Computation, University of Münster, Münster, Germany \*Correspondence: nikos.doltsinis@uni-muenster.de

**Edited by:** Simone Taioli, Bruno Kessler Foundation, Italy

**Reviewed by:** Elise Dumont, Ecole Normale Supérieure de Lyon, France

**Keywords: polymers, excited states,TDDFT, multiscale modeling, exchange-correlation functionals**

## **1. INTRODUCTION**

Predicting molecular properties in excited electronic states has always posed a particular challenge compared to the electronic ground state. Since the breakthrough of time-dependent density functional response theory (TDDFT) in the late 1990s (Casida, 1995; Marques et al., 2006), electronic excitation energies can often be calculated to a fraction of an electron volt at roughly the same computational cost as the ground state. Due to the wide availability of efficient implementations of TDDFT and its quasi "black box" nature, the method became extremely popular, not only among theoreticians but also among experimentalists. Typically, the calculation of an absorption spectrum merely involves a geometry optimization in the ground state followed by a TDDFT calculation. Even the few well-known failures of standard TDDFT, such as charge-transfer and Rydberg excitations, can be easily cured by introducing a long-range correction into the exchange-correlation functional (Iikura et al., 2001; Dreuw et al., 2003, 2010; Dreuw and Head-Gordon, 2004; Yanai et al., 2004; Baer and Neuhauser, 2005; Zhang et al., 2012). However, for very large molecules, it is very difficult to assess the accuracy of TDDFT, because potentially more accurate *ab initio* methods are no longer applicable due to their prohibitive computational cost. The situation becomes more complicated when vibrational effects or even conformational dynamics play an important role, e.g., in spectral broadening. If only a single potential well is involved, one may sample molecular structures from a Wigner distribution (Beutier et al., 2014). If, however, multiple potential minima contribute to the overall spectrum, molecular configurations need to be sampled by molecular dynamics or Monte Carlo simulation.

Moreover, the nature of the molecular environment can crucially impact on the conformational distribution. Hence single molecule calculations are generally insufficient to reproduce experimental spectra. In this article, we will focus on the example of Poly(3-hexylthiophene) (P3HT) – a polymer that is frequently used in organic solar cells and semiconductors and has been studied extensively both experimentally and theoretically (Pope and Swenberg, 1999; Perepichka and Perepichka, 2009; Bhatt et al., 2014; Poelking et al., 2014) (see **Figure 1A**).

# **2. DISCUSSION**

**Figure 1B** shows the absorption spectra of P3HT in an ordered crystal and in solution (Rahimi et al., 2012, 2014). Experimentally, the main absorption maximum in solution is at 450 nm, while the crystal spectrum exhibits a large red-shift to about 700 nm. The spectrum of an amorphous film is situated in between the two extremes (Rahimi et al., 2012, 2014). To understand the origin of the spectral shifts, it is necessary to disentangle the many-body effect on the electronic structure from the differences in the conformational distributions. From a computational point of view, calculating the absorption spectrum of a condensed phase system becomes a rather hopeless task if the intermolecular interactions are too strong to approximate the overall spectrum as the sum of the individual molecular spectra. But even if the molecules can be treated as isolated entities as far as their

electronic structure is concerned, precise knowledge of their conformations is still a prerequisite for the accurate prediction of absorption spectra in condensed phases. This is where it becomes a multiscale problem. The structure of the – disordered – condensed phase system has to be generated by a molecular dynamics (MD) or Monte Carlo (MC) simulation. In the case of polymeric materials, *ab initio* sampling techniques are out of the question because their computational cost is too high for the large time and length scales involved in these systems. Somewhat paradoxically, however, some *ab initio* calculations are generally still necessary to reparametrize the atomistic force field for the specific system under study, since most "out-of-thebox" force fields have been developed for biological molecules. Poor force field parameters can lead to dramatic failures as the calculated absorption spectrum of P3HT in CHCl<sup>3</sup> solution shown in the bottom panel of **Figure 1B** illustrates. The main peak of the theoretical spectrum calculated with TDDFT/LC-BLYP/6-31G\* for geometries sampled from an MD simulation with force field I (see Computational Details) is situated at about 2.1 eV, a staggering 0.7 eV away from the experimental maximum at ca. 2.8 eV. This can mainly be attributed to the fact that the torsional potential between neighboring five-membered rings is too stiff in the standard GROMOS 45a3 force field. We can exclude the possibility that the deviation is due to the TDDFT calculation, since the absorption spectrum of the ordered crystal with planar molecules can be reproduced rather accurately using the LC-BLYP functional and the 6-31+G\* basis set. For the comparison of theoretical and

#### **FIGURE 1 | Continued**

**(A)** Molecular structure of P3HT in the crystal (left) and a snapshot from the MD simulation in CHCl<sup>3</sup> solution with force field II (right). For clarity, the hexyl groups have been replaced by methyl groups. **(B)** Top panel: theoretical LC-BLYP (red) and experimental (black) absorption spectra of P3HT in a perfect crystal. Bottom panel: theoretical absorption spectra of P3HT in CHCl<sup>3</sup> solution obtained with force field I and LC-BLYP (green dotted line), force field II and LC-BLYP (blue dashed line), and force field II and PBE0 (black solid line)

compared to the experimental spectrum (black solid line) (Rahimi et al., 2014). **(C)** Dependence of the position of the main absorption peak of planar P3HT oligomers on the repetition length n. The black circles have been obtained from TDDFT/LC-BLYP/6-31+G\* calculations, while the red line is a Kuhn fit according to equation (1). **(D)** Dependence of theoretical absorption spectra of a P3HT 32mer in CHCl<sup>3</sup> solution on the S–C–C–S torsional angle θ. Each spectrum is for a single geometry with fixed θ: θ = 0°(black), θ = 40°(red), θ = 80°(blue). Lorentz broadening with a FWHM of 50 nm has been applied.

experimental spectra, we need to bear in mind that the typical repetition length *n* in commercial P3HT materials is of the order of 100, which is too large for TDDFT calculations (*n* = 32 was used here). A practical solution that has been applied in the literature is to carry out a series of calculations on small oligomers and extrapolate to long chains (Gierschner et al., 2007). In particular, the Kuhn extrapolation formula for excitation energies (Kuhn, 1948),

$$E = E\_0 \sqrt{1 + 2\frac{k'}{k\_0} \cos\left(\frac{\pi}{N+1}\right)} \quad (1)$$

has been used successfully for a number of different polymer species. In equation (1), *N* is the number of double bonds along the backbone connecting the terminal carbon atoms (*N* = 2*n* in the case of P3HT), *E*<sup>0</sup> is their individual vibrational energy, *k*<sup>0</sup> their force constant, and *k* 0 their coupling constant. A Kuhn fit to the vertical excitation energies of P3HT oligomers calculated with TDDFT/LC-BLYP/6-31+G\* is shown in **Figure 1C**. It demonstrates that the excitation energy of the 32mer used for calculating the spectra in **Figure 1B** is close to the asymptotic limit. It is worth mentioning that typically only a modest basis set can be used for large molecules, introducing a further source of error. In the present case, increase the basis to cc-pVTZ leads to a decrease in excitation energy of the 5mer by 0.02 eV.

Now that we have established that TDDFT in combination with a rangeseparated functional is able to reproduce accurately the spectrum of a planar conformation, let us take a closer look at disordered phases with non-planar polymers. The role of the torsion becomes clear when the MD simulation is repeated with a more accurate – much softer – torsional potential derived from *ab initio* calculations (see Computational Details). As can be seen from the bottom panel of **Figure 1B**, the main absorption maximum obtained with the LC-BLYP functional shifts to a value of about 3.4 eV (blue dashed line) as the torsional angles increase on average. The reason for the shift is thus not the interaction between the solute and the solvent, but rather the temperature dependent conformational change of the P3HT solute. By twisting the monomeric thiophene units of the polymeric chain against each other, the extended π electron system is interrupted and broken down into smaller units as illustrated by the snapshot shown on the rhs of **Figure 1A**. **Figure 1D** illustrates how the absorption spectrum shifts as a function of the torsional angle θ when all thiophene units are twisted by the same angle.

To some extent, the absorption spectrum of P3HT in solution can thus be approximated as that of a distribution of shorter chains. According to equation (1), an excitation energy of 3.4 eV corresponds to a repetition length of *n* = *3* (see also **Figure 1C**). This in turn raises the question as to whether a long-range corrected exchange-correlation functional such as LC-BLYP is still appropriate for a distorted chain. It has been shown that rangeseparated functionals significantly overestimate excitation energies of local transitions (Nguyen et al., 2011). The present observation that the LC-BLYP spectrum is blue-shifted by roughly 0.6 eV with respect to experiment seems to confirm that. A hybrid functional would appear to be a better choice here. Indeed, the PBE0 spectrum seen in **Figure 1B** is in excellent agreement with experiment (B3LYP red-shifts the spectrum by 0.1 eV). This is plausible as the majority of conformations accessed during the dynamics are highly twisted and the minority of near-planar conformations does not significantly influence the average absorption spectrum. For those minority conformations, the PBE0 functional, nevertheless, dramatically underestimates the excitations energies.

The most difficult scenario is the intermediate case between the planar conformation of the crystal and the twisted geometries of the solution. This is true, for example, for a thin, amorphous film of P3HT, whose main absorption peak was measured to lie at about 2.4 eV (Rahimi et al., 2014). Which functional should one choose to describe this situation in which extended π systems are of equal importance as twisted geometries, or many conformations are twisted to a lesser degree than in solution? One might be able to define geometric criteria to assign individual conformations to one or the other extreme, but this would certainly require extensive prior investigation. Furthermore, such a procedure would still not provide a solution for how to deal with intermediate geometries.

## **3. CONCLUSION**

In conclusion, we have shown using the illustrative example of P3HT that the calculation of excited electronic states of polymeric materials is inherently a multiscale problem. Since the electronic properties of each molecule sensitively depend on its conformation, it is imperative to first obtain reliable information about the atomistic structure of the material. For disordered, non-crystalline materials this implies large-scale molecular dynamics or Monte Carlo simulations, often requiring a coarse-grained representation. As an aside, we pointed out that the construction of atomistic and coarse-grained force fields is itself a multiscale task, as the parametrization of each level relies on a suitable set of reference data on the level below.

Once the molecular structure has been established with sufficient accuracy, the electronic structure problem can be addressed. The method with the best accuracy to cost ratio to calculate excited states is certainly TDDFT. However, which functional is the most appropriate crucially depends on the molecular conformation. Range-separated functionals need

to be used when the π delocalization is uninterrupted, while hybrid functionals might be the best choice for highly twisted geometries. What is needed is a functional that is able to yield accurate results over a wide range of distances.

#### **4. COMPUTATIONAL DETAILS**

All absorption spectra were calculated with the Gaussian 09 Rev. D.01 quantum chemistry package (Frisch et al., 2009) using the TDDFT method with the LC-BLYP (Becke, 1988; Lee et al., 1988; Iikura et al., 2001) or PBE0 (Adamo and Barone, 1999) exchange-correlation functional. The default basis set was 6- 31G\*. The crystal spectrum was modeled by single, planar, and regioregular oligomers cut out from the crystal xray structure Dudenko et al. (2012). The lowest energy absorption peak was recalculated with the 6-31+G\* basis set for the 5mer, 10mer, and 20mer; a 6-31+G\* estimate for the 40mer lowest excited state was obtained by subtracting the shift determined for the smaller oligomers (0.04 eV) from the corresponding 6-31G\* value.

For the liquid phase simulation, a P3HT 32mer was embedded in 32653 CHCl<sup>3</sup> molecules in a cubic periodic unit cell of length 15.952 nm. A molecular dynamics run was performed with the Gromacs 4.5.3 program (Van der Spoel et al., 2005; Gromacs, 2015) at 400 K using the Gromos 45a3 force field, a timestep of 1 fs, and a Berendsen thermostat. The point charges were modified according to a ChelpG electrostatic potential fit (Breneman and Wiberg, 1990) for the optimized 4mer with the B3LYP functional (Becke, 1993) and the 6-31G\* basis set. We refer to this force field as force field I. As it turned out, the torsional potential between adjacent thiophene rings is too stiff in force field I. Therefore, a second force field was designed (force field II) whose S–C–C–S dihedral term was reparametrized according to B3LYP/6-31G\* reference calculations [more details will be published separately (Böckmann et al., forthcoming)]. The solution phase absorption spectrum was calculated by averaging over ten P3HT geometries sampled at random from the 2 ns MD trajectory. In each TDDFT calculation, the lowest 30 excited states were computed. Only the P3HT solute was treated

explicitly, while the solvent was taken into account through a PCM model with UFF atomic radii. In the present case, the solvatochromic shift was found to be small; the peak is red-shifted by just 0.03 eV using the PCM compared to calculations in vacuum. Inclusion of explicit solvent molecules might introduce another shift. It should also be noted that sampling a larger number of configurations from the MD run might again slightly change the position and shape of the spectrum, but the high computational effort generally prevents one from doing so for macromolecules.

#### **AUTHOR CONTRIBUTIONS**

MB carried out the calculations and contributed to the manuscript. ND conceived the work and wrote the manuscript.

#### **ACKNOWLEDGMENTS**

We are grateful to Günter Reiter for providing the experimental absorption spectra. *Funding:* This work has been funded by DFG within TRR 61.

#### **REFERENCES**


for long-range charge-transfer excited states: the zincbacteriochlorin-bacteriochlorin and bacteriochlorophyll-spheroidene complexes. *J. Am. Chem. Soc.* 126, 4007–4016. doi:10.1021/ja039556n


processes for generating large single crystals of poly(3-hexylthiophene). *Angew. Chem. Int. Ed.* 51, 11131–11135. doi:10.1002/anie.201205653


functional approach for large-scale simulations. *J. Phys. Condens. Matter* 24, 205801. doi:10.1088/ 0953-8984/24/20/205801

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 02 February 2015; paper pending published: 20 February 2015; accepted: 13 March 2015; published online: 07 April 2015.*

*Citation: Böckmann M and Doltsinis NL (2015) Can excited electronic states of macromolecules with extended* *pi-systems be reliably predicted? A case study on P3HT. Front. Mater. 2:25. doi: 10.3389/fmats.2015.00025 This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials.*

*Copyright © 2015 Böckmann and Doltsinis. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **Surfaces and interfaces: combining electronic structure and electron transport models for describing electron spectra**

*László Kövér\**

*Institute for Nuclear Research, Hungarian Academy of Sciences, Debrecen, Hungary*

**Keywords: surfaces, interfaces, electronic structure, electron transport, electron spectra, models/simulations**

# **Scales of Physical Phenomena Reflected in Electron Spectra Induced from Atoms, Molecules, Nanostructures, and Solids by Photons or Charged Particles**

#### *Edited by:*

*Simone Taioli, Bruno Kessler Foundation, Italy; Charles University, Czech Republic*

#### *Reviewed by:*

*Lucrezia Aversa, National Research Council, Italy Lucia Calliari, Bruno Kessler Foundation, Italy*

> *\*Correspondence: László Kövér lkover@atomki.mta.hu*

#### *Specialty section:*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials*

> *Received: 27 January 2015 Accepted: 02 April 2015 Published: 22 April 2015*

#### *Citation:*

*Kövér L (2015) Surfaces and interfaces: combining electronic structure and electron transport models for describing electron spectra. Front. Mater. 2:35. doi: 10.3389/fmats.2015.00035* There are many different types of physical and chemical phenomena reflected in electron spectra induced from atoms, molecules, nanostructures, and solids by photons or charged particles. **Table 1** compares the spatial, time, and typical excitation energy scales of selected phenomena. The data (Krausz and Ivanov, 2009; Surdutovich and Solov'yov, 2014) in **Table 1** demonstrate that the phenomena influencing electron spectra take place on very wide scales, sometimes encompassing differences of many orders of magnitude. This suggests the need for multiscale approaches when extracting electronic structure or chemical analytical information from electron spectra.

# **Electron Spectra**

Electron spectra induced from solids contain information on (i) the electronic structure (local at atomic, molecular, cluster, and surface/interface level – or delocalized in the bulk of the solid), different single or multiple excitations, and the chemical state resolved chemical composition; (ii) the (ordered, disordered, or amorphous) physical structure of the solid, as well as the lateral or indepth concentration distribution of the constituent chemical species, via effects of electron transport phenomena (including electron diffraction). While core and valence band photoelectron as well as Auger electron spectra reflect the structure of occupied electronic states, shake-up satellites, due to excitations localized at atomic or molecular level, reflect the structure of unoccupied electronic levels. In the case of shake-off excitations, the energy is continuously shared between the two emitted photoelectrons; therefore, the contribution of these phenomena cannot be separated in the photoelectron spectra from the continuous energy distribution background of electrons scattered inelastically in the solid. Shake-off excitations, however, result in energy loss satellite peaks in photoinduced Auger spectra.

Core photoelectron peaks carry localized signatures characteristic of the atomic environment of the emitting atom, while valence band photoelectron spectra reflect properties of delocalized electronic states. On the other hand, information on the site projected density of valence electron states can be obtained from the analysis of core–core valence (CCV) or CVV Auger spectra. The shape of the core photoelectron lines can, however, be influenced by delocalized excitations (e.g., by excitation of electron-hole pairs in the conduction bands of metals or by "intrinsic" type plasmon excitations due to the appearance of the core hole).



From the point of view of energy loss of photoexcited electrons emitted from solids, "intrinsic" type excitations, connected with the suddenly created core holes localized on the core ionized atoms, usually lead to energy loss satellite structures in the electron spectra. Moreover, all electrons emitted from solids are affected by excitations occurring at surface crossings ("surface excitations," localized in two dimensions) and resulting in energy losses, leading to surface plasmon satellite peaks. Finally, collective excitations of weakly bound, free, or near free electrons in the solid by photoinduced electrons during their transport within the solid ("extrinsic" excitations) lead to delocalized volume excitations manifested, e.g., in the presence of bulk plasmon satellite peaks in the spectra. It should be noted that interferences can take place between surface and bulk, extrinsic, and intrinsic collective excitation phenomena.

# **Modeling/Interpreting Electron Spectra Using Classical and Quantum Frame Models**

For describing or simulating the electronic structure and effects of electron transport on electron spectra emitted from solids, combinations of models based on different approximations are used. Here, only some examples are mentioned for illustration. "One step" quantum mechanical (QM) models include full QM calculations accounting for all phenomena, while the "two- or morestep" models assume that different phenomena can be separated. A recent report (Minár et al., 2011) on "one step" models describing angle resolved photoelectron spectra provides examples of applications to metals and alloys.

An example of "one step" QM models is the quantum Landau theory (QLT) developed by Fujikawa et al. (2008). This model explains all plasmon energy loss features occurring due to core level photoemission from solids and fully accounts for intrinsic and extrinsic plasmon excitations and interferences between these processes, for multiple plasmon energy loss features in the spectra, as well as for effects of multiple elastic electron scatterings (including photoelectron diffraction) preceding and following the energy loss of photoelectrons. For high kinetic energy photoelectrons, the QLT performs well, for ~100 eV kinetic energy photoelectrons the model was successfully applied in the case of Al 2s, Na 2s, and Li 1s spectra photoexcited from single crystals (Kazama et al., 2014; Niki et al., 2014).

While "one step" QM models require demanding calculations even in the case of simple ordered systems and it is difficult to include multiple electron scattering, "two- or more-step" models usually assume the separability of the intrinsic and extrinsic processes and use semiclassical or classical approximations for describing the latter phenomena, neglecting interference effects between intrinsic and extrinsic or surface and bulk processes in most cases. This practically means the separation of electron transport and photoionization/excitation processes, focusing on near surface electron transport.

For describing effects of electron transport in solids, different approaches can be followed. A possibility is to deconvolute contributions of electron transport from measured electron spectra in order to obtain spectra carrying information on the (local) electronic structure. An example is the Monte Carlo (MC) simulation-based deconvolution technique for removing contributions due to different kinds (intrinsic, extrinsic, surface, bulk) of multiple plasmon excitations from X-ray photoelectron spectra. A method based on the partial intensity analysis (PIA) approximation (Werner, 2001) using a linearized Boltzmann-type kinetic equation accounts for multiple elastic electron scattering as well, while neglecting core hole and interference effects. The contribution from intrinsic plasmon excitations is assumed to be similar to that from the bulk extrinsic excitations and removed by applying an "adaptive deconvolution." This approach seems to work especially well, e.g., in the case of high energy photoelectron and Auger spectra excited from 3d metals. In this particular case, the asymmetry of the core photoelectron lineshapes (due to the sudden creation of electron-hole pairs in the conduction band upon core level photoionization) can be described using an independent model (Hughes and Scarfe, 1996) based on the joint density of electronic states (that can be obtained from density functional calculations) around the atom with the core hole (Novák et al., in preparation). The MC method based on the parameters of the PIA model, combined with an expert system containing extended databases on physical data characterizing photoionization and electron transport in solids, is applicable for simulation of X-ray photoelectron and Auger spectra excited from solid structures (Werner et al., 2014). It should be noted, however, that accurate simulation of photoelectron spectra is a very demanding task due to the very complex assembly of physical phenomena to be accounted for.

A different approach is first to model/simulate "intrinsic" spectra on the basis of the available knowledge on the atomic and electronic structure, then to simulate the contribution of electron transport to "intrinsic" spectra with the aim to directly compare the resulted synthetic spectra to experimental electron spectra. An example of such an approach is using a mixed QM and MC method for describing Auger spectra photoexcited from SiO<sup>2</sup> nanoclusters of different sizes (Taioli et al., 2009a,b). Within this approach first the "intrinsic" Auger spectra were obtained from cluster *ab initio* QM calculations then the effects of inelastic energy losses of Auger electrons within the nanoclusters on the Auger lineshapes were simulated using the MC method. Synthesized and measured spectra compare fairly well (Taioli et al., 2009a,b). The model assumes fixed nuclei (so phonons are not considered). However, shake-up and shake-off transitions, though neglected in Taioli et al. (2009b), are foreseen by the model.

An example of models describing excitations due to the core hole and to photoelectrons during their transport in solids is the semiclassical dielectric response model (DRM) (Tougaard and Yubero, 2004), which includes the polarization accompanying photoelectrons on their way toward the surface within the solid, near the surface inside and outside of the solid as well as the influence of the stationary hole during the whole photoemission process. The primary excitation spectrum and the dielectric function are input physical quantities to the DRM that, in addition to bulk excitations, describes core hole effects (intrinsic excitations), surface excitations as well as interference effects (between intrinsic and extrinsic or surface and bulk excitations) (Tougaard and Yubero, 2004). On the other hand, if one has detailed knowledge (e.g., from independent experiments or models) on the parametrized primary excitation spectrum, the actual parameter values (energy position, energy widths, and intensity of the component peaks/structures) can be derived from the measured photoelectron spectrum with the help of the DRM and the measured photoelectron spectrum. First, the effective cross section for inelastic electron scattering can be derived using the DRM, then this cross section is convoluted with the parametrized primary excitation spectrum varying the parameters until achieving a good

# **References**


agreement with the experimental XPS spectrum (Pauly et al., 2014). A recent study applying the DRM for determination of the Cu 2p primary excitation spectrum in CuO demonstrates (Pauly et al., 2014) the use of this approach yielding results agreeing with those predicted by the charge transfer multiplet model calculations. The DRM, however, while accounting for multiple inelastic electron scattering, neglects the effects of multiple elastic electron scattering.

# **Challenges in Characterization of Advanced Materials: Need for Multiscale Materials Analysis and Modeling**

Development of advanced materials, e.g., layered structure doped nanocomposites requests 3D monitoring of their chemical, physical, and electronic structure on multiple scales as a function of their mechanical, optical, and electric properties. Such a need for multiscale materials analysis, however, poses great challenges for the available methods of characterization and usually multimethod approaches are needed.

The increasingly popular, non-destructive HArd X-ray PhotoElectron Spectroscopy (HAXPES) with its flexible applicability for chemical state resolved studies of surface (at grazing photon incidence) and deeply (several 10 nm) buried interface structures, is one of the promising candidate methods for multiscale analysis (Kövér, 2010). The information depth of HAXPES using grazing photon incidence (at total X-ray reflection conditions) is set by the depth of penetration of X-rays and it is limited to the uppermost atomic layers. On the contrary, tuning the photon energy in the case of near normal photon incidence makes the deeply buried interfaces accessible. Angular resolved HAXPES and HAX-PES using X-ray standing waves for excitation combined with advanced modeling allow accurate quantitative chemical analysis of concentration depth profiles of deep regions.

# **Summary**

Electron spectra of solids reflect phenomena on different spatial and time scales and events at various levels of localization. Combining models/simulations for local electronic structure and electron transport is necessary for accurate interpretation of experimental electron spectra. Analysis of advanced materials requires novel methods and models for multiscale approaches.


Taioli, S., Simonucci, S., and Dapor, M. (2009a). SURPRISES: when ab initio meets statistics in extended systems. *Comput. Sci. Disc.* 2, 015002. doi:10.1088/ 1749-4699/2/1/015002

Taioli, S., Simonucci, S., Calliari, L., Filippi, M., and Dapor, M. (2009b). Mixed ab initio quantum mechanical and Monte Carlo calculations of secondary emission from SiO2 nanoclusters. *Phys. Rev. B* 79, 085432. doi:10.1103/PhysRevB.79. 085432


*Standard Reference Database 100*. Available at: http://www.nist.gov/srd/ nist100.cfm

**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Kövér. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# A bird's eye view on the concept of multichannel scattering with applications to materials science, condensed matter, and nuclear astrophysics

# **Simone Taioli 1,2\***

<sup>1</sup> European Centre for Theoretical Studies in Nuclear Physics and Related Areas (ECT\*), Trento Institute for Fundamental Physics and Applications (INFN-TIFPA), Bruno Kessler Foundation, Trento, Italy

<sup>2</sup> Faculty of Mathematics and Physics, Charles University, Prague, Czech Republic

\*Correspondence: taioli@ectstar.eu

#### **Edited by:**

Alberto Corigliano, Politecnico di Milano, Italy

#### **Reviewed by:**

Giulio Giuseppe Giusteri, Università Cattolica del Sacro Cuore, Italy

**Keywords: scattering theory, Fano's discrete-continuum interaction, auger effect, ultra-cold fermi gases, stellar nucleosynthesis**

The last few years have witnessed remarkable advances in two branches of condensed matter physics. The first deals with the experimental investigation by electron spectroscopy techniques of the laws governing the behavior of strongly correlated systems. This flourishing experimental area has been fueled most notably by potential applications of high-Tc superconductors (Kordyuk, 2012), graphene (Dresselhaus, 2014), graphenelike structures (Butler et al., 2013), and layered heterostructures (Novoselov et al., 2005) to device engineering. The most remarkable outcome of this analysis is the realization that electronic states other than simple Fermi liquids occur in several materials families. Consequently, a wealth of theoretical and computational methods have been developed to describe accurately the electronic correlation in these low-dimensional many-fermion systems and to investigate unusual physical phenomena that appear in the presence of charge confinement or strong interaction. These methodological developments can be branched in two categories: selfconsistent mean-field approaches, such as Hartree-Fock (HF) and density functional, and approaches beyond mean-field, such as multi-configurational methods and manybody perturbation theory (MBPT) (Szabo and Ostlund, 1996; Onida et al., 2002). In this regard, the new "Fermi-Dirac" physics of graphene and its derivatives represents one of the most significant case studies of interlacing novel theories and experiments so far (Taioli et al., 2009a; Haberer et al., 2010, 2011;Umari et al., 2012;Taioli, 2014).

However, while these methods found extensive use in the specific research field for which they were tailored, mainly related to materials science, it is rather unusual to observe applications crossing the borders between different areas.

At variance, the goal of this opinion article is to provide a "bird's eye view" on a recently developed *ab initio* method (Taioli et al., 2009b,c, 2010), firstly implemented to perform spectral analysis of solids at drastically reduced computational cost. This method has been further extended, as well as successfully used, in apparently distant applications, ranging from materials science to astrophysics. Hereby, after a description of the conceptual framework behind this approach we will prove its ability to span different research fields by interpreting a specific set of physical problems in the meV to MeV range. In particular, we will discuss the Auger effect in SiO<sup>2</sup> (Taioli et al., 2009b), the BEC-BCS crossover in ultracold Fermi gases (UFG) (Simonucci et al., 2011, 2012; Garberoglio et al., 2013), and finally the electron capture in stellar plasma (Simonucci et al., 2013).

The physical picture, under which all the previous phenomena can be framed, is the existence of one (or many) intermediate, quasi-bound state(s) produced by an initial scattering event, embedded in the continuum represented by one or more fragments of the system under investigation.

Depending on the initial kinetic energy of the impinging projectile, the dynamics of the excited target system may undergo a resonant state before decaying. Correlation between system's many constituents in the intermediate state is the driving force to final decay, possibly leading to several outcomes characterized by a particular set of observables ("channels"). Multichannel outcomes are ultimately the fingerprint of the many-body nature of interactions. We stress that, while scattering events may be driven by different forces coupling the initial and final states, be it Coulomb as in the Auger decay or weak as in β-decay, nevertheless the common platform to all these processes is given by the discretecontinuum interaction, named after Fano-Fesbach (Fano, 1961). Scattering theory aims to find positive energy solutions (scattering states) of the following many-body Hamiltonian:

$$(\hat{H} - E)\Psi^-\_{a,\epsilon} = 0 \tag{1}$$

where the scattering wave function |Ψ − <sup>α</sup>,<sup>∈</sup> i describes the motion of a particle asymptotically not interacting with kinetic energy ∈<sup>α</sup> = *E* – *E*α, while the scattering center is in the state classified by the quantum numbers α. The Hamiltonian in Eq. 1 usually contains kinetic and interaction terms. For example in condensed matter systems the electronic Hamiltonian within the Born–Oppenheimer approximation is:

$$H = -\frac{\hbar^2}{2m} \sum\_{i=1}^{N} \nabla\_i^2 + \sum\_{i
$$-\sum\_{i=1}^{N} \sum\_{k=1}^{n} \frac{Z\_k e^2}{|r\_i - R\_k|} \tag{2}$$
$$

where, in order of appearance, we find the electron kinetic energy, the electron–electron, and electron-nuclei Coulomb interaction terms, while *N* and *n* are, respectively, the total number of electrons of mass *m* at positions *r<sup>i</sup>* and of nuclei with charge *Z<sup>k</sup>* at positions *Rk.* Fano's approach basically provides a useful way to write the scattering wavefunction by mixing the discrete quasi-bound spectrum |Φi of the Hamiltonian with *N<sup>c</sup>* continuum channels |χ − β,τ i as follows:

$$\begin{split} \vert \Psi\_{a,\epsilon}^{-} \rangle &= a\_{\alpha}(\in) \vert \Phi \rangle \\ &+ \sum\_{\beta=1}^{N\_{\epsilon}} \int\_{0}^{\infty} \vert \chi\_{\emptyset,\mathfrak{r}}^{-} \rangle C\_{\beta,a}(\mathfrak{r},\mathfrak{e}) d\mathfrak{r} \tag{3} \end{split} (3)$$

where *a*<sup>α</sup> and *C*β*,*<sup>α</sup> are the expansion coefficients. For example, in the case of the Auger process the final channel |χ − β,τ i represents an ionized target in the state β plus an electron in the continuum with kinetic energy τ.

In this regard, to determine the cross section for a given process is more suitable for implementations to write the manybody equation of motion for quantum particles (Eq. 1) in terms of the scattering *T*matrix. The *T*-matrix is defined as (Taioli et al., 2010):

$$\begin{aligned} \hat{T} &= (E - \hat{H})\hat{G}(E)\hat{V} + \hat{V}\hat{G}(E)\hat{V} \\ &= (E - \hat{H} + \hat{V})\hat{G}(E)\hat{V} \\ &= (E - \hat{H}\_0)\hat{G}(E)\hat{V} \end{aligned} \tag{4}$$

where *G*ˆ (*E*) = (*E* − *H*ˆ + *i* ∈) −1 is the Green's function of the full Hamiltonian, *V* is the interaction potential*, H* and *H*<sup>0</sup> are the full and non-interacting Hamiltonian (in Eq. 1 *H*<sup>0</sup> is the sum of kinetic energy and electron-nuclei potential) and *E* is the total energy.

This quantity calculated between initial (*i*) (usually free plane wave) and final (*f*) scattering states (Eq. 3) is a function of the possible ingoing and outgoing channels as well as of the many-body interaction potential within the scattering region. Furthermore, *T* is directly relatable to experiments via the scattering cross section as follows:

$$
\sigma\_{\rm if} \propto |\langle f \vert T \vert i \rangle|^2 \delta(E\_{\rm i} - E\_{\rm f}) \qquad (5)
$$

and to other observables such as partial decay rates into reaction channels, transition life-times, and energy levels of the target (e.g.,in the Korringa-Kohn-Rostoker multiple scattering scheme (Papanikolaou et al., 2002)).

A first example of application of this multichannel scattering approach is concerned with electron spectroscopies. These techniques are based on scattering processes in which the initial state consists of a projectile, typically represented by a photon or an electron, that collides with a target, and final decay states are characterized by the presence of a few fragments recorded by some spectrometers.

Depending on the beam incident energy, the target system may undergo several processes: (i) direct emission, whereby the electron is imparted sufficient energy to be emitted from the solid; (ii) absorption, which consists in the creation of an electron-hole pair to probe the unoccupied density of states; (iii) autoionization and Auger emission, which are the results of non-radiative transitions from neutrally excited (autoionization) or singly ionized (Auger process) systems as a consequence of the impact; (iv) energy loss, when monochromatic electron beams undergo inelastic interactions, change momentum and lose energy. In these processes, the system often goes through intermediate discrete quasi-bound states with lifetimes longer than the collision times, decaying by fragmentation afterward. In this case, the driving force leading to decay is the Coulomb interaction among system's constituents. Observation of fragment types, as well as their angular and energy distribution, provides useful information on the many-particle system's properties.

Historically, the most fruitful approach toward electron spectroscopy has been to understand measured spectra within a one-electron picture (Taioli et al., 2010). This amounts to hide electron–electron interactions by introducing the concept of quasi-particle, defined as a real particle dressed by its interaction with other electrons in the system. Whereby this picture holds, measured electron spectra and related observables (e.g., binding energies) are thus referred to quasi-particles; in this case, mean-field approaches, such as HF and DFT, are generally found in qualitative agreement with experiments.

Alternatively, one could also view electron spectra as a means to learn about electron–electron interactions. Electron correlation effects on spectral lineshapes are observed with respect to the oneelectron picture through the appearance of:


Over the last few years, it has become possible to obtain electron spectra with a total energy resolution (monochromator plus electron spectrometer) considerably smaller than the linewidth of the core/valence-hole level investigated (Taioli et al., 2010). This fact allows accurate determination of the intensity distribution also for those transitions which are split into a large number of sub-levels due to the coupling of the emitted electrons with spectator electrons, nuclei and, collective excitations (plasmons). For a quantitative interpretation of electron spectra recorded on solids, problems mainly arise in the accurate treatment of the electronic correlation, in the localization of the electron emission, in the degree of localization over the valence band, and in the energy losses suffered by the electrons in their way out of the solid.

To take into account both extrinsic and intrinsic features on the same grounds, we proposed a unified method able to calculate electron spectra in solids (Taioli et al., 2009b,c, 2010) using beyond meanfield techniques. This method provides an extension of Fano's resonant multichannel scattering, allowing several intermediate quasi bound states (see Eq. 3) and solving the scattering problem by calculating the *T*-matrix (see Eq. 4) with appropriate boundary conditions. Within this approach, the discrete-continuum interaction modifies the spectral lineshape by including the interchannel interaction between electronic orbitals calculated from beyond mean-field first-principles techniques, while Monte Carlo is used to simulate the effect of inelastic losses on the original line-shape.

With the last step, the theoretical intrinsic spectrum is convoluted with extrinsic losses, so that it can be directly compared to "as-acquired" experimental spectra (Shirley, 1972; Tougaard, 1986). Moreover, shake processes, electron-hole and hole-hole interactions can be considered without semiempirical parameters appearing in Hubbard-like models (Lander, 1953; Cini, 1976; Sawatzky, 1977). An example of application of this method to the calculation of the Si K-LL Auger spectrum in SiO<sup>2</sup> is reported in **Figure 1A**. We stress that the effect of the multichannel nature of the scattering and of the continuum-discrete interaction are both visible respectively in the variety of transitions and in the asymmetric

Related, but in some ways even more puzzling, are the properties of UFG (Giorgini et al., 2008), which represent the ideal playground for studying correlations within a scattering framework. Fermions repel each other due to the Pauli principle. However, recent experimental developments (Regal et al., 2003) have opened the possibility of tuning the interaction between fermionic atoms, typically alkali species, using external magnetic fields. In this way, the scattering length can be varied to values much larger than the mean interatomic distance, by forcing the system to undergo a Feshbach resonance. Under these conditions, resonant scattering between two particles may occur, in which the bound channel becomes energetically close to the continuum spectrum threshold, and may interact via discrete-continuum interaction. The theoretical methods used to describe this "unitary" regime, which is a crossover

lineshape.

**Frontiers in Materials** | Mechanics of Materials April 2015 | Volume 2 | Article 33 | **32**

tion, are mainly based on Quantum Monte Carlo and mean-field approaches (Giorgini et al., 2008), notably the Bogoliubov-de Gennes (BdG) equations. There are two ways in which multi-

fermions in BCS theory and strongly coupled bosons in Bose–Einstein condensa-

channel scattering theory can be applied to study this regime. The first is the calculation of the scattering length describing the interaction strength between alkali atoms. Secondly, we demonstrated (Simonucci et al., 2011, 2012; Garberoglio et al., 2013) that the BdG equations used to describe this system (and thus the pairing function ∆) can be naturally written in term of the *T*-matrix of the underlying potential. Many-body scattering techniques, as for the case of the Auger process, can thus be applied to the assessment of the UFG properties. In particular, using this *T*-matrix based approach one can avoid issues related to the divergence of the interaction arising from the contact potential and, moreover, assess the effect of the actual finite range of the potential on the observables in and out the universality regime (Simonucci et al., 2011, 2012; Garberoglio et al., 2013). In **Figure 1B** we show the application of our approach to the calculation of the main observables in a <sup>6</sup>Li UFG as a function of the interaction effective range appearing in the second order expansion of the *T*-matrix (Simonucci et al., 2011).

To appreciate further more the general nature of our theoretical framework, we applied this multichannel scattering approach to bound-continuum interaction

for determining the abundance of Li in evolved stars. This is a crucial problem to assess Big-Bang nucleosynthesis of Li and still an open issue in astrophysics (Adelberger et al., 2011). Indeed Li production and depletion in stars is not accounted for in details and is particularly hampered by poor knowledge of how β-decay rates are affected in the rapidly varying temperature and density conditions below the envelope of stars (Adelberger et al., 2011). Extrapolations from solar data are extremely uncertain and an accurate estimate of the electron capture on Be is needed. By applying our multichannel approach, we were able to demonstrate that the *e* – capture is proportional to the weak force *T*-matrix, and, assuming the weak force as a contact interaction,to the electron density at the nucleus ρ<sup>0</sup> (Simonucci et al., 2013).

Due to the high temperatures and densities found in stellar plasma (ρ > 150 g/cm<sup>3</sup> and *T* 10<sup>7</sup> K), *e* – capture may occur from both bound and continuum states. Our formalism allows one to calculate accurately the electron density at the nucleus by including bound, partially ionized, and continuum channels into the electron capture estimate on the same grounds and, thus, to go beyond the previous approaches relying on the Debye–Huckel approximation to the electronic screening (Bahcall, 1962).

We stress that the changes in the electron density at the nucleus with respect to previously used approximations (ρ<sup>0</sup> and, thus, half-life as different as 30% with respect to that obtained by Bachall (Bahcall, 1962)

for high density, low temperature, i.e., quantum regime) and an enhanced electron screening, which helps thermonuclear fusion, might have rather wide astrophysical consequences (e.g., for *p*-captures on Li). In **Figure 1C,** we report the new β-decay rate as a function of *T* and ρ.

Summarizing, in this opinion article I discussed a recently developed many-body method, based on the concept of Fano's interaction and the use of the scattering *T*-matrix, able to interpret a variety of physical phenomena, ranging from molecular and solid state physics to astrophysical scenarios.

#### **ACKNOWLEDGMENTS**

This work has received funding from the EU seventh Framework Programme under grant agreement n. 604391 Graphene Flagship. ST acknowledges support by Istituto Nazionale di Fisica Nucleare through the Supercalcolo agreement with Bruno Kessler Foundation, from the European Science Foundation under the INTELBIO-MAT Exchange Grant and FBK for providing economical support through the "research mobility scheme". Finally, ST acknowledges the Institute of Advanced Studies in Bologna for his ISA fellowship.

#### **REFERENCES**


I. The half-life of <sup>7</sup>Be in evolved stars. *Astrophys. J.* 764, 118–129. doi:10.1088/0004-637X/764/2/118


**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 22 January 2015; accepted: 28 March 2015; published online: 20 April 2015.*

*Citation: Taioli S (2015) A bird's eye view on the concept of multichannel scattering with applications to materials science, condensed matter, and nuclear astrophysics. Front. Mater. 2:33. doi: 10.3389/fmats.2015.00033*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials.*

*Copyright © 2015 Taioli. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Mermin differential inverse inelastic mean free path of electrons in polymethylmethacrylate

## **Maurizio Dapor \***

European Centre for Theoretical Studies in Nuclear Physics and Related Areas (ECT\*-FBK), Trento Institute for Fundamental Physics and Applications (TIFPA-INFN), Trento, Italy \*Correspondence: dapor@ectstar.eu

**Edited by:**

Federico Bosia, University of Torino, Italy

#### **Reviewed by:**

Dimitris Emfietzoglou, University of Ioannina, Greece Rafael Garcia-Molina, Universidad de Murcia, Spain

**Keywords: dielectric theory, transport of fast electrons in solids, differential inverse inelastic mean free path, polymethylmethacrylate, polymers**

# **INTRODUCTION**

Electrons continually interact with the matter around us. We use electron beams for our purposes, either on the side of production of materials or on that of their characterization. Let us think to the many applications such as the processing of materials with plasma and to the local melting of materials for joining large components. We use electron beams also in the electron lithography, an important technique utilized for the production of the microelectronics devices. Let us consider the importance of the beams of electrons in the characterization of the materials, performed using techniques such as the electron microscopy and all the electron spectroscopies. Electrons interact with the surfaces of the space-crafts. The plasma– wall interaction in the fusion reactors also involves electron–matter interaction. Electrons play a role also in the cancer proton therapy, where cascades of secondary electrons are produced. These electrons of very low energy are toxic for the human body cells, since they produce damage to the biomolecules due to ionizations/excitations and the resulting break of chemical bonds. Also the secondary electrons which have ultra-low energies – and which, for a long time, were thought to be relatively harmless – are dangerous for the biomolecules due to the so-called "dissociative electron attachment." And, of course, we wish to minimize the effects of the irradiation on the healthy tissues near to the diseased cells.

In all the mentioned cases, the modeling of the interaction of the electrons with the matter is very important, as it can be used to provide a solid theoretical interpretation of the experimental evidences.

The interpretation of the experimental results is often based on the Monte Carlo method (MC). The doping contrast in secondary-electron emission of pnjunctions was investigated by the use of the MC method (Dapor et al., 2008; Rodenburg et al., 2010). The modeling of electron interaction with polymers, and in particular with the polymethylmethacrylate, has been demonstrated to be very important in nano-metrology. Line-scan of resist materials with given geometrical cross-sections deposited on silicon or silicon dioxide, and the corresponding linewidth measurements, obtained with the secondary electrons in the Critical-Dimension Scanning Electron Microscope (CD SEM), require an interpretation that can be performed using MC calculations (Dapor et al., 2010).

The MC method, in turn, requires an accurate description of the differential inverse inelastic mean free path (DIIMFP), in order to calculate the inelastic mean free path (IMFP) and the distribution function for inelastic collisions of electrons causing energy losses less than or equal to given values.

# **MERMIN THEORY**

The Mermin dielectric function (Mermin, 1970) is given by

$$\begin{aligned} \varepsilon\_{\mathsf{M}}\left(\boldsymbol{q},\boldsymbol{\omega}\right) \\ = 1 + \frac{\left(1 + i/\boldsymbol{\omega}\boldsymbol{\tau}\right)}{\left[\boldsymbol{\varepsilon}^{0}\left(\boldsymbol{q},\boldsymbol{\omega} + i/\boldsymbol{\tau}\right) - 1\right]}, \; \text{(1)} \\ = \frac{\left[\boldsymbol{\varepsilon}^{0}\left(\boldsymbol{q},\boldsymbol{\omega} + i/\boldsymbol{\tau}\right) - 1\right]}{\left[\boldsymbol{\varepsilon}^{0}\left(\boldsymbol{q},\boldsymbol{0} + 1\right) - 1\right]}, \; \text{(1)} \end{aligned} , \; \text{(1)}$$

where *q* is the momentum, ω the frequency, τ the relaxation time, and ε 0 (*q*, ω) the Lindhard dielectric constant (Lindhard, 1954)

$$\varepsilon^{0}(\mathbf{q},\omega) = \mathbf{l} + \frac{4\pi e^{2}}{q^{2}}B(\mathbf{q},\omega),\tag{2}$$

*B*(*q*, ω)

$$=\int \frac{d\mathbf{p}}{4\pi^3} \frac{f\_{\mathbf{p}+\mathbf{q}/2} - f\_{\mathbf{p}-\mathbf{q}/2}}{\alpha - (\epsilon\_{\mathbf{p}+\mathbf{q}/2} - \epsilon\_{\mathbf{p}-\mathbf{q}/2})/\hbar}.\tag{3}$$

In these equations *e* is the electron charge, *f<sup>p</sup>* is the Fermi–Dirac distribution, and e*<sup>p</sup>* the free electron energy.

### **ENERGY LOSS FUNCTION**

Let us now consider a superposition of free and bound oscillators. For any oscillator, the energy loss function (ELF) is given by the opposite of the imaginary part of the inverse of the Mermin dielectric function:

$$\operatorname{Im}\left[\frac{-1}{\varepsilon\_{\mathsf{M}}\left(\alpha\_{\mathrm{i}},\chi\_{\mathrm{i}};q,\alpha\right)}\right] = \frac{\varepsilon\_{\mathsf{M}\_{2}}}{\varepsilon\_{\mathsf{M}\_{1}}^{2} + \varepsilon\_{\mathsf{M}\_{2}}^{2}}.\tag{4}$$

where

$$\mathbf{e}\_M = \mathbf{e}\_{M\_1} + i\mathbf{e}\_{M\_2} \tag{5}$$

and ω<sup>i</sup> , and γ<sup>i</sup> are, respectively, the frequency and the damping constant associated to each specific oscillator. A linear combination of Mermin-type ELFs, one per oscillator, allows to calculate the electron ELF for *q* = 0, for any specific material

(Planes et al., 1996; Abril et al., 1998; de Vera et al., 2011)

**kinetic energyT in the range 10–1000 eV**.

$$\begin{aligned} &\operatorname{Im}\left[\frac{-1}{\varepsilon\left(q=0,\omega\right)}\right] \\ &=\sum\_{\mathrm{i}}A\_{\mathrm{i}}\operatorname{Im}\left[\frac{-1}{\varepsilon\_{\mathrm{M}}\left(\omega\_{\mathrm{i}},\gamma\_{\mathrm{i}};q=0,\omega\right)}\right]. \end{aligned} (6)$$

In this equation,*A*<sup>i</sup> , ω<sup>i</sup> , and γ<sup>i</sup> are determined looking for the best fit of the available experimental optical ELF. Actually, as both Mermin and Drude–Lorentz oscillators converge on the same values in the optical limit (i.e., for *q* = 0) (de la Cruz and Yubero, 2007),

$$\begin{aligned} \operatorname{Im}\left[\frac{-1}{\varepsilon\left(q=0,\,\omega\right)}\right] \\ &=\sum\_{\mathrm{i}}A\_{\mathrm{i}}\operatorname{Im}\left[\frac{-1}{\varepsilon\_{\mathrm{M}}\left(\omega\_{\mathrm{i}},\gamma\_{\mathrm{i}};q=0,\,\omega\right)}\right] \\ &=\sum\_{\mathrm{i}}A\_{\mathrm{i}}\operatorname{Im}\left[\frac{-1}{\varepsilon\_{\mathrm{D}}\left(\omega\_{\mathrm{i}},\gamma\_{\mathrm{i}};q=0,\,\omega\right)}\right], \end{aligned} \tag{7}$$

where the Drude–Lorenz functions are given by (Ritchie and Howie, 1977)

$$\begin{split} \operatorname{Im} \left[ \frac{-1}{\varepsilon\_D(\omega\_i, \gamma\_i; q = 0, \alpha)} \right] \\ = \frac{\chi\_i \alpha}{(\omega\_i^2 - \omega^2)^2 + (\gamma\_i \alpha)^2}, \qquad (8) \end{split} $$

the best fit can also be performed using a linear combination of Drude–Lorentz functions, instead of Mermin functions.

For the present work, the parameters provided by de Vera et al. (2011) obtained by calculating the best fit of the Ritsko et al. (1978) experimental optical data, were used.

Once the values of the best fit parameters have been established, the extension out of the optical domain (*q* 6= 0) can be obtained by (Planes et al., 1996; Abril et al., 1998; de Vera et al., 2011)

$$\begin{aligned} &\operatorname{Im}\left[\frac{-1}{\varepsilon\left(q,\,\omega\right)}\right] \\ &=\sum\_{i}A\_{i}\operatorname{Im}\left[\frac{-1}{\varepsilon\_{\mathcal{M}}\left(\boldsymbol{\omega}\_{i},\boldsymbol{\gamma}\_{i};q,\,\omega\right)}\right]. \end{aligned} (9)$$

#### **DIFFERENTIAL INVERSE INELASTIC MEAN FREE PATH AND INELASTIC MEAN FREE PATH**

The ELF allows, in turn, the computation of the DIIMFP, given by

$$\frac{d\lambda\_{\text{inel}}^{-1}}{d\hbar\omega} = \frac{1}{\pi a\_0 T} \int\_{q-}^{q+} \frac{dq}{q} \text{Im} \left[ \frac{-1}{\varepsilon(q, \omega)} \right],\tag{10}$$

where *a*<sup>0</sup> is the Bohr radius,*T* is the kinetic energy of the incident electrons and

$$q\_{\pm} = \sqrt{\frac{2m}{\hbar^2}} \left( \sqrt{T} \pm \sqrt{T - \hbar \alpha} \right). \quad \text{(11)}$$

The Mermin DIIMFP of electrons in PMMA is represented in **Figure 1**, for kinetic energies of the incident electrons in the range from 10 to 1000 eV.

The inverse of the integral of every curve presented in **Figure 1** provides, for each kinetic energy *T*, the IMFP. According to de la Cruz and Yubero (2007), the values of the IMFP calculated using the Tanuma, Powell, and Penn (TPP) empirical predictive formula (Tanuma et al., 1994) are systematically higher than the corresponding values calculated within the Mermin theory. For PMMA, according to our calculation, when *T* = 100 eV, the Mermin IMFP is equal to 6.3 Å while when *T* = 1000 eV, it is equal to 27.6 Å. According to TPP, the IMFP of PMMA is equal to 7.9 Å for *T* = 100 eV and to 33.7 Å for *T* = 1000 eV (Tanuma et al., 1994). Also approaches based on the Drude–Lorentz theory (Emfietzoglou et al., 2013; Dapor, 2014a,b; Dapor et al., 2015) provide values of the IMFP systematically higher than those obtained using the Mermin theory. The IMFP of PMMA calculated according to the Drude–Lorentz theory is equal to 10.1 Å for*T* = 100 eV and to 33.5 Å for *T* = 1000 eV (Dapor, 2014a).

It is not simple to express an opinion about the different approaches used today to calculate the DIIMFP, and hence the IMFP, due to the lack of experimental data and their quite large uncertainty (Dapor et al., 2015). Roberts et al. (1980) provided the experimental values of 29 ± 4 Å for 1196 eV electrons and 33 ± 5 Å for 1328 eV electrons in PMMA. On the one hand, also taking into account of the experimental uncertainty, de Vera et al. (2011) have shown that these values are closer to those predicted by the Mermin theory than to those predicted by the TPP formula,usually taken as a reference and provided, among other IMFP data, by the NIST database. On the other hand, the number of experimental data about electron IMFP in PMMA seems to be definitely too small to judge. Since Mermin theory is much more accurate it is preferable, even if we conclude that, at the moment, both Mermin theory and TTP empirical predictive formula are compatible with the available experimental data about IMFP in PMMA.

#### **CONCLUSION**

Calculations of the DIIMFP of electron in materials are of paramount importance for modeling the transport of electrons in solid targets. In this paper, we presented the calculations of the DIIMFP of electrons in PMMA based on the Mermin theory (Mermin, 1970) and on the de Vera et al. (2011) best fit of the Ritsko et al. (1978) optical data.

#### **ACKNOWLEDGMENTS**

The author wishes to express his gratitude to Isabel Abril, Pablo de Vera, and Rafael Garcia-Molina for their invaluable suggestions. This work was supported by Istituto Nazionale di Fisica Nucleare (INFN) through the Supercalcolo agreement with FBK.

#### **REFERENCES**


glassy carbon, amorphous carbon, and diamond: comparison between two different dispersion laws. *Nucl. Instrum. Methods Phys. Res. B* doi:10.1016/j. nimb.2014.11.115


1 to 300 eV. *J. Chem. Phys.* 69, 3931. doi:10.1063/1. 437131


**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 02 February 2015; accepted: 16 March 2015; published online: 02 April 2015.*

*Citation: Dapor M (2015) Mermin differential inverse inelastic mean free path of electrons in polymethylmethacrylate. Front. Mater. 2:27. doi: 10.3389/fmats.2015.00027*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials.*

*Copyright © 2015 Dapor. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Hydrogen storage in rippled graphene: perspectives from multi-scale simulations

## **Vito Dario Camiola1,2, Riccardo Farchioni 1,2,3,Tommaso Cavallucci 2,3, Antonio Rossi 2,4, Vittorio Pellegrini 1,2,4 and Valentina Tozzini 1,2\***

1 Istituto Nanoscienze, Consiglio Nazionale delle Ricerche, Pisa, Italy

<sup>2</sup> Laboratorio NEST: National Enterprise for nanoScience and nanoTechnology, Scuola Normale Superiore, Pisa, Italy

<sup>3</sup> Physics Department Enrico Fermi, University of Pisa, Pisa, Italy

<sup>4</sup> Graphene Labs, Istituto Italiano di Technologia (IIT), Genova, Italy

\*Correspondence: valentina.tozzini@nano.cnr.it

#### **Edited by:**

Simone Taioli, Bruno Kessler Foundation, Italy

#### **Reviewed by:**

Silvio A. Beccara, Bruno Kessler Foundation, Italy

**Keywords: density functional theory, classical molecular dynamics, rippled graphene, flexural phonons, hydrogen storage**

Exploring new perspectives for green technologies is one of the challenges of the third millennium, in which the need for non-polluting and renewable powering has become primary. In this context, the use of hydrogen as a fuel is promising, since the energy released in its oxidation per unit mass (~142 MJ/kg) is three times that released, on average, by hydrocarbons, and the combustion product is water (Ramage, 1983). Being hydrogen a vector of chemical energy, efficient conservation, and non-dispersive transportation are the main goals. Three issues must be considered to this respect: (i) storage capacity, (ii) storage stability, and (iii) kinetics of loading/release. Commercial technologies are currently based on cryo-compression or liquefaction of H<sup>2</sup> in tanks. These ensure quite a high gravimetric density [GD, point (i)], namely 8–13% in weight of stored hydrogen, and a relatively low cost (Züttel, 2003). However, concerning points (ii) and (iii), these technologies pose problems of safety, mainly due to explosive flammability of hydrogen, and consequent unpractical conditions for transportation and use (Mori and Hirose, 2009). Therefore, research efforts are directed toward solidstate based storage systems (energy.gov, Bonaccorso et al., 2015).

Interactions of hydrogen with materials are classified as physisorption, occurring with H<sup>2</sup> by means of van der Waals (vdW) forces, or chemisorption, i.e., chemical binding of H leading to the formation of hydrides (Mori and Hirose, 2009), requiring dissociative(associative)

chemi(de)sorption of H2. Intermediate nature interactions, sometimes called "phenisorption," can also occur between hydrogen electrons and the electrons of external orbital of metals. Indeed, stable and robust (light) metal hydrides (Sakintuna et al., 2007; Harder et al., 2011) are currently considered an alternative to tanks. Their main drawback is their high chemisorption and chemidesorption barrier, both many times the typical thermal energy, implying slow operational kinetics, which becomes acceptable only at very high temperature. Physisorption, conversely, generally results in barrierless and weak binding. It was considered as a storage mechanism in layered (Zhirko et al., 2007) or porous (Sastre, 2010) materials, and shown to be effective at low temperatures and/or high pressure. Therefore, it generally seems that if storage stability (ii) is improved then the loading/release kinetics (iii) is worsened.

Graphene shows good potential to be an efficient hydrogen-storage medium (Tozzini and Pellegrini, 2013): carbon is among the lightest elements forming layered and porous structures, and graphene is probably the material with the largest surface to mass ratio. These two conditions are in principle optimal to produce high GD [point (i)]. In addition, the chemical versatility of carbon allows it to interact with hydrogen both by physisorption (in sp2 hybridization) and chemisorption (Goler et al., 2013a) (in sp3 hybridization). ["Phenisorption is also obtained in graphene by functionalization with metals (Mashoff et al., 2013)].

On the other hand, concerning points (ii) and (iii), pure graphene does not perform dramatically better than other materials. H<sup>2</sup> easily physisorbs onto graphene layers or within multilayers, but it was theoretically shown (Patchkovskiim et al., 2005) that large GD (6–8%) are reached within multi-layered graphene at cryogenic temperatures, while the room temperature value is at best ~2–3%. This was confirmed by measurements (Klechikov et al., private communication), which also indicate that graphene does not perform better than other carbon based bulk materials, such as nanoporous carbon or carbon nanotubes. In all cases, a key parameter determining GD is the specific surface to volume ratio. Theoretical works also show that stability can be improved (and GD optimized) at specific interlayer spacing (~7–8 Å), due to a cooperative effect of vdW forces (Patchkovskiim et al., 2005). A similar effect is responsible for the accumulation of physisorbed hydrogen within graphene troughs at low temperatures (~100 K) observed in simulations (Tozzini and Pellegrini, 2011). On the other hand, hydrogen chemisorption on graphene produces graphane (Sofo et al., 2007), its completely hydrated alkane counterpart, stable at room temperature [point (ii)] and with 8.2% GD [point (i)]. Graphane, however, shares with other hydrides high chemi(de)sorption barrier (~1.5 eV/atom). As in other materials, physisorption has good kinetics (iii), and

bad storage capacity (i), and stability (ii), while chemisorption has good (i) and (ii), and bad (iii).

Graphene, however, displays extremely peculiar properties: a unique combination of strength and flexibility, and its bidimensionality. Therefore, it can be buckled on different scales down to nanometer (Fasolino et al., 2007; Wang et al., 2011), either statically, i.e., forming stable ripples as an effect of external constraints [i.e., compression or interaction with a substrate (Goler et al., 2013b)], or dynamically, sustaining traveling ripples, i.e., coherent transverse out-of-plane acoustic modes, also called ZA or flexural phonons (Lindsay and Broido, 2010; Xu et al., 2013). We explored the possibility of exploiting these properties in the context of hydrogen storage.

We first evaluated the dependence of chemisorbed H stability on the local curvature by means of a density functional theory (DFT) based study (Tozzini and Pellegrini, 2011). This revealed a linear dependency of the binding energy on the local curvature, leading to an over-stabilization of the adsorbate of up to 1–2 eV on crests and corresponding destabilization within troughs (**Figure 1A**),which is also modified by the presence of already chemisorbed hydrogen (Rossi et al., submitted). Theoretical evaluations were confirmed by Scanning Tunneling Microscopy performed on naturally corrugated graphene grown on SiC (Goler et al., 2013a), which demonstrated the presence of hydrogen prevalently on the convexities. These observations lead to the idea that an inversion of curvature (from concavity to convexity) could detach chemisorbed hydrogen. As a matter of fact, this mechanism was demonstrated in a simulation in which the curvature inversion was realized dynamically by the passage of a ZA coherent phonon of ~2 nm wavelength and ~THz frequency (**Figure 1B**). The simulation shows that hydrogen undergoes associative desorption at room temperature upon curvature inversion, which can be called a "mechanical catalysis": the energy barrier is overcome by means of the energy provided by the traveling phonon.

Once detached, H<sup>2</sup> finds itself in contact with a graphene sheet with traveling ripples. If sheets are at appropriate distances and if ZA phonons are excited with appropriate relative phases in subsequent sheets, the traveling ripples enclose traveling nano-cavities (**Figure 1C**). We showed by means of classical molecular dynamics (MD) simulations using empirical force fields (FF) (Camiola et al., 2014) that these nano-cavities can include and transport H<sup>2</sup> at velocities near to the phase velocity of the phonon for almost micrometric distances before the phonon damps. This effect could be used to transport and pump hydrogen through multilayers, improving the GD of physisorbed hydrogen at room temperature and the kinetics of loading and release.

The existence of ripples and flexural phonons, therefore, offers specific strategies – unique to graphene – to overcome issues (i), (ii), and (iii). The proof of principle of the effectiveness of these strategies was given by means of a multi-scale approach, combining DFT and empirical FF based MD simulations, in order to Camiola et al. Hydrogen in rippled graphene

couple high accuracy in the representation of interactions to the large time and size scales needed to represent the process. However, the realization of these processes relies on the practical possibility of finely manipulating the local curvature and generating and sustaining coherent flexural phonons. In fact, our current efforts are along that route: we are exploring several possibilities for static and dynamic curvature manipulation, including external electric fields (Cavallucci, 2014), functionalization with optically active molecular pillars (Burress et al., 2010), electromechanical pulling, coupling to piezoelectric substrates (Camiola et al., 2014). In this phase, a multi-scale simulation approach will be even more necessary: in the route toward the practical realization of possible devices, one must analyze also macroscopic effects of thermodynamic nature, in addition to the physical phenomena at the sub-nano scale and the dynamical behavior at the nano-to-micro scale. Therefore, we intend to add a continuum representation, where graphene is modeled as a membrane with mechanoelastic properties (Zang et al., 2013) to the DFT and to the MD with empirical FF. Consistency between the representations will be obtained by appropriate bottom-up parameterization and topdown transfer of macroscopic information, and will ensure a complete and accurate representation of the system properties and behavior.

#### **ACKNOWLEDGMENTS**

We gratefully acknowledge financial support by the EU, 7th FP, Graphene Flagship (contract no. NECT-ICT-604391), the CINECA award "ISCRA C" IsC10\_HBG, 2013 and PRACE "Tier0" award Pra07\_1544, and IIT Compunet Platform for computational resources.

#### **REFERENCES**

Bonaccorso, F., Colombo, L., Yu, G., Stoller, M., Tozzini, V., Ferrari, A. C., et al. (2015). Graphene, related two dimensional crystals, and hybrid systems for energy conversion and storage. *Science* 347, 1246501. doi:10.1126/science.1246501


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 15 December 2014; accepted: 08 January 2015; published online: 30 January 2015.*

*Citation: Camiola VD, Farchioni R, Cavallucci T, Rossi A, Pellegrini V and Tozzini V (2015) Hydrogen storage in rippled graphene: perspectives from multi-scale simulations. Front. Mater. 2:3. doi: 10.3389/fmats.2015.00003 This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials.*

*Copyright © 2015 Camiola, Farchioni, Cavallucci, Rossi, Pellegrini and Tozzini. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# MATERIALS

OPINION ARTICLE published: 10 February 2015 doi: 10.3389/fmats.2015.00006

# Emergent models for artificial light-harvesting

## **Celestino Creatore<sup>1</sup>\*, AlexW. Chin<sup>1</sup>\*, Michael A. Parker <sup>1</sup> and Stephen Emmott <sup>2</sup>**

<sup>1</sup> Cavendish Laboratory, University of Cambridge, Cambridge, UK

<sup>2</sup> Computational Science Laboratory, Microsoft Research, Cambridge, UK

\*Correspondence: cc619@cam.ac.uk;ac307@cam.ac.uk

#### **Edited by:**

Simone Taioli, Bruno Kessler Foundation, Italy

#### **Reviewed by:**

Michele Casula, Université Pierre et Marie Curie, France

**Keywords: solar cells, artificial photosynthesis, quantum effects in biology, quantum coherence, artificial nanostructures**

### **INTRODUCTION AND BACKGROUND**

Photovoltaic (PV) solar energy conversion is one of the most promising emerging renewable energy technologies and is likely to play an ever-increasing role in the future of electricity generation. The potential of photovoltaics is considered to be much larger than the current energy demand, which is equivalent to ≈13,000 millions tons of oil (BP, 2014). Together with other renewable energy sources, photovoltaics could thus considerably aid the decarbonization of the energy sector in the near future, an objective, which is (should be) seriously sought after to meet the IPCC 2°C target<sup>1</sup> . Indeed, PV energy conversion is an ever-increasing research field having an impact on other sectors: technological (e.g., dealing with energy storage and distribution issues) and economical (e.g., evaluating the costs of PV devices compared to conventional energy technologies).

The typical working principle of a photovoltaic cell is quite simple, and well understood in terms of doped inorganic semiconductors (Wurfel, 2009). Photovoltaic technology has significantly moved forward since the first silicon (Si) solar cell having an efficiency of ≈4% was invented in 1950s in Bell Laboratories, following the breakthrough of the p–n junction developed by Shockley, Brattain, and Bardeen (Chapin et al., 1954). Since then, research on PV has led to singlejunction and multi-junction solar cells with record efficiencies in controlled (laboratory) environments of ≈29% (Green et al., 2014) and 46% (Fraunhofer-ISE, 2014), respectively. However, the efficiency of the most popular technologies in the commercial market – which belong to the

so called first generation of solar cells, mainly based on crystalline Si – is in the 10–18% range and around 25% in laboratory. A second generation of cells, mainly aimed at reducing the fabrication costs, is based on thin film architectures, where light is absorbed and charge generated in a solid thin layer of semiconductor (e.g., CdTe, CIGS), with efficiencies restrained to around 22%. Beyond these technologies, a variety of other PV designs has emerged, mainly driven by the progress on materials research and the need of low-cost manufacture devices, organic photovoltaic (OPV) being one of the most interesting solutions. The latest scientific developments in this field have taken various directions (using different concepts, materials, and geometries), all having the potential to result in high-efficiency PV devices. Most current approaches, which have been termed third generation PV [see e.g., Green (2003)], explore schemes beyond the guidelines given by Shockley and Queisser (1961), which are based on the assumption that each photon absorbed above the band gap of the semiconductor photovoltaic material generates a single electron-hole pair. These methods include solar cells exploiting multiple-exciton generation in both inorganic and organic materials, hot carriers collection, up- and down-conversions.

Reviewing the extensive literature about the past and present PV schemes is beyond the scope of this *Opinion*, many excellent reviews having already succeeded in this (challenging) task. Here, we will illustrate and highlight our point of view on an emergent strategy for PV and artificial light-harvesting in general. This strategy holds onto an interdisciplinary framework encompassing biology and quantum mechanics. Ideas, which have been developed within this nascent trend, are inspired by the recent discovery of quantum effects in biological systems (Mohseni et al., 2014), and aim to model and reproduce the highly efficient reactions often exhibited by these systems using the tools of quantum mechanics and nanotechnology platforms. Novel theoretical proposals, for instance (Scully, 2010; Blankenship et al., 2011; Scully et al., 2011; Creatore et al., 2013; Zhang et al., 2015), have underlined the relevance of effects of quantum coherence to enhance the performance of photocells whose working principle is inspired by the architecture of biological light-harvesting complexes (LHCs).

#### **EMERGENT SCHEMES FOR ARTIFICIAL LIGHT-HARVESTING**

Mimicking photosynthesis is perhaps the most intuitive scheme for building photovoltaic devices: as photovoltaic solar cells convert sunlight into electricity, biological LHCs use photosynthesis to convert sunlight into chemical energy necessary for sustaining their life. Biological photosynthesis begins with an ultrafast sequence of photophysical events, known as light reactions, which occur via a diverse range of nanoscaled pigment–protein complexes (PPCs) found in photosynthetic plants, algae, and bacteria (Blankenship, 2008). There reactions proceed through the creation of molecular excited states (excitons) by the absorption of solar photons, then followed by intra- and/or inter-PPCexciton energy transfer (EET); finally, the excited state energy reaches a reaction center (RC) where it is used to initiate

<sup>1</sup>This target calls for a maximum global mean temperature rise of 2°C compared to the temperature of pre-industrial times.

charge separation, by releasing an energetic electron for the later stages of the photosynthetic chemistry, i.e., to start chemical reactions to obtain carbohydrates. Intriguingly, the efficiency of these light reactions can often be close to 100%, suggesting a carefully optimization of the deleterious processes (decay and charge recombination, excited state trapping), which plague the attempts to carry out similar processes in artificial solar cells.

Recent experimental and theoretical studies demonstrate that quantummechanical effects (e.g., coherent quantum dynamics) play a previously unexpected and relevant role in biological EET (Mohseni et al., 2014). Quantum effects have been shown to survive even at ambient temperatures and are possibly beneficial to the observed very high efficiencies exhibited by these biological structures. Consequently, there has been an everincreasing interest in understanding and modeling the nanoscaled structures and the related pigment–pigment interactions, which drive the natural light reactions so efficiently. Mimicking the physics of the photosynthetic PPCs could be extremely valuable for the future design of more efficient artificial light-harvesting structures based on the exploitation of quantum effects.

A powerful idea for recognizing the role played by quantum effects in the light reactions, one based on the fundamental principles of quantum mechanics, has recently been proposed by Scully (2010), a pioneer of theoretical quantum optics. By modeling the light reactions (photon absorption/emission, excitonic energy transfer) as a cyclic working process, analogous to a thermodynamic Carnot cycle, Scully has presented a formalism, which allows the performance of the light reactions to be analyzed in a way similar to standard thermodynamics heat engines. Within this framework, the solar photonic energy is described in terms of a "hot" thermal bath that is then converted through photon absorption to low-entropy useful work, i.e., the production of high energy electrons later used for the subsequent photosynthetic chemistry. The relevance of the quantum heat engine (QHE) approach is exemplified by the famous calculation of Shockley and Queisser (1961) who used similar thermodynamical considerations to

determine the fundamental maximum efficiency of the single band-gap (junction) semiconductor solar cells (Wurfel, 2009). More specifically, Scully has shown that in photocells consisting of elements that support quantum coherence, it is possible to suppress the emission process. Spontaneous emission of photons before the absorbed energy can be used to perform any work is indeed the most important as well as unavoidable loss mechanisms in optical energy extraction. Under opencircuit conditions, i.e., when no work current is extracted, the maximum voltage that can be obtained is determined by the detailed balance of absorption and spontaneous emission. Hence, breaking the detailed balance by quenching the spontaneous emission may lead to enhanced power extraction and current generation in photocells (Scully, 2010).

An appealing direction to explore will thus combine the utility and elegance of the quantum heat engine framework and the nanoscale architecture of the highly efficient photosynthetic complexes. As an example developed along these lines, we have recently presented a scheme of a photocell inspired by the molecular elements of biological RCs (Creatore et al., 2013). In photosynthesis, the RC consists of a "special pair" of photon-absorbing pigments (donors), which are coupled through an electron-accepting molecule to a work extraction stage.

The scheme we propose harnesses well understood and experimentally observed quantum effects, such as the formation of delocalized excited (exciton) states resulting from the long-range dipole–dipole interaction (excitonic coupling) between closely spaced pigments. In our RCinspired photocell, such dipolar interaction between the pair or light-absorbing donor molecules leads to the formation of new optically excitable states and quantum interference effects, which can enhance the photocell performance relative to a "classical" scheme in which these effects are not present (**Figure 1**, left panel). In particular, we have shown how dipolar couplings and related quantum interference effects can be tailored to increase the efficiency of charge transfer at the acceptor molecule, thus leading to enhanced light-to-current generation and power extraction at the reaction center. Although not a model of an

actual RC, this proposal shows what could be achieved with the same resources, i.e., a pair of absorbers, an electron acceptor, and the ability to "build" the electronic structure, the quantum states, and their interactions. A recent follow-up (Zhang et al., 2015) (proposed by the Group of S. Kais and one of the groundbreakers of quantum biology, G. S. Engel) of our model has highlighted the originality of the QHE strategy for realizing PV devices based on photosynthesis, finding even greater efficiency gains in coherent trimer structures. The relevance of a new class of bio-inspired quantum nanotechnology, however, is not just restrained to reproduce/enhance efficient operations observed in nature. Higgins et al. (2014), for instance, have shown theoretically that well-established quantum control techniques can be used to achieve an enhanced photon absorption (a phenomenon called "superabsorption") in simple ring nanostructures reminiscent of photosynthetic LHCs (**Figure 1**, right panel), whereas these natural systems themselves are not designed to exploit this effect. Superabsorption, which is the absorptive analog of quantum superradiance, has a range of applications including photon detection in the context of optical sensors and light-based power transmissions (Higgins et al., 2014).

These bio-inspired schemes are potentially realizable in artificial molecular light-harvesting systems through existing organic synthesis routes. A fundamental breakthrough in this direction has been recently achieved by Hayes et al. (2013), who have engineered an artificial synthetic system (consisting of a series of rigid synthetic heterodimers) supporting the persistent electronic coherences observed in biological systems and believed to be beneficial for the process of energy transfer. Another direction for the realization of these models, one based on solid-state nanotechnology architectures, is now viable due to the increasing level of experimental control over nano-optical systems gained in recent years, and involves the use of solid-state chromophores such as quantum dots, superconducting qubits, plasmonic nanostructures, and possibly hybrid combinations of the above. From the theoretical/computational point of view, advanced simulation techniques to complement the Hamiltonian/microscopic

**FIGURE 1 | Left panel**: the photosynthetic reaction center modeled in Creatore et al. (2013), consisting of two donor molecules (red spheres) and one acceptor molecule (green sphere). In **(A)**, the donor D<sup>1</sup> and D<sup>2</sup> are identical and uncoupled and both able to absorb and re-emit light (see the dipole moments µ<sup>1</sup> and µ2). In **(B)**, the dipolar excitonic coupling between D<sup>1</sup> and D<sup>2</sup> leads to new excitable donor states (X<sup>1</sup> and X2) resulting from the quantum superposition of the original donor states: X<sup>1</sup> is a bright state with optical dipole moment µ<sup>x</sup><sup>1</sup> , while X<sup>2</sup> is a dark state, which cannot absorb or emit photons. Under a variety of realistic and experimentally realizable parameters, the presence of a dark state lowers the role of deleterious radiative recombination processes, while increases the efficiency of charge transfer at the acceptor A. This leads to an enhanced current and power generation at the RC, compared to the uncoupled/classical scheme of **(A)**; see Creatore et al. (2013) for details. **Right panel**: a collection of N dipoles arranged in a ring structure reminiscent of the photosynthetic light-harvesting complex LH1. By using environmental quantum control techniques, Higgins et al. (2014) have shown that it is possible to achieve a regime (superabsorption), in which the N dipoles absorb incident photons with a rate proportional to N 2 , this effect being the absorptive analog of quantum superradiance. This model could be relevant to applications to solar energy harvesting, where the excitation induced by the absorption of solar photons are efficiently transferred from the ring structure to a central core absorber and then to a center converting such excitations into stored energy; see Higgins et al. (2014) for details. [Copied with permission from the original paper: Nat. Commun. 5, Article number: 4705 (2014).]

approaches (based on setting up and solving quantum optical master equations) used in the aforementioned proposals, are required. These techniques are useful to precisely assess the optical geometry and the stability of the engineered quantum states, including the role of environmental degrees of freedom (using, e.g., timedependent density matrix renormalization group methods) and details of additional states, e.g., charge-transfer states (using, e.g., tight binding models).

### **CONCLUSION**

Innovative ideas are needed to satisfy the ever-growing interest and expectations on solar energy conversion. In this *Opinion*, we have presented the prospects of an emergent novel strategy to harvest solar energy based on replicating mechanisms developed in nature, and realizable using the theoretical and experimental tools of quantum mechanics. Although at a nascent

stage, this framework is attracting an increasing attention and we think that the proposals here discussed suggest an encouraging trend for biologically inspired PV. Further, the technological capabilities are now fully developed for translating the theoretical schemes in experimental proofs of concept/ prototypes. Besides, we believe that results and solutions established along these new lines will be beneficial to the evolving and interconnected fields of non-equilibrium (Bustamante et al., 2005) and quantum thermodynamics (Gemmer et al., 2004; Kosloff, 2013) as well. Biological processes usually take place at the micro- and nano-scale and evolve in out-of-equilibrium conditions – energy transfer and the motion of molecular motors being paradigmatic examples. The same scenario applies to nanoscale devices, which, through precise quantum-mechanical control, could be employed to unravel the thermodynamic relations between the non-equilibrium evolution of molecular quantum states and the enhanced (or reduced) working efficiency of biological or biologically inspired processes harnessing these states.

## **REFERENCES**


coherence. *Proc. Natl. Acad. Sci. U.S.A.* 108, 15097–15100. doi:10.1073/pnas.1110234108


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 30 December 2014; accepted: 19 January 2015; published online: 10 February 2015.*

*Citation: Creatore C, Chin AW, Parker MA and Emmott S (2015) Emergent models for artificial lightharvesting. Front. Mater. 2:6. doi:10.3389/ fmats.2015. 00006*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials.*

*Copyright © 2015 Creatore, Chin, Parker and Emmott . This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Atomistic model of metal nanocrystals with line defects: contribution to diffraction line profile

# **Alberto Leonardi and Paolo Scardi \***

Department of Civil, Environmental and Mechanical Engineering, University of Trento, Trento, Italy

#### **Edited by:**

Simone Taioli, Bruno Kessler Foundation, Italy

#### **Reviewed by:**

Andrea Piccolroaz, University of Trento, Italy Shangchao Lin, Florida State University, USA

#### **\*Correspondence:**

Paolo Scardi, Department of Civil, Environmental and Mechanical Engineering, University of Trento, Via Mesiano 77, Trento 38123, Italy e-mail: paolo.scardi@unitn.it

Molecular Dynamics (MD) was used to simulate cylindrical Pd and Ir domains with ideal dislocations parallel to the axis. Results show significant discrepancies with respect to predictions of traditional continuum mechanics. When MD atomistic models are used to generate powder diffraction patterns, strong deviations are observed from the usual paradigm of a small crystal perturbed by the strain field of lattice defects.The Krivoglaz–Wilkens model for dislocation effects of diffraction line profiles seems correct for the screw dislocation case if most parameters are known or strongly constrained. Nevertheless the practical implementation of the model, i.e., a free refinement of all microstructural parameters, leads to instability. Possible effects of the experimental practice based on Line Profile Analysis are discussed.

**Keywords: x-ray powder diffraction, line profile analysis, nanocrystalline materials, dislocations, molecular dynamics simulation**

## **INTRODUCTION**

Line broadening has been used since the dawn of X-ray diffraction (XRD) to study the microstructure of crystalline phases. Besides the most basic information on the size and shape of crystalline domains, which in the simplest form is provided by the wellknown Scherrer formula (Scherrer, 1918; Klug and Leroy, 1974), lattice defects have also been extensively studied [e.g., see (Warren, 1990; Krivoglaz, 1996; Snyder et al., 1999; Mittemeijer and Scardi, 2004)].

Analytical techniques to extract information from diffraction line profiles range from simple integral-breadth methods [see Scardi et al. (2004) for a recent review], including accurate studies of isolated peak tails (Wilson, 1955; Groma, 1998; Groma and Székely, 2000) to more complex Fourier analysis (Warren and Averbach, 1950; Warren, 1990), underlying the most recent Whole Powder Pattern Modeling (WPPM) approach (Scardi and Leoni, 2002; Scardi, 2008). All methods for studying dislocations rely on the early studies by Wilson, Krivoglaz, and Wilkens (Wilson, 1955; Wilkens, 1970a,b; Krivoglaz et al., 1983). The last author, in particular, provided an analytical expression for the Fourier Transform (FT) of the diffraction line profile of crystalline domains containing dislocations (Wilkens, 1970a,b)

$$A\_{\{\text{hkl}\}}^D(L) = \exp\left[-2\pi^2 s^2 L^2 < \varepsilon\_{\{\text{hkl}\}}^2(L) > \right]$$

$$= \exp\left[-\frac{\pi b^2}{2} \rho \overline{C}\_{\text{hkl}} s^2 L^2 f^\*\left(\frac{L}{R\_\epsilon}\right)\right] \tag{1}$$

where *L* is the Fourier length (distance between scattering centers) and *s* the scattering vector (*s* = 2 sin(θ)/λ). Main parameters are the average dislocation density (ρ) Burgers vector modulus (*b*) and effective outer cut-off radius (*Re*), which is related to the extension of the effects of the dislocation strain field and, more generally, to dislocation interaction. The *f* \* is a smooth function of *L*/*R<sup>e</sup>* obtained by Wilkens in a heuristic way, to grant integrability of Eq. 1 (Wilkens, 1970b). The anisotropy of the elastic medium and of the dislocation strain field, which depends on the specific dislocation type (e.g., edge or screw) and slip system, is accounted for by the anisotropy or contrast factor *C*hkl. For powder diffraction, an average is used for all equivalent components of a diffraction line profile, which depends on crystal symmetry, and more specifically on the Laue group; given the elastic tensor components, *C*ij and *S*ij, the average contrast factor, *C*hkl, can be calculated for any desired crystalline phase and slip system huvzi {hkl} [e.g., see Martinez-Garcia et al. (2009).

The Krivoglaz–Wilkens approach, using Eq. 1 or similar approximations, is useful and easily implemented in the experimental data analysis; it has been extensively used in materials science (e.g., see work by Klimanek and Kuzel (1988), Kuzel and Klimanek (1988), Kuzel and Klimanek (1989), Ungar et al. (1998), Ungar (2008), Scardi and Leoni (2002), Scardi et al. (2007)] even if, as a matter of fact, it has never been fully validated. A few studies (Kamminga and Delhez, 2000; Kaganer and Sabelfeld, 2011;Kaganer and Sabelfeld, 2014) have tested Eq. 1 against numerical simulations, but the latter were based on the same continuum mechanics expressions for the dislocation strain field underlying Eq. 1, and in any case referred to rather idealized microstructures. An experimental validation, e.g., by Transmission Electron Microscopy (TEM) is not straightforward: quantitative TEM evidence is hard to obtain when the dislocation density is in the range of interest to powder diffraction (typically, above 10<sup>14</sup> m−<sup>2</sup> ), especially after extensive plastic deformation and in small domains.

A useful and quite different point of view can be provided by Molecular Dynamics (MD). MD simulations can realistically describe nanocrystalline domains and clusters with the detail of atomistic models, which can be used to generate powder diffraction patterns from known, designed microstructures (Bulatov et al., 1998; Jacobsen and Schiotz, 2002; Yamakov et al., 2002; Yamakov et al., 2004; Li et al., 2010). The present work is a first step toward this direction, to shed some light on the validity of Eq. 1 and analytical methods relying on it, and the more general meaning of diffraction from defected polycrystalline materials.

In the present paper, the strain field in cylindrical nanocrystals containing screw or edge dislocations is first discussed, moving from the traditional continuous mechanics (stress–strain) description to atomistic simulations; strain effects on XRD patterns from simulated single-crystals and powders are then discussed in terms of broadening of the line profiles. In the last section, a state-of-theart powder diffraction analysis is employed to investigate the simulated patterns, to assess validity of the Krivoglaz–Wilkens approach also considering the additional effects of surface relaxation.

#### **MATERIALS AND METHODS**

Pd and Ir nanocrystals containing line defects were simulated by MD using the LAMMPS code [Large-scale Atomic/Molecular Massively Parallel Simulator – (Plimpton, 1995)], implementing geometrical conditions similar to those underlying the derivation of Eq. 1, i.e., straight dislocation lines in cylindrical regions. First step was the generation of isolated edge or screw dislocations in bulk microstructures, followed by a stabilization of defectcontaining microstructures by the Embedded Atom Method (Daw and Baskes, 1984; Foiles et al., 1986; Sheng et al., 2011). **Figure 1** illustrates geometrical details of the generation process. The starting model of microstructure with screw dislocation was obtained by shifting the atomic coordinates of a perfect single crystal by *u<sup>z</sup>* = *b*θ 2π along the [*hh*0] cylinder axis; to generate edge dislocations, instead, two (110) half-planes were removed along the cylinder axis. Periodic Boundary Conditions (PBCs) for the screw dislocation line were applied along the [*hh*0] axis (**Figure 1C**), whereas they were applied both along the h *hh*2*h* i dislocation axis and along h *hh*0 i for the edge dislocation (**Figure 1A**). The initial microstructures were equilibrated for 1ns using the Langevin thermostat at 300 K combined with a constant Number of atom, Volume, and Energy (constant NVE) integration with 1fs time step. Next to the equilibration, a time trajectory was generated recording sequences of 100 microstructures at 1ps time intervals.

A time-average of the arrangement in space of the atomic positions was computed along the time trajectory, so to cancel the thermal effects out (Leonardi et al., 2011); this Time Averaged Microstructure (TAM) was then used to calculate powder diffraction patterns by the Debye scattering equation (DSE) (Debye, 1915; Gelisio et al., 2010), using the atomic coordinates in the cylindrical regions shown in **Figure 1** (black line wireframe). As in the Krivoglaz–Wilkens approach, dislocations were always straight lines running parallel to the cylinder axis. Effects related to the position of the dislocation line were also considered, randomly displacing the cylinder region of interest and then considering corresponding powder patterns and averages. For comparison, similar procedures were carried out for the same cylindrical regions without any line defect. To assess the role of domain size and surface, few other regions were also considered for different domain shapes (cube and sphere) and for smaller systems. In particular a *D* = 20 nm, *H* = 28.7 nm cylindrical domain was generated,

both as a smaller version of that shown in **Figure 1C** (*D* = 40 nm, *H* = 28.7 nm), and carved out from it.

#### **RESULTS**

**Figure 2** shows the isotropic (volumetric) and deviatoric strain fields for Pd bulk nanocrystals containing edge or screw dislocations. As expected, the edge case gives compressive/tensile (**Figure 2A**) and deviatoric (**Figure 2B**) strains, whereas the screw dislocation gives a predominantly deviatoric strain (**Figure 2D**). In both cases, the dislocation line splits in partials, according to the known reaction: <sup>1</sup> 2 h110i → <sup>1</sup> 6 h211i + 1 6 121 , a feature especially visible for the edge dislocation. The separation distance between edge partials (about 4 nm) is sensibly larger than predicted by simple considerations on surface energy of the faulted region between partials (~1 nm) (Hull and Bacon, 1965), but is well in agreement with other, more recent literature values (Hunter et al., 2011). Also the screw dislocation is not exactly as an ideal straight line, even though splitting and other deviations are less pronounced than in the edge case.

Extended dislocations have a complex effect on diffraction, well beyond a simple broadening of the diffraction line profiles. **Figure 3** shows the XRD intensity distribution on three orthogonal cross-sections of the reciprocal space (RS) for cylindrical nanocrystals containing edge (**Figure 3**, De, Ee, and Fe) or screw (**Figure 3**, D<sup>s</sup> , E<sup>s</sup> , and Fs) dislocations. It is also shown the intensity scattered by the corresponding perfect nanocrystal, i.e., the same

cylinder with no line defect and atoms in perfect, ideal positions (**Figure 3**, Ae, Be, C<sup>e</sup> and A<sup>s</sup> , B<sup>s</sup> , C<sup>s</sup> for edge and screw, respectively). As expected, the strain field gives a growing effect with the distance from the RS origin. Depending on the Miller indices the intensity distribution around points is differently affected, and in some cases splits in two distinct regions. This feature is clearly observed in the edge case, along the h *hh*0 i direction, as an effect of the compressive and tensile strains, respectively above and below the dislocation slip plane (cf. **Figures 1** and **2**). Strain and faulted region in between the two partials affect in a rather complex way the distribution of intensity around all points in RS. Although mediated by the average over different orientations of the cylindrical domains, this complexity is expected to appear also in the corresponding powder patterns, as shown further below.

Another view-point on the atomic displacement effect of extended dislocations is provided by the Pair Distribution Function, shown in **Figure 4**. The same plot for perfect cylinders (i.e., domains with no line defects) gives an array of δ functions at all atom pair distances. These infinitely narrow bars, marked at the bottom of the plot of **Figure 4**, are broadened by the strain field of the dislocation lines, thus producing a sequence of distributions. Distribution widths for edge and screw cases are comparable, although the former gives additional peaks caused by the faulted region between the dislocation partials (**Figures 1** and **2**), which is responsible for a fraction of non-cubic sequence of atomic layer stacking. Effects on the diffraction pattern from a powder of dislocated cylindrical nanocrystals are therefore expected to be quite strong and different for the edge and screw dislocation cases.

**Figure 5** shows the powder patterns simulated by the DSE from the TAM of cylindrical Pd nanocrystals (*D* = *H* = 16 nm) with edge or screw dislocations, and corresponding ideally perfect ("crystallographic") nanocrystals. Despite the different orientation of the cylindrical nanocrystals (cf. **Figure 1**), the nanocrystal shape with equal height and diameter make the crystallographic powder patterns quite similar. Visible differences are caused by the dislocations. The screw case gives a predominant effect of broadening, whereas, the edge case gives broadening and shape effects, the latter caused by the non-cubic atomic layer stacking in the faulted region between edge partials. Such effects are stronger for Pd than Ir (**Figure 6**), as the separation between the partials, and therefore the extension of the stacking fault ribbon in the cylindrical microstructure, is larger for Pd than for Ir.

**Figure 7** shows the effect of position of the dislocation line. To be compatible with assumptions underlying Eq. 1, all dislocations were straight lines parallel to the cylinder axis. The position of the dislocation line has negligible effects in the case of screw dislocations, as a consequence of the mainly deviatoric strain field introduced in the TAM. Differences are quite strong for the edge dislocation case. Features originating from the more complex strain field of the edge type, and the stacking fault region associated to the two partials are always visible. It can only be noted that a pattern obtained by averaging those for different random positions of the edge line is quite similar to the pattern with the edge dislocation along the cylinder axis. However, even such average pattern is clearly affected by strong deviations from the pattern expected for an fcc metal phase.

Molecular Dynamics simulations can easily be extended to different nanocrystal sizes and shapes. **Figure 8** shows the effect of changing the shape of the nanocrystal, following the same procedure described above for cylinders. Powder patterns from nanocrystals containing edge dislocations always show more complex details than the simple broadening provided by Eq. 1. The screw case seems qualitatively closer to the expected effects of Eq. 1, although some peaks present visible splitting [like the (642) line (Wilson, 1955) just above *Q* = 12 in **Figure 8**]. For the screw case in **Figure 9**, we also show the effect of changing the size of the cylindrical nanocrystal: largest effects are caused by changing the diameter, as can be easily explained considering that changing diameter acts on two dimensions, i.e., on the extension of the cylinder base.

#### **DISCUSSION**

The DSE powder patterns described so far can be considered as "experimental" diffraction data and analyzed by a state-of-theart method, like WPPM (Scardi and Leoni, 2002; Scardi, 2008). Besides using Eq. 1 for refining the dislocation parameters (like ρ, *Re* , and*C*hkl),WPPM can model nanocrystalline domains of virtually any size and shape (Leonardi et al., 2012), also considering the presence of stacking faults (Scardi and Leoni, 2002; Scardi, 2008) and other microstructural features responsible for line broadening effects (Scardi, 2008). For example, the complex effect of relaxation of the nanocrystal surface can be described by an additional "strain" profile component, with a corresponding FT given by

$$\begin{split} A\_{\text{\{hkl\}}}^{\text{SR}}(L) &= \exp\left[ -2\pi^2 s^2 L^2 < \mathfrak{e}\_{\text{\{hkl\}}}^2(L) > \right] \\ &= \exp\left[ -2\pi^2 s^2 L^2 \Gamma\_{\text{hkl}}\left(aL + bL^2\right) \right] \end{split} \tag{2}$$

**FIGURE 3 | X-ray diffraction intensity distribution in reciprocal space**. Three different cross-sections are shown for cylindrical Pd nanocrystals (D = H = 16 nm) with edge dislocation (De, Ee, and Fe) and screw dislocation

where Γhkl = 1 + *c h* 2 *k* <sup>2</sup> + *k* 2 *l* <sup>2</sup> + *l* 2*h* 2 / *h* <sup>2</sup> + *k* <sup>2</sup> + *l* 2 <sup>2</sup> = 1 + *c* · *H*<sup>2</sup> accounts for strain anisotropy, referred to the Miller indices of the diffraction peak. Parameters *a*, *b*, *c* can be adjusted to model the strain field and anisotropy.

The FT of the overall diffraction line profile, *A*(*L*), can then be written as a product of all required FTs (Scardi, 2008):

$$A\ (L) = A^{D}\_{\{\text{hkl}\}}\ (L) \cdot A^{S}\_{\{\text{hkl}\}}\ (L) \cdot A^{\text{SR}}\_{\{\text{hkl}\}}\ (L) \dots \tag{3}$$

where*A S* {hkl} (*L*)is the FT of the domain size effect, known in closed analytical form for a cylindrical shape (Langford and Louer, 1982). Therefore, in this specific case, besides refinable parameters of Eqs. 1 and 3 includes height (*H*) and base diameter (*D*) of the cylindrical domain, and adjustable parameters *a*, *b*,*c* of Eq. 2. If necessary, the effect of stacking faults in reasonably low concentration can be described by a *A F* hkl (*L*) component given by Warren's theory for faulted fcc metal structures (Warren, 1990; Scardi, 2008).

(D<sup>s</sup> , E<sup>s</sup> , and F<sup>s</sup> ). XRD intensity from corresponding cylinders without line defects are shown in (Ae, Be, and Ce) and (A<sup>s</sup> , B<sup>s</sup> , and C<sup>s</sup> ). Refer to **Figure 1** for details on the cylindrical domains.

Results of **Figures 5**–**8** and relevant discussion point out the difficulty in modeling the pattern of a powder of cylindrical domains with edge dislocations. The traditional approach, underlying also Eq. 3, considers a perfect cylindrical fcc domain with defects; such a perturbation approach seems little convincing here, as it cannot account in any simple and accurate way for the faulted region between the two partials, with the corresponding hexagonal sequences, and the complex strain field, strongly dependent on split of the dislocation line in partials and their position inside the domain (cf. **Figure 7**). The screw dislocation case seems more "well-behaving," i.e., closer to the assumptions of Krivoglaz–Wilkens theory, mostly involving a line broadening effect. Therefore, in the following we focus on the screw case only, analyzing the corresponding powder patterns by the WPPM approach as described by Eq. 3.

**Figure 10** shows the modeling results for two different powders of cylindrical domains, respectively *D* = 40/*H* = 28.7 nm

**FIGURE 4 | A selected range of the distribution function of atom pair distances in theTime Averaged Microstructure of cylindrical Pd nanocrystals (D** = **H** = **16 nm) with edge (red circular dot) and screw (blue square dot) dislocation lines**. The inset shows the trend across the whole range of distances.

(**Figures 10A–C**) and *D* = 20/*H* = 28.7 nm (**Figures 10D–F**). The result for ideal cylindrical domains (**Figures 10A,D**), i.e., with Pd atoms positioned according to ideal fcc structure, no dislocations, and no surface relaxation, is very good, as expected in case of domain size broadening effects only (Leonardi et al., 2012). Small deviations between DSE pattern and WPPM are expected,owing to the different hypotheses underlying DSE and WPMM, as the former is based on an intrinsically discrete, atomistic model, whereas, the last considers crystalline domains as ideal solid models, i.e.,

**FIGURE 6 | X-ray powder diffraction pattern simulated by DSE from the ("equilibrated")TAM of cylindrical nanocrystals (D** = **H** = **16 nm) with edge dislocations (line), and corresponding ideally perfect ("crystallographic") nanocrystals (dash)**. Results are shown for Pd (red circular dot) and Ir (blue square dot), having respectively larger and smaller

separation distance between partial dislocations.

cylinders with a smooth surface [details can be found in Beyerlein et al. (2011)].

The modeling is still reasonably good when the atomistic model of the same cylindrical domain is equilibrated before, generating

the pattern by the DSE (**Figures 10B,E**). Energy minimization leads to a strain broadening effect caused by the relaxation of the free surface, which adds to the size broadening effect from the finite cylindrical domain. The effect is well represented by the model of Eq. 2. The WPPM result, instead, is much less satisfactory for the cylindrical domain containing a screw dislocation along the axis, showing marked and systematic deviations from the DSE pattern (**Figures 10C,F**): peak width is reasonably matched but details of the peak profile shape are definitely not reproduced. More in particular, the model of Eq. 1 is unstable: even if the contrast factor is fixed to the expected value for a screw dislocation (for the primary slip system of Pd, {111}h110i, *C*hkl = *A* + *B* · *H*<sup>2</sup> = 0.280476 − 0.64335 · *H*<sup>2</sup> ) when ρ and *R<sup>e</sup>* are both allowed to vary, the last diverges (>1012 nm). In the case of the smaller cylinder (*D* = 20/*H* = 28.7 nm), even the cylinder height, when freely refined, tends to wrong (smaller) values. Such instability and drift of the refinement toward wrong values is partly due to the intrinsic correlation between parameters – quite clearly, between ρ and *R<sup>e</sup>* – but is also a result of the complexity of the strain field, which is not fully captured by the model of Eq. 1.

The modeling improves if more parameters are fixed. Indeed, besides the contrast factor we can fix D and H to the model values (for the specific cylinder considered), and also set the dislocation density to the expected value, given by the ratio between dislocation length and cylinder volume: ρ*<sup>t</sup>* = *H*/ - *H* · π(*D*/2 ) 2 = - π(*D*/2 ) 2 −1m−<sup>2</sup> . This gives values of 0.796 × 10−<sup>15</sup> m−<sup>2</sup> and 3.183 × 10−<sup>15</sup> m−<sup>2</sup> , respectively for *D* = 40/*H* = 28.7 nm and *D* = 20/*H* = 28.7 nm cylinders.

Then the only microstructural parameter to be refined, besides unit cell parameters, is the effective outer cut-off radius. This is

how the refinements for the larger cylinder in **Figures 10C,F** were obtained. For a more robust convergence we refined together the patterns of both cylinders, without and with screw dislocation (**Figures 10B,C**, respectively), thus refining the same values of *a*, *b*, and *c* in Eq. 3, together with *R<sup>e</sup>* for the cylinder with screw dislocation. The same procedure was repeated for the smaller cylinder (**Figures 10E,F**).

Best value of *R<sup>e</sup>* for the *D* = 40/*H* = 28.7 nm cylinder is 41.3(1) nm, quite close to twice the cylinder radius. This value is in close agreement with Wilkens model of Eq. 1, where line defects are assumed to be inside the so-called "restrictedly random dislocations regions," whose radius (*Rp*) is related to our definition of *R<sup>e</sup>* as *R<sup>e</sup>* ≈ 2*R<sup>p</sup>* = *D* [see Wilkens (1970b) for definitions, and Armstrong et al. (2006) for a more recent review on the validity of Eq. 1 and relation between *R<sup>e</sup>* and physical lengths of crystalline domains containing dislocations].

For the smaller cylinder, *D* = 20/*H* = 28.7 nm, *R<sup>e</sup>* is proportionally larger, 26.0(1) nm, but still not far from the expected value of 2*Rp*. It is therefore verified the correctness of the hypotheses underlying the Krivoglaz–Wilkens model, although Eq. 1 can be unstable when trying to refine all parameters, especially those which correlate more strongly. The agreement with the model hypotheses increases with the domain size, and we can expect the model to be exact in the limit of very large diameters. Deviations in smaller domains are partly related to the non-ideality of dislocations, but also reflect the effect of the dislocation core region, which in the Krivoglaz–Wilkens model is excluded by an inner cut-off radius (Wilkens, 1970a,b).

Strain effects on the diffraction line profile analysis can be described by an r.m.s. strain (or microstrain) plot, originally proposed by Warren and Averbach (1950), Warren (1990).

This plot, shown in **Figure 11** for both models of cylindrical domains containing a screw dislocation, provides the r.m.s. strain (hε 2 {hkl} (*L*)i 1/2 , the width of the strain distribution) over different distances *L* inside the crystalline domain, taken along the scattering vector direction; as a consequence, the microstrain depends on the hhkli crystallographic direction. **Figure 11** shows trends along h111i, h200i, and h220i for both strain components, respectively due to the dislocation (line) and to the grain boundary relaxation (dash). Strain values increase for smaller domain sizes, following the corresponding increase in dislocation density [from Eq. 1, hε 2 {hkl} (*L*)i <sup>1</sup>/<sup>2</sup> ∝ √ ρ], whereas the effect of the grain boundary strain decreases with increasing diameter.

To better assess the effect of the cylindrical grain boundary, we carved a *D* = 20 nm/*H* = 28.7 nm cylinder from the larger one containing a screw dislocation along its axis, and then generated the powder pattern by the DSE. The WPPM analysis was made considering only the strain effect from the dislocation, with same contrast factor, fixed domain size and dislocation density (3.183 × 10−<sup>15</sup> m−<sup>2</sup> ). The result gave *R<sup>e</sup>* = 26.4(1) nm, nearly the same value refined for the same cylinder with surface relaxation effect (**Figures 10E,F**), thus confirming that *R<sup>e</sup>* approaches the expected 2*R<sup>p</sup>* value if all conditions of the Wilkens model are verified. Also in this case, however, the instability of Eq. (1) is confirmed, as a free refinement of all parameters gives diverging values of *R<sup>e</sup>* , a dislocation density much smaller than expected, and a wrong height of the cylindrical domain.

It is interesting to consider the results shown so far in terms of the main strain component in the studied systems. **Figure 12** shows the deviatoric strain as a function of the radial distance from the cylinder axis, which is also the position of the screw dislocation line (cf. **Figure 1C**).

The strain for an ideal dislocation (Hull and Bacon, 1965) is shown together with the values obtained from MD, for both cylindrical domains (*D* = 20 nm/*H* = 28.7 nm and *D* = 40 nm/*H* = 28.7 nm) with and without screw dislocation.

Apart from the dislocation core region, where the continuum mechanics expression diverges while the MD values stay finite, it is quite evident that the MD trends for domains with screw dislocations result from a combination of the strain from the dislocation with that from the surface relaxation effect. The last steeply decays from the surface toward the inside the domain, and tends to decrease for increasing diameter of the cylindrical domain.

### **CONCLUSION**

Differences between real dislocations and the idealized models of continuum mechanics, become significant when dislocations are confined in nanocrystalline domains. MD simulations show that the traditional continuum mechanics expressions for the strain field only approximately agree with the actual strain field. This was demonstrated for the simple case of straight dislocations in cylindrical Pd nanocrystal, a condition closely matching the hypotheses underlying the Krivoglaz– Wilkens theory on dislocation effects in diffraction. Discrepancies between theory and MD simulations are especially evident for edge dislocations, as an effect of the split into partials and the corresponding stacking faults in between, responsible for a region of non-hexagonal layer stacking. Under these conditions, the traditional perturbation approach, based on a cubic (fcc) Pd phase with deformation caused by lattice defects, seems not appropriate to model XRD patterns obtained from MD simulations. Further discrepancy is observed if the edge dislocation line does not lay along the axis of the cylindrical domain.

The screw case seems more similar and compatible with the Krivoglaz–Wilkens theory. Diffraction patterns generated by MD seem little affected by the position of the dislocation line in the cylindrical domain, and the main effect is a broadening of the diffraction lines, as predicted by the theory. If all parameters of the system – cylindrical domain height and diameter, contrast factor of the screw dislocation in the primary slip system of Pd, dislocation density as given by line defect length divided by the cylinder volume – are fixed, then the dislocation outer cut-off radius is found in good agreement with the Wilkens model for restrictedly random dislocations: *R<sup>e</sup>* is about 10–20% larger than the diameter (2*Rp*) of the restrictedly random regions, the discrepancy decreasing for increasing cylinder diameter.

Despite this positive result, if the powder diffraction pattern generated from the MD simulation is modeled according to the Wilkens theory, the non-ideality of the dislocations and the intrinsic correlations between parameters lead to a strong instability. If all microstructural parameters are allowed to vary without constraints, the outer cut-off radius tends to diverge, leading all other parameters to wrong values. It is therefore likely that significant errors may occur in the experimental practice, when Wilkens model is applied to real materials, e.g., plastically deformed metals; applying restrictions to some parameters, possibly exploiting evidence from other techniques, might significantly help to obtain reliable results from the analysis of the diffraction patterns. Effects are expected to be increasingly significant for smaller domain sizes. Further studies will be required to shed light on this important issue, in the effort to provide a realistic modeling of polycrystalline materials with lattice defects.

#### **REFERENCES**


elastic isotropy applied to hexagonal crystals. *J. Appl. Crystallogr.* 21, 59–66. doi:10.1107/S0021889887009580


Wilson, A. J. C. (1955). The effects of dislocations on X-ray diffraction. *Il Nuovo Cimento* 1, 277–283. doi:10.1007/BF02900634

Yamakov, V., Wolf, D., Phillpot, S. R., Mukherjee, A. K., and Gleiter, H. (2002). Aluminium reveals some surprising behaviour of nanocrystalline aluminium by molecular-dynamics simulation. *Nat. Mater.* 1, 1–4. doi:10.1038/nmat700

Yamakov, V., Wolf, D., Phillpot, S. R., Mukherjee, A. K., and Gleiter, H. (2004). Deformation-mechanism map for nanocrystalline metals by moleculardynamics simulation. *Nat. Mater.* 3, 43–47. doi:10.1038/nmat1035

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 13 November 2014; accepted: 22 December 2014; published online: 03 February 2015.*

*Citation: Leonardi A and Scardi P (2015) Atomistic model of metal nanocrystals with line defects: contribution to diffraction line profile. Front. Mater. 1:37. doi: 10.3389/fmats.2014.00037*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials.*

*Copyright © 2015 Leonardi and Scardi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **Molecular dynamics study on the distributed plasticity of penta-twinned silver nanowires**

*Sangryun Lee and Seunghwa Ryu\**

*Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, Daejeon, South Korea*

#### *Edited by:*

*Federico Bosia, University of Torino, Italy*

#### *Reviewed by:*

*Pratyush Tiwary, ETH Zürich, Switzerland Shan Jiang, University of Missouri, USA Gang Zhang, Agency for Science, Technology and Research (A\*STAR), Singapore*

#### *\*Correspondence:*

*Seunghwa Ryu, Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, KAIST Building N7-4, Room 7107, 291 Daehak-ro, Yuseong-gu, Daejeon 305-701, South Korea ryush@kaist.ac.kr*

#### *Specialty section:*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials*

> *Received: 15 May 2015 Accepted: 27 July 2015 Published: 19 August 2015*

#### *Citation:*

*Lee S and Ryu S (2015) Molecular dynamics study on the distributed plasticity of penta-twinned silver nanowires. Front. Mater. 2:56. doi: 10.3389/fmats.2015.00056* The distributed plasticity of pentatwinned silver nanowires has been revealed in recent computational and experimental studies. However, the molecular dynamics (MD) simulations have not considered the imperfections seen in experiments, such as irregular surface undulations, the high aspect ratio of nanowires, and the stiffness of loading devices. In this work, we report the effect of such inherent imperfections on the distributed plasticity of penta-twinned silver nanowires in MD simulations. We find that the distributed plasticity occurs for nanowires having undulations that are less than 5% of the nanowire diameter. The elastic stress field induced by a stacking fault promotes the nucleation of successive stacking fault decahedrons (SFDs) at long distance, making it hard for necking to occur. By comparing the tensile simulation using the steered molecular dynamics method with the tensile simulation with periodic boundary condition (PBC), we show that a sufficiently long nanowire must be used in the constant strain rate simulations with PBC, because the plastic displacement burst caused by the SFD formation induces compressive stress, promoting the removal of other SFDs. Our finding can serve as a guidance for the MD simulation of crystalline materials with large plastic deformation and in the design of mechanically reliable devices based on silver nanowires.

**Keywords: molecular dynamics, dislocation, stacking fault decahedrons, distributed plasticity**

# **Introduction**

Silver nanowires have been widely used as building blocks in various electronic systems including flexible antennas (Song et al., 2014), highly stretchable strain sensors (Liu and Choi, 2009; Amjadi et al., 2014), transparent conductive films, and flexible conductors (Xu and Zhu, 2012) due to their high conductivity and ease of synthesis (Zhang et al., 2005). Most of those applications are expected to sustain cyclic mechanical loading and be used as components in wearable and flexible devices. For example, silver nanowire–elastomer composite based piezoresistivity sensors experience complex cyclic mechanical loading during their operation (Amjadi et al., 2014) and the nanowire-based touch screen has also been tested under large bending and stretching (Xu and Zhu, 2012). The conductivity of these devices is provided by a percolation network made of silver nanowires and their failure highly depends on the integrity of the network. Therefore, for the reliable design of flexible devices, it is essential to understand the mechanical integrity of the silver nanowire under high strain beyond the plastic limit as well as its recovery under unloading.

In addition, there is a fundamental interest in understanding the plasticity of silver nanowires, for two reasons. First, silver nanowires have a unique penta-twinned structure when chemically synthesized (Zhang et al., 2005; Liu and Choi, 2009; Xu and Zhu, 2012; Amjadi et al., 2014; Song et al., 2014). Nanostructures incorporating nanotwins have been extensively studied with the goal of simultaneously improving their strength and ductility. Second, it is known that the plasticity of a sub-100 nm metal specimen is driven by dislocation nucleation, and thus is very different from bulk counterparts, in which plastic deformation is governed by the motion of existing dislocations. Penta-twinned nanowires can serve as a simple test bed for studying both the role of twin boundaries (TBs) (Zhu and Gao, 2012) and the size effect, together. For example, a few novel properties of penta-twinned nanowires have been revealed including enhance elastic modulus (Zhu et al., 2012), enhanced yield strength (Wu et al., 2006), special multiple conjoint fivefold twins (Jiang et al., 2013a,b), and thermal effects originated from the size effect and the TBs (Wu et al., 2011).

Numerous computational and experimental studies have been devoted to examining the plasticity of penta-twinned nanowires (Cao and Wei, 2006; Wu et al., 2011; Filleter et al., 2012; Sun et al., 2013; Zheng et al., 2014). It was found that the yield of the nanowire initiates with dislocation nucleation from the surface, and plastic deformation is accommodated by a chain of stacking fault decahedrons (SFD) or stacking fault hats (SFH). Multiple distributed plastic zones have been observed for sub-100 nm nanowires, and their mechanism has been explained by the activation and deactivation of SFD chains at stress concentrations due to a surface undulation. Recently, a recoverable plasticity has also been reported by two independent studies (Bernal et al., 2014; Qin et al., 2015) suggesting the crucial role of TBs in blocking the propagation of stacking faults and promoting the recovery under unloading.

However, molecular dynamics (MD) simulation studies have not considered the irregular surface undulation, the high aspect ratio of nanowires, and the stiffness of loading devices that are inherent in experiments. While stress concentrations associated with the irregular surface undulation have been suggested as the origin of distributed plastic zones, there has been no MD simulation that tests the assumptions. The large plastic displacement burst induced by the formation of the SFD chain may significantly affect the stress–strain relationship if nanowires with small aspect ratio are tested in MD simulations. Also, the constant strain rate simulation with a periodic boundary condition (PBC) is equivalent to the use of an infinite stiffness loading device, which may result in a deformation mechanism that is different from the experiments.

In this work, we employ MD simulations to study the effect of the aforementioned factors in the deformation behaviors. We found that the distributed plasticity occurs for nanowires with undulations that are less than 5% of the nanowire diameter. The elastic stress field induced by a stacking fault promotes the delocalization of SFDs, making it hard for necking to occur. For higher diameter undulation, the stress concentration is too high to accommodate a separate plastic zone at a different stressconcentration position. By comparing the tensile simulation using the steered molecular dynamics (SMD) method with the tensile simulation with PBC, we show that a sufficiently long nanowire must be used in the constant strain rate simulations with PBC, because the plastic displacement burst caused by the SFD formation induces compressive stress, which promotes the removal of other SFHs or SFDs.

# **Materials and Methods**

We employed LAMMPS (Plimpton, 1995) and MD++ (MD++; http://micro.stanford.edu/MDpp) MD packages to study the mechanical behavior of penta-twinned silver nanowires. The interactions between atoms were described by embedded-atommethod (EAM) potential. We tested two versions of EAM potentials developed by Sheng et al. (2011) and Williams et al. (2006). The diameter of the nanowires was around 20 nm, and two different aspect ratios (6 and 12) were considered. The penta-twinned nanowire has a *<*110*>* direction along the growth axis and five single face-centered-cubic (FCC) crystals intersect at the {111} TBs (**Figure 1**). To mimic the surface defects and diameter variation observed in experiments (Filleter et al., 2012; Bernal et al., 2014), the nanowire has surface undulations along the axial direction (Agrawal et al., 2009) as shown in **Figure 1B**. At each cross section of the undulation site, we considered an imaginary circle of identical size but centered at a shifted position, and atoms above the circle are removed. The magnitude of shift is a sinusoidal function of the axial position with maximum being 5% of nanowire diameter. We did not observe significant re-arrangements on the undulation sites because there is no sharp surface step. The undulation sites are equally spaced in all nanowires considered in this study. The total number of atoms is about two million for the nanowire with the aspect ratio 6 and four million for the aspect ratio 12. Before the tension, the nanowires are equilibrated at 300 K by utilizing NPT ensemble for 600 ps and NVT ensemble for 600 ps until long period vibration modes decay. Time step and damping parameters are chosen as 1 fs and 0.1, respectively. The

average length of nanowire during the NPT equilibration is used as the reference length for computing engineering strain.

The tensile loading was performed with two methods. One is the conventional tensile simulation where the simulation box is expanded with a constant rate under PBC. The other is the SMD method where one end of the nanowire is pulled by an imaginary spring while the other end is fixed. For both simulations, we used the time step of 1 fs. The tensile loading simulation of the conventional method was performed at 300 K under canonical ensemble (NVT) with an engineering strain rate of 0.1 ns*−*<sup>1</sup> until fracture occurred. We increase the size of the simulation cell in the loading direction every 10 ps with a fixed strain increment (0.001) and run NVT simulations for 10 ps until it reaches equilibrium. We then measured the average stress, and repeated the elongation and NVT simulation sequence multiple times. The virial stress formula was used to calculate the atomic stress of the nanowires (Pendás, 2002; Zimmerman et al., 2004).

For the SMD method, one end of nanowire is fixed and the other end of nanowire is pulled by a virtual spring. On the onset of plasticity, a large plastic strain burst occurs that cannot be captured by the conventional PBC method. Due to the compliance of the virtual spring, the SMD method is suitable for capturing the large strain while it cannot apply a perfect displacement control. When performing SMD method, we chose the spring constant to be comparable to the stiffness of the 120 nm-long nanowires, and the pulling speed of the spring was chosen to achieve the engineering strain rate of 0.1 ns*−*<sup>1</sup> . During SMD simulation, the temperature was kept at 300 K using identical thermostats. To analyze the dislocation structure in the plastic zone, we visualized the dislocation lines with different colors according to their category using OVITO (Alexander, 2010) and Crystal Analysis (Alexander and Karsten, 2010) packages.

# **Results**

First, we compare the deformation behavior of two nanowires with and without diameter undulation (~5% of diameter), as described by the Sheng potential. For the pristine nanowire, the plastic zone is initiated at a random point and extends for a long distance by forming SFD chains, as depicted in **Figure 2A**. In contrast, the deformation of the nanowire with the undulation is accommodated by two distinct plastic zones, activated around the stress concentrated regions, which is similar to the distributed plasticity observed in experiments (Filleter et al., 2012; Bernal et al., 2014) (**Figure 2B**).

Interestingly, for the nanowire with undulations, the plastic zone extends for a large distance before necking occurs, as it does in the pristine nanowire. We found that the toughness moduli, evaluated by integrating the stress–strain curve from zero to fracture strain, are almost identical for the two cases (**Figure 2D**), which indicates that the nanowire accommodates a large plastic deformation despite the presence of stress concentration. Due to this tendency to develop a non-localized plastic zone, the pentatwinned nanowire is tolerant to small surface undulations and can absorb a large amount of energy until failure. This may explain why highly stretchable sensors or touch screens based on the nanowire percolation network can sustain significant cyclic loadings (Xu and Zhu, 2012; Amjadi et al., 2014; Narayanan et al., 2014; Yang et al., 2014; Wang et al., 2015). To analyze the dislocation structures, we drew a double Thompson tetrahedron to visualize the Burgers vector of dislocation of the FCC structure (**Figure 2C**). The most predominant dislocation is the Shockley partial with the Burgers vector 1/6*<*112*>*a (the lines painted green color in the Figure) which represents the leading and trailing part of the stacking fault planes in the SFD chains. The next dominant partial dislocation with Burgers vector 1/9*<*222*>*a (black

tetrahedron painted with different colors by categorizing Burgers vector: red,

tensile simulations of two nanowires.

color in figure) is formed at the intersection between the stacking fault and the TB (Bernal et al., 2014).

To probe whether the observed behavior was independent of the empirical potential choice as well as the small aspect ratio, we compared the tensile simulations of nanowires having four undulation sites with larger aspect ratio (12), using both Sheng and Mishin potentials (**Figures 3A–C**). The different results from the potentials are originated from different choices of fitting properties (Williams et al., 2006; Sheng et al., 2011). Still, regardless of the potential type, we observed the activation of a distinct plastic zone around the undulation sites, and found that some of SFD chains remained until the fracture took place, although the amount of plastic deformation was different between the two potential models.

Nevertheless, the MD simulations indicated that the distributed plasticity is an intrinsic property of the penta-twinned nanowire, independent of length, potentials, and the number of surface undulations. Yet, we note that the abrupt reduction of stress following the SFD chain formation has not been observed in tensile tests in experiments (Filleter et al., 2012; Bernal et al., 2014). This may be attributed to the remarkably high strain rate or the much smaller diameter used in simulation. Further research is necessary to clarify this issue.

To understand the origin of the SFD chain formation, we calculated the energy barrier for the formation of a stacking fault plane in the presence of another stacking fault by using the modified nudged elastic band (NEB) methods (**Figure 4A**) (Ryu et al., 2011). Due to high computational cost, we considered nanowires with smaller diameters (6, 8, 10, and 12 nm). For the diameter of 12 nm, the minimum occurs at a distance of around 3.5 nm (3.31 nm for Sheng potential and 3.61 nm for Mishin potential), as shown in **Figure 4C**. Also, the energy barrier is significantly smaller than the initial stacking fault nucleation energy barrier in the pristine nanowire at the same stress. This explains why, with smaller diameter, the propagation of SFH chains occurs away from the stress-concentration region. We found that the most probable successive nucleation distance δ between the first SFH and second SFH is linearly proportional to the diameter of the nanowire *D* (**Figure 4B**), with the ratio δ ~ 0.28*D*. This linear relationship implies that the observed energy minimum may be caused by an elastic stress field induced by the first SFH.

We illustrated the elastic stress field in the nanowire when a single stacking fault presents at 0 K, as shown in **Figures 4D,E**. The Von-Mises stress of atoms at the surface maximizes around 3 nm, which is close to the most probable nucleation distance δ. Von-Mises stress is calculated from six components of virial stress of each atom, with atomic volume of silver in bulk (10.3 cc/mol) (Wang et al., 1997). We also visualized the Von-Mises stress field of a single crystal nanowire under the same strain (0.04) at 0 K. As shown is **Figures 4D,E**, the stress field induced from a stacking fault in a single crystal nanowire does not have stress concentration, which clearly shows the role of the TBs in the distributed

plasticity of the penta-twinned nanowire. The mismatch between the energy barrier calculation and surface Von-Mises stress plot occurs because the nucleation barrier depends on the integral of the stress field within the dislocation loop at the transition state, not solely on the surface stress. Still, we can conclude that the elastic stress field induced by the stacking fault plane in a penta-twinned nanowire reduces the nucleation barrier of another

stacking fault plane, promoting the formation of a long SFH chain,

the first SFH with Sheng potentials. The diameter of the nanowires is 12 nm and

much like a catalytic reaction. Next, we studied the effect of nanowire length on the deformation behavior. Although the nanowire tested in these experiments was about a few micrometers long, most MD simulations consider a nanowire a few tens of nanometers scale, because of the computational cost (Cao and Wei, 2006; Agrawal et al., 2009; Zheng et al., 2014) with PBC. Nevertheless, that is not enough to fully mimic in simulations the actual experimental conditions. Especially for penta-twinned nanowires, the displacement from a large plastic strain burst is not negligible compared to the nanowire length. For conventional tensile simulations with strain control using PBC, the plastic strain burst in the nanowire induces large stress reduction that effectively causes a significant drop in elastic stress. To overcome the artifacts from the plastic strain burst, we performed the tensile simulation using the SMD method.

of AgNWs as a function of the distance.

We noticed a few differences in deformation behaviors between the two tensile testing methods. First, the extent of the SFH chain after the fracture is much larger in the case of the SMD simulation (**Figures 5A,B**). The large plastic strain burst from SFH chain formation is accommodated well in the SMD method by length extension. However, the strain burst occurs within a fixed nanowire length in the PBC method and gives rise to internal compressive stress within the nanowire, promoting the removal of SFHs.

The stress–strain curves manifest the differences between the two methods (**Figure 5C**). A much sharper stress drop is observed after the initial SFH chain nucleation in PBC due to the reduction of elastic tensile stress caused by the plastic strain burst. In contrast, **Figure 5D** shows that the nanowire length increases abruptly on the verge of plastic deformation under tensile loading by the SMD method. When the nanowires are stretched to the failure strain, the size of the plastic zone remains smaller for the PBC method because a large portion of SFH are removed by the internal compressive stress.

# **Conclusion**

In summary, we studied the distributed plasticity of pentatwinned silver nanowires by performing large scale MD simulations. We observed distributed plasticity when the diameter undulation was *<*5% for two different empirical potentials. A multiple plastic zone formation was explained by dislocation nucleation from local stress concentrations in the undulated region, which leads to SFH chain propagation. The catalytic formation of the SFH chain is caused by the reduction of the nucleation energy barrier due to the concentrated elastic stress field.

The effects of nanowire length and the loading device stiffness were investigated by comparing two tensile testing methods.

# **References**


We found that the plastic strain burst by the SFH chain induces a moderate compressive stress in the simulation cell, causing excessive removal of SFHs. Such an artifact of simulation can be overcome by employing the SMD method, in which the nanowire can accommodate large plastic strain burst by length extension.

# **Acknowledgments**

This work is supported by Basic Science Research Program (2013R1A1A010091) and Global Frontier Program (No. 2013M3A6B1078881) through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT and Future Planning.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationship that could be construed as a potential conflict of interest.

*Copyright © 2015 Lee and Ryu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Frontiers in modeling and design of bio-inspired armors

# **Stefano Signetti <sup>1</sup> and Nicola Maria Pugno1,2,3\***

<sup>1</sup> Laboratory of Bio-Inspired and Graphene Nanomechanics, Department of Civil, Environmental and Mechanical Engineering, University of Trento, Trento, Italy

<sup>2</sup> Centre for Materials and Microsystems, Fondazione Bruno Kessler, Povo (Trento), Italy

<sup>3</sup> School of Engineering and Materials Science, Queen Mary University of London, London, UK

#### \*Correspondence: nicola.pugno@unitn.it **Edited by:**

Giuseppe Saccomandi, Università di Perugia, Italy

#### **Reviewed by:**

Fernando Fraternali, University of Salerno, Italy

**Keywords: bio-inspired armors, impact testing, numerical simulation, theoretical modeling, optimization, toughness**

Impact behavior of materials and structures is of crucial interest. Indeed, every solid may experience collisions in its "mechanical" life. The topic involves a wide range of engineering applications, even in our everyday life: sports protectors [e.g., helmets, Milne et al. (2014)], sensitive portable electronics (Tempelman et al., 2012), bulletproof body armors (Cuniff, 1999), protection systems for buildings and machineries in the civil or defense sectors (NIST, 2005), improvement of crashworthiness in automotive (Schweizerhof et al., 1992), protection of spacecraft and satellite structures from high-velocity micrometeorite or orbital debris impact (NASA, 2015) are some of the most representative.

Parallel to the primary requirement of an effective protection, straightforwardly achievable with a massive armor, research efforts are aimed at weight saving due to essential and binding needs, such us better ergonomics and flexibility (body armors), transportability (vehicles), and in general a more judicious use of materials. Thus, the real goal is the high specific toughness. For some decades, the answer to these tasks has been the adoption of multilayer of textile and composite materials, (Abrate, 1998; Hoog, 2006) based on synthetic fibers (e.g., Kevlar®, Dyneema®) that have allowed to reach protection levels previously unimaginable with metallic targets. Nowadays, in the era of nanomaterials, we are raising the bar to atomistic 2D materials, like graphene, coupling high resistance (Lee et al., 2008) and flaw tolerance (Zhang et al., 2012) at the nanoscale, even for possible application to nanoarmours (Pugno et al., 2007; Lee et al., 2012, 2014). Alternatively, the same goal may be pursued through smart structural solutions to be employed even with traditional

materials, with all the benefits that this option implies. Nature, having worked over the ages for optimizing defense mechanisms against predators attacks or shock loads, is one of the most inspiring sources: as most remarkable examples we mention the coupled hard-soft layers in the *Arapaima gigas* fish's dermis (Yang et al., 2014), the internal undulated walls of the Bombardier Beetle's (*Carabidae*, *Brachinus*) explosion chamber (Lai and Ortiz, 2010), the cross-scale toughening mechanisms in the foam-like structure of dropping fruits (Thielen et al., 2013), dermal armors with scales (Ghosh et al., 2014), and the extreme robustness provided by the spider silk constitutive law (Cranford et al., 2012). On the other side, we could be interested in gaining an efficient strike, like the deadly underwater punch mechanism of the mantis-shrimp (Patek et al., 2005).

Upon impact, several complex physical phenomena take place: elastic–plastic deformation and wave propagation, fracture and fragmentation, heat generation (by yielding and friction), changing of material properties due to strain-rate effects up to phase change. Their occurrence and magnitude depend on the impact velocity that may be very low or up to extreme values (>3 km/s for hypervelocity impact), with increasing challenges for armor resistance as well as for its accurate modeling. The theoretical description of the basic aspects of impact mechanics (Stronge, 2000; Goldsmith, 1999, 2001) has reached a level of advanced maturity but it is in a sort of stalemate due to the severe mathematical complexity in representing the above mentioned phenomena, which also mutually interact. With high speed calculators and the development of computational methods (e.g., finite element method, FEM), simulation has become the favorite design tool, allowing optimization studies. Nonetheless, the advent of nanomaterials and bio-inspiration is further questioning the capabilities of these tools and pushing modeling research.

The traditional stand-alone experimental approach for armor design according to the philosophy "add material until it stops" it is not viable any more. First, multilayer panels can show crosscurrent behavior in relation to material coupling and interface strength, being even non-optimized for increasing areal density (Signetti and Pugno, 2014, **Figure 1**). Technological and economic limits in large scale production of nanomaterials, the difficulties in their manipulation or in their structural arrangement into complex bio-inspired structures require a systematic and reliable design process able to provide a tentative target optimum. With mere experiments is nearly impossible to investigate the whole design space for understanding still unexplained mechanism in order to mimic nature and, why not, do even better.

Going down to nanoarmours at atomic scale, we enter in a new world with completely unexplored scenarios. Deformation, fracture, contact forces are matter of potentials, electronic interactions, affinity, and reactivity of particles of colliding bodies in relation also to their atomic arrangement. Can we still call it only impact *mechanics*? Probably not. Some studies have been published about the protection capabilities of graphene nanoarmours via molecular dynamics [see Ozden et al. (2014) and Shang and Wang (2014)]: many considerations can be done with these simulations, which allow the modeling of systems up to 1 M atoms. However, the uncertainty and insensitivity of the behavior at this

scale would suggest the use *ab initio* methods, inspite of their limitation to a few thousands atoms (e.g., modeling of fewlayer graphite): epitaxy simulations (Verucchi et al., 2012; Taioli et al., 2013) represent a starting knowledge base as well as a collateral application.

Moving up the mesoscale level, one of the main challenges is the modeling of multiple crack nucleation and propagation. Presently, merging more than 3 levels of hierarchy is computationally unfeasible. Some methods have been developed to overcome the problems of the FEM method (erosion mesh-sensitive approach) even if each of them shows known limitations, like the cohesive zone elements (mesh sensitiveness, remeshing required) and the extended finite element method (not applicable for multiple crack interaction). Silling (2000) has proposed a nonlocal reformulation of the standard continuum theory of solid mechanics, called *peridynamics*. Unlike the partial differential equations of the standard theory, the integral equations of peridynamics are

applicable even when cracks and other singularities appear in the deformation field. Thus, continuous and discontinuous media can be modeled with a single set of equations. This theory naturally yields into a meshless method (Silling and Askari, 2005): this has been implemented, for instance, in the acknowledged molecular dynamics code LAMMPS (Sandia, 2015). We believe that this is a promising step toward a real multiscale approach (atomistic to continuum or within the continuum) in the same simulation.

Another advance we believe to be very interesting in the field is the *isogeometric* formulation (Hughes et al., 2005; Cottrell et al., 2009; Temizer et al., 2011), that offers the possibility of integrating finite element analysis with conventional Non-Uniform Rational B-Splinesbased CAD design tools. NURBS work both as geometry descriptors and element basis functions, following the same philosophy of isoparametric elements. Using NURBS it is easy to construct surfaces with C<sup>1</sup> or higher order of continuity and

compared to C<sup>0</sup> finite element geometry. Thus, isogeometric analysis is especially attractive to contact analysis of complex geometries (e.g., undulation, even hierarchical). This should not be limited just to impacts on complex surfaces, let us also think of modeling any problem in which adaptive and tunable wrinkling, e.g., fewlayer graphene as emblematic case (Zang et al., 2013), is exploited to provide superhydrophobicity and self-cleaning properties to surfaces. The integration with CAD would let models to be designed, tested, and adjusted on the go, with a relevant gain in time.

Being still in a relatively primordial phase in the development of these methods it is a gamble to forecast a breakthrough in simulation of bio-inspired and hierarchical nanomaterials for armors; however, it is worth keeping an eye on them since we believe them to be very promising. For sure, a synergistic combination of different and complementary research tools and multidisciplinary expertise (materials science, solid and fluid mechanics, physics) will be

essential to lead in the next years to predictive models and optimization tools. It will be the task of simulation to support good ideas, even futuristic, pushing technology to actually switch ideas to tangible innovation for a new generation of advanced bioinspired (nano)armors with significantly improved specific penetration resistance and energy absorption capability.

#### **ACKNOWLEDGMENTS**

NMP is supported by the European Research Council (ERC StG Ideas 2011 BIHSNAM no. 279985 on "Bio-Inspired Hierarchical Super-Nanomaterials," ERC PoC 2013-I REPLICA2 no. 619448 on "Large area replication of biological antiadhesive nanosurfaces," ERC PoC 2013-II KNOTOUGH no. 632277 on "Supertough knotted fibers"), by the European Commission under the Graphene FET Flagship (WP10 "Nanocomposites," no. 604391) and by the Provincia Autonoma di Trento ("Graphene Nanocomposites," no. S116/2012-242637 and reg. delib. no. 2266). SS acknowledges support from BIH-SNAM. The authors thank Ettore Barbieri for discussion and Ludovic Taxis for the English review.

#### **REFERENCES**


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 02 February 2015; accepted: 17 February 2015; published online: 28 April 2015.*

*Citation: Signetti S and Pugno NM (2015) Frontiers in modeling and design of bio-inspired armors. Front. Mater. 2:17. doi: 10.3389/fmats.2015.00017*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials.*

*Copyright © 2015 Signetti and Pugno. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **A hierarchical lattice spring model to simulate the mechanics of 2-D materials-based composites**

*Lucas Brely <sup>1</sup> , Federico Bosia<sup>1</sup> and Nicola M. Pugno2,3,4 \**

*<sup>1</sup> Department of Physics and "Nanostructured Interfaces and Surfaces" Interdepartmental Centre, Università di Torino, Torino, Italy, <sup>2</sup> Laboratory of Bio-Inspired and Graphene Nanomechanics, Department of Civil, Environmental and Mechanical Engineering, Università di Trento, Trento, Italy, <sup>3</sup> Center for Materials and Microsystems, Fondazione Bruno Kessler, Trento, Italy, <sup>4</sup> School of Engineering and Materials Science, Queen Mary University of London, London, UK*

It is known that structural biological materials such as bone or dentin show unprecedented damage tolerance, toughness, and strength. The common feature of these materials is their hierarchical heterogeneous structure, which contributes to increased energy dissipation before failure occurring at different scale levels. These structural properties are the key to achieve superior nanocomposites. Here, we develop a numerical model in order to simulate the mechanisms involved in damage progression and energy dissipation at different size scales in composites, which depend both on the heterogeneity of the material (defects or reinforcements) and on the type of hierarchical structure. Both these aspects have been incorporated into a 2-D model based on a Hierarchical Lattice Spring Model approach, accounting for geometrical non-linearities and including statistically based fracture phenomena. The model has been validated by comparing numerical results to linear elastic fracture mechanics predictions as well as to finite elements simulations, and then employed to study how hierarchical structural aspects influence composite material (mainly 2d, e.g., graphene based) properties, which is the main novel feature of the approach. Results obtained with the numerical code highlight the dependence of stress distributions (and therefore crack propagation) on matrix properties and reinforcement dispersion, geometry, and properties, and how the redistribution of stresses after the failure of sacrificial elements is directly involved in the damage tolerance of the material.

**Keywords: hierarchical lattice spring model, numerical modeling, fracture mechanics, composite materials, hierarchy**

# **Introduction**

Biological materials often display mechanical properties that differ from traditional engineering materials in that they are capable of simultaneously optimizing competing properties, such as stiffness and density or strength and toughness (Gao et al., 2003; Meyers et al., 2008; Giesa et al., 2011; Ritchie, 2011). The optimization mechanisms found in biomaterials can usually be traced back to their internal structure, which includes various characteristic features, principally heterogeneity and a hierarchical arrangement of microstructural and base components (Zhang et al., 2011; Bosia et al., 2012; Meyers et al., 2013). The challenge in recent years has therefore been to fully understand the mechanisms responsible for such outstanding properties and to replicate them in synthetic materials (Munch et al., 2008; Pugno et al., 2012; Wegst et al., 2015). Composite materials already base their lightweight and directional strengthening properties to the combination of materials with

#### *Edited by:*

*Giuseppe Saccomandi, Università di Perugia, Italy*

#### *Reviewed by:*

*Vincent Legrand, Institut de Recherche en Génie Civil et Mécanique, France Anna Marina Pandolfi, Politecnico di Milano, Italy*

#### *\*Correspondence:*

*Nicola M. Pugno, Dipartimento di Ingegneria Civile, Ambientale e Meccanica, Università di Trento, via Mesiano, 77, 38123 Trento, Italy nicola.pugno@unitn.it*

#### *Specialty section:*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials*

> *Received: 28 February 2015 Accepted: 11 June 2015 Published: 06 July 2015*

#### *Citation:*

*Brely L, Bosia F and Pugno NM (2015) A hierarchical lattice spring model to simulate the mechanics of 2-D materials-based composites. Front. Mater. 2:51. doi: 10.3389/fmats.2015.00051* considerably different properties, but have yet to achieve the simultaneous strength/toughness or stiffness/density combinations obtained in biocomposites (Ashby et al., 1995). New possibilities for super-composites have emerged with the recent introduction of micro and nano reinforcements like carbon nanotube (CNT) or graphene reinforcements (Coleman et al., 2006; Pugno, 2006; Stankovich et al., 2006; Young et al., 2012). The challenge to further develop new materials with radically improved properties is thus to apply biomimetic strategies to synthetic composite materials (Beese et al., 2014) and where possible to implement complex hierarchical structures (Dimas et al., 2013). To do this, numerical approaches are essential, since there are limited possibilities in experimentally exploring different geometries and structures. Therefore, the formulation of reliable numerical models becomes critical. For models to be predictive, they must be able to capture the main relevant aspects, i.e., heterogeneity, complex geometry, scaling, stress concentrations, damage nucleation and evolution, etc.

In the past, we have developed a so-called hierarchical fiber bundle model (HFBM) with these features, with the main aim of modeling fibrous, essentially 1D materials (Pugno et al., 2008). This code was successfully used to model nanocomposites (Bosia et al., 2010), hierarchical organization (Pugno et al., 2012), and fracture and fatigue of self-healing materials (Bosia et al., 2014, 2015). However, this code is insufficient to model more complex composite geometries where shear effects are not negligible, and a generalization to 2-D or 3-D is necessary. Additionally, a 2-D model would be well suited to the simulation of the mechanical behavior of emerging 2-D material systems such as graphene, e.g., in the evaluation of the influence of the presence and type of defects in crystal structure on the overall mechanical properties, which should be considerable given the crystal bi-dimensional structure (Banhart et al., 2010; Lopez-Polin et al., 2015). A 2-D or 3-D formulation can be achieved by adopting the so-called lattice spring model (LSM) approaches, which have been proposed in the past (Buxton and Balazs, 2002; Alava et al., 2006) and are referred to in various manners in the literature, e.g., "Bonded Particle Model" (Potyondy and Cundall, 2004), "Spring Network Model" (Curtin and Scher, 1990), "Random Spring model" (Nukala et al., 2005), and others. Brittle materials with random failure processes can be modeled by introducing random spatial distributions of springs (Beale and Srolovitz, 1988) and/or fixed or random fracture thresholds for the springs (Alava et al., 2006; Zhao et al., 2011). In order to be used for realistic simulations, spring network models need to be verified for consistency with continuum mechanics through homogenization procedures and/or need to be mapped into standard finite elements (Absi and Prager, 1975; Gusev, 2004). More complex versions of the model exploring complicated lattice energy landscapes (Puglisi and Truskinovsky, 2000) or large deformations (Friesecke and Theil, 2002) have also been investigated. A comprehensive review of Lattice models in micromechanics is given in Ostoja-Starzewski (2002).

Here, we develop and validate a 2-D hierarchical lattice spring model (HLSM), implementing the method in a multiscale scheme, and present preliminary results relative to defective reinforced

nanocomposites, elucidating specific mechanisms of damage progressions and crack propagation.

# **Materials and Methods**

To simulate the behavior of heterogeneous and hierarchical complex structures, we adopt and extend a LSM approach based on a 2-D cubic lattice interacting via harmonic springs between nearest and diagonal neighbors, as in Friesecke and Theil (2002). The adopted force–displacement relationship considered here is linear for the sake of simplicity, as discussed below (but this hypothesis can be relaxed). The regular grid consisting of nodes and springs used to discretize the 2-D material portion is shown in **Figure 1**. As mentioned, spring properties need to be assigned appropriately in order to obtain the equivalence of the strain energy of the elementary cell UCell with that of a continuum UContinuum (Absi and Prager, 1975), i.e., UCell = UContinuum. For 2D plane stress problems, we have:

$$\mathbf{U}\_{\text{Continuous}} = \mathbf{V}^\* \left[ \frac{\mathbf{v} \mathbf{E}^\*}{2(1 - \mathbf{v}^2)} (\mathbf{e}\_{\text{xx}} + \mathbf{e}\_{\text{yy}})^2 + \mathbf{G}^\* \left( \mathbf{e}\_{\text{xx}}^{\*2} + \mathbf{e}\_{\text{yy}}^{\*2} \right) \right] \tag{1}$$

$$+ 2\mathbf{G}^\* \mathbf{e}\_{\text{xy}}^{\*2} \tag{1}$$

where V\* is the volume of the elementary cell, E\* its Young's modulus, G\* its shear modulus, *v* its Poisson's coefficient, and εxx, εyy, εxy, are the components of the strain tensor. Using the relation G \* = E \* */*2(1 + ν), Eq. (1) becomes:

$$\mathbf{U\_{Contimum}} = \frac{\mathbf{V^\ast E^\ast}}{2(1+\mathbf{v})} \left[ \frac{\mathbf{v}}{1-\mathbf{v}} (\mathbf{e\_{xx}} + \mathbf{e\_{yy}})^2 + \mathbf{e\_{xx}}^2 + \mathbf{e\_{yy}}^2 + 2\mathbf{e\_{xy}}^2 \right] \tag{2}$$

The stored energy of the unit cell of the considered spring network UCell is obtained by taking the sum of the strain contributions of each spring:

$$\begin{aligned} \mathbf{U\_{Cell}} &= \sum\_{i=1}^{6} \mathbf{U\_{i}} = \frac{1}{2} \text{EV} \left( 2 \mathbf{e\_{xx}}^{2} \right) + \frac{1}{2} \text{EV} \left( 2 \mathbf{e\_{yy}}^{2} \right) \\ &+ \frac{1}{2} \text{EV}' \left( \frac{1}{2} \mathbf{e\_{xx}} + \frac{1}{2} \mathbf{e\_{yy}} + \mathbf{e\_{xy}} \right)^{2} \\ &+ \frac{1}{2} \text{EV}' \left( \frac{1}{2} \mathbf{e\_{xx}} + \frac{1}{2} \mathbf{e\_{yy}} - \mathbf{e\_{xy}} \right)^{2} \end{aligned} \tag{3}$$

where E is Young's modulus of the springs (taking for simplicity springs with the same stiffness), V the volume of the springs along the main axes and V' the volume of the diagonal springs. Equating (2) and (3), and considering for simplicity E = E\* (but the analysis can be generalized to heterogeneous materials), we have:

$$\begin{cases} \mathbf{v} = \frac{1}{3} \\ \mathbf{V}' = 2\mathbf{V} \\ \mathbf{V} = \frac{3}{8}\mathbf{V}' \end{cases} \tag{4}$$

These conditions are therefore enforced for the spring network. The condition on Poisson's ratio of the material is not excessively limiting, since in this study we are interested in discussing general concepts and highlighting qualitative effects, rather than studying specific materials.

As mentioned, linear elastic springs are considered. However, in the considered simulations, non-negligible rotations and/or translations deriving from fracture growth in the modeled material elicit the need for an iterative scheme to avoid large errors in the solution that could be committed in the case of a simple linear analysis. This is due to the occurrence of geometrical non-linearities, due to large displacements and deformations. Therefore, a total Lagrangian formulation is employed, where the governing equations are derived from the Green strain measure (Bathe and Bathe, 1996).

Assuming a quasi-static load case, the equilibrium is based on an iterative procedure to impose the balance between internal (*f* int) and external (*f* ext) forces using a Newton–Raphson iterative scheme:

*f* int *<sup>−</sup> <sup>f</sup>*

i=1

with

*f*

$$f\_{\rm int}^{\rm int} = \sum\_{}^{\rm n\_{\rm i}} f\_{\rm i}^{\rm int} = \sum\_{}^{\rm n\_{\rm i}} \int \langle \mathbf{B}\_{0}^{\rm T} \mathbf{F} \mathbf{E} \ \mathbf{e} \rangle\_{\rm i} \, \mathbf{d} \mathbf{V}\_{\rm i} \tag{6}$$

ext = 0 (5)

where n<sup>s</sup> is the number of springs in the domain, B<sup>0</sup> the transformation matrix, V<sup>i</sup> the spring volume, F the deformation gradient matrix, E the spring's Young's modulus and ε the spring's strain tensor. Considering the discrete governing equation:

i=1

Vi

$$\begin{cases} \mathbf{u}^{\text{h}} \left( \mathbf{X} \right) = \sum\_{\stackrel{\text{h}}{\text{i}}}^{\text{n}\_{\text{n}}} \mathbf{N}\_{\text{i}} \left( \mathbf{X} \right) \mathbf{u}\_{\text{i}} = \mathbf{N} \mathbf{u} \\ \quad \text{e} \left( \mathbf{X} \right) = \sum\_{\stackrel{\text{h}}{\text{i}}}^{\text{n}\_{\text{n}}} \mathbf{N}\_{\text{i},\text{X}} \mathbf{u}\_{\text{i}} = \mathbf{B}\_{0} \mathbf{u} \\ \quad \text{F} \left( \mathbf{X} \right) = \sum\_{\stackrel{\text{h}}{\text{i}}}^{\text{n}\_{\text{n}}} \mathbf{N}\_{\text{i},\text{X}} \mathbf{x}\_{\text{i}} = \mathbf{B}\_{0} \mathbf{x} \end{cases} \tag{7}$$

where n<sup>n</sup> is the number of nodes used for the discretization of the domain, u<sup>h</sup> the trial solution, N a set of linear trial functions, X the domain undeformed configuration, *x* the domain deformed configuration, and u the displacement.

Simulations are carried out in displacement control, as proposed in Batoz and Dhatt (1979). Fracture is implemented based on a failure strain criterion on the single springs. At each load step, a spring is removed from the lattice if its maximum strain ε<sup>i</sup>max is exceeded, simulating local failure:

$$f\_{\mathbf{i}}^{\mathrm{int}} = \begin{cases} \int (\mathbf{B}\_0^T \mathbf{F} \mathbf{E} \mathbf{e})\_\mathbf{i} \mathbf{d} \mathbf{V}\_{\mathbf{i}} & \mathbf{e}\_\mathbf{i} < \mathbf{e}\_{\mathbf{i}\_{\mathrm{max}}}\\ \mathbf{0} & \mathbf{e}\_\mathbf{i} > \mathbf{e}\_{\mathbf{i}\_{\mathrm{max}}} \end{cases} \tag{8}$$

At the smallest considered scale in simulations, ultimate strain values ε<sup>i</sup>max can be assigned deterministically or according to Weibull statistical distributions (Weibull, 1952) to account for material heterogeneity and defectivity occurring at a lower size scale than that considered in the simulation. The cumulative distribution function of a Weibull distribution is given by:

$$\mathcal{W}\left(\mathbf{e}\_{i\_{\text{max}}};\lambda;\mathbf{k}\right) = 1 - \mathbf{e} \begin{pmatrix} \frac{\mathcal{E}i\_{\text{max}}}{\lambda} \end{pmatrix}^{k} \tag{9}$$

Where k is the shape parameter and λ the scale parameter of the distribution. For example, if the lowest simulation scale is of the order of various nanometers in a CNT-reinforced composite, a Weibull distribution to describe the dispersion in CNT ultimate tensile strain values, accounting for their defectivity, can be estimated from experimental or numerical data from the literature.

To model heterogeneous materials, e.g., composites, different mechanical properties are assigned to the springs (e.g., elastic modulus, ultimate strain) as a function of the material type (typically matrix or reinforcement). Since we are primarily interested in evaluating the role on strength of hierarchical and structural effects, brittle fracture is modeled with a linear elastic behavior up to failure. This also allows a reduction in simulation times with respect to a plastic constitutive law; moreover, plasticity at the crack tip is accounted for thanks to the discrete nature of the model.

As discussed in the introduction, structural biomaterials present a hierarchical structure that often differs from one scale level to another. To model these structures extensively using the LSM approach above would lead to huge computational costs, and therefore a multiscale approach is adopted, i.e., subsystems at different size scales in the considered material are modeled. The mechanical properties of a domain obtained as output from a certain hierarchical level become the input for the hierarchical level above. This is illustrated schematically in **Figure 2**. The presence of both defects and reinforcements can be modeled at each hierarchical level, which is simulated separately using representative portions of the domain, starting with the smallest scale or level *n*, supposing that the structure can be modeled with *n* hierarchical levels. Mechanical properties (i.e., Weibull distribution parameters) at level *n*-1 are derived from simulations at level *n* and so on. Given the hierarchical nature of the model, we refer to it as "HLSM." Although in principle the procedure can be extended to a high number of hierarchical levels, and three are shown in **Figure 2**, for simplicity in the following only two levels of hierarchy will be considered.

# **Results**

# **Model Validation**

The main mechanisms the computational model needs to correctly reproduce at microscale are strain concentrations (e.g., at crack tips), load transfer between the matrix and the reinforcement, fracture nucleation, etc. These aspects are addressed in preliminary validation simulations, where HLSM results are compared to finite element method (FEM) calculations. For the latter, the COMSOL Multiphysics commercial package is used, specifically the Structural Mechanics Module. To validate the developed HLSM code, we consider the strain distributions arising in (a) a laterally cracked representative volume element (RVE) and (b) a RVE with a centrally located reinforcement.

## Laterally Cracked RVE

We first consider a pre-cracked RVE subjected to a tensile strain ε<sup>0</sup> = 5%, perpendicular to the direction of the crack (*X* direction), as shown in the schematic in **Figure 3**. The RVE is constituted by a 200 *×* 150 regular node grid in the HLSM and the initial crack length l extends to one-third of the RVE width w. A relatively soft material, simulating a polymeric matrix, is considered, with Young's modulus E<sup>m</sup> = 1 GPa and Poisson's ratio ν<sup>m</sup> = 0.3. The same geometry and boundary conditions are considered in FEM simulations, using a 2-D plane strain approximation, a free triangular element mesh refined in the vicinity of the crack tip, and second-order Lagrange elements. HLSM and FEM results are compared in the plot in **Figure 3** for the vertical strain *YY* taken along the line of the crack, in the *X* direction, with the origin taken as the crack tip. The agreement is very good, showing that the adopted HLSM is capable of simulating correct strain fields in the vicinity of singularities such as crack tips. The small discrepancies can be physically due to the discrete nature of the material and thus can be reduced adopting a finer discretization in the HLSM model. In the plot, we also report the predicted strain profile using quantized fracture mechanics (QFM) (Pugno and Ruoff, 2004),

for which in proximity of the crack tip the strain ε varies as:

$$\mathbf{e}(\mathbf{X}) = \mathbf{e}\_0 \sqrt{\frac{\mathbf{a}}{\mathbf{a} + \mathbf{X}}} \tag{10}$$

where *a* is the fracture quantum, taken as half the discretization length.

Also shown in **Figure 3** is the overall spatial distribution of strain YY in the pre-cracked RVE. The calculated field, through HLSM and FEM, is consistent with results in the literature (Irwin, 1957) and the associated graph shows the strain concentration near and around the crack tip.

### Centrally Reinforced RVE

Reinforcements embedded in a softer matrix behave as stiffeners and take up the strains applied on the surrounding matrix zone. The strain is transferred through the interface between the two materials, which is assumed to be perfect in the present case. The considered RVE is as previously discretized in a 200 by 150 regular node grid, and schematically shown in **Figure 4**. The reinforcement is chosen as linear elastic, with a stiffness 100 times that of the matrix and a high aspect ratio of 80, to model nanoreinforcements such as CNTs or Graphene nanoplatelets, which can achieve such high slenderness ratios. The applied displacement corresponds to a 5% strain in the length direction of the inclusion (*Y* axis). The calculated strain distributions along the directions parallel and perpendicular to the reinforcement (*Y* and *X* axes, respectively) are shown in the plots in **Figure 4**, where HLSM results are compared to FEM calculations. The agreement is again very good. The corresponding full-field spatial distribution of strain in the *Y*-direction around the inclusion is shown in **Figure 5**. It can be seen that deformations inside the matrix decrease around the linear inclusion inside the "reinforced area." The level of load transfer from the matrix to the reinforcement is indicative of the efficiency of the reinforced structure. On the other hand, deformations increase around the tips of the reinforcement, and these regions are therefore typically where crack initiate (as can be observed in the inset of **Figure 5**) and then continue to propagate. These results agree with well-known "shear

**FIGURE 4 | HLSM simulation results for a tensile loading experiment on a centrally reinforced RVE (schematic on the left)**. Strain YY variation along *X* and *Y* axes normalized with respect to width w and length L (bottom right and top right, respectively) compared to that of FEM simulations.

lag" models (Gao and Li, 2005), and also with experimental data in the literature (Bigoni et al., 2008), and also show the interface debonding that can occur near the crack tip. Consequently, in the case of high aspect ratio inclusions, the presence of reinforcements in the matrix provide a stiffening of the material, but also anticipate fracture nucleation and growth due to the stress concentration at the tips. This illustrates one of the practical aspects of the conflict between strength and toughness in the design of composite structured materials.

# **Simulations**

Simulations are first carried out on a model reinforced composite structure to highlight the mechanisms responsible for fracture nucleation and propagation. The considered reinforcements are again much stiffer than the matrix (E<sup>f</sup> = 100Em), have a high aspect ratio (l/w = 50), and are distributed homogeneously in a staggered and aligned arrangement parallel to the loading direction, with a volume fraction corresponding to approximately V<sup>f</sup> = 6.7%. **Figure 6** illustrates the evolution of the simulated damage and crack propagation in such a composite RVE when subjected to tensile loading. Again, strains in the *Y* direction are represented in color scale. The cracks initiate in regions adjacent to the tips of the reinforcement where initial strain concentrations are maximum (**Figure 6A**), propagate horizontally, and are stopped when they reach the closest reinforcement (**Figure 6B**). The load is transferred to the reinforcements through the shear deformation of the matrix, as discussed in Section "Results." At this point, fracture propagates in the direction of the fibers through shear loading, giving matrix-fiber debonding and thus connecting different cracks originated at the tips (**Figure 6C**). This leads to the ultimate failure of the composite (**Figure 6D**). The corresponding stress–strain curves (not shown in this case) display a non-linear behavior which can be assimilated to elasto-plastic

the plastic part. **(C)** debonding leading to the material failure **(D)**.

behavior deriving from crack nucleation and stopping. Thus, the non-linearity is only structure related, due to the linear elastic behavior assumed here for the material constituents. Similar results have been obtained considering various types of bi-material heterogeneous structures, in bone- or nacre-like ("brick/mortar") configurations (Sen and Buehler, 2011).

Next, the influence of porosity at microscopic level has been evaluated by including varying defect percentages in the simulations, in a multi-scale hierarchical approach. The lower hierarchical scale is used to model a matrix containing different volume fractions of uniformly randomly distributed pores or defects, represented here by a percentage of randomly assigned broken springs. Numerical uniaxial tensile tests are performed on these sub-volumes and the corresponding mechanical properties (stiffness and ultimate strain) are obtained. The resulting properties display some statistical dispersion, given the random distribution of the defects which can differ from one simulation to the next, and are thus fitted using a Weibull distribution, which is widely used in modeling the strength of solids (Weibull, 1952). The scale and shape parameters for the derived stiffness and ultimate strain Weibull distributions then become the input for the single springs in the upper hierarchical level, which includes both matrix and reinforcements. The springs in the reinforcements are chosen as previously 100 times stiffer and stronger than those used for the matrix, and reinforcement dimensions and distributions are the same as for those of simulations in **Figure 6**. Results for the upper

hierarchical level (the "composite" level) are shown in **Figure 7** for varying volume fractions of defects in the matrix. The figure illustrates the percentage degradation of stiffness, ultimate strain, strength, and toughness of the resulting composite, normalized with respect to the defect-free material. Simulations predict a toughness reduction of up to 87% and a strength reduction of up to 71% for a 10% defect content, whilst stiffness and ultimate strain reductions are smaller (51 and 27%, respectively) for the same defect volume fraction.

We also study the influence of reinforcement size or shape in the matrix. A single-scale model can be used in this case, and different aspect ratios l/w = 20, 40, 80, 120 for the reinforcements are considered, with the same staggered distribution and material properties as considered previously. Volume fractions of reinforcements between the different cases vary minimally between 5.5% (l/w = 20) and 6.5% (l/w = 120). Results are shown in **Figure 8** and demonstrate how increasing aspect ratios lead to a stiffening and strengthening of the composite. Since volume fractions are virtually the same for the various cases, these effects are due to the non-trivial stress distributions occurring in the matrix for increasing aspect ratios and reinforcement overlap, which lead to a more efficient stress transfer from the matrix to the stiffer reinforcements. Such an effect cannot be captured, e.g., by a widely used simple rule of mixtures, and the example illustrates how structure, as well as material constituents, can control the overall mechanical behavior. The toughness of the composite remains approximately constant for varying reinforcement aspect ratios (a 4% decrease for l/w = 120 with respect to l/w = 20). This is because two competing effects are at play: on the one hand, crack path lengths increase for increasing aspect ratios (schematically shown in **Figure 8**, top), which implies an increase in toughness, on the other hand, the discussed overall material stiffening leads to reduced ultimate strains and a decrease in toughness.

To further highlight this effect, a hierarchical structure is now considered. Results from **Figure 6** show that in a uniaxially loaded heterogeneous structure the matrix sub-regions are subjected to varying loading conditions of tension/shear, depending on the direction of the reinforcements. Therefore, if we consider a hierarchical composite with multiple scale levels each presenting a specific architecture, it is clear that the overall behavior depends on the direction-dependent response of the structure at the lower

scale, which must be derived in different loading conditions and subsequently injected at the next hierarchical level. To illustrate this, we consider a 2-level hierarchical composite constituted

by a directionally defected/porous material at the lowest scale, schematically shown in **Figure 9**, which can appropriately model various systems such as fiber-reinforced foams (Vaikhanski and Nutt, 2003), bioscaffolds (Huang et al., 2013), or the so-called nanolattices or metamaterials (Meza et al., 2014). The hollow structure considered here contains a 20% volume fraction of voids (thus reducing the overall density and increasing the specific stiffness of the matrix, as in Meza et al. (2014), and is composed of overlapping "pores." Simulations show that the representative RVE for this structure displays an anisotropic response in its stress–strain curves, also shown in **Figure 9** (top right). The stiffness and strength as well as toughness of the structure increase with alignment of the load along the direction of the elongated pores, when the beam-like elements are working under tension, as opposed to loading perpendicular to the pores, when the beamlike elements are mainly subjected to bending (as highlighted in **Figure 9**, bottom).

These results obtained for RVEs at the first hierarchical level are then used as the input of for the second hierarchical level, which is a fiber-matrix structure similar to that considered previously (e.g., in **Figure 6**). The fibers are again 100 times stiffer and stronger than the matrix material, not considering voids. The fiber aspect ratio is l/w = 50, with a volume fraction of 6.7% and maximal reinforcement overlap. The results reported in **Figure 10** highlight the impact of anisotropy at the first hierarchical level on overall properties. In the first case (A), the matrix RVEs are oriented so as to obtain maximal strength in the direction of the fibers. Therefore, crack nucleation at the reinforcement tips requires a higher load level. When the matrix in the tip region fails, debonding directly ensues due to the low strength required to damage the matrix in the fiber direction. It is important to notice that only few cracks have opened during the test, and that the evolution to a traversing crack is direct, leading to a brittle fracture with virtually elastic behavior up to failure. In the second case (B), the low strength of the matrix in the fiber direction leads to the opening of various cracks at reinforcement tips. In this case, a higher load can be transferred between fibers before failure occurs, and the overall stress–strain behavior is therefore elastoplastic. Maximum strength is obtained for case A, and maximum ultimate strain and toughness for case B.

Despite its simplicity, this case study highlights the complex interaction between structural architectures at different hierarchical levels. Even considering only two levels of hierarchy, it is possible to simulate and design preferential material orientations depending on the type of expected loading, in order to optimize and calibrate the architectures for improved global material response. This becomes more feasible when increasing the number of hierarchical levels (e.g., seven in the case of bone) and the corresponding range of size scales involved, from nano to macro.

# **Conclusion**

Complex hierarchically built composite structures are the way forward to design materials with exceptional or tailor-made mechanical behavior. Designing and realizing these structures is experimentally challenging, and therefore the need for reliable numerical models to characterize and optimize their behavior is acute. In this study, we have presented a novel numerical HLSM to design structural biomaterials and hierarchical nanocomposites, which correctly accounts for defects/porosity and material heterogeneity at different hierarchical levels, and allows the evaluation of the corresponding stress concentrations, crack nucleation, and propagation. As examples, the influence of reinforcement aspect ratios and defect distributions on crack nucleation and propagation inside the matrix have been studied, using a multiscale approach, in order to model nacre-like composites. Simulations show how structure at lower hierarchical levels affects global properties such as stiffness, ultimate strain, strength, and toughness. Moreover, results show that interaction between these lower level hierarchical properties and macroscopic properties such as the overlapping of the reinforcements in the matrix are key parameters that lead to non-trivial changes in material behavior, e.g., from brittle to ductile. Numerical simulations carried out by Sen and Buehler (2011) show similar results on the behavior of the hierarchical composite structure, i.e., from linear elastic base components, an elasto-plastic stress–strain behavior deriving from crack nucleation and stopping emerges. The additional contribution of the present paper is to further investigate how the composite organization and hierarchical architecture strongly influence the mechanical properties of the material and how it is possible to optimize those using specific criteria, using the same starting matrix/reinforcement materials. This is obtained, for example, by varying the dimension of the reinforcements at a single scale level and considering different orientations in a 2-level hierarchical material.

Thus, the developed HLSM 2-D code has shown reliability in simulating heterogeneity, hierarchy, and fracture propagation inside structured layers, showing promise for the application to more complex hierarchical structures. Thanks to the rapidly developing field of nanocomposite manufacture, it is already possible to artificially create materials with multi-scale hierarchical reinforcements. The developed code could be a valuable support in the design and optimization of these advanced materials, drawing inspiration and going beyond biological materials with exceptional mechanical properties. Further, the model can prove useful in the modeling of 2-D materials such as graphene and graphene-reinforced composites, where defects can play a major role in determining the overall mechanical properties, and their composites. In future works, more complex geometries will be studied by considering 3-dimensional structures and more scale levels (as observed in nature) to further capture optimization aspects in toughening/stiffening

# **References**


strategies in hierarchical materials. Other interesting features that are potentially linked to energy dissipation in biomaterials will also be evaluated, such as non-linear properties of the matrix, non-local bonding, and imperfect matrix/reinforcement interface properties.

# **Acknowledgments**

NP is supported by the European Research Council (ERC StG Ideas 2011 BIHSNAM n. 279985 on Bio-Inspired hierarchical supernanomaterials, ERC PoC 2013-1 REPLICA2 n. 619448 on Large-area replication of biological antiadhesive nanosurfaces, ERC PoC 2013-2 KNOTOUGH n. 632277 on Supertough knotted fibers), by the European Commission under the Graphene Flagship (WP10 "Nanocomposites," n. 604391) and by the Provincia Autonoma di Trento ("Graphene Nanocomposites," n. S116/2012- 242637 and delib. reg. n. 2266). LB and FB acknowledge support from BIHSNAM. Computational resources were provided by HPC@POLITO, a project of Academic Computing within the Department of Control and Computer Engineering at the Politecnico di Torino (http://www.hpc.polito.it), for which we acknowledge the head of the Department, Prof. F. Della Croce.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Brely, Bosia and Pugno. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Computational modeling of discrete mechanical systems and complex networks: where we are and where we are going

# **Andrea Infuso<sup>1</sup>\* and Marco Paggi <sup>2</sup>**

<sup>1</sup> Department of Structural, Geotechnical and Building Engineering, Politecnico di Torino, Torino, Italy

2 IMT Institute for Advanced Studies Lucca, Lucca, Italy

\*Correspondence: andrea.infuso@polito.it

#### **Edited by:**

Maurizio Dapor, Bruno Kessler Foundation, Italy

#### **Reviewed by:** Stefano Mariani, Politecnico di Milano, Italy

**Keywords: discrete systems, molecular dynamics, networks theory, cross-disciplinary analogies**

### **WHERE WE ARE IN MODELING MECHANICAL, TRAFFIC, AND SOCIO-ECONOMICAL DISCRETE SYSTEMS**

Discrete systems have been firstly introduced in physics (Fisher and Wiodm, 1969; Noor and Nemeth, 1980; Adhikari et al., 1996; Kornyak, 2009) to simulate materials at the micro- and nano-scales where a continuum description of matter breaks down. The constituents can be atoms or molecules and their interactions are usually modeled by force fields resulting from chemical potentials or weak van der Waals interactions, depending on the type of bonding. These models have been further exploited in mechanics with the aim of predicting macroscopic properties such as strength and toughness from the non-linear interactions taking place between the constituents at the different scales. Pioneering attempts in mechanics to model discrete mechanical systems are those using lattice models (van Mier et al., 1995; Schlangen and Garboczi, 1997) characterized by a network of nodes connected by links modeled by beams. Although proven to suffer from meshdependency, they have been broadly used for studying the non-linear fracture behavior of concrete at the meso-scale. Another line of research regards the generalized-Born approach (Pellegrini and Field, 2002; Marenich et al., 2013), firstly used in chemistry to model a solute represented as a set of three-dimensional spheres into a continuum medium solvent, then applied in molecular mechanics (called MM/GBSA) to investigate contact and fracture of bodies at the micro-scale.

The high-computational power and the development of powerful open source software allow nowadays the design of wide discrete scalable models composed of up to millions of particles or molecules whose equations of motions and mutual interactions are described by highly nonlinear interatomic potential laws. This is the field of molecular dynamics (MD), which led to the development of specific explicit time integration schemes (like the velocity-verlet integration scheme) to solve large systems of equations with a reduced computational cost. Car and Parrinello (1985) proposed a minimization of the total energy of the system by applying a dynamical simulated annealing based on MD. MD computations can also be coupled with continuum simulations by multi-scale methods. Among the many strategies available in the literature, (Shenoy et al., 1999; Knap and Ortiz, 2001) developed an approach based on the Tadmor's quasi-continuum method (Tadmor et al., 1996) operating on a representative atomistic zone with a reduced number of degrees of freedom. Local minima of the whole system potential energy are determined via the total energy from a cluster of atoms, avoiding the complete calculation of the full atomistic force field. The MD enriched continuum method by Belytschko and Xiao (2003) and Xiao and Belytschko (2003) was also another pioneering approach to couple a potential energy Hamiltonian calculation conducted on a fine scale MD domain with a Lagrangian calculation on a coarse scale continuum domain with an

overlapped bridging domain among the two representations. Recently, an implementation of interatomic potential laws within a displacement-based finite element (FE) formulation has also been proposed in Nasdala et al. (2010), with a rigorous implicit solution scheme, aiming at generating models where non-linear discrete and continuous systems can be suitably combined.

Discrete systems made of nodes and links are also used in other disciplines than mechanics to model transport or socio-economical networks (Caldarelli and Vespignani, 2007;Whrittle, 2012). Based on a continuum version of traffic conservation along a link, Lighthill and Whitham (1955) and Whitham (1974) and independently Richards (1956) proposed the LWR model where the governing equation describing the dependency of the traffic flux on time and on location along a link is a nonlinear hyperbolic partial differential equation (PDE) analogous to that describing the propagation of the front of a wave inside a medium. For a homogenous link where the velocity of traffic is the same at any location and no shocks on the traffic flow are present, the integration of the LWR PDE leads to a non-linear relation between the traffic flow and the density of vehicles, which fully represents the traffic state. Also, in economics, it is of great interest to quantify the effect of breaking a link over the whole network response by simulating the dynamics of flow redistribution. Again a *flow model* can be used as suggested in Zhou et al. (2010) to decode a huge human crowd without distinction between individuals, while small groups are better represented by an *agent-based approach*. Even if computationally demanding, the latter approach can be exploited in economics for competitive resource allocation (Chakraborti et al., 2015) with models characterized by strong heterogeneity and non-equilibrium dynamics among adaptive interacting agents. Similar problems are likewise of paramount importance in socio-economical systems, as, e.g., the prediction of the propagation of false news (trolling behavior) in the web, or the understanding of the socio-economic impact of removing a commercial or a web node (Caldarelli et al., 2008; Quattrociocchi et al., 2013). A graphical representation in **Figure 1** of a mechanical system (**Figure 1A**) and a socio-economic network (**Figure 1B**) emphasizes also the visual similarity between them.

#### **WHERE WE ARE GOING: COMPLEXITY OF DISCRETE SYSTEMS, ANALOGIES, AND DIFFERENCES**

All the previous discrete systems and probably many others analyzed by the scientific Community share the common features of complex systems where the overall properties *emerge* from the non-local nonlinear interactions between their basic network components and can only be predicted by numerically simulating the system response and the dynamic evolution of defects (cracks, node removal, link removal). Failure of the system is generally the result of a percolation of defects at different scales, which leads to complex redistribution of internal forces/flows. Understanding how the system is able to withstand perturbations, i.e., its *resilience* or *flaw-tolerance*, is a problem common to all of the disciplines mentioned in the previous section.

Hence, strong analogies and differences emerge from a deep analysis of network models if examined from a crossdisciplinary perspective. The integral formulation of LWR PDE leads to a non-linear relation between the flux and the vehicle density, which has an ascending branch up to a maximum and then decreased to 0 by further increasing the vehicle density. Interatomic potentials providing the constitutive behavior of a link are also relations between interaction forces and relative displacements of nodes with a similar non-linear trend. The equilibrium of each node is also a common concept: it can be the result of the net sum of the fluxes converging to a node or the net sum of the nodal forces. The removal of a link is typically the result of a propagation of a dislocation inside the material, whereas it represents an interruption of communication between nodes in traffic or in socioeconomic systems. Non-local interactions can be considered both for mechanical and socio-economical systems and can be modeled as parallel links connecting nodes, in addition to serial links connecting only neighboring nodes. Therefore, given the elementary constitutive equations describing the link behavior that can be mechanically modeled as non-linear springs and the state variables characterizing the nodes, the basic dynamics of any discrete system can be in principle simulated according to the numerical techniques proper of non-linear mechanics. A non-locality index can also be used to classify and distinguish between different networks, as shown in Infuso and Paggi (2014) to interpret the response of a discrete system upon removal of nodes in different locations. From numerical simulations in Infuso and Paggi (2014), we observed that the higher the value of the non-locality index, the higher the total force supported by the network, for the same type of node removal.

Significant differences are also present among discrete systems in different disciplines. Links between atoms or molecules are always constructed on a deterministic basis, i.e., all the atoms are linked by a certain amount of links, just with different intensities depending on the spatial distance between the nodes. The situation in socio-economic systems is in reality much more complex, since the existence of links between nodes may depend on the reputation or the economic state of the nodes themselves, to cite two possible influencing state variables. Therefore, links are created dynamically and are probably the result of an optimization problem that should be modeled and solved at the local level. Although these differences may suggest that a unified theory for the simulation of complex networks is still to come and should be the result of a joint cross-disciplinary effort involving mathematicians, engineers, physicists, and economists, they also open interesting

perspectives for mechanics. For instance, the fact that local optimization problems should take place at the local level, in mechanics they could be the result of the mechanical response of a constituent of a metamaterial, whose configuration is such that it can minimize or maximize certain properties depending on the intensity of diffusive external fields (temperature, humidity, electro-chemical potentials, etc.), as it happens in biomechanical applications. Within a bottom-up approach to model this cyber-physical system, given the links formed at the lower scale and their constitutive behavior that may depend on time, the global mechanical response of the network can then be simulated by minimizing the total energy of the system, using numerical methods proper of mechanics.

#### **ACKNOWLEDGMENTS**

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007–2013)/ERC Grant Agreement No. 306622 (ERC Starting Grant "Multifield and Multi-scale Computational Approach to Design and Durability of PhotoVoltaic Modules" – CA2PVM).

## **REFERENCES**


computational aspects. *Eng. Fract. Mech.* 57, 319–332. doi:10.1016/S0013-7944(97)00010-6


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 22 December 2014; accepted: 21 February 2015; published online: 09 March 2015.*

*Citation: Infuso A and Paggi M (2015) Computational modeling of discrete mechanical systems and complex networks: where we are and where we are going. Front. Mater. 2:18. doi: 10.3389/fmats.2015.00018*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials.*

*Copyright © 2015 Infuso and Paggi. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Strain localization and shear band propagation in ductile materials

# **Nicola Bordignon, Andrea Piccolroaz, Francesco Dal Corso and Davide Bigoni \***

Department of Civil, Environmental and Mechanical Engineering (DICAM), University of Trento, Trento, Italy

#### **Edited by:**

Simone Taioli, Bruno Kessler Foundation, Italy

#### **Reviewed by:**

Gianni F. Royer Carfagni, University of Parma, Italy Stefano Vidoli, Sapienza University of Rome, Italy

**\*Correspondence:** Davide Bigoni, via Mesiano 77, Trento I-38123, Italy e-mail: bigoni@ing.unitn.it

A model of a shear band as a zero-thickness non-linear interface is proposed and tested using finite element simulations. An imperfection approach is used in this model where a shear band that is assumed to lie in a ductile matrix material (obeying von Mises plasticity with linear hardening), is present from the beginning of loading and is considered to be a zone in which yielding occurs before the rest of the matrix. This approach is contrasted with a perturbative approach, developed for a J2-deformation theory material, in which the shear band is modeled to emerge at a certain stage of a uniform deformation. Both approaches concur in showing that the shear bands (differently from cracks) propagate rectilinearly under shear loading and that a strong stress concentration should be expected to be present at the tip of the shear band, two key features in the understanding of failure mechanisms of ductile materials.

**Keywords: slip lines, plasticity, failure, stress concentration, rectilinear growth**

#### **1. INTRODUCTION**

When a ductile material is brought to an extreme strain state through a uniform loading process, the deformation may start to localize into thin and planar bands, often arranged in regular lattice patterns. This phenomenon is quite common and occurs in many materials over a broad range of scales: from the kilometric scale in the earth crust (Kirby, 1985), down to the nanoscale in metallic glass (Yang et al., 2005), see the examples, reported in **Figure 1**.

After localization, unloading typically<sup>1</sup> occurs in the material outside the bands, while strain quickly evolves inside, possibly leading to final fracture (as in the examples shown in **Figure 2**, where the crack lattice is the signature of the initial shear band network<sup>2</sup> ) or to a progressive accumulation of deformation bands (as, for instance, in the case of the drinking straws, or of the iron meteorite, or of the uPVC sample shown in **Figure 1**, or in the well-known case of granular materials, where fracture is usually absent and localization bands are made up of material at a different relative density, Gajo et al., 2004).

It follows from the above discussion that as strain localization represents a prelude to failure of ductile materials, its mechanical understanding paves the way to the innovative use of materials in extreme mechanical conditions. For this reason, shear bands have been the focus of a thorough research effort. In particular, research initiated with pioneering works by Hill (1962), Nadai (1950), Mandel (1962), Prager (1954), Rice (1977), Thomas (1961), and developed – from theoretical point of view – into two principal directions, namely, the dissection of the specific constitutive features responsible for strain localization in different materials (for instance, as related to the microstructure, Danas and Ponte Castaneda, 2012; Bacigalupo and Gambarotta, 2013; Tvergaard, 2014) and the struggle for the overcoming of difficulties connected with numerical approaches [reviews have been given by Needleman and Tvergaard (1983) and Petryk (1997)]. Although these problems are still not exhausted, surprisingly, the most important questions have only marginally been approached and are therefore still awaiting explanation. These are as follows: (i) Why are shear bands a preferred mode of failure for ductile materials? (ii) Why do shear bands propagate rectilinearly under mode II, while cracks do not? (iii) How does a shear band interact with a crack or with a rigid inclusion? (iv) Does a stress concentration exist at a shear band tip? (v) How does a shear band behave under dynamic conditions?

The only systematic<sup>3</sup> attempt to solve these problems seems to have been a series of works by Bigoni and co-workers, based on the perturbative approach to shear bands (Bigoni and Capuani, 2002, 2005; Piccolroaz et al., 2006; Argani et al., 2013, 2014). In fact, problems (i), (ii), and (iv) were addressed in Bigoni and Dal Corso (2008) and Dal Corso and Bigoni (2010), problem (iii) in Bigoni et al. (2008), Dal Corso et al. (2008), and Dal Corso and Bigoni (2009), and (v) in Bigoni and Capuani (2005).

The purpose of the present article is to present a model of a shear band as a zero-thickness interface and to rigorously motivate

<sup>1</sup>For granular materials, there are cases in which unloading occurs inside the shear band, as shown by Gajo et al. (2004).

<sup>2</sup>The proposed explanation for the crack patterns shown in **Figure 2** relies on the fact that the fracture network has formed during the plastic evolution of a ductile homogeneously deformed material. Other explanations may be related to bonding of an external layer to a rigid substrate (Peron et al., 2013), or to surface instability (Destrade and Merodio, 2011; Boulogne et al., 2015), or to instabilities occurring during shear (Destrade et al., 2008; Ciarletta et al., 2013).

<sup>3</sup> Special problems of shear band propagation in geological materials have been addressed by Puzrin and Germanovich (2005) and Rice (1973).

#### **FIGURE 1 | Continued**

From left to right, starting from the upper part: A merlon in the Finale Emilia castle failed (during the Emilia earthquake on May 20, 2012) in compression with a typical "X-shaped" deformation band pattern (bricks are to be understood here as the microstructure of a composite material). A sedimentary rock with the signature of an "X-shaped" localization band (infiltrated with a different mineral after formation). A stone axe from a British Island (Museum of Edinburgh) evidencing two parallel localization bands and another at a different orientation. A runestone (Museum of Edinburgh) with several localized deformation bands, forming angles of approximately 45°

between each other. A polished and etched section of an iron meteorite showing several alternate bands of kamacite and taenite. Deformation bands in a strip of unplasticized poly(vinyl chloride) (uPVC) pulled in tension and eventually evolving into a necking. An initially regular hexagonal disposition of drinking straws subject to uniform uniaxial strain has evolved into an "X-shaped" localization pattern. A fracture prevails on a regularly distributed network of cracks in a vault of the Amiens dome. "X-shaped" localization bands in a kaolin sample subject to vertical compression and lateral confining pressure. A thin, isolated localization band in a sedimentary layered rock (Silurian formation near Aberystwyth).

**FIGURE 2 | Regular patterns of localized cracks as the signature of strain localization lattices**. From left to right: Dried mud; Lava cracked during solidification (near Amboy crater); Bark of a maritime pine (Pinus pinaster); Cracks in a detail of a painting by J. Provost ("Saint Jean-Baptiste," Valenciennes, Musée des Beaux Arts).

this as the asymptotic behavior of a thin layer of material, which is extremely compliant in shear (Section 2). Once the shear band model has been developed, it is used (in Section 3) to demonstrate two of the above-listed open problems, namely (ii) that a shear band grows rectilinearly under mode II remote loading in a material deformed near to failure and (iv) to estimate the stress concentration at the shear band tip. In particular, a preexisting shear band is considered to lie in a matrix as a thin zone of material with properties identical to the matrix, but lower yield stress. This is an imperfection, which remains neutral until the yield is reached in the shear band<sup>4</sup> . The present model is based on an imperfection approach and shares similarities to that pursued by Abeyaratne and Triantafyllidis (1981) and Hutchinson and Tvergaard (1981), so that it is essentially different from a perturbative approach, in which the perturbation is imposed at a certain stage of a uniform deformation process<sup>5</sup> .

#### **2. ASYMPTOTIC MODEL FOR A THIN LAYER OF HIGHLY COMPLIANT MATERIAL EMBEDDED IN A SOLID**

A shear band, inside a solid block of dimension *H*, is modeled as a thin layer of material (of semi-thickness *h*, with *h*/*H* 1) yielding at a uniaxial stress σ (*s*) *Y* , which is lower than that of the surrounding matrix σ (*m*) *Y* , **Figure 3**. Except for the yield stress, the material inside and outside the layer is described by the same elastoplastic model, namely, a von Mises plasticity with associated flow rule and linear hardening, defined through the elastic constants, denoted by the Young modulus *E* and Poisson's ratio *v*, and the plastic modulus *Ep*, see **Figure 3B**.

At the initial yielding, the material inside the layer [characterized by a low hardening modulus *Eep* = *EEp*/(*E* + *Ep*)] is much more compliant than the material outside (characterized by an elastic isotropic response *E*).

For *h*/*H* 1, the transmission conditions across the layer imply the continuity of the tractions, *t* = [*t* <sup>21</sup>, *t* <sup>22</sup>] *T* , which can be expressed in the asymptotic form

$$\left\| t\_{21} \right\| = O(h), \quad \left\| t\_{22} \right\| = O(h), \tag{1}$$

where [[·]] denotes the jump operator.

The jump in displacements, [[*u*]] = [[[*u*1]], [[*u*2]]]*<sup>T</sup>* , across the layer is related to the tractions at its boundaries through the asymptotic relations (Mishuris et al., 2013; Sonato et al., 2015)

<sup>4</sup>A different approach to investigate shear band evolution is based on the exploitation of phase-field models (Zheng and Li, 2009), which has been often used for brittle fracture propagation (Miehe et al., 2010).

<sup>5</sup>To highlight the differences and the analogies between the two approaches, the incremental strain field induced by the emergence of a shear band of finite length (modeled as a sliding surface) is determined for a J2-deformation theory material and compared with finite element simulations in which the shear band is modeled as a zero-thickness layer of compliant material.

of highly compliant material (Eep/E 1) embedded in a material block characterized by a dimension H, such that h/H 1; both materials obey the

reported in **(B)**, but having a different yield stress (lower inside than outside the shear band).

$$\mathfrak{a}\_{21}(\|\![\mu\_1]\![\mu\_2]\!) = \frac{E\_{\mathcal{P}}\sqrt{\Im\|\mu\_1\|^2 + 4\|\![\mu\_2]\|^2} + 6h\sigma\_Y^{(s)}}{(3E + 2(1+\nu)E\_{\mathcal{P}})\sqrt{\Im\|\mu\_1\|^2 + 4\|\![\mu\_2]\|^2}} \frac{E\|\![\mu\_1]\!\rangle}{2h} + O(h),\tag{2}$$

$$\mathfrak{a}\_{22}(\|\boldsymbol{\mu}\_{1}\|,\|\boldsymbol{\mu}\_{2}\|)=\frac{(E+2(1-\nu)E\_{\boldsymbol{p}})\sqrt{3\|\boldsymbol{\mu}\_{1}\|^{2}+4\|\boldsymbol{\mu}\_{2}\|^{2}}+8h(1-2\nu)\sigma\_{Y}^{(\boldsymbol{\mu})}}{(1-2\nu)(3E+2(1+\nu)E\_{\boldsymbol{p}})\sqrt{3\|\boldsymbol{\mu}\_{1}\|^{2}+4\|\boldsymbol{\mu}\_{2}\|^{2}}}\frac{E\|\boldsymbol{\mu}\_{2}\|}{2h}+O(h),\tag{3}$$

involving the semi-thickness *h* of the shear band, which enters the formulation as a *constitutive parameter for the zero-thickness interface model* and introduces a *length scale*. Note that, by neglecting the remainder *O(h)*, equations (2) and (3) define non-linear relationships between tractions and jump in displacements.

The time derivative of equations (2) and (3) yields the following asymptotic relation between incremental quantities

$$\dot{\boldsymbol{t}} \sim \left[ \frac{1}{h} \mathbf{K}\_{-1} + \mathbf{K}\_{0} (\|\boldsymbol{\mu}\_{1}\|, \|\boldsymbol{\mu}\_{2}\|) \right] \|\dot{\boldsymbol{u}}\|,\tag{4}$$

where the two stiffness matrices *K* <sup>−</sup><sup>1</sup> and *K*<sup>0</sup> are given by

$$\mathbf{K}\_{-1} = \frac{E}{2(3E + 2(1+\nu)E\_{\mathcal{P}})} \begin{bmatrix} E\_{\mathcal{P}} & 0 \\ 0 & \frac{E + 2(1-\nu)E\_{\mathcal{P}}}{1 - 2\nu} \end{bmatrix}, \qquad (5)$$

$$\mathbf{K}\_{0} = \frac{12E\sigma\_{\mathcal{Y}}^{(\ast)}}{\left(3E + 2(1+\nu)E\_{\mathcal{P}}\right)\left(3\|\|u\_{1}\|\|^{2} + 4\|\|u\_{2}\|\|^{2}\right)^{3/2}}$$

$$\times \begin{bmatrix} \|\|u\_{2}\|\|^{2} & -\|\|u\_{1}\|\|\|u\_{2}\|\|\\ -\|\|u\_{1}\|\|\|u\_{2}\|\| & \|\|u\_{1}\|\|^{2} \end{bmatrix}, \tag{6}$$

Assuming now a perfectly plastic behavior, *E<sup>p</sup>* = 0, in the limit *h*/*H*→0 the condition

$$\|\boldsymbol{\mu}\_2\| = 0 \tag{7}$$

is obtained, so that the incremental transmission conditions equation (4) can be approximated to the leading order as

$$
\dot{\mathbf{t}} \sim \frac{1}{h} \mathbf{K}\_{-1} \| \dot{\mathbf{u}} \|. \tag{8}
$$

Therefore, when the material inside the layer is close to the perfect plasticity condition, the incremental conditions assume the limit value

$$
\dot{\iota}\_{21} = 0, \quad \|\dot{\iota}\_2\| = 0,\tag{9}
$$

which, together with the incremental version of equation (1)2, namely,

$$\|\dot{\mathbf{f}}\|\!\_{22}\| = \mathbf{0},\tag{10}$$

correspond to the incremental boundary conditions proposed in Bigoni and Dal Corso (2008) to define a pre-existing shear band of null thickness.

The limit relations, equations (9) and (10), motivate the use of the imperfect interface approach (Bigoni et al., 1998;Antipov et al., 2001; Mishuris, 2001, 2004; Mishuris and Kuhn, 2001; Mishuris and Ochsner, 2005, 2007) for the modeling of shear band growth in a ductile material. A computational model, in which the shear bands are modeled as interfaces, is presented in the next section.

### **3. NUMERICAL SIMULATIONS**

Two-dimensional plane-strain finite element simulations are presented to show the effectiveness of the above-described asymptotic model for a thin and highly compliant layer in modeling a shear band embedded in a ductile material. Specifically,we will show that the model predicts rectilinear propagation of a shear band under simple shear boundary conditions and it allows the investigation of the stress concentration at the shear band tip.

The geometry and material properties of the model are shown in **Figure 4**, where a rectangular block of edges *H* and *L* ≥ *H* is subject to boundary conditions consistent with a simple shear deformation, so that the lower edge of the square domain is clamped, the vertical displacements are constrained along the vertical edges and along the upper edge, where a constant horizontal displacement *u*<sup>1</sup> is prescribed. The domain is made of a ductile material and contains a thin (*h*/*H* 1) and highly compliant (*Eep*/*E 1*) layer of length *H*/2 and thickness 2*h* = 0.01 mm, which models a shear band. The material constitutive behavior is described by an elastoplastic model based on linear isotropic elasticity (*E* = 200000 MPa, *v* = 0.3) and von Mises plasticity with linear hardening (the plastic modulus is denoted by *Ep*). The uniaxial yield stress σ (*m*) *Y* for the matrix material is equal to 500 MPa, whereas the layer is characterized by a lower yield stress, namely, σ (*s*) *<sup>Y</sup>* = 400 MPa.

The layer remains neutral until yielding, but, starting from that stress level, it becomes a material inhomogeneity, being more compliant (because its response is characterized by *Eep*) than the matrix (still in the elastic regime and thus characterized by *E*). The layer can be representative of a preexisting shear band and can be treated with the zero-thickness interface model, equations (2) and (3). This zero-thickness interface was implemented in the ABAQUS finite element

software<sup>6</sup> through cohesive elements, equipped with the tractionseparation laws, equations (2) and (3), by means of the user subroutine UMAT. An interface, embedded into the cohesive elements, is characterized by two dimensions: a geometrical and a constitutive thickness. The latter, 2*h*, exactly corresponds to the constitutive thickness involved in the model for the interface equations (2) and (3), while the former, denoted by 2*hg*, is related to the mesh dimension in a way that the results become independent of this parameter, in the sense that a mesh refinement yields results converging to a well-defined solution.

We consider two situations. In the first, we assume that the plastic modulus is *E<sup>p</sup>* = 150000 MPa (both inside and outside the shear band), so that the material is in a state far from a shear band instability (represented by loss of ellipticity of the tangent constitutive operator, occurring at *E<sup>p</sup>* = 0) when at yield. In the second, we assume that the material is prone to a shear band instability, though still in the elliptic regime, so that *E<sup>p</sup>* (both inside and outside the shear band) is selected to be "sufficiently small", namely, *E<sup>p</sup>* = 300 MPa. The pre-existing shear band is therefore employed as an imperfection triggering shear strain localization when the material is still inside the region, but close to the boundary, of ellipticity.

#### **3.1. DESCRIPTION OF THE NUMERICAL MODEL**

With reference to a square block (*L* = *H* = 10 mm) containing a pre-existing shear band with constitutive thickness *h* = 0.005 mm, three different meshes were used, differing in the geometrical thickness of the interface representing the pre-existing shear band (see **Figure 5**where the shear band is highlighted with a black line), namely, *h<sup>g</sup>* = {0.05; 0.005; 0.0005} mm corresponding to coarse, fine, and ultra-fine meshes.

The three meshes were generated automatically using the mesh generator available in ABAQUS. In order to have increasing mesh refinement from the exterior (upper and lower parts) to the interior (central part) of the domain, where the shear band is located, and to ensure the appropriate element shape and size according to the geometrical thickness 2*hg*, the domain was partitioned into rectangular subdomains with increasing mesh seeding from the exterior to the interior. Afterwards, the meshes were generated by employing a free meshing technique with quadrilateral elements and the advancing front algorithm.

The interface that models the shear band is discretized using 4-node two-dimensional cohesive elements (COH2D4), while the matrix material is modeled using 4-node bilinear, reduced integration with hourglass control (CPE4R).

It is important to note that the constitutive thickness used for traction-separation response is always equal to the actual size of the shear band *h* = 0.005 mm, whereas the geometric thickness *hg*, defining the height of the cohesive elements, is different for the three different meshes. Consequently, all the three meshes used in the simulations correspond to the same problem in terms of both material properties and geometrical dimensions (although the geometric size of the interface is different), so that

<sup>6</sup>ABAQUS Standard Ver. 6.13 has been used, available on the AMD Opteron cluster Stimulus at UniTN.

**FIGURE 5 | The three meshes used in the analysis to simulate a shear band (highlighted in black) in a square solid block (L** = **H** = **10 mm)**. The shear band is represented in the three cases as an interface with the same constitutive thickness h = 0.005 mm, but with

decreasing geometric thickness hg; **(A)** coarse mesh (1918 nodes, 1874 elements, h<sup>g</sup> = 0.05 mm); **(B)** fine mesh (32,079 nodes, 31,973 elements, h<sup>g</sup> = 0.005 mm); **(C)** ultra-fine mesh (1,488,156 nodes, 1,487,866 elements, h<sup>g</sup> = 0.0005 mm).

corresponds to the material at yielding σ<sup>12</sup> ≥ 500/ 3 = 288.68 MPa. Three different stages of deformation are shown, corresponding to a prescribed

figures are amplified by a deformation scale factor of 25 and the percentages refer to the final displacement.

the results have to be, and indeed will be shown to be, mesh independent<sup>7</sup> .

#### **3.2. NUMERICAL RESULTS**

Results (obtained using the fine mesh, **Figure 5B**) in terms of the shear stress component σ<sup>12</sup> at different stages of a deformation process for the boundary value problem sketched in **Figure 4** are reported in **Figures 6** and **7**.

In particular, **Figure 6** refers to a matrix with high plastic modulus, *E<sup>p</sup>* = 150000 MPa, so that the material is far from the shear band formation threshold. The upper limit of the contour levels was set to the value σ<sup>12</sup> = 500/ √ 3 ' 288.68 MPa, corresponding to the yielding stress of the matrix material. As a result, the gray zone in the figure represents the material at yielding, whereas the material outside the gray zone is still in the elastic regime. Three stages of deformation are shown, corresponding to: the initial yielding of the matrix material (left), the yielding zone occupying approximately one-half of the space between the shear band tip and the right edge of the domain (center), and the yielding completely linking the tip of the shear band to the boundary (right). Note that the shear band, playing the role of a material imperfection, produces a stress concentration at its tip. However, the region of high stress level rapidly grows and diffuses in the whole domain. At the final stage, shown in **Figure 6C**, almost all the matrix material is close to yielding.

**Figure 7** refers to a matrix with low plastic modulus, *E<sup>p</sup>* = 300 MPa, so that the material is close (but still in the elliptic regime) to the shear band formation threshold (*E<sup>p</sup>* = 0). Three stages of deformation are shown, from the condition of initial yielding of the matrix material near the shear band tip (left), to an intermediate condition (center), and finally to the complete yielding of a narrow zone connecting the shear band tip to the right boundary (right). In this case, where the material is prone to shear band localization, the zone of high stress level departs from the shear band tip and propagates toward the right. This propagation occurs in a highly concentrated narrow layer, rectilinear, and parallel to the pre-existing shear band. At the final stage of deformation, shown in **Figure 7C**, the layer of localized shear has reached the boundary of the block.

<sup>7</sup>Note that, in the case of null hardening, mesh dependency may occur in the simulation of shear banding nucleation and propagation (Needleman, 1988; Loret and Prevost, 1991, 1993). This numerical issue can be avoided by improving classical inelastic models through the introduction of characteristic length-scales (Lapovok et al., 2009; Dal Corso and Willis, 2011).

**FIGURE 7 | Contour plots of the shear stress** σ **<sup>12</sup> for the case of material close to shear band instability (E<sup>p</sup>** = **300 MPa)**. The gray region corresponds to the material at yielding σ<sup>12</sup> ≥ 500/ √ 3. Three different stages of deformation are shown, corresponding to a

prescribed displacement at the upper edge of the square domain u<sup>1</sup> = 0.0340 mm **(A)**, u<sup>1</sup> = 0.0351 mm **(B)**, u<sup>1</sup> = 0.03623 mm **(C)**. The displacements in the figures are amplified by a deformation scale factor of 27.

**FIGURE 8 | Contour plots of the shear deformation** γ **<sup>12</sup> for the case of material far from shear band instability (E<sup>p</sup>** = **150000 MPa)**. Three different stages of deformation are shown, corresponding to a prescribed

displacement at the upper edge of the square domain u<sup>1</sup> = 0.037418 mm **(A)**, u<sup>1</sup> = 0.037518 mm **(B)**, u<sup>1</sup> = 0.037522 mm **(C)**. The displacements in the figures are amplified by a deformation scale factor of 25.

Results in terms of the shear strain component γ <sup>12</sup>, for both cases of material far from, and close to shear band instability are reported in **Figures 8** and **9**, respectively. In particular, **Figure 8** shows contour plots of the shear deformation γ <sup>12</sup> for the case of a material far from the shear band instability

(*E<sup>p</sup>* = 150000 MPa) at the same three stages of deformation as those reported in **Figure 6**. Although the tip of the shear band acts as a strain raiser, the contour plots show that the level of shear deformation is high and remains diffused in the whole domain.

**Figure 9** shows contour plots of the shear deformation γ <sup>12</sup> for the case of a material close to the shear band instability (*E<sup>p</sup>* = 300 MPa), at the same three stages of deformation as those reported in **Figure 7**. It is noted that the shear deformation is localized along a rectilinear path ahead of the shear band tip, confirming results that will be reported later with the perturbation approach (Section 4).

The shear deformation γ <sup>12</sup> and the shear stress σ<sup>12</sup> along the *x*-axis containing the pre-existing shear band for the case of a material close to strain localization, *E<sup>p</sup>* = 300 MPa, are shown in **Figure 10**, upper and lower parts, respectively. Results are reported for the three meshes, coarse, fine and ultra-fine (**Figure 5**) and at the same three stages of deformation as those shown in **Figures 7** and **9**. The results appear to be mesh independent, meaning that the solution converges as the mesh is more and more refined.

The deformation process reported in **Figures 7**, **9**, and **10** can be described as follows. After an initial homogeneous elastic deformation (not shown in the figure), in which the shear band remains neutral (since it shares the same elastic properties with the matrix material), the stress level reaches σ<sup>12</sup> = 400/ √ 3 ' 230.9 MPa, corresponding to the yielding of the material inside the shear band. Starting from this point, the pre-existing shear band is activated, which is confirmed by a high shear deformation γ <sup>12</sup> and a stress level above the yield stress inside the layer, −5 mm < *x* < 0 (left part of **Figure 10**). The activated shear band induces a strain localization and a stress concentration at its tip, thus generating a zone of material at yield, which propagates to the right (central part of **Figure 10**) until collapse (right part of **Figure 10**).

In order to appreciate the strain and stress concentration at the shear band tip, a magnification of the results shown in **Figure 10** in the region −0.2 mm < *x* < 0.2 mm is presented in **Figure 11**. Due to the strong localization produced by the shear band, only the ultra-fine mesh is able to capture accurately the strain and stress raising (blue solid curve), whereas the coarse and fine meshes smooth over the strain and stress levels (red dotted and green dashed curves, respectively). The necessity of an ultra-fine mesh to capture details of the stress/strain fields is well-known in computational fracture mechanics,where special elements (quarter-point or extended elements) have been introduced to avoid the use of these ultra-fine meshes at corner points.

For the purpose of a comparison with an independent and fully numerical representation of the shear band, a finite element simulation was also performed, using standard continuum elements (CPE4R) instead of cohesive elements (COH2D4) inside the layer. This simulation is important to assess the validity of the asymptotic model of the layer presented in Section 2. In this simulation, reported in **Figure 12**, the layer representing the shear band is a "true" layer of a given and finite thickness, thus influencing the results (while these are independent of the geometrical thickness 2*h<sup>g</sup>* of the cohesive elements, when the constitutive thickness 2*h* is the same). Therefore, only the fine mesh, shown in **Figure 5B**, was used, as it corresponds to the correct size of the shear band. The coarse mesh (**Figure 5A**) and the ultra-fine mesh (**Figure 5C**) would obviously produce different results, corresponding, respectively, to a thicker or thinner layer. Results pertaining to the asymptotic model, implemented into the traction-separation law for the cohesive elements COH2D4, are also reported in the figure (red solid curve) and are spot-on with the results obtained with a fully numerical solution employing standard continuum elements CPE4R (blue dashed curve).

A mesh of the same size as that previously called "fine" was used to perform a simulation of a rectangular block (*H* = 10 mm,

**FIGURE 11 | Shear and stress concentration at the shear band tip**. Shear deformation γ <sup>12</sup> **(A–C)** and shear stress σ <sup>12</sup> **(D–F)** along the x-axis containing the pre-existing shear band for the case of a material close to a shear band instability E<sup>p</sup> = 300 MPa. Three different stages of deformation are shown, corresponding to a prescribed displacement at the upper edge of the square domain u<sup>1</sup> = 0.0340 mm (left), u<sup>1</sup> = 0.0351 mm (center), u<sup>1</sup> = 0.03623 mm (right).

**idealizations for the shear band: zero-thickness model (discretized with cohesive elements, COH2D4) versus a true layer description (discretized with CPE4R elements)**. Shear deformation γ <sup>12</sup> **(A–C)** and shear stress σ <sup>12</sup> **(D–F)** along the horizontal line y = 0 containing the

instability E<sup>p</sup> = 300 MPa. Three different stages of deformation are shown, corresponding to a prescribed displacement at the upper edge of the square domain u<sup>1</sup> = 0.0340 mm (left), u<sup>1</sup> = 0.0351 mm (center), u<sup>1</sup> = 0.03623 mm (right).

*L* = 4 *H* = 40 mm) made up of a material close to shear band instability (*E<sup>p</sup>* = 300 MPa) and containing a shear band (of length *H*/2 = 5 mm and constitutive thickness 2*h* = 0.01 mm). Results are presented in **Figure 13**. In parts (**Figures 13A,B**) (the latter is a detail of part **Figure 13A**) of this figure the overall response curve is shown of the block in terms of average shear stress σ¯<sup>12</sup> = *T*/*L* (*T* denotes the total shear reaction force at the upper edge of the block) and average shear strain γ¯<sup>12</sup> = *u*1/*H*. In part (**Figure 13C**) of the figure contour plots of the shear deformation γ <sup>12</sup> are reported at different stages of deformation. It is clear that the

**FIGURE 13 | Results for a rectangular domain (L** = **40 mm, H** = **10 mm) of material close to shear band instability (E<sup>p</sup>** = **300 MPa) and containing a pre-existing shear band (of length H/2** = **5 mm and constitutive thickness 2h** = **0.01 mm)**. **(A)** Overall response curve of the block in terms of average shear stress σ¯<sup>12</sup> = T /L, where T is the total shear reaction force at the upper edge of the block, and average shear strain γ¯<sup>12</sup> = u1/H. **(B)** Magnification of

the overall response curve σ¯<sup>12</sup> − ¯γ<sup>12</sup> around the stress level corresponding to the yielding of the shear band. **(C)** Contour plots of the shear deformation γ <sup>12</sup> at different stages of deformation, corresponding to the points along the overall response curve shown in **(B)** of the figure. The deformation is highly focused along a rectilinear path emanating from the shear band tip. The displacements in the figures are amplified by a deformation scale factor of 50.

deformation is highly focused along a rectilinear path emanating from the shear band tip, thus demonstrating the tendency of the shear band toward rectilinear propagation under shear loading.

and referred there as (**Figures 10A,C**). These results, which have been obtained with the fine mesh, show that a strong incremental strain concentration develops at the shear band tip and becomes qualitatively similar to the square-root singularity found in the perturbative approach.

### **4. THE PERTURBATIVE VERSUS THE IMPERFECTION APPROACH**

With the perturbative approach, a perturbing agent acts at a certain stage of uniform strain of an infinite body, while the material is subject to a uniform prestress. Here, the perturbing agent is a pre-existing shear band,modeled as a planar slip surface, emerging at a certain stage of a deformation path (Bigoni and Dal Corso, 2008), in contrast with the imperfection approach in which the imperfection is present from the beginning of the loading.

With reference to a *x*<sup>1</sup> − *x*<sup>2</sup> coordinate system (inclined at 45° with respect to the principal prestress axes *x<sup>I</sup>* – *xII*), where the incremental stress*t*˙*ij* and incremental strain ε˙*ij* are defined (*i*, *j* = 1, 2), the incremental orthotropic response under plane-strain conditions(ε˙*i*<sup>3</sup> = 0)for incompressible materials(ε˙<sup>11</sup> + ˙ε<sup>22</sup> = 0) can be expressed through the following constitutive equations (Bigoni, 2012) 8 .

$$\dot{t}\_{11} = 2\mu \dot{\varepsilon}\_{11} + \dot{p}, \quad \dot{t}\_{22} = -2\mu \dot{\varepsilon}\_{11} + \dot{p}, \quad \dot{t}\_{12} = \mu\_\* \dot{\varphi}\_{12}, \quad (11)$$

where *p*˙ is the incremental in-plane mean stress, while µ and µ<sup>∗</sup> describe the incremental shear stiffness, respectively, parallel and inclined at 45° with respect to prestress axes.

The assumption of a specific constitutive model leads to the definition of the incremental stiffness moduli µ and µ\*. With reference to the J2-deformation theory of plasticity (Bigoni and Dal Corso, 2008), particularly suited to model the plastic branch of the constitutive response of ductile metals, the in-plane deviatoric stress can be written as

$$\mathbf{t}\_I - \mathbf{t}\_{II} = k \varepsilon\_I |\varepsilon\_I|^{(N-1)}.\tag{12}$$

In equation (12), *k* represents a stiffness coefficient and *N* ∈ (0, 1] is the strain hardening exponent, describing perfect plasticity (null hardening) in the limit *N*→0 and linear elasticity in the limit *N*→1. For the J2-deformation theory, the relation between the two incremental shear stiffness moduli can be obtained as

$$
\mu\_\* = N\mu,\tag{13}
$$

so that a very compliant response under shear (µ<sup>∗</sup> µ) is described in the limit of perfect plasticity *N*→0.

The perturbative approach (Bigoni and Dal Corso, 2008) can now be exploited to investigate the growth of a shear band within a solid. To this purpose, an incremental boundary value problem is formulated for an infinite solid, containing a zero-thickness pre-existing shear band of finite length 2*l* parallel to the *x*<sup>1</sup> axis

(see **Figure 15**) and loaded at infinity through a uniform shear deformation γ˙<sup>∞</sup> 12 .

The incremental boundary conditions introduced by the presence of a pre-existing shear band can be described by the following equations:

$$\dot{t}\_{21}(\mathbf{x}\_{\mathrm{l}},0^{\pm}) = 0,\ \|\dot{t}\_{22}(\mathbf{x}\_{\mathrm{l}},0)\| = 0,\ \|\dot{u}\_{2}(\mathbf{x}\_{\mathrm{l}},0)\| = 0,\ \forall |\mathbf{x}\_{\mathrm{l}}| < l.\tag{14}$$

A stream function ψ(*x*1, *x*2) is now introduced, automatically satisfying the incompressibility condition and defining the incremental displacements *u*˙*<sup>j</sup>* as *u*˙<sup>1</sup> = ψ,2, and *u*˙<sup>2</sup> = −ψ,1. The incremental boundary value problem is therefore solved as the sum of ψ°(*x*1, *x*2), solution of the incremental homogeneous problem, and ψ *p* (*x*1, *x*2), solution of the incremental perturbed problem.

The incremental solution is reported in **Figure 16** for a low hardening exponent, *N* = 0.01, as a contour plot (left) and as a graph (along the *x*1-axis, right) of the incremental shear deformation γ˙<sup>12</sup> (divided by the applied remote shear γ˙<sup>∞</sup> <sup>12</sup> ). Note that, similarly to the crack tip fields in fracture mechanics, the incremental stress and deformation display square-root singularities at the tips of the pre-existing shear band. Evaluation of the solution obtained from the perturbative approach analytically confirms the conclusions drawn from the imperfection approach (see the numerical simulations reported in **Figures 9** and **13**), in particular:


We finally remark that, although the tendency toward rectilinear propagation of a shear band has been substantiated through the use of a von Mises plastic material, substantial changes are not expected when a different yield criterion (for instance, pressure-sensitive as Drucker–Prager) is employed.

<sup>8</sup>Note that the notation used here differs from that adopted in Bigoni and Dal Corso (2008), where the principal axes are denoted by *x*<sup>1</sup> and *x*<sup>2</sup> and the system inclined at 45° is denoted by *x*ˆ<sup>1</sup> and *x*ˆ2.

# **5. CONCLUSION**

Two models of shear band have been described, one in which the shear band is an imperfection embedded in a material and another in which the shear band is a perturbation, which emerges during a homogeneous deformation of an infinite material. These two models explain how shear bands tend toward a rectilinear propagation under continuous shear loading, a feature not observed for fracture trajectories in brittle materials. This result can be stated in different words pointing out that, while crack propagation occurs following a maximum tensile stress criterion, a shear band grows according to a maximum Mises stress, a behavior representing a basic micromechanism of failure for ductile materials. The developed models show also a strong stress concentration at the shear band tip, which strongly concur to shear band growth.

### **ACKNOWLEDGMENTS**

DB, NB, and FDC gratefully acknowledge financial support from the ERC Advanced Grant "Instabilities and non-local multiscale modeling of materials" FP7-PEOPLE-IDEAS-ERC-2013- AdG (2014-2019). AP thanks financial support from the FP7- PEOPLE-2013-CIG grant PCIG13-GA-2013-618375-MeMic.

#### **REFERENCES**


Hutchinson, J. W., and Tvergaard, V. (1981). Shear band formation in plane strain. *Int. J. Solids Struct.* 17, 451–470. doi:10.1111/j.1365-2818.2009.03250.x


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 23 January 2015; accepted: 01 March 2015; published online: 23 March 2015. Citation: Bordignon N, Piccolroaz A, Dal Corso F and Bigoni D (2015) Strain localization and shear band propagation in ductile materials. Front. Mater. 2:22. doi: 10.3389/fmats.2015.00022*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials.*

*Copyright © 2015 Bordignon, Piccolroaz, Dal Corso and Bigoni. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Graphene and black holes: novel materials to reach the unreachable

# **Alfredo Iorio\***

Faculty of Mathematics and Physics, Charles University, Prague, Czech Republic \*Correspondence: alfredo.iorio@mff.cuni.cz

**Edited by:** Simone Taioli, Fondazione Bruno Kessler, Italy

#### **Reviewed by:**

Alessio Paris, Fondazione Bruno Kessler, Italy

**Keywords: graphene, quantum field theory, black holes, hawking radiation, quantum gravity**

Physics is an experimental science. Nonetheless, in the last decades, it proved very difficult (if not impossible) to reconcile theoretical investigations of the fundamental laws of Nature with the necessary experimental tests. The divarication between theory and experiments is the central problem of contemporary fundamental physics that, to date, is still unable to furnish a consistent quantum theory of gravity, or to obtain experimental evidences of milestones theories like Supersymmetry [see, e.g., the whole Issue (de Haro et al., 2013) on the status of String Theory].

A widespread view of this problem is that experimental observations of these types of phenomena can only be achieved at energies out of the reach of our laboratories (e.g., the Planck energy 10<sup>19</sup> GeV vs. reachable energies in our laboratories of 10<sup>3</sup> GeV). In our view, due to unprecedented behaviors of certain novel materials, nowadays indirect tests should be considered as a viable alternative to direct observations.

This field of research is not a novelty. It usually goes under the generic name of "analog gravity" (Novello et al., 2002; Volovik, 2003; Barceló et al., 2011), although, perhaps, to call it "bottom-up" approach does it more justice: analog experimental settings on the bottom, fundamental theories of Nature on the top. Nonetheless, for various reasons, it has been seen as little more than a curiosity. It is seen as an amusing and mysterious series of coincidences that (to our knowledge) never are taken as tests of the aspects of the fundamental theories they reproduce. Neither is taken as experimental tests of the fundamental theories, the mysterious and

amusing coincidences of the "top–bottom" approach of the AdS/CFT correspondence [where theoretical constructions of the fundamental world are used to describe experimental results at our energy scale (Maldacena, 2012)].

As well known, with graphene (Novoselov et al.,2004), we have a quantum relativistic-like Dirac massless field available on a nearly perfectly two-dimensional sheet of carbon atoms [see Castro Neto et al. (2009) for a review]. Recent work shows the emergence of gravity-like phenomena on graphene (Iorio, 2011, 2012, 2013; Iorio and Lambiase, 2012, 2014). More precisely, the Hawking effect can take place on graphene membranes shaped as Beltrami pseudospheres of suitably large size, hence even (2 + 1)-dimensional black-hole scenarios (Bañados et al., 1992) are in sight. The Hawking effect here manifests itself through a finite temperature electronic local density of states [for reviews, see Iorio (2014, 2015)].

The predictions of the Hawking effect on graphene are based on the possibility to obtain very specific shapes, e.g., the Beltrami pseudosphere that should recreate, for the pseudoparticles of graphene, conditions related to those of a spacetime with an horizon. What we are first trying (Gabbrielli et al., under preparation) is to obtain a clear picture of what happens to *N* classical particles, interacting via a simple potential, e.g., a Lennard–Jones potential, and constrained on the Beltrami. This will furnish important pieces of information on the actual structure of the membrane. In fact, is well known from similar work with the sphere [that goes under the name of "generalized Thomson problem," see, e.g., Bowick et al. (2000)] that defects will form

more and more, and their spatial arrangements are highly non-trivial, and follow patterns related to the spontaneous breaking of the appropriate symmetry group (Iorio and Sen, 2006). Once the coordinates of the *N* points are found in this way, we need to simulate the behaviors of Carbon atoms arranged in that fashion, hence, we essentially change the potential to the appropriate one, and perform a Density Functional type of computation. The number of atoms that we can describe this way is of the order of 10<sup>3</sup> , highly demanding computer-time wise, but still too small for the reaching of the horizon. Nonetheless, the results obtained will be important to refine various details of the theory. We need to go further, toward a big radius of curvature *r*, when the Hawking effect should be visible.

Any serious attempt to understand Quantum Gravity has to start from the Hawking effect. That is why black holes are at the crossroad of many of the speculations about the physics at the Planck scale. From (Iorio, 2011, 2012, 2013; Iorio and Lambiase, 2012, 2014) a goal that seems in sight is the realization of reliable set-ups, where graphene well reproduces the blackhole thermodynamics scenarios, with the analog gravity of the appropriate kind to emerge from the description of graphene's membrane. The lattice structure, the possibility to move through energy regimes where discrete and continuum descriptions coexist, and the unique features of matter fields whose relativistic structure is induced by the spacetime itself, are all issues related to Quantum Gravity (Loll, 1998; 't Hooft, 2009) that can be explored with graphene.

Many other tantalizing fundamental questions can be addressed with graphene. To mention only two: there are results of Alvarez et al. (2012, 2014) that point toward the use of graphene to have alternative realizations of Supersymmetry, and there are models of the Early Universe, based on (2 + 1)-dimensional gravity (van der Bij, 2007), where graphene might also play a role.

We are lucky that these "wonders" are predicted to be happening on a material that is, in its own right, enormously interesting for applications. Hence, there is expertise worldwide on how to manage a variety of cases. Nonetheless, the standard agenda of a material scientist is of a different kind than testing fundamental laws of Nature. Therefore, let us conclude by invoking the necessity of a dedicated laboratory, where condensed matter and other low energy systems are experimentally studied with the primary goal of reproducing phenomena of the fundamental kind.

For the reasons outlined above, graphene is a very promising material for this purpose.

#### **ACKNOWLEDGMENTS**

The author warmly thanks the LISC laboratories of Fondazione Bruno Kessler and ECT\*, Trento, Italy, for the kind hospitality. Funding: the author acknowledges the Czech Science Foundation (GACR), ˇ Contract No. 14-07983S, for support.

#### **REFERENCES**

de Haro, S., Dieks, D., 't Hooft, G., and Verlinde, E. (eds). (2013). Forty years of string theory: reflecting on the foundations. *Found. Phys.* 43(1). Available at: http://link.springer.com/ journal/10701/43/1/page/1


**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 04 December 2014; accepted: 18 December 2014; published online: 12 January 2015.*

*Citation: Iorio A (2015) Graphene and black holes: novel materials to reach the unreachable. Front. Mater. 1:36. doi: 10.3389/fmats.2014.00036*

*This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials.*

*Copyright © 2015 Iorio. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# ADVANTAGES OF PUBLISHING IN FRONTIERS

FAST PUBLICATION Average 90 days from submission to publication

#### COLLABORATIVE PEER-REVIEW

Designed to be rigorous – yet also collaborative, fair and constructive

RESEARCH NETWORK Our network increases readership for your article

# OPEN ACCESS

Articles are free to read, for greatest visibility

# TRANSPARENT

Editors and reviewers acknowledged by name on published articles

### GLOBAL SPREAD Six million monthly page views worldwide

### COPYRIGHT TO AUTHORS

No limit to article distribution and re-use

IMPACT METRICS Advanced metrics track your article's impact

# SUPPORT By our Swiss-based editorial team

EPFL Innovation Park · Building I · 1015 Lausanne · Switzerland T +41 21 510 17 00 · info@frontiersin.org · frontiersin.org