Phase description of oscillatory convection with a spatially translational mode

We formulate a theory for the phase description of oscillatory convection in a cylindrical Hele-Shaw cell that is laterally periodic. This system possesses spatial translational symmetry in the lateral direction owing to the cylindrical shape as well as temporal translational symmetry. Oscillatory convection in this system is described by a limit-torus solution that possesses two phase modes; one is a spatial phase and the other is a temporal phase. The spatial and temporal phases indicate the position and oscillation of the convection, respectively. The theory developed in this paper can be considered as a phase reduction method for limit-torus solutions in infinite-dimensional dynamical systems, namely, limit-torus solutions to partial differential equations representing oscillatory convection with a spatially translational mode. We derive the phase sensitivity functions for spatial and temporal phases; these functions quantify the phase responses of the oscillatory convection to weak perturbations applied at each spatial point. Using the phase sensitivity functions, we characterize the spatiotemporal phase responses of oscillatory convection to weak spatial stimuli and analyze the spatiotemporal phase synchronization between weakly coupled systems of oscillatory convection.


Highlights:
We develop a phase reduction theory for oscillatory convection with a spatial mode. The theory can be considered as a phase description method for limit-torus solutions. We derive phase sensitivity functions for spatial and temporal phases of convection. We can quantify spatiotemporal phase responses of convection to weak perturbations. We can analyze spatiotemporal phase synchronization between weakly coupled systems.
We recently developed a phase description method for limit-cycle solutions to the following partial differential equations: the nonlinear Fokker-Planck equations that represent the collective dynamics of globally coupled noisy dynamical elements [19], the fluid equations that represent the dynamics of the temperature field in ordinary Hele-Shaw cells [20], and the reaction-diffusion equations that represent rhythmic spatiotemporal patterns in chemical and biological systems [21]. However, there are also examples of spatiotemporal rhythms in systems that further possess spatial translational symmetry; these spatiotemporal rhythms cannot be described by limit-cycle solutions.
For example, rotating annuli and spheres possess continuously rotational symmetry, i.e., continuously translational symmetry in the rotating direction [17,18,22,23]. Consequently, the emergence of spatiotemporal rhythms in such systems brings up two phase modes, i.e., a spatial phase and a temporal phase. Such spatiotemporal rhythms are described by limit-torus solutions. Synchronization of spatiotemporal rhythms with two phase modes has been experimentally investigated using systems of rotating fluid annuli that exhibit traveling and oscillating convection, which is analogous to atmospheric circulation [24,25]. Therefore, a phase description method for limit-torus solutions to partial differential equations is desirable.
In this paper, as the first step, we consider oscillatory convection in a cylindrical Hele-Shaw cell that is laterally periodic. An ordinary Hele-Shaw cell is a rectangular cavity where the gap between two vertical walls is much smaller than the extent of the other two spatial dimensions, and the fluid in the cavity exhibits oscillatory convection under the appropriate parameter conditions (see Refs. [26,27] and references therein). The cylindrical Hele-Shaw cell is a cylindrical version of the ordinary Hele-Shaw cell that possesses spatial translational symmetry in the lateral direction owing to the cylindrical shape. Oscillatory convection in the cylindrical Hele-Shaw cell is therefore described by a limit-torus solution that possesses both spatial and temporal phases.
Here, we formulate a theory for the phase description of oscillatory convection in the cylindrical Hele-Shaw cell. The theory can be considered as a phase reduction method for limit-torus solutions to partial differential equations. The theory can also be considered as a generalization of our phase description method for limit-cycle solutions to partial differential equations such as the nonlinear Fokker-Planck equations [19], fluid equations [20], and reactiondiffusion equations [21]. The phase reduction method for limit-torus solutions enables us to describe the dynamics of the oscillatory convection by two phases (i.e., spatial and temporal phases), and facilitates theoretical analysis of the spatiotemporal phase synchronization properties of the oscillatory convection. On the basis of phase reduction, we characterize the spatiotemporal phase responses of oscillatory convection to weak impulses and analyze the spatiotemporal phase synchronization between weakly coupled systems exhibiting oscillatory convection.
This paper is organized as follows. In Sec. II, we formulate a theory for the phase description of oscillatory convection with a spatially translational mode; supplemental information of the theory is given in App. A and App. B. In Sec. III, we illustrate the theory using a numerical analysis of the oscillatory convection. In Sec. IV, we make a comparison between the theory and direct numerical simulations. Concluding remarks are given in Sec. V.

