Electrochemical Oscillation Model in Electrohydrodynamic Fluid

In electrohydrodynamic phenomena, rhythmic oscillations of current under ac electric fields were observed in dielectric liquid, C4F9OCH3 (Abe et al., Jpn. J. Appl. Phys. (2010)). The oscillations of current and static pressure depended extensively on surface conditions of the electrodes, where an electric double layer on the surface of electrodes is sensitive to its smoothness. The characteristic frequency of current was defined by a crossover point on the power spectrum. This is based on non-synchronization with external ac fields above the critical frequency of current. Simulations are carried out using feedback loop including two-dimensional diffusion process. The observed crossover point on the power spectrum is reproducible qualitatively in the simulations.


Introduction
Synergetic phenomena are seen in physics, chemistry and biology [1,2].Generally, the self-organized oscillations are realized by a balance between activation and suppression in an unstable system.Much attention has been devoted to analyze the above oscillations using nonlinear dissipative dynamics as asymptotic methods.For instance, a simple two-state model, which is simply described by ordinary differential equations, is proposed to explain bacterial behaviors [3].In fluctuating environments, heterogeneous bacterial populations are influenced by the switching rate in the equations.Also, cell population heterogeneity in biochemical oscillating components is represented by feedback loops [4].Mathematical model is formalized to express activation delay in cellular control system.A switching mechanism between two stable steady states corresponds to bistability in response curve, which associates with a biological property such as bacterial virulence.
As a universal phenomenon, electrohydrodynamic (EHD) convection occurs under applied voltage.The electric body force acting on molecules is given by [5], q is charge density.ε and ε 0 are the dielectric constant of the liquid and vacuum, respectively.The first term in Equation ( 1) reveals ion-drag pumping based on Coulomb force.The second and the third terms are polariza-tion pumping as dielectrophoretic force, f DEP .In particular, the third term is electrostriction force, which is induced by non-uniform electric fields.Mainly, induction pumping, ion-drag pumping, polarization pumping, conduction pumping and electro-osmotic pumping are considered as EHD pumping mechanism.Serious problem of the ion-drag pumping is that dielectric liquid is easily damaged at ultimate high voltage because of including a charge injection mechanism [6].Accompanying a strong increment of electric current, pumping becomes unstable and the pressure is not reproducible.For an industrial application, it is pointed out that low durability in EHD module is not sufficient.Recently, the EHD conduction pumping is investigated experimentally and theoretically [7][8][9].The conduction pumping is based on the heterocharge layers near the electrodes.Under extremely high voltage, an electric double layer (EDL) in the vicinity of surface of electrodes is formed in the dielectric liquid, although dielectric liquid as non-electrolyte has no EDL under conventional voltage.In electrochemistry, an EDL consists of the Helmholtz layer and diffuse layer.In the dielectric liquid C 4 F 9 OCH 3 , self-organized oscillations of electric current in EHD are induced by ac electric fields [10].The rhythmic current oscillations depend on surface conditions of electrodes.In the low pressure region, an EDL contributes to pumping mechanism described by the second term in Equation ( 2), since formation of the EDL is sensitive to the surface conditions.Furthermore, under ac electric fields, two kinds of relaxation frequency (f 1 and f 2 ) of current were observed separately.Particularly at around the first relaxation frequency (f 1 ), which is smaller than f 2 , the oscillating components of current are divided into two regions.One is the synchronized region with the external ac frequency (f < f 1 ).The other is the non-synchronization region above f 1 .
In this study, we simulate rhythmic oscillations using the idea of feedback loops between the thickness of EDL and polarization.A combination of simple differential equations having an activator and an inhibitor can reproduce the observed cyclic current.By extending to twodimensional (2D) diffusion, the observed synchronization/non-synchronization between the external fields as an input and current as a response is simulated in spite of the complicated system.

Experimental
Dielectric liquid is C 4 F 9 OCH 3 (Sumitomo 3M Co., HFE-7100).The high-voltage power source is model APL-20K05PBX (Max-Elec Co.).ac frequencies (0.1 -30 Hz) are generated by a functional generator (Iwatsu Test Instruments Co, SG-4115).Model 6485 (Keithley Instruments) was utilized for monitoring electric current.For static pressure measurements, the samples and the electrodes were assembled inside a dry box in order to reduce atmospheric H 2 O and O 2 .Inside the dry box, helium gas was flowing to reduce humidity below 10% RH [10].
13 C-NMR experiments were carried out to investigate high voltage damage effect to molecules.500 MHz NMR spectrometer (Bruker Co., DMX500) was used for the measurements.Samples were equilibrated with additive of methyl-d3 alcohol-d (Acros Organics Co.) and dissolved in this solvent for analysis.The 13 C-NMR resonance was referenced to the signal (δ = 49.0 ppm) of the additive.

