cazeaux@vt.edu
Office: McBryde Room 570
(540) 231-3446

Paul Cazeaux
460 McBryde Hall
225 Stanger Street
Virginia Tech
Blacksburg, VA 24061-1026, USA

cazeaux@vt.edu
Office: McBryde Room 570
(540) 231-3446

Paul Cazeaux
460 McBryde Hall
225 Stanger Street
Virginia Tech
Blacksburg, VA 24061-1026, USA

List of publications

Preprints

Non-Abelian topological defects and strain mapping in 2D moiré materials.
(with R. Engelke, H. Yoo, S. Carr, K. Xu, R. Allen, A.M. Valdivia, M. Luskin, E. Kaxiras, M. Kim, J.H. Han and P. Kim)
arXiv:2207.05276 (15 pages)

Abstract. We present a general method to analyze the topological nature of the domain boundary connectivity that appeared in relaxed moiré superlattice patterns at the interface of 2-dimensional (2D) van der Waals (vdW) materials. At large enough moiré lengths, all moiré systems relax into commensurated 2D domains separated by networks of dislocation lines. The nodes of the 2D dislocation line network can be considered as vortex-like topological defects. We find that a simple analogy to common topological systems with an S1 order parameter, such as a superconductor or planar ferromagnet, cannot correctly capture the topological nature of these defects. For example, in twisted bilayer graphene, the order parameter space for the relaxed moiré system is homotopy equivalent to a punctured torus. Here, the nodes of the 2D dislocation network can be characterized as elements of the fundamental group of the punctured torus, the free group on two generators, endowing these network nodes with non-Abelian properties. Extending this analysis to consider moiré patterns generated from any relative strain, we find that antivortices occur in the presence of anisotropic heterostrain, such as shear or anisotropic expansion, while arrays of vortices appear under twist or isotropic expansion between vdW materials. Experimentally, utilizing the dark field imaging capability of transmission electron microscopy (TEM), we demonstrate the existence of vortex and antivortex pair formation in a moiré system, caused by competition between different types of heterostrains in the vdW interfaces. We also present a methodology for mapping the underlying heterostrain of a moiré structure from experimental TEM data, which provides a quantitative relation between the various components of heterostrain and vortex-antivortex density in moiré systems.

arXiv:2206.05580 (36 pages)

Abstract. Twisted bilayer graphene gives rise to large moiré patterns that form a triangular network upon mechanical relaxation. If gating is included, each triangular region has gapped electronic Dirac points that behave as bulk topological insulators with topological indices depending on valley index and the type of stacking. Since each triangle has two oppositely charged valleys, they remain topologically trivial. In this work, we address several questions related to the edge currents of this system by analysis and computation of continuum PDE models. Firstly, we derive the bulk invariants corresponding to a single valley, and then apply a bulk-interface correspondence to quantify asymmetric transport along the interface. Secondly, we introduce a valley-coupled continuum model to show how valleys are approximately decoupled in the presence of small perturbations using a multiscale expansion, and how valleys couple for larger defects. Thirdly, we present a method to prove for a large class of continuum (pseudo-)differential models that a quantized asymmetric current is preserved through a junction such as a triangular network vertex. We support all of these arguments with numerical simulations using spectral methods to compute relevant currents and wavepacket propagation.

Relaxation and domain wall structure of bilayer moire systems.
(with D. Clark, R. Engelke, P. Kim and M. Luskin.)
arXiv:2211.12274 (36 pages)

Abstract. Moire patterns result from setting a 2D material such as graphene on another 2D material with a small twist angle or from the lattice mismatch of 2D heterostructures. We present a continuum model for the elastic energy of these bilayer moire structures that includes an intralayer elastic energy and an interlayer misfit energy that is minimized at two stackings (disregistries). We show by theory and computation that the displacement field that minimizes the global elastic energy subject to a global boundary constraint gives large alternating regions of one of the two energy-minimizing stackings separated by domain walls. We derive a model for the domain wall structure from the continuum bilayer energy and give a rigorous asymptotic estimate for the structure. We also give an improved estimate for the L2-norm of the gradient on the moire unit cell for twisted bilayers that scales at most inversely linearly with the twist angle, a result which is consistent with the formation of one-dimensional domain walls with a fixed width around triangular domains at very small twist angles.

