Electro-Oxidation of p-Silicon in Fluoride-Containing Electrolyte: A Physical Model for the Regime of Negative Differential Resistance

TL;DR

This paper develops a physical model for the negative differential resistance observed during the electro-oxidation of p-silicon in fluoride electrolytes, aligning well with experimental data on current-voltage behavior and oxide layer properties.

Contribution

It introduces a novel physical model accounting for partial and full oxidation states of silicon and their influence on etch rate, explaining the negative differential resistance phenomenon.

Findings

Model reproduces negative slope in current-voltage characteristics

Simulations match measured oxide thickness and dissolution valence

Impedance spectra are consistent with experimental data

Abstract

When Si is anodically oxidized in a fluoride containing electrolyte, an oxide layer is grown. Simultaneously, the layer is etched by the fluoride containing electrolyte. The resulting stationary state exhibits a negative slope of the current-voltage characteristics in a certain range of applied voltage. We propose a physical model that reproduces this negative slope. In particular, our model assumes that the oxide layer consists of both partially and fully oxidized Si and that the etch rate depends on the effective degree of oxidation. Finally, we show that our simulations are in good agreement with measurements of the current-voltage characteristics, the oxide layer thickness, the dissolution valence, and the impedance spectra of the electrochemical system.

Figures

Equations

$$({\mathrm{H}}_2\mathrm{O}{)}_{\text{solution}}$$

$${\mathrm{Si}+2h}^{+}+{\mathrm{O}}^{2-}$$

$$\mathrm{SiO}+\alpha {\mathrm{h}}^{+}+{\mathrm{O}}^{2-}$$

$${\mathrm{Si}+4h}^{+}+2{\mathrm{O}}^{2-}$$

$${\mathrm{SiO}}_x+6\mathrm{HF}$$

$${\partial }_t{x}_{\mathrm{a}}(t)$$

$${\partial }_t{x}_{\mathrm{b}}(t)$$

$$r_1(t)$$

$$r_3(t)$$

$$r_{\mathrm{etch}}(t)$$

$${\partial }_t{n}_{\mathrm{ion}}(x,t)$$

$${\partial }_t{n}_{\mathrm{SiO}}(x,t)$$

$${\nabla }^2\phi (x,t)=-\frac{z_{\mathrm{ion}}F}{{\epsilon }_{\mathrm{ox}}}{n}_{\mathrm{ion}}(x,t)$$

$$J_{\mathrm{ion}}(x,t)$$

$$J_{\mathrm{SiO}}(x,t)$$

$$r_2(x,t)$$

$$J_{\mathrm{ion}}(x_{\mathrm{a}}(t),t)$$

$$J_{\mathrm{SiO}}(x_{\mathrm{a}}(t),t)$$

$$n_{\mathrm{ion}}(x_{\mathrm{b}}(t),t)$$

$$J_{\mathrm{SiO}}(x_{\mathrm{b}}(t),t)$$

$$\phi (x_{\mathrm{a}}(t),t)$$

$$\phi (x_{\mathrm{b}}(t),t)$$

$$I(t)=I_{\mathrm{reac}}(t)+I_{\mathrm{ch}}(t)$$

$$I_{\mathrm{reac}}(t)=2F\left[r_1(t)+2\cdot r_3(t)+\underset{x_{\mathrm{a}}(t)}{\overset{x_{\mathrm{b}}(t)}{\int }}\mathrm{d}x{r}_2(x,t)\right]$$

$$I_{\mathrm{ch}}(t)=-{\epsilon }_{\mathrm{ox}}{\partial }_t\left[\underset{h\searrow 0}{\lim }{\partial }_x\phi (x_{\mathrm{a}}(t)+h,t)\right]$$

$$\nu (t)=\frac{I_{\mathrm{reac}}(t)}{F{\partial }_t{x}_{\mathrm{a}}(t)n_{\mathrm{ox}}}.$$

$$U_{\omega }(t)=1\mathrm{mV}\cdot \sin (\omega t)+U_0$$

$$Z_{\omega }(t)=\frac{{\tilde{U}}_{\omega }(t)}{{\tilde{I}}_{\omega }(t)}$$

Full text

1 Introduction

