Experimental study on the bifurcation of a density oscillator depending on density difference

Limit-cycle oscillations have been investigated intensively since they are related to a wide variety of dynamical phenomena both in natural and artificial systems. There are mainly two scenarios for the realization of oscillations. One is a harmonic-like oscillation realizing in an energy conservative system, and the other is a limit-cycle oscillation in which energy is alternately supplied and dissipated in one cycle. Thus, the observation of oscillation in open systems is the mark of the existence of a limit-cycle oscillation. Since energy dissipation is unavoidable in macroscopic and mesoscopic systems, many oscillatory phenomena observed in such systems can be regarded as limit-cycle oscillations.

One of the most famous experimental systems that exhibit limit-cycle oscillation is the chemical oscillation represented by the Belousov-Zhabotinsky (BZ) reaction, in which oxidation and reduction processes occur alternately and the solution color changes according to such oxidation/reduction reactions [1,2]. Other types of oscillators have also been reported; e.g., Briggs-Rauscher (BR) reaction, glycolysis, and electrolysis [3]. Another class is hydrodynamical systems such as the Bénard convection system and the density oscillator. These hydrodynamical systems are often discussed in association with the atmospheric circulation and thermohaline circulation [4,5]. Martin firstly reported the density oscillator, and he discussed it as a simplified experimental system of thermohaline circulation [6].

The limit-cycle oscillations in the hydrodynamical systems are important not only for the relation to such atmospheric and thermohaline circulations but also as the suitable model systems to investigate the bifurcation phenomena seen in dynamical systems. The hydrodynamical systems have many degrees of freedom, but they often show successive bifurcation structure from a trivial state to oscillatory states, quasi-periodic oscillation states, and chaotic states that can be described by a model with a small number of variables [7].

One of the simple and appropriate experimental systems that exhibit hydrodynamic limit-cycle oscillations is a density oscillator. In this system, the gravitational instability originating from an upset density profile of higher- and lower-density fluids induces upstream and downstream alternation, which can be regarded as a limit-cycle oscillation. From Martin's first report in 1970, there have been experimental and theoretical papers on the density oscillators. Some are on the theoretical analysis on the mechanism of oscillation [8–15] and others are on the coupling between two or more oscillators [16–21]. For the bifurcation phenomenon, some theoretical studies have been performed based on the reduced ordinary differential equations, which have predicted various types of bifurcations according to various bifurcation parameters [10,11]. However, the experimental analyses on the bifurcation structure with physical parameters between the resting and oscillatory states are missing, and the detailed behaviors close to the bifurcation point remain unclear.

In the present paper, we report the experimental results on the bifurcation between the resting and oscillatory states by changing the density of the higher-density fluid as a bifurcation parameter. We measured the oscillation amplitude and period depending on the density of the higher-density fluid, which indicate the bifurcation structure between the resting and oscillatory states.

All aqueous solutions and distilled water were prepared with the Elix water purification system Elix UV 3 (Merck, Darmstadt, Germany). A sodium chloride (NaCl) aqueous solution as a higher-density fluid was prepared by dissolving NaCl (Wako Pure Chemicals, Tokyo, Japan) into distilled water with various weight/volume concentrations c. All the aqueous solutions and pure water ($c = 0\ \text{g/L}$, for control) were degassed for 1 hour in a vacuum chamber before use. The observation chamber of the density oscillator was made of acrylic plates, which is composed of a smaller container surrounded by a larger container for higher- and lower-concentration solutions, respectively. The acrylic plates were 10 mm in thickness except for the bottom of the smaller container with 2 mm thickness. The center of the bottom of the smaller container had a small cylindrical hole with a diameter of 1 mm and a length of 2 mm to connect these two containers. They were kept isolated by blocking the hole with a needle before the measurements. The dimensions of the observation chamber and a photograph of the experimental setup are shown in fig. 1.

Zoom In Zoom Out Reset image size

We put 245 mL of the NaCl aqueous solution in the smaller container, and 1400 mL of pure water in the larger container to reach the same water levels in these two containers. After removing the blocking of the small hole, the flow was induced through the hole if the concentration difference was greater than a critical value. The time series of the surface height of the solution in the outer container was measured with a laser displacement meter (32.77 fps, LT9010M, Keyence, Japan) [12]. We measured the surface height for 4000 s. As the surface height adequately relaxed to a steady resting or oscillatory state in 2000 s, we analyzed the data from 2000 s to 4000 s. It should be noted that the steady oscillatory state with the repetitive macroscopic flow is characteristic for limit-cycle oscillations. To visualize the profile of the density difference between higher and lower NaCl concentration, we utilized the optical-index difference in a similar manner with the shadowgraph method [22] using a light emitting diode (LED) as a light source and a plastic sheet (TAMIYA Inc., Shizuoka, Japan) as a screen. The obtained time series of the surface height were processed as follows: The surface height linearly decreased due to the evaporation, and thus the linear trend was reduced by subtracting the least-squares fitted linear function. Then, the data were smoothed by applying a band-pass filter between 0.003 and 0.019 Hz. Using the smoothed data and their time derivatives, we detected the local maximum and minimum points. The time series of the time derivative were calculated as the slope of adjacent $\pm10$ points (local data points for 0.64 s) obtained by the least-squares fitting. We obtained the mean peak-to-peak distances and the mean time intervals of the local minima as the oscillation amplitude and the period, respectively, for each condition of NaCl concentration c.

