Self-organized Multi-Frequency Clusters in an Oscillating Electrochemical System with Strong Nonlinear Coupling

We study the spatio-temporal dynamics of the oscillatory photo-electrodissolution of n-type Si in a fluoride-containing electrolyte under anodic potentials using in-situ ellipsometric imaging. When lowering the illumination intensity step-wise, we successively observe uniform oscillations, modulated amplitude clusters, and the coexistence of multi-frequency clusters i.e., regions with different frequencies, with a stationary domain. We argue that the multi-frequency clusters emerge due to an adaptive, nonlinear, and nonlocal coupling, similar to those found in the context of neural dynamics.

In this Letter, we report the emergence of selforganized multi-frequency clusters from a uniform oscillatory state during the photo-electrodissolution of an n-Si wafer when reducing the illumination intensity. Through the illumination, valence-band holes are created. Their movement parallel to the surface constitutes a nonlocal spatial coupling [23]. In addition, an external resistance in series with the electrode acts as a global synchronizing force on the dynamics [24]. Hence, there are two dominant types of coupling, a global synchronizing coupling and a long-range coupling through diffusion and migration of valence-band holes. Below we will argue that their interaction creates a nonlinear coupling that, in turn, promotes the formation of multi-frequency clusters.
The oscillatory photo-electrodissolution of Si in fluoride-containing electrolytes involves the electrochemical oxidation of Si to SiO 2 according to (1) and the chemical etching of SiO 2 via where 1 ≤ λ VB ≤ 4 is the amount of valence-band holes h + [25]. The experiments were conducted with an n-doped (1-10 Ωcm) Si (111) sample as the working electrode in a three-electrode setup. The electrolyte was an aqueous solution containing 0.06 M NH 4 F and 142 mM H 2 SO 4 . The uniformity of the electrode surface was monitored in-situ with an ellipsometric imaging setup which probes the change in optical path length at the electrochemical interface. The changes in optical path length are converted into an intensity signal ξ (x, t) (x denoting space and t time) and recorded with a CCD-camera (JAI CV-A50). We will present the ellipsometric intensity ξ (x, t) as a percentage of the saturation threshold of the CCD camera. To allow for the oxidation of n-type Si, the electrode was illuminated with a linearly polarized He-Ne laser (HNL150L-EC, Thorlabs), the intensity of which was adjusted with a linear polarization filter. Further experimental details can be found in in the appendix and in Ref. [23].
In the measurement presented below, the illumination intensity is the bifurcation parameter. We initialized the electrode by applying a constant voltage at a high illumination intensity and then decreased the illumination step by step. At each step, we waited until transients had died out and then measured the dynamics for 10 3 s. In order to characterize the dynamics of our system, we define the amplitude A (x, t) and phase φ (x, t) of the ellipsometric intensity signal ξ (x, t) at each pixel by calculating the analytic signal ζ (x, t) via the Hilbert transform H (ξ (x, t)) (for details see Ref. [26]): Having determined the time-series of the amplitude and of the phase, we extracted the dominant frequency ν (x) at each point from a linear fit to φ vs. t. Exemplary states from a measurement series can be seen in Fig. 1 where the temporally averaged amplitude A (x, t), the dominant frequency ν (x), and a snapshot of the phase at an arbitrary instant in time, φ (x, t = 807 s), are shown in the first, second, and third row respectively. The four columns depict measurements at four different illumination intensities. The initial, highly illuminated state is shown in column I. Here, the system oscillates uniformly with the same amplitude, frequency, and phase at each point in space, cf. Ref. [9,27].
Upon lowering the illumination intensity ( Fig. 1 column II), the electrode splits into a region with higher and a region with lower amplitude. These two regions still oscillate with the same average frequency, but the oscillation phase differs between points in the higher-and lower-amplitude regions. In other words, amplitude clusters have formed.
However, the data shown in Fig. 1 column II do not give the full picture of the dynamics. This can be seen if we look at Fig. 2. Here, the temporal evolution of the phase at an exemplary point, marked by a cross in the phase plots in Fig. 1 column I and II, is depicted in a frame rotating uniformly with the dominant frequency of the point in question. Starting with the higher illumination ( Fig. 2(a)) we observe only a simple modulation with the same dominant frequency as the one of the rotating frame. This is in fact the second harmonic of the dominant frequency of the original time series φ (t) and thus stems from its slight relaxational character. In the case with the lower illumination ( Fig. 2(b)), when the amplitude clusters have formed, we observe a further slow modulation of the phase evolution. This suggests that the system not only underwent a pitchfork bifurcation leading to amplitude clusters but also a secondary Hopf bifurcation creating the modulated oscillations.
When lowering the illumination further, two drastic changes are observed ( Fig. 1 column III). First, the mean amplitude differentiates further in space, suppressing the oscillations nearly completely on a part of the electrode. In this region, the very small amplitude combined with experimental noise leads to apparent discontinuities in the phase, rendering the determination of the dominant frequency impossible. Therefore, in the second and the third row of Fig. 1, we depict points with A (x, t) < 0.35%, in grey. Second, and perhaps even more astonishing, focusing our attention on the region that exhibits well defined oscillations, A (x, t) > 0.35%, we observe that the dominant frequency is not uniform anymore. Rather, the frequencies appear to accumulate around three plateau values, as apparent from the turquoise, red, and yellow patches in Fig. 1 column III, whereby the higher frequencies are found in the regions with higher mean amplitude.
In the last state ( Fig. 1 column IV) the features that appeared in column III become more pronounced; on a part of the electrode the amplitude is practically completely suppressed. In other words, on this part of the electrode we observe amplitude death [28]. Likewise, the frequency differences across the oscillating part of the electrode become more pronounced. Equal, or at least very similar frequencies now appear in connected regions, whereby the frequency distributions of the two outer orange and blue regions are very narrow, and the frequency distribution of the middle, 'mediating region' is somewhat broader, ranging from light-blue to yellow. Indeed, we witness the self-organized formation of multi-frequency clusters in a homogeneous oscillatory medium. Considering the snapshot of the phase distribution, we observe that the faster region on the right oscillates nearly uniformly whereas the more slowly oscillating region on the left exhibits a continuous distribution of the phases over 2π rad. This travelling-wave-type feature can be seen as the continuum version of a splay state in networks of coupled oscillators. Interestingly, the existence of mixed-type multi-frequency clusters consisting of a splay-type cluster and a phase-synchronized cluster, as we observe it here, has also been found in simulations of networks of phase oscillators with adaptive coupling, yet with the difference that the simulated phase-synchronized clusters occurred in antipodal pairs [16][17][18].
A key to understanding the changes in the dynamics is to realize that our bifurcation parameter controls the effective number of degrees of freedom in the system. At high illumination intensity µ there are more than sufficient valence-band holes for the oxidation process to take place equally everywhere on the electrode surface. Hence, the concentration of holes n h is effectively constant, and does not impact the uniform oscillation. The oscillations are synchronized by a global coupling arising from the presence of the external resistor and the potentiostatic control: Here, ϕ SC (x, t) and ϕ OX (x, t) are the potential drops across the space charge layer of Si, and the SiO 2 oxide layer respectively, U is the applied voltage, R is the external resistance, A is the electrode area, and i is the local current. The last term in Eq. (4) describes the potential drop across the external resistor which depends on the total current. Since at high illuminations ϕ SC (x, t) remains constant, the oscillating total current causes oscillations in ϕ OX (x, t), which in turn influence the reaction rate and thus the oscillations. Hence, our electrochemical oscillator creates a mean field I = x ∈A i (t) dx that feeds back into the dynamics of the oscillations, However, as we lower the illumination intensity µ, n h becomes so low that, at some point, it starts to limit the reaction current. To compensate for the lower reaction rate, ϕ SC (x, t) increases and thus becomes time dependent as well. Hence, we face a situation where a physical quantity, namely n h , starts to change in time when the value of a parameter crosses a threshold: The dynamics of n h now depend on the oscillating mean field I and also feeds back to said mean field as well as to the dynamics of the original 'base' oscillator. Our oscillating medium is thus nonlinearly coupled as soon as n h becomes a degree of freedom of the dynamics. If we consider our spatially continuous system as being composed of infinitesimally small base oscillators w k , one realizes that the nonlinear coupling is of the same type as the general physical setting for nonlinearly coupled oscillators formulated by Rosenblum and Pikovsky [5]: Here, w k forms a base oscillator, the ensemble of all oscillators produces some mean fields g, and v is a coupling variable that modulates the global coupling in a nonlinear way; µ is a bifurcation parameter. In our system, µ is the illumination intensity, g is the total current I, and v can be identified with n h . Nonetheless, in contrast to the global feedback variable v in Eq. (8), n h is not strictly global, but rather a nonlocal variable which depends on space [23]. Since n h influences the coupling strength of the base oscillator, cf. Eq. (5), the coupling becomes not only nonlinear but also space dependent. In this view, it is similar to the adaptive coupling discussed in Ref. [16][17][18], where multifrequency clusters were found. We thus attribute the occurrence of our multi-frequency clusters to the combination of the nonlinear and nonlocal coupling, which allows for a self-organized adaptation of the coupling strength: At parameter values at which multi-frequency clusters are found, the intra-cluster coupling strengths as well as all mutual inter-cluster coupling strengths differ.
Adopting a different perspective, one realizes that our dynamics also contain features that have been discussed in connection with certain types of chimera states. Provata considered a birhythmic model [29]. When coupling these oscillators nonlocally, synchronized domains oscillating in either of the two bistable limit cycles could be stabilized. The interfacial regions mediating between the domains with different frequencies oscillated asynchronously with frequency components of both adjacent regions. Provata interpreted her two frequency-domains separated by a 'more frequency' incoherent region as a chimera state. Our multi-frequency cluster (Fig. 1, column IV) exhibits the same features. This can be seen in Fig. 3(a) where we present the absolute value of the Fourier coefficients of the three main frequencies, 24 mHz, 27 mHz, and 32 mHz, of the local Fourier spectra along a the dashed line in the ν plot in Fig. 1 collumn IV. While in the left low-and the right high-frequency regions the contribution of the other two frequencies are very small, in the middle region we find not only the third, dominant frequency at 27 mHz, but also a significant contribution of the frequencies of the two adjacent regions, just as in Provata's model system. In our case, the dynamics of the mediating region appears phenomenologically rather coherent. On this basis, we would not classify our state as a typical chimera state, in contrast to the chimera states observed previously during silicon photoelectrodissolution [9,10].
However, the classification of chimera states is a very intricate matter. This becomes clear when considering the relation of multi-frequency clusters of different spa-tial extensions in continuous systems to weak chimera states in ensembles of four coupled discrete oscillators [30]. Loosely speaking, a chimera state in such a system is characterized by two oscillators being synchronized, with the frequency ν 1 , while the other two oscillators possess frequencies which differ from each other as well as from ν 1 [30]. If, to a first approximation, we neglect the interfacial regions, a multi-frequency cluster state with three different regions in a spatially extended system can likely be reduced to a low-dimensional system of three coupled oscillators where the coupling is weighted by the size of the domains. Thus, the dynamics of multi-frequency clusters in systems with many degrees of freedom is in a sense equivalent to the one of a weak chimera state. With this in mind it appears worthwhile to differentiate between chimera states where the number of incoherent oscillators scales with the system size and states where it does not. Kemeth et al. have presented considerations along these lines and coined the first type of chimera state extensive which would suggest that the three frequency cluster state could be classified as an intensive chimera state [31]. In this respect, an important question to be investigated in the future is whether one can identify general dynamical properties that determine whether a weak chimera state behaves as an intensive or extensive chimera, in the sense defined here when successively increasing the number of oscillators.
Another issue concerning the classification of the dynamics arises when regarding only the global picture of the dominant frequencies, neglecting any spatial information, as shown in Fig. 3(b). Here, the dominant frequencies as found in state IV of Fig. 1 are sorted in ascending order. The first about thirty thousand entries with the value 0 Hz arise from the region where we observe amplitude death. For higher indices, we clearly observe three plateaus. These reflect our three frequency domains. However, the transitions between these plateaus are not sharp but instead occur continuously in a finite index range. As such, this graph is reminiscent of the distribution of dominant frequencies in 2-frequency chimera states [32][33][34], and to some extent also of the ones in hybrid chimera states [14], which are composed of a fully synchronized, a nearly-synchronized, and an incoherent part. Hence, Fig. 3(b) shows again that there might be aspects in the dynamics of multi-frequency clusters in continuous media that are related to those of chimera states.
In conclusion, our experimental observation of multifrequency clusters is not only an exceptional example where a self-organized adaptive coupling was observed in a non-living system, but also reveals important open problems concerning the properties of multi-frequency states in continuous systems, such as their relation to chimera states or requirements on the adaptive coupling for their existence.
The authors would like to thank Felix P. Kemeth, We used a three-electrode setup with an n-doped (1-10 Ωcm) single crystalline (111) Si sample as the working electrode. Before the experiments, the Si sample was equipped with a back contact by thermally evaporating aluminum on to the back and then annealing it at 250 • C for 15 min. Then the front side of the electrode was treated with an oxygen plasma in order to rid it of any organic contamination. The sample was then mounted on a custom-made PTFE sample holder using a conductive silver paste and sealed using red silicone rubber (Scrintex 901, Ralicks GmbH); 15-25 mm 2 of the Si sample were left exposed, forming the active area of the working electrode. The active area was cleaned by wiping the electrode with acetone-drenched precision wipes and sequentially immersing it in acetone, ethanol, methanol, and ultrapure water (18.2 MΩcm) for 10 min each.
The mounted electrode was then placed in the center of the cell with the Hg|Hg 2 SO 4 reference electrode placed behind it. For the counter electrode, we bent a Pt wire (99.99 % Chempur) into a circle and placed it symmetrically in front of the working electrode. In order to control the voltage between working and reference electrode, we used a FHI-2740 potentiostat (electronics laboratory of the Fritz-Haber-Institut, Berlin, Germany).
The aqueous electrolyte had a total volume of 500 ml and contained 0.06 M NH 4 F and 142 mM H 2 SO 4 , yielding a pH of 1 in accordance with the dissociation constants found in Ref. [35]. The electrolyte was purged with argon for 30 min before the experiment and an argon overpressure was kept throughout all measurements via a gas inlet above the electrolyte. The electrolyte was also stirred using a magnetic stirrer at 20 Hz throughout all measurements.
All glassware was cleaned in HNO 3 and subsequently in a 1 M aqueous KOH solution and stored in ultrapure water. Platinum and PTFE parts were cleaned in Piranha solution. All organic cleaning solvents were AnalaR NORMAPUR grade (VWR Chemicals). All electrolyte components were Suprapur grade (Merck).

