{'fraction_non_alphanumeric': 0.037903746221512755, 'fraction_numerical': 0.017865027340963895, 'mean_word_length': 4.495334079880553, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bifurcation structure of the flame oscillation

A flame exhibits a limit-cycle oscillation, which is called"flame flickering"or"puffing", in a certain condition. We investigated the bifurcation structure of the flame oscillation in both simulation and experiment. We performed a two-dimensional hydrodynamic simulation by employing the flame sheet model. We reproduced the flame oscillation and investigated the parameter dependences of the amplitude and frequency on the fuel-inlet diameter. We also constructed an experimental system, in which we could finely vary the fuel-inlet diameter, and we investigated the diameter-dependences of the amplitude and frequency. In our simulation, we observed the hysteresis and bistability of the stationary and oscillatory states. In our experiments, we observed the switching between the stationary and oscillatory states. As fluctuations can induce the switching in the bistable system, switching observed in our experiments suggested the bistability of the two states. Therefore, we concluded that the oscillatory state appeared from the stationary state through the subcritical Andronov-Hopf bifurcation in both the simulation and experiment. The amplitude was increased and the frequency was decreased as the fuel-inlet diameter was increased. In addition, we visualized the vortex structure in our simulation and discussed the effect of the vortex on the flame dynamics.


I. INTRODUCTION
Periodic behaviors can be observed in biological [1,2], chemical [3,4], and hydrodynamical [5,6] systems. They have been understood as nonlinear oscillations and have been studied in relation to the bifurcation, synchronization, and pattern formation. One of the simple phenomena that are understood as limit-cycle oscillations is the flame oscillation, which is often called "flame flickering" or "puffing" [7,8]. One of the authors reported that a single candle burns in a stationary manner, whereas the merged flame on a bundle of three-or-more candles oscillates at a typical frequency of around 10 Hz [9]. The study reported that the flame oscillation emerges with a greater amount of fuel supply than a threshold. Flame oscillations are generally observed in two types of diffusion flames: a pool flame, where the fuel vapor is introduced into the system through the evaporation from a liquid pool [10], and a jet flame, where the fuel stream is supplied via a gaseous jet [11]. Note that candle flames are classified as a pool flame. In both systems, the frequency in the oscillatory state is universally proportional to d −1/2 over a wide range of the fuel-inlet diameter d [12]. The flame oscillation has attracted the attention of many researchers not only of fundamental sciences but also in the industrial fields, since the oscillatory combustion leads to more emission of pollutants, such as unburned hydrocarbons and soot, than stationary combustion.
From the viewpoint of dynamical systems, the emergence of the flame oscillation was considered to be a transition from a stationary state to an oscillatory state [22]. Many other types of transitions with regard to flame dynamics have also been investigated: transition between the flame oscillations with and without "pinchoff" [14], transition between periodic, quasiperiodic, and chaotic oscillations with a rotating cylindrical burner [23], and switching between the axisymmetric and asymmetric shapes of oscillating flames [18,24]. It is also reported that two-or-more oscillating flames can synchronize according to the spatial arrangement of fuel inlets [9]. The transitions between various synchronization modes have been observed [9,[17][18][19][25][26][27][28][29][30][31][32][33] and have also been investigated using mathematical models [9,25].
When we consider transitions in terms of bifurcation phenomena, investigation of bifurcation structures allows us to evaluate the validity of mathematical models and can provide insights into the mechanism of the phenomena. However, most studies of flame oscillations have focused on reproducing the phenomena, and only a few studies have focused on the bifurcation structures. Moreno-Boza et al. theoretically investigated the chemohydrodynamic instability of the stationary state [20] with the flame sheet model, in which the reaction rate of the combustion is sufficiently high and the Lewis number, the ratio of the thermal diffusion coefficient to the mass diffusion coefficient, is unity [34]. They compared the result of linear stability analysis with that from hydrodynamic simulation under the assumption that the flame oscillation occurs through a supercritical Andronov-Hopf bifurcation from the stationary state, and reported some discrepancies, which were left for future investigation [20]. The Andronov-Hopf bifurcation is one of the well-known bifurcation structures where the oscillatory state with a finite frequency appears from the stationary state. This bifurcation is classified into two types: the supercritical and subcritical Andronov-Hopf bifurcations [22,35]. Considering that only the finite frequency has been reported for the flame oscillation [10][11][12]19], the bifurcation structure of the flame oscillation should be classified into the supercritical or subcritical Andronov-Hopf bifurcation. To further understand the flame oscillation, it is essential to identify the bifurcation structure embedded in this phenomenon.
Thus, in the present study, we investigated the bifurcation structure of the flame oscillation in both simulation and experiment. We performed a two-dimensional hydrodynamic simulation employing the flame sheet model, and investigated the dependences of the amplitude and frequency on the fuel-inlet diameter. We considered that our two-dimensional simulation model includes the essential factors, such as the combustion reaction, buoyancy, and convective heat transfer, to reproduce the flame behavior. In order to confirm the simulation result, we also constructed an experimental system in which we could finely vary the fuel-inlet diameter, and investigated the dependences of the amplitude and frequency on the diameter. We then identified the bifurcation structure from the results of the simulation and experiment. We also visualized the vortex structure in our simulation and discussed the effect of the vortex on the flame dynamics based on the phase description of the flame oscillation.

II. SIMULATION SETUP
To reproduce the flame oscillation, we considered the hydrodynamics and the combustion dynamics for the gas, which is composed of ethanol, oxygen, and other nonreactive compounds. We performed two-dimensional hydrodynamic simulation assuming that the gas around the flame is a compressible Newtonian fluid. We used the equation of continuity and the Navier-Stokes equation, which are nondimensionalized by the characteristic velocity, length, and density scales. Here, t is the time, ∂ t is the partial differential operator with respect to t, x α (α = 1, 2) is the spatial coordinate, ∂ α is the partial differential operator with respect to x α , and δ αβ is the Kronecker delta. Here and henceforth, we use the subscripts α and β to denote the index of the spatial coordinate and use the summation convention. ρ, v α , and P are the density, velocity, and pressure field, respectively.
Re is the Reynolds number, the ratio of the inertial force to the viscous force, and F r is the Froude number, the ratio of the inertial force to the gravitational force. ρ 0 is the density of the gas with no combustion. In order to consider the ethanol combustion simply, we adopted the flame sheet model [34]. In this model, the temperature and the mass fractions of the ethanol, oxygen, and other compounds commonly depend on a single variable, "mixture fraction" Z [36], and it follows the advection-diffusion equation, Here, P r is the Prandtl number, the ratio of the kinematic viscosity to the thermal diffusivity. The temperature field T depends on Z, as where Z st is the stoichiometric mixture fraction, T 0 is the temperature with no combustion, and T 1 is the highest temperature in the system. T 1 and Z st should depend on the actual fuel type and the initial and/or boundary condition of the mass fraction. In order to ignore the sound wave, we adopted the low Mach number approximation, assuming that the Mach number, which is the ratio of the flow velocity to the sound velocity, in the flame oscillation was negligibly smaller than unity [37]. We also assumed that the gas around the flame is the ideal gas, where the dependence of the density on the mass fraction was ignored. Under these approximations, ρ and T satisfy the equation of state for the ideal gas, For numerical simulation, we considered a rectangular region, −X 1 ≤ x 1 ≤ X 1 , 0 ≤ x 2 < X 2 , and defined the floor as x 2 = 0. We also defined the region for the wick as −d/2 ≤ x 1 ≤ d/2, 0 ≤ x 2 ≤ h, where d and h are the width and height of the wick, respectively. As we assumed a bilateral symmetry of the system, we set the symmetry axis to x 1 = 0. The actual calculation region, half of the system, was thus 0 ≤ x 1 ≤ X 1 , 0 ≤ x 2 ≤ X 2 . The nonslip boundary condition for the velocity, v ⊥ = v = 0, was adopted for the floor and surfaces of the wick. Here, v ⊥ and v are the velocity components in the directions perpendicular and parallel to the boundary, respectively. A slip boundary condition for the velocity, v ⊥ = ∂ ⊥ v = 0, was adopted for the symmetry axis. Here, ∂ ⊥ and ∂ are the partial differential operators with respect to x ⊥ and x , which are the spatial coodinates in the directions perpendicular and parallel to the boundary, respectively. The velocity at x 1 = X 1 and x 2 = X 2 follows the Neumann boundary condition, ∂ ⊥ v ⊥ = ∂ ⊥ v = 0. The pressure at the symmetry axis, the floor, and the surfaces of the wick was determined so that it satisfies the Navier-Stokes equation in Eq. (2). The pressure at x 1 = X 1 , and x 2 = X 2 follows the Dirichlet boundary condition, P = 0. The Neumann boundary condition for the mixture fraction, ∂ ⊥ Z = 0, was adopted for the symmetry axis, the floor, and the side of the wick. Considering that the fuel evaporates from the upper surface of the wick, the Dirichlet boundary condition for the mixture fraction, Z = 1, was adopted for the upper surface of the wick. The mixture fraction at x 1 = X 1 and x 2 = X 2 follows the Dirichlet boundary condition, Z = 0. To solve the equation of continuity and the Navier-Stokes equation, we used the fractional step method [38] and calculated the variables on a staggered grid. We used an explicit method for the advection-diffusion equation to calculate the time evolution of Z. The convective terms were handled by the upwind scheme and the other terms were handled by the central difference scheme. We set the time step dt = 0.01 and the spatial mesh dx 1 = dx 2 = 1. The size of the calculation region was set as X 1 = 50 and X 2 = 100. h was set constant at 10 and d was varied, which could only be an even number due to the bilateral symmetry in a discrete grid. Dimensionless parameters were set as Re = ̺U L/µ = 21.58, P r = µc/κ = 0.7099, and F r = U 2 /gL = 2.041, where L = 0.002 m, U = 0.2 m/s, and ̺ = 1.0 kg/m 3 were the characteristic length, velocity, and density, respectively. µ = 1.853 × 10 −5 Pa · s, κ = 0.026 W/(m · K), and c = 1000 J/(kg · K) were the viscosity, thermal conductivity, and specific heat capacity of the air, respectively, and g = 9.8 m/s 2 was the gravitational acceleration. We set P r = 0.1789 for the upper surface of the wick so that the diffusion may be great enough to provide sufficient amount of ethanol from the upper surface of the wick. The other simulation parameters were set as follows: ρ 0 = 1.18, T 0 = 3, and Z st = 0.2199. T 1 was varied finely as a control parameter. We performed the simulations for d = 8, 10, 12, and 14 in the following procedure. First, we calculated 1,000 steps for T 1 = T 0 = 3 with dt = 0.1, which corresponded to the diffusion process of ethanol into the air before ignition. The initial conditions for v α and Z were v 1 = v 2 = Z = 0 in the whole region. In this situation, buoyancy should not work and convective flow should not occur since we assumed that the density depends only on the temperature as in Eq. (5). Next, we reset T 1 = 12, which corresponded to the ignition, and set t = 0 at this instant. We calculated till t = 3, 000 with dt = 0.01, and the variables in the final state of this calculation were set as an initial condition of the following calculation with different T 1 . We increased T 1 from 12 to 24 in increments of 0.5, and then decreased from 24 to 12 in decrements of 0.5. We reset t = 0 when we varied T 1 , and calculated till t = 3, 000 with dt = 0.01 for each T 1 .

III. SIMULATION RESULTS
The sequential snapshots of the temperature field T and the velocity field v α after a sufficiently long time are shown in Fig. 1. The stationary state, in which the flame height was constant, is shown in Fig. 1(a) for d = 14 and T 1 = 14, while the oscillatory state, in which the flame height oscillated, is shown in Fig. 1 flame, which is indicated by an arrow in each snapshot in Fig. 1. As shown in the velocity fields, we observed a strong upward flow above the wick in both states. The time series of the flame height for the parameters in Fig. 1 Fig. 2(d) was observed depending on the initial condition. The stationary and oscillatory states were observed when T 1 was increased and decreased, respectively, which indicates the bistability of the two states.
In the oscillatory state, the difference between the maximum and minimum values of the flame height was employed for the amplitude, and the frequency at which the Fourier spectrum took the maximum within a range of 0.05-0.20 was employed for the frequency of the oscillation. In the stationary state, the amplitude should be 0. T 1 -dependences of the amplitude and frequency for each d are shown in Fig. 3(a) and (b), respectively. The phase diagram on the d-T 1 plane is shown in Fig. 3(c). The stationary state, in which the amplitude was zero, was observed for small T 1 and d. In contrast, the oscillatory state, in which the amplitude was finite, was observed for large T 1 and d. In the oscillatory state, the amplitude was increased and the frequency was decreased as d was increased for fixed T 1 . In order to investigate the bifurcation structure, we varied T 1 finely around the bifurcation point for each d. . T 1 -dependence of the amplitude is shown in Fig. 3(d), (e), and (f) for d = 10, 12, and 14, respectively. The system was in the stationary or oscillatory state depending on how the parameter had been varied before reaching the same T 1 , which clearly indicated the existence of hysteresis and bistability. These results suggested that the oscillatory state appeared from the stationary state through the subcritical Andronov-Hopf bifurcation.

IV. EXPERIMENTAL SETUP
To compare with the parameter dependences observed in the simulation, we performed experimental observations of the flame for various wick diameters. Figure 4(a) shows the experimental setup. We prepared wicks by cutting diatomite (B161WH, Soil, Japan) into the shape depicted in Fig. 4(b) with a milling machine (monoFab SRM-20, Roland, Japan). The wick diameter d was varied from 8 mm to 28 mm in increments of 2 mm as a control parameter. In order to prevent evaporation of the fuel from the side surface of the wick, we covered the side surface with an aluminum tape (width: 10 mm, thickness: 0.1 mm, B07SF5GCHT, TeenitorJP, China), placed it in a glass Petri dish (diameter: 45.5 mm, height: 18.8 mm) (Fig. 4(c)), and covered it with a lid (D50S, Thorlabs, USA) ( Fig. 4(d)). We performed the experiments in the following procedure. First, we immersed the wick covered with the tape in ethanol (Wako Pure Chemicals, Japan), and degassed it for 10 min in a vacuum chamber in order to let it absorb the ethanol. Next, we prepared the setup shown in Fig. 4(a), and ignited the wick in a darkroom. At the same time, we started recording videos of the flame for 90 s using a high-speed camera (300.16 fps, 540 × 256 pixels, STC-MBS43U3V, Omron Sentech, Japan) equipped with an objective lens (L-600-12, Hozan, Japan). We performed the experiments six times for each d with the above procedure. The wicks were reused up to two times.

V. EXPERIMENTAL RESULTS
The sequential snapshots of the flame taken every 0.02 s are shown in Fig. 5. The stationary state, in which the flame height was constant (Fig. 5(a)), and the oscillatory state, in which the flame height oscillated (Fig. 5(b)), were observed for d = 12 mm and 28 mm, respectively.
In order to quantify the dynamics of the flame shape, we performed the binarization on the brightness of each pixel in the recorded videos, and defined the vertical component of the top position of the flame as the flame height. The time series of the flame height are shown in Fig. 6. We observed the stationary state, the oscillatory state, and the quasiperiodic state. Here, the quasiperiodic oscillation was composed of two-or-more incommensurate frequencies. We also observed the switching among these states.  In order to determine the frequency of the flame oscillation, we performed the Fourier transform on the time series of the flame height. The power spectrum corresponding to each time series shown in Fig. 6 is shown in Fig. 7. In the oscillatory state, there was a peak around 10 Hz and its multiples. In the stationary state, there was no peak (Fig. 7(a)(i)) or a much smaller peak than that in the oscillatory state ( Fig. 7(a)(iii) and (b)(ii)). In the quasiperiodic state, there were some peaks at equal intervals around the highest peak (Fig. 7(c)(ii) and (d)(ii)).
Considering that the state of the system could vary over time, we measured the amplitude and frequency of the flame oscillation every 1 s in the time course. We performed the Fourier transform on the time series of the flame height during the 4 s window around the considered time, and the frequency at which the Fourier spectrum took the maximum in the frequency range between 5 Hz and 30 Hz was employed for the frequency of the oscillation of the flame height. The amplitude of the flame oscillation was defined as the difference between the local maximum and minimum values in the time course of the flame height. Since the frequency of the flame oscillation was about 10 Hz, we detected the local maximum values in the time course of the flame height so that the time interval between them was greater than 0.05 s. The minimum values between the local maximum values were employed for the local minimum values. We averaged the amplitude over the time window of 1 s at every 1 s in its time course. Note that the quasiperiodic state is classified into the oscillatory state through this procedure.
We detected 60 pairs of the amplitude and frequency for 20-79 s from the ignition time in each experiment. The scatter plots of the amplitude and frequency for each d are shown in Fig. 8(a). In the stationary state, the amplitude was about 0.0 mm and the frequency was distributed over a wide range, and in the oscillatory state, the amplitude was finite and the frequency was almost constant at one or two values. Since both the stationary and oscillatory states were observed for the same d, they may be bistable. Since the frequency typically took one of the two values in the oscillatory state for 16 ≤ d ≤ 28 mm, two-types of oscillatory states exist: oscillatory state A, where the frequency is lower (Fig. 7(b)(i), (c)(iii), and (d)(i)), and oscillatory state B, where the frequency is higher (Fig. 7(b)(iii), (c)(i), and (d)(iii)). We also observed the switching between oscillatory states A and B through the stationary state or quasiperiodic state as shown in Figs. 6 and 7.
In order to distinguish the stationary and oscillatory states, we defined a common threshold value of the amplitude as 5.5 mm for all d as shown in Fig. 8(a) and (b). The histograms of the frequency in the oscillatory state for each d are shown in Fig. 8(c). In order to distinguish oscillatory states A and B, we defined a threshold value of the frequency where the count was minimum between the two peaks for each d as shown in Fig. 8(a) and (c). In case there were multiple candidates where the counts were minimum between the two peaks, we employed the value closest to the peak of oscillatory state A. The classification by the thresholds is shown with different markers in Fig. 8(a), and the classification was appropriate.
The results of the mean and standard deviation of the amplitude and frequency for each state are summarized in Fig. 9(a) and (b). As d was increased, the amplitudes increased in both oscillatory states A and B, and the frequency decreased in oscillatory state A. In contrast, the frequency did not clearly decrease in oscillatory state B. The fraction of time duration for which each state was observed is shown in Fig. 9(c). As d was increased, the fraction of the stationary state decreased and that of the oscillatory state increased. The fraction of oscillatory state B was maximum at d = 24 mm.

VI. DISCUSSION
Our experimental results shown in Fig. 8 suggested that the stationary and oscillatory states were bistable, and the noise, such as gentle convective flows and fluctuation of fuel supply, may induce the switching between these states. Therefore, We checked whether a fluctuation induced the switching in the simulation. As the fluctuation of fuel supply, we changed the Prandtl number P r at the upper surface of the wick from P r = 0.1789 into P r = 0.2236 at t = 3000 and changed back into P r = 0.1789 at t = 3000 + ∆t. We show ∆t-dependence of the amplitude for d = 10, T 1 = 23.0 after a sufficiently long time in Fig. 10. When the time duration ∆t of the perturbation was greater than ∆t = 0.5, the system switched from the stationary state to the oscillatory state ( Fig. 10(a)). The system switched from the oscillatory state to the stationary state when 1.6 ≤ ∆t ≤ 2.4 ( Fig. 10(b)). Therefore the switching between the stationary and oscillatory states could be reproduced in the simulation.
As d was increased, the states observed in our experi- ments typically changed in the following order: stationary state, bistable state, and oscillatory state. Thus, the oscillatory state may appear from the stationary state through the subcritical Andronov-Hopf bifurcation, which was in agreement with our simulation results. In the simulation (Fig. 3), the bistable region and the oscillatory amplitude therein were so small that the parameter dependence of the amplitude may be similar to the one under the assumption of the supercritical Andronov-Hopf bifurcation. Therefore, we suggest that the hysteresis and bistability should be adequately evaluated to determine the bifurcation points.
d-dependences of the amplitude and frequency obtained in the simulation were qualitatively consistent with those for oscillatory state A obtained in the experiments. In a previous study, Xia et al. reproduced the power law, in which the frequency was proportional to g/d, with the theoretical analysis based on vortex dynamics [12]. Following their results, the frequency for d ≈ 10 1 mm and g ≈ 10 1 m/s 2 , typical parameters in our experiments, was almost 10 1 Hz, which was at an order consistent with that of our experimental results. Although this power law was observed over three orders of magnitude of g/d, the range of g/d in our exper-iments, which was limited near the bifurcation point as 18.7 s −1 ≤ g/d ≤ 35 s −1 , was so small that we could not check whether the power law held. Chen et al. investigated the parameter dependence of the frequency by varying the number of candles, which could correspond to d in the viewpoint of the fuel supply; the frequency was decreased with an increase in the number of bundled candles [19]. The results in our simulation and experiment were qualitatively consistent with their results.
We experimentally observed the two types of oscillations with different amplitudes and frequencies: oscillatory states A and B. Previously, Cetegen et al. and Bunkwang et al. reported bistability of the axisymmetric and asymmetric states [18,24]. Although the frequencies in oscillatory states A and B were quantitatively consistent with those in the reported axisymmetric and asymmetric states, respectively, the difference in the flameshape symmetry between oscillatory states A and B was not observed in our experiment. Therefore, whether oscillatory state B corresponded to the asymmetric state is still unclear, and further investigation is needed in both three-dimensional hydrodynamic simulation and experiments by focusing on the axial symmetry breaking.
Following some previous studies [13][14][15][16][17][18]21], we checked the flow structure in our simulation to understand the role of the large-scale hydrodynamic structure on the flame oscillation. In order to focus on the time variation of the flow in the oscillatory state, we obtained the velocity fieldv α , which is defined as the deviation from the time-averaged flow field aŝ where τ is the period of the oscillation. The sequential snapshots ofv α after a sufficiently long time are shown in Fig. 11(a) and (b). The vorticity field, is also plotted with the color map. We observed the vortex pair with the upward velocity along the center axis (VPU) and that with the downward velocity along the center axis (VPD) rising alternately. Such vortex dynamics was observed in a previous experimental study, in which the velocity field in the flame oscillation was measured with particle image velocimetry and "invalid velocity vector analysis" based on the proper orthogonal decomposition [13]. Note that the vortex pair in the two-dimensional systems corresponds to the vortex ring in the three-dimensional systems. In order to clarify the effect of the vortex pairs on the oscillation of the flame height, we investigated the timing at which the vortex pair rose. The top position of the flame and the center positions of the vortices, at whicĥ v 1 =v 2 = 0 holds, are indicated in Fig. 11(b), and the time series of the x 2 -components are shown in Fig. 11(c). respectively. We quantified this instant by the phase θ of the oscillation of the flame height. Here, one period corresponds to 0 ≤ θ < 1, and the instant when the flame top is the highest corresponds to θ = 0. The phases when VPU and VPD pass through the mean of the flame height are set as θ = θ U and θ = θ D , respectively. The dependences of θ U and θ D on d and T 1 are shown in Fig. 11(d).
We found that θ U ≃ 1/4 and θ D ≃ 3/4 irrespective of the parameters d and T 1 . Here, we discuss the mechanism of the flame oscillation by focusing on the horizontal flow. An illustration of the behavior during one period of the flame oscillation is shown in Fig. 11(e). The horizon-  tal inflow toward the flame decreased the flame height at θ = θ U , whereas the horizontal outflow increased the flame height at θ = θ D . Thus, we considered that the horizontal outflow and inflow may supply the hotter and colder gases, which may expand and shrink the flame, respectively, under the flame sheet model. It is notable that the vertical flow did not explain the phase relationship, since the upward and downward flow may supply more and less fuel, respectively. The above discussion supports that VPU and VPD cause the decrease and increase in the flame height, respectively, and may enhance the instability of the stationary flame, which leads to the oscillation. The oscillation of flame height and the vortex pair periodically rising may enhance each other with time delay, which can lead to the hysteresis and bistability.

