## Abstract

In 2008, researchers at the Hewlett–Packard (HP) laboratories published a paper in *Nature* reporting the development of a new basic circuit element that completes the missing link between charge and flux linkage, which was postulated by Chua in 1971 (Chua 1971 *IEEE Trans. Circuit Theory* **18**, 507–519 (doi:10.1109/TCT.1971.1083337)). The HP memristor is based on a nanometre scale TiO_{2} thin film, containing a— doped region and an undoped region. Further to proposed applications of memristors in artificial biological systems and non-volatile RAM, they also enable reconfigurable nanoelectronics. Moreover, memristors provide new paradigms in application-specific integrated circuits and field programmable gate arrays. A significant reduction in area with an unprecedented memory capacity and device density are the potential advantages of memristors for integrated circuits. This work reviews the memristor and provides mathematical and SPICE models for memristors. Insight into the memristor device is given via recalling the quasi-static expansion of Maxwell’s equations. We also review Chua’s arguments based on electromagnetic theory.

## 1. Introduction

Based on the International Technology Roadmap for Semiconductors (ITRS) report (ITRS 2007), it is predicted that, by 2019, 16 nm half-pitch dynamic random access memory (DRAM) cells will provide a capacity of around 46 GB cm^{−2}, assuming 100 per cent area efficiency. Interestingly, memristors (MRs) promise extremely high capacity of more than 110 and 460 GB cm^{−2} for 10 and 5 nm half-pitch devices, respectively (Williams 2008; Lewis & Lee 2009). In contrast to DRAM memory, MRs provide a non-volatile operation, as is the case for flash memories. Hence, such devices can continue the legacy of Moore’s Law for another decade. Furthermore, inclusion of molecular electronics and computing as an alternative to complementary metal–oxide–semiconductor (CMOS) technologies in the recent ITRS report emphasizes the significant challenges of device scaling (Jones 2009). Moreover, Swaroop *et al.* (1998) demonstrated that the complexity of a synapse, in an analogue very large-scale integration (VLSI) neural network implementation, is minimized by using a device called the programmable metallization cell. This is an ionic programmable resistive device and an MR can be employed to play the same role.

Research on MR applications in various areas of circuit design, alternative materials and spintronic MRs, and especially MR device/circuit modelling, have appeared in the recent literature: (i) SPICE macro-modelling using linear and non-linear drift models (Benderli & Wey 2009; Biolek *et al.* 2009*b*; Chen & Wang 2009; Kavehei *et al.* 2009; Wang, D. *et al.* 2009*a*; Zhang *et al.* 2009), (ii) application of MRs in different circuit configurations and their dynamic behaviour (Li *et al.* 2009; Sun *et al.* 2009*a*,*b*; Wang, D. *et al.* 2009*b*), (iii) application of MR-based dynamic systems to image encryption in Lin & Wang (2009), (iv) a fine resolution programmable resistor using an MR in Shin *et al.* (2009), (v) an MR-based op-amp circuit and filter characteristics of MRs by Yu *et al.* (2009) and Wang, W. *et al.* (2009), respectively, (vi) an MR receiver (MRX) structure for ultra-wide band (UWB) wireless systems (Itoh & Chua 2008; Witrisal 2009), (vii) a memristive system I/O nonlinearity cancellation in Delgado (2009), (viii) the number of required MRs to compute a function by Lehtonen & Laiho (2009), and its digital logic implementation using an MR-based crossbar architecture in Raja & Mourad (2009), (ix) different physical mechanisms to store information in MRs (Driscoll *et al.* 2009; Wang, X. *et al.* 2009), (x) interesting fabricated non-volatile memory using a *flexible MR*, which is an inexpensive and low-power device solution, is also recently reported (Gergel-Hackett *et al.* 2009) and (xi) using the *memductance* concept to develop an equivalent circuit diagram of a transmission line has been carried out by Mallégol (2003, appendix 2).

There are still many problems associated with the MR device level. For instance, it is not clear which leakage mechanisms are associated with the device. The flexible MR is able to retain its non-volatile characteristic for up to 14 days or, equivalently, up to 4000 flexes (Gergel-Hackett *et al.* 2009). Thus, the non-volatility feature eventually vanishes after a short period. Strukov & Williams (2009) also investigated this particular feature as a ratio of volatility to switching time.

One disadvantage of using MRs is switching speed. The volatility-to-switching speed ratio for MR cells in the Hewlett–Packard (HP) crossbar structure is around 10^{3} (Strukov & Williams 2009), which is much lower than the ratio for DRAM cells, 10^{6} (ITRS 2007); therefore, the switching speed of MRs is far behind DRAM. However, unlike DRAM, resistive random access memory (RRAM) is non-volatile. In terms of yield, DRAM and RRAM are almost equal (Lewis & Lee 2009). Generally, endurance becomes very important once we note that the DRAM cells must be refreshed at least every 16 ms, which means at least 10^{10} write cycles in their lifetime (Lewis & Lee 2009). Unfortunately, MRs are far behind DRAM in terms of endurance, but the HP team is confident that there is no functional limitation against improvement of MRs (Williams 2008). Finally, another advantage of RRAM is readability. Readability refers to the ability of a memory cell to report its state. This noise immunity in DRAM is weak because each DRAM cell stores a very small amount of charge particularly at nanometer dimensions, while in RRAM cells, e.g. MRs, it depends on the difference between the on and off state resistances (Lewis & Lee 2009). This difference was reported up to one order of magnitude (Williams 2008).

It is interesting to note that there are devices with similar behaviour to an MR (e.g. Buot & Rajagopal 1994; Beck *et al.* 2000; Krieger *et al.* 2001; Liu *et al.* 2006; Waser & Masakazu 2007; Ignatiev *et al.* 2008; Tulina *et al.* 2008), but the HP scientists were the only group that found the link between their work and the missing MR postulated by Chua. Moreover, it should be noted that physically realized MRs must meet the mathematical requirements of MR devices or memristive systems that are discussed by Kang (1975) and Chua & Kang (1976). For instance, the hysteresis loop must have a double-valued bow-tie trajectory. However, for example, in Krieger *et al.* (2001), the hysteresis loop shows more than two values for some applied voltage values. Beck *et al.* (2000) demonstrated one of the perfect examples of MR devices. Kim, K. H. *et al.* (2009) recently introduced a multi-level one-time programmable (OTP) oxide diode for crossbar memories. They focused on an OTP structure that basically uses one diode and one resistor, 1D–1R, as obtaining a stable device for handling multiple programming and erasing processes is much more difficult than an OTP device. In terms of functionality, OTP devices are very similar to memristive elements, but in terms of flexibility MRs are able to handle multiple programming and erasing processes.