Published or accepted papers

Randomized Algorithms for Rounding in the Tensor-Train Format. SISC
(with H. Al Daas, G. Ballard, E. Hallman, A. Międlar, M. Pasha, T.W. Reid, and A.K. Saibaba)
SIAM Journal on Scientific Computing, to appear. doi.org/10.1137/21M1451191 (22 pages)

Abstract. The tensor-train (TT) format is a highly compact low-rank representation for high-dimensional tensors. TT is particularly useful when representing approximations to the solutions of certain types of parametrized partial differential equations. For many of these problems, computing the solution explicitly would require an infeasible amount of memory and computational time. While the TT format makes these problems tractable, iterative techniques for solving the PDEs must be adapted to perform arithmetic while maintaining the implicit structure. The fundamental operation used to maintain feasible memory and computational time is called rounding, which truncates the internal ranks of a tensor already in TT format. We propose several randomized algorithms for this task that are generalizations of randomized low-rank matrix approximation algorithms and provide significant reduction in computation compared to deterministic TT-rounding algorithms. Randomization is particularly effective in the case of rounding a sum of TT-tensors (where we observe 20x speedup), which is the bottleneck computation in the adaptation of GMRES to vectors in TT format. We present the randomized algorithms and compare their empirical accuracy and computational time with deterministic alternatives.

Physical Review B, 101 (224107), 2020. doi.org/10.1103/PhysRevB.101.224107 (14 pages)

Abstract. The incommensurate stacking of multilayered two-dimensional materials is a challenging problem from a theoretical perspective and an intriguing avenue for manipulating their physical properties. Here we present a multiscale model to obtain the mechanical relaxation pattern of twisted trilayer van der Waals (vdW) heterostructures with two independent twist angles, a generally incommensurate system without a supercell description. We adopt the configuration space as a natural description of such incommensurate layered materials, based on the local environment of atomic positions, bypassing the need for commensurate approximations. To obtain the relaxation pattern, we perform energy minimization with respect to the relaxation displacement vectors. We use a continuum model in combination with the generalized stacking fault energy to describe the interlayer coupling, obtained from first-principles calculations based on density functional theory. We show that the relaxation patterns of twisted trilayer graphene and WSe2 are “moiré of moiré,” as a result of the incommensurate coupling two bilayer moiré patterns. We also show that, in contrast to the symmetry-preserving in-plane relaxation in twisted bilayers, trilayer relaxation can break the two fold rotational symmetry about the xy plane when the two twist angles are equal.

Energy minimization of 2D incommensurate heterostructures. Arch. Rational Mech. Anal.
(with M. Luskin, D. Massatt)
Archive for Rational Mechnics and Analysis, 235, 2020, pp 1289-1325. doi.org/10.1007/s00205-019-01444-y (37 pages)

Abstract. We derive and analyze a novel approach for modeling and computing the mechanical relaxation of incommensurate 2D heterostructures. Our approach parametrizes the relaxation pattern by the compact local configuration space rather than real space, thus bypassing the need for the standard supercell approximation and giving a true aperiodic atomistic configuration. Our model extends the computationally accessible regime of weakly coupled bilayers with similar orientations or lattice spacing, for example materials with a small relative twist where the widely studied large-scale moire patterns arise. Our model also makes possible the simulation of multi-layers for which no interlayer empirical atomistic potential exists, such as those composed of MoS2 layers, and more generally makes possible the simulation of the relaxation of multi-layer heterostructures for which a planar moire pattern does not exist.

Atomic reconstruction at van der Waals interface in twisted bilayer graphene. Nature Materials
(with H. Yoo, R. Engelke, S. Carr, S. Fang, K. Zhang, S.H. Sung, R. Hovden, A.W. Tsen, T. Taniguchi, K. Watanabe, G-C. Yi, M. Kim, M. Luskin, E.B. Tadmor, E. Kaxiras and P. Kim)
Nature Materials, 18, 2019, pp 448-453. doi.org/10.1038/s41563-019-0346-z (6 pages)