II. PHASE DESCRIPTION OF OSCILLATORY CONVECTION
In this section, we formulate a theory for the phase description of oscillatory convection in a cylindrical Hele-Shaw cell that is laterally periodic. The theory can be considered as an extension of our phase description method for oscillatory convection in the ordinary Hele-Shaw cell [20] to that in the cylindrical Hele-Shaw cell.

A. Dimensionless form of the governing equations
The dynamics of the temperature field T (x, y, t) in the cylindrical Hele-Shaw cell is described by the following dimensionless form (see Ref. [26] and references therein): The Laplacian and Jacobian are respectively given by where we assumed that the curvature effects due to the cylindrical shape are negligible (see Refs. [28,29] for curvature effects, although the subject of these references is not thermal convection but viscous fingering). The first and second terms on the right-hand side of Eq. (1) represent diffusion and advection, respectively. The stream function ψ(x, y, t) is determined from the temperature field T (x, y, t) as follows: where the Rayleigh number is denoted by Ra. The stream function ψ(x, y, t) also gives the fluid velocity field v(x, y, t), i.e., v(x, y, t) = ∂ψ ∂y , − ∂ψ ∂x .
(5) Figure 1 shows a schematic diagram of the cylindrical Hele-Shaw cell. The system is defined in the following rectangular region: x ∈ [0, 2] and y ∈ [0, 1]. Because this Hele-Shaw cell has a cylindrical shape, the system possesses a 2-periodicity with respect to x. The boundary conditions for the temperature field T (x, y, t) are given by T (x, y, t) where the temperature at the bottom (y = 0) is higher than that at the top (y = 1). The stream function ψ(x, y, t) satisfies the periodic boundary condition on x and the Dirichlet zero boundary condition on y, i.e., Owing to the homogeneity of Eqs. (1)(4) and the periodic boundary condition on x, given in Eqs. (6)(8), this system possesses continuous spatial translational symmetry with respect to x. We also note that no conserved quantity exists in this system.