Figure 2 shows the upstream and downstream series through the small hole connecting the higher- and lower-concentration solutions. Characteristic oscillatory behavior of the density oscillator is observed. The concentration of the NaCl aqueous solution was $c = 30\ \text{g/L}$. In this setup, the upstream of the lower-concentration solution and the downstream of the higher-concentration solution were visualized as dark and light images, respectively, since the cylindrical flow shape acts as concave or convex optical lenses. The upstream and downstream series were repeated with a typical period of ${\sim}50\ \text{s}$ in this condition. In this timescale, the cylindrical flows did not mix with the surrounding solutions, i.e., diffusion of solutes is negligible. Figures 2(b) and (c) show the details for the switching of these opposite flows, where the arrow-like shaped flow enters the other solution at an almost constant velocity.

Zoom In Zoom Out Reset image size

For the quantitative measurement of the oscillation, we used a time series of the surface height of the solution in the outer container obtained by the laser displacement meter. The characteristic time series of the surface height is shown in fig. 3. Figure 3(a) shows the result for $c = 0\ \text{g/L}$ as a control condition. The unprocessed time series (blue line) shows a constantly-decreasing trend due to the slow evaporation of water, even though a total volume change of the aqueous phase is negligibly small. By subtracting the linear trend (black line), which is obtained by the least-squares method, time evolution of the surface height h(t) is obtained. h(t) directly reflects the fluid flow through the small hole between the inner and outer containers. While the time series for $c = 0\ \text{g/L}$ is almost constant with small fluctuation (fig. 3(a)), that for $c = 30\ \text{g/L}$ exhibits steady oscillation (fig. 3(b)), which corresponds to the oscillatory flow through the hole. Note that the linear trends for both resting $(c = 0\ \text{g/L})$ and oscillatory $(c = 30\ \text{g/L})$ states were measured in 3 hours under the condition with similar humidity, and showed a similar decreasing rate with the slopes of $-0.00746\ \mu \text{m}\cdot \text{s}^{-1}$ and $-0.00772\ \mu \text{m}\cdot \text{s}^{-1}$, respectively, which are compatible with the typical rate for water evaporation. For further quantification of the amplitude and period of the oscillation, we applied a band-pass filter in a frequency range 0.003–0.019 Hz (red-dashed line).

Zoom In Zoom Out Reset image size

Figure 4 shows the concentration dependence of the typical time series of the surface height for $c = 0\ \text{g/L}$ (control), $c = 1.5\ \text{g/L}$, $c = 3\ \text{g/L}$, $c = 10\ \text{g/L}$, and $c = 30\ \text{g/L}$. In addition to the original trend-subtracted data (gray points), we plotted smoothed data obtained by a moving average with adjacent $\pm64$ points (3.91 s) (black line) and the data with a band-pass filter (red-dashed line) (figs. 4(a)–(e)). The corresponding whole dynamics including the initial transient states are shown in fig. S2 in the SM. While the surface heights for small c, e.g., $c = 0\ \text{g/L}$, are in a resting state, those for $c \geq 1.5\ \text{g/L}$ exhibit characteristic oscillation. According to the increase in c, the amplitude becomes larger and the period becomes shorter. The detailed wave profiles, shown in each inset, also change from smooth sinusoidal-like oscillations to relaxation oscillations. Figures 4(f)–(j) show the corresponding trajectories of the trend-subtracted smoothed data in a phase space $(h,\mathrm{d}h/\mathrm{d}t)$. According to the increase in c, the trajectories corresponding to the limit-cycle orbits with larger amplitudes appear. Also in the trajectories for $c \geq 10\ \text{g/L}$, we can clearly see a characteristic shape of relaxation oscillation. Though the smoothed plot for small oscillation in $c = 1.5\ \text{g/L}$ shows no clear cyclic trajectory due to the long-term fluctuation of the baseline, the plot after a band-pass filter shows sufficiently large oscillatory amplitude, which will be discussed later.