This paper focuses on the MR device and reviews its device-level properties. Although the MR as a device is new, it was conceptually postulated by Chua (1971). Chua predicted that an MR could be realized in the form of a purely dissipative device as a fourth fundamental circuit element. Thirty-seven years later, Stan Williams and his group in the Information and Quantum Systems Laboratory at HP developed the MR in device form (Strukov *et al.* 2008).

In §2, we review the MR and its characteristics as a nano-switch, which was developed by HP (Strukov *et al.* 2008), and we review its properties based on the early mathematical models. We introduce a new model using a parameter we call the *resistance modulation index* (RMI). Owing to the significance of ionic drift that plays the most important role in the memristive effect, this section is divided into two: (i) a linear and (ii) a nonlinear drift model. Section 3 presents a preliminary SPICE macro-model of the MR and different types of circuit elements in combination with a proposed MR macro-model. Section 4 describes an interpretation of the MR based on electromagnetic theory by recalling Maxwell’s equations. Finally, we summarize this review in §5.

## 2. MR device properties

Traditionally, there are only three fundamental passive circuit elements: capacitors, resistors and inductors, discovered in 1745, 1827 and 1831, respectively. However, one can set up five different mathematical relations between the four fundamental circuit variables: electric current *i*, voltage *v*, electric charge *q* and magnetic flux *φ*. For *linear* elements, *f*(*v*,*i*)=0, *f*(*v*,*q*)=0 where *i*=(*dq*/*dt*) (*q*=*Cv*), and *f*(*i*,*φ*)=0, where *v*=(*dφ*/*dt*) (*φ*=Li), indicate linear resistors, capacitors and inductors, respectively.

Chua (1971) proposed that there should be a fourth fundamental passive circuit element to set up a mathematical relationship between *q* and *φ*, which he named the *MR* (a portmanteau of *memory* and *resistor*). Chua predicted that a class of MRs might be realizable in the form of a pure solid-state device without an internal power supply.

Williams (2008), at HP, announced the first fabricated MR device (Strukov *et al.* 2008). However, a resistor with memory is not a new thing. If taking the example of non-volatile memory, it dates back to 1960 when Bernard Widrow introduced a new circuit element named the *memistor* (Widrow *et al.* 1960). The reason for choosing the name of memistor is exactly the same as the MR, a resistor with memory. The memistor has three terminals and its resistance is controlled by the time integral of a control current signal. This means that the resistance of the memistor is controlled by charge. Widrow devised the memistor as an electrolytic memory element to form a basic structure for a neural circuit architecture called ADALINE (ADAptive LInear NEuron), which was introduced by him and his postgraduate student, Marcian Edward ‘Ted’ Hoff (Widrow *et al.* 1960). However, the memistor is not exactly what researchers were seeking at the nanoscale. It is just a charge-controlled three-terminal (transistor) device. In addition, a two-terminal nano-device can be fabricated without nanoscale alignment, which is an advantage over three-terminal nano-devices (Lehtonen & Laiho 2009). Furthermore, the electrochemical memistors could not meet the requirement for the emerging trend of solid-state integrated circuitry.

Thirty years later, Thakoor *et al.* (1990) introduced a solid-state thin-film tungsten-oxide-based, non-volatile memory. The concept is almost the same as the HP MRs. Their solid-state memistor is electrically re-programmable, it has variable resistance and it is an analogue synaptic memory connection that can be used in electrical neural networks. They claimed that the resistance of the device could be stabilized at any value, between 100 kΩ and 1 GΩ. This solid-state memistor has a thick, 60–80 nm, layer of silicon dioxide that achieves non-volatility over a period of several months. This thick electron blocking layer, however, causes a considerable reduction in the applied electric field. As a consequence, they reported very high programming voltages (±25–30) and very long (minutes to hours) switching times.

In the 1960s, the very first report on the hysteresis behaviour of a current–voltage curve was published by Simmons (1963), which is known as the Simmons tunnelling theory. The Simmons theory generally characterizes the tunnelling current in a metal–insulator–metal (MIM) system. A variable resistance with hysteresis was also published by Simmons & Verderber (1967). They introduced a thin-film (20–300 nm) silicon dioxide doped with gold ions sandwiched between two 200 nm metal contacts. Thus, overall, the device is a MIM system, using aluminium metal contacts. Interestingly, their system demonstrates a sinh function behaviour as also reported in Yang *et al.* (2008) for the HP MR. It is also reported that the switching from high to low impedance takes about 100 ns. As their device is modelled as an energy-storage element, it cannot be an MR, because MRs remember the total charge passing through the port and do not store charge (Oster 1974).

The sinh resistance behaviour of MRs can be used to compensate the linearity of analogue circuits. In Varghese & Gandhi (2009), an MR-based amplifier was proposed using the behaviour of a sinh resistor (Tavakoli & Sarpeshkar 2005) as an MR element. In Shin *et al.* (2009), such a structure was also introduced as an example without mentioning the sinh resistance behaviour of MRs. The idea in Tavakoli & Sarpeshkar (2005) is to characterize a sinh function-type circuit that can be used to linearize a tanh function-type circuit behaviour, e.g. a CMOS differential amplifier. A theoretical analysis shows that a sinh function significantly reduces the third harmonic coefficient and as a consequence reduces nonlinearity of the circuit.

Chua (1971) mathematically predicted that there is a solid-state device which defines the missing relationship between four basic variables. Recall that a resistor establishes a relation between voltage and current, a capacitor establishes a charge–voltage relation and an inductor establishes a current–flux relationship, as illustrated in figure 1. Notice that we are specifically discussing nonlinear circuit elements here.

