Reduced Bloch mode expansion for periodic media band structure calculations

Mahmoud I. Hussein


Reduced Bloch mode expansion (RBME) is presented for fast periodic media band structure calculations. The expansion employs a natural basis composed of a selected reduced set of Bloch eigenfunctions. The reduced basis is selected within the irreducible Brillouin zone at high symmetry points determined by the medium’s crystal structure and group theory (and possibly at additional related points). At each of the reciprocal lattice selection points, a number of Bloch eigenfunctions are selected up to the frequency/energy range of interest for the band structure calculations. As it is common to initially discretize the periodic unit cell and solution field using some choice of basis, RBME is practically a secondary expansion that uses a selected set of Bloch eigenvectors. Such expansion therefore keeps, and builds on, any favourable attributes a primary expansion approach might exhibit. Being in line with the well-known concept of modal analysis, the proposed approach maintains accuracy while reducing the computation time by up to two orders of magnitudes or more depending on the size and extent of the calculations. Results are presented for phononic, photonic and electronic band structures.

1. Introduction

In the early years following the birth of quantum mechanics through the work of de Broglie, Schrödinger, Heisenburg and others (see Messiah (1964) for a brief history), the nature of non-interacting electron motion in a perfect lattice was revealed by Felix Bloch through the formulation of Bloch theory (Bloch 1928). His theory has rigorously, yet simply, formulated the concept of electron bands in crystals, hence the synonymous term band theory. Subsequent work based on band theory and the Pauli exclusion principle (Pauli 1925) has led to the formal classification of all crystals into metals, semiconductors and insulators (Peierls 1929, 1930; Wilson 1931a,b). Band theory has since played a major role in elucidating the electrical properties of crystalline solids. Band theory has also provided a basis for the study of optical, thermal and magnetic properties of crystals in association with the motion of other quantum particles: photons, phonons and magnons. In the late 1980s, yet another type of periodic system emerged only to ignite further research interest in Bloch wave propagation, namely, photonic and phononic crystals (John 1987; Yablonovitch 1987; Ho et al. 1990). These are macroscale periodic materials that can be designed, via band engineering, to classically control the flow of light, sound and/or mechanical vibrations in a predetermined manner within the solid. In doing so, these modern materials have opened up a new technological frontier in optical, acoustic and elastic devices (e.g. Kushwaha 1996; Joannopoulos et al. 1997; Sigalas et al. 2005; Prather et al. 2006 and references therein).

Central to the application of band theory in all of the above disciplines is the process of band structure calculations. Owing to crystallographic symmetry, the Bloch wave solution needs to be applied only to a single unit cell in the reciprocal lattice space covering the first Brillouin zone (BZ) (Brillouin 1953). Further utilization of symmetry reduces the solution domain, even more, to the irreducible BZ (IBZ). The target of the calculations typically is to generate a plot of the frequency or energy versus wave vector (along various directions inside or along the boundaries of the IBZ). The density of states, as a function of frequency or energy, is also commonly computed.

(a) Classical calculations

Depending on the problem of interest, periodic systems can be modelled in a discrete or a continuous fashion. An important discrete problem in materials science is the calculation of band structure for lattice dynamics systems based on classical interatomic potentials (Maradudin 1974). Discrete problems and relevant modelling methods are also encountered in the field of structural dynamics (Mead 1973; Langley 1996). There are several techniques for band structure calculations for classical continuous systems such as photonic and phononic crystals. These include the plane-wave method (Ho et al. 1990; Leung & Liu 1990; Zhang & Satpathy 1990; Meade et al. 1993; Johnson & Joannopoulos 2001), the transfer matrix method (Pendry & MacKinnon 1992), the multiple scattering method (Leung & Qiu 1993; Wang et al. 1993), the finite-difference time-domain method (Chan et al. 1995), the finite-difference method (Yang 1996), the finite-element (FE) method (Axmann & Kuchment 1999; Dobson 1999; Pask et al. 2001), the meshless method (Jun et al. 2003), the multiple multipole method (Moreno et al. 2002), the wavelet method (Checoury & Lourtioz 2006; Yan & Wang 2006) and the pseudospectral method (Chiang et al. 2007), among others (Botten et al. 2001; Marrone et al. 2002; Yuan & Lu 2006) (see Busch et al. (2007) for a review). Some of the methods involve expanding the periodic domain (or potential) and the wave field using a truncated basis. This provides a means of classification in terms of the type of basis, e.g. the plane-wave method involves a Fourier basis expansion and the FE method involves a real space basis expansion. The pros and cons of the various methods are discussed in depth in the literature, e.g. Busch et al. (2007).