B. Convective components of the temperature field
To simplify the boundary conditions in Eq. (7), we consider the convective component X(x, y, t) of the temperature field T (x, y, t) as follows: Substituting Eq. (10) into Eq. (1) and Eq. (4), we derive the following equations: and Applying Eq. (10) to Eqs. (6)(7), we obtain the following boundary conditions for the convective component X(x, y, t): That is, the convective component X(x, y, t) satisfies the periodic boundary condition on x and the Dirichlet zero boundary condition on y.
In the derivation below, it should be noted that Eq. (12) can also be written in the following form: where the Green's function G(x, y, x , y ) is the solution to under the periodic boundary condition on x and the Dirichlet zero boundary condition on y. From the translational symmetry with respect to x, the Green's function G(x, y, x , y ) possesses the following property: G(x, y, x , y ) = G(x − x , y, 0, y ). In the following two subsections, we analyze the dynamical equation (11) with Eq. (12) or Eq. (15) under the boundary conditions given by Eqs. (13) (14) and Eqs. (8)(9).
C. Limit-torus solution and its Floquet-type system In general, a stable limit-torus solution to Eq. (11), which represents oscillatory convection in the cylindrical Hele-Shaw cell, can be described by (e.g., see Fig. 4 in Sec. III) The spatial phase and traveling velocity are denoted by Φ and c, respectively; the temporal phase and oscillation frequency are denoted by Θ and ω, respectively. The spatial and temporal phases indicate the "position" and "oscillation" of the convection, respectively. Figure 2 shows a schematic diagram of the limit-torus solution X 0 (x − Φ, y, Θ). The limit-torus solution X 0 (x − Φ, y, Θ) satisfies the 2-periodicity with respect to Φ and the 2π-periodicity with respect to Θ, i.e., Substituting Eq. (17) into Eqs. (11)(12), we find that the limit-torus solution X 0 (x − Φ, y, Θ) satisfies the following equation: where the stream function ψ 0 (x − Φ, y, Θ) is determined by (e.g., see Fig. 5 in Sec. III) From Eq. (10), the corresponding temperature field T 0 (x, y, Θ) is given by (e.g., see Fig. 5 in Sec. III) Let u(x − Φ, y, Θ, t) represent a small disturbance to the limit-torus solution X 0 (x − Φ, y, Θ), and consider a slightly perturbed solution Equation (11) is then linearized with respect to u(x − Φ, y, Θ, t) as follows: Here, the linear operator L(x − Φ, y, Θ) is explicitly given by where Similarly to the limit-torus solution X 0 (x − Φ, y, Θ), the function u(x − Φ, y, Θ) satisfies the periodic boundary condition on x, the Dirichlet zero boundary condition on y, and the 2π-periodicity with respect to Θ. In Eq. (26), the function ψ u (x − Φ, y, Θ) is the solution to under the periodic boundary condition on x and the Dirichlet zero boundary condition on y. Note that the linear operator L(x − Φ, y, Θ) is periodic with respect to both Φ and Θ. Therefore, Eq. (24) is a Floquet-type system with two zero-eigenvalues; one is associated with spatial translational symmetry breaking and the other is associated with temporal translational symmetry breaking.
Defining the inner product of two functions as we introduce the adjoint operator of the linear operator L(x − Φ, y, Θ) by By partial integration, the adjoint operator L * (x − Φ, y, Θ) is explicitly given by where Similarly to u(x − Φ, y, Θ), the function u * (x − Φ, y, Θ) also satisfies the periodic boundary condition on x, the Dirichlet zero boundary condition on y, and the 2π-periodicity with respect to Θ. In Eq. (31), the two functions, i.e., ψ * u,x (x, y, Θ) and ψ * u,y (x, y, Θ), are the solutions to under the periodic boundary condition on x and the Dirichlet zero boundary condition on y, respectively. Details of the derivation of the adjoint operator L * (x − Φ, y, Θ) are given in App. A.

D. Floquet zero eigenfunctions
In the calculation below, we utilize the Floquet eigenfunctions associated with the two zero-eigenvalues, i.e., Note that the two right zero eigenfunctions, i.e., U s (x − Φ, y, Θ) for the spatial phase Φ and U t (x − Φ, y, Θ) for the temporal phase Θ, can be chosen as (e.g., see Fig. 7 in Sec. III) which are confirmed by differentiating Eq. (20) with respect to Φ and Θ, respectively. For the inner product (28) with the two right zero eigenfunctions (38) (39), the corresponding two left zero eigenfunctions, i.e., U * s (x − Φ, y, Θ) and U * t (x − Φ, y, Θ), are orthonormalized as for p, q = s, t. Here, we note that the following equation holds (see also Refs. [6,19,20]): Therefore, the following orthonormalization condition is independently satisfied for each Θ: Here, we describe a numerical method for obtaining the left zero eigenfunctions. From Eqs. (36)(37), the left zero eigenfunctions, for p = s, t, which can be transformed into by substituting Θ = −ωs. A relaxation method using Eq. (44), which can also be called the adjoint method ( see Refs. [6][7][8][9][10] for limit-cycle solutions to ordinary differential equations and Refs. [19][20][21] for limit-cycle solutions to partial differential equations ), is convenient to obtain the left zero eigenfunctions for the limit-torus solution. In the following two subsections, we derive a set of phase equations for oscillatory convection in the cylindrical Hele-Shaw cell using the limit-torus solution and its Floquet zero eigenfunctions.

