My current main research interest concerns developing a mathematically rigorous theory of fracture at the atomistic scale. The term fracture refers to a separation of a material into two or more pieces by means of a thin opening, a crack, propagating through it due to the accumulation of stress. Cracks are abundant in humanbuilt structures and are, for instance, a common sight in the walls of ageing buildings, as seen to the right in the a photo of a crack in the wall of my student accomodation at Warwick.
The atomistic scale is known to lie at the heart of fracture, with cracks propagating via atomistic mechanisms involving breaking interaction bonds between individual atoms. If one zooms in on a crack tip, firstly it leads to a simplified geometry and then the atomistic structure is revealed, as illustrated below.
Atomistic brittle fracture
In brittle fracture, with no plastic deformation occurring around the crack tip, the fundamental insight of the standard continuum fracture mechanics lies in identifying a single parameter \(K\), known as the stress intensity factor and showing that if it exceeds a materialdependent critical value, \(K_{\rm c}\) then it is energetically favourable for the crack to propagate.
The key atomistic phenomenon not captured in the continuum theory is known as lattice trapping, a term referring to there existing a range of values for the stress intensity factor \((K_,K_+)\) for which the crack remains locally stable, i.e. latticetrapped, thus conceptually invalidating the key insight of the continuum theory. The critical stress intensity factor \(K_{\rm c}\) computed directly from finitesize atomistic models lies somewhere inbetween and is known to deviate from the continuum theory prediction.
My work in this area is at the junction of mathematical analysis of the underlying models and the development of novel algorithms leading to quantifiable improvements in the accuracy of computer simulations.
The basic premise is to consider deformations of some crystalline lattice \(\Lambda \subset \mathbb{R}^n\) of the form $$ y\,:\Lambda \to \mathbb{R}^k,\quad y(m) = K\hat{u}(m\alpha) + u(m), $$ where \(K \in \mathbb{R}\) is the stress intensity factor of the mode of fracture considered, \(\hat{u}\) is the continuum theory predictor , known to give rise to a strain field \(D\hat{u}(m) \sim m^{1/2}\), \(\alpha \in \mathbb{R}^n\) is the crack tip position and \(u\) is the atomistic correction, which is constrained to belong to function space \(\dot{\mathcal{H}^1}\), which is a discrete variant of the canonical Sobolev space \(\dot{H}^1(\R^n)\). The key achievement of [2] is proving for certain models that at an equilibrium the atomistic correction strain field decays \(Du(m)~\sim~m^{3/2}\), as confirmed in numerical tests shown to the right.
A finite domain approximation restricting $$ u(m) = 0\;\forall m \not\in\Omega_R := \Lambda \cap B_R, $$ gives rise to the analysis of cell size effects. Under the reasonable assumption that there exists a continuous path of equilibrium configurations in the infinite problem \(s~\mapsto~(u_s,K_s,\alpha_s)\), it is proven in [3] that for \(R\) large enough there exists an approximate path of equlibrium configurations \(s~\mapsto~(u^R_s,K^R_s,\alpha^R_s)\) such that $$ \Du_sDu^R_s\_{\ell^2} + K_s  K^R_s + \alpha_s  \alpha_s^R \lesssim R^{1/2}. $$ While the convergence rate follows directly from the decay properties of \(Du\), a recourse to bifurcation theory for Banach spaces is needed for this to to apply to continuous paths of equilibria.
The work in [5] looks at practical ways of computing approximate equilibrium paths. In particular we have proposed the novel numericalcontinuationenhanced flexible boundary (NCFlex) scheme, which appears to improve the convergence rate to \(\mathcal{O}(R^{1})\) and is also capable of identifying index1 saddle equilibria with ease, thus constituting an indespensible tool for the study of energy barriers, as shown to the right.
It is hoped that more accurate NCFlextype schemes can be systematically derived by marrying this approach with the idea of development of solutions first proposed in [1], which entails prescribing a more accurate continuum predictor.
The NCFlex scheme applied to two realistic models of silicon yields the paths of equilibrium configurations shown below, which reveals subtle qualitative difference at the unstable segments of the equilibrium paths that would be near impossible to capture with methods known to date.
Atomistic and mesoscale nearcracktip plasticity
In a recent independent endeavour [4], I have introduced a novel geometric and functional framework of the lattice manifold complex, in which a cracked crystalline material is described by a CW complex structure defined on the Riemann surface $$ \mathcal{S} = \lbrace (z,w) \in \Complex^2 \,  \, w^2 = z \rbrace $$ corresponding to the complex square root mapping. The essence of the construction is presented in the following figure.
Considering deformations of a cracked crystal in this framework has led to establishing existence of nearcracktip plasticity equilibria, in which a crack opening and a finite number of dislocations are present. This has unlocked a number of further studies about nearcracktip plasticity, most notably asking about a rigorous passage from the atomistic scale to the mesoscopic scale, as well the study of equilibria in the upscaled representation, which is the most relevant scale to the problem. More on that will appear here as the work progresses.
Stochastic approaches to atomistic material modelling
In my capacity as a Research Associate at Cardiff University, I am working with Dr Angela Mihai and Dr Thomas Woolley on a new informationtheoretic stochastic approach to atomistic material modelling, based around the idea of invoking Maximum Entropy Principle to infer the probability distribution of parameters specifying the interatomic potential. We first employ this new framework to study atomistic fracture from a stochastic point of view.
As an illustration, consider the simplified case of Mode III fracture under antiplane displacements \(y(m)~=~(0,0, y_3(m)),\,\forall m \in \Lambda\) with energy (morally) of the form $$ \mathcal{E}(y) = \sum_{m \in \Lambda} \sum_{\rho \in \mathcal{R}} \phi(y_3(m+\rho)  y_3(m))), $$ where \(\mathcal{R}\) is the set of interaction stencils and \(\phi\) is the pairpotential.
If \(\phi\) is of the form $$ \phi(r) = \alpha\left(1\exp(\beta r^2)\right), $$ then it is of interest to understand how an inherent uncertainty in the choice of \(\alpha\) and \(\beta\) propagates to uncertaintities in the key quantities describing atomistic crack propagation, namely \(K_{\rm c}\) and \( (K_,K_+)\) (as discussed in the section on atomistic brittle fracture above). A plot of functions \(\phi\) for different parameter choices (with mean highlighted) can be seen to the right.
The choice of parameters has to be such that the notion of lattice stability is satisfied, which, when fed into the Maximum Entropy Principle framework, implies that both parameters should follow Gamma distribution. A direct calculation possible in this simple scenario then reveals that \(K_{\rm c} = C \beta^{1/2}\) for some crystallinestructuredependent known constant \(C\), thus revealing the propagation of uncertainties, as can be seen in a histogram plot below.
The ongoing work is about formalising a general maximumentropyprinciplebased framework for uncertainty quantification in atomistic material modelling and, in the particular case of fracture, accounting for the phenomenon of lattice trapping.
Collaborators:

 Cardiff University, UK:
 Angela Mihai (postdoctoral mentor), Thomas Woolley.

 HeriotWatt University, UK:
 Heiko Gimperlein, Lucia Scardia.

 University of British Columbia, Canada:
 Christoph Ortner (PhD advisor).

 Kanazawa University, Japan:
 Masato Kimura, Patrick van Meurs.

 University of Warwick, UK:
 Julian Braun, Thomas Hudson (PhD advisor), James Kermode.

 CINaM, Marseille, France:
 Petr Grigorev.