text:
# Interaction-driven spin precession in quantum-dot spin valves

## Abstract

We analyze spin-dependent transport through spin valves composed of an interacting quantum dot coupled to two ferromagnetic leads. The spin on the quantum dot and the linear conductance as a function of the relative angle $`\theta `$ of the leads’ magnetization directions is derived to lowest order in the dot-lead coupling strength. Due to the applied bias voltage spin accumulates on the quantum dot, which for finite charging energy experiences a torque, resulting in spin precession. The latter leads to a non-trivial, interaction-dependent, $`\theta `$-dependence of the conductance. In particular, we find that the spin-valve effect is reduced for all $`\theta \pi `$.

Introduction. – The field of spin- or magnetoelectronics has attracted much interest, for both its beautiful fundamental physics and its potential applications. One famous spin-dependent transport phenomenon is the tunnel magnetoresistance (TMR) in a spin-valve geometry, in which two ferromagnetic metals are separated by an insulating layer serving as a tunnel barrier julliere . The transmission through the barrier decreases as the relative angle $`\theta `$ between the magnetizations of two ferromagnets is increased from $`0`$ to $`\pi `$. Within a single-particle picture the $`\theta `$-dependent part of the transmission can be shown slonczewski ; angular  to be proportional to $`\mathrm{cos}\theta `$.

Transport based on tunneling has also been extensively studied in nanostructured devices such as semiconductor quantum dots (QDs) or metallic single-electron transistors. Recently, magnetotransport through those devices has attracted much interest. This includes normal or ferromagnetic metallic islands coupled to ferromagnetic leads ono ; metal-theory  as well as spin-dependent transport from ferromagnets through QDs QD-theory ; Kondo ; QD-exp . Precession of a single magnetic atom spin in an external magnetic field has been detected, but only in the power spectrum of the tunneling current precession\_e ; precession\_t .

In this Letter, we study the effect of strong Coulomb interaction in a single-level QD (or a magnetic impurity impurity ) attached to ferromagnetic leads on the average dot spin and the linear conductance in the weak dot-lead coupling limit ($`\mathrm{\Gamma }k_BT`$, where $`\mathrm{\Gamma }`$ is the intrinsic line width of the dot levels) comment\_0 . We find an interaction-driven spin precession, even in the absence of an external magnetic field. This spin precession is predicted to be clearly visible in the linear conductance as a reduction of the spin-valve effect and a nontrivial $`\theta `$-dependence. For any non-parallel configuration, transport is reduced as compared to the parallel one (spin-valve effect). In the absence of Coulomb interaction, the $`\theta `$-dependence follows simply $`\mathrm{cos}\theta `$ slonczewski ; angular . The presence of a finite charging energy, however, leads to a reduction of the spin-valve effect. This can be understood by the interplay of spin accumulation caused by the bias voltage and an interaction-dependent spin torque due to an effective exchange interaction between the spin in the dot and the leads, which in turn generates spin precession, detectable in the conductance.

The model. – We consider a small QD with one spin-degenerate energy level $`ϵ`$ participating in transport. The left and right lead are magnetized along $`\widehat{𝐧}_\mathrm{L}`$ and $`\widehat{𝐧}_\mathrm{R}`$ (see Fig. 1), with relative angle $`\theta `$. The total Hamiltonian is $`H=H_{\mathrm{dot}}+H_\mathrm{L}+H_\mathrm{R}+H_{\mathrm{T},\mathrm{L}}+H_{\mathrm{T},\mathrm{R}}`$. The first part, $`H_{\mathrm{dot}}=ϵ_\sigma c_\sigma ^{}c_\sigma +Un_{}n_{}`$, describes the dot energy level plus the charging energy $`U`$ for double occupation of the dot. The leads are modeled by $`H_r=_kϵ_ka_{rk}^{}a_{rk}`$ with $`r=\mathrm{L},\mathrm{R}`$. For simplicity, we assume them to be half metallic, i.e., only majority spins have a finite density of states. Tunneling between lead and dot is described by $`H_{\mathrm{T},\mathrm{L}}=t_k(a_{\mathrm{L}k}^{}c_++h.c.)`$, where $`c_+`$ is the Fermi operator for an electron on the QD with spin along $`\widehat{𝐧}_\mathrm{L}`$. It is convenient to quantize the dot spin along the $`z`$-direction in the coordinate system defined by $`\widehat{𝐞}_x=(\widehat{𝐧}_\mathrm{L}+\widehat{𝐧}_\mathrm{R})/|\widehat{𝐧}_\mathrm{L}+\widehat{𝐧}_\mathrm{R}|`$, $`\widehat{𝐞}_y=(\widehat{𝐧}_\mathrm{L}\widehat{𝐧}_\mathrm{R})/|\widehat{𝐧}_\mathrm{L}\widehat{𝐧}_\mathrm{R}|`$, and $`\widehat{𝐞}_z=(\widehat{𝐧}_\mathrm{R}\times \widehat{𝐧}_\mathrm{L})/|\widehat{𝐧}_\mathrm{R}\times \widehat{𝐧}_\mathrm{L}|`$. The tunnel Hamiltonian, then, is

$$H_{\mathrm{T},\mathrm{L}}=\frac{t}{\sqrt{2}}\underset{k}{}(e^{i\theta /4}a_{\mathrm{L}k}^{}c_{}+e^{i\theta /4}a_{\mathrm{L}k}^{}c_{}+h.c.),$$
(1)
and $`H_{\mathrm{T},\mathrm{R}}`$ is the same but with $`\mathrm{L}\mathrm{R}`$ and $`\theta \theta `$. Due to tunneling the dot level acquires a finite width $`\mathrm{\Gamma }=2\pi |t|^2N`$, where $`N=N_\mathrm{L}=N_\mathrm{R}`$ is the density of states of the majority spins of the left and right lead.

With this choice of the quantization axis, the model studied here appears similar to those for Aharonov-Bohm interferometers which contain a single-level QD in each arm AB1 ; AB2 ; AB3 . In the present case, the two “arms” are labeled by the spin $`\sigma =,`$ along the $`z`$-direction, and the Aharonov-Bohm phase due to an enclosed magnetic flux corresponds to the angle $`\theta `$ between the leads’ magnetizations. The limit $`U\mathrm{}`$ of the present model is, in fact, equivalent to the one studied in Sec. IV. C. of Ref. AB2, .

Linear conductance and average spin. – We make use of the analogy between the quantum-dot spin valve and the Aharonov-Bohm interferometer and express the current in terms of Green’s functions of the dot electrons, as shown in Eq. (4.3) of Ref. AB2, . Here, we are interested in first-order transport in $`\mathrm{\Gamma }`$, for which the linear-response conductance $`G^{\mathrm{lin}}=(I/V)|_{V=0}`$ simplifies to comment\_1

$`G^{\mathrm{lin}}`$ $`=`$ $`{\displaystyle \frac{e^2}{h}}\mathrm{\Gamma }{\displaystyle }d\omega \{\mathrm{Im}G_{}^{\mathrm{ret}}(\omega )f^{}(\omega )`$ (2)
$`+\mathrm{sin}{\displaystyle \frac{\theta }{2}}f(\omega ){\displaystyle \frac{G_{}^>(\omega )}{(eV)}}`$
$`+\mathrm{sin}{\displaystyle \frac{\theta }{2}}[1f(\omega )]{\displaystyle \frac{G_{}^<(\omega )}{(eV)}}\}.`$
Here, $`f(\omega )`$ is the Fermi function, $`G_{\sigma \sigma ^{}}(\omega )`$ are the Fourier transforms of the Green’s functions $`G_{\sigma \sigma ^{}}^>(t)=ic_\sigma (t)c_\sigma ^{}^{}(0)`$, $`G_{\sigma \sigma ^{}}^<(t)=ic_\sigma ^{}^{}(0)c_\sigma (t)`$, and $`G_{\sigma \sigma ^{}}^{\mathrm{ret}}`$ is the usual retarded Green’s function. Contributions involving $`\mathrm{Im}G_{}^{\mathrm{ret}}(\omega )`$, $`G_{}^>(\omega )`$, and $`G_{}^<(\omega )`$ are accounted for in the prefactor 2. Since $`\mathrm{\Gamma }`$ already appears explicitly in front of the integral, all Green’s functions are to be taken to zeroth order in $`\mathrm{\Gamma }`$. In this limit, we find $`(1/\pi )\mathrm{Im}G_{}^{\mathrm{ret}}(\omega )=(P_0^0+P_{}^{})\delta (\omega ϵ)+(P_{}^{}+P_d^d)\delta (\omega ϵU)`$, $`G_{}^>(\omega )=2\pi iP_{}^{}\delta (\omega ϵU)`$, and $`G_{}^<(\omega )=2\pi iP_{}^{}\delta (\omega ϵ)`$, where $`P_\chi ^{}^\chi =|\chi ^{}\chi |`$ are elements of the stationary density matrix (to zeroth order in $`\mathrm{\Gamma }`$) of the quantum dot subsystem, with $`\chi ,\chi ^{}=0`$ (empty dot), $`,`$ (singly-occupied dot), and $`d`$ (doubly-occupied dot).

The main task is now to determine these density matrix elements to zeroth order in $`\mathrm{\Gamma }`$. They contain as well the information about the average occupation and spin on the QD. The diagonal matrix elements, $`P_\chi ^\chi `$, are nothing but the probabilities to find the QD in state $`\chi `$, i.e., the dot is empty with probability $`P_0^0`$, singly occupied with $`P_1P_{}^{}+P_{}^{}`$, and doubly occupied with $`P_d^d`$. A finite spin can only emerge for single occupancy. The average spin $`𝐒=(S_x,S_y,S_z)`$ is related to the off-diagonal matrix element $`P_{}^{}`$ via $`S_x=\mathrm{Re}P_{}^{}`$, $`S_y=\mathrm{Im}P_{}^{}`$, and $`S_z=(1/2)(P_{}^{}P_{}^{})`$.

