Mathematical Research in Materials Science: Opportunities and Perspectives (1993)

The properties of all materials ultimately depend on their component atoms. Thus, a fundamental goal of theoretical work in materials science is the description of properties of bulk materials based only on the properties of and interactions between their constituent atoms and molecules. During the last 40 years, tremendous progress has been made toward this goal of understanding and predicting material properties from such a microscopic viewpoint. Research of this nature falls into several broad categories as follows:

• Crystalline materials are characterized by having all the constituent atoms arranged about fixed positions in a periodic lattice. Thus, given the known atomic arrangements, as obtained, say, from x-ray diffraction experiments, theory seeks both to predict and to explain experimental observations concerning the electronic, optical, magnetic, transport, and other properties of these crystalline materials.

• While a number of other materials have their constituent atoms and molecules arranged about fixed positions, these positions are not known in advance and are the desired objects of theory. Such systems involve surfaces of materials, where the surface atoms often rearrange themselves from the bulk configurations, defects in materials (such as dopants in semiconductors), where the host atoms generally alter their arrangement in response to the guest, and folded proteins, where the molecule is found in a complicated equilibrium globular structure that is related to its biological activity. Theories of these materials must first be capable of determining the equilibrium arrangements of atoms before researchers can study a host of other properties of these materials. Liquids and glasses, on the other hand, have more complicated arrangements of the constituent atoms and molecules. Liquids owe their fluidity to the fact that the constituent atoms are not centered about fixed positions but are slowly able to diffuse through the whole system. Thus, theories of liquids seek to determine the statistical properties of the equilibrium atomic and molecular configurations.

• Many systems are not in thermal equilibrium, and the object of theory becomes a desire to calculate the dynamical evolution of a nonequilibrium assembly of atoms and molecules.

The ultimate theory in all three cases would proceed completely from first principles, using quantum mechanics to describe electrons, spins, and nuclear motions. This has been the object of solid-state physics theories for crystalline materials, where the nuclei undergo only small-amplitude vibrations about equilibrium positions. However, when the nuclei are permitted to move over more substantial distances, such as in the examples described in 2 above, it is often adequate to treat the nuclear motion in terms of classical equations of motion, with the forces between atoms and molecules calculated using quantum mechanics or, when unavailable, from empirical potential functions. The latter is the realm of conventional molecular dynamics (Abraham, 1986), which has been applied to a wide variety of systems, including surface reconstruction, liquid and glass structures and dynamics, and others. In many other systems, the interatomic forces depend too intimately on nuclear

positions, necessitating the simultaneous use of quantum mechanics to compute intermolecular forces and classical mechanics to compute the nuclear motions (Car and Parrinello, 1985).

The quantum mechanical description of a bulk material is given by the Schrödinger equation, which is a partial differential equation (PDE) involving on the order of 10²³ degrees of freedom. Cohesive energies, of solids, for instance, are obtained from differences in eigenvalues of this equation for the solid and the separated atoms, but the accuracy of solution required is enormous because this cohesive energy represents only on the order of 10-4 of the total eigenvalue. Thus, mathematical approximations are necessary, and this has led to the development of a wide variety of methods, including model Schrödinger equations that contain within them a set of empirical interactions, pseudopotential methods that likewise rely on empirical information or (more recently) on accurate quantum mechanical calculations for atomic systems upon which curves/surfaces are fitted, and the nonempirical density functional and ab initio quantum mechanical methods.

Density functional theory (Hohenberg and Kohn, 1964; Kohn and Sham, 1965; Callaway and March, 1984; see also Interfaces in Polymer Systems in Chapter 3) provides a framework for the calculation of the electronic structure and total energy of any solid-state or molecular system. However, the underlying Hohenberg-Kohn theorem only proves the existence of the density functional, but does not describe its exact form. Thus, various approximate and heuristic density functionals are used in practice. Rather little has been done to compute these density functionals from first principles theory (Levy, 1991; Freed and Levy, 1982), and mathematical assistance would be welcomed (see Chapters 3 and 8) in proving theorems concerning properties of density functionals (Levy, 1991) and in developing new approximation schemes as means for systematically improving available and perforce approximate density functional methods.