When a silicon electrode is electrochemically dissolved in a fluoride containing electrolyte, it is covered by a layer of silica which can change its thickness in peculiar spatio-temporal patterns chazalviel1992ionic ; lewerenz1993origin ; lehmann1996origin ; foll2009self ; schonleber2012high ; miethe2009irregular ; miethe2012ellipsomicroscopic ; proost2014origin . (For a comprehensive overview of the work until 2003 see Chapter 5 in zhang2001electrochemistry .) These spatio-temporal patterns make silicon in hydrofluoric solution one of the most relevant electrochemical systems for basic studies of non-linear dynamics. The most prominent example is the chimera state in which one part of the ensemble oscillates in synchrony, while the rest is turbulent abrams2004chimera ; schmidt2014coexistence ; schonleber2014pattern ; schmidt2015chimeras ; patzauer2017autonomous . Although it is possible to find some attempt of modeling schonleber2016comparison ; zensen2014capacitance ; chazalviel1992theory , not even the simplest, periodic, spatially homogeneous oscillations have so far been consistently and convincingly ascribed to a concrete physical mechanism.

When trying to model the system one has to consider the following processes. Applying a postive voltage to a silicon electrode in an aqueous electrolyte results in the formation of silicon dioxide zhang2001electrochemistry . As the oxidation reaction proceeds, the resulting layer of oxide eventually passivates the electrode, inhibiting further reaction. However, if the solution contains HF or HF${}_{2}^{-}$ species, the oxide layer is etched away uhlir1956electrolytic ; osseo1996dissolution . The oxide layer thickness may, thus, attain a steady state at which the rates by which the oxide layer grows into the silicon electrode and is dissolved by the fluoride containing solution balance exactly. The current voltage characteristics of this equilibrium is very similar for a wide range of parameters chazalviel1991voltammetric and can exemplarily be seen in Figure 1 for voltages positive to the porous oxide formation. The peak marks the voltage above which an oxide layer is formed eddowes1990anodic ; blackwood1992electrochemical .

This equilibrium can become unstable at certain experimental parameter values which then leads to spatio-temporal variations of the oxide layer thickness.

In this article we follow the mindset that complicated dynamics bifurcate from simple dynamics, and hence that the first step towards a comprehensive understanding must be to find out which physical processes lead to the equilibrium properties, for example the stationary current voltage characteristics and the electrochemical impedance. So far there is a model for the current voltage characteristics in the electropolishing regime at voltages below the current peak in Figure 1 cheggou2009theoretical . There is also a more abstract model for the electrochemical impedance in the resonant voltage regime chazalviel1992theory . We present a physical model that explains both the current voltage characteristics and the impedance spectrum in the regime of negative differential resistance. Moreover, it reproduces the dependence of the oxide layer thickness and the dissolution valence on the applied voltage.

The central idea of our model is that the oxide is etched with a rate that depends on its composition. This assumption was motivated by measurements of dissolution valence eddowes1990anodic ; blackwood1992electrochemical ; schonleber2012high , measurements of the current voltage characteristics, and by the fact that the etching is likely to occur purely chemically, i.e. without charge transfer memming1966anodic ; eddowes1990anodic ; cattarin2000electrodissolution . Then, naïvely, the applied voltage should not influence the etching of the oxide layer. However, the rates of the purely chemical etching and the voltage dependent oxidation have to be the same at the steady state. Hence, there has to be an indirect dependence of the etch rate on the applied voltage. We suggest that this indirect dependence can be explained by the change in oxide composition as the voltage is varied. This change could arise due to the accumulation/depletion of ions or defects moving through the oxide layer. In this article we propose a physical model incorporating this idea.

First we explain our applied methods in Section 2. The model is described and defined in Section 3. In Section 4 the predictions of the model are compared to the experimental data which we obtained from cyclic voltammmetry, ellipsometry, and dynamic multi-frequency analysis as well as to values of the dissolution valence from the literature eddowes1990anodic ; blackwood1992electrochemical ; schonleber2012high . Finally, we conclude the paper in Section 5.

2 Methods

2.1 Experimental methods