Regardless of the type of system and type of method used for band structure calculations, the computational effort is usually high because it involves solving a complex eigenvalue problem and doing so numerous times as the value of the wave vector, k, is varied. The size of the problem, and hence the computational load, is particularly high for the following cases. The first is when the unit cell configuration requires a large number of degrees of freedom to be adequately described. This, for example, could be due to a complex unit cell material phase topology or due to a complex interatomic potential, both necessitating a finely resolved description. The second case is when the presence of defects is incorporated in the calculations. Defects are known to have a physical influence extending over very long ranges in space. This, in turn, requires choosing correspondingly large unit cells, known as supercells, for the band structure calculations. Consequently, large cells imply large numbers of degrees of freedom. The third case is concerned with an emerging area of research, namely, band structure optimization (Cox & Dobson 1999 for classical problems; Franceschetti & Zunger 1999 for semi-classical problems; and Burger et al. 2004 for a review). By varying certain parameters, the band structure can be tailored to the desired form. Needless to say, optimization studies inherently involve numerous repeated forward calculations.

Some techniques have been developed to expedite classical band structure calculations, examples include utilization of the multigrid concept (Chern et al. 2003), development of fast iterative solvers for the Bloch eigenvalue problem (Axmann & Kuchment 1999; Dobson 1999; Johnson & Joannopoulos 2001) and extension of homogenization methods to capture dispersion (McDevitt et al. 2001; Wang & Sun 2002; Hussein 2004; Nagai et al. 2004; Hussein & Hulbert 2006; Gonella & Ruzzene 2008). In this paper, we provide a fundamentally different approach for fast band structure calculations. The approach, which we refer to as reduced Bloch mode expansion(RBME), is based on an expansion employing a natural basis composed of a selected reduced set of Bloch eigenfunctions. This reduced basis is selected within the IBZ at high symmetry points determined by the crystal structure and group theory (and possibly at additional related points). At each of these high symmetry points, a number of Bloch eigenfunctions are selected up to the frequency range of interest for the band structure calculations. As mentioned above, it is common to initially discretize the problem at hand using some choice of basis. In this manner, RBME constitutes a secondary expansion using a set of Bloch eigenvectors, and hence keeps and builds on any favourable attributes a primary expansion approach might exhibit. The underlying idea of Bloch mode reduction utilizing eigenfunctions/vectors selected at high symmetry points in the IBZ was presented earlier in Hussein (2004) and Hussein & Hulbert (2006). In that work, however, the reduction process was embedded within a multi-scale two-field variational formulation that masked the simplicity and broad applicability of the approach.

The proposed RBME approach is in line with the well-known concept of modal analysis, which is widely used in various fields in the physical sciences and engineering. The concept of modal analysis is rooted in the idea of extracting a reduced set of representative information on the dynamical nature of a complex system. This practice is believed to have originated by the Egyptians in around 4700 BC in their quest to find effective ways to track the flooding of the Nile and predict celestial events. Over time, spectral methods emerged as a tool for the study of cosmology and optics until modal expansion emerged, as presented by Clebsch (1862) and Rayleigh (1877) in his study of elastic vibrations, as a modelling paradigm in its own right. For a thorough discussion on the historical development of modal analysis, the reader is referred to Døssing (1995).

(b) Ab initio calculations

Our RBME approach applies to ab initio band structure calculations as well. There are numerous techniques for this category of periodic media and their classification is beyond the scope of this paper (see Martin (2004) and Nemoshkalenko & Antonov (1998) for reviews and descriptions of available methods). Most relevant to the approach we propose is a range of methods involving mathematical operations in the reciprocal lattice space. These methods focus on predicting physical quantities, including band structure eigenvalues, at any value of wave vector k (i.e. the k-point) in the BZ through a reduced cost calculation that depends on the solution at a specific k-point or a reduced set of k-points. Common among these is the k·p method that describes band dispersion using perturbation theory around regions that are identified to be of physical importance (Luttinger & Kohn 1955; Kane 1956). The perturbative nature of this method constrains its application to local regions in the BZ unless a large set of expansion functions are used (Wood & Zunger 1996). Another method is the Fourier interpolation scheme (Giannozzi et al. 1991; Gonze & Lee 1997), which is based on the well-known special k-point integration schemes (Baldereschi 1973; Chadi & Cohen 1973; Monkhorst & Pack 1976). This method is associated with the application of density functional theory in defining the electronic potential. It involves evaluating the dynamical matrix at a coarse grid of carefully selected special k-points, and then using Fourier interpolation to yield the force constant matrix in real space. The force constants are then inversely transformed in order to economically obtain the dynamical matrix and band solution at any point in the BZ. While well suited for a wide range of solids, there are cases where a considerably refined integration mesh is still required, such as when the real-space force constants have an extended range (Baroni et al. 2001). Another method uses maximally localized Wannier functions (Souza et al. 2001) and the Slater–Koster (1954) interpolation scheme for efficient computation of dispersion curves. This is a non-variational method and hence does not produce interpolated solutions that satisfy the governing equations (Pau 2007). Other methods have been proposed that utilize physically significant basis functions, such as pure spin basis functions in relativistic calculations (MacDonald et al. 1980). A recent study (Pau 2007) presents an approach that involves a projection operation onto a reduced basis consisting of solutions obtained at certain k-points. These k-points are determined by a posteriori error estimator. The RBME approach provides an alternative to the ab initio model reduction methods mentioned above and has the following qualities: (i) it covers the entire BZ, (ii) it is suitable for media with extended interatomic force ranges, (iii) it sustains the variational character of a problem when it exists, (iv) it does not require a built-in posteriori error estimator, and (v) it is simple and formulated in a manner that allows it to apply ‘on top’ of a chosen primary expansion technique.