It is remarkable that on the r.h.s of Eq. (2) derivatives of Green’s function with respect to bias voltage $`V`$ appear. As a consequence, not only the equilibrium density matrix elements enter the linear conductance, but also linear corrections in $`V`$ are involved. In equilibrium, $`V=0`$, the density matrix is diagonal with $`P_0^0=1/Z`$, $`P_{}^{}=P_{}^{}=\mathrm{exp}(\beta ϵ)/Z`$, $`P_d^d=\mathrm{exp}[\beta (2ϵ+U)]/Z`$, with $`Z=1+2\mathrm{exp}(\beta ϵ)+\mathrm{exp}[\beta (2ϵ+U)]`$ and $`\beta =1/(k_\mathrm{B}T)`$. As a consequence, the average spin on the QD vanishes at $`V=0`$ comment\_2 . With applied bias voltage, though, a finite spin can accumulate, which yields a finite $`(d𝐒/dV)|_{V=0}`$.

We now determine the needed density matrix elements by using the real-time transport theory developed in Ref. diagrams, . The starting point is the generalized stationary master equation in Liouville space,

$$\left(ϵ_{\chi _1}ϵ_{\chi _2}\right)P_{\chi _2}^{\chi _1}+\underset{\chi _1^{},\chi _2^{}}{}P_{\chi _2^{}}^{\chi _1^{}}\mathrm{\Sigma }_{\chi _2^{},\chi _2}^{\chi _1^{},\chi _1}=0,$$
(3)
where $`\chi _1`$ and $`\chi _2`$ label the QD states, and $`ϵ_{\chi _1}`$ and $`ϵ_{\chi _2}`$ are the corresponding energies. The matrix elements $`P_{\chi _2}^{\chi _1}`$ of the density matrix are connected to each other in Eq. (3) by terms $`\mathrm{\Sigma }_{\chi _2^{},\chi _2}^{\chi _1^{},\chi _1}`$ which can be viewed as generalized transition rates in Liouville space. They are defined as irreducible self-energy parts of the propagation in Liouville space and are represented as diagram blocks on a Keldysh contour. For a detailed derivation of this diagrammatic language, the generalized master equation, and the rules on how to calculate a diagram we refer to Ref. diagrams, .

In the following, we write $`P=\overline{P}+\widehat{P}+\mathrm{}`$ and $`\mathrm{\Sigma }=\overline{\mathrm{\Sigma }}+\widehat{\mathrm{\Sigma }}+\mathrm{}`$, where $`\overline{P}`$ and $`\overline{\mathrm{\Sigma }}`$ denote the equilibrium limit, and $`\widehat{P}`$ and $`\widehat{\mathrm{\Sigma }}`$ are the linear corrections in $`V`$. Using the symmetry of the model, we find the relations $`\overline{P}_\chi ^\chi ^{}=\overline{P}_{\stackrel{~}{\chi }}^{\stackrel{~}{\chi }^{}}`$ and $`\overline{\mathrm{\Sigma }}_{\chi ,\chi ^{\prime \prime }}^{\chi ^{},\chi ^{\prime \prime \prime }}=\overline{\mathrm{\Sigma }}_{\stackrel{~}{\chi },\stackrel{~}{\chi }^{\prime \prime }}^{\stackrel{~}{\chi }^{},\stackrel{~}{\chi }^{\prime \prime \prime }}`$ and $`\widehat{P}_\chi ^\chi ^{}=\widehat{P}_{\stackrel{~}{\chi }}^{\stackrel{~}{\chi }^{}}`$ and $`\widehat{\mathrm{\Sigma }}_{\chi ,\chi ^{\prime \prime }}^{\chi ^{},\chi ^{\prime \prime \prime }}=\widehat{\mathrm{\Sigma }}_{\stackrel{~}{\chi },\stackrel{~}{\chi }^{\prime \prime }}^{\stackrel{~}{\chi }^{},\stackrel{~}{\chi }^{\prime \prime \prime }}`$, where $`\stackrel{~}{\chi }`$ is obtained from $`\chi `$ by the transformation $``$. For transitions from diagonal states in Liouville space to diagonal ones we find $`\widehat{\mathrm{\Sigma }}_{\chi ,\chi ^{}}^{\chi ,\chi ^{}}=0`$. Finally, we drop all $`\mathrm{\Sigma }`$ terms which connect states in Liouville space that are not compatible, at least to lowest order in $`\mathrm{\Gamma }`$. It turns out that it is sufficient to specify Eq. (3) for $`\chi _1=\chi _2=`$ as well as for $`\chi _1=`$, $`\chi _2=`$. For the linear correction in $`V`$ we get

$`0`$ $`=`$ $`\widehat{P}_{}^{}\overline{\mathrm{\Sigma }}_,^,+\widehat{P}_{}^{}\left(\overline{\mathrm{\Sigma }}_,^,\overline{\mathrm{\Sigma }}_,^,\right)`$ (4)
$`0`$ $`=`$ $`\widehat{P}_{}^{}\overline{\mathrm{\Sigma }}_,^,+\widehat{P}_{}^{}\left(\overline{\mathrm{\Sigma }}_,^,\overline{\mathrm{\Sigma }}_,^,\right)`$ (5)
$`+\overline{P}_0^0\widehat{\mathrm{\Sigma }}_{0,}^{0,}+\overline{P}_{}^{}\left(\widehat{\mathrm{\Sigma }}_,^,+\widehat{\mathrm{\Sigma }}_,^,\right)+\overline{P}_d^d\widehat{\mathrm{\Sigma }}_{d,}^{d,}.`$
We evaluate the all necessary diagrams $`\mathrm{\Sigma }_{\chi _2^{},\chi _2}^{\chi _1^{},\chi _1}`$ explicitly. Eventually, we find the solution

$$\widehat{P}_{}^{}=\frac{i}{4}\frac{eV}{k_BT}\overline{P}_1\mathrm{cos}^2\alpha (\theta )\mathrm{sin}\frac{\theta }{2}$$
(6)
and $`\widehat{P}_{}^{}=\widehat{P}_{}^{}=i\widehat{P}_{}^{}\mathrm{tan}\alpha (\theta )`$. We used the definition

$$\mathrm{tan}\alpha (\theta )=\frac{A}{1f(ϵ)+f(ϵ+U)}\mathrm{cos}\frac{\theta }{2},$$
(7)
where $`A=\frac{1}{\pi }\mathrm{Re}\left[\mathrm{\Psi }\left(\frac{1}{2}+i\frac{\beta ϵ}{2\pi }\right)\mathrm{\Psi }\left(\frac{1}{2}+i\frac{\beta (ϵ+U)}{2\pi }\right)\right]`$, and $`\mathrm{\Psi }(x)`$ denotes the digamma function. This means that the spin accumulated in the QD is (for $`eVk_BT`$)

$$|𝐒|=\sqrt{S_y^2+S_z^2}=\frac{eV}{4k_BT}\overline{P}_1\mathrm{cos}\alpha (\theta )\mathrm{sin}\frac{\theta }{2},$$
(8)
with $`\alpha (\theta )`$ being the angle enclosed by the quantum-dot spin and the $`y`$-axis, $`\mathrm{tan}\alpha (\theta )=S_z/S_y`$.

With the result for $`\widehat{P}_{}^{}`$ we are able to obtain the linear conductance. It can be written in the compact form

$$G^{\mathrm{lin}}=G^{\mathrm{lin},\mathrm{max}}\left(1\mathrm{cos}^2\alpha (\theta )\mathrm{sin}^2\frac{\theta }{2}\right).$$
(9)
This equation together with the condition for $`\alpha (\theta )`$, Eq. (7), is the central result of this paper. The conductance is maximal for parallel magnetization, $`\theta =0`$. Its value is $`G^{\mathrm{lin},\mathrm{max}}=(\pi e^2/h)(\mathrm{\Gamma }/k_BT)[1f(ϵ+U)]f(ϵ)[1f(ϵ)+f(ϵ+U)]/[f(ϵ)+1f(ϵ+U)]`$.

We can straightforwardly generalize our theory to allow for arbitrary spin polarization $`p(N_{\mathrm{maj}}N_{\mathrm{min}})/(N_{\mathrm{maj}}+N_{\mathrm{min}})`$ in the leads. In this case, we get $`G^{\mathrm{lin}}/G^{\mathrm{lin},\mathrm{max}}=1p^2\mathrm{cos}^2\alpha _p(\theta )\mathrm{sin}^2(\theta /2)`$ with $`\mathrm{tan}\alpha _p(\theta )=p\mathrm{tan}\alpha _{p=1}(\theta )`$.

A reduction of the spin-valve effect has also been found for hybrid systems of ferromagnets with Luttinger liquids egger  or normal metals brataas . Both its physical origin and the $`\theta `$-dependence of the conductance are different from our proposal, though.

Results and discussion. – It is interesting to analyze how Coulomb interaction in the dot affects the spin accumulation in our model. In the absence of charging energy, $`U=0`$, we find $`A=0`$ and $`\alpha =0`$, and the accumulated spin is along the $`y`$-direction, i.e., along $`\widehat{𝐧}_L\widehat{𝐧}_R`$. A finite charging energy, however, yields a rotation of the quantum-dot spin within the $`yz`$-plane by an angle $`\alpha `$, accompanied by a reduction of the total accumulated spin by $`\mathrm{cos}\alpha `$. The origin of the torque responsible for this spin rotation is an interaction-dependent effective spin exchange of the quantum-dot spin with the lead spins. In the subspace of single dot occupation an effective spin Hamiltonian can be derived from the full model by means of a Schrieffer-Wolff transformation. If we now employ a mean-field picture and replace the lead spin operators by their average value, we end up with the simple effective Hamiltonian $`H^{\mathrm{eff}}=A\mathrm{\Gamma }\mathrm{cos}(\theta /2)S_x`$ for the subspace under consideration. In this effective model, a spin in the $`yz`$-plane experiences a torque and starts to precess, as described by classical Bloch equations. Together with the rate $`\mathrm{\Gamma }_{10}+\mathrm{\Gamma }_{1d}=\mathrm{\Gamma }[1f(ϵ)+f(ϵ+U)]`$ to diminish the $`z`$-component of the quantum-dot spin by changing the dot state from single occupation to an empty or doubly-occupied dot via tunneling, we get the Bloch equation $`dS_z/dt=A\mathrm{\Gamma }\mathrm{cos}(\theta /2)S_y(\mathrm{\Gamma }_{10}+\mathrm{\Gamma }_{1d})S_z`$, from which we can extract the angle $`\alpha (\theta )`$ of the rotated spin in the stationary limit $`dS_z/dt=0`$, and we recover Eq. (7).