E. Oscillatory convection with weak perturbations
In this subsection, we consider oscillatory cylindrical-Hele-Shaw convection with a weak perturbation to the temperature field T (x, y, t) described by the following equation: The weak perturbation is denoted by p(x, y, t). Substituting Eq. (10) into Eq. (45), we obtain the following equation for the convective component X(x, y, t): Using the idea of phase reduction [2], we derive a set of phase equations from the perturbed equation (46). That is, using the left zero eigenfunctions, i.e., U * s (x − Φ, y, Θ) and U * t (x − Φ, y, Θ), we project the dynamics of the perturbed equation (46) onto the unperturbed limit-torus solution with respect to the spatial and temporal phases as follows: where we approximated X(x, y, t) by the unperturbed limit-torus solution X 0 (x − Φ, y, Θ) in Eqs. (47) (48). Therefore, the two phase equations describing the oscillatory cylindrical-Hele-Shaw convection with weak perturbation are approximately obtained in the following forms: where the phase sensitivity functions for the spatial and temporal phases are defined as (e.g., see Fig. 8 in Sec. III) The phase equations (49)(50) are the main results of this paper. These equations can also be considered as an extension of that describing oscillatory convection in the ordinary Hele-Shaw cell [20]. As found from Eqs. (49)(50), the spatial and temporal phases are coupled; therefore, nontrivial spatiotemporal phase dynamics are revealed. Furthermore, we consider the case of the perturbation described by a product of two functions as follows: That is, the space-dependence and time-dependence of the perturbation are separated. In this case, the phase equations (49)(50) are written in the following forms: where the effective phase sensitivity functions for the spatial and temporal phases are given by (e.g., see Fig. 10 in Sec. III) We note that the forms of Eqs. (54)(55) are essentially the same as those of the phase equations which are derived for a perturbed limit-torus oscillator described by a finite-dimensional dynamical system (see Refs. [30,31]). Finally, it should be noted that we can also consider oscillatory cylindrical-Hele-Shaw convection with weak boundary forcing as mentioned in App. B.

F. Weakly coupled systems of oscillatory convection
In this subsection, we consider weakly coupled systems of oscillatory cylindrical-Hele-Shaw convection described by the following equation: for (σ, τ ) = (1, 2) or (2, 1). Two identical systems of oscillatory cylindrical-Hele-Shaw convection are mutually coupled through corresponding temperatures at each spatial point [56], where the coupling parameter is denoted by . Substituting Eq. (10) into Eq. (58) for each σ, we obtain the following equation for the convective component As in Eq. (12), the stream function ψ σ (x, y, t) of each system is determined by Here, we assume that unperturbed oscillatory cylindrical-Hele-Shaw convection is a stable limit-torus solution and that the coupling between the systems of oscillatory convection is sufficiently weak. Under this assumption, as in the preceding subsection, we obtain a set of phase equations from Eq. (59) as follows: whereΓ These two functions, i.e., , depend on the spatial phase difference between the systems and the temporal phases of both systems.
Introducing the slow phase variables as we rewrite Eqs. (61)(62) asφ By applying the averaging method with respect to the temporal phases, Eqs. (67)(68) are written in the following forms:φ where the phase coupling functions for the spatial and temporal phases are given by (e.g., see Fig. 12 in Sec. III) Therefore, we obtain the following phase equations: where the phase coupling functions are explicitly described as The phase coupling functions for the spatial and temporal phases, i.e., depend only on the spatial and temporal phase differences. Let the spatial and temporal phase differences respectively be defined as From Eqs. (73)(74), we obtain the following equations by subtraction: where These two functions, i.e., Γ which represent the anti-symmetry with respect to the origin, i.e., ∆Φ = ∆Θ = 0. Finally, we note that the forms of Eqs. (73)(74) are the same as that of the phase equations which are derived from weakly coupled limit-torus oscillators described by finite-dimensional dynamical systems (see Refs. [30,31]). That is, a system of oscillatory convection with a spatially translational mode can be reduced to a set of phase equations, similarly to an ordinary limit-torus oscillator.