In the next sections, a description of direct, full-model-based calculations of band structures is given, followed by a description of the RBME process and its application in a discretized setting using the FE method. Numerical examples are given for band structure, Bloch mode shape and density of state calculations. Some results are then presented to demonstrate the accuracy and efficiency of the calculations following the proposed approach. Linear elastic phononic media are chosen as a platform for presenting the formulations and results in these sections. Other sections cover the implementation of RBME for photonic and electronic band structure calculations. The application of the proposed approach to the plane-wave method is also presented, in the context of photonic crystals. Finally, a summary and conclusions section is provided.

2. Direct implementation of Bloch theory

Bloch theory (Bloch 1928) describes the wave function of a particle in an infinite periodic medium, such as a crystal, in terms of wave functions at the reciprocal space vectors of a Bravais lattice (the one-dimensional temporal version of Bloch theory is known as Floquet theory; Floquet (1883)). The Bloch response, f, is a product of a Bloch periodic function, Embedded Image, and a phase multiplier Embedded Image 2.1 where x={x,y,z) is the position vector, k={kx,ky,kz) is the wave vector, Embedded Image and ω and t denote frequency and time, respectively. Bloch theory can also be expressed as Embedded Image 2.2 where R denotes primitive lattice vectors.

We now develop the application of Bloch theory in the context of three-dimensional phononic crystals modelled as a continuous linear elastodynamic periodic medium. Actual numerical implementation is covered in §5a for a two-dimensional plane-strain model of a phononic crystal. The application of Bloch theory to other types of periodic media, such as for photonic and electronic band structure calculations, are briefly presented in §5b,c.

The governing continuum equation of motion for a heterogeneous medium is Embedded Image 2.3 where σ is the stress tensor, u is the displacement field, ρ is the density and a superposed dot denotes differentiation with respect to time. For an elastic medium, Embedded Image 2.4 where C is the elasticity tensor and ∇S denotes the symmetric gradient operator, that is, Embedded Image 2.5 In equation (2.5), (·)T denotes a transpose operation. Substituting equation (2.4) into equation (2.3) yields Embedded Image 2.6 which is the strong form of the general elastodynamic problem. Using operator notation, equation (2.6) can be expressed as Embedded Image 2.7 where Embedded Image and L[u]=∇·C:∇Su. A three-dimensional periodic medium is characterized by an infinitely repeated unit cell ω. In a phononic crystal, the unit cell is composed of two, or more, material phases. We will assume the material-to-material interfaces to be ideal. For the periodic unit cell domain ω, equation (2.6) has a Bloch solution of the form Embedded Image 2.8 where Embedded Image is the displacement Bloch function that is periodic in the domain ω. Using equation (2.8), the displacement gradient is Embedded Image 2.9 where the symbol ⊗ denotes outer product. Substitution of equations (2.5), (2.8) and (2.9) into equation (2.6) gives the strong form of the Bloch eigenvalue problem Embedded Image 2.10 that is, Embedded Image 2.11 where λ=−ω2, Embedded Image and Embedded Image is the Bloch operator.

Using variational principles, the weak form for equation (2.6) can be expressed as Embedded Image 2.12 where w is a weighting function. Upon Bloch transformation using equation (2.8), the weak form for the Bloch eigenvalue problem is (Hussein 2004; Hussein & Hulbert 2006) Embedded Image 2.13 where Embedded Image is the Bloch-transformed weighting function. The operator on the l.h.s. of equation (2.13) arising from the undergone transformation is Hermitian. In general, the r.h.s. operator can also be Hermitian, as this depends on the form of Embedded Image.

With periodic boundary conditions applied on ω, equation (2.13) is solved for wave vector k spanning the first BZ over the unit cell. For a simple cubic crystal, for example, k∈[−π/d,π/d]3, where d is the side length of the cell. When symmetry allows (which is the usual case for crystals), it is sufficient for k to cover only the IBZ whose vertices are the high symmetry points determined by the crystal’s structure and group theory (points shown in figure 1 for a two-dimensional square lattice and a three-dimensional simple cubic lattice, as examples). Following common practice, the k-space for band diagram calculations is narrowed even further to only considering wave vectors pointing to positions along the border of this zone (circuit lines shown in figure 1). Considering a two-dimensional square lattice, the considered wave vector endpoints are along the paths Γ→X, X→M and M→Γ. For a three-dimensional simple cubic lattice, the IBZ circuit consists of the paths Γ→X, X→M, M→R and R→Γ. For a unit cell with side length d, k=[0,0,0] at Γ,k=[π/d,0,0] at X, k=[π/d,π/d,0] at M and k=[π/d,π/d,π/d] at R. This process generates a representative band structure (frequency versus wave vector relation) for the periodic medium. For computational purposes frequencies ωp(k), where p denotes branch number, are only evaluated at a discrete (uniformly distributed) set of nk points along the entire IBZ border. The discrete wave vector set is denoted kj=[k1,k2,k3]j, j=1,…,nk. The number of k-points on a single straight line along the IBZ border connecting two high symmetry points is denoted as lk, i.e. nk=3(lk−1)+1 and nk=4(lk−1)+1 for the two-dimensional and three-dimensional cases, respectively, shown in figure 1. For density of states calculations, the k-space sampling covers the entire IBZ, not only the border.