Ellipsometric imaging
We used the ellipsometric imaging setup sketched in Fig. 4 to monitor changes in the optical path through the oxide layer in-situ. The non-polarized light coming from the LED (Linos, HiLED, 470 nm) becomes elliptically polarized once it passes through the Glan-Thompson prism and the zero-order λ/4 plate. The beam is then reflected off the working electrode at an angle close to the Brewster angle of water and Si α = 70 • . Depending on the length of the optical path at the electrochemical interface, the ratio between the s-and the p-polarized components of the light change. The polarization is then converted into an intensity signal by letting the light pass through a second Glan-Thompson filter. The intensity was measured using a CCD-camera (JAI CV-A50) and digitized using a frame grabber card (PCI-1405, National Instruments). The spatial average of the frame was sampled at 10 Hz and one spatially resolved frame was saved each second. The CCD gives a linear response to the intensity of the incoming illumination, up to a saturation threshold; we present the ellipsometric intensity as a percentage of this threshold.
In general, the light intensity from the LED varies slightly across the electrode. This leads to a variation of the raw ellipsometric intensity ξ (x, t) raw depending on the position x on the electrode. To adjust for this variation, a background correction was applied by subtracting the temporal average of the raw data ξ (x) raw individually at every point. In addition to this background variation of the intensity, the contrast positively correlates with the absolute value of the LED illumination intensity. Hence, a point on the sample that is illuminated with a high background intensity will have a higher contrast. To counter this, we correct each individual pixel by dividing its value by its temporal average. This correction factor is then normalized by multiplying with the spatial average of the temporal average of the entire image.
In total, the correction suppresses the signal from pixels with high temporal average and enhances the signal from pixels with low temporal average. The complete background correction is summarized in Eq. (9): with ξ (x, t) being the corrected local time series and ξ (x, t) raw the spatial average of the temporal average of the raw data. To reduce the noise, we smoothed the data in the temporal direction by using a Savitzky-Golay filter with a 2nd degree polynomial and a 15 point window. In addition, the data was binned into 5x5 pixels bins.

Illumination
Since n-type Si mainly interacts with the electrolyte through valence-band processes (see Eq. (1) in the manuscript), the sample had to be illuminated to allow for anodic oxidation. For this purpose, a linearly polarized He-Ne laser (HNL150L-EC, Thorlabs) was used. The fact that the laser was linearly polarized allowed us to adjust the illumination intensity with a polarization filter mounted on a motorized rotation mount (KPRM1E/M, Thorlabs) and placed directly after the laser. After the polarizer, the beam was widened using a beam expander and then passed through an iris diaphragm. This allowed for the illumination of the entire sample with the central, more uniform, part of the beam. A sketch of the illumination setup can be seen in Fig. 4. Note that, the intensity of the laser was much higher than the intensity of the LED used for the ellipsometric imaging.
At the beginning of the measurement series, we initialized the electrode in a uniformly oscillating state via a potential step from OCP to U app = 4.15 V vs SHE, and then held the voltage constant whilst illuminating with a high illumination intensity.