Consequently, *φ*=*f*_{M}(*q*) or *q*=*g*_{M}(*φ*) defines a charge-controlled (flux-controlled) MR. Then, (*dφ*/*dt*)=(*df*_{M}(*q*)/*dq*)(*dq*/*dt*) or (*dq*/*dt*)=(*dg*_{M}(*φ*)/*dφ*)(*dφ*/*dt*), which implies *v*(*t*)=(*df*_{M}(*q*)/*dq*)*i*(*t*) or *i*(*t*)=(*dg*_{M}(*φ*)/*dφ*)*v*(*t*). Note that *M*(*q*)=(*df*_{M}(*q*)/*dq*) for a charge-controlled MR and *W*(*φ*)=(*dg*_{M}(*φ*)/*dφ*) for a flux-controlled memristor, where *M*(*q*) is the incremental memristance and *W*(*φ*) is the incremental memductance, where *M*(*q*) is in the units of ohms and the units of *W*(*φ*) is in siemens (Chua 1971).

Memristance, *M*(*q*), is the slope of the *φ*–*q* curve, which for a simple case is a piecewise *φ*–*q* curve with two different slopes. Thus, there are two different values for *M*(*q*), which is exactly what is needed in binary logic. For detailed information regarding typical *φ*–*q* curves, the reader is referred to Chua (1971).

It is also obvious that if *M*(*q*)≥0, then instantaneous power dissipated by an MR, *p*(*i*)=*M*(*q*)(*i*(*t*))^{2}, is always positive so the MR is a passive device. Thus, the *φ*–*q* curve is a monotonically increasing function. This feature is exactly what is observed in the HP MR device (Strukov *et al.* 2008). Some other properties of the MR such as zero crossing between current and voltage signals can be found in Chua (1971) and Chua & Kang (1976). The most important feature of an MR is its pinched hysteresis loop *v*–*i* characteristic (Chua & Kang 1976). A very simple consequence of this property and *M*(*q*)≥0 is that such a device is purely dissipative like a resistor.

Another important property of an MR is its excitation frequency. It has been shown that the pinched hysteresis loop is shrunk by increasing the excitation frequency (Chua & Kang 1976). In fact, when the frequency increases towards infinity, the MR acts as a linear resistor (Chua & Kang 1976).

Interestingly enough, an attractive property of the HP MR (Strukov *et al.* 2008), which is exclusively based on its fabrication process, can be deduced from the HP MR simple mathematical model (Strukov *et al.* 2008) and is given by
2.1
where *β* has the dimensions of magnetic flux *φ*(*t*). Here, *β*=*D*^{2}/*μ*_{D} in units of *sV* ≡*Wb*, where *μ*_{D} is the average drift mobility in units of cm^{2} s^{−1} V^{−1} and *D* is the film (titanium dioxide; TiO_{2}) thickness. Note that *R*_{OFF} and *R*_{ON} are simply the ‘on’ and ‘off’ state resistances. Also, *q*(*t*) defines the total charge passing through the MR device in a time window, *t*−*t*_{0}. Notice that the MR has an internal state (Chua & Kang 1976). Furthermore, as stated in Oster (1974), , the state variable in a charge-controlled MR gives the charge passing through the device and does not behave as a storage charge as in a capacitor, as incorrectly reported in some works, e.g. Simmons & Verderber (1967). This concept is very important from two points of view. First of all, an MR is not an energy-storage element. Second, this shows that the MR is not merely a nonlinear resistor, but is a nonlinear resistor with charge as a state variable (Oster 1974).

Five years after Chua’s paper on the MR (Chua 1971), he and his graduate student, Kang, published a paper defining a much broader class of systems, named *memristive systems*. From the memristive systems viewpoint, a generalized definition of an MR is *v*(*t*)=*R*(*w*)*i*(*t*), where *w* defines the internal state of the system and (*dw*/*dt*)=*f*(*w*,*i*) (Chua & Kang 1976). Based on this definition, an MR is a special case of a memristive system.

The HP MR (Strukov *et al.* 2008) can be defined in terms of memristive systems. It exploits a very thin-film TiO_{2} sandwiched between two platinum (Pt) contacts and one side of the TiO_{2} is doped with oxygen vacancies, which are positively charged ions. Therefore, there is a TiO_{2} junction where one side is doped and the other is undoped. Such a doping process results in two different resistances: one is a high resistance (undoped) and the other is a low resistance (doped). Hence, HP intentionally established a device that is illustrated in figure 2. The internal state variable, *w*, is also the variable length of the doped region. Roughly speaking, when we have a non-conductive channel and when we have a conductive channel. The HP MR switching mechanism is further discussed by Yang *et al.* (2008).

### (a) Linear drift model

The MR’s state equation is at the heart of the HP memristive system mechanism (Wang 2008; Joglekar & Wolf 2009). Let us assume a uniform electric field across the device; thus there is a linear relationship between drift–diffusion velocity and the net electric field (Blanc & Staebler 1971). Therefore, the state equation is 2.2