All electrochemical measurements were performed with p-doped silicon electrodes ((111) orientation, 5-25 $\mathrm{\Omega}\,\text{cm}$). Ohmic back contact preparation, surface pretreatment and cleaning were done as described in patzauer2017autonomous . The aqueous electrolyte was prepared from ultra pure water (ELGA Purelab Ultra, 18.2 M$\mathrm{\Omega}\,$cm, TOC $<\,$3ppb), NH4F (Merck, p.a.) and H2SO4 (Merck, Suprapur). For all experiments shown here, the electrolyte contained 0.05 M NH4F and the pH was adjusted to 2.3, according to dissociation constants from the literature cattarin2000electrodissolution . The counter electrode was a Pt wire and a sat. Hg/Hg2SO4 electrode served as reference elecrode.

The in-situ ellipsometric measurements in Figure 4 and the cyclic voltammogram in Figure 1 and 3 were obtained with the spectroelectrochemical set-up described in miethe2012ellipsomicroscopic . In short, the stationary working electrode is illuminated with elliptically polarized, blue light from an LED with 470 nm wavelength through a window in the cell wall. After reflection at an incident angle of approximately 70 ∘ the polarization of the light has changed, depending on the optical path through the oxide. The change in polarization is converted to an intensity signal by letting the light pass through a polarizer and is measured using a CCD camera. The experiments were controlled with an FHI potentiostat (Fritz-Haber-Institut, Elab), and the electrolyte was stirred with a magnetic stirrer.

The dynamic impedance measurements battistel2016analysis ; koster2017dynamic in Figures 7b, 8b, 9 were done using a rotating disk electrode set up (AFMSRCE, Pine Research Instrument, Inc) with a disc diameter of 5 mm. A multisine perturbation was generated with a 2-channel waveform generator (Keysight) and added to the potential applied by the potentiostat (BioLogic Science Instruments, SP-300) to the working electrode. For the evaluation of the data we employed dynamic multi-frequency analysis as outlined in battistel2016analysis ; koster2017dynamic . Further experimental details of the impedance measurements and the cell used can be found in koster2018measurement .

2.2 Numerical methods

The numerical results were obtained with the proprietary finite element software COMSOL 5.2 comsol2016comsol . We solved partial differential equations numerically on a one dimensional spatial domain perpendicularly to the electrode, which represents the oxide layer of finite thickness. From the available COMSOL-physics-components we used the “General Form PDE” component for chemical species, the “Poisson’s Equation” component for the electrostatic potential, and “Moving Mesh” component to simulate the growth and dissolution of the oxide layer. The spatial domain was split into 2000 finite elements of the same length, which were adapted continuously as the boundaries moved. All simulations were run with the COMSOL-study “Time Dependent” and all solver configurations were left at their default except for the “Relative tolerance” which we reduced to $10^{-7}$ for the electric impedance calculations and to $10^{-5}$ in all other cases. All the parameters, variables and physical constants used in the simulations are listed in Tables 1, 2, 3.

3 Model description

In this Section we model the oxide layer, SiOx, that covers the Si anode in hydrofluoric solution. We aim to reproduce the negative differential resistance, i.e. $\mathrm{d}I/\mathrm{d}\varphi_{\text{WE}}<0$, where $I$ is the electric current and $\varphi_{\text{WE}}$ is the potential of the working electrode.

3.1 Concept

We assume the electrodissolution of Si to occur in three steps which are sketched in Figure 2. The initial step is silicon oxidation at the interface between the silicon and the oxide layer. In this step, the oxide layer grows into the silicon electrode. The pure silicon is, thus, replaced by a partially oxidized silicon species, e.g. SiO.

The second step happens inside the layer volume and transforms the partially oxidized silicon, SiO, into the final oxidation state, SiO2. The transformation from Si to SiO2 is likely to involve several sub-steps and further intermediate oxidation states, correspondingly. Our model merges these sub-steps for simplicity. In the third and final step the layer is etched away purely chemically. The etch rate is assumed to depend on the oxide composition and to be faster if the fraction of partial oxide is higher.

One possible way to realize the mechanism we stated above is given by the following set of reactions. Ions O2- enter the oxide layer by a reaction like:

[TABLE]

Under the present high electric fields (V/nm) the ions travel through the layer and oxidize the silicon in two steps