We emphasize that both the rates for changing the number of dot electrons and the spin precession are of the same order in $`\mathrm{\Gamma }`$, which is why the the angle $`\alpha `$ is $`\mathrm{\Gamma }`$-independent. The two types of processes correspond to two different kinds of diagrams $`\mathrm{\Sigma }_{\chi _2^{},\chi _2}^{\chi _1^{},\chi _1}`$. In spin-precession terms, all four labels $`\chi _1,\chi _1^{},\chi _2^{},\chi _2`$ represent single occupation, $`,`$. They describe first-order virtual charge fluctuations during which the spin is rotated. In contrast, tunneling rates which change the number of dot electrons are described by diagrams with $`\chi _1=\chi _1^{}=0`$ or $`\chi _1=\chi _1^{}=d`$ or $`\chi _2=\chi _2^{}=0`$, or $`\chi _2=\chi _2^{}=d`$. It is crucial for a consistent theory of first-order transport to include both types of diagrams in Eqs. (4) and (5).

The linear conductance as a function of the level energy $`ϵ`$ is plotted in the upper panel of Fig. 2 for interaction strength $`U/k_BT=10`$ and different values of the angle $`\theta `$. For parallel magnetization, $`\theta =0`$, there are two conductance peaks located near $`ϵ=0`$ and $`ϵ=U`$, respectively. With increasing angle $`\theta `$, transport is more and more suppressed due to the spin-valve effect. However, this suppression is not uniform, as would be in the absence of charging energy. In contrast, the spin-valve effect is less pronounced in the valley between the two peaks, where the dot is dominantly singly occupied, and spin accumulation can occur. As a consequence, the two peaks move towards each other with increasing $`\theta `$.

The differential spin accumulation $`dS/dV`$ in units of $`k_BT`$ is illustrated in the middle panel of Fig. 2. It is clear that single occupation of the dot is required for spin accumulation, i.e., the plotted signal is high in the valley between the two conductance peaks.

As explained above, an effective exchange interaction between quantum-dot spin and spin of the leads yield a rotation of the accumulated spin in the $`yz`$-plane by an angle $`\alpha `$. The lower panel of Fig. 2 depicts the evolution of $`\alpha `$ as a function of the level energy $`ϵ`$. This angle is large in the valley between the conductance peaks, getting close to $`\pm \pi /2`$. A special point is $`ϵ=U/2`$, at which, due to particle-hole symmetry, the effective exchange interaction vanishes. As a consequence, $`\alpha `$ shows a sharp transition from positive to negative values, accompanied with a peak in the accumulated spin.

In the upper panel of Fig. 3 we show the linear conductance as a function of $`\theta `$ for four values of $`ϵ`$. For $`ϵ/k_BT=3`$ and $`1`$, the dot is predominantly empty, and the $`\theta `$-dependence of the conductance is almost harmonic. For $`ϵ/k_BT=1`$ and $`3`$, however, the spin-valve effect is strongly reduced, and conductance is enhanced, except in the regime close to antiparallel magnetization, $`\theta =\pi `$. The enhancement of conductance is related to the fact that the spin precession reduces the angle between the accumulated spin and the magnetization direction of the drain electrode. This is even better illustrated in the lower panel of Fig. 3 which shows the same curves but normalized to conductance at parallel magnetization ($`\theta =0`$). For $`ϵ/k_BT=3`$, the conductance stays almost flat over a broad range, and then establishes the spin-valve effect only in a small region around $`\theta =\pi `$.

Finally, we comment that a finite spin-flip relaxation time $`\tau _{\mathrm{sf}}`$ will reduce the spin-valve effect and limit its observability to $`\mathrm{\Gamma }>\tau _{\mathrm{sf}}^1`$ fujisawa . The main prediction of our theory, the deviation from the $`\mathrm{cos}\theta `$-law, will not be affected by $`\tau _{\mathrm{sf}}`$, as long as a $`\theta `$-dependence is visible.

To summarize, there are two pronounced features in the linear conductance which proves the spin precession proposed in this Letter: (i) the shift of two adjacent conductance peaks towards each other with increasing angle $`\theta `$ (Fig. 2), and (ii) strong deviation from the $`\mathrm{cos}\theta `$ law for the spin-valve effect (Fig. 3).

Acknowledgements. – We thank J. Barnaś, A. Brataas, M. Braun, L. Glazman, S. Maekawa, and G. Schön for discussions. This work was supported by the Deutsche Forschungsgemeinschaft under the Center for Functional Nanostructures and the Emmy-Noether program.


format: text

text:
# Differential Complexes and Numerical Stability

## 1. Introduction

During the twentieth century chain complexes, their exactness properties, and commutative diagrams involving them pervaded many branches of mathematics, most notably algebraic topology and differential geometry. Recently such homological techniques have come to play an important role in a branch of mathematics often thought quite distant from these, numerical analysis. Their most significant applications have been to the design and analysis of numerical methods for the solution of partial differential equations.

Let us consider a general problem, such as a boundary value problem in partial differential equations, as an operator equation: given data $`f`$ in some space $`Y`$ find the solution $`u`$ in some space $`X`$ to the problem $`Lu=f`$. A numerical method discretizes this problem through the construction of an operator $`L_h:X_hY_h`$ and data $`f_hY_h`$ and defines an approximate solution $`u_hX_h`$ by the equation $`L_hu_h=f_h`$. Of course the numerical method is not likely to be of value unless it is consistent which means that $`L_h`$ and $`f_h`$ should be close to $`L`$ and $`f`$ in an appropriate sense.

Before we speak of solving the original problem, numerically or otherwise, we should first confront the question of whether it is well-posed. That is, given $`fY`$, does a unique $`uX`$ exist, and, if so, do small changes of $`f`$ induce small changes in $`Y`$? The analogous questions for the numerical method, whether given $`f_hY_h`$ a unique $`u_hX_h`$ is determined by the discrete equation $`L_hu_h=f_h`$, and whether small changes in $`f_h`$ induce small changes in $`u_h`$, is the question of *stability* of the numerical method. A common paradigm, which can be formalized in many contexts of numerical analysis, is that a method which is consistent and stable is convergent.

Well-posedness is a central issue in the theory of partial differential equations. Of course, we do not expect just any PDE problem to be well-posed. Well-posedness hinges on structure of the problem which may be elusive or delicate. Superficially small changes, for example to the sign of a coefficient or the type of boundary conditions, can certainly destroy well-posedness. The same is true for the stability of numerical methods: it often depends on subtle or elusive properties of the numerical scheme. Usually stability reflects some portion of the structure of the original problem that is captured by the numerical scheme. However in many contexts it is not enough that the numerical scheme be close to the original problem in a quantitative sense for it to inherit stability. That is, it may well happen that a consistent method for a well-posed problem is unstable. In this paper we shall see several examples where the exactness properties of discrete differential complexes and their relation to differential complexes associated with the PDE are crucial tools in establishing the stability of numerical methods. In some cases the homological arguments have served to elucidate or validate methods that had been developed over the preceding decades. In others they have pointed the way to stable discretizations of problems for which none were previously known. They will very likely play a similar role in the eventual solution of some formidable open problems in numerical PDE, especially for problems with significant geometric content, such as in numerical general relativity. As in other branches of mathematics, in numerical analysis differential complexes serve both to encode key structure concisely and to unify considerations from seemingly very different contexts.

In this paper we shall discuss only finite element methods since, of the major classes of numerical methods for PDE, they are the most amenable to rigorous analysis, and have seen the greatest use of differential complexes. But complexes have recently arisen in the study of finite differences, finite volumes, and spectral methods as well.

## 2. Finite element spaces

A finite element space on a domain $`\mathrm{\Omega }`$ is a function space defined piecewise by a certain assembly procedure which we now recall; cf. . For simplicity, here we shall restrict to spaces of piecewise polynomials with respect to a triangulation of an $`n`$-dimensional domain by $`n`$-simplices with $`n=2`$ or $`3`$ (so implicitly we are assuming that $`\mathrm{\Omega }^2`$ is polygonal or $`\mathrm{\Omega }^3`$ is polyhedral). On each simplex $`T`$ we require that there be given a function space of *shape function* $`W_T`$ and a set of *degrees of freedom*, i.e., a set of linear functionals on $`W_T`$ which form a basis for the dual space. Moreover, each degree of freedom is supposed to be associated with some subsimplex of some dimension, i.e., in three dimensions with a vertex, an edge, a face, or the tetrahedron itself. For a subsimplex which is shared by two simplices in the triangulation, we assume that the corresponding functionals are in one-to-one-correspondence. Then the finite element space $`W_h`$ is defined as those functions on $`\mathrm{\Omega }`$ whose restriction to each simplex $`T`$ of the triangulation belongs to $`W_T`$ and for which the corresponding degrees of freedom agree whenever a subsimplex is shared by two simplices.

The simplest example is obtained by choosing $`W_T`$ to be the constant functions and taking as the only degree of freedom on $`T`$ the $`0`$th order moment $`\varphi _T\varphi (x)𝑑x`$ (which we associate with $`T`$ itself). The resulting finite element space is simply the space of piecewise constant functions with respect to the given triangulation. Similarly we could choose $`W_T=_1(T)`$ (by $`_p(T)`$ we denote the space of polynomial functions on $`T`$ of degree at most $`p`$), and take as degrees of freedom the moments of degrees $`0`$ and also those of degree $`1`$, $`\varphi _T\varphi (x)x_i𝑑x`$. Again all the degrees of freedom are associated to $`T`$ itself. This time the finite element space consists of all piecewise linear functions. Of course, the construction extends to higher degrees.