Integrating equation (2.2) gives
where *w*(*t*_{0}) is the initial length of *w*. Hence, the speed of drift under a uniform electric field across the device is given by *v*_{D}=(*dw*(*t*)/*dt*). In a uniform field, we have *D*=*v*_{D}×*t*. In this case, *Q*_{D}=*i*×*t* also defines the amount of required charge to move the boundary from *w*(*t*_{0}), where , to distance *w*(*t*_{D}), where . Therefore, *Q*_{D}=(*β*/*R*_{ON}), so
2.3
If *x*(*t*)=*w*(*t*)/*D*, then
2.4
where *q*(*t*)/*Q*_{D} describes the amount of charge that is passed through the channel over the required charge for a conductive channel.

Using Strukov *et al.* (2008), we have
2.5

By inserting *x*(*t*)=*w*(*t*)/*D*, equation (2.5) can be rewritten as
2.6

Now assume that *q*(*t*_{0})=0, then *w*(*t*)=*w*(*t*_{0})≠0, and, from equation (2.6),
2.7
where *r*=*R*_{OFF}/*R*_{ON} and *M*_{0} is the memristance value at *t*_{0}. Consequently, the following equation gives the memristance at time *t*:
2.8
where Δ*R*=*R*_{OFF}−*R*_{ON}. When *R*_{OFF}≫*R*_{ON}, *M*_{0}≈*R*_{OFF} and equation (2.1) can be derived from equation (2.8).

Substituting equation (2.8) into *v*(*t*)=*M*(*q*)*i*(*t*), when *i*(*t*)=*dq*(*t*)/*dt*, we have
2.9

Recalling that *M*(*q*)=*dφ*(*q*)/*dq*, the solution is
2.10

Using Δ*R*≈*M*_{0}≈*R*_{OFF}, equation (2.10) becomes
2.11

Consequently, using equation (2.4) if *Q*_{D}=*D*^{2}/*μ*_{D}*R*_{ON}, so the internal state of the MR is
2.12

The current–voltage relationship in this case is 2.13

In equation (2.13), the inverse square relationship between memristance and TiO_{2} thickness, *D*, shows that, for smaller values of *D*, the memristance shows improved characteristics, and for this reason the MR imposes a small value of *D*.

In equations (2.10)–(2.13), the only term that significantly increases the role of *φ*(*t*) is lower *Q*_{D}. This shows that at the micrometre scale 1/*R*_{OFF}*Q*_{D}=1/*rβ*=*μ*_{D}/*rD*^{2} is negligible and the relationship between current and voltage reduces to a resistor equation.

Substituting *β*=*D*^{2}/*μ*_{D}, which has the same units as magnetic flux, into equation (2.13), and considering *c*(*t*)=*φ*(*t*)/*β*=*μ*_{D}*φ*(*t*)/*D*^{2} as a normalized variable, we obtain
2.14
where is what we call the (RMI) (Kavehei *et al.* 2009).

Owing to the extremely uncertain nature of nanotechnologies, a variability-aware modelling approach should always be considered. Two well-known solutions to analyse a memristive system were investigated by Chen & Wang (2009): (i) a Monte Carlo simulation for evaluating the (almost) complete statistical behaviour of a device and (ii) Corner analysis. Considering the trade-off between time complexity and accuracy between these two approaches shows the importance of using a simple and reasonably accurate model, because finding the real corners is highly dependent on the model accuracy. The RMI could be one of the device parameters in the model extraction phase, so it would help to provide a simple model with fewer parameters.

Joglekar & Wolf (2009) clarified the behaviour of two MRs in series. As shown in figure 3, they labelled the polarity of an MR as *η*=±1, where *η*=+1 signifies that *w*(*t*) increases with positive voltage or, in other words, the doped region in MR is expanding. If the doped region, *w*(*t*), is shrinking with positive voltage, then *η*=−1. In other words, reversing the voltage source terminals implies MR polarity switching. In figure 3*a*, the doped regions are simultaneously shrunk so the overall memristive effect is retained. In figure 3*b*, however, the overall memristive effect is suppressed (Joglekar & Wolf 2009).

Using the MR polarity effect and equation (2.2), we thus obtain 2.15 Then, with a similar approach, we have 2.16 There is also no phase shift between current and voltage signals, which implies that the hysteresis loop always crosses the origin as demonstrated in figure 4. For further investigation, if a voltage, , is applied across the device, the magnetic flux would be . The inverse relation between flux and frequency shows that, at very high frequencies, there is only a linear resistor.

Figure 4 demonstrates equation (2.16) in MATLAB for five different frequencies, where *ω*_{0}=2*πv*_{0}/*β*, using the Strukov *et al.* (2008) parameter values, *D*=10 nm, *μ*_{D}=10^{−10} cm^{2} s^{−1} V^{−1}, *R*_{ON}=100 Ω, *R*_{OFF}=16 kΩ, *v*_{0}=1 V and *η*=−1. The simulation with the same parameter values is shown in figure 5.

Using different parameters causes a large difference in the hysteresis and MR characteristics. In Witrisal (2009), a new solution for UWB signals using MR devices was introduced. Applying the parameter values given by Witrisal (2009) results in a significant difference in the value of *ω*_{0}. Substituting parameter values given in the paper gives *ω*_{0}≈4 GHz instead of *ω*_{0}≈50 kHz, using *t*_{0}≈0.1 ms and the parameter values from Strukov *et al.* (2008). The new parameter values are *D*=3 nm, , *R*_{ON}=100 Ω, *R*_{OFF}=10 kΩ and *v*_{0}=0.2 V. Using these parameters shows that, even though the highest and lowest memristance ratio in the last case, Figure 4, is around 2, here the ratio is approximately equal to 120.

### (b) Nonlinear drift model

The electrical behaviour of the MR is directly related to the shift in the boundary between doped and undoped regions, and the effectively variable total resistance of the device when an electrical field is applied. Basically, a few volts supply voltage across a very thin film, e.g. 10 nm, causing a large electric field. For instance, it could be more than 10^{6} V cm^{−1}, which results in a fast and significant reduction in the energy barrier (Blanc & Staebler 1971). Therefore, it is reasonable to expect a high nonlinearity in ionic drift–diffusion (Waser *et al.* 2009).

One significant problem with the linear drift assumption is the boundaries. The linear drift model, equation (2.2), suffers from problematic boundary effects. Basically, at the nanoscale, a few volts cause a large electric field that leads to a significant nonlinearity in ionic transport (Strukov *et al.* 2008). A few attempts have been carried out so far to consider this nonlinearity in the state equation (Strukov *et al.* 2008; Benderli & Wey 2009; Biolek *et al.* 2009*b*; Strukov & Williams 2009). All of them proposed a simple *window function*, *F*(*ξ*), which is multiplied by the right-hand side of equation (2.2). In general, *ξ* could be a variable vector, e.g. *ξ*=(*w*,*i*) where *w* and *i* are the MR’s state variable and current, respectively.

In general, the window function can be multiplied by the right-hand side of the state variable equation, equation (2.2),
2.17
where *x*(*t*)=*w*(*t*)/*D* is the normalized form of the state variable. The window function makes a large difference between the model with linear and that with nonlinear drift assumptions at the boundaries. Figure 6 shows such a condition considering a nonlinear drift assumption at the critical, or boundary, states. In other words, when a linear drift model is used, simulations should consider the boundaries and all constraints on initial current, initial voltage, maximum and minimum *w*, etc. These constraints cause a large difference in output between linear and nonlinear drift assumptions. For example, it is impossible to achieve such a realistic curve, as in figure 6, using the linear drift modelling approach.

In Strukov *et al.* (2008), the window function is a function of the state variable and it is defined as *F*(*w*)=*w*(1−*w*)/*D*^{2}. The boundary conditions in this case are *F*(0)=0 and *F*(*D*)=((1−*D*)/*D*)≈0. It meets the essential boundary condition , except there is a slight difference when . The problem with this boundary assumption is when an MR is driven to the terminal states, and , then , so no external field can change the MR state (Biolek *et al.* 2009*b*). This is a fundamental problem of the window function. The second problem of the window function is it assumes that the MR remembers the amount of charge that passed through the device. Basically, this is a direct result of the state equation, equation (2.2) (Biolek *et al.* 2009*b*). However, it seems that the device remembers the position of the boundary between the two regions, and not the amount of charge.

In Benderli & Wey (2009), another window function has been proposed that is slightly different from that in Strukov *et al.* (2008). This window function, *F*(*w*)=*w*(*D*−*w*)/*D*^{2}, approaches zero when , and when then . Therefore, this window function meets both the boundary conditions. In fact, the second window function imitates the first function when we consider *x*=*w*/*D* instead of *w*. In addition, there is another problem associated with these two window functions, namely the modelling of approximate linear ionic drift when 0<*w*<*D*. Both of the window functions approximate nonlinear behaviour when the MR is not in its terminal states, *w*=0 or *w*=*D*. This problem is addressed by Joglekar & Wolf (2009), who propose an interesting window function to address the nonlinear and approximately linear ionic drift behaviour at the boundaries and when 0<*w*<*D*, respectively. Nonlinearity (or linearity) of their function can be controlled with a second parameter, which we call the *control parameter*, *p*. Their window function is *F*(*x*)=1−(2*x*−1)^{2p}, where *x*=*w*/*D* and *p* is a positive integer. Figure 7*a* demonstrates the function behaviour for different 2≤*p*≤10 values. This model considers a simple boundary condition, *F*(0)=*F*(1)=0. As demonstrated, when *p*≥4, the state variable equation is an approximation of the linear drift assumption, *F*(0<*x*<1)≈1.

The most important problem associated with this model is revealed at the boundaries. Based on this model, when an MR is at the terminal states, no external stimulus can change its state. Biolek *et al.* (2009*b*) tackle this problem with a new window function that depends on *x*, *p* and MR current, *i*. Basically, *x* and *p* are playing the same role in their model and the only new idea is using current as an extra parameter. The window function is *F*(*x*)=1−(*x*−*sgn*(−*i*))^{2p}, where *i* is MR current and *sgn*(*i*)=1 when *i*≥0, and *sgn*(*i*)=0 when *i*<0. As a matter of fact, when the current is positive, the doped region length, *w*, is expanding. Figure 7*b* illustrates the window function behaviour.

All window functions suffer from a serious problem. They are only dependent on the state variable, *x*, which implies that the MR remembers the entire charge that is passing through it. Moreover, based on the general definition of the time-invariant MR’s state equation and current–voltage relation, or , *f* and *g* must be *continuous* *n*-dimensional vector functions (Kang 1975, ch. 2). However, the last window function, *F*(*x*)=1−(*x*−*sgn*(−*i*))^{2p}, does not provide the continuity condition at the boundaries, or . Biolek *et al.* (2009*b*) did not use the window function in their recent publication (Biolek *et al.* 2009*a*).

In Strukov & Williams (2009), the overall drift velocity is identified with one linear equation and one highly nonlinear equation, *υ*=*μE*, when *E*≪*E*_{0} and , for *E*∼*E*_{0}, where *υ* is the average drift velocity, *E* is an applied electric field, *μ* is the mobility and *E*_{0} is the characteristic field for a particular mobile ion. The value of *E*_{0} is typically about 1 MV cm^{−1} at room temperature (Strukov & Williams 2009). This equation shows that a very high electric field is needed for exponential ion transport. They also showed that the high electric field is lower than the critical field for dielectric breakdown of the oxide. Reviewing the available window functions indicates that there is a need for an appropriate model that can define MR states for strongly nonlinear behaviour where *E*∼*E*_{0}.

In addition to the weakness of nonlinear modelling in the original HP model, there are some other drawbacks that show that the connection between physics and electronic behaviour was not well established. Moreover, the currently available electronic models for MRs followed the exact pathway of the HP modelling, which is mostly due to the fact that the underlying physical conduction mechanism is not fully clear yet (Kim, T. H. *et al.* 2009). One weakness is that the HP model does not deliver any insight about capacitive effects, which are naturally associated with MRs. These capacitive effects will later be explained in terms of a memcapacitive effect in a class of circuit elements with memory. In Kim, T. H. *et al.* (2009), the MR behaviour was realized using an infinite number of crystalline magnetite (Fe_{3}O_{4}) nanoparticles. The device behaviour combines both memristive (time-varying resistance) and memcapacitive (time-varying capacitance) effects, which deliver a better model for the nonlinear properties. Their model description for the current–voltage relationship is given as
2.18
where *R*(*x*,*t*) and *C*(*x*,*t*) are the time-varying resistance and capacitance effects, respectively. The time-dependent capacitor is a function of the state variable, *x*(*t*)=*w*(*t*)/*D* and Δ*C*(*t*), where Δ*C*(*t*) is defined as the additional capacitance caused by changing the value of capacitance (Kim, T. H. *et al.* 2009), *C*(*x*,*t*)=((*C*_{ON}−Δ*C*(*t*))/(1−*x*(*t*))), where *C*_{ON} is the capacitance at *x*=0. The state variable is also given by
2.19

Kim, T. H. *et al.* (2009) also investigated the impact of temperature variation on their Fe_{3}O_{4} nanoparticle MR assemblies of *D* equal to 9, 12 and 15 nm. It was reported that the change in electrical resistivity (specific electrical resistance), *ρ*_{r}, as an explicit function of temperature can be defined by , which means that there is a significantly increasing behaviour as temperature decreases. As a consequence, for example, there is no hysteresis loop signature at room temperature, *T*=295 K (*D*=12 nm), while at *T*=210 K it shows a nice bow-tie trajectory. As they claimed, the first room temperature reversible switching behaviour was observed in their nanoparticle memristive system (Kim, T. H. *et al.* 2009).

It is worth noting that there are also two other elements with memory named the *memcapacitor* (memcapacitive, MC, systems) and *meminductor* (meminductive, MI, systems). Di Ventra *et al.* (2009) postulated that these two elements could also be developed someday in device form. The main difference between these three elements, the MR, MC and MI, is that the MR is not a lossless memory device and dissipates energy as heat. However, at least in theory, the MC and MI are lossless devices because they do not have resistance. Di Ventra *et al.* (2009) also investigated some examples of using different nanoparticle-based thin-film materials. MRs (memresistive systems) are identified by a hysteresis current–voltage characteristic, whereas MC and MI systems introduce Lissajous curves for charge voltage and flux current, respectively. Similar to MRs, there are two types of these elements; therefore, the three circuit elements with memory (mem systems/devices) can be summarized as follows.

*Memristors* (MR systems). An MR is a one-port element whose instantaneous electric charge and flux linkage, denoted by *q*_{mr} and *φ*_{mr}, respectively, satisfy the relation *F*_{mr}(*q*_{mr},*φ*_{mr})=0. It has been proven that these devices are passive elements (Strukov *et al.* 2008). As discussed, they cannot store energy, so *v*_{mr}(*t*)=0 whenever *i*_{mr}(*t*)=0 and there is a pinched hysteresis loop between current and voltage. Thus, the charge–flux curve is a monotonically increasing function. An MR acts as a linear resistor when frequency goes towards infinity and as a nonlinear resistor at low frequencies. Owing to the nonlinear resistance effect, *dv*_{mr}/*dt*=*R*(*t*)(*di*_{mr}/*dt*)+*i*_{mr}(*t*)(*dR*/*dt*) should be used instead of *dv*_{mr}/*dt*=*R*(*t*)(*di*_{mr}/*dt*). There are two types of control process:

—

*nth*order current-controlled MR systems,^{1}:.

—

*nth*order voltage-controlled MR systems, :.

*Memcapacitors* (MC systems). An MC is a one-port element whose instantaneous flux linkage and time integral of electric charge, denoted by *φ*_{mc} and *σ*_{mc}, respectively, satisfy the relation *F*_{mc}(*φ*_{mc},*σ*_{mc})=0. The total added/removed energy to/from a voltage-controlled MC system, , is equal to the linear summation of areas between *q*_{mc} and *v*_{mc} curves in the first and third quadrants with opposite signs. Owing to the nonlinear capacitance effect, *dq*_{mc}/*dt*=*i*_{mc}(*t*)=*C*(*t*)(*dv*_{mc}/*dt*)+*v*_{mc}(*t*)(*dC*/*dt*) should be used instead of *dq*_{mc}/*dt*=*C*(*t*)(*dv*_{mc}/*dt*), so . In principle, an MC can be a passive, an active and even a dissipative^{2} element (Di Ventra *et al.* 2009). If and capacitance is varying between two constant values, *C*_{ON} and *C*_{OFF}, then the MC is a passive element. It is also important to note that, assuming a zero-charged initial state for a passive MC, the amount of removed energy cannot exceed the amount of previously added energy (Di Ventra *et al.* 2009). An MC acts as a linear capacitor when frequency tends to infinity and as a nonlinear capacitor at low frequencies. There are two types of control process:

*— nth*order voltage-controlled MC systems, :.

—

*nth*order charge-controlled MC systems, :.

*Meminductors* (ML systems). An MI is a one-port element whose instantaneous electric charge and time integral of flux linkage, denoted by *q*_{ml} and *ϱ*_{ml}, respectively, satisfy the relation *F*_{ml}(*ϱ*_{ml},*q*_{ml})=0. In the total stored energy equation in an ML system, , the nonlinear inductive effect, *dφ*_{ml}/*dt*=*v*_{ml}(*t*)=*L*(*t*)(*di*_{ml}/*dt*)+*i*_{ml}(*t*)(*dL*/*dt*) should be taken into account. Thus, . Similar to MC systems, in principle an ML system can be a passive, an active and even a dissipative element and using the same approach, under some assumptions, they behave like passive elements (Di Ventra *et al.* 2009). There are two types of control process:

