- Email: [email protected]
Monte Carlo simulation of current–voltage characteristics in metal–insulator–metal thin film structures
Thin Solid Films 433 (2003) 321–325
Monte Carlo simulation of current–voltage characteristics in metal– insulator–metal thin film structures S. Gravanoa,*, E. Amrb, R.D. Goulda, M. Abu Samrab a Thin Films Laboratory, Department of Physics, School of Chemistry and Physics, Keele University, Keele, Staffs, ST5 5BG, UK Department of Physics, College of Science and Technology, Al Quds University, Abu Deis, P.O. Box 20002, Jerusalem, West Bank
b
Abstract Many metal–insulator–metal (MIM) thin film sandwich structures exhibit an electroforming process, following which the resulting direct current (DC) current–voltage I–V characteristic exhibits ohmic conduction and voltage-controlled differential negative resistance (VCNR) at low and high voltages, respectively. One of the more successful models of these phenomena is the filamentary model, in which it is assumed that the electroforming process establishes a population of ohmic filaments within the insulating matrix, which span the metal contacts. The VCNR behaviour results from the progressive cessation of conduction in individual filaments owing to Joule heating effects. Early work showed that by postulating plausible probability distributions of filament resistances, I–V characteristics typical of those obtained experimentally could be derived. More recently more realistic normal distributions of filament resistances, radii and cross-sectional areas have been considered, and these too have yielded characteristics in accordance with experiment. In the present work normal distributions of filament resistances and radii are considered, but in contrast to previous work a novel method using a Monte Carlo simulation was employed. Furthermore, the effects of using a combination of two different filament resistivity values were explored. Variations in conductance and shifts in the current peak were explicable using existing theory. The approach adopted in this work lends itself to the exploration of more complex filament distributions, and therefore to a better representation of the underlying filament population, including those having a multivariate distribution. 䊚 2003 Elsevier Science B.V. All rights reserved. Keywords: Monte Carlo simulation; Electrical properties; Metal-insulator-metal structures (MIM)
1. Introduction The current–voltage (I–V) characteristics shown by metal–insulator–metal (MIM) thin film sandwich structures yield useful information concerning the dominant conduction process through the insulator. Under certain circumstances a phenomenon known as electroforming occurs. Essentially this consists of applying a direct current (DC) voltage, in excess of a few volts, between the metal electrodes. Providing this voltage exceeds a level known as the forming voltage VF, electroforming occurs and the sample shows various characteristic features. These include a considerably increased conductance, voltage-controlled differential negative resistance (VCNR) behaviour, electron emission and various memory effects. Early models of the conduction process in electroformed structures were based on the establishment *Corresponding author.
of a modified band structure during electroforming w1,2x. Following these a novel filamentary model was proposed by Dearnaley et al. w3x in which the conductivity is dominated by a matrix consisting of a large number of parallel conducting paths, or filaments. At low voltage the conduction is essentially ohmic, but as the voltage is increased Joule heating within the filaments causes them to cease to conduct or rupture, and the current level falls in the negative resistance region of the VCNR characteristic. Fig. 1 w4x shows a typical current–voltage characteristic for a Au–CaBr2 –Au sandwich structure. Ic represents the current circulating through the sample and Ie the current emitted as functions of the applied voltage Ub. We are not concerned with emission current in this simulation work. This characteristic is representative of those for many different insulator and electrode combinations. In common with other characteristics presented in the literature, the circulating current initially increases with voltage faster than expected from Ohm’s
0040-6090/03/$ - see front matter 䊚 2003 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 0 - 6 0 9 0 Ž 0 3 . 0 0 3 8 7 - 0
322
S. Gravano et al. / Thin Solid Films 433 (2003) 321–325
tribution of resistances is then given by
PŽr.s
Fig. 1. A typical experimental current–voltage characteristic for an electroformed Au–CaBr2 –Au sample. Ic represents the current circulating through the sample, Ie the current emitted and Ub the applied voltage. Reproduced from w4x. Reprinted with permission from R.D. Gould and C.A. Hogarth, Phys. Stat. Sol. (A) 23, 531. Copyright Wiley-VCH, 1974.
law, following typically either an IAI0sinh(VyV0 ) or an IAI0(VyV0)2 dependence, where I is the current through the sample, V is the applied voltage, and I0 and V0 are constants w3x. The conducting filaments were initially characterised in terms of a probability density distribution of filament resistances P(r) which was triangular, where r represents the filament resistance w3x. Later work on the form of the function P(r) suggested that a parabolic distribution equally well described the experimental results w5x, while it was also argued that a normal distribution of filament resistances was more fundamental w6x, and that the earlier suggested distributions were merely approximations to this. More recently the present authors have suggested that there is no a priori reason for suggesting that the filament resistances are normally distributed, and that more basic attributes such as the filament radius w7x or cross-sectional area w8x may be normally distributed. Indeed, both of these assumptions have been shown to yield computed characteristics in accord with experiment. Ray and Hogarth w6,9x have also suggested that more than one variable may be normally distributed, for instance the resistivity and the filament radius. The degree of mathematical difficulty for analytical solutions, even for a single distribution of filament crosssectional areas, becomes increasingly complex, and numerical methods are normally required w8x. For multivariate distributions a more easily applicable method is clearly required, and the present work is an attempt to develop such a method. Before presenting the general results obtained in such initial simulations, it is necessary to outline the filamentary model for normal distributions of filament resistances and radii, and to briefly describe the simulation. 2. The filamentary conduction model and simulation We consider first the situation when the filament resistances r are normally distributed. The normal dis-
w
1
y2psr2
expxy y
Žryrm.2 z
|
2s2r
~
(1)
where rm is the mean and sr is the standard deviation. A filament ruptures when the power dissipated in the filament reaches a particular value b2, which is determined by the heat flow symmetry from the filament, and depends primarily on the filament dimensions and density, the filament and bulk temperatures and the insulator thermal conductivity w10x. b2 is effectively a factor indicating how the power is dissipated, and it follows from this that the voltage Vp at which a filament of resistance r will rupture is given by Vpsbr1y2
(2)
The conductance G(s1yR), where R is the resistance, of an assembly of N filaments when a voltage V is applied is given by ru
1 sN R
|
PŽr. rv
ry
(3)
dr
where ru represents the maximum filament resistance and rv represents the value of the lowest resistance in the distribution which is still conducting at voltage V, and may be calculated from Eq. (2). The lowest resistance in the distribution is denoted by rt. When VFbr1y2 t , all the filaments in the distribution are conducting, and the lower limit of the integral rv is given 1y2 by rt. For br1y2 t FVFbru , where some filaments are conducting the lower limit rvsV 2 y b2, and for VGbr1y2 where no filaments conduct rvsru. u If we then consider the filament radii r to be normally distributed, the normal distribution of radii is given by
PŽr.s
1
y2psr2
w
expxy y
Žryrm.2 z 2s2r
| ~
(4)
where rm is the mean and sr is the standard deviation. Although expressions analogous to Eq. (2) and Eq. (3) also exist for the distribution of filament radii w7x, it was not necessary to use these in the simulation. All that was necessary was to calculate the resistance of each filament, using rs
r0d pr2
(5)
where r0 is the filament resistivity and d is the filament length.
S. Gravano et al. / Thin Solid Films 433 (2003) 321–325
In the simulation the current I was calculated for different values of voltage V. For a specified voltage V filaments with resistance r)V 2 y b2 exhibit ohmic conduction and contribute to the overall current, whereas those with rFV 2 y b2 do not conduct. Since we have made the approximation of ohmic conduction through individual filaments, it is clear that the initial I–V dependence will be linear where all filaments are conducting. The simulation was based on a Monte Carlo approach in which for each filament a value of either resistance or radius was randomly generated using the ¨ Box–Muller method w11x. The resistance or radius value for the filament was determined by specifying the mean and standard deviation of P(r) or P(r), as appropriate. For simulations involving the radius, the resistance was calculated from Eq. (5). If r)V 2 y b2 the filament contributes a current Vy r to the overall current I. If rFV 2 y b2 the filament does not conduct, and makes no contribution to the overall current. The total current flowing with voltage V applied was obtained by summing the currents through the individual filaments. This was then repeated over a range of voltages to generate the complete I–V characteristic. 3. Results and discussion The results of a selection of the simulations are shown in the following figures. Two fixed filament resistivity values r0 were used. The lower of these is that adopted by Dearnaley et al. w3x of 0.002 V m, which assumes that the resistivity is determined largely by that of metallic gold atoms derived from the electrodes, which are injected into the insulating matrix in the filamentary regions. A higher resistivity value of 0.0034 V m was also used in some simulations. This value was chosen to reflect an alternative view of filamentary conduction, in which carbon is considered to be an important component of the filamentary material w12x. Note that these resistivity values do not correspond directly to those of gold and carbon, but are chosen (i) to yield realistic current values for the assumed number of filaments N, and (ii) to reflect different filament compositions, depending on the different origins of a major filament component. Following previous work on filamentary conduction, other parameters of the model were as proposed by Dearnaley et al. w3x. In particular these were the number of filaments Ns5=106, filament length ds4=10y8 m and heat flow constant bs 2.1=10y3 W1y2. These values were also used in earlier work concerned with normally distributed filament radii w7x, and therefore allow direct comparison with this work. In order to demonstrate the feasibility of the method, initial work was performed using a normal distribution of filament resistances, while the bulk of the work used a more realistic normal distribution of filament radii.
323
Fig. 2. Current–voltage characteristics derived from normal distributions of filament resistances with mean resistance 8=107 V, and standard deviations 4=106 V (A), 8=106 V (B) and 16=106 V (C). The lower resistivity value of 0.002 V m was used in all cases.
Fig. 2 shows the I–V characteristics obtained using a normal distribution of filament resistances P(r), having a mean resistance rms8=107 V, and with three different standard deviations increasing from 4=106 V (A), through 8=106 V (B) to 16=106 V (C); these respectively represent rm y20, rm y10 and rm y5. The resistivity was 0.002 V m in all cases. At low voltages ohmic conduction is observed, as would be expected if no filaments are ruptured, since we assume ohmic conductivity in each individual filament. The distribution with the smallest standard deviation (A) yields the sharpest VCNR region, since in this case nearly all the filaments rupture over a very small voltage range. As the standard deviation increases the I–V characteristic becomes more realistic, until for curve C there is a fairly wide region of falling current between 5 and 9 V. Thus, fairly good agreement is obtained with experimental measurements, and with characteristics obtained using a normal distribution of resistances w6x, but without the use of a Monte Carlo simulation. Differences between curves presented in ref. w6x and the present work are attributable to differences in the parameters of the normal distribution adopted. Fig. 3 illustrates characteristics obtained for normally distributed filament radii, having a mean value rms 15.3=10y10 m and standard deviations of 7.5=10y11 m (A), 1.03=10y10 m (B) and 3.8=10y10 m (C), corresponding approximately to rm y20 (A), rm y15 (B) and rm y4 (C). Again the resistivity value adopted was 0.002 V m. Similar behaviour to that obtained in Fig. 2 was obtained, with low voltage ohmic behaviour and more developed VCNR for the distribution with the largest standard deviation. More importantly, these curves exactly replicate those obtained analytically for the same filament parameters w7x, demonstrating that the present simulation method is both viable and accurate.
324
S. Gravano et al. / Thin Solid Films 433 (2003) 321–325
Fig. 3. Current–voltage characteristics derived from normal distributions of filament radii, with mean radius 15.3=10y10 m and standard deviations of 7.5=10y11 m (A), 1.03=10y10 m (B) and 3.8=10y10 m (C). The lower resistivity value of 0.002 V m was used in all cases.
Fig. 4 shows characteristics obtained using identical normal filament radii distributions to those of Fig. 3, but with the exception that the assumed filament resistivity is 0.0034 V m. This higher resistivity filament material would be appropriate, for instance, in the case of filaments containing carbonaceous material rather than the electrode metal. As expected, the currents for each of the three characteristics are lower than for the corresponding characteristic of Fig. 3, as a result of the higher resistivity material. There also appears to be a shift to higher voltages of the current peak, e.g. the peak of curve C shifts from approximately 5 to 6 V. This is consistent with the following expression for the voltage at which a filament of radius r will cease to conduct w7 x : b B r0d E1y2 F Vps C rD p G
Fig. 4. Current–voltage characteristics derived from normal distributions of filament radii, with mean radius 15.3=10y10 m and standard deviations of 7.5=10y11 m (A), 1.03=10y10 m (B) and 3.8=10y10 m (C). The higher resistivity value of 0.0034 V m was used in all cases.
of the lower resistivity, the main effect is an extension of the VCNR region to higher voltages, especially for the smaller standard deviations. This clearly reflects the importance of the higher resistance filaments, which as discussed earlier, continue conducting to higher applied voltage levels. There is also a slight reduction in current levels at lower voltages, reflecting the higher resistance of some of the filaments. For distributions containing a smaller proportion of lower conductivity filaments (i.e. 50 and 10%) it was observed that for the smaller standard deviations turning points were present in the I–V characteristics. Although this may be a possible cause of the instabilities in the negative resistance region, such characteristics are not generally observed experimentally. It is clear from the work of Ray and Hogarth w13x that realistic I–V characteristics are gen-
(6)
This predicts that the voltage Vp at which a filament will cease to conduct is proportional to yr0, where r0 is the filament resistivity. Rupturing of filaments at higher voltages for Fig. 4 than for Fig. 3, inevitably results in a shift of the current peak to higher voltages. In order to determine the effects on the I–V characteristics of having two different filament resistivities, some based on gold electrode material and some based on carbonaceous material, a series of simulations was performed in which the proportion of filaments present varied from 90% lower resistivityy10% higher resistivity to 10% lower resistivityy90% higher resistivity. An example is shown in Fig. 5, in which 90% of the filaments were of the lower resistivity and 10% of the higher resistivity. The mean and standard deviations of the filament radii distributions were as in Figs. 3 and 4. Comparing Fig. 5 to Fig. 3, in which all filaments are
Fig. 5. Current–voltage characteristics derived from normal distributions of filament radii, with mean radius 15.3=10y10 m and standard deviations of 7.5=10y11 m (A), 1.03=10y10 m (B) and 3.8=10y10 m (C). 90% of the filaments had the lower resistivity value of 0.002 V m and 10% the higher resistivity value of 0.0034 V m.
S. Gravano et al. / Thin Solid Films 433 (2003) 321–325
325
method, requiring the consideration of two filament distributions. The utility of the approach in considering even more complex situations with additional filament distributions, and particularly with multivariate distributions, is obvious. 4. Conclusions
Fig. 6. Current–voltage characteristics derived from a normal distribution of filament radii, with mean radius 15.3=10y10 m and standard deviation 3.8=10y10 m. The captions indicate the proportion of lower resistivity (0.002 V m) to higher resistivity (0.0034 V m) filaments.
erated for the higher standard deviations, and we do not therefore consider these effects to be significant. In Fig. 6 three I–V characteristics are shown which are based on the usual mean radius rms15.3=10y10 m and with the highest standard deviation of 3.8=10y10 m (rm y4). As discussed above, these are more likely to yield realistic characteristics than those with smaller standard deviations. The proportions of lower to higher resistivity filaments are indicated on the figure; thus the curve labelled 90%: 10% corresponds to curve C of Fig. 5. Fig. 6 confirms that in moving from a lower to a higher resistivity value the peak current level decreases, and the current peak shifts to higher voltages. It is worth mentioning that all the computed characteristics initially show ohmic behaviour, notwithstanding the faster than ohmic experimental behaviour illustrated in Fig. 1. This follows from the approximation made that the conductivity through each individual filament is ohmic. In order to resolve this discrepancy it is clear that non-linear conductivity, of the type mentioned in Section 1 w3x needs to be considered in future simulations. This approach has been adopted in the work of Sutherland w14x for ZnS samples and clearly lends itself to the Monte Carlo approach considered here. The above examples demonstrate that the Monte Carlo simulation method gives similar results to analytic calculations, while being simpler to use. In particular, the results presented in Figs. 5 and 6 would be more cumbersome to replicate without using a simulation
It has been demonstrated that a Monte Carlo simulation allows the generation of I–V characteristics for electroformed structures which are in agreement with previous analytical work and with experiment, with the exception of the initial linear behaviour. The method has been tested for normal probability density distributions of filament resistances and radii, the latter also for filaments having two different resistivity values. Provided the filament distributions are suitably defined, the method is also applicable to non-normal probability density distributions. Further work aims to utilise the method for the investigation of samples having multivariate probability density distributions, and to apply it to the case of non-linear conduction through individual filaments. Acknowledgments The authors (particularly E. Amr) would like to acknowledge the financial support of the British Council in East Jerusalem towards a joint Masters’ programme between Al Quds and Keele Universities. References w1x T.W. Hickmott, J. Appl. Phys. 35 (1964) 2679. w2x J.G. Simmons, R.R. Verderber, Proc. Roy. Soc. (A) 301 (1967) 77. w3x G. Dearnaley, D.V. Morgan, A.M. Stoneham, J. Non-Cryst. Solids 4 (1970) 593. w4x R.D. Gould, C.A. Hogarth, Phys. Stat. Sol. (A) 23 (1974) 531. w5x R.D. Gould, Thin Solid Films 57 (1979) 33. w6x A.K. Ray, C.A. Hogarth, Thin Solid Films 127 (1985) 69. w7x S. Gravano, R.D. Gould, Int. J. Electron. 73 (1992) 315. w8x S. Gravano, R.D. Gould, Thin Solid Films 248 (1994) 263. w9x A.K. Ray, C.A. Hogarth, Int. J. Electron. 69 (1990) 97. w10x R.D. Gould, J. Non-Cryst. Solids 55 (1983) 363. w11x H. Gould, J. Tobochnik, An Introduction to Computer Simulation Methods, Part 2, Addison-Wesley, New York, 1988. w12x H. Pagnia, N. Sotnik, Phys. Stat. Sol. (A) 108 (1988) 11. w13x A.K. Ray, C.A. Hogarth, Thin Solid Films 141 (1986) 201. w14x R.R. Sutherland, J. Phys. D: Appl. Phys. 4 (1971) 468.
Monte Carlo simulation of current–voltage characteristics in metal– insulator–metal thin film structures S. Gravanoa,*, E. Amrb, R.D. Goulda, M. Abu Samrab a Thin Films Laboratory, Department of Physics, School of Chemistry and Physics, Keele University, Keele, Staffs, ST5 5BG, UK Department of Physics, College of Science and Technology, Al Quds University, Abu Deis, P.O. Box 20002, Jerusalem, West Bank
b
Abstract Many metal–insulator–metal (MIM) thin film sandwich structures exhibit an electroforming process, following which the resulting direct current (DC) current–voltage I–V characteristic exhibits ohmic conduction and voltage-controlled differential negative resistance (VCNR) at low and high voltages, respectively. One of the more successful models of these phenomena is the filamentary model, in which it is assumed that the electroforming process establishes a population of ohmic filaments within the insulating matrix, which span the metal contacts. The VCNR behaviour results from the progressive cessation of conduction in individual filaments owing to Joule heating effects. Early work showed that by postulating plausible probability distributions of filament resistances, I–V characteristics typical of those obtained experimentally could be derived. More recently more realistic normal distributions of filament resistances, radii and cross-sectional areas have been considered, and these too have yielded characteristics in accordance with experiment. In the present work normal distributions of filament resistances and radii are considered, but in contrast to previous work a novel method using a Monte Carlo simulation was employed. Furthermore, the effects of using a combination of two different filament resistivity values were explored. Variations in conductance and shifts in the current peak were explicable using existing theory. The approach adopted in this work lends itself to the exploration of more complex filament distributions, and therefore to a better representation of the underlying filament population, including those having a multivariate distribution. 䊚 2003 Elsevier Science B.V. All rights reserved. Keywords: Monte Carlo simulation; Electrical properties; Metal-insulator-metal structures (MIM)
1. Introduction The current–voltage (I–V) characteristics shown by metal–insulator–metal (MIM) thin film sandwich structures yield useful information concerning the dominant conduction process through the insulator. Under certain circumstances a phenomenon known as electroforming occurs. Essentially this consists of applying a direct current (DC) voltage, in excess of a few volts, between the metal electrodes. Providing this voltage exceeds a level known as the forming voltage VF, electroforming occurs and the sample shows various characteristic features. These include a considerably increased conductance, voltage-controlled differential negative resistance (VCNR) behaviour, electron emission and various memory effects. Early models of the conduction process in electroformed structures were based on the establishment *Corresponding author.
of a modified band structure during electroforming w1,2x. Following these a novel filamentary model was proposed by Dearnaley et al. w3x in which the conductivity is dominated by a matrix consisting of a large number of parallel conducting paths, or filaments. At low voltage the conduction is essentially ohmic, but as the voltage is increased Joule heating within the filaments causes them to cease to conduct or rupture, and the current level falls in the negative resistance region of the VCNR characteristic. Fig. 1 w4x shows a typical current–voltage characteristic for a Au–CaBr2 –Au sandwich structure. Ic represents the current circulating through the sample and Ie the current emitted as functions of the applied voltage Ub. We are not concerned with emission current in this simulation work. This characteristic is representative of those for many different insulator and electrode combinations. In common with other characteristics presented in the literature, the circulating current initially increases with voltage faster than expected from Ohm’s
0040-6090/03/$ - see front matter 䊚 2003 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 0 - 6 0 9 0 Ž 0 3 . 0 0 3 8 7 - 0
322
S. Gravano et al. / Thin Solid Films 433 (2003) 321–325
tribution of resistances is then given by
PŽr.s
Fig. 1. A typical experimental current–voltage characteristic for an electroformed Au–CaBr2 –Au sample. Ic represents the current circulating through the sample, Ie the current emitted and Ub the applied voltage. Reproduced from w4x. Reprinted with permission from R.D. Gould and C.A. Hogarth, Phys. Stat. Sol. (A) 23, 531. Copyright Wiley-VCH, 1974.
law, following typically either an IAI0sinh(VyV0 ) or an IAI0(VyV0)2 dependence, where I is the current through the sample, V is the applied voltage, and I0 and V0 are constants w3x. The conducting filaments were initially characterised in terms of a probability density distribution of filament resistances P(r) which was triangular, where r represents the filament resistance w3x. Later work on the form of the function P(r) suggested that a parabolic distribution equally well described the experimental results w5x, while it was also argued that a normal distribution of filament resistances was more fundamental w6x, and that the earlier suggested distributions were merely approximations to this. More recently the present authors have suggested that there is no a priori reason for suggesting that the filament resistances are normally distributed, and that more basic attributes such as the filament radius w7x or cross-sectional area w8x may be normally distributed. Indeed, both of these assumptions have been shown to yield computed characteristics in accord with experiment. Ray and Hogarth w6,9x have also suggested that more than one variable may be normally distributed, for instance the resistivity and the filament radius. The degree of mathematical difficulty for analytical solutions, even for a single distribution of filament crosssectional areas, becomes increasingly complex, and numerical methods are normally required w8x. For multivariate distributions a more easily applicable method is clearly required, and the present work is an attempt to develop such a method. Before presenting the general results obtained in such initial simulations, it is necessary to outline the filamentary model for normal distributions of filament resistances and radii, and to briefly describe the simulation. 2. The filamentary conduction model and simulation We consider first the situation when the filament resistances r are normally distributed. The normal dis-
w
1
y2psr2
expxy y
Žryrm.2 z
|
2s2r
~
(1)
where rm is the mean and sr is the standard deviation. A filament ruptures when the power dissipated in the filament reaches a particular value b2, which is determined by the heat flow symmetry from the filament, and depends primarily on the filament dimensions and density, the filament and bulk temperatures and the insulator thermal conductivity w10x. b2 is effectively a factor indicating how the power is dissipated, and it follows from this that the voltage Vp at which a filament of resistance r will rupture is given by Vpsbr1y2
(2)
The conductance G(s1yR), where R is the resistance, of an assembly of N filaments when a voltage V is applied is given by ru
1 sN R
|
PŽr. rv
ry
(3)
dr
where ru represents the maximum filament resistance and rv represents the value of the lowest resistance in the distribution which is still conducting at voltage V, and may be calculated from Eq. (2). The lowest resistance in the distribution is denoted by rt. When VFbr1y2 t , all the filaments in the distribution are conducting, and the lower limit of the integral rv is given 1y2 by rt. For br1y2 t FVFbru , where some filaments are conducting the lower limit rvsV 2 y b2, and for VGbr1y2 where no filaments conduct rvsru. u If we then consider the filament radii r to be normally distributed, the normal distribution of radii is given by
PŽr.s
1
y2psr2
w
expxy y
Žryrm.2 z 2s2r
| ~
(4)
where rm is the mean and sr is the standard deviation. Although expressions analogous to Eq. (2) and Eq. (3) also exist for the distribution of filament radii w7x, it was not necessary to use these in the simulation. All that was necessary was to calculate the resistance of each filament, using rs
r0d pr2
(5)
where r0 is the filament resistivity and d is the filament length.
S. Gravano et al. / Thin Solid Films 433 (2003) 321–325
In the simulation the current I was calculated for different values of voltage V. For a specified voltage V filaments with resistance r)V 2 y b2 exhibit ohmic conduction and contribute to the overall current, whereas those with rFV 2 y b2 do not conduct. Since we have made the approximation of ohmic conduction through individual filaments, it is clear that the initial I–V dependence will be linear where all filaments are conducting. The simulation was based on a Monte Carlo approach in which for each filament a value of either resistance or radius was randomly generated using the ¨ Box–Muller method w11x. The resistance or radius value for the filament was determined by specifying the mean and standard deviation of P(r) or P(r), as appropriate. For simulations involving the radius, the resistance was calculated from Eq. (5). If r)V 2 y b2 the filament contributes a current Vy r to the overall current I. If rFV 2 y b2 the filament does not conduct, and makes no contribution to the overall current. The total current flowing with voltage V applied was obtained by summing the currents through the individual filaments. This was then repeated over a range of voltages to generate the complete I–V characteristic. 3. Results and discussion The results of a selection of the simulations are shown in the following figures. Two fixed filament resistivity values r0 were used. The lower of these is that adopted by Dearnaley et al. w3x of 0.002 V m, which assumes that the resistivity is determined largely by that of metallic gold atoms derived from the electrodes, which are injected into the insulating matrix in the filamentary regions. A higher resistivity value of 0.0034 V m was also used in some simulations. This value was chosen to reflect an alternative view of filamentary conduction, in which carbon is considered to be an important component of the filamentary material w12x. Note that these resistivity values do not correspond directly to those of gold and carbon, but are chosen (i) to yield realistic current values for the assumed number of filaments N, and (ii) to reflect different filament compositions, depending on the different origins of a major filament component. Following previous work on filamentary conduction, other parameters of the model were as proposed by Dearnaley et al. w3x. In particular these were the number of filaments Ns5=106, filament length ds4=10y8 m and heat flow constant bs 2.1=10y3 W1y2. These values were also used in earlier work concerned with normally distributed filament radii w7x, and therefore allow direct comparison with this work. In order to demonstrate the feasibility of the method, initial work was performed using a normal distribution of filament resistances, while the bulk of the work used a more realistic normal distribution of filament radii.
323
Fig. 2. Current–voltage characteristics derived from normal distributions of filament resistances with mean resistance 8=107 V, and standard deviations 4=106 V (A), 8=106 V (B) and 16=106 V (C). The lower resistivity value of 0.002 V m was used in all cases.
Fig. 2 shows the I–V characteristics obtained using a normal distribution of filament resistances P(r), having a mean resistance rms8=107 V, and with three different standard deviations increasing from 4=106 V (A), through 8=106 V (B) to 16=106 V (C); these respectively represent rm y20, rm y10 and rm y5. The resistivity was 0.002 V m in all cases. At low voltages ohmic conduction is observed, as would be expected if no filaments are ruptured, since we assume ohmic conductivity in each individual filament. The distribution with the smallest standard deviation (A) yields the sharpest VCNR region, since in this case nearly all the filaments rupture over a very small voltage range. As the standard deviation increases the I–V characteristic becomes more realistic, until for curve C there is a fairly wide region of falling current between 5 and 9 V. Thus, fairly good agreement is obtained with experimental measurements, and with characteristics obtained using a normal distribution of resistances w6x, but without the use of a Monte Carlo simulation. Differences between curves presented in ref. w6x and the present work are attributable to differences in the parameters of the normal distribution adopted. Fig. 3 illustrates characteristics obtained for normally distributed filament radii, having a mean value rms 15.3=10y10 m and standard deviations of 7.5=10y11 m (A), 1.03=10y10 m (B) and 3.8=10y10 m (C), corresponding approximately to rm y20 (A), rm y15 (B) and rm y4 (C). Again the resistivity value adopted was 0.002 V m. Similar behaviour to that obtained in Fig. 2 was obtained, with low voltage ohmic behaviour and more developed VCNR for the distribution with the largest standard deviation. More importantly, these curves exactly replicate those obtained analytically for the same filament parameters w7x, demonstrating that the present simulation method is both viable and accurate.
324
S. Gravano et al. / Thin Solid Films 433 (2003) 321–325
Fig. 3. Current–voltage characteristics derived from normal distributions of filament radii, with mean radius 15.3=10y10 m and standard deviations of 7.5=10y11 m (A), 1.03=10y10 m (B) and 3.8=10y10 m (C). The lower resistivity value of 0.002 V m was used in all cases.
Fig. 4 shows characteristics obtained using identical normal filament radii distributions to those of Fig. 3, but with the exception that the assumed filament resistivity is 0.0034 V m. This higher resistivity filament material would be appropriate, for instance, in the case of filaments containing carbonaceous material rather than the electrode metal. As expected, the currents for each of the three characteristics are lower than for the corresponding characteristic of Fig. 3, as a result of the higher resistivity material. There also appears to be a shift to higher voltages of the current peak, e.g. the peak of curve C shifts from approximately 5 to 6 V. This is consistent with the following expression for the voltage at which a filament of radius r will cease to conduct w7 x : b B r0d E1y2 F Vps C rD p G
Fig. 4. Current–voltage characteristics derived from normal distributions of filament radii, with mean radius 15.3=10y10 m and standard deviations of 7.5=10y11 m (A), 1.03=10y10 m (B) and 3.8=10y10 m (C). The higher resistivity value of 0.0034 V m was used in all cases.
of the lower resistivity, the main effect is an extension of the VCNR region to higher voltages, especially for the smaller standard deviations. This clearly reflects the importance of the higher resistance filaments, which as discussed earlier, continue conducting to higher applied voltage levels. There is also a slight reduction in current levels at lower voltages, reflecting the higher resistance of some of the filaments. For distributions containing a smaller proportion of lower conductivity filaments (i.e. 50 and 10%) it was observed that for the smaller standard deviations turning points were present in the I–V characteristics. Although this may be a possible cause of the instabilities in the negative resistance region, such characteristics are not generally observed experimentally. It is clear from the work of Ray and Hogarth w13x that realistic I–V characteristics are gen-
(6)
This predicts that the voltage Vp at which a filament will cease to conduct is proportional to yr0, where r0 is the filament resistivity. Rupturing of filaments at higher voltages for Fig. 4 than for Fig. 3, inevitably results in a shift of the current peak to higher voltages. In order to determine the effects on the I–V characteristics of having two different filament resistivities, some based on gold electrode material and some based on carbonaceous material, a series of simulations was performed in which the proportion of filaments present varied from 90% lower resistivityy10% higher resistivity to 10% lower resistivityy90% higher resistivity. An example is shown in Fig. 5, in which 90% of the filaments were of the lower resistivity and 10% of the higher resistivity. The mean and standard deviations of the filament radii distributions were as in Figs. 3 and 4. Comparing Fig. 5 to Fig. 3, in which all filaments are
Fig. 5. Current–voltage characteristics derived from normal distributions of filament radii, with mean radius 15.3=10y10 m and standard deviations of 7.5=10y11 m (A), 1.03=10y10 m (B) and 3.8=10y10 m (C). 90% of the filaments had the lower resistivity value of 0.002 V m and 10% the higher resistivity value of 0.0034 V m.
S. Gravano et al. / Thin Solid Films 433 (2003) 321–325
325
method, requiring the consideration of two filament distributions. The utility of the approach in considering even more complex situations with additional filament distributions, and particularly with multivariate distributions, is obvious. 4. Conclusions
Fig. 6. Current–voltage characteristics derived from a normal distribution of filament radii, with mean radius 15.3=10y10 m and standard deviation 3.8=10y10 m. The captions indicate the proportion of lower resistivity (0.002 V m) to higher resistivity (0.0034 V m) filaments.
erated for the higher standard deviations, and we do not therefore consider these effects to be significant. In Fig. 6 three I–V characteristics are shown which are based on the usual mean radius rms15.3=10y10 m and with the highest standard deviation of 3.8=10y10 m (rm y4). As discussed above, these are more likely to yield realistic characteristics than those with smaller standard deviations. The proportions of lower to higher resistivity filaments are indicated on the figure; thus the curve labelled 90%: 10% corresponds to curve C of Fig. 5. Fig. 6 confirms that in moving from a lower to a higher resistivity value the peak current level decreases, and the current peak shifts to higher voltages. It is worth mentioning that all the computed characteristics initially show ohmic behaviour, notwithstanding the faster than ohmic experimental behaviour illustrated in Fig. 1. This follows from the approximation made that the conductivity through each individual filament is ohmic. In order to resolve this discrepancy it is clear that non-linear conductivity, of the type mentioned in Section 1 w3x needs to be considered in future simulations. This approach has been adopted in the work of Sutherland w14x for ZnS samples and clearly lends itself to the Monte Carlo approach considered here. The above examples demonstrate that the Monte Carlo simulation method gives similar results to analytic calculations, while being simpler to use. In particular, the results presented in Figs. 5 and 6 would be more cumbersome to replicate without using a simulation
It has been demonstrated that a Monte Carlo simulation allows the generation of I–V characteristics for electroformed structures which are in agreement with previous analytical work and with experiment, with the exception of the initial linear behaviour. The method has been tested for normal probability density distributions of filament resistances and radii, the latter also for filaments having two different resistivity values. Provided the filament distributions are suitably defined, the method is also applicable to non-normal probability density distributions. Further work aims to utilise the method for the investigation of samples having multivariate probability density distributions, and to apply it to the case of non-linear conduction through individual filaments. Acknowledgments The authors (particularly E. Amr) would like to acknowledge the financial support of the British Council in East Jerusalem towards a joint Masters’ programme between Al Quds and Keele Universities. References w1x T.W. Hickmott, J. Appl. Phys. 35 (1964) 2679. w2x J.G. Simmons, R.R. Verderber, Proc. Roy. Soc. (A) 301 (1967) 77. w3x G. Dearnaley, D.V. Morgan, A.M. Stoneham, J. Non-Cryst. Solids 4 (1970) 593. w4x R.D. Gould, C.A. Hogarth, Phys. Stat. Sol. (A) 23 (1974) 531. w5x R.D. Gould, Thin Solid Films 57 (1979) 33. w6x A.K. Ray, C.A. Hogarth, Thin Solid Films 127 (1985) 69. w7x S. Gravano, R.D. Gould, Int. J. Electron. 73 (1992) 315. w8x S. Gravano, R.D. Gould, Thin Solid Films 248 (1994) 263. w9x A.K. Ray, C.A. Hogarth, Int. J. Electron. 69 (1990) 97. w10x R.D. Gould, J. Non-Cryst. Solids 55 (1983) 363. w11x H. Gould, J. Tobochnik, An Introduction to Computer Simulation Methods, Part 2, Addison-Wesley, New York, 1988. w12x H. Pagnia, N. Sotnik, Phys. Stat. Sol. (A) 108 (1988) 11. w13x A.K. Ray, C.A. Hogarth, Thin Solid Films 141 (1986) 201. w14x R.R. Sutherland, J. Phys. D: Appl. Phys. 4 (1971) 468.