Phase Transitions and Duality in Adiabatic Memristive Networks

0.75 MB·

by  Forrest C. Sheldon

Phase Transitions and Duality in Adiabatic Memristive Networks
Forrest C. Sheldon∗ and Massimiliano Di Ventra†
Department of Physics, University of California San Diego, La Jolla, California 92093, USA
(Dated: August 12, 2016)
Thedevelopmentofneuromorphicsystemsbasedonmemristiveelements-resistorswithmemory
- requires a fundamental understanding of their collective dynamics when organized in networks.
Here,westudyanexperimentallyinspiredmodeloftwo-dimensionaldisorderedmemristivenetworks
subjecttoaslowlyrampedvoltageandshowthattheyundergoafirst-orderphasetransitioninthe
conductivityforsufficientlyhighvaluesofmemory,asquantifiedbythememristiveON/OFFratio.
6
1 Weinvestigatetheconsequencesofthistransitionforthememristivecurrent-voltagecharacteristics
0 both through simulation and theory, and uncover a duality between forward and reverse switching
2 processes that has also been observed in several experimental systems of this sort. Our work sheds
considerable light on the statistical properties of memristive networks that are presently studied
g
both for unconventional computing and as models of neural networks.
u
A
INTRODUCTION namics and scale-invariant properties, but have not clar-
1
1 ified whether further collective behavior might arise in
higher dimensions [15].
Although systems that display resistive switching -
] To fill this gap, we study the statistical properties
h also referred to as “memristive elements” (resistors with
c memory) - have been actively studied since the 1960s, of two-dimensional disordered memristive networks sub-
e ject to a slowly ramped voltage. In this ‘adiabatic’
they have recently received renewed interest in view of
m
regime, where the applied voltage/current varies much
theirpossibleuseincomputation,bothaslogicandmem-
- more slowly than the memristance change of individual
ory components [1]. Of particular note is the tendency
t
a elements, thereisastronganalogybetweenthebehavior
of some to display a history-dependent decay constant,
st allowing a transition between a volatile and non-volatile of the network and an equilibrium thermal system. Our
. aimistounderstandthedynamicsofboththedisordered
t regime of memory [2, 3]. The resulting dynamics bear a
a devices being assembled in experiment, and the ordered
m closeresemblancetotheshort-termandlong-termpoten-
(butstronglyheterogeneous)networksbeingproposedas
tiationobservedinbiologicalsynapseswhicharethought
- novel computational architectures. To this end, we for-
d to be of central importance to learning and plasticity in
n the brain [4]. This resemblance has inspired the realiza- mulate a general model for networks in this limit which
o tion of experimental systems that seek to combine the is similar to those employed to describe metal-insulator
c transitions and electrical failure. We thus posit that the
memory features of biological synapses with the struc-
[
tural complexity of neural tissue [5, 6]. In fact, research transitions identified in these fields should also occur in

3 memristive networks and summarize work done to de-
isbeingperformedtoassesstheadvantageofusingmem-
v scribe the dynamics in these fields. Through simulations
ristive elements in non-von Neumann architectures and
2
we obtain current-voltage (I-V) curves for various values
7 is already showing that their networks dynamically or-
oftheG /G ratioanddiscussthefeaturesimpliedby
7 ganizeintothesolutionsofcomplexcomputationalprob- on off
5 lems,therebyperformingthecomputationdirectlyinthe theadiabaticlimit. TheI-Vcurvesfoundshowaduality
0 betweenforwardandreverseswitchingprocessesthathas
memory and avoiding the separation between logic and
. been observed in several experimental systems [16, 17]
1 memory units [7–10].
0 and clarifies the role of boundary conditions in the net-
All of these studies however, still lack a fundamental
6 work. Thesefeaturesarecapturedbyamean-fieldtheory
understanding of the role of time non-locality in the dy-
1 and cluster approximation which clarify the internal dy-
: namicsofmemristivenetworksandtheirstatisticalprop-
v namics and account for the features of the I-V curves.
erties. For instance, high density (∼ 109 elements/cm2)
i These results shed considerable light on the statistical
X disordered networks of memristive Ag/Ag S/Ag atomic
2 and collective properties of networks with memristive el-
r switches have been fabricated showing collective switch-
a ements that are being currently explored for neuromor-
ing behavior between a low-resistance (G ) and high-
on phic applications.
resistance (G ) state, and possible critical states po-
off
tentially useful in neuromorphic computation [4, 11, 12].
Some theoretical work has attempted to reproduce sev-
eral features observed in the experiments by performing MODEL
simulations in relatively small networks but an under-
standingofthedynamicsofsuchsystemsisfarfromclear As inspiration and as a test bed for disordered mem-
[13, 14]. Theoretical investigations of one-dimensional ristive networks we consider the Ag/Ag S/Ag gapless
2
memristive networks have shown complex temporal dy- atomic switches experimentally fabricated in [5, 6, 18].
2
FIG. 1. The network conductance for several values of G are plotted against the applied voltage for both (A) G →G
on off on
and (B) G → G . The network conductance G has been scaled to vary from 0 to 1 and the range of the applied voltage
on off
shortened to focus on the point of transition. The insets show a typical network biased by a voltage V just following the
+
transition where the formation of a (A) conducting backbone and (B) crack may be observed.
The model we consider though, is quite general and can conductivity increases, more current is diverted through
beappliedtoavarietyofothermaterialsandsystems[1]. the element from neighboring bonds making the transi-
A bipolar memristor subject to an external bias will tion from Goff to Gon unstable. On the timescale of
transition to a conductive state (with conductance G ) the slowly varying applied voltage, elements of the net-
on
in one direction and to an insulating state (with con- workwillappeartoswitchdiscretelybetweentheinsulat-
ductanceG )whenbiasedintheother. Thischangeis ing and conducting state. Therefore, we can effectively
off
generallyinducedbytherearrangementofionsintheap- model the elements as switching discretely from a con-
pliedelectricfield, asinthecaseofsilversulfidebetween ductance Goff to a conductance Gon, when a threshold
two silver electrodes, where drifting ions form a filament current, It >0, is crossed
structure eventually bridging the insulator [19].
G(I)=G +(G −G )θ(I−I ), (1)
These switching dynamics are typically subject to a off on off t
threshold in the applied bias, below which the conduc-
whereI isthecurrentthroughthedevice, andθ(·)isthe
tancewillnotvary,orwillvaryonlyslowly. Whenconsid-
Heaviside function of the argument.
eringthedynamicsofasinglememristor,suchthresholds
Similar considerations would lead us to conclude that
may be described interchangeably in the applied voltage
the reverse direction is ‘stable’, as passing an element’s
or current. However when embedded in a network, the
threshold leads to a decrease in the conductance and a
two can lead to quite different dynamics. When the con-
corresponding decrease in the current, bringing it back
ductivity of an element within a network is increased,
below the threshold. The devices would thus explore the
its current increases while its voltage decreases, and vice
full continuum of memristance between G and G
versa when its conductivity is decreased. For the elec- on off
leading to a gradual RESET-like behavior on the net-
tric field driven migration of vacancies or ions [20], we
work level. However, we find evidence both in exper-
expect that in order to simulate the behavior of actual
imental data and simulations that this effect does not
devices, all elements must be subject to a current rather
occur or is not significant in describing the behavior of
thanvoltagethreshold. Theroleoftemperaturehasbeen
thenetworkforawiderangeofparameters. Forinstance,
emphasized by several studies, especially those focussing
in the atomic switch networks produced by Stieg et al.
on the dissolution of the filament. Such an effect may
[12] sharp fluctuations in the network conductivity are
be taken into account by the explicit inclusion of a tem-
observed in both directions, thought to be due to the
perature in the model, or by taking a threshold in the
switchingofindividualelements(seeFig. 3cinRef. [12]).
power dissipated in an element P = GI2 which gives a
t t Suchsharpbehaviormaybeaccountedforbyassuminga
current threshold that changes with the device conduc-
(cid:113) nonlineardependenceoftheconductanceonthefilament
tance It = PGt. In the discretely switching model we lengthdueto,forexample,thetransitionfromtunneling
consider here, the two choices are identical, but when to ballistic conductance.
the full spectrum of memristances is allowed for, as in Disorder at the network level can similarly render the
non-adiabatic regimes, such effects may be important. continuumofmemristivevaluesirrelevanttothenetwork
Thepresenceofacurrentthresholdhasimmediatecon- dynamics. If the current diverted from a switching ele-
sequences on the dynamics of a network. As the thresh- ment does not cross another’s threshold, the increasing
old of a device in the insulating state is crossed and its current from the boundaries will continue to transition
3
that element to G . As the conductance and external fully occupied network. It is also noted in [22] that this
off
voltage range in which this occurs is very small relative is obtained upon coarse-graining a randomly pruned lat-
to the network scale, this is the same as if that element tice, and thus may be a better approximation to the
had switched discretely from G to G . Simulations thermodynamic limit than simulations on structurally
on off
of networks in which the full range of memristance was disordered lattices. The use of a disorder distribution
accessiblehavenotshownasignificantchangeindynam- also affords us more direct control over the relationship
ics and we believe the discretely switching model to be between the disorder and the dynamics, simplifying the
appropriate for a wide range of networks. We thus make searchoveralargeparameterspaceofpossiblestructural
a similar assumption for the reverse direction, obtaining disorders.
the equations for an element by exchanging G with It is worth noting that this type of model has been
on
G , and I with I in Eq. (1) for a different threshold arrived at in several contexts involving the interplay be-
off t
I <0. tween conduction and disorder in 2D systems. For in-
t
Memristive elements are also generally polar. The stance, a similar model was first applied to the study
Ag/Ag S/Ag atomic switches formed in atomic switch of the random fuse model for electrical failure (G →
2 on
networks [5, 6] are gapless type devices (see Hasegawa G = 0) [22]. A uniform distribution of thresholds on
off
et al. [21] for a review of types of atomic switches and the interval [1−w,1+w] was considered and behavior
theirswitchingprocesses). Theirsymmetricmetalliccon- examined as a function of network size L and w. Brit-
figuration (typically gapless switches have two differing tle and ductile regimes of behavior were identified, both
metallic electrodes, e.g., Ag/Ag S/Pt) suggests that at of which culminated in the formation of a lateral ‘crack’
2
the point of formation within the network, no preferred severing the network. In the brittle, or narrow disorder
direction within the switch has been selected. Polarity regime,thisoccurredasanavalancheuponthefirstbond
is instead instilled through a ‘formation’ step in which failing, while for larger w there was a regime of diffuse
a bias is applied to the network causing filament struc- failure, causing the networks to progressively deform be-
tures to form in the switches throughout. After a joule- fore the formation of a crack. In the thermodynamic
heatingassisteddissolutionofthethinnestpartofthefil- limit, only the brittle regime survived, except for the
ament, the junctions display bipolar resistive switching. case w → 1 when the disorder distribution extended to
The polarity of the internal switches is thus determined 0. More recently, in metal-insulator transitions (MIT)
by the direction of current propagation from the bound- [23] the G → G transition was examined for its
off on
aries. Thatthismustbetrueintheexperimentalsystems G /G dependence, finding that a transition occurs
on off
isevidentfromthefactthatthenetworkundergoesresis- in which the network conductivity exhibits a discontin-
tive switching as a whole. Without a majority of switch- uous jump, corresponding to the formation of a ’bolt’.
ing polarities coinciding with the direction of currents Similarly,aconductingbackbonehasbeenfoundtoform
from the boundaries, the network would switch between along the direction of current flow in the case of MIT
identical states with half the switches in the G state and dielectric breakdown [23] for sufficiently large values
off
and half in G and not display the pinched hysteresis of the G /G ratio.
on on off
observed in experiment. We thus assign the polarity of Interest in the dynamics of individual resistive switch-
elements in the network to coincide with the direction of ing (RS) devices has also led to new models such as the
currents flowing from the boundaries. Random Circuit Breaker (RCB) model [24] for unipolar
We now turn to a network of these elements. We con- devices, capable of reproducing the conductivity dynam-
sider an architecture as depicted in the insets of Fig. 1, ics of a unipolar device in the SET and RESET opera-
wheretheupperandlowerboundariesofthenetworkare tions. IntheRCBmodel, elementsofalatticetransition
held to some constant voltage or total current. The dia- toaconductivestatewhenavoltagethresholdiscrossed,
mond lattice is chosen such that all elements participate and back to an insulating state when another threshold
equally in conduction and networks are periodic in the is crossed in the same direction.
direction transverse to the current flow to mitigate finite In view of these previous results, we thus expect the
size effects. transitions observed in MIT and electrical failure stud-
Including disorder is most directly done through ran- ies to occur in memristive networks in the adiabatic
dom pruning of the lattice, however the networks pro- limit. The G → G transition will correspond to
off on
duced are subject to strong finite size effects, requiring theformationofaconductivebackboneor‘bolt’through
thesimulationofmanyinstancesoflargenetworkstoob- the network along the direction of current flow and the
tain regular results. The random pruning of the lattice G →G transition will correspond to the formation
on off
imposes a current distribution over the elements that is of a ‘crack’ transverse to the direction of current flow
simply scaled by the external boundary conditions. This severing the network. In both directions we anticipate
distributiondeterminestheorderinwhichelementscross a trivial brittle, or narrow disorder regime in which the
their thresholds and is thus equivalent, at a mean-field transitionoccursuponthefirstelementswitchingandall
level, to assuming a distribution in the thresholds of a elements transition within a narrow range of the applied
4
voltage,aswouldbethecaseifallelementshadthesame G . Each element was assigned a current threshold I
on t
threshold. Forbroaddisorderweexpectaductileregime from Uniform(0,1) in each direction, beyond which they
in which there will be diffuse switching leading up to the transition from G → G , or vice versa. Network di-
off on
transitionandactivityoverabroadrangeofappliedvolt- mensions were chosen so that the total network conduc-
ages. While the ductile regime does not survive in the tance varied between G and G (N = N = 128).
off on x y
thermodynamic limit for electrical failure, the modest Theinitialvoltageissettothevaluerequiredtocrossthe
size of memristive networks experimentally realized [1] lowest threshold in the network. Once that element has
suggests the ductile regime is still significant in their dy- switched,voltagesandcurrentsarerecalculatedthrough-
namics. Of particular interest to us is the influence of outthenetworkwiththeexternalvoltageheldfixed,and
such transitions on the I-V curves of the network. From allotherelementswhosecurrentsexceedtheirthresholds
the perspective of the experimenter such effects may be are switched. This is repeated until no currents exceed
desirable, such as providing strong sensitivity across a the thresholds of their elements, at which point the volt-
small voltage range or signalingthe solution of a compu- age is raised until another threshold is crossed and this
tationalproblem,orundesirable,byreducingthenumber process repeated. The forward and reverse protocols are
of internal states accessible to the network. identical aside from the initial state and switching direc-
Investigations of the behavior of memristive networks tion of the elements.
have been limited to date. In [13], the full time in- InFig.1weshowthenetworkconductancesasafunc-
tegration of a memristive network was undertaken for tion of applied voltage for various values of G (setting
on
networks of moderate size. Elements did not include a G = 1), in both forward G → G and backward
off off on
thresholdintheirdynamicsandthenetworkswereinves- G → G transitions, and for threshold distribution
on off
tigatedfortheirdependenceonthefractionofmemristive p(t)=Uniform(0,1). The displayed networks have a lin-
to ohmic conductors p, and their AC response. It was ear size of N =128 memristors which we found large
x/y
found that for poor ohmic conductors (G = G ) the enough to achieve regular results over multiple realiza-
off
networks exhibited pinched hysteresis curves only when tions of the disorder. Network conductances have been
p>0.5, the percolation threshold. For good ohmic con- scaled to lie on the interval [0,1]. We note that for small
ductors (G=G ), a strongly memristive phase was ob- values of G in both directions, the conductance is a
on on
servedforp>p inwhichthenetworksswitchedabruptly smooth function of the voltage. As G is increased, a
c on
and for p<p a weakly memristive phase was observed, discontinuity forms in the slope which sharpens, appear-
c
similar to that for poor ohmic conductors and p>p . ingalmostcontinuousuntiladiscontinuousjumpappears
c
Modelingperformedby[14]ofthe“atomicswitchnet- for large G analogous to a first-order phase transition.
on
works”simulatedsmallnetworksincludingvolatilityand Similar behavior was seen for a variety of other distri-
noise in their memristor model, and showed an opening butions of sufficient breadth (not shown) with the point
oftheI-Vcurvesasthenoisetermwasreducedand1/fγ of transition, however, being distribution dependent. In
(0≤γ ≤2) scaling of the power spectral density for the the insulating transition, we have scaled the voltage by
small networks simulated (L ≈ 10). While such studies G such that the current densities of all networks are
on
reproducephenomenaseeninexperiment,littleworkhas initially equal, bringing the transitions to the same scale
been done to analyze the behavior of these models and in both polarities. While the dependence on G /G
on off
understand how memristors in networks interact. Here, in the forward direction has been shown in MIT [23], we
we instead focus on building an understanding of mem- are not aware of a similar demonstration in the reverse
ristive networks in the adiabatic regime, where analogies direction,possiblyasmostworkhasfocusedonelectrical
with thermodynamic systems are strongest. By moder- breakdown in fuse networks in which G =0.
off
ating the strength of interactions through the G /G The insets in Fig. 1 each show the configuration of
on off
ratio, we examine the transition that occurs in each di- a network just following a transition. In the forward
rection and it’s effects on the I-V curves of the network. transition (A), this corresponds to the emergence of a
The features of this transition and the resulting hystere- conducting backbone spanning the network. Similarly,
sis curves are well captured by a cluster approximation inthereversetransition(B)acrackformsseparatingthe
thatwellapproximatesthebehaviorofthenetworkabout network transverse to the direction of current flow.
the point of transition. The corresponding I-V curves are shown in Fig. 2.
While the voltage scale depends only on the geometry of
thenetwork,thecurrentscaledependsonG /G and
on off
SIMULATIONS has been rescaled so that the axes coincide. The strong
asymmetry of the curves is due to the use of the same
Simulations were carried out for a square lattice at a threshold distribution for G →G and G →G
on off off on
variety of sizes, G /G ratios, and threshold distri- processes. In the G state, an equivalent current will
on off on
butions p(t). The network dimensions were chosen such be reached for a voltage which is a factor 1/G lower.
on
thattheconductivityofthenetworkvariedfromG to As there is not an obvious physical choice for how to
off
5
FIG. 2. Simulated hysteresis curves of memristive networks for G = 4, 100, and 1000. The discrete jump in conductance
on
becomesevidentintheI-VcurvesforlargevaluesofG Theasymmetryofthecurvesistheresultofelementthresholdsbeing
on
in terms of the current, and thus the transition occurring at a factor of 1/G lower voltage than the corresponding forward
on
transition. TheinsetsshowanalyticallyreachedI-Vcurveswhichdemonstratesimilarasymmetryandtheemergenceofajump
inthecurrent. Thereversebranchhasbeenrescaledandplottedinthesecondrow. ForlargeG ,thetransitionappearsasa
on
noisyareaneartothey-axis. Herethejumpsinconductivityatafixedvoltagegiverisetosharpdecreasesinthecurrentthat
are opposed by the subsequent increase in voltage.
connect the distributions for forward and reverse switch- gion about the transition however, the jump of the for-
ing, we have displayed the negative voltage sections of ward branch appears as a fluctuating region in the re-
the curves separately, on their own current scales. The verse branch. Here, avalanches of transitioning elements
discrete jump in the conductance becomes evident for sharply reduce the current, which is then opposed by a
sufficiently large values of G , but is less apparent than subsequent increase of the externally applied voltage. If
on
in the conductivity plots due to the long tail following weinsteadrunthereversebranchwithcurrent-controlled
thetransition. Pastthetransition,voltagestepsbetween boundary conditions, this fluctuating region becomes a
thresholds become longer as current is diverted into the discrete jump as seen in the first row of Figure 3. In
conducting backbone. The transition thus has an in- the second row, upon the mapping (2), the curves align
hibitory effect on the remaining bonds, opening the hys- nearlyexactly. Thiscorrespondencebetweenthetwopro-
teresiscurvesandspreadingthememristivestatesovera cesses is just the familiar I-V duality of electrical cir-
wider voltage range. Thus, while the transition reduces cuit theory [25]: the diamond lattice is dual to itself
thenumberofaccessiblestates,itincreasestheresolution and taking G→ 1 in all links takes a voltage-controlled
G
between those remaining. insulator-to-metaltransitiontoacurrent-controlledelec-
InthesecondrowofFig. 2,thescaledreversebranches trical failure process.
oftheI-Vcurvesappearsremarkablysimilartotheposi-
In the context of memristor networks, this indicates
tivebranch,butwiththerolesofthecurrentandvoltage
that the voltage-controlled I-V curves are dual to the
exchanged. This suggests the following mapping of the
current-controlled I-V curves upon exchanging the roles
reverse branch to the forward branch,
of voltage and current and the direction of the switching
IN VG N process. Runningthemodelinthecurrentcontrolledset-
V → y , I → off x, (2)
tingthusgivesnearlyidenticalresults,butexchangesthe
N G N
x on y
fluctuatingregioninthereversedirectionforthediscrete
where N and N are the lattice dimensions. In the re- jump in the forward direction.
x y
6
FIG. 3. For a current-controlled network switching from G → G the transition gives a discrete jump, as in the voltage-
on off
controlled G →G case. After the mapping in Eqn. 2, the reverse branch of the hysteresis curves (red) have been plotted
off on
over the forward branch for several values of G demonstrating the duality between the two processes.
on
It is important to note that this connection between ical quantity considered is the power dissipated by the
the forward and reversed hysteresis loops has been ob- network,
served experimentally in individual memristive systems
[16, 17] as well. While observing the duality for a net- P =(cid:88)g v2 =(cid:88) i2j (3)
work of resistors is clear, the extension to the dynamic j j gj
j j
non-linear elements considered here is surprising, espe-
whereg istheconductanceofanelementinthenetwork
cially when considering the asymmetry in the switching j
and v (i ) is its voltage drop (current). We require that
processbetweentheforwardandreverseddirections,one j j
theaveragepowerdissipatedmatchthepowerdissipated
of which corresponds to a backbone forming along the
direction of current propagation, and the other corre- by the network GnetV2 = GIn2et and assume that all el-
sponding to a crack forming transverse to the current ements experience a mean-field voltage VMF or current
flow. IMF leading to the equations,
In the remainder of the paper, we analytically investi- G V2 (cid:40)
gate our model with the aim of understanding the major net  N(cid:104)g(cid:105)V2
I2 = MF (4)
featuresofthesimulatedI-Vcurves,namelytheexistence N(cid:104)1(cid:105)I2 .
of a transition for sufficiently large Gon/Goff, the long Gnet g MF
tailfollowingthetransition, andthedualitybetweenthe
The choice of the LHS is determined by the bound-
forward and reverse switching processes, leading to the
ary conditions applied to the network but the choice of
I-V curves displayed in the insets of Figs. 2 and 3.
the RHS is not constrained. An interesting form is the
‘voltage-voltage’choice,leadingtothemean-fieldvoltage
(cid:115)
MEAN-FIELD THEORY
G V
V = net√ (5)
MF (cid:104)g(cid:105) N
As a first step towards an analysis of the model, we
develop its mean-field theory. The method followed is which unlike other choices displays a transition in both
similartothatofZapperiet al. [26]employedtoanalyze directions. Tomakeprogresswerequirethenetworkcon-
random fuse networks. In this form, the central phys- ductance G . Below the transition, where switching of
net
7
FIG.5. Thelhs(blue)andrhs(red)ofEq.(6)areplottedfor
several values of the applied voltage. For low values of G
on
their intersection gives a solution that is smooth function of
the voltage (panel A). A transition develops for intermediate
(cid:113)
FIG. 4. h(f) = G(f) is plotted for several values of values that can appear continuous (panel B). For sufficiently
(cid:104)g(cid:105)(f) large values of G , a transition occurs where the solution
G . In the inset, (cid:104)g(cid:105)(f) and G(f) are plotted for G = on
on on jumps discontinuously. The transition voltage is highlighted
100. Notethatwhentheaverageconductanceincreasesmore
in panels B) and C).
quickly than the network conductance as for small f, h(f) is
decreasing and vice versa for large f.
the conductance ratio G /G and is independent of
on off
elementsisprimarilydrivenbythethresholddistribution the length scale of the disorder, both amounting only to
andnotbytheinfluenceofnearbyswitchedelements,the a rescaling of the applied local field v. The lhs and rhs
conductivityofthenetworkmaybewellapproximatedby of Eq. (6) are plotted for several values of the voltage
an effective medium theory, giving Gnet ≈ Geff(f) as a in Figure 5. For the chosen distribution and Gon/Goff
function only of the fraction of the devices in the ON ratios a transition is evident at the point
state.
The functions (cid:104)g(cid:105)(f), G (f), and h(f) = (cid:113)Geff(f) 1=p(h(f)v)h(cid:48)(f)v 0≤f ≤1. (7)
eff (cid:104)g(cid:105)(f)
are plotted in Figure 4. We can understand the non-
We also note that the inflection of the curves from
monotonic form of h(f) as arising from competition
the RHS of Eq. (6) shows a trend that looks almost like
between switching elements concentrating current away
a continuous transition, corresponding to the behavior
from other elements, and the increasing conductance of
seen in simulations for intermediate values of G (see
on
the network pulling more current in at the boundaries.
Figure 1).
Forsmallf,theaverageconductanceofanelementisin-
An exactly analogous treatment may be undertaken
creasing faster than the network conductivity, indicating
for the transition from G → G . We regard f =
on off R
that current is concentrated away from other elements
1−f as the fraction of devices in their G state and
off
more quickly than it increases at the boundaries, and v = √V as the positive voltage applied in the reverse
the mean-field voltage decreases. For larger f, the net- N
direction. Theeffectivemediumconductivitymaybeob-
work conductance begins to increase more quickly than
tained from the substitution f → 1−f , G (f ) =
the average conductance, pulling in current faster than R eff,R R
G (1−f ) and similarly with the average conductiv-
switchingelementscanconcentrateit,andthemean-field eff R
ity (cid:104)g(cid:105) (f ) = (cid:104)g(cid:105)(1−f ), giving a mean field voltage
voltageincreases. Theincreasingregimeatlargef allows R R (cid:113) R
for a phase transition to occur from Goff → Gon, and VMF = hR(fR)v = G(cid:104)egf(cid:105)fR,R(f(RfR))v. As the mechanisms
the decreasing portion at low f allows for the possibility for turning ON and OFF within the individual atomic
of the reverse transition when the voltage is reversed. switches are not the same, we take a possibly different
To determine whether a transition occurs for a partic- probability distribution p (t) for the reverse switching
R
ular disorder distribution, we derive a self-consistency thresholds. With these definitions we obtain the self-
equation, ensemble-averaging over the number of ele- consistency equation as before,
ments that have switched for a given mean-field voltage.
ForthetransitionfromG toG ,thefractionthathas
off on (cid:90) hR(fR)v
switched will approach the average fraction of elements
f = p (t)dt. (8)
R R
with thresholds below the mean-field voltage, 0
(cid:90) h(f)v In Fig. 5 this quantity has been plotted for several val-
f = p(t)dt. (6) ues of G and the applied voltage. We observe a first-
on
0 order phase transition similar to that observed in the
Because the applied field enters multiplicatively, the dy- G →G branch. However, this transition occurs for
off on
namics given by the mean-field theory depend only on much lower values of G and near the limiting value of
on
8
1D MODELS
The mean-field assumption, that all elements experi-
ence either a voltage V or current I , is equivalent
MF MF
to replacing the network with a 1D parallel or series ar-
rangementofmemristorswhoseboundaryconditionsare
then matched (through an effective medium theory) to
the behavior of the original network. In fact, the quanti-
ties
FIG.6. Solvingtheselfconsistencyequations(6)and(8)fora n
uniformdistributionleadstothehysteresiscurvesabove. For (cid:104)g(cid:105)=Goff + N(Gon−Goff) (9)
small values of G /G , the networks are smooth, but as
on off
the ratio is increased a discrete jump emerges in the forward
direction. A similar jump near the end of the reverse branch (cid:28)1(cid:29) G + n(G −G )
is only barely discernible. = on N off on (10)
g G G
on off
whichappearinthemean-fieldequations(4)arealsothe
conductancesG /N foraseriesandparallelnetworkof
net
N memristors with n in the ON state. Having seen that
the conductance. As we proceed from f = 1 to 0, in
the transition is restricted to a small subset of the net-
the vicinity of f = 0 the average conductance (cid:104)g(cid:105)(f) is
work,wedonotexpectthatincludingtheentirenetwork
decreasing more rapidly than the network conductance
in the backbone will capture the behavior in the vicin-
G (f). This leads the internal switches to redistribute
eff
ity of the transition (where homogeneity, and thus the
currenttotheirneighborsfasterthanthedecreaseofthe
effective medium theory, fails). We thus first explore the
total current at the boundary. This increases the mean-
opposite extreme by ignoring the presence of the rest of
field voltage overall, and hence promotes a transition.
thenetworkandconsideringonlythoseelementsinvolved
Solving the self-consistency equations gives the mean- in the conducting backbone or crack.
field hysteresis curves plotted in Fig. 6. While these Intheforwarddirection,thetransitioningelementsare
curves display a qualitative similarity to simulation, sev- a collection of memristors in series of length Ny held at
eral features of the mean-field theory are lacking. As a voltage V. Such an arrangement with a fraction f in
we have already noted, the choice for the form of the the ON state admits a current,
mean-field theory (voltage-voltage, current-voltage, etc.
G G V
) are not prescribed a priori and each choice will give a I(f)= on off (11)
G +f(G −G )N
slightly different account of the dynamics. All of these on off on
share the feature that transitions occur due to competi-
and the fraction of elements in the ON state may be
tion between a changing current at the boundaries and
determined self-consistently,
the internal sharing of currents within the network, as
summarized by the function h(f) which is some ratio (cid:90) I(f)
f = ρ(t)dt. (12)
between the network conductance and an average con-
ductance (cid:104)g(cid:105), (cid:104)1(cid:105). In both directions, the transition will 0
g
eventuallyproceed(asG →∞)fromsomecriticalfrac- The distribution ρ(t) is the distribution of thresholds
on
tion f to a fully switched network f = 1 as opposed to in the conducting backbone which should be related to
c
the finite jump and long tail seen in simulations. the distribution of thresholds across the network. While
an explicit calculation of ρ is difficult, a reasonable ap-
Having observed the internal form of the transition in proximation on the diamond lattice should be ρ(t) =
simulation, we can see that in contrast to thermal tran- 2p(t)(1 − F(t)) where F(t) is the cumulative distribu-
sitions, where the phase transition occurs homogenously tion function of the threshold distribution, such that the
throughoutthesystem, theconductivitytransitionscon- currentalwaysselectsthepathwithlowerthreshold. We
sist only of a d=1 dimensional conducting backbone in notethatthisconcentratesthethresholddistributionto-
the forward direction and a d=D−1 dimensional crack wards its lowest values but does not strongly alter the
in the reverse. In D =2 both of these correspond to one behavior of the theory. In the interest of simplicity, we
dimensionalsubsystemsofthenetworkandsothemean- thusmaintainouruseoftheUniform(0,1)distributionin
field theory which considers all elements equally cannot illustrating the features of each approach. Equation (12)
modelitaccurately,especiallyintheregimefollowingthe is plotted in Fig. 7 for several values of G in which we
on
formation of the backbone. In the following, we consider observe the transition first occurring at G = 2, and
on
methods for modeling the formation of the backbone. then progressing to a jump to f =1 for larger values.
9
FIG. 7. The lhs (blue) and rhs (red) of Eq. (12) are plotted
for several values of the applied voltage for the distribution
Uniform(0,1). For a collection of memristors in series, the
transitionisthecompletionoftheconductingbackbone,with
all elements in the G state.
on
In the reverse direction, we consider a collection of FIG. 8. The systems and sub-units considered in the clus-
memristors in parallel of length Nx corresponding to the ter approximations. As a model of the conducting backbone
crack that will eventually sever the network. As the embedded in the network we consider (A) a series ‘chain’ of
‘crack’ is separated from the boundaries, the boundary memristorsinparallelsubjecttoaslowlyrampedvoltageV+
composed of sub-units of (B) pairs of memristors in parallel
conditions are instead supplied by the network. As ev-
subjecttoaslowlyrampedcurrent. Asamodelofthecrack,
ery strip of memristors perpendicular to the direction of
weconsider(C)aparallelstripofmemristorsinseriessubject
current propagation will have a current G V passing
net toaslowlyrampedcurrentconsistingofsub-unitsof(D)pairs
throughthem,thereverseswitchingprocessofavoltage- of memristors in series subject to a slowly ramped voltage.
controllednetworkshouldbebestdescribedbyacurrent-
controlled strip of memristors in parallel. For the mo-
ment,weagainignorethepresenceoftherestofthenet- CLUSTER MODELS
work and consider an isolated set of elements. The con-
ductance of the 1D strip of memristors in parallel with a In order to include the influence of the conducting
fraction fR in the Goff state is backbone or crack on the rest of the network and thus
the behavior of the network past the transition, we use
N(G +f (G −G )) (13) a cluster approach similar to the Bethe-Kikuchi approx-
on R off on
imation in equilibrium thermodynamics [27, 28]. To this
which gives the current through an element in the ON end we replace each memristor in the series chain with
state, two memristors in parallel as in Figure 8 A and subject
theentirechaintoaslowlyrampedvoltage. Eachpairin
G G V thechainisthussubjecttoacurrentI slowlyraisedfrom
I(f )= on net . (14)
R G +f (G −G )N zero (See Fig. 8 B). If the threshold of each memristor,
on R off on
t is drawn independently from a distribution p(t), the
i
f may be similarly found self-consistently probability distribution for the conductance of the pair
R
G is,
(cid:107)
(cid:90) I(fR) (cid:90) ∞ (cid:90) ∞
f = ρ(t)dt. (15)
R p(2G ;I)=2 dt dt p(t )p(t )
0 off 1 2 1 2
I/2 t1
Theresultingmean-fieldtheoryisjustthereverseofthat (cid:90) I/2 (cid:90) ∞
p(G +G ;I)=2 dt dt p(t )p(t )
fbourttwhiethfoarwvaorldtagsweistccahliengsmparlolecresbsy(atafkaicntgorGonfeGt =.Gon) off on GoGffof+fGIon 1 t1 2 1 2
on
Here, physical considerations from the switching pro- +2(cid:90) GoGffof+fGIon dt (cid:90) ∞ dt p(t )p(t )
cesses have led us to two dual structures: a series chain 0 1 GoffI 2 1 2
of memristors transitioning from G →G subject to Goff+Gon
off on
a ramped voltage as a model of the conducting back- p(2G ;I)=2(cid:90) GoGffof+fGIon dt (cid:90) t2dt p(t )p(t ).
bone, and a parallel strip of memristors transitioning on 2 1 1 2
0 0
from G → G subject to a ramped current for the
on off (16)
crack severing the network. Each of these demonstrates
a transition in which the networks proceed from some A long chain of such pairs in series (Fig. 8 A), each with
f =f to f =1 at a critical voltage or current. conductance G will possess a total conductance close
c (cid:107),i
10
to The two structures considered above are again dual
and while the self-consistency equation may be solved as
(cid:28) (cid:29) (cid:28) (cid:29)
1 1
=N (17) before, we instead note that the exchange
G G
chain (cid:107)
IN VG N
V → y , I → off x (21)
andwemaydeterminethecurrentthroughthechainself- N G N
x on y
consistently as the smallest solution of the equation
takes the above self-consistency equation, to that of the
V/N (cid:28) 1 (cid:29) forward switching process (we consider a square network
x = (18)
I G(cid:107) Nx =Ny =N to avoid dimensional factors)
(cid:28) (cid:29)
V/N 1
where the I dependence of (cid:104) 1 (cid:105) has entered through the N(cid:104)G (cid:105)V =I → = (22)
averaging. Using again theG(cid:107)distribution Uniform(0,1), series I G(cid:107)
thesolutiontothisequationhasbeenplottedintheinsets The reversed switching process thus maps exactly to the
of Figs. 2 and 3 for several values of Gon. forward switching process upon exchanging the roles of
Here, we see the finite jump in conductance observed thevoltageandcurrentandscalingappropriately,asseen
in simulations followed by the gradual switching of the in the simulated I-V curves of Figure 3.
remaining memristors. This is due to the memristor
in the ON state diverting current away from its neigh-
bor. While a current of only I = 2 (in units of G V) MEAN-FIELD DYNAMICS
off
is required to switch the first of the pair, a current
of I ≈ (G + G )/G is required to guarantee Although the avalanche dynamics of mean-field mod-
on off off
the switching of the second, which for large values of els akin to the random-field Ising model are well known
G /G is considerable. [29], we include here a brief discussion in the interest of
on off
Ananalogoustreatmentofthereversedswitchingpro- completeness.
cess requires replacing each element of the parallel strip The mean-field theories considered lead to self-
of memristors with two elements in series (Fig. 8 C) consistency equations of the form
subject to a slowly ramped current. Each pair is then
(cid:90) h(f)v
subject to a slowly ramped voltage V (Fig. 8 D). Again
f = p(t)dt (23)
drawingthethresholdsindependentlyfromadistribution
0
p(t), the distribution for the conductance of the pair is
where v is an external field and h(f) is some function
G (cid:90) ∞ (cid:90) ∞ of the fraction of switched elements f. We consider this
p( on;V)=2 dt dt p(t )p(t )
2 1 2 1 2 relationtobeinitiallysatisfiedandraisethevoltageuntil
GonV/2 t1 the next threshold is passed. This takes f → f + 1,
p( GonGoff ;V)=2(cid:90) GonV/2 dt (cid:90) ∞dt p(t )p(t ) increasing the limit of (23) by h(cid:48)(f)v. The probabiliNty
Gon+Goff GGoonffG+ofGfoVn 1 t1 2 1 2 that n memristors are switched ONN by this increase is
(cid:90) GonGoffV (cid:90) ∞ given by a Poisson distribution
+2 Goff+Gon dt1 dt2p(t1)p(t2) µn
0 GGoonffG+ofGfoVn pn = n!e−µ, µ=p(h(f)v)h(cid:48)(f)v. (24)
p(Goff;V)=2(cid:90) GGoonffG+ofGfoVn dt (cid:90) t2dt p(t )p(t ). Each of the n memristors will cause a similar increase in
2 2 1 1 2 themean-fieldvoltage,andthuswillgiverisetothesame
0 0
(19) distribution. Therefore, by switching a single memristor
gives rise to a poissonian branching process. This may
AstripofNy ofthese,eachwithconductanceGseries,i, bebroughtintoamoreusefulformbycalculatingtheto-
in parallel (Fig. 8 C) will have conductance Ny(cid:104)Gseries(cid:105). tal number of memristors switched in a single branching
As discussed above, such a strip embedded within a net- process, or the avalanche size distribution. This leads to
work will be subject to a current I and thus satisfy the Borel distribution
Ny(cid:104)Gseries(cid:105)V =I, (20) p = (µS)S−1e−µS, S =1,2,... (25)
S S!
where the voltage across the strip is determined self-
with µ = p(h(f)v)h(cid:48)(f)v. This distribution has mean
consistently as the smallest solution of the above equa-
and variance,
tion and the voltage dependence of (cid:104)G (cid:105) has entered
series
through the averaging over the voltage dependent distri- 1 µ
(cid:104)S(cid:105)= , σ2 = . (26)
bution above. 1−µ S (1−µ)3