—

*nth*order current-controlled ML systems, :.

—

*nth*order flux-controlled ML systems, :.

Biolek *et al.* (2009*a*) introduced a generic SPICE model for mem devices. Their MR model was discussed in an earlier section. Unfortunately, there are no simulation results available in Biolek *et al.* (2009*a*) for MC and ML systems.

It is worth noting that there is no equivalent mechanical element for the MR. In Chen *et al.* (2009), a new mechanical suspension component, named a *J-damper*, has been studied. This new mechanical component has been introduced and tested in Formula 1 racing, delivering significant performance gains in handling and grip (Chen *et al.* 2009). In that paper, the authors attempted to show that the J-damper, which was invented and used by the McLaren team, is in fact an *inerter* (Smith 2002). The inerter is a one-port mechanical device where the equal and opposite applied force at the terminals is proportional to the relative acceleration between them (Smith 2002). Despite the fact that the *missing mechanical circuit element* has been chosen because of a ‘spy scandal’ in the 2007 Formula 1 race (Formula1 2007), it may be possible that a ‘real’ missing mechanical equivalent to an MR may someday be found, as its mechanical model is described by Oster & Auslander (1972) and Oster (1974).

## 3. SPICE macro-model of MR

Basically, there are three different ways to model the electrical characteristics of the MRs—SPICE macro-models, hardware description language (HDL) and C programming. The first approach, SPICE macro-models, is more appropriate as it is more readable for most of the readers and available in all SPICE versions. There is also another reason for choosing the SPICE modelling approach. Regardless of common convergence problems in SPICE modelling, we believe that it is a more appropriate way to describe a real device operation. Moreover, using the model as a subcircuit can highly guarantee reasonably high flexibility and scalability features.