[TABLE]

or in just a single step

[TABLE]

Reaction (R1), occurs at rate $r_{1}(t)$ at the boundary between silicon and a mixed layer of SiO and SiO2. This mixture is referred to as SiOx in the following. Inside the mixed oxide layer the SiO created in Reaction (R1) can be further oxidized to SiO2 in Reaction (R2) at rate $r_{2}(x,t)$. Alternatively, Si can be directly oxidized to SiO2 in Reaction (R3) at rate $r_{3}(t)$.

The last step is the etching of oxide at the boundary between oxide and solution, named SiOx/solution interface, with a rate $r_{\mathrm{etch}}(t)$ depending on the oxide composition right at the boundary. The corresponding reaction reads

[TABLE]

with some $x$ between 1 and 2. The dissolution of SiO2 in hydrofluoric solution is well investigated osseo1996dissolution ; cattarin2000electrodissolution and more complicated than Reaction (R4), but our simplification will serve the purpose.

Although the reactions that we assumed above lead to good predictions, the mechanism could probably also be implemented with other reactions. One could for example consider OH- ions and Si(OH)x instead of O2- ions and SiOx.

3.2 Kinetic description

We transcribed the ideas from Section 3.1, which are sketched in Figure 2, to a set of coupled partial differential equations of time $t$ and of one spatial dimension $x$ which points perpendicularly away from the electrode surface. These partial differential equations describe the transport and the reactions inside the oxide layer and at its boundaries. As etching (R4), reaction (R1), and reaction (R3) proceed, the corresponding interfaces move. Let ${x_{\mathrm{a}}}(t)$ and ${x_{\mathrm{b}}}(t)$ be the time-dependent position of the Si/SiOx and SiOx/solution interfaces, respectively. Assuming a constant molar oxide density ${n_{\mathrm{ox}}}$ for simplicity, the velocities can be written in terms of the respective reaction rates:

[TABLE]

where the rates per unit area $r_{1}(t)$, $r_{3}(t)$ and $r_{\mathrm{etch}}(t)$ are defined below in Equation (3), (4) and (5).

Since we exclusively consider situations where the silicon surface is covered with some (though possibly thin) oxide layer, the rates $r_{1}(t)$ and $r_{3}(t)$ are assumed to be limited solely by a lack of O2- ions. Silicon and h+, the other reactants involved, are present in high concentrations because we consider p-type silicon under anodic bias. The etch rate $r_{\mathrm{etch}}(t)$ for a certain oxide mixture is linearly interpolated between the etch rate of pure SiO2 and the etch rate of pure SiO with proportionality coefficients $\eta_{\mathrm{SiO}_{2}}$, $\eta_{\mathrm{SiO}}$:

[TABLE]

where ${n_{\mathrm{ox}}}$ is the molar density of silicon in the oxide, and $n_{\mathrm{SiO}}(x,t)$, $n_{\mathrm{ion}}(x,t)$ are the molar densities of SiO and O2- respectively. The concentration of SiO2 is expressed as ${n_{\mathrm{ox}}}-n_{\mathrm{SiO}}(x,t)$. Hereby we assume SiO2 to be the non-SiO fraction of the entire molar oxide density ${n_{\mathrm{ox}}}$. The rate coefficients $k_{1}$ and $k_{3}$ depend on the potential drop across the space charge layer which is assumed to be constant in our model, because the differential capacity of the space charge layer is large under forward bias sze2006physics . Thus, $k_{1}$ and $k_{3}$ are constant as well.

Unlike the reactions (R1) and (R3) which are interface reactions, reaction (R2) occurs inside the oxide volume. From the participating species we expect electrons/holes to leave/enter the oxide quickly via direct tunneling or trap hopping which has been observed in the electroluminescence spectra of electrically biased oxide films hasegawa1988anodic . The remaining degrees of freedom in the oxide layer are then the concentrations of SiO and O2-. The dynamics of these two species is dictated by their respective mass balance equations which contain transport and reaction terms:

[TABLE]

where $J_{\mathrm{ion}}(x,t)$ and $J_{\mathrm{SiO}}(x,t)$ are the molar fluxes of O2- and SiO respectively.