When the atomic positions are assumed known, density functional methods reduce the quantum mechanical description to a nonlinear integrodifferential eigenvalue problem for an effective single-particle Hamiltonian operator whose potential energy depends on the operator's own eigenfunctions. The problem is usually solved iteratively until self-consistency is achieved between the eigenfunctions and the potential energy. The development of fast Fourier transform methods has been a key mathematical contribution to improving the efficiency of this process. The effective single-particle Hamiltonian operator's eigenvalues are identified with the electron energy levels and the eigenfunctions with the electron wave functions. The eigenfunctions and eigenvalues are then used to calculate densities of states, optical spectra, carrier life times, and so forth. Most implementations of density functional methods employ what is called the local-density approximation to incorporate exchange and correlation in the effective single-particle potential energy. In recent years, a number of methods, including quantum Monte Carlo approaches (Zhu and Louie, 1991), have been developed to go beyond this approximation, resulting in significant improvement in the calculated electron energy levels. There is growing use of electronic structure calculations based on approximate solution of the Schrödinger equation for a cluster of atoms (see Chapter 8). The main theoretical problem here lies in extrapolating solutions for small clusters to the bulk or in extracting parameters for model Hamiltonian treatments of, for instance, high-temperature superconductivity.

Density functional theory also yields the self-consistent ground-state electron density and the total energy. A calculation of the total energy as a function of small deviations from perfect crystallinity enables the computation of the vibrational spectra (phonons) of the solid. When the lowest-energy atomic arrangement is unknown in advance and is to be calculated, computation of the electronic energy for this configuration of the atoms enables determination of the total forces on each atom in the system. These forces are then used to calculate the configuration of minimum energy either by employing energy-minimization techniques or by seeking a configuration in which forces on all atoms vanish. Earlier approaches involved the repeated solution of the self-consistent electronic eigenvalue problem for each atomic arrangement encountered in the energy minimization process. Recently, a number of techniques have been developed in which the electronic and nuclear degrees of freedom are treated simultaneously (Car and Parrinello, 1985; Payne et al., 1986, 1992). These techniques have resulted in vast improvements in efficiency for the computation of both equilibrium configurations and dynamics. Other advances enable the treatment of complex materials (Ellis et al., 1990) and consider the use of alternate quantum mechanical descriptions from ab initio quantum computations or ab initio pseudopotential methods.

For perfect periodic crystals, the size of the problem is determined by the number of atoms in the primitive unit cell. When periodicity is broken by the presence of surfaces, interfaces, and defects, the problem becomes mathematically more complex. A number of methods have been developed to handle this reduced symmetry. Some methods are based on scattering theory or Green's functions (Appelbaum et al., 1975; Baraff and Schlüter, 1978; Bernholc et al., 1978). With the advent of very powerful computers, however, these techniques have lost ground to a simpler approach based on supercells. For example, a large cell containing a point defect is repeated periodically. The cells must be large enough so that interactions between defects in neighboring cells are insignificant. Similarly, studies of surface properties or of interfaces use periodically repeated slabs. Solutions are typically constructed by expanding the eigenfunctions in terms of a basis set. A variety of basis sets have proven useful in different cases, but the power of modern computers favors the use of plane waves because the size of the basis set can be increased systematically until convergence is achieved (Ihm et al., 1979). The treatment of dynamical properties, such as conductivity, is generally carried out with semiempirical Hamiltonians, for example, a tight-binding Hamiltonian that is constructed with a localized basis set of atomic orbitals (Khan and Broughton, 1991). The Hamiltonian is defined in terms of its matrix elements, but there are major conceptual and mathematical problems associated with the theoretical determination of these Hamiltonian matrix elements; these offer research opportunities.