VII. CONCLUSION
We investigated the bifurcation structure of the flame oscillation in a 2D hydrodynamical simulation and experiment. In our simulation, we reproduced the flame oscillation and observed the hysteresis and bistability between the stationary and oscillatory states. In our experiments, we constructed a system in which we could finely vary the wick diameter d as a parameter, and we observed the switching between the stationary and oscillatory states. We concluded that the oscillatory state appeared from the stationary state through the subcritical Andronov-Hopf bifurcation in both the simulation and experiment. In addition, the amplitude was increased and the frequency was decreased as the wick diameter d was increased in both the simulation and experiment.
Following the previous studies, we analyzed the flow structure in our simulation and observed the vortex pairs rising periodically; this provides insights into the mechanisms of the flame oscillation. In addition, we also experimentally observed the quasiperiodic oscillation, which has never been reported before, and another oscillatory state, where the frequency was higher. Understanding of these oscillations may contribute to controlling the flame, and thus, they need to be investigated in detail in the future. Identification of the bifurcation structure of the flame oscillation may also facilitate the investigation of synchronization of coupled flame oscillations using mathematical models.

VIII. ACKNOWLEDGMENTS
This work was supported by JSPS KAKENHI Grant Nos. JP19H00749 and JP21K13891 (HI), and Nos.