A more common piecewise linear finite element space occurs if we again choose $`W_T=_1(T)`$, but take as degrees of freedom the maps $`\varphi \varphi (v)`$, one associated to each vertex $`v`$. In this case the assembled finite element space consists of all *continuous* piecewise linear functions. More generally we can choose $`W_T=_p(T)`$ for $`p1`$, and associate to each vertex the evaluation degrees of freedom just mentioned, to each edge the moments on the edge of degree at most $`p2`$, to each face the moments on the face of degree at most $`p3`$, and to each tetrahedron the moments of degree at most $`p4`$. The resulting finite element space, called the *Lagrange finite element* of degree $`p`$, consists of all continuous piecewise polynomials of degree at most $`p`$. Figure 1 shows a mesh of a two dimensional domain and a typical function in the space of Lagrange finite elements of degree $`2`$ with respect to this mesh.

Mnemonic diagrams as in Figure 2 are often associated to finite element spaces, depicting a single element $`T`$ and a marker for each degree of freedom.

Next we describe some finite element spaces that can be used to approximate vector-valued functions. For brevity we limit the descriptions to the $`3`$-dimensional case, but supply diagrams in both $`2`$ and $`3`$ dimensions. Of course we may simply take the Cartesian product of three copies of one of the previous spaces. For example, the element diagrams shown on the left of Figure 3 refer to continuous piecewise linear vector fields in two and three dimensions. More interesting spaces are the *face elements* and *edge elements* essentially conceived by Raviart and Thomas  in two dimensions and by Nedelec  in three dimensions. In the lowest order case, the face elements take as shape functions polynomial vector fields of the form $`p(x)=a+bx`$ where $`a^3`$, $`b`$ and $`x=(x_1,x_2,x_3)`$, a 4-dimensional subspace of the 12-dimensional space $`_1(T,^3)`$ of polynomial vector fields of degree at most $`1`$. The degrees of freedom are taken to be the $`0`$th order moments of the normal components on the faces of codimension $`1`$, $`p_fp(x)n_f𝑑x`$ where $`f`$ is a face and $`n_f`$ the unit normal to the face. The element diagram is shown in the middle column of Figure 3. In the lowest order case the edge elements shape functions are polynomial vector fields of the form $`p(x)=a+b\times x`$ where $`a,b^3`$, which form a 6-dimensional subspace of $`_1(T,^3)`$. The degrees of freedom are the $`0`$th order moments over the edges of the component tangent to the edge, $`p_ep(x)t_e𝑑x`$, as indicated on the right of Figure 3.

Each of these spaces can be generalized to arbitrarily high order. For the next higher order face space, the shape functions take the form $`p(x)=a(x)+b(x)x`$ where $`a_1(T,^3)`$ and $`b_1(T)`$ a linear scalar-valued polynomial. This gives a subspace of $`_2(T,^3)`$ of dimension $`15`$, and the degrees of freedom are the moments of degree at most $`1`$ of the normal components on the faces and the moments of degree $`0`$ of all components on the tetrahedron. This element is indicated on the left of Figure 4. For the second lowest order edge space, the shape functions take the form $`p(x)=a(x)+b(x)\times x`$ with $`a,b_1(T,^3)`$, giving a $`20`$-dimensional space. The degrees of freedom are the tangential moments of degree at most $`1`$ on the edges (two per edge) and the tangential moments of degree $`0`$ on the faces (two per face). This element is indicated on the right of Figure 4.

The choice of the shape functions and the degrees of freedom determine the smoothness of the functions belonging to the assembled finite element space. For example, the Lagrange finite element spaces of any degree belong to the Sobolev space $`H^1(\mathrm{\Omega })`$ of $`L^2(\mathrm{\Omega })`$ functions whose distributial first partial derivatives also belong to $`L^2(\mathrm{\Omega })`$ (and even to $`L^infty(\mathrm{\Omega })`$). In fact, the distributional first partial derivative of a continuous piecewise smooth function coincides with its derivative taken piecewise and so belongs to $`L^2`$. Thus the degrees of freedom we imposed in constructing the Lagrange finite elements are sufficient to insure that the assembled finite element space $`W_hH^1(\mathrm{\Omega })`$. In fact more is true: for the Lagrange finite element space with shape function spaces $`W_T=_p(T)`$, we have

$$W_h=\{uH^1(\mathrm{\Omega })|u|_TW_T\text{ for all simplices }T\text{ of the triangulation}\}.$$

This says that, in a sense, the degrees of freedom impose exactly the continuity required to belong to $`H^1`$, no less and no more.

In contrast, the discontinuous piecewise polynomial spaces are subsets of $`L^2(\mathrm{\Omega })`$ but not of $`H^1(\mathrm{\Omega })`$, since their distributional first derivatives involve distributions supported on the interelement boundaries, and so do not belong to $`L^2(\mathrm{\Omega })`$.

For the vector-valued finite elements there are more possibilities. The face and edge spaces contain discontinuous functions, and so are not contained in $`H^1(\mathrm{\Omega },^3)`$. However, for vector fields belonging to one of the face spaces the normal component of the vector field does not jump across interelement boundaries, and this implies, via integration by parts, that the distributional divergence of the function coincides with the divergence taken piecewise. Thus the face spaces belong to $`H(\mathrm{div},\mathrm{\Omega })`$, the space of $`L^2`$ vector fields on $`\mathrm{\Omega }`$ whose divergence belongs to $`L^2`$. Indeed, for these spaces the degrees of freedom impose exactly the continuity of $`H(\mathrm{div})`$, no less or more. For the edge spaces it can be shown that the tangential components of a vector field do not jump across element boundaries, and this implies that the edge functions belong to $`H(\mathrm{curl},\mathrm{\Omega })`$, the space of $`L^2`$ vector fields whose curl belongs to $`L^2`$. Again the degrees of freedom impose exactly the continuity needed for inclusion in $`H(\mathrm{curl})`$.

## 3. Discrete differential complexes

The de Rham complex

$$\begin{array}{ccccccc}^0(\mathrm{\Omega })& \stackrel{d}{}& ^1(\mathrm{\Omega })& \stackrel{d}{}& \mathrm{}& \stackrel{d}{}& ^n(\mathrm{\Omega })0\end{array}$$

is defined for an arbitrary smooth $`n`$-manifold $`\mathrm{\Omega }`$. Here $`^k(\mathrm{\Omega })`$ denotes the space of differential $`k`$-forms on $`\mathrm{\Omega }`$, i.e., for $`\omega ^k(\mathrm{\Omega })`$ and $`x\mathrm{\Omega }`$, $`\omega (x)`$ is an alternating $`k`$-linear map on the tangent space $`T_x\mathrm{\Omega }`$. The operators $`d:^k(\mathrm{\Omega })^{k+1}(\mathrm{\Omega })`$ denote exterior differentiation. This is is a complex in that the composition of two exterior differentiations always vanishes. Moreover, and if the manifold is topologically trivial, then it is exact.

If $`\mathrm{\Omega }`$ is a domain in $`^3`$, then we may identify its tangent space at any point with $`^3`$. Using the Euclidean inner product, the space of linear maps on $`^3`$ may be identified by $`^3`$ as usual, so $`^1(\mathrm{\Omega })`$ may be identified with the space $`C^{\mathrm{}}(\mathrm{\Omega },^3)`$ of smooth vector fields on $`\mathrm{\Omega }`$. Moreover, the space of alternating bilinear maps on $`^3`$ may be identified with $`^3`$ by associating to a vector $`u`$ the alternating bilinear map $`(v,w)det(u|v|w)`$. Thus we have an identification of $`^2(\mathrm{\Omega })`$ with $`^3`$ as well. Finally the only alternating trilinear maps on $`^3`$ are given by multiples of the determinant map $`(u,v,w)cdet(u|v|w)`$, and so we may identify $`^3(\mathrm{\Omega })`$ with $`C^{\mathrm{}}(\mathrm{\Omega })`$. In terms of such *proxy fields*, the de Rham complex becomes

$$\begin{array}{ccccccc}C^{\mathrm{}}(\mathrm{\Omega })& \stackrel{\mathrm{grad}}{}& C^{\mathrm{}}(\mathrm{\Omega },^3)& \stackrel{\mathrm{curl}}{}& C^{\mathrm{}}(\mathrm{\Omega },^3)& \stackrel{\mathrm{div}}{}& C^{\mathrm{}}(\mathrm{\Omega },)0.\end{array}$$
(3.1)
Alternatively we may consider $`L^2`$-based forms and the sequence becomes

$$\begin{array}{ccccccc}H^1(\mathrm{\Omega })& \stackrel{\mathrm{grad}}{}& H(\mathrm{curl},\mathrm{\Omega })& \stackrel{\mathrm{curl}}{}& H(\mathrm{div},\mathrm{\Omega })& \stackrel{\mathrm{div}}{}& L^2(\mathrm{\Omega },)0.\end{array}$$

The finite element spaces constructed above allow us to form discrete analogues of the de Rham complex. Given some triangulation of $`\mathrm{\Omega }^3`$, let $`W_h`$ denote the space of continuous piecewise linear finite elements, $`Q_h`$ the lowest order edge element space, $`S_h`$ the lowest order face element space, and $`V_h`$ the space of piecewise constants. Then $`\mathrm{grad}W_hQ_h`$ (since $`Q_h`$ contains all piecewise constant vector fields belonging to $`H(\mathrm{curl})`$ and the gradient of a continuous piecewise linear is certainly such a function), $`\mathrm{curl}Q_hS_h`$ (since $`S_h`$ contains all piecewise constant vector fields belonging to $`H(\mathrm{curl})`$), and $`\mathrm{div}S_hV_h`$. Thus we have the discrete differential complex