When material properties may adequately be computed by first averaging out the electrons, the difficult task still remains of constructing suitable interatomic potentials (Baskes et al., 1992; Smith and Srolovitz, 1992). Analytical expressions are typically assumed that incorporate pairwise interactions plus additional many-body terms, depending on the nature of the material. The parameters in the potentials are tuned to reproduce either experimental data (for example, lattice constants, elastic constants, and so on) or theoretical results calculated by first-principles techniques. Once constructed, the interatomic potentials are used to predict a variety of other equilibrium and dynamical properties of both the

materials whose properties are used to construct the potentials and other materials containing the same constituent atoms. Similar approaches are applied to the description of polymers, proteins, and glasses.

The construction of accurate and reliable interatomic potentials remains a pressing problem that can benefit from new mathematical sciences research. First-principles calculations provide an enormous database concerning the interaction potentials. However, the construction of computationally convenient interatomic potentials involves a many-parameter, highly nontrivial, nonlinear fitting problem.

Currently, the frontiers in atomic-scale theories concern the computation of equilibrium atomic configurations and the evolution of nonequilibrium configurations at nonzero temperatures. The former is usually carried out with direct energy minimization techniques or with Monte Carlo simulation methods (Hayes, 1993). The latter is carried out using molecular dynamics, which can be implemented with a first-principles Hamiltonian, a semiempirical tight-binding Hamiltonian, or with interatomic potentials. The method is in effect an implementation of statistical mechanics. In order to calculate the evolution of a system, externally controlled variables such as pressure or volume, temperature, and so on must be specified. Such specification defines a statistical mechanical ensemble, and the calculation simulates the evolution of a single member of the ensemble in terms of the "trajectories" of all the particle coordinates and momenta in phase space. Equilibrium properties are calculated as time integrals. In addition, elastic constants and linear transport coefficients are obtained using formulas derived from fluctuation theory. Finally, the detailed dynamical evolution allows one to unravel the atomistic mechanisms that underlie a variety of processes and phenomena. For example, studies have been made of dislocation formation and motion, plastic flow, grain boundary sliding and strengthening, crack propagation, tribological phenomena, chemical reactivity, and so forth (Landman et al., 1992; Stillinger and Weber, 1987).

The mathematical sciences have made valuable contributions to the developments described above, especially to the numerical implementations. Examples are fast Fourier transforms, multidimensional integration, curved-space description (Riemannian metric), solutions of nonlinear equations, nonlinear regression, conjugate-gradient methods, and eigenvalue methods for large or sparse matrices.

The main limiting factors in computational atomic-scale research today are the number of atoms that can be treated and the length of the simulation. Parallel computers are becoming an important ingredient in the impetus to study larger systems. The current limits are a few hundred atoms for first-principles calculations and several hundred million atoms for conventional molecular dynamics using interatomic potentials. In time scales, first-principles calculations are limited to picoseconds, whereas conventional molecular dynamics is limited to nanoseconds. Mathematical challenges lie in the development of algorithms to handle large numbers of atoms, in the analysis of instabilities that may arise in the governing equations, and in the development of general asymptotic methods to pass to larger scales.

The above does not exhaust interesting areas involving electronic structure and other atomic-scale theories for properties of materials. One other area is electron transport. In most cases, the classical Boltzmann transport equation is adequate, with the electrons (and

holes in the case of semiconductors) having the energy dispersions dictated by the material's quantum energy bands. Boltzmann's equation can be solved in various degrees of approximation. The drift-diffusion approximation has been the workhorse of the microelectronics industry in modeling devices (Cole et al., 1990). At the next level of sophistication, the hydrodynamic approximation, more accurate solutions are possible (Blotekjaer, 1970; Rudan and Odeh, 1986). Finally, exact solutions of the Boltzmann equation are possible by Monte Carlo techniques (Jacoboni and Reggiani, 1983; Fischetti and Laux, 1988, 1991). Mathematicians played key roles in many of these developments. Computer implementations of the Monte Carlo techniques are, however, extremely time consuming, calling for the development of more efficient algorithms (see Chapter 8 ).