Figure 1.

Unit cell in reciprocal lattice space with the IBZ, high symmetry k-points (solid circles) and intermediate k-points (hollow circles) shown. (a) Two-dimensional square unit cell and (b) three-dimensional simple cubic unit cell.

3. Reduced Bloch mode expansion

In this section, the formulation for RBME is given, as well as two eigenfunction selection schemes for the reduced basis.

(a) Formulation

Recall the Bloch eigenvalue problem given in equations (2.10) or (2.11). The displacement Bloch function will now be expressed as a superposition of Bloch mode eigenfunctions Embedded Image, Embedded Image, Embedded Image 3.1 where βj is the coordinate for mode j. Equation (3.1) represents a linear transformation to modal coordinates. Substituting into equation (2.11) gives Embedded Image 3.2 Equation (3.2) is then multiplied by Embedded Image and integrated over the unit cell domain ω, Embedded Image 3.3 To achieve the targeted reduction in the size of the problem, only a small truncated set of m modes is retained for the Bloch mode expansion, Embedded Image 3.4 If we utilize all the Bloch eigenfunctions corresponding to the k-point eigenvalue problem for which the Bloch mode expansion is implemented, then this is essentially a unitary similarity transformation. Using the same eigenfunctions, or a broader selection of eigenfunctions as described in §3b, for transforming the eigenvalue problem at other k-points is a deviation from this condition.

(b) Choice of reduced expansion basis

The proposed RBME approach will only produce accurate results if a suitable choice of reduced basis, Embedded Image, j=1,…,m, is made. For this selection, we propose the two schemes that were presented in Hussein (2004) and Hussein & Hulbert (2006). The first consists of choosing eigenfunctions corresponding to the high symmetry points that characterize the periodic lattice, as determined by the crystal structure and group theory. Figure 1 shows these points as solid circles for a two-dimensional square lattice and a three-dimensional simple cubic lattice, i.e. the Γ, X and M points and the Γ, X, M and R points, respectively. The second scheme involves choosing eigenfunctions at all the high symmetry points mentioned above, as well as at intermediate points centrally intersecting the straight lines joining these high symmetry points along the border of the IBZ. In practice, these would be points along the particular IBZ circuit lines that are intended for the band structure calculations, i.e. all the circles (solid and hollow) representing the Γ, Δ, X, Z, M and Σ points in two dimensions and the Γ, Δ, X, Z, M, T, R and Λ points in three dimensions, as shown in figure 1. This scheme is therefore richer and slightly more computationally expensive than the previous one. We refer to the first selection scheme as a 2-point expansion scheme, because eigenfunctions corresponding to two k-points on each of the straight lines bordering the IBZ are used. Similarly, we refer to the second selection scheme as a 3-point expansion scheme. In principle, other schemes may be developed and can be tailored to the particular application and purpose for the band structure calculations. Furthermore, resourceful utilization of the Bloch eigenfunctions can be made by only using those eigenfunctions that correspond to selection k-points belonging to the current IBZ circuit line that is being evaluated. An equally important consideration in setting up the RBME basis is the number and choice of eigenfunctions to retain at each of the k-points selected. In order for the band structure calculations to be accurate up to the frequency range of interest, this number should be equal to, or slightly higher than, the number of dispersion branches that are to be computed.

In this manner, both spatial and temporal dynamical characteristics are considered in the selection of the RBME basis. The spatial characteristics are related to the selection of the k-points and the temporal characteristics are related to the choice of the eigenfunctions (in frequency space) to retain at each of the selected k-points. The approach is, in principle, applicable to any type of lattice symmetry, not just square and cubic lattices. The utilization of periodic symmetry has been used previously in photonic crystals analysis, namely in the context of cutting the number of transfer operations in the transfer matrix method (Li & Ho 2003; Matias et al. 2003) and cutting the number of plane waves employed in the plane-wave method (Zhang et al. 2004).

4. Finite-element discretization

The FE method is one of the various methods for discretizing a periodic medium to facilitate numerical calculation of the band structure. A brief description of FE implementation for phononic band structure calculations is now provided.

The Galerkin FE method (explained in detail in Hughes (1987)) involves introducing a finite-dimensional approximation of the solution domain, its boundaries and the governing variational eigenvalue problem (2.13) that is to be solved. Upon discretization, the unit cell domain ω is divided into nel finite elements, Embedded Image 4.1 where ωe denotes the domain for each element, e=1,…,nel. The solution field, Embedded Image, is discretized locally in domain ωe using well-chosen local piecewise polynomials, known as shape functions Na, Embedded Image 4.2 where da is the discrete displacement vector at node a. Upon application to equation (2.13) with incorporation of periodic boundary conditions, the following algebraic eigenvalue problem emerges: Embedded Image 4.3 where M and K are the global mass and stiffness matrices, respectively, and Embedded Image is the discrete Bloch vector that is periodic in the domain ω. If the chosen weighting function in equation (2.13) is not a function of k, then M is independent of k (as is assumed in equation (4.3)). The global matrices M and K are assembled from the element-level matrices, as shown in appendix A. Also provided in appendix A is the definition of what we refer to as a Bloch element, that is, the form of the FE matrices associated with the Bloch transformation.