Zoom In Zoom Out Reset image size

From the band-pass filtered data, we obtained the amplitudes and the periods of the oscillations for various concentrations c, as shown in figs. 5(a) and (b), respectively. The results in fig. 5(a) show that the amplitude shows a steep increase above a critical concentration $c_\mathrm{c}\simeq 1\ \text{g/L}$, and it monotonically increases for larger concentrations. Here, we determined one standard deviation of the amplitude in the control condition (c = 0 g/L) as the threshold amplitude between the resting and oscillatory states. By this standard, $c \leq 1\ \text{g/L}$ is in a resting state, and $c \geq 1.125\ \text{g/L}$ shows a limit-cycle oscillation. For the conditions with limit-cycle oscillations $c \geq 1.125\ \text{g/L}$, the oscillation period T is plotted against various c in fig. 5(b). The period T monotonically decreases with the increase in the concentration c, and converges to $T \simeq 50\ \text{s}$ for larger c.

Zoom In Zoom Out Reset image size

The driving force of the density oscillator is the gravitational energy. In fact, the density profile is upset in the initial stage; in other words, the higher-density fluid is located above the lower-density fluid. This instability can induce the limit-cycle oscillation through the dissipation of the gravitational energy. From the experimental results, the limit-cycle oscillation was observed when the concentration of NaCl aqueous solution, c, was above a critical value $c_\textrm{c} \simeq 1 \ \text{g/L}$. Considering that the gravitational energy is proportional to the density, the results suggest that the oscillatory flow is induced by the upset density profile but it should suffer from the resistance due to the fluid viscosity. The gravitational energy due to the upset profile can overcome the viscosity for $c > c_\textrm{c}$, while it cannot for $c < c_\textrm{c}$.

The existence of the critical value $c_\textrm{c}$ for the instability is also interpreted from the transport phenomenon. Since the solute transportation is dominated by diffusion in the resting state and by the advection in the oscillatory state, the Rayleigh number $Ra = (gL^3\Delta\rho/\rho)/(\nu D)$ provides a critical value for the density difference Δρ to drive the oscillation. Here, $g\,({\sim}10\,\mathrm{m\cdot s^{-2}})$ is the gravitational acceleration, $L\,({\sim}0.5\,\mathrm{mm})$ is the hole radius as a characteristic length, $\Delta\rho/\rho\,({\sim}10^{-2})$ is the ratio of the density difference, $\nu\,({\sim}1.0\times10^{-6}\ \text{m}^{2}\cdot \text{s}^{-1})$ is the kinematic viscosity of water at room temperature, and $D\,({\sim}10^{-9}\ \text{m}^{2}\cdot \text{s}^{-1})$ is the diffusion constant of solutes in water at room temperature. If we take a critical Rayleigh number $Ra\,({\sim}10^3)$ typical for Rayleigh-Bénard instability, the critical density difference $\Delta\rho_\textrm{c}$, namely, the critical solute concentration $c_\textrm{c}$ can be estimated as $c_\textrm{c} \sim 0.8\ \text{g/L}$, which agrees well with the bifurcation point ${\sim}1.0\ \text{g/L}$ obtained in our experiment. Note that this value is much smaller compared to the parameter ranges in the previous experiments [8,12]. Actually, the decreases in the density of salt water per a cycle are roughly estimated as ${\sim}{-}3\%$ in Steinbock et al. [8] and ${\sim}{-}0.2\%$ in Ueno et al. [12], which is much larger than our experiment, ${\sim}{-}0.003\%$. These previous experiments showed continuous shifts of the water surface in mm–cm length scale during the measurements due to the decrease in the density difference by mixing of the solutions between the two containers. In contrast, our experiment measured the motion of surface height in μm length scale, and therefore, the bifurcation phenomenon at the small density difference as well as the water evaporation could be detected.

From the viewpoint of the dynamical systems, this transition can be discussed in terms of bifurcation. There are several scenarios known for the transition from a resting state to an oscillatory state; supercritical Hopf bifurcation, saddle-node bifurcation of cycles (subcritical Hopf bifurcation), infinite-period bifurcation (saddle-node homoclinic bifurcation, saddle-node bifurcation on an invariant circle (SNIC)), homoclinic bifurcation, and so on [23–26]. For a supercritical Hopf bifurcation, the amplitude of oscillation raises from a bifurcation point in proportion to the square root of the distance from the bifurcation point, while the amplitude has a finite value at the bifurcation point for the other three types.