III. NUMERICAL ANALYSIS OF OSCILLATORY CONVECTION
In this section, we illustrate the theory developed in Sec. II using a numerical analysis of the oscillatory convection in the cylindrical Hele-Shaw cell.

A. Spectral transformation
For numerical simulations using the pseudospectral method performed in Sec. III B, we first describe a spectral transformation. Considering the boundary conditions for X(x, y, t), given in Eqs. (13)(14), we introduce the following spectral decomposition: where the spectral transformation of X(x, y, t) is given by Similarly, the limit-torus solution X 0 (x − Φ, y, Θ) introduced in Eq. (17) is decomposed as The corresponding spectral transformation of the limit-torus solution is described by where we assign the origin of the spatial phase, i.e., Φ = 0, to the spatial pattern X 0 (x, y, Θ) that satisfies the following property: which is unique when a pair of vortices exist in the system, as considered below (e.g., see Fig. 4 and Fig. 5 in Sec. III B). When visualizing the temporal phase of the limit-torus in the infinite-dimensional state space, we project the limit-torus solution X 0 (x − Φ, y, Θ) onto the H 0,2 -H 0,4 plane as which are real numbers and depend only on the temporal phase Θ. When determining the spatial phase of the limit-torus solution, we introduce the following complex order parameter: which corresponds to Eq. (88) with j = −1 and k = 1. From Eq. (89), the spatial phase Φ is determined by As in X(x, y, t), considering the boundary conditions for Z s (x − Φ, y, Θ) and Z t (x − Φ, y, Θ), i.e., the periodic boundary condition on x and the Dirichlet zero boundary condition on y as described in Eqs. (A20)(A21), we also introduce the following spectral decompositions: where the spectral transformations are given by The spatial power spectra of the phase sensitivity functions averaged over the temporal phase are defined as (99)

B. Limit-torus solution and phase sensitivity functions
We first summarize our numerical simulations in this paper. As mentioned in the preceding subsection, we applied the pseudospectral method, which is composed of a Fourier expansion with 256 modes for the periodic boundary condition on x and a sine expansion with 128 modes for the Dirichlet zero boundary condition on y. The initial values were prepared such that the system exhibits oscillatory convection with a pair of vortices. Because the Rayleigh number was fixed at Ra = 400, the traveling velocity and oscillation frequency were c = 0 and ω 532, respectively. As mentioned below, this limit-torus solution possesses reflection symmetry because the traveling velocity is exactly zero, i.e., c = 0. Here, we note that the theory developed in Sec. II is applicable for the case of non-zero traveling velocity. Figure 3 shows the limit-torus orbit projected onto the H 0,2 -H 0,4 plane, which was obtained from our numerical simulations of the dynamical equation (11). Snapshots of the limit-torus solution X 0 (x − Φ, y, Θ) are shown in Fig. 4. In addition, several quantities associated with the limit-torus solution are shown as follows: snapshots of the temperature field T 0 (x − Φ, y, Θ) and the stream function ψ 0 (x − Φ, y, Θ) are shown in Fig. 5; snapshots of the fluid velocity v 0 (x − Φ, y, Θ) = (∂ y ψ 0 (x − Φ, y, Θ), −∂ x ψ 0 (x − Φ, y, Θ)) are shown in Fig. 6; snapshots of the right zero eigenfunctions, U s (x − Φ, y, Θ) and U t (x − Φ, y, Θ), are shown in Fig. 7; snapshots of the phase sensitivity functions, Z s (x − Φ, y, Θ) and Z t (x − Φ, y, Θ), are shown in Fig. 8. We note that the phase sensitivity functions, Z s (x − Φ, y, Θ) and Z t (x − Φ, y, Θ), were obtained using the adjoint method, i.e., the relaxation method for Eq. (44) with the orthonormalization condition given by Eq. (42).
As seen in Fig. 5, the spatial phase Φ can be considered as the "position" of the hot plume in the convection, whereas the temporal phase Θ represents the "oscillation" of the convection. We also note that the phase sensitivity functions, Z s (x − Φ, y, Θ) and Z t (x − Φ, y, Θ), are spatially localized as seen in Fig. 8. Namely, when the spatial phase is Φ = 1, the amplitudes of Z s (x − Φ, y, Θ) and Z t (x − Φ, y, Θ) with respect to the temporal phase Θ in the bottom-left, bottom-right, and top-center regions of the system are much larger than those in the other regions; this fact reflects the dynamics of the spatial pattern of the convective component X 0 (x − Φ, y, Θ) shown in Fig. 4.
The limit-torus solution representing the oscillatory convection with Eq. (100) possesses the following reflection symmetry (see Fig. 4): From Eq. (38) and Eq. (39), the right zero eigenfunctions, U s (x − Φ, y, Θ) and U t (x − Φ, y, Θ), respectively possess the following reflection anti-symmetry and reflection symmetry (see Fig. 7): Therefore, the phase sensitivity functions, or the left zero eigenfunctions, also possess the following properties (see Fig. 8): Moreover, in this case, the limit-torus solution X 0 (x − Φ, y, Θ) and the phase sensitivity functions, Z s (x − Φ, y, Θ) and Z t (x − Φ, y, Θ), also possess point symmetry. For each Θ, the limit-torus solution and phase sensitivity functions possess the following properties with respect to a certain point of the system (see Fig. 4 and Fig. 8): where x δ = x − Φ/2 and y δ = y − 1/2.