Solution of equation (4.3) at selected k-points provides the eigenvectors from which a reduced Bloch modal matrix, denoted Ψ, is formed column-wise. This matrix is used to expand Embedded Image in analogy to the continuous version of RBME expressed in equation (3.4). In discrete form, the expansion follows the matrix multiplication Embedded Image 4.4 where Embedded Image is a vector of modal coordinates for the selected unit cell Bloch mode shapes. In equation (4.4), n and m refer to the number of rows and number of columns for the matrix equation and, respectively, denote the total number of degrees of freedom and the total number of retained Bloch modes. For the k-point selection schemes described in §3b, m<<n. Substituting equation (4.4) into equation (4.3), and pre-multiplying by Ψ* (where (·)* denotes a complex conjugate transpose operation), Embedded Image 4.5 yields a reduced eigenvalue problem of size m×m, Embedded Image 4.6 where Embedded Image and Embedded Image are the generalized mass and stiffness matrices. From a coding point of view, Ψ is constructed ahead of the main band structure calculations and so is the Embedded Image matrix. The Embedded Image matrix is set up before the calculations and is evaluated inside the loop spanning the wave vector space. Such features are favourable for both speed and memory utilization.

Remark 4.1

The fact that the operator associated with Bloch transformation is Hermitian is a favourable property for discretizing and solving the full-model eigenvalue problem that emerges. The process of strict diagonalization via mode expansion requires the selected eigenvectors to be orthogonal with respect to the system matrices, e.g. M and K in equation (4.3) for the elastodynamic problem. Owing to the use of Bloch eigenvectors across different k-points in the proposed approach, this orthogonality condition is not held exactly in the RBME methodology. Nevertheless, if the system matrices are Hermitian, the reduction will be successful. While the Bloch operator is Hermitian, it is essential that this property is sustained through the discretization, such as the case with the FE method (and the plane-wave method as demonstrated in §7).

5. Numerical examples

(a) Phononic band structure calculations

We consider a linear elastic, isotropic, continuum model for a two-dimensional phononic crystal under plane-strain conditions. Applying Newton’s equations, coupled in-plane P (longitudinal) and SV (shear vertical) wave propagation modes can be predicted. The algebraic Bloch eigenvalue problem is given in equation (4.3); details on the mass and stiffness matrices are provided in appendix A. The reduced elasticity matrix, denoted D, for this model has the form Embedded Image 5.1 where λ and μ denote Lamé’s constants.

As an example, a square lattice is considered with a bi-material unit cell. One material phase is chosen to be stiff and dense (denoted by subscript ‘2’) and the other compliant and light (denoted by subscript ‘1’). In particular, a ratio of Young’s moduli of E2/E1=16 and a ratio of densities of ρ21=8 are chosen for all calculations presented. A Poisson ratio of 0.34 is assumed for both phases. The topology of the material phase distribution in the unit cell is shown in figure 2. The unit cell is discretized into 45×45 uniformly sized 4-node bilinear quadrilateral finite elements, i.e. nel=2025 (mesh shown in figure 2). With the application of periodic boundary conditions, the number of degrees of freedom is n=4050. The k-space is discretized such that lk=49; thus, a total of nk=145 k-points are evaluated to generate the band structure. Figure 3 shows the calculated band structure using 2-point and 3-point expansions. For generalization, k in the figure, and in the rest of the paper, is used to describe the non-dimensional wave vector. Density of state calculations are also provided in figure 3 based on a total of nk=561 k-space sampling across the IBZ. In both sets of calculations, eight modes (p=1,…,mp; mp=8) were selected at each of the selection points in k-space, i.e. a total of 24 eigenvectors (m=24) were used to form the Bloch modal matrix for the 2-point expansion calculations and a total of 48 eigenvectors (m=48) for the 3-point expansion calculations. The results for the full model are overlaid for comparison, indicating excellent agreement (despite the reduction in model size from 4050×4050 to 24×24 and 48×48, respectively). The slightly more expensive 3-point expansion scheme gives better results at k-space regions far away from the high symmetry points Γ, X and M. The convergence of the results as a function of mp is shown in figure 4 for two specific points in the band diagram: the third branch at the Γ point (point A in figure 3a,b) and the first branch at the Z point (point B in figure 3a,b). The error plotted here is that of the RBME results when compared with the ‘benchmark’ full model. It is observed that, for point A, a significant drop in error occurs when three or more modes per selection point are retained, highlighting the importance of retaining modes up to at least the frequency of interest. The results show that the error for point A with mp=8 is 0.2741 per cent based on a 2-point expansion and 0.2981 per cent based on a 3-point expansion. For point B, with mp=8, the error is 1.7683 per cent based on a 2-point expansion and 4.4913×10−11 per cent based on a 3-point expansion. The latter value is very small because the 3-point expansion employs modes associated with the Z-point that coincides with point B. In general, modes (eigenvectors) associated with more intermediate selection points in the k-space, and more frequencies per k-space selection point, can be added to reduce error as needed. The RBME approach also performs well in predicting Bloch mode shapes, as demonstrated in figure 5. The displacement field (based on 3-point expansion) in figure 5a,b represents, respectively, the Bloch mode shapes corresponding to the same points A and B. Calculations using the same parameters were also carried out for internal segments in the IBZ, as illustrated in the insets of figure 6. The band structures (based on 2-point expansion) for the Γ→XZ, Γ→Z and Γ→ZM paths are shown in the figure. Excellent agreement is again observed in comparison with the full-model results. Finally, calculations for unit cells with elasticity contrast ratio of up to E2/E1=2000 (while keeping the density ratio the same) were carried out (not shown), and the reduced model results accurately matched that of the full model.