Some currently available mathematical results may provide new inspiration to materials scientists if this well-known mathematics is brought to their attention. A recent example is the study of graphite-like carbon structure with negative Gaussian curvature (Mackay and Terrones, 1991). Construction of such structures depends on the mathematics of periodic minimal surfaces that was developed more than a century ago (Schwarz, 1890; Weierstrass, 1866). These minimal surfaces are also relevant to the study of microphase separation in block copolymer systems (see Chapter 3).

Strong electronic correlations play a dominant role in phenomena such as the metal-nonmetal transition and (probably) high-temperature superconductivity. Standard density functional theory in the local density approximation does not accommodate these strong correlations, and even ab initio quantum mechanical methods encounter difficulties when many electrons must be correlated. Typically, one has to construct a model theory that captures the essence of the strong correlations. Quite often, a model theory, though simple in appearance, may present significant mathematical challenges to obtaining practical solutions to specific problems. Performing direct numerical calculations for the many-body problem of strongly correlated electrons is another frontier. The main limitation is again on the number of electrons that can be accommodated. There are efforts to construct model Hamiltonians with reduced degrees of freedom by introducing a mapping from the true Hamiltonian (McMahan et al., 1990), but these efforts remain heuristic in nature. Some progress, however, is available for small-molecule systems (Freed, 1983, 1989). Further mathematically rigorous developments in this area would be very valuable.

Many workers have explored methods beyond density functional methods in order to describe quantum effects for systems in which material properties involve excited states, be they optical or high thermal excitations. Wiener path integral methods have been extensively used to model (Chandler and Wolynes, 1981; Berne and Thirumalai, 1986) the thermal properties of the solvated electron (Rossky and Schnitker, 1988), the role of quantized rotational and librational motion in liquid water (Kuharski and Rossky, 1985), superfluid transitions in liquid helium (Ceperley and Pollock, 1984), melting in rare-gas clusters (Doll et al., 1990), and the thermal behavior of both nuclei and electrons in sodium clusters (Hays and Hall, 1991; Hall and Prince, 1991; Hall, 1992). Simulated annealing techniques have enabled the computation of electronically excited potential surfaces in conjunction with molecular dynamics simulations of the nuclear motion. The first thermal simulations combining the full electronic Schrödinger equation and molecular dynamics involve simulations of both ground and electronically excited states of metal atoms solvated in rare-

gas clusters, and yield both structural and optical properties of clusters exhibiting isomerization and phase transitions (Tsoo et al., 1990, 1992; Estrin et al., 1992). Simulated annealing has been applied to the Na₄ cluster using a generalized valence bond electronic wavefunction, but only a 500-femtosecond test run has been reported so far (Hartke and Carter, 1992). Ab initio quantum electronic structure methods for simulations would benefit if mathematical scientists improved the computer algorithms or advanced current theory.

Atomic-scale theories have proved to be very valuable in elucidating the properties of all sorts of materials. Electronic structure calculations have been essential in understanding electronic, optical, and transport properties, including elastic constants, conductivities, magnetic properties, and many more. Molecular dynamics is an explanatory mechanism and contributes to determining structural and mechanical properties of complex materials. It should be borne in mind, however, that atomic-scale theories are limited to length scales of a maximum of a few hundred angstroms. Most properties of complex materials, however, are determined by their collective microstructure (grain boundaries, precipitates, inclusions, microvoids, and so on) where the length scale is measured in microns and beyond (Baskes et al., 1992; Pantelides, 1992). Similarly, time scales of atomic simulations are measured in picoseconds or nanoseconds, whereas many diffusive and slip phenomena that underlie plasticity and other processes occur over time scales of microseconds, milliseconds, and beyond. It is therefore necessary to develop theories relating results obtained on atomic scales to desired material properties on more macroscopic length and time scales. A wide-open frontier for new mathematical sciences research is that of building connecting links between materials science theories developed for different length scales.