C. Effective phase sensitivity functions
In this subsection, we calculate the effective phase sensitivity functions obtained in Sec. II E. Before illustrating the effective phase sensitivity functions, the spatial power spectra of the phase sensitivity functions averaged over the temporal phase as defined in Eq. (98) and Eq. (99), i.e., P s (j, k) and P t (j, k), are shown in Fig. 9(a) and Fig. 9(b), respectively. Owing to the point-symmetry of the phase sensitivity functions, given in Eqs. (107)(108), both P s (j, k) and P t (j, k) exhibit checkerboard patterns. The power of the mode (j, k) = (7, 3) is the largest for P s (j, k) and P t (j, k). To illustrate the effective phase sensitivity functions, we consider a corresponding spatial pattern a(x, y) = cos(πjx) sin(πky), where j = 7 and k = 3. Figure 9(c) shows the spatial pattern a(x, y), for which the effective phase sensitivity functions are shown in Fig. 10 with respect to Φ and Θ. In addition, the effective phase sensitivity functions with Φ = 0.75 are shown in Fig. 11 as a function of Θ. As seen in Fig. 10 and Fig. 11, the effective phase sensitivity functions, ζ s (Φ, Θ) and ζ t (Φ, Θ), exhibit both positive and negative values. Namely, when the effective phase sensitivity function for the spatial phase, ζ s (Φ, Θ), is positive (negative), the spatial phase Φ is advanced (delayed) by applying a positive perturbation. Similarly, when the effective phase sensitivity function for the temporal phase, ζ t (Φ, Θ), is positive (negative), the temporal phase Θ is advanced (delayed) by applying a positive perturbation.

IV. COMPARISONS WITH DIRECT NUMERICAL SIMULATIONS
In this section, we compare the theoretical values obtained in Sec. III with direct numerical simulations of oscillatory convection.

A. Spatiotemporal phase responses of oscillatory convection to weak impulses
In this subsection, we make a comparison of the effective phase sensitivity functions between the theoretical values obtained in Sec. III C, i.e., Fig. 11, and direct numerical simulations of oscillatory convection with weak impulses described by Eq. (46) using Eqs. (53)(109). The comparison of the effective phase sensitivity functions, ζ s (Φ = 0.75, Θ) and ζ t (Φ = 0.75, Θ), between the direct numerical simulations with impulse intensity and the theoretical curves are shown in Fig. 13. The simulation results agree quantitatively with the theory.