In the experiments, the amplitude of the oscillation seems to increase from zero with an increase in c for $c > c_{\rm c}$ as shown in fig. 5(a). If we assume that the system undergoes the supercritical Hopf bifurcation, the dependence of the amplitude on the concentration c is suggested to be proportional to $(c-{c_\mathrm{c}^\mathrm{H}})^{1/2}$ when the system is close to the bifurcation point. Here, $c_\mathrm{c}^\mathrm{H}$ is the critical concentration for the supercritical Hopf bifurcation. We performed power-law fittings of the experimental data to determine the bifurcation point and the power exponent. First, we checked the power spectra for each concentration c to determine the range of the "close-to-bifurcation-point region", where the critical behavior is expected. Figure 6 shows the power spectra corresponding to the conditions shown in fig. 4. Above $c_\mathrm{c}$, a clear first peak for the characteristic limit-cycle oscillation appears within a frequency range of the band-pass filter, 0.003–0.019 Hz. For higher concentrations, $c\geq3\ \text{g/L}$, harmonics originating from the nonlinear waveform of relaxation oscillation appear, i.e., the system is in the "far-from-bifurcation-point region". Thus, we performed linear fitting of the squared amplitude in the "close-to-bifurcation-point region" as shown in fig. 7(a), and obtained the critical concentration as ${c_\mathrm{c}^\mathrm{H}} = 1.03\ \text{g/L}$. Here, the closest oscillatory condition $c = 1.125\,\mathrm{g/L}$ was eliminated from the fitting, because the amplitude was close to the noise level. Then, the double logarithmic plot of the amplitude vs. $c-{c_\mathrm{c}^\mathrm{H}}$ shown in fig. 7(b) provides the power exponent close to the bifurcation point, resulting in 0.482, which is close to $1/2$. It suggests that the bifurcation would be classified into the supercritical Hopf bifurcation.

Zoom In Zoom Out Reset image size

The supercritical Hopf bifurcation is also characterized by a finite angular velocity in the phase space at the bifurcation point. In fig. 5(b), the period increased as c approached the bifurcation point from higher concentrations $c > c_\mathrm{c}$. This behavior recalls infinite-period bifurcation, in which the amplitude is finite while the period diverges at the bifurcation point $c_\mathrm{c}^\mathrm{I}$ with the scaling of $(c-c_\mathrm{c}^\mathrm{I})^{-1/2}$. We further checked the critical behavior of the period vs. $c-c_\mathrm{c}^\mathrm{H}$ in the double logarithmic plot shown in fig. 7(c), and confirmed that the period does not seem to diverge. Instead, it remains constant around the bifurcation point, which is also consistent with the supercritical Hopf bifurcation. Moreover, we also estimated another critical concentration $c_\mathrm{c}^\mathrm{I}$ by assuming the infinite-period bifurcation through the critical behavior of the inversed square period. This test resulted in $c_\mathrm{c}^\mathrm{I} = 4.87\times10^{-2}\ \text{g/L}$ as shown in fig. 7(d), which is apparently inconsistent with the experimentally indicated $c_\mathrm{c}\simeq1.0\ \text{g/L}$ as observed in the amplitude shown in fig. 5(a). While our experimental results and analyses strongly support the supercritical Hopf bifurcation, the precise asymptotic behaviors are still not clear due to the experimental limitations. From the experimental observation only, it is difficult to completely exclude the possibilities for other types of bifurcations. From the present results, the saddle-node bifurcation of cycles was dismissed since the bistability between the resting and oscillatory states was not observed. Further study is needed to identify the bifurcation class of the density oscillator depending on the density difference.

Zoom In Zoom Out Reset image size

We investigated the transition from the resting state to the oscillatory state depending on the density difference in a density oscillator. The limit-cycle oscillation, where the upstream and downstream alternations occur, was observed for the higher density difference, while no hydrodynamic flow was observed for the lower density difference. The detailed data close to the bifurcation point provide a critical exponent close to $1/2$ and a finite period around the bifurcation point, which is consistent with the supercritical Hopf bifurcation. Further experimental and theoretical studies should be performed to exactly identify the bifurcation class since the experimental results still showed the ambiguity.

The density oscillator is an excellent experimental system for the limit-cycle oscillation with hydrodynamic instability since it can be treated from the standpoint of physics; complex processes such as chemical reaction and phase transition are not included, and thus only hydrodynamics is taken into consideration. The present experimental study will give fundamental knowledge on the nonlinear oscillation with hydrodynamic instability from the viewpoint of bifurcation theory in dynamical systems, and help further studies using density oscillators.

This work was supported by JSPS KAKENHI Grant Nos. JP19K14675, JP16H03949. It was also supported by the Japan-Poland Research Cooperative Program "Spatio-temporal patterns of elements driven by self-generated, geometrically constrained flows" and the Cooperative Research Program of "Network Joint Research Center for Materials and Devices" (No. 20191030).