An MR can be established by connecting an appropriate nonlinear resistor, inductor or capacitor across port 2 of an MR mutator, an ML mutator or an MC mutator, respectively^{3} (Chua 1971). These mutators are nonlinear circuit elements that may be described by a SPICE macro-model (i.e. an analogue behavioural model of SPICE). The macro-circuit model realization of a type 1 MR mutator based on the first realization of the MR (Chua 1971) is shown in figure 8.

In this model, the MR mutator consists of an integrator, a current-controlled voltage source (CCVS) ‘H’, a differentiator and a voltage-controlled current source (VCCS) ‘G’. The nonlinear resistor is also realized with resistors *R*_{1}, *R*_{2} and a switch. Therefore, the branch resistance is 1 kΩ for *V* <2 V and 2 kΩ for *V* ≥2 V. The input voltage of port 1, *V*_{1}, is integrated and connected to port 2 and the nonlinear resistor current, *I*_{2}, is sensed with the CCVS ‘H’ and differentiated and converted into current with the VCCS ‘G’.

SPICE simulations with the macro-model of the MR are shown in figures 9 and 10. In this particular simulation, a monotonically increasing and piecewise-linear *q*–*φ* function is assumed as the MR characteristic. This function is shown in figure 9*b*. The simulated MR has a value of 1 kΩ when the flux is less than 2 Wb, but it becomes 2 kΩ when the flux is equal to or higher than 2 Wb. The critical flux (*φ*_{c}) can be varied with the turn-on voltage of the switch in the macro-model. Figure 9*a* shows the pinched hysteresis characteristic of the MR. The input voltage to the MR is a ramp with a slope of ±1 V s^{−1}. When the input voltage ramps up with a slope of ±1 V s^{−1}, the memristance is 1 kΩ and the slope of the current–voltage characteristics is 1 mA V^{−1} before the flux reaches the *φ*_{c}. But when the flux becomes 2 Wb, the memristance value is changed to 2 kΩ and the slope is now 0.5 mA V^{−1}. After the input voltage reaches the maximum point, it ramps down and the slope is maintained, because the memristance is still 2 kΩ. Figure 10 shows the MR characteristics when a step input voltage is applied. Initially, the memristance is 1 kΩ, so the input current is 1 mA. When the flux reaches 2 Wb (1 V×2 s), the memristance is 2 kΩ and so the input current is now 0.5 mA as predicted. The developed macro-model can be used to understand and predict the characteristics of an MR.