Abstract. Interfaces between crystalline materials have been an essential engineering platform for modern electronics. At the interfaces in two-dimensional (2D) van der Waals (vdW) heterostructures, the twist-tunability offered by vdW crystals allows the construction of a quasiperiodic moir\'e superlattice of tunable length scale, leading to unprecedented access to exotic physical phenomena. However, these interfaces exhibit more intriguing structures than the simple moir\'e pattern. The vdW interaction that favors interlayer commensurability competes against the intralayer elastic lattice distortion, causing interfacial reconstruction with significant modification to the electronic structure. Here we demonstrate engineered atomic-scale reconstruction at the vdW interface between two graphene layers by controlling the twist angle. Employing transmission electron microscopy (TEM), we find local commensuration of Bernal stacked graphene within each domain, separated by incommensurate structural solitons. We observe electronic transport along the triangular network of one-dimensional (1D) topological channels as the electronic bands in the alternating domains are gapped out by a transverse electric field. The atomic scale reconstruction in a twisted vdW interface further enables engineering 2D heterostructures with continuous tunability.

Relaxation and Domain Formation in Incommensurate 2D Heterostructures. Phys. Rev. B
(S. Carr, D. Massatt, S.B. Torrisi, P.C, M. Luskin, E. Kaxiras)
Physical Review B, 98 224102 (2018). doi.org/10.1103/PhysRevB.98.224102 (12 pages)

Abstract. We introduce configuration space as a natural representation for calculating the mechanical relaxation patterns of incommensurate two-dimensional (2D) bilayers, bypassing supercell approximations to encompass aperiodic relaxation patterns. The approach can be applied to a wide variety of 2D materials through the use of a continuum model in combination with a generalized stacking fault energy for interlayer interactions. We present computational results for small-angle twisted bilayer graphene and molybdenum disulfide (MoS2), a representative material of the transition metal dichalcogenide (TMDC) family of 2D semiconductors. We calculate accurate relaxations for MoS2 even at small twist-angle values, enabled by the fact that our approach does not rely on empirical atomistic potentials for interlayer coupling. The results demonstrate the efficiency of the configuration space method by computing relaxations with minimal computational cost for twist angles down to 0.05º, which is smaller than what can be explored by any available real space techniques. We also outline a general explanation of domain formation in 2D bilayers with nearly-aligned lattices, taking advantage of the relationship between real space and configuration space.

Compression of Wannier functions into Gaussian-type orbitals, Comp. Phys. Comm.
(with A. Bakhta, E. Cancès, S. Fang, E. Kaxiras)
Computer Physics Communications, 230, 2018, pp 27-37. doi.org/10.1016/j.cpc.2018.04.011 (12 pages)

Abstract. We propose a greedy algorithm for the compression of Wannier functions into Gaussian-polynomials orbitals. The so-obtained compressed Wannier functions can be stored in a very compact form, and can be used to efficiently parameterize effective tight-binding Hamiltonians for multilayer 2D materials for instance. The compression method preserves the symmetries (if any) of the original Wannier function. We provide algorithmic details, and illustrate the performance of our implementation on several examples, including graphene, hexagonal boron-nitride, single-layer FeSe, and bulk silicon in the diamond cubic structure.

Quantum plasmons with optical-range frequencies in doped few-layer graphene, Phys. Rev. B
(S.N. Shirodkar, M. Mattheakis, P.C, P. Narang, M. Soljačić, E. Kaxiras)
Physical Review B. 97, 195435 (2018). doi.org/10.1103/PhysRevB.97.195435 (6 pages)

Abstract. Although plasmon modes exist in doped graphene, the limited range of doping achieved by gating restricts the plasmon frequencies to a range that does not include the visible and infrared. Here we show, through the use of first-principles calculations, that the high levels of doping achieved by lithium intercalation in bilayer and trilayer graphene shift the plasmon frequencies into the visible range. To obtain physically meaningful results, we introduce a correction of the effect of plasmon interaction across the vacuum separating periodic images of the doped graphene layers, consisting of transparent boundary conditions in the direction perpendicular to the layers; this represents a significant improvement over the exact Coulomb cutoff technique employed in earlier works. The resulting plasmon modes are due to local field effects and the nonlocal response of the material to external electromagnetic fields, requiring a fully quantum mechanical treatment. We describe the features of these quantum plasmons, including the dispersion relation, losses, and field localization. Our findings point to a strategy for fine-tuning the plasmon frequencies in graphene and other two-dimensional materials.