$$\begin{array}{ccccccc}W_h& \stackrel{\mathrm{grad}}{}& Q_h& \stackrel{\mathrm{curl}}{}& S_h& \stackrel{\mathrm{div}}{}& V_h0.\end{array}$$
(3.2)
This differential complex captures the topology of the domain to the same extent as the de Rham complex. In particular, if the domain is topologically trivial, then the sequence is exact.

It is convenient to abbreviate the above statement using the element diagrams introduced earlier. Thus we will say that the following complex is exact:

$$\begin{array}{ccccccccccc}& & \text{}& \stackrel{\mathrm{grad}}{}& \text{}& \stackrel{\mathrm{curl}}{}& \text{}& \stackrel{\mathrm{div}}{}& \text{}& & 0\end{array}$$

By this we mean that if we assemble finite element spaces $`W_h`$, $`Q_h`$, $`S_h`$, and $`V_h`$ using the indicated finite elements and a triangulation of a topologically trivial domain, then the corresponding discrete differential complex (3.2) is exact.

There is another important relationship between the de Rham complex (3.1) and the discrete complex (3.2). The defining degrees of freedom determine projections $`\mathrm{\Pi }_h^W:C^{\mathrm{}}(\mathrm{\Omega })W_h`$, $`\mathrm{\Pi }_h^Q:C^{\mathrm{}}(\mathrm{\Omega },^3)Q_h`$, and so on. In fact $`\mathrm{\Pi }_h^W`$ is just the usual interpolant, $`\mathrm{\Pi }_h^V`$ is the $`L^2`$-projection into the piecewise constants, and the projections $`\mathrm{\Pi }_h^Q`$ and $`\mathrm{\Pi }_h^S`$ onto the edge and face elements are determined by the maintenance of the appropriate moments. It can be checked, based on Stokes theorem, that the following diagram commutes.

$$\begin{array}{ccccccccccc}& & C^{\mathrm{}}(\mathrm{\Omega },)& \stackrel{\mathrm{grad}}{}& C^{\mathrm{}}(\mathrm{\Omega },^3)& \stackrel{\mathrm{curl}}{}& C^{\mathrm{}}(\mathrm{\Omega },^3)& \stackrel{\mathrm{div}}{}& C^{\mathrm{}}(\mathrm{\Omega },)& & 0\\ & & \mathrm{\Pi }_h^W& & \mathrm{\Pi }_h^Q& & \mathrm{\Pi }_h^S& & \mathrm{\Pi }_h^V& & \\ & & W_h& \stackrel{\mathrm{grad}}{}& Q_h& \stackrel{\mathrm{curl}}{}& S_h& \stackrel{\mathrm{div}}{}& V_h& & 0\end{array}$$
(3.3)
The finite element spaces appearing in this diagram, with one degree of freedom for each vertex for $`W_h`$, for each edge for $`Q_h`$, for each face for $`S_h`$, and for each simplex for $`V_h`$, are highly geometrical. In fact, recalling the identifications between fields and differential forms, we may view these spaces as spaces of piecewise smooth differential forms. They were in fact first constructed in this context, without any thought of finite elements or numerical methods, by Whitney . The spaces were reinvented, one-by-one, as finite element spaces in response to the needs of various numerical problems, and the properties which are summarized in the commutative diagram above were slowly rediscovered as needed to analyze the resulting numerical methods. The connection between low order edge and face finite elements and Whitney forms was first realized by Bossavit .

Analogous statements hold for higher order Lagrange, edge, face, and discontinuous finite elements. For example, the following diagram commutes and has exact rows:

$$\begin{array}{ccccccccccc}& & C^{\mathrm{}}(\mathrm{\Omega },)& \stackrel{\mathrm{grad}}{}& C^{\mathrm{}}(\mathrm{\Omega },^3)& \stackrel{\mathrm{curl}}{}& C^{\mathrm{}}(\mathrm{\Omega },^3)& \stackrel{\mathrm{div}}{}& C^{\mathrm{}}(\mathrm{\Omega },)& & 0\\ & & & & & & & & & & \\ & & \text{}& \stackrel{\mathrm{grad}}{}& \text{}& \stackrel{\mathrm{curl}}{}& \text{}& \stackrel{\mathrm{div}}{}& \text{}& & 0\end{array}$$

We shall see many other discrete differential complexes below.

## 4. Stability of Galerkin methods

Consider first the solution of the Dirichlet problem for Poisson’s equation on a domain in $`^n`$:

$$\mathrm{\Delta }u=f\text{ in }\mathrm{\Omega },u=0\text{ on }\mathrm{\Omega }.$$

The solution can be characterized as the minimizer of the energy functional

$$(u):=\frac{1}{2}_\mathrm{\Omega }|\mathrm{grad}u(x)|^2𝑑x_\mathrm{\Omega }f(x)u(x)𝑑x$$

over the Sobolev space $`\stackrel{˚}{H}^1(\mathrm{\Omega })`$ (consisting of $`H^1(\mathrm{\Omega })`$ functions vanishing on $`\mathrm{\Omega }`$), or as the solution of the weak problem: find $`u\stackrel{˚}{H}^1(\mathrm{\Omega })`$ such that

$$_\mathrm{\Omega }\mathrm{grad}u(x)\mathrm{grad}v(x)𝑑x=_\mathrm{\Omega }f(x)v(x)𝑑x\text{for all }v\stackrel{˚}{H}^1(\mathrm{\Omega }).$$

We may define an approximate solution $`u_h`$ by minimizing the Dirichlet integral over a finite dimensional subspace $`W_h`$ of $`\stackrel{˚}{H}^1(\mathrm{\Omega })`$; this is the classical Ritz method. Equivalently, we may use the Galerkin method, in which $`u_hW_h`$ is determined by the equations

$$_\mathrm{\Omega }\mathrm{grad}u_h(x)\mathrm{grad}v(x)𝑑x=_\mathrm{\Omega }f(x)v(x)𝑑x\text{for all }vW_h.$$

After choice of a basis in $`W_h`$ this leads to a system of linear algebraic equations, and $`u_h`$ is computable.

Let $`T_h`$ denote the discrete solution operator $`fu_h`$. Then it is easy to check that $`T_h`$ is bounded as a linear operator from $`H^1(\mathrm{\Omega }):=\stackrel{˚}{H}^1(\mathrm{\Omega })^{}`$ to $`\stackrel{˚}{H}^1(\mathrm{\Omega })`$ by a constant that depends only on the domain $`\mathrm{\Omega }`$ (and, in particular, doesn’t increase if the space $`W_h`$ is enriched). This says that the Galerkin method is *stable*. A consequence is the *quasioptimality estimate*

$$uu_h_{H^1}c\underset{vW_h}{inf}uv_{H^1},$$
(4.4)
for some constant $`c`$ depending only on the domain $`\mathrm{\Omega }`$. Note that there is no restriction on the subspace $`W_h`$ to obtain this estimate. Galerkin’s method for a coercive elliptic problem is always stable and convergence depends only on the approximation properties of the subspace. A natural choice for $`W_h`$ is the Lagrange finite element space of some degree $`p`$ with respect to some regular simplicial mesh of maximal element size $`h`$, in which case Galerkin’s method is a standard finite element method. In this case the right hand side of (4.4) is $`O(h^p)`$ provided that $`u`$ is sufficiently smooth.

Next consider the related eigenvalue problem, which arises in the determination of the fundamental frequencies of a drum. That is, we seek standing wave solutions $`w(x,t)`$ to the wave equation on some bounded domain $`\mathrm{\Omega }^2`$ which vanish on $`\mathrm{\Omega }`$. Assuming that the tension and density of the drum membrane are unity, these solutions have the form $`w(x,t)=\alpha \mathrm{cos}(\sqrt{\lambda }t)u(x)+\beta \mathrm{sin}(\sqrt{\lambda }t)u(x)`$ where $`\alpha `$ and $`\beta `$ are constants and $`u`$ and $`\lambda `$ satisfy the eigenvalue problem

$$\mathrm{\Delta }u=\lambda u\text{ in }\mathrm{\Omega },u=0\text{ on }\mathrm{\Omega }.$$

The eigenvalues $`\lambda `$ form a sequence of positive numbers tending to infinity. The numbers $`\sqrt{\lambda }/(2\pi )`$ are the fundamental frequencies of the drum and the functions $`u`$ give the corresponding fundamental modes.

The eigenvalues and eigenfunctions are characterized variationally as the critical values and critical points of the Rayleigh quotient

$$(u)=\frac{_\mathrm{\Omega }|\mathrm{grad}u(x)|^2𝑑x}{_\mathrm{\Omega }|u(x)|^2𝑑x},$$

defined for nonzero $`u`$ belonging to the Sobolev space $`\stackrel{˚}{H}^1(\mathrm{\Omega })`$. The classical Rayleigh-Ritz method for the approximation of eigenvalue problems determines approximate eigenvalues $`\lambda _h`$ and eigenfunctions $`u_h`$ as the critical values and points of the restriction of $``$ to the nonzero elements of some finite dimensional subspace $`W_h`$ of $`\stackrel{˚}{H}^1(\mathrm{\Omega })`$. Equivalently, we can write the eigenvalue problem in weak form: find $`\lambda `$ and nonzero $`u\stackrel{˚}{H}^1(\mathrm{\Omega })`$ such that

$$_\mathrm{\Omega }\mathrm{grad}u(x)\mathrm{grad}v(x)𝑑x=\lambda _\mathrm{\Omega }u(x)v(x)𝑑x\text{for all }v\stackrel{˚}{H}^1(\mathrm{\Omega }).$$
(4.5)
The Galerkin approximation of the eigenvalue problem, which is equivalent to the Rayleigh-Ritz method, seeks $`\lambda _h`$ and nonzero $`u_hW_h`$ such that