B. Spatiotemporal phase synchronization between weakly coupled systems of oscillatory convection
In this subsection, we make a comparison on time evolution of phase differences between the theoretical values obtained in Sec. III D, i.e., Fig. 12(c), and direct numerical simulations of two weakly coupled systems of oscillatory convection described by Eq. (59). The comparison of the time evolution of the phase differences between the direct numerical simulations with coupling intensity and the theoretical curves are shown in Fig. 14. The simulation results agree quantitatively with the theory.

V. CONCLUDING REMARKS
Our investigations in this paper are summarized as follows. In Sec. II, we formulated the theory for the phase description of oscillatory convection with a spatially translational mode. In particular, we derived the phase sensitivity functions for the spatial and temporal phases. Details of the derivation of the adjoint operator, which provides the phase sensitivity functions, are given in App. A. Treatments of the boundary forcing are described in App. B. We also derived the phase coupling functions for the spatial and temporal phases. In Sec. III, we illustrated the theory using the numerical analysis of the oscillatory convection. In particular, we obtained the phase sensitivity functions and phase coupling functions. In Sec. IV, we made comparisons between the theory and direct numerical simulations: the spatiotemporal phase responses of oscillatory convection to weak impulses and the spatiotemporal phase synchronization between two weakly coupled systems of oscillatory convection. The theoretical predictions were successfully confirmed by the direct numerical simulations.
Second, we note the spatial reflection symmetry of spatial patterns. Consider the pattern formation in a system of spatial reflection symmetry: when a spatial pattern does not break the reflection symmetry, the traveling velocity is zero, c = 0; meanwhile, when a spatial pattern does break the reflection symmetry, the traveling velocity becomes non-zero, c = 0. The traveling velocity of the oscillatory cylindrical-Hele-Shaw convection is zero; however, the phase description method itself is applicable for the case of non-zero traveling velocity. If the reflection symmetry of limittorus solution, i.e., Eq. Third, we note ubiquitousness of the limit-torus solutions to partial differential equations. A limit-torus solution to an ordinary differential equation represents a quasi-periodic oscillator, and a phase description method for the quasiperiodic oscillator has also been developed [30,31]. However, there have been few studies on the synchronization of quasi-periodic oscillators; this may be due to the fact that limit-cycle or chaotic oscillations are ubiquitous, but quasiperiodic oscillations are rather rare in ordinary differential equations. In contrast, a limit-torus solution is ubiquitous in a partial differential equation that possesses some spatial translational symmetry. From this point of view, a systematic analysis of a set of phase equations, such as Eqs. (54)(55), Eqs. (73)(74), and Eqs. (79)(80), is meaningful and important. It should also be noted that these phase equations are universal and invariant for limit-torus solutions.
Finally, we remark the broad applicability of our phase description approach, which is not restricted to oscillatory cylindrical-Hele-Shaw convection. The phase description method can be generalized to traveling and oscillating localized convection (i.e., a traveling breather) in a binary fluid system (see, e.g., Ref. [49]) or traveling and oscillating convection in a rotating fluid annulus system (see, e.g., Ref. [50]). Similar phase description methods can also be developed for oscillating spots of zero traveling velocity (see, e.g., Ref. [51]) or traveling breathers of non-zero traveling velocity (see, e.g., Ref. [52]) in reaction-diffusion systems.
In this appendix, we describe the details of the derivation of the adjoint operator L * (x−Φ, y, Θ) given in Eqs. (30)(31) (see also, e.g., Refs. [53,54] for mathematical terms). The derivation procedure is similar to that performed in Ref. [20]. From Eqs. (25) (26), the linear operator L(x − Φ, y, Θ) is given by the following form: Using the Green's function G(x, y, x , y ) given in Eq. (16), the function ψ u (x − Φ, y, Θ) given in Eq. (27) can also be written in the following form: In Eqs. (A6)(A7), we perform the following manipulations: where we used the following abbreviations: and also defined the following functions: Here, we note that Eqs. (32)(33) can be derived by applying the Laplacian to Eqs. (A14)(A15), respectively. In this way, the adjoint operator L * (x − Φ, y, Θ), defined in Eq. (29), is obtained as Similarly to the boundary conditions for u(x − Φ, y, Θ) as the adjoint boundary conditions are given by Each term of the bilinear concomitant S[u * (x − Φ, y, Θ), u(x − Φ, y, Θ)] vanishes for the following reasons: the 1st, 2nd, 5th, 7th, 9th, 10th, and 11th terms in Eq. (A23) become zero owing to the 2-periodicity with respect to x for all the functions; the 3rd, 4th, 6th, and 8th terms, the Dirichlet zero boundary condition on y for u and u * ; the last term (i.e., the 12th term), the 2π-periodicity with respect to Θ for all the functions.