Figure 2.

Finite-element mesh for a two-dimensional square unit cell. The stiff/dense material phase is in black, and the compliant/light material phase is in white.

Figure 3.

Phononic band structure and density of states (DOS) calculated using the full model (matrix size: 4050×4050) and the RBME model following (a) 2-point expansion (matrix size: 24×24, black solid line, full model; grey dashed line, RBME model) and (b) 3-point expansion (matrix size: 48×48, black solid line, full model; grey dashed line, RBME model). The two-dimensional unit cell, IBZ and eigenvector selection points are shown in the insets. The FE method was used for the primary expansion.

Figure 4.

Percentage error, e, in calculated frequency for points A and B in the band diagram (see figure 3) versus number of retained Bloch modes per selection point, mp. Results for both 2-point (black lines) and 3-point (grey lines) expansions are shown.

Figure 5.

Mode shapes for points (a) A and (b) B in the frequency spectrum (see figure 3). Vector fields from the full model (black arrows) and an RBME model following 3-point expansion (grey arrows) are shown.

Figure 6.

Phononic band structure for different k-point paths in the interior of the IBZ. Finite-element calculations from the full model (black solid line) and an RBME model following 2-point expansion (grey dashed line) are shown. The two-dimensional unit cell, IBZ, k-point path and eigenvector selection points are shown in the insets.

(b) Photonic band structure calculations

The RBME approach is now applied to photonic band structure calculations. A two-dimensional model is considered for lossless electromagnetic waves propagating in-plane (i.e. the xy plane). In this case, TE-(H field in the z direction) and TM-(E field in the z direction) polarized waves are, respectively, described by the following form of Maxwell’s equations: Embedded Image 5.2 and Embedded Image 5.3 where ϵr is the dielectric constant and c is the speed of light. The photonic crystal medium is periodic, i.e. ϵr(x+Rxy)=ϵr(x), where Rxy are primitive lattice vectors with zero z-components. The geometry and material phase topology of the unit cell shown in figure 2 are again considered for an example problem. The two materials chosen are GaAs (represented in black) and air (represented in white). The ratio of the dielectric constants is 11.4. Figure 7a,b shows the photonic band structures for TE polarization and TM polarization, respectively, using the same number of finite elements and the same k-space sampling parameters as in §5a. Unlike the phononic crystal case above, in which the wave field is a vector field, the present models are governed by scalar equations (5.2) and (5.3). The size of the problem in each is therefore smaller, n=2025. The RBME results, which, in this example, are based on 2-point expansion and the use of eight Bloch modes for every selected k-point (i.e. mp=8 and m=24), are in excellent agreement with those of the full model. The error with respect to the full model for point A (location marked in the band diagrams) is 1.4361×10−12 per cent for TE polarization and 0.095 per cent for TM polarization. For point B (also shown in the band diagrams), the error is 0.0013 per cent and 5.0592×10−4 per cent for TE and TM polarization, respectively.

Figure 7.

Photonic band structure calculated using the full model (matrix size: 2025×2025, black solid line) and RBME model following 2-point expansion (matrix size: 24×24, grey dashed line). (a) TE polarization and (b) TM polarization. The two-dimensional unit cell, IBZ and eigenvector selection points are shown in the insets. The FE method was used for the primary expansion.

(c) Electronic band structure calculations

In this section, we demonstrate the applicability of RBME to electronic band structure calculations for a three-dimensional model. We start with the time-independent single electron Schrödinger equation Embedded Image 5.4 where Embedded Image is the Planck’s constant and m, E(k) and ψ(k,x) are, respectively, the electron effective mass, energy and wave function. The electronic potential V (x) in a perfect crystal satisfies the relation V (x)=V (x+R). We consider a three-dimensional generalized Kronig–Penney model potential (Kronig & Penney 1931) for a simple cubic lattice, Embedded Image 5.5