$$_\mathrm{\Omega }\mathrm{grad}u_h(x)\mathrm{grad}v(x)𝑑x=\lambda _h_\mathrm{\Omega }u_h(x)v(x)𝑑x\text{for all }vW_h.$$
(4.6)
We now discuss the convergence of this method. Let $`\lambda `$ denote the $`j`$th eigenvalue of the problem (4.5). In the interest of simplicity we assume that $`\lambda `$ is a simple eigenvalue, so the corresponding eigenfunction $`u`$ is uniquely determined up to sign by the normalization $`u_{H^1}=1`$. Similarly let $`\lambda _h`$ and $`u_h`$ denote the $`j`$th eigenvalue of (4.6). It can then be proved (see, e.g.,  for much more general results) that there exists a constant $`c`$ such that

$$uu_h_{H^1}c\underset{vW_h}{inf}uv_{H^1},|\lambda \lambda _h|cuu_h_{H^1}^2.$$
(4.7)
In short, the eigenfunction approximation is quasioptimal and the eigenvalue error is bounded by the square. Again there is no restriction on the space $`W_h`$.

Figure 5 reports on the computation of the eigenvalues of the Laplacian on an elliptical domain of aspect ratio $`3`$ using Lagrange finite elements of degree $`1`$.

Now consider an analogous problem, the computation of the resonant frequencies of an electromagnetic cavity occupying a region $`\mathrm{\Omega }^3`$. In this case we wish to find standing wave solutions of Maxwell’s equations. If we take the electric permittivity and the magnetic permeability to be unity and assume a lossless cavity with perfectly conducting boundary, we are led to the following eigenvalue problem for the electric field: find nonzero $`E:\mathrm{\Omega }^3`$, $`\lambda `$ such that

$$\mathrm{curl}\mathrm{curl}E=\lambda E,\mathrm{div}E=0\text{ in }\mathrm{\Omega },E\times n=0\text{ on }\mathrm{\Omega }.$$
(4.8)
This is again an elliptic eigenvalue problem and the eigenvalues form a sequence of positive numbers tending to infinity. The divergence constraint is nearly redundant in this eigenvalue problem. Indeed if $`\mathrm{curl}\mathrm{curl}E=\lambda E`$ for $`\lambda >0`$, then $`\mathrm{div}E=\lambda ^1\mathrm{div}\mathrm{curl}\mathrm{curl}E=0`$ since the divergence of a curl vanishes. Thus the eigenvalue problem

$$\mathrm{curl}\mathrm{curl}E=\lambda E\text{ in }\mathrm{\Omega },E\times n=0\text{ on }\mathrm{\Omega },$$
(4.9)
has the same eigenvalues and eigenfunctions as (4.8) except that it also admits $`\lambda =0`$ as an eigenvalue, and the corresponding eigenspace is infinite-dimensional (it contains the gradients of all smooth functions vanishing on the boundary of $`\mathrm{\Omega }`$). The eigenvalues and eigenfunctions are now critical points and values of the Rayleigh quotient

$$(E)=\frac{_\mathrm{\Omega }|\mathrm{curl}E(x)|^2𝑑x}{_\mathrm{\Omega }|E(x)|^2𝑑x},$$

over the space of nonzero fields $`E`$ in $`\stackrel{˚}{H}(\mathrm{curl},\mathrm{\Omega })`$, which is defined to be the space of functions for which both the above integrals exist and are finite and which have vanishing tangential component on the boundary (i.e., $`E\times n=0`$ on $`\mathrm{\Omega }`$).

In Figure 6 we show the result of approximating a two-dimensional version of this eigenvalue problem using the Rayleigh-Ritz method or, equivalently, the Galerkin method with continuous piecewise linear vector fields on $`\mathrm{\Omega }`$ whose tangential components vanish on the boundary (the first element depicted in Figure 3). For $`\mathrm{\Omega }`$ we take a square of side length $`\pi `$, in which case the nonzero eigenvalues are known to be all numbers of the form $`\lambda =m^2+n^2`$ with $`0m,n`$ not both zero, and the corresponding eigenfunctions are $`E=(\mathrm{sin}my,\mathrm{sin}nx)`$. For the mesh pictured, the finite element space has dimension $`290`$. We find that $`73`$ of the $`290`$ computed eigenvalues are between $`0`$ and $`10`$ and that they have no tendency to cluster near the integers $`1,1,2,4,4,5,5,8,9,9`$ which are the exact eigenvalues between $`0`$ and $`10`$. Thus this numerical method is useless: the computed eigenvalues bear no relation to the true eigenvalues! The analogue of (4.7) is surely not true.

If instead we choose the lowest order edge elements as the finite element space (Figure 3, top right), we get very different results. Using the same mesh, the edge finite element space has dimension $`472`$. It turns out that $`145`$ of the computed eigenvalues are zero (to within round-off), and the subsequent eigenvalues are $`0.9998`$, $`0.9999`$, $`2.0023`$, $`3.9968`$, $`4.0013`$, …, i.e., excellent approximations of the exact eigenvalues. See Figure 7.

The striking difference between the behavior of the continuous piecewise linear finite elements and the edge elements for the resonant cavity problem is a question of stability. We shall return to this below, after examining stability in a simpler context.

## 5. Stability of mixed formulations

Consider now the Dirichlet problem

$$\mathrm{div}C\mathrm{grad}u=f\text{ in }\mathrm{\Omega },u=0\text{ on }\mathrm{\Omega },$$

where $`\mathrm{\Omega }`$ is a domain in $`^3`$ and the coefficient $`C`$ is a symmetric positive definite matrix at each point. We may again characterize $`u`$ as a minimizer of the energy functional

$$u\frac{1}{2}C\mathrm{grad}u\mathrm{grad}udxfu𝑑x$$

and use the Ritz method. This procedure is always stable.

However, for some purposes it is preferable to work with the equivalent first order system

$$\sigma =C\mathrm{grad}u,\mathrm{div}\sigma =f.$$
(5.10)
The pair $`(\sigma ,u)`$ is then characterized variationally as the unique critical point of the functional

$$(\sigma ,u)=_\mathrm{\Omega }(\frac{1}{2}C^1\sigma \sigma +u\mathrm{div}\sigma )𝑑x_\mathrm{\Omega }fu𝑑x$$
(5.11)
over $`H(\mathrm{div},\mathrm{\Omega })\times L^2(\mathrm{\Omega })`$. Note that $`(\sigma ,u)`$ is a saddle-point of $``$, not an extremum. Numerical discretizations based on such saddle-point variational principles are called *mixed methods*.

It is worth interpreting the system (5.10) in the language of differential forms, because this brings some insight. The function $`u`$ is a $`0`$-form, and the operation $`u\mathrm{grad}u`$ is just exterior differentiation. The vector field $`\sigma `$ is a proxy for a $`2`$-form and the operation $`\sigma \mathrm{div}\sigma `$ is again exterior differentiation. The loading function $`f`$ is the proxy for a $`3`$-form. Since $`\mathrm{grad}u`$ is the proxy for a $`1`$-form, it must be that the operation on differential forms that corresponds to multiplication by $`C`$ takes $`1`$-forms to $`2`$-forms. In fact, if we untangle the identifications, we find that multiplication by $`C`$ is a Hodge star operation. A Hodge star operator defines an isomorphism of $`^k(\mathrm{\Omega })`$ onto $`^{3k}(\mathrm{\Omega })`$. To determine a particular such operator, we must define an inner product on the tangent space $`^3`$ at each point of $`\mathrm{\Omega }`$. The positive definite matrix $`C`$ does exactly that. Many of the partial differential equations of mathematical physics admit similar interpretations in terms of differential forms. For a discussion of this in the context of discretization, see .

A natural approach to discretization of the mixed variational principle is to choose subspaces $`S_hH(\mathrm{div},\mathrm{\Omega })`$, $`V_hL^2(\mathrm{\Omega })`$ and seek a critical point $`(\sigma _h,u_h)S_h\times V_h`$. This is of course equivalent to a Galerkin method and leads to a system of linear algebraic equations. However in this case, *stability is not automatic*. It can happen that the discrete system is singular, or more commonly, that the norm of the discrete solution operator grows unboundedly as the mesh is refined.

In a fundamental paper, Brezzi  established two conditions that together are sufficient (and essentially necessary) for stability. Brezzi’s theorem applied to a wide class of saddle-point problems, but for simplicity we will state the stability conditions for the saddle-point problem associated to the functional (5.11).

* There exists $`\gamma _1>0`$ such that

$$_\mathrm{\Omega }C^1\tau \tau 𝑑x\gamma _1\tau _{H(\mathrm{div})}^2,$$

for all $`\tau S_h`$ such that $`\mathrm{div}\tau vdx=0`$ for all $`vV_h`$.
* There exists $`\gamma _2>0`$ such that for all $`vV_h`$ there exists nonzero $`\tau S_h`$ satisfying

$$_\mathrm{\Omega }v\mathrm{div}\tau dx\gamma _2v_{L^2}\tau _{H(\mathrm{div})}.$$

Theorem (Brezzi) If the stability conditions (S1) and (S2) are satisfied, then $``$ admits a unique critical point $`(\sigma _h,u_h)`$ over $`S_h\times V_h`$, the solution operator $`f(\sigma _h,u_h)`$ is bounded $`L^2(\mathrm{\Omega })H(\mathrm{div},\mathrm{\Omega })\times L^2(\mathrm{\Omega })`$, and the quasioptimal estimate

$$\sigma \sigma _h_{H(\mathrm{div})}+uu_h_{L^2}c\underset{(\tau ,v)S_h\times V_h}{inf}(\sigma \tau _{H(\mathrm{div})}+uv_{L^2})$$

holds with $`c`$ depending on $`\gamma _1`$ and $`\gamma _2`$.