Appendix B: Oscillatory convection with weak boundary forcing
In this appendix, we consider oscillatory cylindrical-Hele-Shaw convection with weak boundary forcing described by the following equation: where the stream function ψ(x, y, t) is determined from the temperature field T (x, y, t) as The above two equations are a reproduction of Eqs. (1)(4) for readability. The boundary conditions for the temperature field T (x, y, t) are now given by T (x, y, t) where the weak boundary forcing applied at the bottom (y = 0) and the top (y = 1) is described by p B (x, t) and p T (x, t), respectively. The stream function ψ(x, y, t) satisfies the Dirichlet zero boundary condition on y, i.e., Extending the transformation given by Eq. (10), we consider the following transformation of the temperature field (see also, e.g., Ref. [55] for this type of transformation): where the function P (x, y, t) is defined as We note that P (x, y, t)| y=0 = p B (x, t) and P (x, y, t)| y=1 = p T (x, t). By considering the boundary forcing, the stream function ψ(x, y, t) is also decomposed as Substituting Eqs. (B6)(B8) into Eq. (B1), we derive the following equation: where the first-order terms associated with the boundary forcing are given by Applying Eqs. (B6)(B8) to Eq. (B2) and also considering the linearity of Eq. (B2), we obtain the following equations: We note that Eq. (B11) corresponds to Eq. The two stream functions, ψ X (x, y, t) and ψ P (x, y, t), also satisfy the Dirichlet zero boundary condition on y, i.e., In this way, oscillatory convection with boundary forcing is found to be exactly described by Eq. (B9), which possesses the form of Eq. (46) from the viewpoint of the phase reduction. When the boundary forcing is absent, i.e., = 0, the system is assumed to exhibit oscillatory cylindrical-Hele-Shaw convection described by the following limit-torus solution: which is a reproduction of Eq. (17) for readability. Under the assumption that the boundary forcing is sufficiently weak, as in Sec. II E, using the phase sensitivity functions given by Eqs. (51)(52), we derive a set of phase equations from Eq. (B9) as follows:Φ where the effective boundary forcing function is given by We note that the two functions, X 0 and P 0 , are given by Eq. (20) and Eq. (21), respectively. We also note that the second-order term of Eq. (B9), i.e., 2 J(ψ P , P ), is negligible owing to the smallness of . Furthermore, we consider the case in which the weak forcing is spatially homogeneous as follows: For this case, the effective boundary forcing function is simplified as In summary, the phase description method is applicable to the oscillatory cylindrical-Hele-Shaw convection with weak forcing applied to the boundary given in Eqs. (B3)(B4) as well as to the bulk given in Eq. (45). Because the spatial extent of the cylindrical-Hele-Shaw cell is quasi-two-dimensional, external forcing can be easily applied not only to the boundary but also to the bulk in experiments; however, in general, although external forcing can be easily applied to boundary, it may be difficult to apply external forcing to bulk of fluid systems in experiments. In fact, in the rotating fluid annulus experiments [24,25], perturbations are applied to the boundary. Therefore, the treatments of boundary forcing as developed in this appendix are required to apply the phase description method to such fluid systems.       ). The spatial pattern a(x, y) of the perturbation is given by Eq. (109) and is shown in Fig. 9(c).