15832
J. Phys. Chem. C 2007, 111, 15832-15838
Alternate Current Nonequilibrium Molecular Dynamics Simulations of Yttria-Stabilized Zirconia† Qingyin Zhang and Kwong-Yu Chan* Department of Chemistry, the UniVersity of Hong Kong, Pokfulam Road, Hong Kong S.A.R., China ReceiVed: May 28, 2007; In Final Form: August 13, 2007
An alternating current nonequilibrium molecular dynamics (AC-NEMD) simulation technique is applied to study the ionic conduction behavior of 8 mol % yttria-stabilized zirconia (YSZ). A maximum conductivity is observed at a frequency of 3 × 1013 Hz when the current response is in phase with the applied field. The current response lags the applied field at lower frequencies but leads the applied field at higher frequencies, in both cases having lower conductivity than when current response is in phase. Represented by a Cole-Cole plot of complex admittance representation, the AC-NEMD data indicate a resistance-capacitance (RC) circuit behavior at low frequency and resistance-inductance (RL) circuit behavior at high frequency. The frequency dependent ionic conductivity may be correlated to the length scale and time scale of oxygen transport between vacancy sites. In the AC-NEMD, the generated ohmic heat is effectively removed, and a constant temperature is maintained by an Anderson thermostat. Separate NEMD simulations performed with a Nose-Hoover thermostat give mostly identical results, within statistical error bars.
1. Introduction Yttria-stabilized zirconia (YSZ) is an important electrolyte material used in solid oxide fuel cells (SOFC) and oxygen sensors. The cubic fluorite phase of ZrO2 appears only above 2370 °C1 but can be stabilized at room temperature by doping the solid electrolyte with lower valency oxides, such as Y2O3. The high conductivity of YSZ is due to some Y3+ ions replacing Zr4+ in the fluorite structure of zirconia. The lower valency of Y3+ leads to vacancy sites in the cubic structure, which are originally occupied by oxygen ions. These vacancies allow oxygen ions to hop through and facilitate oxygen transport. The ionic conductivities of YSZ are 0.02 Ω-1 cm-1 at 800 °C and 0.1 Ω-1cm-1 at 1000 °C.2 In recent years, YSZ has been intensively studied by molecular dynamics simulations.3-21 Simulation research has focused on the self-diffusion coefficients (or conductivities) of oxygen ions, the activation energy of diffusion, and the structures of oxygen vacancies. Shimojo et al.3,4 used an equilibrium molecular dynamics (EMD) method to investigate the oxygen diffusions of 4.85, 10.2, and 22.7 mol % YSZ and found that 10.2 mol % YSZ had the maximum oxygen diffusion coefficient. Okazaki et al.5 proposed an ordered YSZ structure in which the oxygen diffusion coefficient was larger than that in a random YSZ structure. Yamamura et al.10 calculated the ionic conductivity of YSZ at concentration between 5.9 and 12.5 mol % Y2O3. They observed that 8 mol % YSZ exhibited a maximum conductivity. Fisher and Matsubara11,13 investigated the oxygen diffusion through grain boundaries of YSZ by EMD simulations. In the scale of simulations studied, they found that grain boundaries were detrimental to ionic conductivity. Sawaguchi et al.14 studied the diffusion coefficients of oxygen ions and the activation energy of diffusion. Their results showed that the activation energy of oxygen ions increased linearly with yttria concentrations. Kilo †
Part of the “Keith E. Gubbins Festschrift”. * Corresponding author. E-mail:
[email protected]. Fax: (852) 2859 7919.
et al.16 used both experimental technique and molecular dynamics simulations to study the oxygen diffusion in YSZ. The temperature was 650-1200 K, and the yttria concentration was 8-24 mol %. Their experimental and simulation results all showed that 10 mol % YSZ had the highest diffusion coefficient, but the diffusion results obtained by EMD simulations were higher than those of experiments. Arima et al.17 studied the transport properties of YSZ between 300∼2000 K. Except for diffusion properties, the constant pressure heat capacity was also obtained through the EMD simulations. Hayashi et al.20 used EMD simulations to study the thermal expansion coefficients of YSZ with the Y3+ content. Devanathan et al.21 investigated the diffusion coefficient of oxygen ions at 1125∼2500 K. They found that 8 mol % YSZ had the highest coefficient, and the activation energy of diffusion linearly increased with increasing Y2O3 content. Most simulations of solid oxide electrolytes reported were performed with the equilibrium molecular dynamics method. In EMD simulations, there is no net flux of oxygen ions, as opposed to the situation of a current flow in the operation of a solid oxide electrolyte fuel cell. Only the self-diffusion of ions is calculated and ionic conductivity is evaluated by the NernstEinstein relationship between diffusion coefficient and conductivity. To simulate the situation of current flow, Tang et al.19 first employed the nonequilibrium molecular dynamics (NEMD) method with a direct current (DC) to study the conductivity of YSZ. In their results, the Arrhenius activation energy obtained by NEMD was in better agreement with experimental data compared with the results of EMD simulations. In their other simulation study of aqueous electrolyte confined in a nanopore, Tang et al.22 observed a large discrepancy between EMD and NEMD results showing the Nernst-Einstein relationship not to be appropriate in a nanoconfined system. The DC nonequilibrium molecular dynamics (DC-NEMD) method is therefore more realistic than the EMD method in simulating current flow in a solid oxide electrolyte. In experiments, it is also very common to use alternating current (AC) and impedance
10.1021/jp0741152 CCC: $37.00 © 2007 American Chemical Society Published on Web 10/03/2007
MD Simulations of Yttria-Stabilized Zirconia
J. Phys. Chem. C, Vol. 111, No. 43, 2007 15833
techniques to investigate electrode-electrolyte interfaces and determine properties of electrolyte. Relaxation dynamics in multiple time scales can be probed by the application of AC and the analyses in frequency domain. The use of alternating current nonequilibrium molecular dynamics (AC-NEMD) method in simulating oxygen ion transport will be attractive and offers a new angle to probe the properties of the solid oxide electrolyte. The use of AC-NEMD has been applied to simulate an aqueous electrolyte confined in a nanopore.23,24 The modes of ionic transport are markedly different in a solid oxide electrolyte compared with aqueous electrolyte in which more than one ion species are mobile in the liquid phase. To further understand the behavior of YSZ and the intricate dynamics of ionic transport, it will be interesting to apply the AC-NEMD method. We report here, to our best knowledge, the first AC-NEMD simulations of a solid electrolyte. 2. Simulation Method
Figure 1. Unit cell structure of YSZ (gray spheres, Zr4+; black spheres, Y3+; shaded circles, O2-; crossed circle, oxygen vacancy).
In the AC-NEMD simulations, a frequency-dependent alternating electric field E(ω) is applied along the Z direction of the simulation box. The equation of motion can be expressed as
the real parts of the electric field and current density, respectively. Hence,
q3 i )
pi m
p3 i ) Fi + qiEz(ω)kˆ
1 V
N
qiVi,z ∑ i)1
(2)
where Vi,z is the axial component of the velocity of ith oxygen ion, qi is the valence charge of the ion, and V is the volume of the simulation box. The relationship of frequency dependent conductivity, σ(ω), current density, J(ω,), and electrical field, E(ω), can be expressed as
σ(ω) ) J(ω)/E(ω)
(3)
Here, the bold face symbols denote complex numbers representing current density, electric field, and conductivity. In simulations, the external AC electric field was applied in the Z direction. Without loss of generality, we used Ez(ω) and Jz(ω) to denote the z components of the electric field and current density, respectively. In complex space, the electrical field can be expressed as
Ez(ω) ) |Ez|ei(ωt-π/2)
(4)
where |Ez| is the modulus of the electric field at the particular frequency ω. In the linear response range, the current density can be similarly represented as
Jz(ω) ) |Jz|ei(ωt-π/2+φ)
(6a)
Jz ) Jzo sin(ωt + φ)
(6b)
and
(1b)
where the bold letters pi and qi are momentum and position vectors of ion i, m is its mass, Fi is the interaction force on the ion, qi is the charge of the ion, Ez is the strength of the applied field, ω is the frequency of external alternate current, and k is the unit vector in the same direction as the applied field. The current density of oxygen ions, Jz(ω), is calculated by
Jz(ω) )
Ez ) Ezo sin(ωt)
(1a)
(5)
where φ represents the phase lag between the electric field and the current field. For convenience, Ez and Jz are to represent
In eq 6, Ezo and Jzo denote the amplitudes of the electric field and the current density, respectively. It can be seen from eqs 4 and 5 that they are also equal to the moduli |Ez| and |Jz| at the given frequency. The real functions Ez and Jz are the physical quantities implemented and measured in the NEMD simulations reported in this study. In the linear response range, the complex conductivity can be obtained by eq 3 and resolved into real and imaginary parts, denoted by σ′(ω) and σ′′(ω), respectively. The real and imaginary parts of the electrical conductivity are computed in the equations.23
σ′(ω) ) Re[σ(ω)] ) lim Re
Ez0f0
[
]
Jz0 ei(ωt+φ-π/2) i(ωt-π/2)
Ez0 e
) lim
Ez0f0
Jz0 cos φ (7a) Ez0
and
σ′′(ω) ) lim
Ez0f0
Jz0 sin φ Ez0
(7b)
In eq 7a,b, the phase difference between external field and response current can be measured by fitting the current density with a sine function and comparing it with the electric field function. Here, it is assumed that the current response is monochromatic with the same frequency as the applied field. The 8 mol % YSZ has a cubic fluorite type structure with the cations forming a face-centered-cubic (FCC) arrangement. The simulation cell is composed of 3 × 3 × 3 unit cells, containing 92 Zr4+, 16 Y3+, and 208 O2- ions. The lattice of zirconium(IV) oxide was first constructed and the positions of yttrium ions were generated by random replacements of zirconium ions. To keep the neutrality of the simulation box, a corresponding number of oxygen ions were removed accordingly. A unit structure of 8 mol % YSZ is shown in Figure 1.
15834 J. Phys. Chem. C, Vol. 111, No. 43, 2007
Zhang and Chan
TABLE 1: Interatomic Potential Parameters i-j
Aij (kJ mol-1)
Bij (m-1)
Cij (kJ nm6 mol-1)
Zr-Zr Zr-Y Y-Y Zr-O Y-O O-O
0 0 0 9.512 × 104 12.978 × 104 2.196 × 106
1 1 1 2.660 × 1010 2.865 × 1010 6.711 × 1010
0 0 0 0 0 2.691 × 10-3
In the AC-NEMD simulations, a sinusoidal electric field with a fixed frequency was imposed in the Z direction. The current response was monitored by tracing the velocity of oxygen ions with time. The interaction between ions consists of short-range Born-Meyer-Buckingham (BMB) potential and the long-range Coulomb interaction potential. The BMB potential is expressed as
U(rij) ) Aij e-Bijrij -
Cij rij6
(8)
where rij is the distance between the ions i and j; Aij, Bij, and Cij are the parameters of the interaction potential. The parameters used are identical to previously adopted values in refs 25 and 26 and are given in Table 1. The interaction between cations is assumed to be negligible with Aij identically equal to zero. The corresponding choices of Bij are arbitrarily assigned to be unity to be consistent in other presentations of eq 8 where the reciprocal coefficients are used in the exponents. The Coulomb interaction has the form
Uc(rij) )
q i qj 4π0rij
(9)
where 0 is the permeability of free space, and qi ) zie and qj ) zie are the charges of the ions. The charge valencies zi, of zirconium, yttrium, and oxide ions are +4, +3, and -2, respectively. Ewald summations27 are used to calculate the longrange Coulomb interactions between ions in the simulation cell. The total potential contains short repulsive potential between ionic cores, van der Vaals interaction, and coulomb interaction. Thus, this potential is based on the rigid ion approach, as in almost all previously reported MD work. In principle, the use of a shell model and allowing polarization of oxygen ions is a better approach. In a previous shell model MD study of CaF2, however, it was found that the polarizability had little effect on diffusion dynamics in the superionic state.28 Temperatures of the simulations here were at 900, 1273, and 1759 K. The thermal expansion was accounted for by using the same number of ions but different sizes of the simulation cell. The length of the simulation cell, a(T), increases with temperature and is given by the equation
a(T) ) a(T ) 300 K)[1 + R(T - 300)]
(10)
where a(T ) 300 K) ) 5.136 Å and R ) 1 × 10-5 K-1 is the thermal expansion coefficient for all of the compositions of yttria-stabilized zirconium.8,29 The use of the isobaricisothermal ensemble (NPT) will be a better approach and avoid the inconsistency of potential model and such an artificial volume adjustment. The previous EMD simulations using the same BMB potential model16 showed agreement in the expansion coefficient in eq 10 with experimental results.28 They have, however, small adjustments of less than 1% in the cationanion interaction parameters Aij and Bij compared with those in Table 1.
Figure 2. Radial distribution functions of three different atom pairs in 8 mol % YSZ at 1759 K in an AC electric field of strength E* ) 3 × 10-6 (lines, ω* ) 0.07193; solid symbols, ω* ) 0.04316; open symbols, ω* ) 0.01726).
TABLE 2: Phase Shifts of Current Response in 8 mol % YSZ with Different Electric Field Strengths at 1759 K Simulated with an Anderson Thermostata ω*
NT (fs)
φ1
∆φ1
0.01726 0.02158 0.02877 0.03453 0.04316 0.04795 0.05395 0.06166 0.07193
500 400 300 250 200 180 160 140 120
-47.5 -52.1 -30.0 -17.5 0.900 20.5 47.0 61.3 83.6
(21.2 (25.2 (31.0 (13.7 (10.3 (5.63 (3.05 (4.56 (5.03
φ2
∆φ2
-61.1 (24.9 -53.0 (22.6 -33.6 (16.5 -20.3 (8.50 5.85 (6.44 24.8 (3.53 45.3 (2.23 63.9 (1.77 82.7 (2.93
φ3
∆φ3
-54.8 (20.8 -46.0 (17.3 -30.9 (9.51 -21.6 (7.89 4.28 (2.34 25.8 (2.25 44.2 (2.39 66.5 (0.972 85.3 (2.36
a Subscripts 1, 2, and 3 represents results of E1* ) 1.0 × 10-6, E2* ) 2.0 × 10-6, E3* ) 3.0 × 10-6, respectively.
Because an AC field was imposed, the heat generated by the external field must be removed efficiently. Both Anderson30 and Nose-Hoover31,32 thermostats were separately used. In the Anderson thermostat, the collision frequency ν was 0.00045 in reduced unit (ref 19). In the Nose-Hoover thermostat, the parameter Q ) 0.5 (reduce unit) was used. The time step of simulation was set to 1 fs. Both thermostats controlled the desired temperature well. 3. Results and Discussion The first sets of NEMD simulations are performed at 900, 1273, and 1759 K under an Anderson thermostat control. The three different choices of amplitude for the external sinusoidal electric field are 4.47, 8.94, and 13.4 × 107 V m-1. These correspond to reduced units of E* ) 1.0 × 10-6, 2.0 × 10-6, and 3.0 × 10-6, respectively. In experiments, a voltage of 1.0 can be applied over a thin film solid-electrolyte layer of 10-7m so the electric field strengths simulated are not unrealistic values. The total run time of a simulation in each case is 1 ns, corresponding to 1 million time steps. The frequency of the sine electrical field ranges from 1.26 × 1013 Hz to 5.24 × 1013 Hz corresponding to reduced units of ω* ) 0.01726 to 0.07193. These values are much larger than reported experimental frequencies for YSZ impedance measurements. At a given temperature, YSZ maintains a stable solid structure at all magnitudes and frequencies of electric fields applied, as shown in the pair distribution curves in Figure 2. The Zr-Zr and Zr-O pair distributions show pronounced isolated peaks, characteristic of a solid like structure. The oxygen-oxygen pair has a lower first peak at about 0.26 nm and a broader and more
MD Simulations of Yttria-Stabilized Zirconia
J. Phys. Chem. C, Vol. 111, No. 43, 2007 15835
Figure 3. Current density response to electric field at different frequencies and magnitudes under an Anderson thermostat (solid line is current response and dotted line is waveform of original field applied; increasing frequency from left to right, ω* ) 0.1726, 0.4316, 0.7193; increasing electric field strengths from top to bottom, E* ) 1 × 10-6, 2 × 10-6, 3 × 10-6).
Figure 5. Phase shift of currents in 8 mol % YSZ versus frequency of electric field at different temperatures (solid symbols, E* ) 1 × 106; open symbols, E* ) 3 × 106; error bars are statistical deviations in a simulation run). Figure 4. Comparison of current response between 4 × 4 × 4 and 3 × 3 × 3 cell under an Anderson thermostat. (Increasing electric field strengths from top to bottom: E* ) 1 × 10-6, 2 × 10-6, 3 × 10-6.)
connected second peak at about 0.37 nm. The oxygen distribution is more continuous, because of its larger number and the existence of vacancy sites. This continuous distribution also favors oxygen transport, and almost all of the current is carried by the oxygen anions. Figure 3 shows the current responses at
different combinations of electrical fields and frequencies at 1759 K. Comparing the series of subfigures in Figure 3, some general trends of the effect of frequency and electric field can be observed. The panels from right to left show increasing frequency from ω ) 0.01726 to 0.07193 (whereas from top to bottom the electric field strength increases from E*) 1 × 10-6 to 3 × 10-6). At low frequency (left panel), the current response lags the electric field as the solid lines have peaks to the left of the potential peaks. The locations of the current response peaks
15836 J. Phys. Chem. C, Vol. 111, No. 43, 2007
Zhang and Chan
Figure 6. Field dependent conductivities of 8 mol % YSZ simulated with an Anderson thermostat at (a) 900, (b) 1273, and (c) 1759 K.
are shown as labeled numbers in the subfigures. At the intermediate frequency (middle panel), the current response is in phase with the applied sinusoidal electric field, whereas at higher frequency (right panel), the current response leads the electric waveform. A more detailed and complete comparison of the phase shifts are tabulated in Table 2 and shown graphically in Figure 5 for results over a range of temperature 900 to 1759 K. In checking the size effect of simulation cell, we employ additional simulations with a cell containing 4 × 4 × 4 unit cells for the case of 1759 K with ω* ) 0.4316. The comparison of current responses in simulations with a 4 × 4 × 4 cell and a 3 × 3 × 3 cell is shown in Figure 4. Both sets of simulations give very similar results. In Table 2, NT is defined as the number of time steps (or in femtoseconds) for one period of the applied sinusoidal electric field. It is interesting to observe the general trend of phase shifts from negative values (phase lag) at low frequencies to positive values (phase lead) at high frequencies. This result is similar to the earlier result observed in liquid electrolyte confined in a nanopore. The negative phase shift corresponds to a capacitance behavior whereas a phase lead corresponds to an inductance behavior. For a confined electrolyte, the one-dimensional capacitor was explained as accumulation of ion pairs at low frequencies when oppositely traveling ions cannot get past each other in a severely confined nanopore.23,24 YSZ is a bulk electrolyte and the confinement effect is less obvious. We can, however postulate that current is mainly carried by oxygen anions along a path of oxygen vacancies. The passage between vacancies is restricted with the presence of an energy barrier. Depending on length scale and time scale, this can lead to accumulation of charges when ions are driven by the external field over a longer time scale at a low frequency. The time scale for passage between vacancies can be roughly estimated. According to eq 10, the simulation box length is 5.167, 5.186, and 5.211 Å at 900, 1273, and 1759 K, respectively. The distance between oxygen vacancies can be estimated if one oxygen vacancy exists in a unit cell, as shown in Figure. 1. The time required to travel between oxygen vacancies without barrier can be estimated from the average thermal velocity of oxygen ion at different temperatures, without considering migration due to the electric field. These estimates are 381, 319, and 272 fs, for 900, 1273, and 1759 K, respectively. At 1759 K, the 272 fs corresponds roughly to the value of Nτ where zero phase-shift was observed. The regularity of the current waveform also depends on the frequency. In the left panel of a lower frequency, the current response curves show larger relative fluctuations with multiple peaks and shoulders appearing. The large error bars are observed at low frequency in Figure 5. This may be explained by the tortuous passage
Figure 7. Frequency dependence conductivity of 8 mol % YSZ simulated with an Anderson thermostat.
between oxygen vacancies, which is more evident at a low frequency when ions travel a longer distance before reverting direction at the turn of electric field. The relative fluctuations and noises, however, can be reduced with a stronger field applied as observed by comparing top and bottom subfigures in Figure 3. These noises are partially due to thermal fluctuations, which are random in three directions but less apparent with a stronger electric field applied in the z direction. In Figure 5, the phase lag appears to be smaller at further lower frequencies and tends toward the origin (direct current) at the zero frequency limitation. Unfortunately, this cannot be concluded because of the large scattering of data and insufficient CPU time to simulate low frequency behavior. Conductivity of 8 mol % YSZ can be determined by eq 3 using the data of current responses at different field strengths. The conductivity determined assuming a linear response to a given field strength at a fixed frequency is shown in Figure 6. Fairly consistent values of conductivities are determined at different electric field strengths, indicating good linear response behavior, and the results can be extrapolated to obtain the zero-field conductivity. The imaginary part of admittance is calculated in the same manner. The real part of admittance (conductivity) thus obtained is plotted in Figure 7 as a function of frequency. The electrical conductivity of YSZ has a maximum value at a certain frequency. This frequency corresponds to that of zero phase-shift in Figure 5. This maximum conductivity may be closely related to the optimum mode of oxygen transport between vacancy sites. In a given cycle of an applied waveform at a low frequency, an oxygen anion travels a large distance, beyond the next vacancy site. This can result in an ionic arrangement with increased polarity and undesirable configuration during field reversal. On the other
MD Simulations of Yttria-Stabilized Zirconia
Figure 8. Cole-Cole plot of frequency dependent admittance in 8 mol % YSZ simulated with an Anderson thermostat.
Figure 9. Cole-Cole plots of frequency dependent conductivities in 8 mol % YSZ simulated with a Nose-Hoover thermostat.
hand, at a high frequency, an oxygen anion vibrates about its equilibrium position without being connected to the vacancy sites. Both of these situations do not carry much current. The maximum AC conductivity observed here is 2 orders of magnitude higher than the DC conductivity observed experimentally and in previous DC-NEMD results.19 The lowfrequency conductivity in Figure 7 does tend toward a very small number. It appears that the phenomenon of high-frequency AC ionic conduction can be markedly different from the DC conduction. The effect of frequency on conductivity can also be interpreted by representation in the complex domain, indicating capacitance and inductance effects with negative and positive phase shifts, respectively. The complex representation is called admittance, and Figure 8 plots the imaginary part versus the real part of admittance for the 8 mol % YSZ simulated at three different temperatures. This σ′-σ′′ plot or Cole-Cole plot is analogous to the Z′-Z′′ plot of impedance in the complex domain, when impedance is the reciprocal of admittance. The plot of complex impedance is a popular experimental technique called impedance spectroscopy (IS).33 In IS, an equivalent electrical circuit with inductance-resistance-capacitance (LRC) components is often proposed to fit the resulting Z′-Z′′ plot and interpret the corresponding physical phenomena that occurred. In Figure 8, admittance starts at low frequency in the fourth quadrant with a negative complex component indicating the presence of a capacitance. At an intermediate frequency, the maximum conductivity (real part) is reached with a zero phase shift when the curve crosses the x axis. At the highest frequencies, the admittance curve resides in the first quadrant,
J. Phys. Chem. C, Vol. 111, No. 43, 2007 15837 and an inductance component is indicated. The admittance approaches the origin at the infinite frequency limit. Most experimental impedance spectroscopy data show only a semicircle in the first quadrant, which is equivalent to a semicircle of admittance in the fourth quadrant of the Cole-Cole plot. The equivalent circuit of YSZ and other solid-state electrolytes interpreted from experimental IS are dominated by resistancecapacitance (RC) components without any inductance-resistance (LR) components. A number of factors can contribute this difference between the AC-NEMD results here and the experimental results. Earlier impedance spectroscopy data was reported by Bauerle33,34 for stabilized zirconia interfaced with a platinum electrode, and the range of frequency studied was below 100 kHz. In recent years, the frequency range increased to 1 MHz but is still much below the frequency studied here in NEMD. There is a big gap between the lowest AC frequency that can be simulated and the highest frequency applied in experiments. Recent experiments also focused more on the performance of the anode or cathode35-39 and the electrode-electrolyte interface40-42 whereas our present NEMD study is focused on a uniform bulk YSZ material without electrode-electrolyte interfaces. Another important difference is the presence of grain boundaries in YSZ bulk material studied experimentally. The relaxation frequency of bulk lattice is 1 or 2 orders of magnitude higher that in between grain boundaries.43 These differences make it difficult to have a direct comparison of the present NEMD results and experimental data. Nevertheless, the ACNEMD simulation presented here is a useful step toward the comprehensive modeling of solid oxide electrolyte at the atomic scale. The higher frequency results can be useful in understanding relaxation phenomena at shorter length scales and time scales. Setting up domain structures and including grain boundaries in AC-NEMD are also possible. Previous equilibrium molecular dynamics simulations have included artificially setting up lattice arrangements5 and grain boundaries11,13 and have found higher conductivity along certain paths. On the other hand, the modeling of electrolyte-electrode interfaces to match experimental impedance results have only been made at the continuum level.44 Applying a multi-scale modeling approach may be needed as an intermediate step to close the gap between theory and experiments. The Anderson thermostat successfully maintains the temperature of the NEMD simulations. To investigate possible effects of the thermostat, separate NEMD simulations are carried out using the alternative Nose-Hoover thermostat. For all temperatures, frequencies, and electric fields, the resulting current density responses are quite similar in separate NEMD simulations using the two thermostats. Figure 9 shows the Cole-Cole plot of admittance data obtained from NEMD simulations using a Nose-Hoover thermostat. Within the statistical error bars, the results are mostly the same as those obtained with the Anderson thermostat in Figure 8. 4. Conclusions An AC nonequilibrium molecular dynamics (AC-NEMD) simulation technique is successfully applied for the first time to study the ionic conduction behavior of 8 mol % yttriastabilized zirconia (YSZ). A maximum conductivity is observed at a frequency of 3 × 1013 Hz when the current response is in phase with the applied field. The current response lags the applied field at lower frequencies and leads the applied field at higher frequencies and, in both cases, have a lower conductivity than when current response is in phase. The frequency dependent ionic conductivity may be correlated to the length scale and
15838 J. Phys. Chem. C, Vol. 111, No. 43, 2007 time scale of relaxation between oxygen vacancy sites. Because of the differences in frequency range and the presence of electrode-electrolyte interface and grain boundaries in experiments, a direct comparison between NEMD results and experiments cannot be made. Ohmic heat is effectively removed, and a constant temperature is maintained by both an Anderson thermostat and a Nose-Hoover thermostat. NEMD simulations with inclusion of grain boundaries and electrode electrolyte interfaces may be made in the future to bridge the gap between experiments and atomistic simulations. Acknowledgment. This research has been supported by a Hong Kong Research Grants Council grant (HKU 7006/05P). References and Notes (1) Smith, D. K.; Cline, C. F. J. Am. Ceram. Soc. 1962, 45, 249. (2) Larminie, J.; Dicks, A. Fuel Cell Systems Explained; John Wiley & Sons Ltd: Chichester, England, 2003. (3) Shimojo, F.; Okabe, T.; Tachibana, F.; Kobayashi, M.; Okazaki, H. J. Phys. Soc. Jpn. 1992, 61, 2848. (4) Shimojo, F.; Okazaki, H. J. Phys. Soc. Jpn. 1992, 61, 4106. (5) Okazaki, H.; Suzuki, H.; Ihata, K. Phys. Lett. A 1994, 188, 2910. (6) Okazaki, H.; Suzuki, H.; Ihata, K. J. Phys. Soc. Jpn. 1994, 63, 3556. (7) Brinkman, H. W.; Briels, W. J.; Verweij, H. Chem. Phys. Lett. 1995, 247, 386. (8) Li, X. Y.; Hafskjold, B. J. Phys.: Condens. Matter 1995, 7, 1255. (9) Khan, M. S.; Islam, M. S.; Bates, D. R. J. Mater. Chem. 1998, 8, 2299. (10) Yamamura, Y.; Kawasaki, S.; Sakai, H. Solid State Ionics 1999, 126, 181. (11) Fisher, C. A. J.; Matsubara, H. Solid State Ionics 1998, 115, 311. (12) Suzuki, K.; Kubo, M.; Oumi, Y.; Miura, R.; Takaba, H.; Fahmi, A.; Chatterjee, A.; Teraishi, K.; Miyamoto, A. Appl. Phys. Lett. 1998, 73, 1502. (13) Fisher, C. A. J.; Matsubara, H. Molecular dynamics investigations of grain boundary phenomena in cubic zirconia. Comput. Mater. Sci. 1999, 14, 177. (14) Sawaguchi, N.; Ogawa, H. Simulated diffusion of oxide ions in YO1.5-ZrO2 at high temperature. Solid State Ionics 2000, 128, 183. (15) Perumal, T. P.; Sridhar, V.; Murthy, K. P. N.; Easwarakumar, K. S.; Ramasamy, S. Physica A 2002, 309, 35.
Zhang and Chan (16) Kilo, M.; Argirusis, C.; Borchardt, G.; Jackson, R. A. Phys. Chem. Chem. Phys. 2003, 5, 2219. (17) Arima, T.; Fukuyo, K.; Idemitsu, K.; Inagaki, Y. J. Mol. Liq. 2004, 113, 67. (18) Kilo, M.; Taylor, M. A.; Argirusis, C.; Borchardt, G.; Jackson, R. A.; Schulz, O.; Martin, M.; Weller, M. Solid State Ionics 2004, 175, 823. (19) Tang, Y. W.; Zhang, Q. Y.; Chan K. Y. Chem. Phys. Lett. 2004, 385, 202. (20) Hayashi, H.; Saitou, T.; Maruyama, N.; Inaba, H.; Kawamura, K.; Mori, M. Solid State Ionics 2005, 176, 613. (21) Devanathan, R.; Weber, W. J.; Singhal, S. C.; Gale, J. D. Solid State Ionics 2006, 177, 1251. (22) Tang, Y W.; Szalai, I.; Chan, K. Y. J. Phys. Chem. A 2001, 105, 9616. (23) Tang, Y. W.; Chan, K. Y.; Szalai, I. Nano. Lett. 2003, 3, 217 (24) Tang, Y. W.; Szalai, I.; Chan, K. Y. Mol. Phys 2002, 100, 1497. (25) Lewis, G. V.; Catlow, C. R. A. J. Phys. C: Solid State Phys. 1985, 18, 1149. (26) Dwivedi, A.; Cormack, A. N. Philos. Mag. A 1990, 61, 1. (27) Ewald, P. P. Ann. Phys. 1921, 64, 253. (28) Lindan, P. J. D.; Gillan, M. J. J. Phys.: Condens. Matter 1993, 5, 1019. (29) Ingel, R. P.; Lewis, D. J. Am. Ceram. Soc. 1986, 69, 325. (30) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384. (31) Nose, S. J. Chem. Phys. 1984, 81, 511. (32) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (33) Barsoukov, E.; Macdonald, J. R. Impedance Spectroscopy, Theory, Experiment, and Applications; John Wiley & Sons, Inc.: Hoboken, NJ, 2005. (34) Bauerle, J. E. J. Phys. Chem. Solids 1969, 30, 2657. (35) Bebelis, S.; Neophytides, S. Solid State Ionics 2002, 152-153, 447. (36) Lee, C. H.; Lee, C. H.; Oh, S. M. Solid State Ionics 1997, 98, 39. (37) Kek, D.; Panjan, P.; Wanzenberg, E.; Jamnik, J. J. Euro. Ceram. Soc. 2001, 21, 1861. (38) Kim, J. D.; Kim, G. D.; Moon. J. W.; Park. Y. I.; Lee, W. H.; Kobayashi, K.; Nagai, M.; Kim, C. E. Solid State Ionics 2001, 143, 379. (39) Barbucci, A.; Bozzo, R.; Cerisola, G.; Costamagna, P. Electrochim. Acta 2002, 47, 2183. (40) Sridhar, S.; Stancovski, V.; Pal, U. B. Solid State Ionics 1997, 100, 17. (41) Bay, L.; Jacobsen, T. Solid State Ionics 1997, 93, 201. (42) Wanzenberg, E.; Tietz, F.; Kek, D.; Panjan, P.; Stover, D. Solid State Ionics 2003, 164, 121. (43) Morales, J. C. R.; Marrero-Lopez, D.; Irvine, J. T. S.; Nunez, P. Mater. Res. Bull. 2004, 39, 1299. (44) Bessler, W. G. Solid State Ionics 2005, 176, 997.