Experimental Results
Figure 1 shows NMR spectra before and after applying high voltage, respectively.Conditions of damaging test were 15.8 kV and 20 μA.The dc high voltage was continuously applied to the dielectric liquid for 2 hours.In spite of high voltage, there is no difference between applied and non-applied liquids.
Under ac electric fields, the power spectrum, G(ν), was calculated in order to analyze the oscillating components of the observed current [10].ν is the frequency of current oscillation and G(ν) is the Fourier transform of autocorrelation function.The inherent ν is easily determined on the G(ν). Figure 2 shows the inherent ν as a function of external ac frequency, f. ν has a linearity below f 1 with satisfying ν = f.This means that the oscillating current synchronized with f.Whereas, non-synchronization with  f appeared above f 1 .

Simulations
In dielectric liquid under high voltage, an EDL and the heterocharge layers are formed as shown in Figure 3 [7][8][9].The EDL consists of the compact layer (the Helmholtz layer) and the diffuse layer.Thickness of the compact layer is around 0.5 -1 nm.On the surface, a few molecular layers are fixed with some orientational order on the surface.The diffuse layer locates between the compact layer and the bulk liquid.The thickness of the diffuse layer depends on a property of dielectric liquids.Under high voltage, the heterocharge layers distribute additionally in outer part of the diffuse layer [7][8][9].A non-equilibrium state between the ionic dissociation and recombination is induced under high voltage.By the dissociated ions, the heterocharge layers, whose thickness is a few mm, are formed near the electrodes.Since Coulomb interaction as attractive force are affected on the heterocharge layers, pressure difference occurs between the electrodes and the heterocharge layers.This is a picture of the conduction pumping in EHD.
In the simulations, we assume that the relaxations of the relative hard "diffuse layer" and the soft "heterocharge layers" correspond to f 1 and f 2 , respectively.The assumption satisfies that surface roughness influences directly on the static pressure.Low-pressure regime at around f 1 is explained by the surface sensitive diffuse layer, while high pressure one above f 2 might be connected with relaxation of the heterocharge layers.Here, in the simulations, we focus only on the low-pressure region, since the high-pressure regime at around f 2 becomes quite complicated in the current oscillations.By applying the feedback loops to the EDL model, relaxation process of the diffuse layer is simulated systematically.