ESAIM: Mathematical Modelling and Numerical Analysis 52(2), 2018, pp 729 - 749. doi.org/10.1051/m2an/2017057 (21 pages)

Abstract. The recent fabrication of weakly interacting incommensurate two-dimensional lattices requires an extension of the classical notion of the Cauchy-Born strain energy density since these atomistic systems are not periodic. In this paper, we rigorously formulate and analyze a Cauchy-Born strain energy density for weakly interacting incommensurate one-dimensional lattices (chains) as a large body limit and we give error estimates for its approximation by the popular supercell method.

Twistronics: Manipulating the Electronic Properties of Two-dimensional Layered Structures through their Twist Angle, Phys. Rev. B
(S. Carr, D. Massatt, S. Fang, P.C., M. Luskin and E. Kaxiras).
Physical Review B 95, 075420 (2017). doi:10.1103/PhysRevB.95.075420 (6 pages)

Abstract. The ability in experiments to control the relative twist angle between successive layers in two- dimensional (2D) materials offers a new approach to manipulating their electronic properties; we refer to this approach as "twistronics". A major challenge to theory is that, for arbitrary twist angles, the resulting structure involves incommensurate (aperiodic) 2D lattices. Here, we present a general method for the calculation of the electronic density of states of aperiodic 2D layered materials, using parameter-free hamiltonians derived from ab initio density-functional theory. We use graphene, a semimetal, and MoS2, a representative of the transition metal dichalcogenide (TMDC) family of 2D semiconductors, to illustrate the application of our method, which enables fast and efficient simulation of multi-layered stacks in the presence of local disorder and external fields. We comment on the interesting features of their Density of States (DoS) as a function of twist-angle and local configuration and on how these features can be experimentally observed.

Journal of Mathematical Physics 58, 063502 (2017). doi:10.1063/1.4984041 (23 pages)

Abstract. We give an exact formulation for the transport coefficients of incommensurate two-dimensional atomic multilayer systems in the tight-binding approximation. This formulation is based upon the C*-algebra framework introduced by Bellissard and collaborators to study aperiodic solids (disordered crystals, quasicrystals, and amorphous materials), notably in the presence of magnetic fields (quantum Hall effect). We also present numerical approximations and test our methods on a one-dimensional incommensurate bilayer system.

Analysis of rippling in incommensurate one-dimensional coupled chains, Multiscale. Model. and Simul.
(with M. Luskin and E. Tadmor)
Multiscale Modeling and Simulation 15(1), 2017, pp. 56-73. doi:10.1137/16M1076198 (18 pages)

Abstract. Graphene and other recently developed 2D materials exhibit exceptionally strong in-plane stiffness. Relaxation of few-layer structures, either free-standing or on slightly mismatched substrates occurs mostly through out-of-plane bending and the creation of large-scale ripples. In this work, we present a novel double chain model, where we allow relaxation to occur by bending of the incommensurate coupled system of chains. As we will see, this model can be seen as a new application of the well-known Frenkel-Kontorova model for a one-dimensional atomic chain lying in a periodic potential. We focus in particular on modeling and analyzing ripples occurring in ground state configurations, as well as their numerical simulation.

Perturbation theory for weakly coupled two-dimensional layers, J. Mater. Res.
(G.Tritsaris, S. Shirodkar, T. Kaxiras, P. C., M. Luskin, P. Plechac and E. Cancès.)
Journal of Material Research 31(7), 2016, pp. 959-966. doi:10.1557/jmr.2016.99 (8 pages)

Abstract. A key issue in two-dimensional structures composed of atom-thick sheets of electronic materials is the dependence of the properties of the combined system on the features of its parts. Here, we introduce a simple framework for the study of the electronic structure of layered assemblies based on perturbation theory. Within this framework, we calculate the band structure of commensurate and twisted bilayers of graphene (Gr) and hexagonal boron nitride (h-BN), and of a Gr/h-BN heterostructure, which we compare with reference full-scale density functional theory calculations. This study presents a general methodology for computationally efficient calculations of two-dimensional materials and also demonstrates that for relatively large twist in the graphene bilayer, the perturbation of electronic states near the Fermi level is negligible.