Now, if a 1 kHz sinusoidal voltage source is connected across the MR model, the flux does not reach 2 Wb, so *M*=*M*_{1}=1 kΩ and *i*=10 mA. Figure 11*a*–*c*(i) shows the MR characteristics when a sinusoidal input voltage is applied. As shown in figure 11*a*–*c*(ii), when the voltage source frequency reduces to 10 Hz, the flux linkage reaches 2 Wb within 30 ms. Based on this result, there are two levels of memristance, *M*=*M*_{1}=1 kΩ and then it changes to *M*_{2}=2 kΩ.

Another interesting study is needed to verify that the model is working properly when there is a parallel, series, RM (resistor–memristor), LM (inductor–memristor) or CM (capacitor–memristor) network. First of all, let us assume that there are two MRs with the same characteristic as shown in figure 9. Analysing parallel and series configurations of these MRs is demonstrated in figure 12*a*,*b*, respectively. In both parts of figure, the left *I*–*V* curve shows a single MR.

The simulation results show that the series and parallel configurations of MRs are the same as their resistor counterparts. It means the equivalent memristances, *M*_{eq}, of two MRs in series and parallel are *M*_{eq}=2*M* and *M*_{eq}=*M*/2, respectively, where *M* is memristance of a single MR. The second step is RM, LM and CM networks. In these cases, a 10 V step input voltage has been applied to circuits. As mentioned before for a single MR based on the proposed model, while the flux linkage is less than or equal to the critical flux, *φ*≤*φ*_{c}, *M*=*M*_{1}=1 kΩ, and when the flux is more than the critical flux value, *φ*>*φ*_{c}, *M*=*M*_{2}=2 kΩ. Recall that the critical flux value based on the *q*–*φ* curve is *φ*_{c}=2 Wb. Figure 13*a*–*c*(i–iii) illustrates RM, CM and LM circuits and their response to the input step voltage, respectively.