One Dimensional Activator-Inhibitor Model
Dissipative structure is demonstrated using feedback theory [4].Time lag in system is obtained by a combination of positive and negative feedback.The simplest feedback is given by, d d where α is constant.If α is positive, a positive feedback loop affects in the system.Schematic positive and negative feedback loops are illustrated in Figures 4(a) and (b).Generally, the negative feedback loops as inhibitor stabilize the system.For instance, numerical calculations converge to the constant value at any initial values.In contrast, the positive feedback loop (α > 0) enables it to diverge in the system.Furthermore, self-organization is described by a competition between an activator, α, and an inhibitor, h.The differential equation in one-dimensional (1D) model [1] is given by, D a and D h are diffusion constants.k 1 , k 2 , k 3 and k 4 are adequate constants, while ϕ is a variable parameter.We apply the idea to the observed current oscillations in EHD.Here, we focus on the relaxation of the diffuse layer denoted by f 1 .The equations without considering diffusion process are simplified by, δ is thickness of the diffuse layer as an activator, while P is polarization of the diffuse layer as an inhibitor.ϕ(t) reveals applied voltage.The feedback loops in Equations ( 5) and ( 6) are illustrated schematically in Figure 4(c).
The simplified model is based on the following assumptions: 1) the diffuse layer develops with increasing applied voltage.k 1 is responsible for promoting factor of the diffuse layer; 2) Oppositely, large enough δ is compressed by the outer heterocharge layers.Then, k 2 means an influence of the heterocharge layers; 3) In Equation ( 6), proportionally to δ, P increases.Self-suppressing term as an inhibitor is indispensable to prevent a divergence of P, since P is regarded as conservative quantity.
In Equation ( 5), if molecules in the diffuse layer become well ordered, the orientated molecules belong to the compact layer.P reduces the thickness of the diffuse layer.Then, P affects on the first term in Equation ( 5) as an inhibitor.
At the first step, the current oscillations are demonstrated under dc electric fields.Numerically, simulations are carried out by the finite difference method.At k 2 = 5, k 3 = 1 and ϕ(t) = 5, δ and P are simulated at each k 1 values.At k 1 = 2.5, both δ and P is dumping with a lapse of time (Figure 5(a)).At the condition, non-EDL state appears.Surprisingly, at k 1 = 3.5, limit cycle oscillation (LCO) [2] is induced as shown in Figure 5(b).The periodic oscillation in time shows the closed attractor in the δ-P plane.In fact, only using smooth surface of electrodes, rhythmic current was observed even under dc electric fields.Furthermore, at k 1 = 5, the EDL is formed steadily (Figure 5(c)).As a transition process from the non-EDL to EDL state, LCO occurs under dc electric fields.In Figures 5(a)-(c), a common feature of a relation between δ and P is that P varies gradually after development of the diffuse layer.Here, k 1 and k 2 are regarded as promoting and suppressing coefficients, respectively.In order to clarify the k 1 and k 2 dependence of LCO, three states are plotted in the k 1 -k 2 plane (Figure 6).The trivial balance between k 1 and k 2 governs the current oscillations.
Next, simulations under ac fields are carried out using ϕ(t) = ϕ 0 sin(2 πft).The current density, J, is connected  with polarization in electromagnetics, whose notation is A serious problem of the above model without diffusion is that LCO turns to be a non-EDL state at the high frequency region.Hence, the non-synchronization behavior such as LCO is not seen in the simulations at any parameters.

One-and Two-Dimensional Diffusion Model
As another approach to the non-synchronized behavior at f 1 < f, we add realistic diffusion process to the EDL model with the following assumptions: 1) In a flat part of surface, an EDL develops on the smooth surface (Figure 7); 2) In a rough part, the EDL is not generated on its own; 3) But, in the rough surface, the EDL is supplied by diffusion from the EDL on the flat one.At the same time, smoothness factor of the surface (ξ) is defined to be a ratio of flat part in the simulation size.
The one-dimensional (1D) model is simulated on 100 pixels along x direction (Figure 7).The EDL part is centered in the 1D lattice.For instance, at ξ = 10%, flat part is assigned in 45 55 i   90 i .Periodic boundary conditions are employed in the simulation box.The simulating conditions are determined to be k 1 = 3.5, k 2 = 5, k 3 = 1 and ϕ 0 = 5 with satisfying the LCO state under dc fields.Diffusion constant, D, is selected to be D = 5.By averaging each pixel along the x direction, total current densities are calculated.Figure 8(a) shows ac frequency dependence of current oscillations at fixing the smoothness factor (ξ = 80%, 10   ).Below 0.2 Hz, ν increases with increasing f.However, the simulated ν-f relation does not coincide with the synchronization line described as ν = f.On the other hand, above 0.2 Hz, constant ν  region appears.The 1D diffusion model, however, is insufficient to represent the fully relaxed diffuse layer.
To improve the 1D model, we extend it to two-dimensional (2D) diffusion process.In the 2D model, 60 × 60 pixels are prepared in the x -y plane of the simulation box, where periodic boundary condition is also employed.A part of an EDL is located in the center of 2D square lattice.On the 2D model, we must consider a shape of the EDL additionally.Square and circular shapes of the EDL are verified in the simulations.Consequently, it is found that the shapes are not so significant for the current oscillations.Then, we performed systematic simulations using the circle EDL. Figure 8(b) reveals the current oscillations as a function of ac frequency using the same parameters of the 1D model (ξ = 78.5%).At low frequency region, the complete accordance with the synchronization line (ν = f) at 0.4 < f is obtained on the 2D model.Furthermore, relative small ν is the seen at nonsynchronization region.This supports that feedback loops including diffusion can simulate the relaxed diffuse layer.It should be notices that, at 0.35 < f < 4.0 Hz (green closed circles in Figure 8(b)), the simulated oscillations are divided into two components, that is, ν H and ν L , as a transition state from synchronized to nonsynchronized regions.Apart from coexistence of rhythmic currents, the simulation results are very similar to the observed one.The 2D model is found to be more realistic in EHD.The above features such as synchronization/ non-synchronization is well simulated at around ξ = 80%.
As another approach to verify the 1D and 2D models, current oscillations in real time are plotted in Figures 9(a)-(c) compared with the observed one (f = 0.3 Hz).Not only time periodicity but also current profile in the 2D model coincides with the observed one in real time (Figure 9(a)).In addition, the weak modulation of the current profile is well described in the 2D model.On a contrary, the 1D model is insufficient both for time periodicity and the current profile.