Projective multiscale time-integration for electrostatic particle-in-cell methods, Accepted for publication in Comp. Phys. Comm.
(with J. Hesthaven)
Accepted for publication in Computer Physics Communications. doi.org/10.1016/j.cpc.2018.10.012 (26 pages)

Abstract. The simulation of problems in kinetic plasma physics are often challenging due to strongly coupled phenomena across multiple scales. In this work, we propose a wavelet-based coarse-grained numerical scheme, based on the framework of Equation-Free Projective Integration, for a kinetic plasma system modeled by the Vlasov-Poisson equations. A kinetic particle-in-cell (PIC) code is used to simulate the meso scale dynamics for short time intervals. This allows the extrapolation over long time-steps of the behavior of a coarse wavelet-based discretization of the system. To validate the approach and the underlying concepts, we perform two 1D1V numerical experiments: nonlinear propagation and steepening of an ion wave, and the expansion of a plasma slab in vacuum. The direct comparisons to resolved PIC simulations show good agreement. We show that the speedup of the projective integration scheme over the full particle scheme scales linearly with the system size, demonstrating efficiency while taking into account fully kinetic, non-Maxwellian effects. This suggests that the approach is potentially interesting for kinetic plasma problems with a large separation of scales.

Mathematical Models and Methods in Applied Sciences 25(6), 2015, pp 1125-1177. doi:10.1142/S0218202515500293 (53 pages)

Abstract. We are interested in the mathematical modeling of the deformation of the human lung tissue, called the lung parenchyma, during the respiration process. The parenchyma is a foam-like elastic material containing millions of air-filled alveoli connected by a tree-shaped network of airways. In this study, the parenchyma is governed by the linearized elasticity equations and the air movement in the tree by the Poiseuille law in each airway. The geometric arrangement of the alveoli is assumed to be periodic with a small period ε > 0. We use the two-scale convergence theory to study the asymptotic behavior as ε goes to zero. The effect of the network of airways is described by a nonlocal operator and we propose a simple geometrical setting for which we show that this operator converges as ε goes to zero. We identify in the limit the equations modeling the homogenized behavior under an abstract convergence condition on this nonlocal operator. We derive some mechanical properties of the limit material by studying the homogenized equations: the limit model is nonlocal both in space and time if the parenchyma material is considered compressible, but only in space if it is incompressible. Finally, we propose a numerical method to solve the homogenized equations and we study numerically a few properties of the homogenized parenchyma model.

ESAIM: Proceedings and Surveys, 48, 2015, pp 156-168. doi:10.1051/proc/201448006 13 pages)

Abstract. Our work is motivated by numerical homogenization of materials such as concrete, modeled as composites structured as randomly distributed inclusions imbedded in a matrix. In this paper, we propose a method for the approximation of the periodic corrector problem based on boundary integral equations. The fully populated matrices obtained by the discretization of the integral operators are successfully dealt with using the ℋ-matrix format.

Homogenization of a model for the propagation of sound in the lungs, Multiscale. Model. Simul.
(with C. Grandmont and Y. Maday)
Multiscale Modeling and Simulation 13(1), 2015, pp 43-71. doi:10.1137/130916576 (29 pages)

Abstract. In this paper, we are interested in the mathematical modeling of the propagation of sound waves in the lung parenchyma, which is a foam-like elastic material containing millions of air-filled alveoli. In this study, the parenchyma is governed by the linearized elasticity equations, and the air by the acoustic wave equations. The geometric arrangement of the alveoli is assumed to be periodic with a small period $\varepsilon>0$. We consider the time-harmonic regime forced by vibrations induced by volumic forces. We use the two-scale convergence theory to study the asymptotic behavior as $\varepsilon$ goes to zero and prove the convergence of the solutions of the coupled fluid-structure problem to the solution of a linear-elasticity boundary value problem.

Mathematical Modelling and Numerical Analysis, 48(1), 2013, pp 27-52. doi:10.1051/m2an/2013093 (23 pages)