Let $\varphi(x,t)$ be the electrostatic potential. Since O2- is the only charged species considered inside the oxide, the total charge density in the layer is $z_{\mathrm{ion}}Fn_{\mathrm{ion}}(x,t)$, where $z_{\mathrm{ion}}$ is the number elementary charges per ion, i.e. -2, and $F$ is the Faraday constant. Poisson’s equation thus reads:

[TABLE]

Introducing the diffusion constants $D_{\mathrm{ion}},D_{\mathrm{SiO}}$ and the rate coefficient $k_{2}$, the transport and reaction terms in Equation (6) and (7) can be modeled as:

[TABLE]

where diffusive transport was assumed to be Fickian in the Equations (9),(10) and the Einstein-Smoluchowski relation was used to derive the ionic migration term in Equation (9). $R$ and $T$ are the gas constant and temperature.

The diffusion of SiO reflects the thermal restructuring of the solid oxide which we expect to be relevant on the considered length scales of some nanometers, after trying various diffusion coefficients for SiO.

3.3 Boundary conditions

Continuity yields Equation (12) and (13) as boundary conditions for the Si/SiOx interface at ${x_{\mathrm{a}}}$:

[TABLE]

For the SiOx/solution interface at ${x_{\mathrm{b}}}(t)$ we assume a certain fixed concentration $n_{\mathrm{ion}}^{0}$ of O2- that is tantamount to assuming that Reaction (R0) is always in equilibrium. The oxide dissolution is accounted for by the movement of the SiOx/solution interface in Equation (2). Thus, the flux of SiO at the boundary should be set to zero:

[TABLE]

To come up with the boundary conditions for the electrostatic potential $\phi(x,t)$ we made the following simplifying assumption about the space charge layer. The differential capacity of the space charge layer of p-type silicon under strong anodic bias is very large, like in the case of a metal-insulator-semiconductor capacitor sze2006physics . The same is true for the Helmholtz layer which also has a large capacity. Therefore, any change in the total applied voltage will drop on the much smaller capacitor formed by the oxide layer. Thus, the oxide sees the externally applied voltage with an offset that is fixed. Under these assumptions we end up with the following boundary conditions for the electrostatic potential:

[TABLE]

where $U(t)$ equals the externally applied voltage plus/minus a constant offset.

3.4 Total electric current

The total electric current $I(t)$ that flows into the working electrode is the sum of the reaction current $I_{\mathrm{reac}}(t)$ and capacitive charging current $I_{\mathrm{ch}}(t)$.

[TABLE]

The current $I(t)$ is used in Section 4 to predict a cyclic voltammogram and an impedance spectrum and compare them to electrochemically measured data.

The reaction current $I_{\mathrm{reac}}(t)$ is calculated from the reaction rates $r_{1}(t)$, $r_{3}(t)$ and $r_{2}(x,t)$, where $r_{1}(t)$,$r_{3}(t)$ were defined as rates per surface and $r_{2}(x,t)$ as a rate per volume.

[TABLE]

Remember that we assumed that the electrons/holes created/consumed in the bulk at rate $r_{2}(x,t)$ leave/enter the oxide via trap hopping or direct tunneling, more or less immediately. Thus, the current resulting from $r_{2}(x,t)$ is contributing to $I_{\mathrm{reac}}(t)$ without any further delay.

The capacitive charging current $I_{\mathrm{ch}}(t)$ is calculated from the electric potential $\varphi(x,t)$ as

[TABLE]

which is simply the time derivative of the integrated space charge in the semiconductor. Equation (20) was obtained by integrating Poisson’s equation from $-\infty$ to ${x_{\mathrm{a}}}(t)+h$ and differentiating both sides of the equation by $t$.

3.5 Discussion on free parameters

There are quite some tunable parameters in the model. The electrochemical potential of ions in equilibrium and the transport and reaction rates were chosen somewhat arbitrarily, but for the remaining parameters one can make an educated guess. For the molar density $n_{\mathrm{ox}}$ and the permittivity $\varepsilon_{\mathrm{ox}}$ one can approximately take known values of thermal oxide. The etch rates of partial and stoichiometric oxide $\eta_{\mathrm{SiO}},\,\eta_{\mathrm{SiO_{2}}}$ can be roughly estimated from the experimentally observed etch rates. The offset between the model voltage parameter $U$ and the actual voltage across the interface can be estimated considering the dependence of the space charge layer’s differential capacity on the voltage drop across it.