Discussion
In NMR experiments, it was found that dielectric molecules are not damaged even at 15.8 kV.Typically, high static pressure and low current regime means that the conduction pumping without charge injection can prevent not only dielectric breakdown but also a molecular damage.At the same time, it is predicted that slightly dissociated ions under high voltage might recover completely to normal dielectric molecules when the applied voltage is set back to 0 kV.
The activator-inhibitor simulations are demonstrated for rhythmic oscillations of currents in EHD.The model is constructed considering dynamics of the diffuse layer.Consequently, we found that additional 2D diffusion process is indispensable to explain the experimentally obtained time periodicity and profile of current at low frequency region.In the 2D diffusion model, 2D periodic arrays of flat parts have a geometrical adjustment on diffusion process.For instance, the first nearest neighbor distance between the flat parts is different from the second nearest one.A variety of a pair distance between flat parts can correspond to wide range of the synchronization.Therefore, the current using the 1D diffusion model cannot synchronize with the external electric field due to a lack of degree of freedom.
Even in the 2D model, discrepancy between the experiment and simulation is seen at a transition state from synchronized to non-synchronized regions: In the ex-periment, a gradual decay of ν was observed at the transient frequency region, while coexistence of rhythmic currents represented by ν L and ν H is simulated.We deduce that actual flat part on the surface of electrodes is not arranged periodically.The other reason is that the flat part has a size distribution on the surface.Thus, diffusion process at some EDL is not so correlated with others.We interpret that the simulated two components are obscured in real EHD pumping.

Summary
Based on the experimental results, the EDL model including 2D diffusion process is proposed to express the rhythmic oscillations of currents.At low frequency region, the synchronization with ac fields in the simulations is in good agreement with the experimental result.In addition, at high frequency region, the diffuse layer is relaxed because of inherent size of the diffuse layer.The crossover point suggests that synchronized current oscillations change to non-synchronized one.The relaxation process is provided by the simple feedback loops in an open system.
For better understanding the EHD pumping mechanism, the dynamic property of the EDL is necessary in addition to the conventional fluid dynamics.We confirm that stable formation of the EDL including molecular conformations is a key to realize a high efficient EHD module.

Figure 2 .
Figure 2. Current frequency, ν, under ac electric fields, f.The frequency of oscillating components is calculated by Fourier transform of the observed current.ν changes at f 1 = 0.5 Hz drastically.

Figure 5 .
Figure 5.Time dependence of thickness of the diffuse layer and its polarization on the one-dimensional model under dc electric fields.Solid and dotted lines reveal thickness of the diffuse layer (δ) and polarization (P), respectively.Promoting constant, k 1 , can determine three different states.(a) Non-electric double layer (EDL) state (k 1 = 2.5); (b) limit cycle oscillation (LCO) (k 1 = 3.5) and (c) EDL state (k 1 = 5) are formed in the simulations.

Figure 6 .
Figure 6.k 1 -k 2 plot.k 1 and k 2 reveal promoting and suppressing of an electric double layer (EDL), respectively.Limit cycle oscillations (LCO) appears as a transient process from non-EDL to EDL state.

Figure 7 .
Figure 7.The extended electric double layer (EDL) model including one-dimensional diffusion process.In a rough part of the surface of electrodes, an EDL is not generated on the rough surface.By diffusion process from the flat part, the EDL is supplied.

Figure 8 .
Figure 8. Simulated current frequency, ν, as a function of external ac frequency, f, using (a) one-dimensional diffusion model and (b) two-dimensional one.

Figure 9 .
Figure 9.Time dependences of (a) the observed current; (b) the simulated current using the 2D model and (c) the simulated current using the 1D model.