We choose, for an example problem, Rydberg (Ry) atomic units (a.u.) and the following values for the potential parameters: a=2 a.u., b=3 a.u. and V0=6.5 Ry (Pask et al. 2001). The unit cell is discretized into 18×18×18 uniformly sized 8-node trilinear hexahedral finite elements, i.e. nel=5832. With the application of periodic boundary conditions, the number of degrees of freedom is n=5832 (noting that equation (5.4) is a scalar equation). The k-space is discretized such that lk=49; thus, a total of nk=193 k-points are evaluated to generate the band structure. Figure 8 shows the electronic band structure following the k-space paths Γ→X→M→R→Γ. As in the previous cases, the RBME results (based here on 2-point expansion and use of eight Bloch modes for every selected k-point, i.e. mp=8 and m=32) provide an excellent approximation to the results of the full model, noting that the reduction in the size of the matrix problem here is from 5832×5832 to 32×32. The error with respect to the full model for points A and B (locations shown in figure 8) is 4.9525×10−11 per cent and 0.0019 per cent, respectively.

Figure 8.

Electronic band structure calculated using full model (matrix size: 5832×5832, black solid line) and RBME model following 2-point expansion (matrix size: 32×32, grey dashed line). The three-dimensional unit cell, IBZ and eigenvector selection points are shown in the inset. The FE method was used for the primary expansion.

6. Computational efficiency

In addition to the quality of the primary expansion method, a key component in efficient band structure calculations is the technique and algorithm utilized for the solution of the full Bloch eigenvalue problem. Several techniques and algorithms have been developed in the literature for fast eigenvalue problem calculations, e.g. Axmann & Kuchment (1999), Dobson (1999) and Johnson & Joannopoulos (2001), with the computational complexity for formulations yielding sparse Hermitian matrices roughly being around Embedded Image for a single k-point calculation. For models with many degrees of freedom (large value of n), and for calculations with a large number of k-points, nk, to be evaluated, the total cost of generating a band structure or density of states plot can be prohibitive. The RBME approach has the advantage that it is compatible with available primary expansion approaches (as demonstrated for real-space FE expansion in the above examples and for Fourier space plane-wave expansion in §7). It also utilizes efficient techniques and algorithms for solving the Bloch eigenvalue problem of the full model at the selection points (for setting up the modal matrix). Therefore, any positive attributes pertaining to employed methods from the point of view of efficiency are, in principle, still kept and utilized in the RBME model.

A series of calculations were performed on a single computer processor involving two-dimensional plane-strain phononic band structure calculations in order to assess the computational efficiency of the RBME approach as a function of the size of the numerical problem. First, a series of calculations with increasing number of finite elements, nel, are considered with the value of the k-space sampling rate fixed to lk=49 (nk=145). As an inverse measure of efficiency, the ratio r is defined as Embedded Image 6.1 that is, the lower the value of r, the higher the efficiency. The time taken to obtain the modal matrix, Ψ, is included in the numerator. Figure 9a shows the value of r as a function of nel for both the 2-point and 3-point expansion schemes. It is observed that, beyond 500 finite elements, the value of r almost levels off at values ranging from 0.05 to 0.15, with the 3-point expansion being slightly more expensive than the 2-point expansion as expected. In figure 9b, the results for two specific finite-element resolutions were considered, a coarse resolution with nel=729 and a fine resolution with nel=3969. A series of calculations were carried out to generate r as a function of lk, again for both the 2-point and 3-point expansion schemes. A significant decrease in the value of r is observed as lk is increased, and, in particular, it is shown that an increase in the number of degrees of freedom (as represented by the finite-element resolution) adds substantially to the savings brought about by RBME. For the nel=3969 model with 2-point expansion, two orders of magnitude reduction in computational cost is recorded in these calculations, as shown in figure 9b. Further reduction in computational cost is expected with more FE model refinement and higher k-space sampling rates.

Figure 9.

Computational efficiency: ratio of RBME-model to full-model calculation times, r, versus (a) number of finite elements, nel (for lk=17, nk=49) and (b) number of sampled k-points along the border of the IBZ, nk (for two two-dimensional finite-element meshes). Black line, 2-point RBME; grey line, 3-point RBME.

7. Application of reduced Bloch mode expansion to the plane-wave method

In this section, the RBME approach is applied to the plane-wave method. As an example, we will again consider the case study of TE-polarized wave propagation in a photonic crystal with squared elements (the same problem covered in §5b using the FE method). Following Sakoda (2004), the reciprocal of the dielectric constant, 1/ϵr(x), in equation (5.2) is expanded in a Fourier series, Embedded Image 7.1 where Gxy denotes reciprocal lattice vectors with zero z-component, and so is the field variable Hz(x), Embedded Image 7.2 Substituting equations (7.1) and (7.2) into equation (5.2) results in a system matrix A(k) that defines the eigenvalue problem which is to be solved to obtain the band structure. Full-model implementation details for unit cells with rectangular dielectric elements are provided in Plattner (2003). One important consideration for using plane waves as the primary bases is the orthogonality of Bloch eigenvectors with respect to each other between sets obtained from different selection k-points. While this apparently was not a necessity for the FE case, it is here. We therefore orthonormalize (e.g. using a Gram–Schmidt scheme) the sets of retained Bloch eigenvectors in the modal matrix Ψ ahead of the full-model reduction. In an analogous manner as in §6, once the matrix A obtained from the plane-wave expansion is in place, pre- and post-multiplication by Ψ* and Ψ are, respectively, carried out to generate the reduced matrix Embedded Image, Embedded Image 7.3 Results for the same unit cell geometry, material topology, dielectric constant and k-space sampling rate considered in §5b are shown in figure 10. The band structure is plotted for a full model computed using npw=2601 plane waves and for a 2-point RBME model constructed using m=24 Bloch modes (i.e. a reduction in matrix problem size from 2601×2601 to 24×24). The dispersion curves are again in excellent agreement and are, in fact, practically identical to those obtained by the FE method and shown in figure 7a. Figure 10b shows the percentage error as a function of mp for the two points A and B, marked in figure 10a, with values for mp=8 as low as 3.0157×10−11 per cent for point A and 0.0013 per cent for point B. Finally, computational efficiency is demonstrated in figure 10c, which shows the value of r as a function of npw. When over 6000 plane waves are used, the ratio of calculation times reaches approximately 0.1, indicating an order of magnitude in savings. For a three-dimensional photonic crystal model, with a larger number of plane waves and a finer k-space sampling rate, the value of r is expected to further decrease considerably.