4 Results and discussion

4.1 Current-voltage characteristics

The model is supposed to reproduce the negative slope in the current-voltage characteristics. In an electrochemical set up, this negative slope is measured between the broad current peak and the resonant current plateau in Figure 1, if the applied voltage is cycled. We simulated such cyclic voltammetry by making the model parameter $U$ a periodic triangular function of time, since $U$ represents the potential difference between the silicon electrode and the electrolyte. The parameter $U$ differs from the actual applied voltage by an offset due to the space charge layer and the Helmholtz layer, which are discussed in Section 3.3, but are not quantified in this model.

After an initial transient the resulting voltammogram is robust for a wide range of initial conditions. Figure 3 shows a simulated voltammogram with a scan rate of 1 mV s*-1* and the corresponding electrochemical measurement. In both curves the current decreases with higher electrode potential. With faster scan rates we observed that the hysteresis of both the simulated and the electrochemically measured cyclic voltammogram become larger and that the measured clockwise direction of the hysteresis is matched in the simulations. If the scan rate is low enough, the simulated cyclic voltammogram converges to a single line, corresponding to the stationary current voltage curve, again in agreement with experiments.

4.2 Oxide layer properties

The thickness of the oxide layer as obtained during the simulated cyclic voltammetry is plotted in Figure 4, together with the corresponding electrochemical measurement which shows an in-situ ellipsometric signal. The ellipsometric signal is approximately proportional to the oxide layer thickness but also depends on the refractory index inside the layer and in front of it. It can be used as a qualitative estimation of the layer thickness miethe2012ellipsomicroscopic . In both plots the thickness increases with the applied voltage.

The main idea of the model is that the etch rate of the oxide depends on its composition. Thus, let us have a look at Figure 5a which shows the stationary oxide composition for some voltages $U$. With increasing voltage $U$ the stationary layer thickness $x_{\text{b}}(t)-x_{\text{a}}(t)$ increases. Consider the four cases in Figure 5a. The concentration of partial oxide $n_{\mathrm{SiO}}(x,t)$ is the highest at the Si/SiOx interface $x_{\text{a}}(t)$, where the partial oxide is created. Closer to the SiOx/solution interface $x_{\text{b}}(t)$ the fraction of partial oxide becomes smaller as it is used up in Reaction (R2). The remaining fraction right at the SiOx/solution interface $x_{\text{b}}(t)$ determines the etch rate by Equation (5). Note that this fraction of partial oxide is much smaller at the SiOx/solution interface if the layer is thicker and $U$ is larger. This comes from an increased volume that is available for the ions $n_{\mathrm{ion}}(x,t)$ to participate in the bulk Reaction (R2), which consumes partial oxide. Consequently, if $U$ is larger, the layer is etched more slowly, according to Equation (5).

The concentration of ions $n_{\mathrm{ion}}(x,t)$ is the highest at the SiOx/solution interface $x_{\text{b}}(t)$ where the electrochemical potential of O2- in the oxide and the electrochemical potentials of protons and water in the electrolyte are assumed to be in equilibrium. With increasing distance from the SiOx/solution interface $x_{\text{b}}(t)$, which corresponds to going left in Figure 5a, the concentration of ions decreases. This is because the ions enter the oxide at the SiOx/solution interface $x_{\text{b}}(t)$ and are then used up in the oxide volume by Reaction (R2) and at the Si/SiOx interface $x_{\text{a}}(t)$ by Reaction (R1) and (R3). The concentration of ions in the layer determines the curvature of the electrostatic potential $\varphi(x,t)$ in Figure  5b by Poisson’s equation. Its derivative gives us the electric field $-\partial_{x}\varphi(x,t)$, which monotonically increases from about 0.2 V/nm at the SiOx/electrolyte interface to 2-3 V/nm at the Si/SiOx interface.

4.3 Dissolution valence