The stability conditions of Brezzi strongly limit the choice of the mixed finite element spaces $`S_h`$ and $`V_h`$. Condition (S1) is satisfied if the indicated functions $`\tau S_h`$, those whose divergence is orthogonal to $`V_h`$, are in fact divergence-free. (In practice, this is nearly the only way it is satisfied.) This certainly holds if $`\mathrm{div}S_hV_h`$, and so such as inclusion is a common design principle of mixed finite element spaces. On the other hand, condition (S2) is most easily satisfied if $`\mathrm{div}S_hV_h`$, because in this case, given $`vV_h`$, we can choose $`\tau S_h`$ with $`\mathrm{div}\tau =v`$, so $`_\mathrm{\Omega }v\mathrm{div}\tau dx=v_{L^2}^2`$, and the second condition will be satisfied as long as we can insure that $`\tau _{H(\mathrm{div})}\gamma _2^1v_{L^2}`$. In short, we need to know that $`\mathrm{div}`$ maps $`S_h`$ onto $`V_h`$ and that $`\mathrm{div}|_{S_h}`$ admits a bounded one-sided inverse.

The face elements of Raviart-Thomas and Nedelec were designed to satisfy both these conditions. Specifically, let $`S_h`$ again denote the space of face elements of lowest degree (whose element diagram is shown in the middle of the second row of Figure 3), and $`V_h`$ the space of piecewise constants.<sup>1</sup><sup>1</sup>1It may seem odd to seek $`u_h`$ in $`V_h`$, a space of discrete $`3`$-forms, rather than in a space of $`0`$-forms, since $`u`$ is a $`0`$-form. The resolution is through a Hodge star operator, this time formed with respect to the Euclidean inner product on $`^3`$. In the mixed method $`u_h`$ is a discrete $`3`$-form, approximating the image of $`u`$ under this star operator. We know that $`S_hH(\mathrm{div},\mathrm{\Omega })`$ so these elements are admissable for the mixed variational principle. Moreover, we have $`\mathrm{div}S_hV_h`$, so (S1) holds.

To verify (S2), we refer to the commutative diagram (3.3). Given $`vV_h`$, we can solve the Poisson equation $`\mathrm{\Delta }\varphi =v`$ and take $`\sigma =\mathrm{grad}\varphi `$ to obtain a function with $`\mathrm{div}\sigma =v`$ and $`\sigma _{H^1}Cv_{L^2}`$. Now let $`\tau =\mathrm{\Pi }_h^S\sigma S_h`$. Then

$$\mathrm{div}\tau =\mathrm{div}\mathrm{\Pi }_h^S\sigma =\mathrm{\Pi }_h^V\mathrm{div}\sigma =\mathrm{\Pi }_h^Vv=v,$$

where we have used the commutativity and the fact that $`vV_h`$. Moreover $`\tau _{H(\mathrm{div})}c\sigma _{H^1}c^{}v_{L^2}`$, where we used the boundedness of $`\mathrm{\Pi }_h^S`$ on $`H^1(\mathrm{\Omega },^3)`$. This shows that $`\mathrm{div}V_h=S_h`$ and establishes a bound on the one-sided inverse, and so verifies (S2). Of course, the same argument shows the stability of a mixed method based on higher order face elements as well.

Thus we see that the stability of the mixed finite element method depends on the properties of the spaces $`V_h`$ and $`S_h`$ encoded in the rightmost square of the commutative diagram (3.3).

Now let us return to the resonant cavity eigenvalue problem (4.9) for which we explored the Galerkin method: find $`\lambda _h`$, $`0E_hQ_h`$ such that

$$_\mathrm{\Omega }\mathrm{curl}E_h\mathrm{curl}Fdx=\lambda _h_\mathrm{\Omega }E_hF𝑑x\text{for all }FQ_h.$$
(5.12)
We saw that if $`Q_h\stackrel{˚}{H}(\mathrm{curl},\mathrm{\Omega })`$ is taken to be a space of edge elements this method gives good results in that the positive eigenvalues of the discrete problem are good approximations for the positive eigenvalues of the continuous problem. However, the simple choice of Lagrange finite elements did not give good results. We now explain the good performance of the edge elements based on the middle square of the commutative diagram (3.3). Following Boffi et. al  we set $`P_h=\mathrm{curl}Q_h`$ and introduce the following mixed discrete eigenvalue problem: find $`\lambda _h`$, $`0(E_h,p_h)Q_h\times P_h`$ such that

$$_\mathrm{\Omega }E_hF𝑑x+_\mathrm{\Omega }\mathrm{curl}Fp_hdx=0\text{for all }FQ_h,$$
(5.13)
$$_\mathrm{\Omega }\mathrm{curl}E_hqdx=\lambda _h_\mathrm{\Omega }p_hq𝑑x\text{for all }qP_h.$$
(5.14)
It is then easy to verify that if $`\lambda _h`$, $`E_h`$ is a solution to (5.12) with $`\lambda _h>0`$, then $`\lambda _h`$, $`(E_h,\lambda _h^1\mathrm{curl}E_h)`$ is a solution to (5.13), and if $`\lambda _h`$, $`(E_h,p_h)`$ is a solution to (5.13) then $`\lambda _h>0`$ and $`\lambda _h`$, $`E_h`$ is a solution to (5.12). In short, the two problems are equivalent except that the former admits a zero eigenspace which the mixed formulation suppresses. As explained in , the accuracy of the mixed eigenvalue problem (5.13) hinges on the stability of the corresponding mixed source problem. This is a saddle-point problem of the sort studied by Brezzi, and so stability depends on conditions analogous to (S1) and (S2). The proof of these conditions in case $`Q_h`$ is the space of edge elements follows, as in the preceding stability verification, from surjectivity and commutativity properties encoded in the diagram (3.3).

The diagram can also be used to explain the zero eigenspace computed with edge elements. Recall that in the case of the mesh shown in Figure 6, this space had dimension $`145`$. In fact, this eigenspace is simply the null space of the curl operator restricted to $`Q_h`$. Referring again to the commutative diagram (3.3), this is the gradient of the space $`W_h`$ of linear Lagrange elements vanishing on the boundary. Its dimension is therefore exactly the number of interior nodes of the mesh.

## 6. The elasticity complex

Let $`𝕊`$ denote the space of $`3\times 3`$ symmetric matrices. Given a volumetric loading density $`f:\mathrm{\Omega }^3`$, the system of linearized elasticity determines the displacement field $`u:\mathrm{\Omega }^3`$ and the stress field $`\sigma :\mathrm{\Omega }𝕊`$ induced in the elastic domain $`\mathrm{\Omega }`$ by the equations

$$\sigma =Cϵu,\mathrm{div}\sigma =f,$$

together with boundary conditions such as $`u=0`$ on $`\mathrm{\Omega }`$. Here $`ϵu`$ is the symmetric part of the matrix $`\mathrm{grad}u`$, and the elasticity tensor $`C:𝕊𝕊`$ is a symmetric positive definite linear operator describing the particular elastic material, possibly varying from point to point.

The solution $`(\sigma ,u)`$ may be characterized variationally as a saddle-point of the Hellinger-Reissner functional

$$(\sigma ,u)=_\mathrm{\Omega }(\frac{1}{2}C^1\sigma :\sigma +u\mathrm{div}\sigma )dx_\mathrm{\Omega }fudx$$
(6.15)
over $`H(\mathrm{div},\mathrm{\Omega },𝕊)\times L^2(\mathrm{\Omega },^2)`$ (i.e., $`\sigma `$ is sought in the space of square-integrable symmetric-matrix-valued functions whose divergence by rows is square-integrable, and $`u`$ is sought among all square-integrable vector fields).

For a mixed finite element method, we need to specify finite element subspaces $`S_hH(\mathrm{div},\mathrm{\Omega },𝕊)`$ and $`V_hL^2(\mathrm{\Omega },^2)`$ and restrict the domain of the variational problem. Of course the spaces must be carefully designed if the mixed method is to be stable: the analogues of the stability conditions (S1) and (S2) must be satisfied. The functional (6.15) is quite similar in appearance to (5.11) and so it might be expected that the mixed finite elements developed for the latter (the face elements for $`\sigma `$ and discontinuous elements for $`u`$) could be adapted to the case of elasticity. In fact, the requirement of symmetry of the stress tensor and, correspondingly, the replacement of the gradient by the symmetric gradient, changes the structure significantly. Four decades of searching for mixed finite elements for elasticity beginning in the 1960s did not yield any stable elements with polynomial shape functions.

Using discrete differential complexes, R. Winther and the author recently developed the first such elements for elasticity problems in two dimensions . (The three-dimensional case remains open.) For elasticity, the displacement and stress fields cannot be naturally interpreted as differential forms and the relevant differential complex is not the de Rham complex. In three dimensions it is instead the *elasticity complex*:

$$\begin{array}{ccccccccccc}𝕋& & C^{\mathrm{}}(\mathrm{\Omega },^3)& \stackrel{ϵ}{}& C^{\mathrm{}}(\mathrm{\Omega },𝕊)& \stackrel{J}{}& C^{\mathrm{}}(\mathrm{\Omega },𝕊)& \stackrel{\mathrm{div}}{}& C^{\mathrm{}}(\mathrm{\Omega },^3)& & 0.\end{array}$$

Here the operator $`J`$ is a *second order differential operator* which acts on a symmetric matrix field by first replacing each row with its curl and then replacing each column with its curl to obtain another symmetric matrix field. The resolved space $`𝕋`$ is the six-dimensional space of infinitesimal rigid motions, i.e., the same space of linear polynomials $`a+b\times x`$ which arose as the shape functions for the lowest order edge elements. If the domain $`\mathrm{\Omega }`$ is topologically trivial, this complex is exact. Although it involves a second order differential operator, and so looks quite different from the de Rham complex, Eastwood  recently pointed out that it can be derived from the de Rham complex via a general construction known as the Bernstein-Gelfand-Gelfand resolution.

In two dimensions the elasticity complex takes the form

$$\begin{array}{ccccccccc}_1& & C^{\mathrm{}}(\mathrm{\Omega })& \stackrel{J}{}& C^{\mathrm{}}(\mathrm{\Omega },𝕊)& \stackrel{\mathrm{div}}{}& C^{\mathrm{}}(\mathrm{\Omega },^2)& & 0,\end{array}$$