If we assume that the memristance value switches at time *T*_{d}, then for 0≤*t*≤*T*_{d}, *φ*≤*φ*_{c}, and *M*=*M*_{1}=1 kΩ. Therefore, in the RM circuit we have *V*_{M}=*V* (*M*_{1}/(*R*+*M*_{1})). For *R*=1 kΩ, *V*_{M}=5 V and then *i*_{M}=5 mA, so *T*_{d}=*φ*_{c}/*V*_{M}=0.4 s. Likewise, when *t*>*T*_{d}, *φ*>*φ*_{c}, we have *V*_{M}=*V* (*M*_{2}/(*R*+*M*_{2}))=6.7 V and *i*_{M}=3.3 mA. Both cases have been verified by the simulation results.

In the LM circuit, we have the same situation, so for 0≤*t*≤*T*_{d}, *φ*≤*φ*_{c}, *V*_{M}=*V* =10 V, *i*_{M}=10 mA (*R*=1 kΩ) and *T*_{d}=0.2 s. For *t*>*T*_{d}, MR current is *i*_{M}=5 mA. MR’s current changing is clearly shown in figure 13*a*–*c*(iii). The CM circuit simulation also verifies a change in time constant from *M*_{1}*C* to *M*_{2}*C*.

As another circuit example of using the new MR model, an op-amp integrator has been chosen. The model of an op-amp and circuit configuration is shown in figure 14. If *C*=100 μC, then we have 0.1 and 0.2 s as the time constant of the circuit at *t*≤*T*_{d} and *t*>*T*_{d}, respectively.

It is worth mentioning that, recently, a few simple SPICE macro-models have been proposed by Benderli & Wey (2009), Biolek *et al.* (2009*b*) and Zhang *et al.* (2009) but none of them consider the model response with different circuit element types, which is an important step to verify the model correctness.

## 4. Interpreting MR in electromagnetic theory

In his original paper, Chua (1971) presented an argument based on electromagnetic field theory for the existence of the MR. His motivation was to interpret the MR in terms of the so-called *quasi-static expansion* of Maxwell’s equations. This expansion is usually used to give an explanation for the elements of circuit theory within the electromagnetic field theory.

Chua’s argument hints that an MR might exist, although it never proved that this device can in fact be realized physically. In the following, we describe how Chua argued for an MR from a consideration of quasi-static expansion of Maxwell’s equations. To consider this expansion, we use Maxwell’s equations in their differential form
4.1
4.2
4.3and
4.4
where ** E** is electric field intensity (V m

^{−1}),

**is magnetic flux density (Wb m**

*B*^{−2}),

**is electric current density (A m**

*J*^{−2}),

*ρ*is electric charge density (C m

^{−3}),

**and**

*H***are magnetic field intensity (A m**

*D*^{−1}) and electric flux density (C m

^{−2}), ∇⋅ and ∇× are divergence and curl operators.

The idea of a quasi-static expansion involves using a process of successive approximations for time-varying fields. The process allows us to study electric circuits in which time variations of electromagnetic field are slow, which is the case for electric circuits.

Consider an entire family of electromagnetic fields for which the time rate of change is variable. The family of fields can be described by a time-rate parameter *α*, which is time rate of change of charge density *ρ*. We can express Maxwell’s equations in terms of the *family time*, *τ*=*αt*, and the time derivative of ** B** can be written as
4.5

Other time derivatives can be expressed similarly. In terms of the family time, Maxwell’s equations take the form
4.6
which allow us to consider different values of the family time *τ*, corresponding to different time scales of excitation. Note that in equation (4.6) ** E**,

**,**

*H***,**

*D***and**

*J***are also functions of**

*B**α*and

*τ*, along with the position (

*x*,

*y*,

*z*). This allows us (Fano

*et al.*1960) to express, for example,

**(**

*E**x*,

*y*,

*z*,

*α*,

*τ*) as power series in

*α*, 4.7 where the zero and first orders are 4.8

Along with these, a similar series of expressions for ** B**,

**,**

*H***and**

*J***are obtained and can be inserted into equation (4.6), with the assumption that every term in these series is differentiable with respect to**

*D**x*,

*y*,

*z*and

*τ*. This assumption permits us to write, for example, 4.9 and, when all terms are collected together on one side, this makes equation (4.6) take the form of a power series in

*α*that is equated to zero. For example, the first equation in equation (4.6) becomes 4.10 which must hold for all

*α*. This can be true if the coefficients of all powers of

*α*are separately equal to zero. The same applies to the second equation in equation (4.6) and one then obtains the so-called

*nth*-order Maxwell equations, where

*n*=0,1,2,… . For instance, the

*zero-order*Maxwell equations are 4.11 and 4.12 and the

*first-order*Maxwell equations are 4.13 and 4.14

The quasi-static fields are obtained from only the first two terms of the power series equation (4.10), while ignoring all the remaining terms and by taking *α*=1. In this case, we can approximate ** E**≈

**E**

_{0}+

*E*_{1},

**≈**

*D***D**

_{0}+

*D*_{1},

**≈**

*H***H**

_{0}+

*H*_{1},

**≈**

*B***B**

_{0}+

*B*_{1}and

**≈**

*J***J**

_{0}+

*J*_{1}.

Circuit theory, along with many other electromagnetic systems, can be explained by the zero-order and first-order Maxwell equations for which one obtains *quasi-static fields* as the solutions. The three classical circuit elements *resistor*, *inductor* and *capacitor* can then be explained as electromagnetic systems whose quasi-static solutions correspond to certain combinations of the zero-order and the first-order solutions of equations (4.11)–(4.14)

However, in this quasi-static explanation of circuit elements, an interesting possibility was unfortunately dismissed (Fano *et al.* 1960) as it was thought not to have any correspondence with an imaginable situation in circuit theory. This is the case when both the first-order electric and the first-order magnetic fields are *not* negligible. Chua argued that it is precisely this possibility that provides a hint towards the existence of a fourth basic circuit device.

Chua’s argument goes as follows. Assume there exists a two-terminal device in which *D*_{1} is related to *B*_{1}, where these quantities are evaluated instantaneously. If this is the case, then this device has the following two properties:

— zero-order electric and magnetic fields are negligible when compared with the first-order fields, i.e.

*E*_{0},*D*_{0},*B*_{0}and*J*_{0}can be ignored;— the device is made from

*nonlinear*material for which the first-order fields become related.

Assume that the relationships between the first-order fields are expressed as
4.15
4.16and
4.17
where and are one-to-one continuous functions defined over space coordinates only. Combining equation (4.14), in which we have now *D*_{0}≈0, with equation (4.15) gives
4.18
As the curl operator does not involve time derivatives, and is defined over space coordinates, equation (4.18) says that the first-order fields *H*_{1} and *E*_{1} are related. This relation can be expressed by assuming a function and we can write
4.19
Now, equation (4.17) can be re-expressed by using equation (4.19) as
4.20
where the ° operator is the composition of two (or more) functions. Also, as is a one-to-one continuous function, equation (4.16) can be re-expressed as
4.21
Inserting from equation (4.21) into equation (4.20) then gives
4.22

Equation (4.22) predicts that an instantaneous relationship can be established between *D*_{1} and *B*_{1} that is realizable in an MR. This completes Chua’s argument using Maxwell’s equations for a quasi-static representation of the electromagnetic field quantities of an MR.

## 5. Conclusion

In this paper, we surveyed key aspects of the MR as a promising nano-device. We also introduced a behavioural and SPICE macro-model for the MR and reviewed Chua’s argument for the MR by performing a quasi-static expansion of Maxwell’s equations. The SPICE macro-model has been simulated in PSpice and shows agreement with the actual MR presented by Strukov *et al.* (2008). The model shows expected results in combination with a resistor, capacitor or inductor. A new op-amp-based MR is also presented and tested.

Nanoelectronics not only deals with the nanometre scale, materials and devices but also implies a revolutionary change even in computing algorithms. The von Neumann architecture is the base architecture of all current computer systems. This architecture will need revision for carrying out computation with nano-devices and materials. There are many of these different components, such as processors, memories, drivers, actuators and so on, but they are poor at mimicking the human brain. Therefore, for the next generation of computing, choosing a suitable architecture is the first step and requires deep understanding of relevant nano-device capabilities. Obviously, different capabilities might create many opportunities as well as challenges. At present, industry has pushed nanoelectronics research for the highest possible compatibility with current devices and fabrication processes. However, the MR motivates future work in nanoelectronics and nanocomputing based on its capabilities.

In this paper, we addressed some possible research gaps in the area of the MR and demonstrated that further device and circuit modelling are urgently needed. The current approach to device modelling is to introduce a physical circuit model with a number of curve fitting parameters. However, such an approach has the limitation of requiring a large number of parameters. Using a nonlinear drift model results in more accurate simulation at the cost of a much more complicated set of mathematical equations. Initially, behavioural modelling can be used; nonetheless, a greater modelling effort is needed to accommodate both the defect and process variation issues. An interesting follow-up would be the development of mapping models based on the MR to neuronmorphic systems that deal with architectural-level challenges, such as defect tolerance and integration into current integrated circuit technologies.

## Acknowledgements

We would like to thank Leon Chua, at UC Berkeley, for useful discussion and correspondence. This work was partially supported by grant no. R33-2008-000-1040-0 from the World Class University (WCU) project of MEST and KOSEF through CBNU, and is gratefully acknowledged.

## Footnotes

↵1 Current (and voltage) control is a better definition for MRs because they do not store any charge (Oster 1974). In Di Ventra

*et al.*(2009), it is specified as current (and voltage) control instead of charge (and flux) control as specified by Chua (1971).↵2 Adding energy to the system.

↵3 For further details about the mutator, the reader is referred to Chua (1968).

- Received October 21, 2009.
- Accepted February 15, 2010.

- © 2010 The Royal Society