Authors: Yukiteru Murakami, Katharina Krischer
Categories: Article, electrochemistry, applied sciences, semiconductor manufacturing
Source: iScience
Authors: Yukiteru Murakami, Katharina Krischer
Silicon is a key semiconductor electrode and is known for exhibiting oscillations during anodic oxidation in fluoride-containing solutions. This article introduces a mathematical model that captures these oscillations. It describes the formation and dissolution of an oxide layer, determining the oxide composition and the electrostatic potential perpendicular to the electrode. Oscillations occur if the following three conditions are the etching speed increases with defects in the oxide layer, the defect density decreases with increasing electric field at the Si-oxide interface, and a sufficient time delay exists between production and etching of the oxide. Numerical simulations reproduce experimental results well. Based on these results, we derive a simplified time-delay model. Using linear stability analysis, we confirm the essential role of the time delay for the oscillations. The basic mechanisms are universal and in line with the point defect model for growth and dissolution of passive films on metal electrodes.
Because of its importance in the electronics industry and its application in electrochemical devices, such as sensors or photoelectrochemical cells, silicon is arguably the semiconductor material whose electrochemistry has been most thoroughly studied among all semiconductors.^1^ Yet, there are many fundamental processes occurring at the electrified Si/electrolyte interface that are still only poorly understood. Here the situation seems little changed from twenty years ago, when Macdonald, in the preface to the book “Electrochemistry of Si and its oxides,”^1^ pointed out that “… the lack of a comprehensive account of the electrochemistry of silicon in aqueous solution at the fundamental level is surprising and troubling […] some of the electrochemical properties of this element are not as well known as might be warranted by its importance in a modern technological society.” One of these still puzzling, yet intensively studied phenomena is self-sustained oscillations in current or potential in a fluoride-containing electrolyte, first reported nearly 70 years ago by Turner.^2^ The last decades revealed that these oscillations are connected to a variety of further nonlinear phenomena, such as birhythmicity^3^ and even bichaoticity,^4^ i.e., the coexistence of different types of stable oscillations or chaotic attractors, at the same parameter values, or the formation of spatiotemporal patterns.^5^^,^^6^ Two notable examples of the latter are chimera states, i.e., synchronization patterns in ensembles of identical oscillators or uniform oscillating media that are composed of a synchronously and an incoherently oscillating part,^7^ and frequency clusters,^8^ that are composed of a few regions that oscillate with distinct frequencies. In both cases, Si electrooxidation is one of the few experimental systems where the spontaneous formation of such patterns has been observed, making the “Si oscillator” a prototypical system for the experimental exploration of pattern formation or the validation of theories. Understanding the mechanism of current oscillations during the electrodissolution of Si is therefore important both for our fundamental interest in the dynamics of nonlinear systems and, furthermore, for controlling or designing the uniformity and quality in silicon oxidation and etching processes used in the nanofabrication industry.
As mentioned above, the study of Si oscillations has a long history. Turner already pointed out that, during the oscillations, a passivating silicon oxide film covers the electrode surface. The oxide film is formed electrochemically and chemically etched by fluoride species according to reactions (R0a) and (R0b), respectively.(R0a)Si+2H2O+λVBh+→SiO2+4H++(4−λVB)e−,(R0b)SiO2+6HF→SiF62−+2H2O+2H+,where 1≤λVB≤4 is the number of the valence band holes taking part in the oxidation reaction.^9^ Since holes are required for the oxidation process, the reaction (R0a) proceeds with either p-type silicon or illuminated n-type silicon. Under potentiostatic conditions, the oscillatory oxidation current is observed over a wide range of parameters.
From the outset, it was speculated that the oscillations are caused by some kind of interplay between oxide film formation and dissolution. This triggered a large number of investigations into the film properties during oscillations using various in situ and ex situ techniques, such as IR spectroscopy,^10^^,^^11^ ellipsometry^12^ atomic force microscopy,^13^ X-ray reflectometry,^13^ microwave reflectivity,^11^ and transmission electron microscopy.^12^ An overview of experimental studies on current oscillations up to 2001 can be found in chapter five of Electrochemistry of Silicon and its Oxide.^1^ The studies showed that many properties of the oxide film oscillated along with the current or potential, most notably the thickness of the oxide layer.^10^ Furthermore, it was observed that pores in the oxide opened and closed during oscillations.^12^
Most theoretical studies were motivated by the latter observation of opening and closing pores during oscillations. Many of them dealt with the question of how local, microscopically small oscillating domains synchronize to give rise to macroscopic oscillations in current density, not addressing the physicochemical mechanism of how the local oscillations occur in the first place.^14^^,^^15^ An exception hereof is the current-burst model suggested by Föll and coworkers. The model assumes that the ionic current through the oxide layer flows through pores that open at a certain strength of the electric field across the oxide layer, allowing for oxide growth. The pores close again only at a considerably lower field strength and only reopen when the oxide has been sufficiently etched such that the critical field strength for pore opening is reached again. Monte-Carlo simulations of the current-burst model agree in many respects with experimentally observed oscillations. However, today we know that this picture is incomplete. The system exhibits birhythmicity, and thus, there are two distinct types of oscillation patterns, so-called high-amplitude oscillation (HAO) and low-amplitude oscillation (LAO).^4^^,^^16^ HAOs are observed under relatively high voltages and typically have a relaxation-type character, while LAOs are sinusoidally shaped and occur at lower voltages. Their coexistence in certain parameter intervals underlines that their oscillation mechanisms also differ. The current-burst model seems to describe HAOs.
The importance of a delay for oscillations was hypothesized by Smith and Collins, who assumed the formation of “hard” and “soft” oxides that are etched at different rates.^17^ While this idea already captured the intuitive role of delayed feedback on the system, it remained conceptual and did not identify the specific physicochemical processes responsible for generating the delay, nor did it provide evidence. Duportal et al.^18^ demonstrated that the experimentally observed impedance spectra are consistent with a defect-mediated picture of oxide development and dissolution, suggesting that the Si oscillations can be captured with some type of a point-defect model. The present work builds directly on these ideas and experimental findings by developing a minimal physicochemical model that explicitly resolves defect generation, transport, and etching, thereby establishing a direct link between point-defect dynamics, time delay, and oscillatory behavior.
In this article, we present a mathematical model for LAOs during Si electrodissolution that involves just a few electrochemical and chemical reaction steps, yet predicts the important dependencies on the experimental parameters. The model inherently contains a well-defined, physically grounded time delay as the key ingredient that destabilizes the steady state and gives rise to sustained oscillations. More precisely, this basic physicochemical model predicts oscillations when the following three conditions are (1) the dissolution rate of the oxide film increases with an increasing number of partially oxidized silicon species, below also called defects, in the oxide; (2) the fraction of partially oxidized species decreases under a high electric field at the Si/oxide interface because they become fully oxidized; and (3) there is a sufficient time delay between the instant at which the defect is produced at the Si/oxide interface and the instant at which it is exposed to the electrolyte, and thus, etched. In other words, oscillations require a sufficiently thick oxide film. The interplay between etch rate and defect formation rate constitutes a negative feedback. However, the time delay suppresses the negative feedback for a certain time, causing an effective positive feedback that promotes stable oscillation. The fundamental oscillation mechanism is illustrated in Figure 1. It does not seem to be specific to Si electrodissolution but may also apply to other semiconductor or metal electrodes in which anodic formation of oxide layers and their dissolution are accompanied by self-sustained oscillations.Figure 1A conceptual diagram of the silicon-oscillation mechanismThe vertical axis represents the coordinate perpendicular to the oxide layer, and the horizontal axis represents time. A change in electric field strength at the Si/oxide interface alters the average oxidation state of the Si oxide (the number of defects). After a time delay, the altered composition is exposed to the electrolyte interfaces, resulting in a different etching rate. This, in turn, prompts a change in thickness and, consequently, a change in electric field strength that causes again an altered composition of the oxide. This delayed response weakens the stabilizing negative feedback, providing the physical origin of oscillatory behavior.
Delay-induced oscillations, in which a stable state becomes oscillatory once a negative feedback is sufficiently delayed, are well known in biological systems. To give two examples besides the famous Mackey-Glass equation,^19^ we mention population growth as an early study,^20^ or gene expression as a more recent one.^21^ This places the present findings in a general framework of time-delayed nonlinear dynamics, rather than being relevant only to silicon electrochemistry.
The article is organized as In Model derivation, we derive the mathematical model that describes the dynamics of the oxide film. We first introduce the governing equations for the potential and the different Si-oxide species in the oxide layer (Governing equations within the Si-oxide layer). In Boundary conditions for the chemical species at the Si/oxide interface and Boundary conditions for the chemical species at the oxide/electrolyte interface, we specify the boundary conditions for the Si-oxide species at the Si/oxide and the oxide/electrolyte interfaces, respectively. Then, we explain how we define a moving frame that accounts for the time-dependent motions of the two interfaces (Interface motion and moving frame), before we derive the boundary conditions for the electrostatic potential in Boundary conditions for the electrostatic potential. This section is concluded by a discussion of the numerical implementation of the model. In Simulations, numerical simulations of the model are presented together with a comparison to experimental findings under Potentiostatic conditions and Potentiodynamic conditions. A detailed analysis of the temporal evolutions of the various simulated quantities allows us then to deduce the oscillation mechanism in Discussion. Based on the insights of the oscillation mechanism, we derive a simplified model that consists of two coupled delay differential equations (Simplified model). Simulation results of this simplified model together with a linear stability analysis are presented in Behavior of the simplified model. In the closing paragraph of the Discussion, a summary of the results and a short outlook are provided. At the end, we mention Limitations of the study. Experimental conditions and parameter values used in the simulations can be found in STAR Methods.
In the following, we derive the mathematical model for Si electrodissolution in fluoride containing electrolytes. The model calculates the temporal evolution of the state of the oxide layer bounded by two moving interfaces, the oxide/electrolyte and the Si/oxide interfaces. For this purpose, we consider the following charge-transfer reactions at the Si/oxide interface, diffusion and migration of defects within the oxide layer, and ion transfer as well as chemical etching reactions at the oxide/electrolyte interface. In mathematical terms, we have to solve continuity equations for the chemical species in the oxide layer together with Poisson’s equation, which determines the electrostatic potential within the oxide. The reactions at the two interfaces are part of the boundary conditions, as is the potentiostatic control of the experiment.
The oxidation of Si at the Si/oxide interface occurs in many reaction steps and involves various oxide species with Si oxidation states from I to IV. Furthermore, substoichiometric oxides may also be further oxidized within the oxide layer. For the sake of simplicity, we consider only three species in the model and neglect charge-transfer processes in the oxide. The species considered are fully oxidized silicon oxide, SiO2, one substoichiometric oxide species, SiO, and a charged defect in the oxide lattice, SiO2+, in other words, an O2− vacancy. These species are transported through the oxide and react further at the oxide/electrolyte interface.
A positively charged oxygen vacancy moves by exchanging its position with neighboring O2− ions, thereby carrying the current through the oxide and simultaneously leading to the buildup of SiO2 in the opposite direction. Owing to the electric field and the number density gradient of SiO2+, oxygen vacancies are transported toward the oxide/electrolyte interface, which is equivalent to a transport of O2− species to the Si/oxide interface, and amounts as well to a net transport of SiO2 toward the Si/oxide interface. This is depicted schematically in Figure 2. Finally, at the oxide/electrolyte interface, oxygen vacancies are filled through rapid reaction with water, and SiO and SiO2 are chemically etched. The processes are summarized in Figure 3, and will be further motivated in the following subsections. The scheme can be seen as a type of point defect model, as originally developed for metal/oxide/electrolyte interfaces.^22^ In the following subsections, we will derive the corresponding mathematical model.Figure 2Schematic illustration of oxygen-vacancy transport through the SiO2 latticeSiO2+ transport is equivalent to a transport of SiO2 in the opposite direction.Figure 3Schematic illustration of the considered reaction steps at the two interfacesAt the Si/oxide interface, electrochemical oxidation (R1) and (R2). At the oxide/electrolyte interface, vacancy filling (R3) and chemical etching, (R4) and (R5), occur.
We assume that the oxide layer is laterally uniform, and thus, consider only variations along the spatial coordinate z perpendicular to the oxide layer. The continuity equations for the three considered species read as (Equation 1)∂[SiO2+]∂t=DSiO2+∂2[SiO2+]∂z2+μSiO2+∂∂z[SiO2+]∂ϕ∂z(Equation 2)∂[SiO]∂t=DSiO∂2[SiO]∂z2(Equation 3)[SiO2]=ρatom−[SiO]−[SiO2+]The positive z direction points from the electrolyte into the direction of the Si bulk. The number density of species Xi is indicated by the square bracket notation [Xi], and ϕ is the electrostatic potential. DXi is the diffusion coefficients of species Xi in the oxide, and μSiO2+ the mobility of oxygen vacancies. ρatom is the total Si atomic density of the SiO2 lattice.
The first term of Equation 1 describes diffusion, and the second one migration of oxygen vacancies according to the Nernst-Planck approximation. As a neutral species, SiO is only subject to diffusion (Equation 2). The number density of SiO2 follows the algebraic mass balance Equation 3. For simplicity and clarity, the volume expansion that occurs when silicon is oxidized to silicon oxide is not taken into account.
The electrostatic potential ϕ obeys Poisson’s (Equation 4)−∂2ϕ∂z2=2q[SiO2+]ϵ0ϵoxwhere q, ϵ0, and ϵox are the elementary charge, the vacuum permittivity, and the relative permittivity of silicon oxide, respectively. The relative permittivity of silicon oxide is assumed to be independent of its composition.
In the following sections, we specify the boundary conditions for the chemical species and the electrostatic potential at the two interfaces that describe Si electrodissolution at moderate voltages.
At the Si/oxide interface, the oxide grows into the Si bulk by the generation of substoichiometric Si oxide through electric field-driven “hopping” of oxygen ions from the Si oxide lattice into the space-charge layer of the Si, leaving positively charged oxygen vacancies behind. The substoichiometric oxides are, in part, further oxidized. As a result, the SiO2 lattice expands toward the silicon bulk, with a high number density of positively charged oxygen vacancies at the Si/oxide interface. In our minimal reaction scheme, we describe these processes by the following two (R1)Si+SiO2+2h+→SiO2++SiO(R2)SiO→SiO2++2e−.It is known that the initial charge-transfer step involves valence-band holes, h+, while the further oxidation steps can occur through hole capture or electron injection.^23^ Here, we consider an electron-injection step in (R2).
The kinetic equations for the rates of formation of SiO and oxygen vacancies according to (R1) and (R2), r1 and r2, respectively, are derived according to the following The rate of the first oxidation step of a Si atom according to (R1) is proportional to the number density of SiO2 species and of valence-band holes at the Si/oxide interface, which is located at zb. The number density of valence-band holes is proportional to the electric field at the interface, Einter= ∂ϕ/∂z|b.^24^ Furthermore, we assume that the transfer of oxide ions into the space-charge layer is an activated process that depends exponentially on Einter. We therefore express the oxidation rate r1 as (Equation 5)r1=k1Einter[SiO2]|zbexp(α1Einter)with the rate constant k1 and the temperature dependent constant α1.
The further oxidation of SiO through electron injection (R2) depends on the surface-state density and its occupation level. Electron injection through the surface states occurs only when the energy of occupied surface states reaches the bottom edge of the conduction band of the silicon bulk. The occupation level of a surface state increases as −qϕ increases, which happens when the oxide becomes thinner, and therefore, the electric field at the interface becomes larger. We chose to express the occupation level of the surface states through the potential drop across the oxide monolayer adjacent to the Si bulk, uSiO2|zb. This is schematically depicted in Figure 4. Furthermore, we assume that the number of occupied surface states is limited toward higher energies, which we described by a sigmoidal dependence of the oxidation rate r2 on uSiO2|zb. We therefore approximate the rate r2 by(Equation 6)r2=k2[SiO]|zb1+exp(α2(α3−uSiO2|zb)),where k2 is a rate constant and the constants α2,α3 determine the shape of the sigmoidal function. k1 and k2 are defined so that the dimension of r1 and r2 are 1/(m2s).Figure 4Schematic comparison of the degree of filling of surface states at different values of uSiO2|zbThe filled and unfilled dots represent electrons and holes, respectively. In case (A), where uSiO2|zb is small, the highest level of occupied surface states is lower than the conduction band energy of silicon bulk, and electron injection through the surface states is slow. On the other hand, in case (B), where uSiO2|zb is large, the high energy of the occupied surface states causes (R2) to proceed rapidly.
The boundary conditions for the chemical variables [SiO] and [SiO2+] at the Si/oxide interface, i.e., at z=zb, are given by the fluxes of these species through the interface, which result from the reaction rates Equation 5 and Equation 6.(Equation 7)∂[SiO2+]∂zzb=−r1+r2DSiO2+(Equation 8)∂[SiO]∂zzb=−r1−r2DSiO
At the oxide/electrolyte interface, the oxygen vacancies are filled through the reaction with water according to (R3). In this way, oxygen ions enter the oxide and SiO2 is formed.(R3)SiO2++H2O→fastSiO2+2H+The overall electrooxidation reaction of Si, which is usually written as (R0a) thus consists of charge-transfer reactions at the Si/oxide interface and ion-transfer reactions at the oxide/electrolyte interface, in our case (R0a) = (R1) + (R2) + 2×(R3) with λVB=2.
The ion-transfer reaction (R3) is assumed to be so fast that the number density of oxygen vacancies in the oxide layer adjacent to the electrolyte, hereafter after the top oxide layer, described with the subscript zt, is assumed to be (Equation 9)[SiO2+]|zt=0
The boundary condition for [SiO] is easiest expressed in a moving frame (see Interface motion and moving frame). We denote the spatial coordinate in the moving frame by z˜. As will become clear below, in the moving frame, the flux of SiO species at the oxide/electrolyte interface is zero. The boundary condition for the number density of [SiO] thus reads as (Equation 10)∂[SiO]∂zzt˜=0
So far, we have specified the boundary conditions at the interfaces Si/oxide and oxide/electrolyte, or, expressed more formally, at z=zb and z=zt. However, both interfaces, and consequently the positions of zb and zt, move in time, which makes solving the model numerically very complicated. This challenge can be overcome when using a moving frame.
The velocity vox with which the Si/oxide interface moves toward the Si-bulk (in the positive z-direction in the laboratory frame) is determined by the rate with which SiO is formed (R1), and thus, reads(Equation 11)vox=r1ρatom.
At the oxide/electrolyte interface, the Si oxide is etched chemically by fluoride species from the electrolyte. We do not formulate the stoichiometry of individual etching steps here since they are irrelevant for the model. Rather, we assume that substoichiometric Si oxides are etched faster than SiO2^17^(R4)SiO2|zt→HF,HF2ketchSiO2SiF62−aq(R5)SiOzt→HF,HF2−ketchSiO>ketchSiO2SiF62−aq.Here, the ketchi are rate constants, and (SiF62−)aq indicates silicon hexafluoride ions dissolved in the electrolyte. Note that according to Equation 9, we do not have to take into account etching of SiO2+ species. Thus, the overall etch rate, vetch, measured in m/s, can be written as(Equation 12)vetch=kvetch([SiO2]zt+kSiO[SiO]|zt)kvetch is proportional to ketchSiO2. kSiO>1 defines how much faster SiO is etched than SiO2. In the simulation, a representative value of kSiO consistent with a previous report on Si electrodissolution kinetics was chosen.^25^
We define the moving coordinate z˜ so that, independent of time, the oxide/electrolyte interface is always at z˜=0 and the Si/oxide interface at z˜=1, as illustrated in Figure 5. Note that z˜ is a unit-less quantity. Hence, as long as we calculate the physical quantities in this moving frame, the positions of the two boundaries of the system are fixed. If we denote the initial oxide thickness at t=0 by L0, we can write for the oxide thickness at t=t(Equation 13)L(t)=L0+δL(t),where(Equation 14)δL(t)=∫0t(vox(t′)−vetch(t′))dt′.Analogously, we define the oxide/electrolyte interface displacement, δLetch(t) between t=0 and t=t by(Equation 15)δLetch(t)=∫0tvetch(t′)dt′.The conversion between z and z˜ can then be written as(Equation 16)z˜=1L(t)z−δLetch(t).Here, δLetch(t) represents the offset of the two frames, and 1/L(t) defines the stretching ratio of the two different coordinates.Figure 5Illustration of the connection between the laboratory frame and the moving frameLaboratory frame z is on the left and moving frame z˜ is on the right. In the moving frame, the oxide/electrolyte and Si/oxide interfaces are always at z˜=0 and 1, respectively, independent of time.
In the moving coordinate system, the time evolution of a number density at each position z˜ is caused by the local time evolution in the laboratory frame and the moving frame’s apparent effect. The relation between the two different time derivatives is expressed as follows(Equation 17)DDt=∂∂t−∂z∼∂t∂∂z∼.Here, ∂/∂t represents the derivative at the static (laboratory) z coordinate, and D/Dt that at the z˜ coordinate. It takes into account how a certain quantity is associated with changes over time due to the motion of the coordinate system. ∂z˜/∂t is given by(Equation 18)∂z˜∂t=−1L(t)(vox(t)−vetch(t))z˜+vetch(t),whereby the first term proportional to z˜ takes into account that the oxide-layer thickness changes with time, and the second one is the velocity with which the oxide/electrolyte interface moves, which also depends on time.
To complete the model and to be able to solve Poisson’s Equation 4, we need to specify the boundary conditions of the electrostatic potential. A schematic potential distribution across the Si/oxide/electrolyte layers is shown in Figure 6. The potential drops across the space-charge layer, usc, and the total oxide layer, uox, are much larger than those across the silicon bulk and the electrolyte (including the double layer). Therefore, we neglect the latter two and assume that the applied potential, uapp, is distributed over the oxide layer and the space-charge layer(Equation 19)uapp=uox+usc.Choosing the oxide/electrolyte interface as the origin of the potential scale, the boundary conditions at the oxide/electrolyte and Si/oxide interfaces are given by(Equation 20)ϕzt=0,(Equation 21)ϕzb=uapp−usc, respectively.Figure 6Schematic representation of the potential distribution across the Si/oxide/electrolyte interfaceThe applied potential uapp is distributed over the oxide layer, uox, and the space-charge layer, usc.
Next, we determine the temporal evolution of usc. Therefore, we model the Si/oxide interface by the equivalent circuit depicted in Figure 7. According to Kirchhoff’s current law,(Equation 22)izb=ire+Cscduscdt,where izb, ire, and Csc are the migration current density into the oxide, the electrochemical reaction current density, and the capacitance per unit area at the space-charge layer, respectively. The migration current density at the interface is given by(Equation 23)izb=σzbEzb,where σzb and Ezb are the conductivity and the electric field at the bottom of the oxide, z=zb. Note that we here take the derivative when approaching the boundary from the oxide side, i.e., from the +z direction, whereas the electric field at the boundary approached from the Si side (the −z direction) is denoted by Einter. The resistivity of the oxide at zb is calculated according to(Equation 24)σzb=2q([SiO2+]|zb+kleak)μSiO2+.The small constant kleak is added to avoid numerical difficulties. The reaction current density ire is defined by(Equation 25)ire=2q(r1+r2)The time evolution of usc is calculated using Equation 22. Finally, the electric field at the space-charge layer, Einter, which is necessary to calculate r1 and r2 is a function of usc and given by(Equation 26)Einter=2kBTqLDFϕ(usc)(Equation 27)Fϕ(usc)=exp−qkBTusc+qkBTusc−1,where kB, T, and LD are the Boltzmann constant, temperature, and Debye length, respectively. This relation is calculated by approximating the Si/oxide/electrolyte interface as a semiconductor-insulator-conductor connection.^24^Figure 7Equivalent circuit of the Si/oxide interfaceThe migration current density into the oxide is the sum of the electrochemical reaction current density and the space-charge layer capacitive current density.
An illustration of our spatial discretization scheme is shown in Figure 8. We discretized the oxide into n=25 grid points. Because the thickness of the oxide is fixed at 1 on the z˜ coordinate, the physical thickness of each discretized length element, lox, changes over time. For easier implementation of the boundary conditions, we define additional top and bottom layers which sandwich the discretized oxide and have a fixed length, lm.(Equation 28)lox|zt=lox|zb=lm=constFigure 8Discretization of the oxideThe physical thickness of the top and bottom layers is fixed, while that of the middle layers changes over time.
In this section, the numerical simulation results of the complete model are presented and compared with the experimental results.
The model exhibits self-sustained current and oxide-layer thickness oscillations under potentiostatic conditions, i.e., at a constant value of the applied voltage uapp. More precisely, with the parameter values used, oscillations were obtained for uapp≥4.6V while for uapp<4.6V, the etching rate exceeded the oxidation rate and no oxide layer formed.
Simulated time series of the current density A and the oxide-layer thickness B at various values of uapp in the oscillatory region are shown in Figure 9. The oscillation amplitudes and periods obviously change with the applied voltage. Furthermore, the average oxide-layer thickness increases with uapp, while the mean current density changes only marginally.Figure 9Simulation results of the full model(A) Current time series, (B) time series of oxide-layer thickness, (C) oscillation period, (D) current amplitude, and (E) average oxide-layer thickness vs applied voltage. The square root-like increase of the current amplitude is a fingerprint of a Hopf bifurcation.
The dependence of the oscillation period, the amplitude of the current oscillations, and the time-averaged oxide depth on uapp is shown in plates C, D, and E. The period of the oscillation increases linearly with uapp while the current amplitude follows a square root relationship. The nearly linear increase of the mean oxide-layer thickness with the applied potential at an almost constant average current density reflects the fact that the time-averaged electric field strength adjusts to a similar value regardless of the applied potential in a stable oscillatory state. This picture agrees with the ideas for film growth in the point defect model.^22^
For comparison, experimental results are shown in Figure 10. The measurements were conducted under potentiostatic conditions in the range between 3 and 4VvsMSE with a series external resistance of 1.77kΩcm2. Further details of the experiments are summarized in the STAR Methods. Obviously, all the characteristics of the oscillations from the simulations agree with the experimental observations. In both cases, the shape of the current oscillations is nearly harmonic while the mean current hardly changes with the applied voltage (Figure 10 A). Furthermore, the oscillation period (Figure 10 B) and the amplitudes of the current (Figure 10 C) increase linearly and square root-like with uapp, respectively. In addition, current oscillations only occur when the applied potential exceeds a threshold value of 3.2VvsMSE, and the average current level stays nearly the same. Finally, the linear increase of the average oxide thickness with the applied voltage agrees with the experimental result presented by Duportal et al.^18^Figure 10Experimental results(A) Current time series for various values of the applied voltage, (B) oscillation period, (C) current current vs applied voltage. In all experiments, a series of external resistance of 1.77kΩcm2 was inserted between the working electrode and the potentiostat. The plots show qualitatively similar results to Figure 9.
Having seen that the simulations at fixed parameter values agree well with the experimental results, we next compare experiments and simulations at cyclic voltage sweeps. Typical experimental current-voltage curves obtained at different scan speeds are depicted in Figure 11. They can be roughly divided into three regions. Negative to the current maximum at uapp=1.3VvsMSE is the electropolishing region (I) in which a surface oxide forms and dissolves at equal rates such that no stable oxide builds up on the Si surface.^1^ In the negative differential resistance region (II) between 1.3 and 2.5VvsMSE, a stable oxide film develops on the electrode surface. At the lower potential edge of that region, the oxide composition is dominated by suboxides. They are increasingly oxidized to SiO2 with increasing voltage.^25^ Region (III) in the voltage interval between 2.5 and 6.5VvsMSE is the oscillatory region that we describe in our model. Here, the current density is limited by the etching speed, and sustained current oscillations are observed at any constant value of uapp. However, when the potential is swept, in the forward scan toward higher potentials, the oscillations are suppressed or have a considerably smaller amplitude than in the backward scan. Moreover, the average current density is higher in the forward scan than in the backward scan. This hysteresis increases with the scan speed.Figure 11Experimental current voltage curves with different scan speedsThe red, green, and blue curves represent 10mV/s, 5mV/s, and 2mV/s scans, respectively. The arrows indicate the scan directions. The experimental details are given in Experimental model and study participant details. Regions (I), (II), and (III) represent electropolishing, negative differential resistance, and oscillatory regions, respectively.
Simulations of current-voltage curves during cyclic voltage sweeps in region (III) covered by the model are shown for different scan rates in Figure 12. The forward scans are shown as solid lines, and the backward scans as dashed lines. At all scan rates, the oscillation amplitudes are larger for the backward scan than for the forward scan. Furthermore, comparing the oscillation amplitudes at the different scan rates, they increase as the scan rate increases. These features show good agreement with the measurements. However, different from the experimental data, a hysteresis in the average current density between the forward and reverse scans is not clearly recognizable in the simulated curves. Nevertheless, the qualitative agreement of the trends of the oscillation amplitudes with both scan rate and scan direction can be regarded as a validation of the model.Figure 12Simulated current potential curves with different scan speedsThe red, green, and blue curves represent 28mV/s, 14mV/s, and 7mV/s scans, respectively. The forward scans are shown in solid lines, and the backward scans in dashed lines.
So far, we have discussed simulation results that can be compared to experimental observations. In this subsection, we will look at the temporal evolution of simulated quantities that provide further insight into the oscillation mechanism. Figure 13 depicts the spatiotemporal evolution of [SiO] in the oxide layer during oscillations, where the spatial coordinate z is the static laboratory coordinate. The blue upper and gray lower triangles show where the electrolyte and the Si bulk are located. It can be seen that the growth of the oxide layer into the Si bulk and the etching of the oxide by the electrolyte occur with a nearly constant velocity, while the space-time plot of the [SiO] forms approximately a horizontal stripe pattern. The latter indicates that transport of SiO through diffusion is negligible. Rather, the suboxide defects remain very close to the z position of their formation and only reach the electrolyte interface as the overlying oxide is etched. We denote the time delay between formation and etching of oxide species by τ. In other words, τ is the time required to etch the oxide layer once, and thus, τ is proportional to the oxide thickness.Figure 13Heatmap of the simulated spatiotemporal evolution of the SiO fraction in the oxide layerThe SiO fraction is normalized to the number density of Si and given in the right color bar. The upper right and lower left triangular areas are filled with electrolyte and p-Si bulk, respectively. uapp=5V.
The role that the time delay τ plays in the oscillation mechanism becomes evident when zooming in on the development of the oxide layer over slightly more than one oscillation period and analyzing the minor changes in etching rate and oxide-layer thickness in more detail. Figure 14 displays the section of Figure 13 between 330 and 414s in plate C together with the etching velocity of the oxide, the oxide-layer thickness, and the temporal evolution of average electric field strength across the oxide in plates A, B, and D, respectively. The oscillation mechanism can be deduced from the changes at the four points in time, which are labeled (i)–(iv).Figure 14Simulated time series(A) Etching speed, (B) oxide-layer thickness, (C) fraction of partially oxidized defects, and (D) averaged electric field intensity across the oxide. uapp=5V. The color bar gives the normalized SiO fraction in the oxide layer.
At time (i), the electric field just starts to increase, which causes the reaction rate of (R2) to increase more strongly than that of (R1). In other words, the relative concentration of SiO at the bottom layer of the oxide decreases, as evidenced by the transformation from red to orange at the bottom of the oxide layer in plate C. (ii) marks the time when this lowered SiO defect number density reaches the oxide/electrolyte interface. Hence, time (ii) is τ later than time (i). The lower defect number density causes the etch rate to decrease such that the thickness of the oxide increases, which, in turn, reduces the electric field intensity. At (iii), the electric field intensity is so low that the SiO number density increases again. 2τ later than (i), at (iv), this increased number density of SiO reaches the oxide/electrolyte interface and causes the etching velocity to increase again, resulting in a decrease of the oxide thickness and an increase of the electric field, and so the cycle starts anew. These results validate the conjectured mechanism summarized in Figure 1.
The oscillation mechanism entails that approximately two oxide layers are etched during one oscillation period. Direct measurements of τ are not available, but this model prediction agrees well with estimations of the relation between the oscillation period and τ from literature data, more precisely from measured oxide thicknesses, the average oxidation rate, which are equal to the average etching speed, and the period of the oscillations.^18^^,^^26^ Furthermore, the model predicts that oscillations occur only above a threshold voltage and that the period of the oscillation increases with increasing potential. Also, these two model predictions agree well with experimental results,^18^ see also Figure 10.
The above considerations revealed that the time delay τ between the change in the ratio of the SiO and SiO2 number densities at the Si/oxide interface and its effect on the etch rate resumes the role of a positive feedback and is essential for the occurrence of the oscillations. The relations between (i) and (iii) (E↗ → [SiO]b↘→vetch↘→L↗→E↘) as well as between (ii) and (iv) form negative feedback loops that damp perturbations. The strength of this negative feedback is controlled by the sensitivity of defect production rate to the electric field (Equation 6) and by the dependence of the etching rate on the defect concentration, parameterized by kSiO. If this negative feedback acted instantaneously, perturbations would be damped, and no sustained oscillations would occur. Only when the delay is sufficiently large, i.e., the oxide sufficiently thick, can oscillations occur. Thus, the time delay effectively weakens the immediate negative feedback and destabilizes the steady state. On the other hand, any specific reaction or transport steps within the oxide seem to be irrelevant for the oscillation mechanism. This suggests that the model can be considerably simplified, which we will do in the next section.
As discussed below, the oxide layer takes a passive role in the oscillation mechanism. Its only role is to provide a time delay between the formation and the etching of Si-oxide species. In this section, we show that indeed the electrochemical oscillation can be qualitatively described using a time-delay model with only two the oxide-layer thickness, l and the number density of substoichiometric oxide species, [SiO]. With this simplified model, we only aim to capture the qualitative behavior. All variables and parameters are dimensionless.
The time evolution of the oxide-layer thickness l is determined by the difference in the rates with which bulk Si is oxidized, i.e., the velocity with which the Si/oxide interface moves, vox, and the etching rate, that determines the velocity of the oxide/electrolyte interface, vetch.(Equation 29)dldt=vox−vetch.As above, we define the velocity of the Si/oxide interface by the reaction speed of (R1). When Si and SiO2 are present in large quantities at the Si/oxide interface, the reaction rate is determined by the number of valence-band holes at the interface, which depends on the electric field at the interface. For the sake of simplicity, we assume here that vox is proportional to the electric field, E(t), and that the electric field at the interface can be approximated by the average field across the oxide layer(Equation 30)vox(t)=γ1E(t),where γ1 is a constant. Furthermore, since the potential drop across the space-charge layer is negligibly small compared to the one across the oxide layer, we approximate E(t) by(Equation 31)E(t)=uappl(t).The etching rate is given by Equation 12.
In the simplified model, we express the number density of the partially oxidized defects at the oxide/electrolyte interface through its number density at the bottom layer. As seen in Simulations, SiO species hardly diffuse in the oxide. Rather, in the stationary (laboratory) frame, they practically remain at the z position where they were formed. Thus, they reach the top layer when the oxide layer above them is etched, i.e., τ time units later. When [SiO(t)] denotes the number density of SiO at the bottom layer at time t, the defect number density at the top layer is, thus, well approximated by [SiO(t−τ)], and, consequently, the number density of SiO2 at the top layer at time t is equal to 1−[SiO(t−τ)]. Furthermore, the results of the full model show that τ does not change much during an oscillation. In the simplified model, we approximate τ as a constant. The assumed changes of the oxide position and composition during a time interval of τ are sketched in Figure 15. Accordingly, the etching rate is expressed as(Equation 32)vetch(t)=γ2(1−[SiO(t−τ)]+kSiO[SiO(t−τ)]),where γ2 and kSiO are constants.Figure 15Illustration of the Si/oxide/electrolyte interface in the simplified modelDuring the time τ, the two interfaces change from the left state to the right state.
The number density of partially oxidized defects at the bottom layer, [SiO], results from the oxidation reactions (R1) and (R2). We consider a differential equation that approximately describes the time evolution of [SiO] at the Si/oxide interface. In the full model, the effects of reactions (R1) and (R2) were considered in the moving frame. In this moving coordinate system, the interface moves at a velocity vox=r1/ρatom. When we consider an infinitesimally thin volume at this moving interface, the influx and outflux of SiO due to reaction (R1) are equal, meaning that (R1) does not change [SiO]. On the other hand, the influence of (R2) must satisfy the second oscillation condition, namely that the number of defects decreases with increasing electric field. In simulations of the full model, the value of [SiO]|zb stabilizes near a constant value. Therefore, in this simplified model, we assume that [SiO] relaxes to a given value (the first term in Equation 33), and we define a sigmoid function (the second term) that decreases the defect fraction with increasing E, as does (R2) in the full model.(Equation 33)d[SiO]dt=γ4(γ5−[SiO(t)])−γ6[SiO(t)]1+exp(γ7(γ8−E)).
In this section, we present the simulations obtained using the simplified model and examine the oscillation mechanism in more detail through linear stability analysis.
Figure 16 shows simulation results obtained with the simplified model Equation 29 to Equation 33 under fixed potential conditions between uapp=0.5 and 1.75 with τ=50. The time series of the number density of the defects at the oxide/electrolyte interface and the oxide-layer thickness are plotted in plates A and B, respectively. Obviously, this simplified model also captures the oscillations. The dependence of the oscillation period and the average oxide thickness on the applied voltage is plotted in plates C and D, respectively. As in the full model, both the period of the oscillations and the oxide depth increase linearly with the potential. However, in the simplified model, the period increases by less than 20% while the thickness becomes nearly four times as thick. In contrast, the full model predicts an increase in oscillation period by a factor of 2 in a voltage interval in which the average thickness does not fully double (cf. Figure 9). This discrepancy arises from the fact that in the simulations, we used the same value of τ, whereas in the real system, the greater thickness at higher values of uapp also implies a larger value of the delay τ. In positive terms, the period is again about twice as long as the delay time, τ, namely (2.4±0.2)×τ, which is an inherent feature of the mechanism and, as already mentioned above, in good agreement with experimental data.^18^^,^^26^Figure 16Simulation results of the simplified model(A) Time series of the SiO number density, (B) time series of the oxide-layer thickness, (C) oscillation period, and (D) time-averaged oxide-layer thickness as a function of the applied voltage. The different colors in (A) and (B) represent different values of the applied voltage between 0.5 and 1.75. In (A), only three curves are shown to avoid the overlap of the data.
To support this interpretation, in Figure 17, we show simulation results in the same potential range as in Figure 16 but with adjusted values of τ, which were taken to be proportional to the average oxide depth as obtained in simulations with a constant value of τ (see Figure 16D). The selected values of τ in ascending order of potential are 50, 75, 100, 126, 151, and 177. Again, the defect number density and the oxide thickness time series are plotted in plates A and B, respectively. Here, the oscillation amplitude of the defect number density stays at a similar value in the potential range, and the oscillation amplitude of the oxide thickness increases with increasing potential. This relation is plotted in C.Figure 17Simulation results of the simplified model Equations 29 and 33 with adjusted ττ was adjusted to the thickness of the oxide layer as obtained with a fixed value of τ (see Figure 16) for each value of uapp. (A) Time series of the SiO number density, (B) time series of the oxide-layer thickness, and (C) time-averaged oxide-layer thickness as a function of the applied voltage. The different colors in (A) and (B) represent different applied potentials between 0.5 and 1.75. In ascending order of the potential, τ is set to 50, 75, 100, 126, 151, and 177. In (A), only three curves are shown to avoid the overlap of the data.
Having seen that the simplified model qualitatively captures the Si oscillations, in this section, we substantiate the interpretation of the oscillation mechanism, in particular, the essential role of the delay time, τ, by means of a linear stability analysis.
For simplicity, we rename the two variables as follows.(Equation 34)l(t)[SiO(t)]≡x(t)y(t)≡x(t)The governing equations are then given by(Equation 35)x(t)˙=γ1uappx(t)−γ2(1+γ3y(t−τ))(Equation 36)y(t)˙=γ4(γ5−y(t))−γ6y(t)1+expγ7γ8−uappx(t).At a fixed point, y(t) and y(t−τ) are equal. Therefore, it is straightforward to numerically solve for the fixed point, (x∗,y∗), determined by x(t)˙=0 and y(t)˙=0. We linearize the system around the fixed point and analyze the linear stability for small perturbations. With the following definitions of f(x,y) and g(x,y)(Equation 37)f(x,y)≡γ1uappx−γ2(1+γ3y)(Equation 38)g(x,y)≡γ4(γ5−y)−γ6y1+expγ7γ8−uappx,the partial derivatives are given by(Equation 39)∂f∂x=−γ1Ueffx2(Equation 40)∂f∂y=−γ2γ3(Equation 41)∂g∂x=γ6yγ7uappx2expγ7γ8−Ueffx1+expγ7γ8−uappx2(Equation 42)∂g∂y=−γ4−γ61+expγ7γ8−uappx.The linearized system can then be expressed as(Equation 43)x(t)=Ax(t)+Bx(t−τ),where A and B are defined by(Equation 44)A≡∂f∂x0∂g∂x∂g∂y(x∗,y∗),B≡0∂f∂y00(x∗,y∗).
If the oscillation mechanism discussed above is correct, the stability of the fixed point should depend on the delay time, τ, losing its stability for increasing τ. This can be verified by means of a Laplace transform. With the following definition of the functions p(s,±1)(Equation 45)p(s,±1)≡det(sI−(A±B)),the fixed point is stable if (a) the real parts of the roots of p(s+,1)=0 are both negative and (b) at least one of the real parts of the roots of p(s−,−1)=0 is positive.^27^ The numerical solution of Equation 45 shows that this system satisfies the delay-dependent stability conditions(Equation 46)Re(s+1)<0∧Re(s+2)<0(Equation 47)Re(s−1)>0∨Re(s−2)>0.For a certain value of τ, Re(s1,2)=0, and the roots si are purely imaginary. This means that a Hopf bifurcation occurs and a limit cycle is born. Correspondingly, the amplitude exhibits a square root dependence as a function of the voltage, as seen in Figure 9 D. The existence of a Hopf bifurcation in this system is also confirmed in the experiments.^28^
In summary, starting from ideas of the point defect model that describes the motion of point defects in passive films on metal surfaces, we derived a mathematical model for the formation and dissolution of an oxide layer during the anodic polarization of silicon in fluoride-containing electrolytes. Our minimal model incorporates just two defects, substoichiometric SiO and oxygen valencies, SiO2+, whose temporal evolution is described by mass balance equations coupled to Poisson’s equation in the oxide layer. Oscillations arise when three conditions are Firstly, the more substoichiometric SiO the oxide layer contains, the faster it is etched. Secondly, the fraction of substoichiometric oxide formed at the Si/oxide interface is higher, the lower the electric field at the interface is. In other words, a large electric field yields more fully oxidized SiO2. Thirdly, the time delay between the formation of oxide at the Si/oxide interface and its etching at the oxide/electrolyte interface must be sufficiently long. This third condition creates positive feedback, which, in conjunction with the negative feedback loops of the first two conditions, enables oscillations. The numerical results of the model agree well with experimentally observed data. In particular, the dependence of current, oxide-layer thickness, and period on the applied voltage is captured under potentiostatic and potentiodynamic conditions. Furthermore, the oxide layer is etched approximately twice during one oscillatory cycle in simulations and experiments.
The oscillation mechanism was validated using a simplified two-variable time-delay differential model that contained only the essential features of the oscillation mechanisms. The simplified model allowed us to perform a linear stability analysis that further confirmed the interpretation of the oscillation mechanism.
We believe that the model provides a solid starting point for further investigations in various directions. First, oscillations were observed during anodic polarization of many semiconductor and metal electrodes, most of which are not yet understood. While our current work focuses on Si, the oscillation mechanism discovered here is not related to specific properties of silicon electrochemistry. Rather, it seems to be likely that the three components of the mechanism are also present in other electrode materials and could therefore explain the occurrence of oscillations in various systems. The proposed mechanism suggests clear experimental strategies for testing whether the delay-induced mechanism causes oscillations in in other electrochemical systems. For example, systematic variation of film layer thickness or electrolyte composition would provide a direct means of probing the role of transport-induced delay in these systems. Within the framework of the point defect model, the dynamic formation and annihilation of point defects at the metal/oxide and oxide/electrolyte interfaces can introduce intrinsic delays between defect generation, migration, and oxide growth/dissolution. Such time delays may give rise to oscillatory behaviors analogous to those observed in anodic silicon electrodissolution. Second, the model appears to provide the basic properties of the electrochemistry of oxide-covered silicon electrodes. Extensions of the model seem to be straightforward and should uncover further fundamental properties of Si electrochemistry. This includes a better description of the composition of oxide,e.g., by taking into account electrolyte species that penetrate the oxide, the formation of pores at high electric field strengths, as considered in the current-burst model and most likely responsible for the emergence of HAOs, or charge-transfer reactions within the oxide, which seem to dominate the dynamics in the negative differential resistance region (II) in cyclic voltammograms (cf. Figure 11). Thirdly, extending the model to include spatial directions parallel to the electrode surface should provide information about properties of the spatial coupling of different locations of the electrode through valence-band holes, and their role in spatial pattern formation. It has been speculated that this coupling is responsible for the formation of unusual patterns such as chimera states or frequency clusters. Therefore, expanding this research will establish an important link between the numerous theoretical studies on these patterns and their experimental realization.
Experimental studies indicate that increasing the Ohmic resistance in series to the electrode promotes oscillations, and in particular stabilizes oscillations which are damped in the absence of the resistance. In our simulations, this promoting feature of an additional series resistance was not seen. The reason is not apparent. There seem to be different possibilities, e.g., that it results from the spatial extension of the system, or from charge-transfer reactions within the oxide that have been neglected here but, as described above, are necessary for capturing the negative differential resistance in the CV (region II in Figure 11). When extending the model, particular attention should be paid to this aspect.
Requests for further information and resources should be directed to and will be fulfilled by the lead contact, Katharina Krischer (krischer@tum.de).
This study did not generate new unique reagents.
•The raw experimental data supporting the findings of this study are available at Mendeley Data (https://doi.org/10.17632/dkd92jj7xk.1).•The code is accessible via GitHub (https://github.com/Yukiteru-san/Time-delay-as-the-origin-of-oscillations-in-anodic-Si-electrodissolution).•Any additional information required to reanalyze the data reported in this article is available from the lead contact upon request.
We thank M. Duportal for fruitful discussions. This project was funded by the 10.13039/501100001659Deutsche Forschungsgemeinschaft (DFG, 10.13039/501100001659German Research Foundation, project KR1189/18-2).
Y.M.: conceptualization, methodology, investigation, and writing – original draft; K.K.: conceptualization, writing – review and editing, funding acquisition, resources, and supervision.
The authors declare no competing interests.
During the preparation of the original draft, the authors used Grammarly in order to polish the English. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.
REAGENT or RESOURCESOURCEIDENTIFIERDeposited dataRaw experimental dataThis paperhttps://doi.org/10.17632/dkd92jj7xk.1Software and algorithmsSimulation code for full and simplified modelsThis paperhttps://github.com/Yukiteru-san/Time-delay-as-the-origin-of-oscillations-in-anodic-Si-electrodissolutionOtherAdditional information required to reanalyze the datalead contactAvailable upon request
The electrochemical experiment data shown in this paper were measured with a three-electrode cell. The working electrode was made of p-doped 1-10Ωcm2 silicon with a back contact by thermally evaporated aluminum on the back side and annealed at 250°C for 15 minutes. The reference electrode was a saturated Hg|Hg2SO4, and the counter electrode a platinum wire. The electrolyte comprised 60mMNH4F and 142mMH2SO4 at room temperature. The electrolyte was deaerated by bubbling with argon for 20 minutes. During the measurement, the electrolyte was stirred by a magnetic stirrer.
Numerical simulations were conducted using MATLAB R2024b. Time integration was performed using the explicit Euler method using a time step Δt=1×10−5. We verified numerical convergence by performing simulations with different discretization steps in time and space. All parameters used in the model are compiled below. Parameters that are experimentally inaccessible were chosen within physically reasonable ranges to ensure that the model behavior remains consistent with experimental characteristics.
The following parameters were used for the calculation presented in Simulations.
q=1.602×10−19[C]kBT=1.381×10−23×300[J]ϵ0=8.85×10−12[F/m]ϵsc=11.7[−]ϵox=3.8[−]ρ0=1×1021[1/m3]ρatom=2.66×1028[1/m3]lm=3.135×10−10[m]
LD=kBTϵ0ϵsc/ρ0q[m]kvox=1ρatom[m3]DSiO2+=μSiO2+kBT2q[m2/s]
L0=4×10−9[m]μSiO2+=4×10−17[m2/(Vs)]DSiO=1×10−20[m2/s]Csc=2.85×10−2[F/m2]kvetch=4.8×10−11/ρatom[m4/s]kSiO=6[−]k1=6×10−6/(2qρatom)[m2/(Vs)]α1=1/7×10−3[m/V]k2=5.3/(2qρatom)[1/(ms)]α2=66[1/V]α3=5/11[V]kleak=1×10−5×ρatom[1/m3]
n=25[−]dt=1×10−5[s]
γ1=0.5[−]γ2=0.2[−]γ3=5[−]γ4=0.2[−]γ5=0.6[−]γ6=10[−]γ7=50[−]γ8=0.45[−]
Figures in this paper are plotted using Matlab, Microsoft PowerPoint, and Origin.