MATHEMATICAL AND NUMERICAL METHODS
There are three hierarchical levels of understanding, corresponding to three levels of spatial detail, into which models of materials can be organized. The microscopic level is concerned with the properties and evolution of assemblages of large numbers of atoms and molecules, often arranged as a crystal. The macroscopic level deals with bulk averages across microstructure and is the domain of continuum mechanics and thermodynamics. The mesoscopic level is intermediate to the other two. Each of these levels has its own set of concepts and materials parameters. In principle, these parameters can be derived from properties of a lower level or measured experimentally. In either case, they provide an appropriate basis for the quantitative description of matter at the level in question, while omitting details of behavior important only at a lower level. A theory of materials consists in the calculation of material behavior within each level and the application of parameters calculated on one level to a more macroscopic one.
Within each level, the challenge is to predict phenomena from appropriate equations of motion or conditions of equilibrium. On the microscopic level, the dynamics of atoms and molecules are described by Newton's law or Schrödinger's equation. On the macroscopic level, materials properties are generally represented by a set of partial differential equations that express energy, mass, and momentum conservation and are formulated to represent the symmetry of the material to which they are applied. On the mesoscopic level, theoretical approaches are less well defined. Ginzburg-Landau and Langevin equations, however, are being used increasingly. In the Ginzburg-Landau approach, the free energy of the material is expanded as a power series of a quantity, called the order parameter, and its gradients. The value of the order parameter shows the local state (or phase) of the material (Landau and Lifshitz, 1969). In the Langevin approach, particle dynamics is modified by an interaction with a random force that represents an interaction, the details of which are undeterminable or nonessential (Parisi, 1988). Frequently, there is no sharp boundary between the mesoscopic and macroscopic levels, except the degree to which an application may demand further constraints in modeling.
The complete formulation of a materials science problem is often too complicated, involving the entire atomic or molecular structure and focus. Mathematical modeling is needed to simplify the problem while preserving the relevant physics; the model can then be studied via mathematical and computational tools. At each level, models typically lead to nonlinear equations; solutions of these equations can exhibit coherent structures, chaos, or complex patterns. Understanding the models and their solutions presents formidable analytical and numerical challenges that require extending existing methods, developing new approaches, and incorporating techniques from other fields.
Some mathematical research areas with materials science applications are summarized below. Sometimes, an opportunity is specific to a material, property, or phenomenon. Often a problem area is generic, spanning the hierarchy of materials modeling.
In analyzing and predicting the phases of a material and the nature of transitions between these phases, an electronic structure calculation (Callaway, 1991) is an important tool. Typically, electronic structure calculations are used to develop models of interatomic forces and to predict the zero-temperature configuration, that is, the crystal structure, of a large number of atoms and molecules. In these areas, algorithms are needed that can handle complicated materials efficiently and effectively.
Improved electronic structure methods and approximations are needed to treat systems with a large number of atoms or molecules. This is particularly true for such complicated structures as polymers, alloys, ceramics, and materials in which defects such as surfaces and dislocations are present. Improving the calculations of alloy properties will probably require the development of radically new methods because disorder is generally present in these materials. Determining microscopic forces between defects is an example of an especially fundamental and challenging calculation.
Improving many electronic structure methods requires going beyond the local density approximation to density functional theory. Density functional methods are based on a theorem that says the ground-state electronic energy is a universal functional of the electronic density, but that theorem provides no information on the nature of the functional. Mathematical analysis expresses the functional in pieces that are exactly known and lumps the unknown parts into a piece that is approximated, known as the local density approximation. This latter piece contains the dominant part of the electronic correlation energy, which is intrinsically a quantum many-body effect. Accurate computation of the electron structure of many materials requires improved approximations of this energy. How to do this has been the object of research for a number of years. Can bounds or inequalities be established for this part of the energy? and How must the assumptions be changed for different classes of systems? are natural questions to ask.
A general issue in electronic structure calculations is the effective utilization of parallel computers. The use of parallel computers for electronic structure calculations has been developing slowly because of the size of the computational codes, uncertainty in possible computational gains, and the difficulty in adjusting well-developed procedures to novel computational environments. Most parallel computing in this area has involved the joint efforts of a materials theorist, an applied mathematician, and a computer scientist.
Ab initio electronic structure calculations for clusters, which are common activities of the quantum chemist, are an example of a situation where strong coupling appears between treating the electron correlation energy properly and using parallel computers. The major difficulty with these methods, as with their counterparts used in condensed matter physics (mentioned in Chapter 2), is the need to correlate many electrons. Knowing how to effectively exploit parallel computers for such intensive calculations would be beneficial.
In recent years, electronic structure calculations and the study of quantum many-body phenomena have been attempted by Monte Carlo methods (Doll and Gubernatis, 1990; De Raedt and von der Linden, 1992). The quantum mechanics imposes significant constraints on the nature and utility of the algorithms used. In almost all cases, these algorithms possess what is called the sign problem, in which the transition probabilities needed in the Monte
Carlo step become negative (Schmidt and Kalos, 1984; Loh et al., 1990). In some algorithms, approximate procedures allow useful calculations to be performed; in others, the sign problem is so bad that the study of the models believed to be basic for understanding the materials, such as high-temperature superconductors, are essentially prohibited. The sign problem is a significant factor limiting the advance of theoretical many-body physics.
The study of phase transitions and their dynamics sits in a branch of physics called statistical mechanics. This branch plays an important and useful role in many areas of materials research. From the specification of the interactions between many particles (obtained from electronic structure calculations or from assumed models), statistical mechanics attempts to predict the macroscopic behavior associated with these particles and their interactions. While many exactly solvable models and exact results exist, most predictions are based on approximate theories, and there is always a need to develop better approximations and methods.
In quantum mechanical problems, where the above mentioned density functional theory is used, the Kac-Feynman path integral is also often used as the basis for analytic and numerical work (Feynman, 1972). In classical statistical mechanics, analogous mathematical constructs called Wiener integrals and density functionals exist. A Wiener integral represents the motion of particles by the trajectories of these particles as they undergo Brownian motion. Density functionals describe the spatial distribution of matter when the location of one or more particles is given (Feynman, 1972). These classical statistical mechanics constructs are used, for instance, in polymer research and in the theory of melting of solids and freezing of liquids. What is lacking is a complete understanding of how to perform systematic approximation of these density functionals and how to gauge the effects of approximations on predicted properties.
One of the most remarkable developments in the modern theory of statistical mechanics has been the renormalization group method. This method is a prescription for the systematic removal of irrelevant length or energy scales from a problem (Wilson, 1975). Successfully applied to many second-order phase transition problems, some electronic structure problems, and some polymer problems, its applicability has been limited by its computational needs, an incomplete understanding of its mathematical structure, and the lack of a prescription for applying its concepts to a wider spectrum of problems. The increased use of renormalization group methods in some areas has brought with it the need for more insight into dealing with asymptotic series and their resummation.
For equilibrium behavior, computer simulations have provided important insights, and, with the use of parallel computers, many of these simulations now involve tens of millions of particles. The physics of a phase transition, however, leads to difficulties in performing such simulations, difficulties that become more severe as the transition is approached. For Monte Carlo simulations of second-order phase transitions, recent algorithmic breakthroughs involving the use of Monte Carlo moves of clusters, as opposed to moves of single particles, have reduced the computing time for such simulations by orders of magnitude (Swendsen and Wang, 1987; Wolff, 1989; Hayes, 1993).
These algorithms were first developed for Ising models, the physical variables (spins) of which assume discrete values, and were later extended to other models, including the XY and sigma models, the variables of which assume continuous values. To date, however, these
algorithms fail for a class of problems characterized by ''frustrated'' interactions. In a frustrated system, the minimal energy configuration, say for a local arrangement of spins, becomes impossible to satisfy because of the constraints imposed by the minimal energy configuration of a neighboring region. Frustration leads to the existence of many local minima of the energy with values very close to the global minimum. This condition radically alters the dynamics of the system, renders its analysis difficult, and makes simulations increasingly difficult as the temperature is lowered.
A comparably challenging class of problems, involving random systems, is often called the spin-glass problem after physical systems in which low concentrations of atoms with magnetic moments are randomly positioned in a material and induce an interaction between each other that varies in sign. The problem is one of optimizing a nonlinear function of a large number of variables, which is an issue throughout materials research. Polymer problems and protein folding problems are solved by optimizing nonlinear functions (see Chapter 3). How does one know that the global minimum has been reached? and What techniques can one use to reach it efficiently? are major questions awaiting answers. One recent approach uses a random walk in the function to be minimized; from a series of extrema that are generated, the global minimum can be determined with high probability as the number of walks increases (Berg, 1993).
Theoretical work on the microscopic level has concentrated on predicting the properties of a system at or not far from equilibrium. Yet, the synthesis, processing, and use of many materials often occurs under conditions far from equilibrium. Both the analytical and the computational treatments of such problems are underdeveloped. Problem areas include the kinetics of chemical reactions. The kinetics of chemical reactions is an inherently nonlinear area and therefore presents a considerable challenge.
Molecular dynamics is increasingly used to study the kinetics of phase transitions and other nonequilibrium problems. A limitation of this method is that often the physical time simulated by these calculations is short compared to the time scale of the phenomenon of interest. Can this method be modified to extend the simulation to longer physical times? Part of the problem is that many time scales may be present in the problem. Can the longer scales be studied without computing time being taken up by the shorter scales?
Applied mathematicians have had a large impact on materials science at the continuum level. For example, applied mathematicians did the asymptotic analysis of the stress field near a crack-tip, which has yielded catalogs of figures of merit (stress intensity factors) for fracture prediction based on simple crack models (Babuška and Miller, 1984; Vasilopoulos, 1988; Sumaratna and Ting, 1986; Reitich, 1991). Surface-stability analysis of moving interfaces has led to an increased understanding of microstructural growth during solidification; see two subsections, Phase Transformations and Pattern Formation and Dendritic Growth, in Chapter 4. Finite-element and finite-difference methods are routinely used to map stress concentrations in many load problems. Much research in this area is already in progress, yet much remains to be done. In particular, improvements are needed
in both the analysis and numerics of partial differential equations, especially nonlinear ones.
Improving ways to generate adaptive grids is important in many problems associated with the growth and processing of materials, where the boundaries are often complicated. In these problems, a moving interface (marking the separation of a liquid and solid) or a free-moving boundary (marking the extent of solidification) can develop complex spatial structures and become unstable (Huang and Glicksman, 1981a, b; Langer, 1987; Gollub, 1991). These interface developments define the material's mesoscale structure, which in turn determines the material's macroscopic parameters and behavior. These problems often are associated with nonlinear diffusion processes and place great demands on finite-difference and finite-element methods because of the need to track the fronts and simultaneously to capture structure and patterns developing at increasingly reduced length scales. Improved numerical methods for such problems are needed.
A fascinating feature of many interface and growth problems is the occurrence of self-similar structures (Voorhees, 1985; Voorhees and Schaefer, 1987; Hardy et al., 1991). Analysis and computation would benefit from a clear understanding of the conditions under which self-similar structures occur and knowledge of the geometric details of those that can occur. Many materials problems involve phenomena on disparate time and spatial scales, resulting in stiff partial differential equations. To further complicate the situation, in some problems the partial differential equations change type, for example, from elliptic to hyperbolic, as time evolves (Keyfitz and Shearer, 1990). Much work on the theory and numerical solution of partial differential equations under these challenging conditions awaits attention.
Often in materials research, the scientist is interested in inferring from measurements made on a macroscopic level information about the mesoscale. For instance, acoustic measurements may lead to inferences about the location and geometry of a crack in a material, which in turn may lead to a prediction of whether the crack is one likely to initiate fracturing of the material. This and related situations are inverse problems and are inherently ill posed. Ill-posed problems are a class of problems of mathematical interest for which improved methods would significantly influence present understanding and characterization of materials. Ill-posed problems also are basic to such bread-and-butter tasks as parameter fitting of models.
The complete solution of a problem is not always essential or desired. Of interest may be only what happens at long or short times, long or short wave numbers, and so on. This simplification leads to the use of asymptotic expansion methods. Of course, once the terms in an asymptotic expansion have been produced, the next natural thing is to sum this series to extend the solution, say, from long-time behavior to the earlier behavior leading to it.
Characteristic of materials modeling on the macroscopic level is the use of constitutive relations. Familiar cases are the relations between current and electric field and between stress and strain, in which quantities such as conductivity and elastic stiffness, respectively, appear. These quantities are a way of specifying the properties and behavior of the material on the macroscopic scale. With new material types being developed and studied, extending current constitutive models and developing new constitutive models become important. A context in which this is the case is a new family of viscoelastic
materials (Renardy et al., 1987). In general, the viscoelastic response of materials, especially of polymers, needs more attention. Plasticity has been and continues to be an important area of investigation. Today, nonlinear behavior in materials, including nonlinear optical behavior, is at the forefront of investigation.
Between the microscopic and macroscopic scales is the mesoscale, with many significant materials problems. On this scale, the materials scientist is concerned primarily with characterizing and controlling the microstructure of a material in ways that allow one to predict macroscopic properties. The mesoscale structure is usually a long-lived state of the material that is far from equilibrium. At the mesoscopic level, most materials are inhomogeneous, and macroscopic descriptions usually require an average over, or an explicit treatment of, the microstructural detail.
The technique of "homogenization" of the microstructure-level partial differential equation appears useful for such mesoscale problems. Further development of homogenization procedures appears warranted. In some respects, homogenization methods build on an older body of work, effective medium approximations, that for a long time has been the subject of considerable analytical efforts. Lord Rayleigh himself pondered such problems (Rayleigh, 1892).
A material's mesostructure exhibits polycrystalline, phase-separated, and other variations in properties over regions that are too large to be treated microscopically but are sufficiently large that they have the same properties and parameters as the macroscopic bulk material. The basic question in effective medium theory is, Given the properties of the constituents, what are the properties of the aggregate? On the microscopic level, the same question is asked in the theory of alloys: Given the properties of the individual atoms or molecules, what are the properties of the combination? Not surprisingly, these two areas share approximation techniques (Gubernatis, 1978), but considerable attention has been given to the development of bounds and variational principles with regard to effective medium approximations, issues that have not been focused on as much at the microscopic level. Today, many applications require that one consider the effective nonlinear response (elastic and optical) of a material to external forces. Analysis of extreme responses and the development of statistical theories of fracture are also needed. The question of what structures produce optimal response properties is also of current interest.
As a phase transition evolves in time, phases nucleate, separate, grow, and disappear. As a material interface or free-moving surface becomes unstable, rich complex structures evolve into columnar grains, lamellar eutectics, dendritic growths of fractal structures, and other structural features. Both views of microstructural evolution, that of evolving phase transitions and that of instability development in interfaces and surfaces, point to the often exceptionally complex mesoscopic morphology of materials. See Chapters 4 and 5 for examples.
In Ginzburg-Landau (GL) theory, the free energy is expanded as a power series of a quantity called the order parameter. If the material is inhomogeneous, terms with the
gradient of the order parameter are included. The various terms in this expansion are adjusted to be consistent with the symmetry of the material. The resulting expressions can be used to investigate the phases of the material. From some microscopic models, the leading terms in this expansion can be explicitly obtained and the associated parameters identified with quantities that are calculable or measurable. This remarkable situation is true for the Bardeen-Cooper-Schrieffer (BCS) theory of superconductivity (Fetter and Walceka, 1971; see also the section Superconductivity in Chapter 4) and the theory of lattice vibrations (Maradudian and Fein, 1962; Cowley, 1964). The nature of the martensitic transformation has recently been successfully investigated in this manner (Krumhansel and Gooding, 1989; see also the Chapter 4 section, Martensite and Shape-Memory Materials). It is expected that similar methods of analysis will be fruitful in other areas of materials science in the future.
One research need presented by the GL approach is for a way to avoid redeveloping it on a case-by-case basis. Another research challenge arises when minimizing the GL energy functional to obtain the order parameter's power series structure, and therefore the structure of the material; this may present nonconvex minimization problems with no or non-unique solutions, about which more understanding is needed.
One approach to the kinetics of phase transitions is to make the order parameter time dependent and allow it to evolve in time under the action of forces taken to be proportional to gradients of the GL energy functional. These forces may be supplemented by viscous terms and random forces. There are also related approaches, expressing time-dependent position and momentum and that include viscous and random forcing terms, based on Langevin equations. The theorist is now faced with solving sets of nonlinear dissipative stochastic partial differential equations that are first-order in time. Such partial differential equations have been the subject of considerable analysis. However, the analysis needs to be made accessible, and more efficient and effective numerical algorithms need to be developed to integrate these equations.
A number of interesting materials problems exist for which a GL approach is ineffective. An example of such a situation is microphase separation in block copolymers in the strong segregation limit (see Chapter 3 sections, Block Copolymers and Interfaces in Polymer Systems). There is much work to be done in developing effective methods to treat mesoscale phenomena.
POTENTIALLY APPLICABLE MATHEMATICAL SCIENCES DEVELOPMENTS
The preceding three sections have focused on mathematical opportunities in materials science that require the extension of existing methods or the development of new approaches. In the present section, developments in the mathematical sciences are identified that might be applicable to materials problems but are not yet widely exploited.
As already mentioned, parallel computing provides important opportunities for electronic structure calculations. Such applications are only one example of the need and potential for parallel computing in materials science research. Opportunities exist for
parallel computation of material properties and behavior in nearly every area of computational materials science. From the algorithmic point of view, there are two issues. The first is that mathematicians have developed and implemented parallel versions of many basic linear algebra methods, finite-element and finite-difference methods, fast Fourier transforms, sparse-matrix methods, and so forth. Much of this off-the-shelf technology could be used in materials science calculations, but little of it is known to the materials scientist. Expanding the awareness and use of state-of-the-art parallel computation capabilities and expertise in the materials science community would be an important contribution. The second issue is that efficient use of parallel architectures generally requires rethinking the algorithm on the global scale. Directly porting a serial or vector algorithm to a parallel computer typically leads to disappointing performance gains. Many of the equations, boundary conditions, geometries, initial conditions, and so on that are germane to materials problems differ from those traditionally addressed by applied mathematicians in such applications as hydrodynamics. These issues underscore an obvious area of opportunity for interdisciplinary collaboration between materials scientists and mathematical scientists.
Wavelets (Chui, 1992), neural networks (Hertz et al., 1991), and cellular automata (Wolfram, 1983) are topics of mathematical research in which there is currently great activity. These areas are finding applications is such fields as image enhancement, signal processing, and DNA sequencing. While they are still relatively unknown to the materials community, they could have a place in materials science. Some ways in which these areas might he able to assist the materials scientist are examined below.
A natural use of wavelets would be the compression of the vast amount of data produced by experiments and simulations. One can also ask to what extent wavelet transforms (or their generalizations) could be used in place of the Fourier transform that has been used routinely in materials research for over a century. Can analysis and computations based on wavelets yield information that is hard to obtain by Fourier analysis? In electronic structure analysis, could wavelets be an attractive alternative to plane-wave and Wannier bases, which emphasize, respectively, the delocalized or localized nature of the electron? Electron density often exhibits spatial characteristics intermediate between the two (delocalized and localized).
The fields of neural networks, cellular automata, and statistical mechanics have many conceptual and model interrelationships with potential benefits for materials science. Spin-glass-like problems (Binder and Young, 1986; Fischer and Hertz, 1990) involve the optimization issues that arise in neural networks that can "learn"; an algorithm that can learn how to solve a problem while in the process of solving it may open possibilities in materials processing and synthesis. If these processes were adequately modeled, could one specify desired properties and have the network adjust the experimental controls to produce them? Can neural networks be used to determine the optimal processing conditions?
Monte Carlo simulations of Ising models and related models are examples of stochastic cellular automata. With the specification of very simple rules for going from one stage of the calculations to the next, both stochastic and deterministic automata can produce rich patterns reminiscent of mesoscale structure. To what extent can these structures be better understood by simple automata models than by the more complex and difficult-to-analyze models normally used? Statistical mechanicians already have considerable
experience with this type of problem.
A form of deterministic cellular automata, the lattice gas method for solving partial differential equations, is successful in solving hydrodynamics problems (Doolen, 1991). With lattice gas methods, simulations of fluid turbulence are available that provide direct checks of many approximate theories. Calculations of fluid flowing through a porous material, a problem important for oil recovery from sandstone, have also been successful. These methods have been applied to the study of phase-separation and chemical-reaction kinetics.
A strength of these methods is the ease with which they can handle complicated boundaries and boundary conditions. Because of the ability to use millions of cells, these methods have the potential of more easily tracking fronts and thus may be a considerable asset in solving moving interface and free-boundary problems associated with solidification and grain growth. The methods also adapt well to parallel computers; specialized computers that perform just this type of calculation could be designed.
A basic problem with lattice gas methods is determining whether the rules specified are sufficient to ensure that the calculated phenomena correspond to physical phenomena. To what level of accuracy do lattice gas methods solve the Navier-Stokes equations, the diffusion equation, and so forth? For fluid problems, the answer appears to be on relatively firm ground. This correspondence follows from analyses analogous to those used in statistical mechanics for determining macroscopic behavior in the kinetic theory of gases, from which the Navier-Stokes equations for a collection of interacting particles can be derived. For other physical systems, more analysis and understanding are needed.
Many of the concepts and methods of probability and statistics used in materials science are products of mathematical sciences research done much earlier in this century or previous centuries. Yet, probability and statistics sport many newer concepts and methods that may influence the analysis and design of materials experiments; see, for example, Chapters 9 and 10 in National Research Council (1991d). What impact might methods of nonparametric function estimation, modeling, and simulation have on the interpretation and analysis of materials experiments and simulations? Recently, Bayesian methods of statistical analysis have been used successfully in several areas of materials science, including optimization of the design of neutron scattering interferometers (Sivia et al., 1990), analysis of reflectivity data (Sivia et al., 1991), and solving ill-posed problems associated with the use of quantum Monte Carlo calculations (Gubernatis et al., 1991). The saying "With enough parameters one can fit anything" recedes into the background as these methods estimate not only what fit to the parameters is most probable but also what is the most probable number of parameters needed (Sivia and Carlile, 1992).
The ubiquitous need for nonlinear optimization has been highlighted several times in this report. Nonlinear optimization is also used in comparing theory with experiment (see, for example, the Modeling Protein Structure and Dynamics section of Chapter 3). Such optimization problems are especially challenging because of noise and the incomplete nature of the data.