We calculated the reaction valence $\nu(t)$ in the presented model in terms of the reaction current $I_{\mathrm{reac}}(t)$, the Faraday constant $F$, the Si/SiOx interface position ${x_{\mathrm{a}}}(t)$ and the molar density of Si-atoms in the SiOx layer ${n_{\mathrm{ox}}}:$

[TABLE]

The calculated reaction valence during simulated cyclic voltammetry is shown in Figure 6. Due to the slow scan rate there is no noticeable hysteresis, if there is one at all. The simulated valence increases monotonically with the voltage from a value of about 3.4 close to the current peak to 4 at the plateau where it saturates. The same behavior has been confirmed by gravimetric measurements eddowes1990anodic and by rotating ring disc measurements blackwood1992electrochemical ; schonleber2012high The exact value of the valence in the limit $U\rightarrow 0$ can be adjusted by changing the rate constant of Reaction (R3).

4.4 Impedance spectra

To obtain the impedance spectra for our model described in Section 3 we set the parameter $U$ to a different function for each investigated angular frequency $\omega$, denoted by $U_{\omega}(t)$:

[TABLE]

The model equations were solved for many periods of $U_{\omega}(t)$ such that the resulting electric current $I_{\omega}(t)$ converged and became sine-shaped. We calculated the corresponding analytical signals $\widetilde{U}_{\omega}(t)$ and $\widetilde{I}_{\omega}(t)$ with the Matlab R2017b function hilbert(…) MATLAB:2017b . This means that $\widetilde{U}_{\omega}(t)$ and $\widetilde{I}_{\omega}(t)$ are obtained from ${U}_{\omega}(t)$ and ${I}_{\omega}(t)$ by removing the negative and zero frequency components. The impedance $Z_{\omega}(t)$ is then calculated as

[TABLE]

which would be constant in time if the current response was perfectly linear. $Z_{\omega}(t)$ is, of course, not completely constant in our non-linear small amplitude case. Thus, we determined the value for the impedance $Z_{\omega}$ at a certain angular frequency $\omega$ by time-averaging $Z_{\omega}(t)$ over some periods. Alternatively one could have calculated the impedance from the values of the Fourier transformed signals at $\omega$, which is equivalent for a linear response.

The impedance spectra obtained are plotted in Figure 7a and 8a and the corresponding electrochemically measured impedance spectra in Figure 7b and 8b. The plots in Figure 9 exemplarily compare one simulated impedance spectrum to one measured impedance spectrum. It can be seen that the simulated spectrum looks similar to the measured spectrum around 2 V vs SHE.

The model allows us to assign a physical meaning to some features of the measured impedance spectra which are illustrated in the equivalent circuit in Figure 10. For example, when omitting the charging term $I_{\mathrm{ch}}(t)$ in the numerical calculation (see Equation (20)), the small semi-circle at high frequencies in the Nyquist plot disappears. Thus, it is clear that the small semi-circle emerges due to interfacial charging, which means that opposing charges accumulate at the Si/SiOx interface and at the SiOx/solution interface. Hence, any equivalent circuit would require a capacitor $C_{1}$, see Figure 10.

At an intermediate frequency at which the small and the large semi-circle merge, the impedance is almost perfectly ohmic. This means that this frequency is sufficiently low for the layer capacity to be negligibly small but high enough so that the oxide thickness is not affected. Thus, the impedance at this frequency corresponds to the resistivity of the layer for a fixed thickness. This is represented in the equivalent circuit by a resistor $R_{1}$. Since the simulated impedance spectra at low frequencies resemble a semi-circle which ends at a real and negative value when $\omega\rightarrow 0$, we fit a capacitor $C_{2}\ll C_{1}$ in parallel with a negative resistance $R_{2}$. This approximation seems to be valid for the electrochemically measured impedance spectra in Figure 7b, at $U\approx 2$ V vs SHE.

There are some measurable features that the model does not reproduce. It can be seen in the Bode plot in Figure 8b that at high frequencies the measured spectra approach a phase that differs from $-\pi/2$, where $-\pi/2$ would correspond to an ideal capacitor. This is probably caused by an uncompensated resistance in series, which leads to a complex phase of zero at $\omega\rightarrow\infty$.