where now the second order differential operator is

$$J=\left(\begin{array}{cc}\frac{^2}{x_2^2}& \frac{^2}{x_1x_2}\\ \\ \frac{^2}{x_1x_2}& \frac{^2}{x_1^2}\end{array}\right).$$

In the lowest order case, the finite elements we introduced in , for which the element diagrams can be seen in Figure 8, use discontinuous piecewise linear vector fields for the displacement field and a piecewise polynomial space which we shall now describe for the stress field. The shape functions on an arbitrary triangle $`T`$ are given by

$$S_T=\{\tau _3(T,𝕊)|\mathrm{div}\tau _1(T,^2)\},$$

which is a $`24`$-dimensional space consisting of all quadratic symmetric matrix fields on $`T`$ together with the divergence-free cubic fields. The degrees of freedom are

* the values of three components of $`\tau (x)`$ at each vertex $`x`$ of $`T`$ (9 degrees of freedom)
* the values of the moments of degree $`0`$ and $`1`$ of the two components of $`\tau n`$ on each edge $`e`$ of $`T`$ (12 degrees of freedom)
* the value of the three components of the moment of degree $`0`$ of $`\tau `$ on $`T`$ (3 degrees of freedom)

Note that these degrees of freedom are enough to ensure continuity of $`\tau n`$ across element faces, and so will furnish a finite element subspace of $`H(\mathrm{div},\mathrm{\Omega },𝕊)`$. The continuity is not however, the minimal needed for inclusion in $`H(\mathrm{div})`$. The degrees of freedom also enforce continuity at the vertices, which is not required for membership in $`H(\mathrm{div})`$. For various reasons, it would be useful to have a mixed finite element for elasticity that does not use vertex degrees of freedom. But, as we remark below, this is not possible if we restrict to polynomial shape functions.

In order to have a well-defined finite element, we must verify that the $`24`$ degrees of freedom form a basis for the dual space of $`S_T`$. We include this verification since it illustrates an aspect of the role of the elasticity complex. Since $`dimS_T=24`$, we need only show that if all the degrees of freedom vanish for some $`\tau S_T`$, then $`\tau =0`$. Now $`\tau n`$ varies cubically along each edge, vanishes at the endpoints, and has vanishing moments of degree $`0`$ and $`1`$. Therefore $`\tau n0`$. Letting $`v=\mathrm{div}\tau `$, a linear vector field on $`T`$, we get by integration by parts that

$$_Tv^2𝑑x=_T\tau :ϵvdx+_T\tau nv𝑑s=0$$

since the integral of $`\tau `$ vanishes as well as $`\tau n`$. Thus $`\tau `$ is divergence-free. In view of the exactness of the elasticity complex, $`\tau =Jq`$ for some smooth function $`q`$. Since all the second partial derivatives of $`q`$ belong to $`_3(T)`$, $`q_5(T)`$. Adjusting by an element of $`_1(T)`$ (the null space of $`J`$), we may take $`q`$ to vanish at the vertices. Now $`^2q/s^2=\tau nn=0`$ on each edge, whence $`q`$ is identically zero on $`T`$. This implies that the gradient of $`q`$ vanishes at the vertices. Since $`^2q/sn=\tau nt=0`$ on each edge (with $`t`$ a unit vector tangent to the edge), we conclude that $`q/n`$ vanishes identically on $`T`$ as well. Since $`q`$ has degree at most $`5`$, it must vanish identically.

Let $`\mathrm{\Pi }_h^S:C^{\mathrm{}}(\mathrm{\Omega },𝕊)S_h`$ denote the projection associated with the supplied degrees of freedom, and $`\mathrm{\Pi }_h^V:C^{\mathrm{}}(\mathrm{\Omega },^2)V_h`$ the $`L^2`$-projection. For any triangle $`T`$, $`\tau C^{\mathrm{}}(\mathrm{\Omega },𝕊)`$, and $`v_1(\mathrm{\Omega },^2)`$, we have

$$_T\mathrm{div}(\tau \mathrm{\Pi }_h^S\tau )v𝑑x=_T(\tau \mathrm{\Pi }_h^S\tau ):ϵvdx+_T(\tau \mathrm{\Pi }_h^S\tau )nv𝑑s.$$

The degrees of freedom entering the definition of $`\mathrm{\Pi }_h^S`$ ensure that the right hand side vanishes, and from this we obtain the commutativity $`\mathrm{div}\mathrm{\Pi }_h^S\tau =\mathrm{\Pi }_h^V\mathrm{div}\tau `$ which is essential for stability. (Actually a technical difficulty arises here, since $`\mathrm{\Pi }_h^S`$ as given is not bounded on $`H^1(\mathrm{\Omega },𝕊)`$. See  for the resolution.) Note that, by their definitions, $`\mathrm{div}S_hV_h`$ and, using the commutativity, we have $`\mathrm{div}S_h=V_h`$, i.e., $`\begin{array}{ccc}S_h& \stackrel{\mathrm{div}}{}& V_h0\end{array}`$ is exact. To complete this to a discrete analogue of the elasticity complex, we define $`Y_h`$ to be the inverse image of $`S_h`$ under $`J`$. Then $`Y_h`$ is exactly the space of $`C^1`$ piecewise quintic polynomials which are $`C^2`$ at the vertices of the meshes. This is in fact a well-known finite element space, called the Hermite quintic or Argyris space, developed for solving $`4`$th order partial differential equations (for which the inclusion in $`H^2(\mathrm{\Omega })`$ and therefore $`C^1`$ continuity is required). The shape functions are $`_5(T)`$ and the $`21`$ degrees of freedom are the values of the function and all its first and second partial derivatives at the vertices and the integrals of the normal derivatives along edges. We then have a *discrete elasticity complex*

$$\begin{array}{ccccccccc}_1& & Y_h& \stackrel{J}{}& S_h& \stackrel{\mathrm{div}}{}& V_h& & 0,\end{array}$$

or, diagrammatically,

$$\begin{array}{ccccccccc}_1& & \text{}& \stackrel{J}{}& \text{}& \stackrel{\mathrm{div}}{}& \text{}& & 0.\end{array}$$

Moreover this sequence is exact and is coupled to the two-dimensional elasticity sequence via a commuting diagram:

$$\begin{array}{ccccccccc}_1& & C^{\mathrm{}}(\mathrm{\Omega })& \stackrel{J}{}& C^{\mathrm{}}(\mathrm{\Omega },𝕊)& \stackrel{\mathrm{div}}{}& C^{\mathrm{}}(\mathrm{\Omega },^3)& & 0\\ & & \mathrm{\Pi }_h^Y& & \mathrm{\Pi }_h^S& & \mathrm{\Pi }_h^V& & \\ _1& & Y_h& \stackrel{J}{}& S_h& \stackrel{\mathrm{div}}{}& V_h& & 0\end{array}$$

The right half of this diagram encodes the information necessary to establish the stability of our mixed finite element method.

The Hermite quintic finite elements arose naturally from our mixed finite elements to complete the commutative diagram. Had they not been long known, we could have used this procedure to devise a finite element space contained in $`H^2(\mathrm{\Omega })`$. In fact, on close scrutiny we can see that any stable mixed finite elements for elasticity with polynomial shape functions will give rise to a finite element space with polynomial shape functions contained in $`H^2(\mathrm{\Omega })`$. However, it is known that such spaces are difficult to construct and complicated. In fact, it can be proved that an $`H^2`$ finite element space must utilize shape functions of degree at least $`5`$ and the first and second partial derivatives at the vertices must be among the degrees of freedom . This helps explain why mixed finite elements for elasticity have proven so hard to devise. In particular, we can rigorously establish the stress elements must involve polynomials of degree $`3`$, and that vertex degrees of freedom are unavoidable.

In addition to the element just described, elements of all greater orders are also introduced in . The elements of next higher order can be seen as the final two elements in this discrete elasticity complex.

$$\begin{array}{ccccccccc}_1& & \text{}& \stackrel{J}{}& \text{}& \stackrel{\mathrm{div}}{}& \text{}& & 0.\end{array}$$

It is also possible to simplify the lowest order element slightly. To do this we reduce the displacement space from piecewise linear vector fields to piecewise rigid motions, and we replace the stress space with the inverse image under the divergence of the reduced displacement space. This leads to a stable element shown in this exact sequence:

$$\begin{array}{ccccccccc}_1& & \text{}& \stackrel{J}{}& \text{}& \stackrel{\mathrm{div}}{}& \text{}& & 0.\end{array}$$

Because of the unavoidable complexity of $`H^2`$ finite elements, practitioners solving $`4`$th order equations often resort to *nonconforming* finite element approximations of $`H^2`$. This means that the finite element space does not belong to $`H^2`$ in that the function or the normal derivative may jump across element boundaries, but the spaces are designed so that jumps are small enough in some sense (e.g., on average). The error analysis is more complicated for nonconforming elements, since in addition to stability and approximation properties of the finite element space, one must analyze the *consistency error* arising from the jumps in the finite elements. In  Winther and the author investigated the the possibility of nonconforming mixed finite elements for elasticity, which, however are stable and convergent, and developed two such elements. These are related to nonconforming $`H^2`$ elements via nonconforming discrete elasticity complexes, two of which are pictured here:

$$\begin{array}{ccccccccc}_1& & \text{}& \stackrel{J}{}& \text{}& \stackrel{\mathrm{div}}{}& \text{}& & 0.\\ \\ _1& & \text{}& \stackrel{J}{}& \text{}& \stackrel{\mathrm{div}}{}& \text{}& & 0.\end{array}$$

In both cases the shape function space for the stress is contained between $`_1(T,𝕊)`$ and $`_2(T,𝕊)`$. The nonconforming $`H^2`$ finite element depicted in these diagrams was developed for certain $`4`$th order problems in . Note the nonconforming mixed elasticity elements are significantly simpler than the conforming ones (and, in particular, don’t require vertex degrees of freedom).