Abstract. In this paper we develop and study numerically a model to describe some aspects of sound propagation in the human lung, considered as a deformable and viscoelastic porous medium (the parenchyma) with millions of alveoli filled with air. Transmission of sound through the lung above 1 kHz is known to be highly frequency–dependent. We pursue the key idea that the viscoelastic parenchyma structure is highly heterogeneous on the small scale $\varepsilon$ and use two–scale homogenization techniques to derive effective acoustic equations for asymptotically small $\varepsilon$. This process turns out to introduce new memory effects. The effective material parameters are determined from the solution of frequency–dependent micro–structure cell problems. We propose a numerical approach to investigate the sound propagation in the homogenized parenchyma using a Discontinuous Galerkin formulation. Numerical examples are presented.

Thesis manuscripts

Abstract. We present macroscopic models of the mechanical behavior of the human lung’s parenchyma obtained by the two–scale homogenization method. The parenchyma is a porous material, with a huge number of air–filled alveoli connected to the exterior air by the bronchial tree. This complex microscopic structure defines the macroscopic behavior. We propose to study two problems in particular : modeling the deformation of the parenchyma while taking into account the ventilation by the bronchial tree, and the sound or ultrasound propagation through the parenchyma.
The first part focuses on the coupling between parenchyma and bronchial tree. We begin by describing a model for the parenchyma deformation. We model (i) the parenchyma as a linear elastic material, (ii) the alveoli as periodically distributed cavities in the macroscopic parenchyma domain and (iii) the bronchial tree as a dyadic resistive tree. We write the equations of the model as a coupled fluid–structure system modeling the three–dimensional parenchyma’s displacement and depending on a parameter $\varepsilon$ which corresponds to the size of the periodicity cell. We study the two–scale convergence of the solutions of this system under an abstract hypothesis that describes the convergence of the action of the tree on the parenchyma. We obtain a macroscopic description of the parenchyma as a viscoelastic material where the tree induces a spatially non–local dissipation. In this part, we also study the abstract condition we have introduced. We propose two models for the irrigation of the domain by the tree inspired by the lung’s structure and for which the abstract condition can be verified. Finally, we describe a numerical method for the macroscopic problem and we illustrate the previous work by numerical simulations in two dimensions.
The second part focuses on the sound wave propagation in the parenchyma. We do not take into account the effect of the bronchial tree in this case. We homogenize in the frequency domain a first model coupling the linearized elasticity equations in the parenchyma and the acoustic equation in the air. We rigorously obtain the Rice model which describes sound propagation at low frequencies. We encounter a difficulty because the problem we investigate, of Helmholtz type, is not well–posed for all values of the frequency. To show the result, we use an argument by contradiction based on the Fredholm alternative. Then, we homogenize a second model which takes into account the viscoelastic and heterogeneous nature of the parenchyma at the microscopic level. The macroscopic viscoelastic coefficients depend on frequency. The material exhibits some new memory effects compared to its individual components. We propose a numerical method based on discontinuous Galerkin finite elements to solve the homogenized problem we have obtained. The numerical results obtained in a two–dimensional test case show that this model enables us to recover some physiological observations on the propagation of low–frequency ultrasound.

Autour de la modélisation du poumon. (in French) Master's thesis.
Abstract. Mon stage de Master, et ma thèse, ont pour objet la description en termes mathématique du poumon. La modélisation mathématique et numérique des divers aspects du corps humain sont aujourd’hui de plus en plus sollicités par le monde médical. Ces modélisations vont de l’étude de l’évolution de certaines molécules, de bactéries, à la modélisation du comportement d’organes du corps humain. Ainsi, l’analyse de modèles couplés fluide - structure a pu conduire à la modélisation de flux sanguins dans les artères, autour d’anévrismes par exemple, ouà la modélisation du coeur. Je me suis plus particulièrement intéressé à la modélisation du poumon, qui vise à une meilleure compréhension des mécanismes susceptibles d’être rencontrés par les médecins lors des mesures effectuées en milieu hospitalier, par exemple dans les cas de patients asthmatiques.
Ces problèmes ont donné lieu à un certain nombre de modèles simplifiés, visant à prendre en compte tel ou tel aspect, ou tel ou tel paramètre, jusqu'à des modèles tridimensionnels en géométrie réelle. Pour les expliquer, nous allons d’abord présenter rapidement la physiologie du poumon. Le modèle expérimental décrit est celui présenté par Weibel [1] .