Figure 10.

Results obtained by the plane-wave method for the TE-polarization photonic crystal problem considered in §5b. (a) Band structure based on full model (matrix size: 2601×2601, black solid line) and based on 2-point RBME model (matrix size: 24×24, grey dashed line); unit cell, IBZ, and selection points are shown in the inset. (b) Percentage error, e, in calculated frequency for points A and B (locations shown in (a)) versus number of retained Bloch modes per selection point, mp. (c) Computational efficiency: ratio of RBME model to full-model calculation times, r, versus number of plane waves used, npw.

8. Summary and conclusions

RBME was presented as an approach for efficient and accurate calculation of band structures for periodic media with any type of lattice symmetry. This infinite media modal analysis approach involves expanding the Bloch solution at all calculation k-points using, in its discrete/discretized form, a selected reduced set of Bloch eigenvectors to form the expansion basis. This basis is selected a priori within the IBZ at high symmetry points determined by the crystal structure and group theory, i.e. the Γ, X, M and R points for the three-dimensional simple cubic lattice. Intermediate k-points along the lines connecting the high symmetry points can be selected in addition, i.e. a set that would now consist of the Γ, Δ, X, Z, M, T, R and Λ points for the three-dimensional simple cubic lattice. The former selection is referred to as a 2-point expansion, and the latter as a 3-point expansion. At each of the reciprocal lattice selection points, a number of Bloch eigenvectors are selected up to the frequency/energy range (or dispersion branch number) of interest for the band structure calculations. Much like modal analysis of finite systems in structural dynamics and other disciplines, the analyst, in general, can decide, as needed, on the number and type of modes to retain. Here, however, the decision is on the number and locations of selection points in k-space and on the number and identity of eigenvectors to retain, in frequency space, for each k-space selection point. As it is common to initially discretize the periodic unit cell and solution field using some choice of basis, e.g. using finite-element or plane-wave bases, RBME is practically a secondary expansion that keeps, and builds on, any favourable attributes a primary expansion approach might exhibit.

Numerical examples were presented using both finite elements and plane waves for the primary expansion. Results presented for two-dimensional plane strain phononic band structure calculations within and along the boundaries of the IBZ are in excellent agreement with those obtained from the full model. Bloch mode shapes and density of state predictions also agree. While 2-point expansion gives accurate band structures, the slightly more expensive 3-point expansion scheme provides higher accuracy, especially at k-points far away from the high symmetry points. Two orders of magnitude in reduction of computation time has been recorded for the proposed approach, and more savings are expected for models of larger size and for calculations based on more refined sampling of the reciprocal lattice space (or IBZ circuit lines). The success of the approach was also demonstrated for two-dimensional photonic and three-dimensional electronic band structure calculations.

Appendix A. Bloch element

An FE treatment of equation (2.13) yields the algebraic eigenvalue problem given in equation (4.3). The global mass and stiffness matrices M and K are assembled from the local matrices me and ke representing a Bloch element, Embedded Image A 1 where δij is the Kronecker delta function, defined as Embedded Image Embedded Image A 2

In equations (A 1) and (A 2), Embedded Image refers to the FE assembly operation in which the local matrices of elements e=1,…,nel are assembled into one global matrix, M and K, respectively. Furthermore, i,j=1,…,nsd, where nsd denotes the number of spatial dimensions, and a,b=1,…,nen, where nen denotes the number of element nodes. Hence, p,q=1,…,nee, where nee denotes the number of element equations and satisfies nee=nsd+nen. With this notation, the following formulae hold: p=nsd(a−1)+i and q=nsd(b−1)+j (Hughes 1987). Recall that D is the reduced elasticity matrix. The matrix Embedded Image is the Bloch-transformed discrete strain operator Embedded Image A 3 For nsd=2 (i.e. a two-dimensional domain), Embedded Image A 4 and for nsd=3 (i.e. a three-dimensional domain), Embedded Image A 5

This formalism allows straightforward implementation, as it only requires application of periodic boundary conditions and a simple modification to the B matrix in existing FE codes.


    • Received November 26, 2008.
    • Accepted June 2, 2009.


View Abstract