Another feature in the measured spectra, which is not reflected by the model, is the kink at frequencies between 0.1 Hz and 1 Hz. This kink becomes a loop that grows for higher voltage until at some voltage it possibly causes the resonance that is known to occur in the system ozanam1992resonant , by touching the origin in the Nyquist plot. If the impedance becomes zero at a certain frequency, it corresponds to a Hopf-bifurcation koper1994far , which could explain the emergence of oscillations that are observed in the system turner1958electropolishing ; schonleber2012high . However, the loop is not explained by the model and the scenario in the electrochemical set up is probably even more complicated, as can be guessed from earlier impedance measurements at other parameters ozanam1992resonant . Unfortunately, measuring the impedance of the steady state in the resonant regime appears to be challenging.

A third discrepant feature is measured at smaller anodic voltages. As the voltage is decreased in the electrochemical set up, one eventually passes the current peak and arrives at the electropolishing regime where the differential resistance is positive. Thus, the presumed contact point between the impedance curve and the real axis at zero frequency in the Nyquist plot moves from the left half plane to the right, as show in Figure 7b. The model, in contrast, is not valid for voltages this low, as there is no oxide layer of relevant thickness.

4.5 Origin of the negative differential resistance

From our model we can extract the necessary physical conditions leading to the negative differential resistance. First we note that in order to obtain the negative slope in the current voltage characteristics only the reactions (R0), (R1), (R2), and (R4) are required. Reaction (R3), in contrast, was only introduced to adjust the valence to above 3 at very small layer thicknesses, where volume Reaction (R2) hardly occurs.

Let us assume that the system is in a steady state on the negative differential branch. An increase in the voltage $U$ leads to an increase in the amount of O2- ions pulled into the oxide layer by reaction (R0) which increases the rate of the first oxidation step (R1), as the latter requires O2- ions. The resulting thickening of the oxide layer works against (R1) because the path that the O2- ions have to travel in order to reach the Si/SiOx interface becomes longer. Therefore, the thickening of the oxide layer partly compensates the rate increasing effect of the larger voltage $U$ on the first oxidation step (R1). The the second oxidation step (R2), however, which happens inside the layer volume, is supported by the thicker layer, or at least it is less suppressed compared to the first oxidation step (R1), which happens at the Si/SiOx interface. Thus, the oxide layer now contains a higher fraction of fully oxidized silicon, SiO2, as compared to a smaller applied voltage $U$. Consequently, the oxide layer is etched more slowly by Reaction (R4). This thickens the oxide layer even further such that the rate increasing effect of the larger voltage $U$ on the first step (R1) is overcompensated.

Note that without considering reaction (R2), any increase in voltage is compensated by a thicker oxide layer in such a way that the electric field inside the oxide layer as well as the current remain constant, independently of the ratios or the absolute vales of the rate constants of SiO (R1) and SiO2 (R3) formation.

In short, the negative differential resistance is the result of the interplay between the further oxidation of partially oxidized SiOx inside the oxide layer and the dependence of the etch rate on the stoichiometry of the oxide.

5 Conclusion

We presented a physical model that reproduces the cyclic voltammogram in the region of negative differential resistance for voltages below the resonant current plateau, the monotone increase in the valence, and the monotone increase in the oxide layer thickness. Moreover, our simulations reproduce our measured impedance spectra, which were obtained via dynamic multi-frequency analysis.

The model is build on principles of non-equilibrium thermodynamics, solid state chemistry and semi-conductors physics. It is constructed in a way to capture the behavior qualitatively and semi-quantitatively with minimal complexity. To this end, we could reproduce all measured trends assuming that the reaction is completely limited by migration of ions in the oxide layer. Hence, potential dependences of the electrochemical reaction rates or changes in the potential drops across the space charge layer in the Si or the Helmholtz layer could be neglected. On the other hand, a dependence of the chemical etch rate of SiOx on the oxide stoichiometry and on the (slow) further oxidation of SiOx within the oxide layer are crucial elements for reproducing the negative differential branch in the current-voltage characteristic.

This model presents a basis on which one can build to further uncover nonlinear features of the Si electrodissolution dynamics, such as current oscillations or the patterns observed with n-type Si electrodes. For this it is likely that some of the simplifying assumptions entering our base model have to be relaxed.