Spatial localisation beyond steady states in the neighbourhood of the Takens--Bogdanov bifurcation

The coincidence of a pitchfork and Hopf bifurcation at a Takens-Bogdanov (TB) bifurcation occurs in many physical systems such as double-diffusive convection, binary convection and magnetoconvection. Analysis of the associated normal form, in one dimension with periodic boundary condition, shows the existence of steady patterns, standing waves, modulated waves and travelling waves, where the values of coefficients of the terms in the normal form classify all possible different bifurcation scenarios in the neighbourhood of the TB bifurcation (Dangelmayr&Knobloch, 1987). In this work we develop a new and simple pattern-forming PDE model, based on the Swift-Hohenberg equation, adapted to have the TB normal form at onset, which allows us to explore the dynamics in a wide range of bifurcation scenarios, including in domains much wider than the lengthscale of the pattern. We identify two bifurcation scenarios in which coexistence between different types of solutions is indicated from the analysis of the normal form equation. In these scenarios, we look for spatially localised solutions by examining pattern formation in wide domains. We recover two types of localised states, that of a localised steady state in the background of the trivial state and that of a spatially localised travelling wave in the background of the trivial state which have previously been observed in other systems. Additionally, we identify two new types of spatially localised states: that of a localised steady state in a modulated wave background and that of a localised travelling wave in a steady state background. The PDE model is easy to solve numerically in large domains and so will allow further investigation of pattern formation with a TB bifurcation in one or more dimensions and the exploration of a range of background and foreground pattern combinations beyond steady states.


Introduction
The Takens-Bogdanov bifurcation exhibits a variety of dynamical behaviours. The bifurcation occurs when a Hopf and pitchfork bifurcation coincide, and arises in physical problems such as double-diffusive convection in a horizontal layer of fluid heated from below (Rucklidge, 1992). This system has two 2 of 26 ALRIHIELI, RUCKLIDGE AND SUBRAMANIAN competing gradients that drive motion in the fluid: the temperature gradient and the solute gradient. With low solute gradient, the first bifurcation from the resting (trivial) state as the temperature gradient is increased is a pitchfork bifurcation leading to steady convection. With larger solute gradient, the bifurcation changes to a Hopf bifurcation leading to oscillatory convection. In double-diffusive convection with idealised boundary conditions, these two forms of convection set in with the same horizontal wavelength (Veronis, 1965). In two dimensions with periodic boundary conditions the point where the pitchfork and Hopf bifurcation coincide is called the Takens-Bogdanov point. The normal form that describes the amplitude close to the Takens-Bogdanov (TB) point for a system with O(2) symmetry is (Dangelmayr & Knobloch, 1987) z = µz + A|z| 2 z + ε νż +C (żz + zż) z + D|z| 2ż + O(ε 2 ), ε 1 (1.1) where z is the complex amplitude of the pattern, µ and ν are the unfolding parameters, A,C and D are constants, the dot denotes differentiation with respect to time, and ε controls how close the system is to the TB point. Different bifurcation scenarios obtained by the analysis of the amplitude equation are found close to onset (Dangelmayr & Knobloch, 1987;Knobloch, 1986), with several different types of patterns: steady states (SS), travelling waves (TW), standing waves (SW) and modulated waves (MW). In domains that are many times wider than the preferred wavelength, extended TW, SW, and MW solutions have been found in numerical investigations of the partial differential equations (PDEs) for thermosolutal convection (Deane et al., 1988;Spina et al., 1998;Turton et al., 2015). In a similar scenario, in binary convection, the application of a thermal gradient to a mixture sets up a competing concentration gradient due to the Soret effect. In this system, a transition from SS to TW has been observed in numerical simulations (Zhao & Tian, 2015;Barten et al., 1995a), while a nonlinear SW solution has been numerically obtained by Matura et al. (2004) and Jung et al. (2004).
In addition to patterned states that are spatially extended, i.e., span the entire domain, parameter values where both the trivial state and a periodic SS state are both dynamically stable allow for the existence of spatially localised states. In the subcritical regime with coexistence between the trivial state and the periodic SS branches, spatially localised steady states (LSS) in a background of the trivial state undergoing homoclinic snaking have been obtained in numerical investigations of thermosolutal convection (Beaume et al., 2011) and binary convection (Batiste et al., 2006;Mercader et al., 2009). The snaking branches behave like those familiar from the Swift-Hohenberg equation (Burke & Knobloch, 2007). At a given Rayleigh number, odd and even branch solutions with different number of rolls can be found.
For binary convection, the system undergoes a subcritical Hopf bifurcation to TW for negative separation ratio (Zhao & Tian, 2015). In the parameter regime where the TW bifurcate subcritically from the conduction state, localized travelling waves (LTW) have also been obtained. The LTW solution refers to the spatially localized cells whose envelope moves with a characteristic speed in a background of the trivial state. In contrast to LSS, the LTW have fixed and uniquely selected width, which was discovered in experimental (Kolodner, 1994(Kolodner, , 1991cNiemela et al., 1990;Kolodner, 1991a) and numerical (Barten et al., 1991;Taraut et al., 2012;Barten et al., 1995b) studies of binary convection, with a negative separation ratio = −0.08. This was also observed later in numerical simulations of the full system of binary convection with different but still small negative separation ratios of −0.123 (Watanabe et al., 2012) and −0.1 (Zhao & Tian, 2015).
In this paper, we develop a new and simple model as a useful description of the qualitative features of double-diffusive convection. Our model is a PDE based on the Swift-Hohenberg equation but adapted to SPATIAL  have the TB normal form at onset. This allows for an exploration of the dynamics of localised steady and time dependent patterns in very wide domains. Our model can access most of the bifurcation scenarios that occur in the TB normal form and so it is relevant to other pattern-forming systems whose dynamics can be reduced to a TB normal form.The model recovers LSS and LTW as documented above, as well as two new localised patterns: LSS in an oscillating background and LTW moving through a background of SS.
In Section 2, we develop the linear part of the model by reproducing the dynamics of doublediffusive convection. In Section 3, we discuss the nonlinearities which we can add to the model, taking into account Lyapunov stability. In Section 4, the model is reduced to the Taken-Bogdanov normal form by applying a weakly nonlinear analysis. In Section 5, we identify parameter combinations in the model at which we can observe different dynamical behaviours close to the TB bifurcation. In Section 6, we obtain localized SS with trivial state background and localized SS with SW background using time stepping and observe snaking in the branch of localized SS with trivial state background using continuation. Localized TW with trivial state and SS background are discussed in Section 7 using time stepping. We conclude in Section 8.

Designing the linear dynamics near the TB bifurcation
The first part of designing a model system that has a Takens-Bogdanov (TB) bifurcation requires the possibility for both a pitchfork and a Hopf bifurcation. We build a minimal model for the TB bifurcation by reproducing the dynamics of double-diffusive convection, starting with the design of the linear dynamics. Two different density gradients drive the dynamics in double-diffusive convection: thermal gradients quantified with a thermal Rayleigh number Ra and solutal gradients quantified with a solutal Rayleigh number Rs. The stable quiescent state in the system becomes unstable with increasing thermal gradients and starts to convect. When Rs is less than a critical value Rs c , the quiescent state undergoes a pitchfork bifurcation leading to steady convection as the temperature gradient Ra increases. At large solutal gradients with Rs > Rs c , this behaviour changes and the quiescent state loses stability via a Hopf bifurcation leading to oscillatory convection. The point where the primary bifurcation changes from pitchfork to Hopf bifurcation with (Ra, Rs) = (Ra c , Rs c ) is called the Takens-Bogdanov point, as shown in Fig. 1(a).
In order to replicate this behaviour, we use two control parameters ν and µ in the new model, where the change of sign in ν corresponds to the loci of pitchfork bifurcations and the change of sign in µ (with ν < 0) corresponds to the loci of Hopf bifurcations, as shown in Fig. 1. Such an identification allows us to decompose the behaviour close to the TB point as follows. Starting from large negative values, variations of parameters above the diagonal µ = ν, as shown in Fig. 1(b), allow for the occurrence of pitchfork bifurcation, while parameter variations below this diagonal causes a Hopf bifurcation to occur, followed by a pitchfork bifurcation. In this way we are able to replicate the different bifurcation scenarios from double-diffusive convection.
The second factor to include is the variation of the linear stability threshold with wavenumber k. In the case of double-diffusive convection, the linear stability thresholds for pitchfork and Hopf bifurcations are as shown in Fig. 2. In all three cases, we observe that both the pitchfork and Hopf thresholds vary with the wavenumber k. For example, the threshold for pitchfork bifurcations (shown in pink) has a minimum Ra = Ra 0 at k = k c = π √ 2 . Near this critical wavenumber k c , the critical Rayleigh number varies with the square of the distance of the current wavenumber k from k c . This means that close to k = k c , Ra 0 can be expanded as where Ra 0 = d 2 Ra 0 /dk 2 . In order to reproduce this behaviour, we incorporate such a variation into the two linear stability thresholds by writing µ and ν in terms of (k 2 c − k 2 ) 2 as (2.2) Here the constant b is used to change the rate at which Hopf bifurcation varies with (k 2 c − k 2 ) 2 . From the case of double diffusion shown in Fig. 2, we see that we require b > 1. Other choices can be made depending on the specific system of interest.
Close to the TB point, the PDEs that govern small amplitude thermosolutal convection can be reduced to the linear second order Van der Pol-Duffing equation (Veronis, 1965) where u is the amplitude of the lowest-order mode of the stream function, κ, λ are unfolding parameters and the dot indicates the derivative with respect to time. We start with this dynamical equation and look to fit in our two parameters in terms of the unfolding parameters in this model below. This is equivalent to an eigenvalue problem for σ of the form (2.5) A pitchfork bifurcation (with one zero eigenvalue) occurs when the determinant is zero, i.e., λ = 0, and a Hopf bifurcation (with purely imaginary eigenvalues) occurs when the trace is zero and the determinant is positive, i.e., κ = 0 and λ < 0. At the TB point, the system has two zero eigenvalues, i.e. κ = λ = 0. We can relate µ and ν from Eqn. (2.2) to this two dimensional dynamical system by setting (2.6) This gives the relations that λ = 0 when µ = (k 2 cPF − k 2 ) 2 , and κ = 0 when ν = b(k 2 cHopf − k 2 ) 2 . This implies that the linear operator in terms of Eqn. (2.6) is then given as (2.7) Consequently, we can write the linear equation (2.3) as where u(t) is now the amplitude at lowest order of the mode e ikx . This equation in Fourier space can be converted into a PDE by replacing k 2 with −∂ 2 /∂ x 2 and considering u to be a function of x and t. The ODE (2.8) converted to a PDE is where u tt = ∂ 2 u/∂t 2 and u t = ∂ u/∂t are the derivatives with respect to time, and the parameters k cPF and k cHopf are the wavenumbers for the pitchfork and Hopf bifurcations. The general choice of k cPF = k cHopf , is relevant to problems where the pitchfork and Hopf bifurcation have different critical wavenumbers, e.g., magnetoconvection (Chandrasekhar, 1961;Arter, 1983;Weiss, 1981;Proctor & Weiss, 1982) and rotating convection (Clune & Knobloch, 1993;Veronis, 1966Veronis, , 1968Zimmermann et al., 1988). Since we are interested in thermosolutal convection, where the pitchfork and Hopf bifurcations have the same critical wavenumbers, for this paper we let k cPF = k cHopf = 1. Then the linear second order partial differential equation that models the dynamics at small amplitudes near a TB point takes the following from: The dispersion relation can be determined by studying the eigenvalues of the model given by (2.11) The growth rate σ is a function of the wavenumber k, so some modes could decay and while others can grow. If the eigenvalues for all k have negative real parts the evolution decays and the zero solution is linearly stable. If any eigenvalue for any k has positive real part the evolution grows and the zero solution is linearly unstable.

Selection of the nonlinear terms
The model given in Eqn. (2.10) is a second order linear partial differential equation in time which has been designed to reproduce the linear stability results of double-diffusive convection. This section deals with determining the choice of nonlinear terms in the model, with two aims: first, to have a globally stable nonlinear dynamical system and second, to be able to access the variety of dynamical behaviours that can occur in the neighbourhood of the TB bifurcation. We consider only those nonlinearities which are invariant with respect to reflection about the x-axis (to replicate this symmetry in double-diffusive convection) and include nonlinearities both in the field u and its time derivative u t . We can additionally classify the nonlinear terms in the order at which they contribute to the dynamics. For example, quadratic nonlinear terms can include while cubic nonlinearities can include terms such as Similar terms (without the u t contributions) have been incorporated into various generalisation of the Swift-Hohenberg equation (Burke & Dawes, 2012;Kozyreff & Tlidi, 2007;Crawford & Riecke, 1999). Including all the nonlinearities mentioned would make the model very complicated. We consider the two criteria identified at the beginning of this section and proceed to create a candidate nonlinear model of the form, In this model we have chosen simple polynomial nonlinearities in the form of u 2 , u 3 and u 3 t along with the additional mixed nonlinear term u 2 u t , each with a constant coefficient.
The first step is to ensure global stability during evolution. In order to test this, we look to determine a Lyapunov function for the dynamics in the spirit of the Swift-Hohenberg equation. With this in mind, we consider the Lyapunov functional where L is the length of the periodic domain. Our aim is to show that this function is bounded below and decreases with time for large u and u t . Having C 1 < 0 is sufficient for F (t) to be bounded from below. Differentiating F (t) with respect to time we have If u t = 0 for all x, then we see that the above relation vanishes with dF /dt = 0 and we have an equilibrium. Now, we consider a non-equilibrium point in the dynamics with large non-zero values of u and u t and want show that the rate of change of F is still negative. In order to do this, we can simplify the expression above a bit further. At large u and u t , the last two quartic terms dominate, so we write them as follows: (3.7) We can re-cast the states into an amplitude-phase form with where R(x,t) is a large radius and φ (x,t) is an angle. Substituting (3.8) into (3.7) and simplifying, we get If trajectories are to remain bounded for any choice of ν, dF /dt must be negative for any (u, u t ) large enough as long as u t is not zero for all x. This is guaranteed if T < 0 for any φ as long as sin φ is not zero for all x, which requires Therefore the model second order partial differential equation takes the following form where µ and ν are the control parameters, b > 1, Q 1 ,C 1 ,C 2 and C 3 are constants coefficients with C 1 < 0 to make F bounded below, and we choose C 2 < 0 and C 3 < 0 to make F decrease with time at large values of u or u t . Now that we have identified conditions on the coefficients of the nonlinear terms in order to ensure global stability, we look at our second goal: to access all the possible dynamical behaviours close to a TB bifurcation. The process to check this in detail is discussed in section 5. In order to increase the number of scenarios that are possible, we add one quadratic nonlinear term and two cubic nonlinear terms as shown below.
We do not calculate the Lyapunov functional for this updated model, but rather rely on the conditions that we have derived previously in terms of C 1 , C 2 and C 3 to provide global stability.

Reduction to the Takens-Bogdanov normal form
In this section, we use weakly nonlinear theory to reduce our model PDE from Eqn. (3.12) to the TB normal form in Eqn. (1.1). We consider ε 1 to be a small parameter that parameterises the small amplitude solution such that u = O(ε). This solution can be found in the neighbourhood of the TB bifurcation point where µ and ν vary with values of the order of O(ε 2 ). This scaling satisfies both the normal forms of a pitchfork and a Hopf bifurcation and has been used to analyse the TB problem previously (Knobloch & Proctor, 1981;Dangelmayr & Knobloch, 1987). Therefore, we scale the field, time and the parameters µ and ν as follows: By substituting the scalings into the governing equation (3.12), we get the following equation, where we write the terms up to O(ε 4 ) explicitly. Note that C 3 is not present in the above equation since it contributes only at the sixth order of ε. Matching terms of order ε and ε 2 we get where L is the linear operator defined as (4.5) The linear Eqn. (4.3) obtained at O(ε) can be solved for u 1 by assuming exponential solutions of the form (4.6) Here F 1 and its complex conjugate are functions of time. Substituting this and its derivatives into Eqn. (4.4), we can solve for u 2 as where G 1 andḠ 1 are functions of time.
At O(ε 3 ), we have We can solve for u 3 by assuming it to be of the form The equation (4.8) has contributions into the e ix lengthscale. However, this lengthscale has already been accounted for in the solution for u 1 . Therefore, we need to enforce the condition that the coefficient of e ix terms are zero as a solvability condition. Substituting (4.6), (4.7) and (4.9) into (4.8) and collecting the terms that contribute at the e ix lengthscale, we get We can then solve for the unknowns H 0,2,3 by collecting the coefficients of the constant, e 2ix and e 3ix respectively as (4.12) (4.14) At O(ε 4 ), we have (4.15) The solvability condition from this order of equations is the e ix component of (4.15) as below. (4.20) In order to ensure that both the solvability conditions arising from both the O(ε 3 ) and O(ε 4 ) equations are satisfied, we use a reconstitution procedure (Rucklidge & Silber, 2009) to combine equations (4.10) and (4.16) into a single PDE by defining a new variable By unscaling time and the parameters according to (4.23) Substituting F 1 = z ε − εG 1 in (4.23) and collecting terms of O(1) in the above equation, we get the TB normal form as In this section, we identify parameter combinations in the model at which we can observe different dynamical behaviours close to the TB bifurcation. In order to do this, we relate the parameters in the current model with those investigated in detail by Dangelmayr & Knobloch (1987). In the rest of this paper, we refer to this paper as DK. DK identify different bifurcation scenarios close to a TB bifurcation and classify them in terms of the value of coefficient A as well as the ratio D/(2C + D) where 2C + D is defined as M. The ratio can be related to the coefficients in our model as below.
We consider a few special cases below to illustrate the need for certain nonlinear terms in the model. First, let us consider the case with Q 2 = C 4 = C 5 = 0. Then using expressions in Eqns. (4.11) and (4.20), the values of A and D will be negative and the fraction D/M is always bounded between the values In this range of ratios, we can only access one type of bifurcation behaviour near the TB bifurcation identified from the classification given in DK as A < 0, case II−. Secondly, we consider the case with only C 4 = C 5 = 0, which gives the expression for the ratio D/M as Considering that Lyapunov stability requires C 1 ,C 2 and C 3 to be negative, we consider the case with C 1 = C 2 = C 3 = −1 and b = 2. By plotting contours of A and D M = c, where c = 1 5 , 1 2 , 0.7, 0.74, 3 4 , 4 5 , 1, we can obtain a range of normal form cases. Figure 3 shows the regions of different cases accessible in this case as a function of the two quadratic coefficients Q 1 and Q 2 . From this figure, we see that we can access all bifurcation scenarios with A < 0, while in the case with A > 0, we are unable to access cases I− and IV−. For other choices of parameters, it is possible to access all cases on the normal form with M < 0, while still satisfying the Lyapunov stability requirement. However, cases of the normal form with M > 0 still remain inaccessible with C 4 = C 5 = 0.
Thirdly, when we allow for all nonlinear terms in the model, we can access all cases of the normal form listed in DK, i.e., 18 cases with A < 0 and 8 cases with A > 0 both with M < 0 and M > 0. As a first reference, we give a list of parameters that allow us to reach all the cases with M < 0 in Table 1. Each of these cases has a different bifurcation scenario close to the TB bifurcation.
Of all these cases, we choose two to look at in detail. Our choice is guided by the predicted stability diagrams from DK, which indicate the possibility of interesting coexistence regions of different states, e.g., coexistence of the trivial and a patterned steady state (SS), coexistence of SS and standing wave (SW), etc. Cases of particular interest are IV− with A > 0 and I− with A < 0 in DK and marked in blue in Table 1.

Localized SS
The first case that we discuss is the one labelled as case IV− with A > 0 in DK. Figure 4 small unstable branch of SS that lies in the third quadrant (where the trivial state is stable). Further, we observe a stable branch of SW in the fourth quadrant along with the unstable branch of SS. We identify this bifurcation scenario to be interesting as we expect the possibility for localised steady states in a background of trivial state (LSS-TS) in the third quadrant and the possibility for localised steady state in a background of standing waves (LSS-SW) in the fourth quadrant.
Figure 4(b) shows the predicted bifurcation diagrams (in black lines) from the TB normal form equation. The two cases correspond to the cases of variation of parameters along the diagonal above and below the line µ = ν respectively. The bifurcation above the diagonal in (ν, µ)-plane (see Figure  4(b)(i) ) has only a subcritical SS branch from a pitchfork bifurcation at µ = 0. The trivial state is stable in the region where µ < 0 and ν < 0. The bifurcation below the diagonal in (ν, µ)-plane (see Fig. 4 b, right panel) has a subcritical SS branch from a pitchfork bifurcation at µ = 0. Stable SW and unstable TW bifurcate from the trivial state at a Hopf bifurcation where ν = 0, µ < 0. The stable SW branch terminates on the subcritical SS branch with the formation of a heteroclinic orbit at a global bifurcation SL s connecting two saddles (the names of these bifurcations follows DK). The unstable TW branch terminates at the subcritical SS at L m . This scenario has been investigated analytically and numerically in thermosolutal convection (Knobloch & Proctor, 1981;Huppert & Moore, 1976;Da Costa et al., 1981;Knobloch et al., 1986) and is important because it was an early example of the discovery and analysis of how a Shil'nikov heteroclinic orbit (Knobloch et al., 1986) can lead to chaotic dynamics, though this is in a different parameter regime from that which we will investigate. To the bifurcation diagrams in Figure 4(b), we add predictions (from global stability requirements) of a large amplitude stable branch of steady patterned state (in red lines). This figure now illustrates the possible coexistence regions that could allow for the localised states detailed in the previous paragraph.
In order to confirm the existence of a stable large amplitude branch of patterned state (SS), we first SPATIAL LOCALISATION BEYOND STEADY STATES

of 26
Case Coefficients of the nonlinearities in DK look at asymptotic states accessible via time stepping. We treat the model numerically by discretizing the PDE both in time and space. In space, we discretize the model using spectral methods (Cox & Matthews, 2002;Kassam & Trefethen, 2005;Canuto et al., 2012;Hussaini & Zang, 1987) and fast Fourier transform (FFT) with 16 grid points per wavelength. In time, we discretize the model using the exponential time differencing method (ETD) (Cox & Matthews, 2002;Kassam & Trefethen, 2005;Canuto et al., 2012;Hussaini & Zang, 1987). The ETD method solves the linear parts of the PDEs exactly followed by a second order approximation of the nonlinear parts. We compare the stability region from solving the model with the stability region obtained from the normal form (Dangelmayr & Knobloch, 1987). To do this we solve the PDE and plot the type of solutions in (ν, µ)-plane, using polar coordinates defined as ν = r cos(θ ), µ = r sin(θ ), (6.1) where r is the radius that controls how far ν and µ are from the TB point, and θ is the angle that controls the position of ν and µ in the (ν, µ)-plane. Note that the Hopf bifurcation occurs at ν = 0 and µ < 0 which correspond to θ = 270 • . The pitchfork bifurcation occurs at µ = 0 which correspond to θ = 0 • and θ = 180 • . Figure 4(c) shows the results of starting time stepping from different initial conditions for a variety of system parameters at two different radii r = 0.7 and r = 0.9. At radius r = 0.7, we start from initial conditions close to a pitchfork bifurcation at θ = 180 • . We are able to obtain large amplitude  (Dangelmayr & Knobloch, 1987), where the panel (i) represents the bifurcation above the diagonal in the (ν, µ)-plane and the panel (ii) represents the bifurcation below the diagonal in the (ν, µ)-plane. In this bifurcation diagrams, the black lines (solid for stable solutions, dashed for unstable solutions) are from DK, while the red lines are inferred from our system given Lyapunov stability. (c) Plot of the solution from solving the model (3.12) by time stepping with parameters values Q 1 = 0.9, Q 2 = −0.2,C 1 = −0.2,C 2 = C 3 = −1,C 4 = −0.5,C 5 = 6 and b = 2 for radius r = 0.7, 0.9. A Hopf bifurcation occurs at θ = 270 • and a pitchfork bifurcation occurs at θ = 180 • and θ = 0 • . The red x and green square refer to extended SS and SW solutions, respectively. The half line SL s is the line of heteroclinic bifurcations where SW joins the small-amplitude unstable SS. SS solutions (shown as red crosses) as well as recover the trivial state (not shown with markers). The trivial state is stable when µ < 0 and ν < 0 until we reach a Hopf bifurcation close to θ = 270 • . At this bifurcation, the trivial state loses stability and a new branch of SW are created (shown as green crosses). The amplitude of the SW branch increases with increasing 270 • < θ < 297 • which is close to the prediction of an SL s bifurcation at an angle of θ = 308 • from the normal form analysis. The complete circle of red crosses observed for the large amplitude SS in Fig. 4(c) at r = 0.7 indicates that the solution branch of large amplitude SS is disconnected from the low amplitude SS solution branch when viewed as a function of θ .
At a slightly larger radius with r = 0.9, time stepping identifies similar coexistence of trivial state and large amplitude SS solutions in the third quadrant. However, the large amplitude SS solutions exist only till θ = 245 • . At θ = 270 • , we encounter the Hopf bifurcation, as before, resulting in the branch of SW solutions. The branch of SW exists in the range of 270 • < θ < 309 • and terminate close to the prediction of the SL s bifurcation. We find that we are able to recover the large amplitude SS branch again from θ = 297 • . The fact that we are able to identify two θ values where the large amplitude SS solution branch terminates indicates that they are the locations of saddle-node bifurcations where the large amplitude SS branch connects with the low amplitude SS branch. Figure 4(c) has identified two bistable regions: bistability between the trivial state and large-amplitude SS when µ < 0, ν < 0 and bistability between small-amplitude SW and large-amplitude SS in the region between ν = 0, µ < 0 and the half line SL s . We now look to obtain localized states in these regions. To do this we follow the process. First, we increase the domain size to allow 64 wavelengths and perform time stepping to find an extended SS. Then, we use a sech-envelope with different widths to construct several initial conditions and perform time stepping again to obtain nearby dynamically stable localized states. Using this method we are able to get localized steady states, LSS in both of the bistable regions with two different backgrounds. First, we find LSS-TS, which has a localised steady state with a trivial state (TS) as a background in the region where the trivial state and large-amplitude SS are both stable (µ < 0 and ν < 0). Figure 5 (a) shows one example of LSS with TS background for radius 0.7 and θ = 200 • (ν = −0.54, µ = −0.45). There are other LSS-TS with different widths, depending on the initial conditions, and to investigate these in detail we perform numerical continuation in next section.
Second, we find LSS with an MW background in the region where the large-amplitude SS and the small-amplitude SW are both stable. The bistability occurs in the region between the Hopf bifurcation at θ = 270 • (ν = 0, µ < 0) to the half line SL s at θ ≈ 308.4 • , as mentioned above. The small-amplitude MW background solutions move as a function of time. Initially, the small-amplitude solutions are SW with large-amplitude SS in the middle of the domain. As time increases the SW turn in to MW. This change occurs because the left-right symmetry of the SW solutions is broken by the SS solutions in the middle. Figure 5 (b) shows one example of LSS with MW background for radius 0.7 and θ = 280 • (ν = 0.12, µ = −0.69). Many widths of this class of solutions can be obtained by altering the initial width of the sech-envelope.
Unlike the LSS with the trivial state background, in this case two patterns (large-amplitude SS and small-amplitude MW) coexist, which suggests the presence of a spatial heteroclinic orbit between the SS and MW states.

Numerical continuation
In the following we use numerical continuation to compute steady numerical solutions of the model for both the extended SS and the LSS. The method we use is based on Newton iteration and pseudo arclength continuation (Doedel et al., 1991). Seting the time-derivative terms to zero in (3.12) results in the steady Swift-Hohenberg equation: We note that even when the solutions depend only on µ, the stability depends on both µ and ν. The initial guesses for the branches of extended SS and LSS are obtained from time stepping. Since θ controls the values of µ and ν in the (ν, µ)-plane, we will use θ as the bifurcation parameter with fixed radius. We show the bifurcation diagram of all solutions as functions of θ and µ, where µ = r sin(θ ), against the norm as a measure of the solutions, where We will perform continuation at two different radii r = 0.7 and r = 0.9 in a one dimensional domain allowing 16 wavelengths with 16 grid points per wavelength.
6.1.1 Continuation for radius 0.9 Starting from initial guesses obtained from time stepping in numerical continuation with radius r = 0.9, we perform continuation to obtain the extended SS and the LSS-TS solutions. Figure 6 shows the solutions from numerical continuation with θ as the control parameter in (a) for radius r = 0.9 and (b) for radius r = 0.7. In Fig. 6, we represent the full branch of extended SS in black, the LSS branch with even numbers of peaks L e in orange and the LSS branch with odd numbers of peaks L o in green. Along the odd branch L o the midpoint (x = Lx/2) of the localized state is always a global maximum (with an odd number of maxima), while along the even branch L e the midpoint (x = Lx/2) is a global minimum (with an even number of maxima). Figure 7 shows Fig. 6 (a) for µ where µ = d sin(θ ) in more detail with examples of LSS with different widths. The extended SS branch emerges subcritically from the trivial state at θ = 180 • (the pitchfork bifurcation µ = 0) and undergoes a saddle-node bifurcation at θ = 245.6 • (µ = −0.81). The branch changes at the saddle-node to a large-amplitude stable state. It reaches the maximum amplitude when θ = 90 • (µ = 0.9) and decreases until reaches to the second saddle-node at θ = −65.6 • (µ = −0.81). The brach then decrease further until terminates back to a pitchfork bifurcation θ = 0 • (µ = 0). Because the solutions depend on µ only there is a symmetry θ → 180 • − θ .
Similarly, we use time stepping to obtain the initial guess for the even L e and odd L o branches for localized state and then do the continuation. Both branches emerge subcritically close to the pitchfork bifurcation with small-amplitude and undergo a series of saddle-node bifurcations producing the homoclinic snaking. Each branch adds an oscillation on each side as it snakes back and forth, until they reach the width of the domain where they terminate on the SS branch, at the saddle-node point. The snaking region occurs between θ = 209.5 • (µ = −0.44) and θ = 240.5 • (µ = −0.78). The Hopf bifurcation occurs when θ = 270 • (where ν = 0 and µ = −0.9). Note that, for this choice of radius the homoclinic region occurs away from the Hopf bifurcation. The trivial solution is stable in the region where µ < 0 and ν < 0 and unstable anywhere else.
6.1.2 Continuation for radius 0.7 In this section, we will examine the numerical continuation for smaller radius r = 0.7. The difference here from the previous case that we discussed for r = 0.9, is that here the saddle-node bifurcation of the SS branch at µ = −0.81 is beyond the minimum value of µ possible, i.e., µ = −0.7. This implies that the left side of the snaking branch cannot be reached. To obtain the numerical branch for LSS as well as the extended SS branch we use the same process as for radius 0.9. We show the results in Fig. 6 (b) for θ against the norm.
Starting from large extended SS as initial guesses in continuation, we find that the branch of extended SS has large-amplitude with maximum at θ = 90 • (µ = 0.7) and minimum amplitude at θ = 270 • (µ = −0.7). This branch continues around the circle from θ = 0 • to θ = 360 • and does not connect to the pitchfork bifurcation, which indicates that there is no saddle-node bifurcation.
To find the localized branch we start continuation from initial guesses with various numbers of peaks obtained from time stepping. The continuation shows that the odd and even branches with at least 15 peaks form isolas. They start from the region where the trivial state stable (µ < 0 and ν < 0) and continue to the region where the trivial state not stable after passing the Hopf bifurcation (ν = 0 and µ < 0). Each isolated branch reach the minimum amplitude at θ = 270 • (µ = −0.7). The snaking region of isolated branch occurs between θ = 219.5 • (µ = −0.44) and θ = 320.5 • (µ = −0.44). The two lower branches with one and two peaks bifurcate close the pitchfork bifurcation θ = 180 • (µ = 0) with smallamplitude, reach maximum at θ = 270 • (µ = −0.7), and then terminate with small-amplitude close to the pitchfork bifurcation θ = 360 • (µ = 0).
From time stepping, the LSS with the trivial state background are stable in the region where µ < 0 and ν < 0 and become unstable after passing the Hopf bifurcation in the region where µ < 0 and ν > 0. After the Hopf bifurcation, LSS-TS exist, but the stable solutions are the LSS with an MW background where there is bistability between SS and SW. Figure 8 summarizes the results obtained from numerical continuation for radius r = 0.7 and r = 0.9 in the (ν, µ)-plane. The black curve refers to the extended SS and the blue curve refers to the region where LSS with the trivial state background exist at both radii. The region between the red dashed lines identifies the snaking region for radius r = 0.9 and the isolas for radius r = 0.7. As we take smaller radius the snaking region becomes narrower. This also indicates that for any radius r < 0.44 no stable LSS can occur since the bistability between two stable states lies beyond the snaking region. In Fig. 8 18  8: (ν, µ)-plane from solving the PDE using numerical continuation for radii 0.9 and 0.7. The pink and red cures refer to the pitchfork bifurcation and Hopf bifurcation, respectively. The black curve refers to the extended SS solutions and the blue curve refers to the snaking region for radius r = 0.9, as shown in Fig. 7 (a) and the region of isolas for radius r = 0.7 as shown in Fig. 7(b). At r = 0.9, this homoclinic snaking occurs between µ = −0.44, ν = −0.78 and µ = −0.78, ν = −0.44. The SN point occurs at µ = −0.81 for all radii. Red dashed lines mark the snaking region as shown previously in Fig. 7(a). The localized solutions stable in the region where the trivial states stable (µ < 0, and ν < 0) and unstable in the region where the trivial states become unstable (µ < 0, and ν > 0). The green star markers refer to the LSS with an MW background and are obtained from time stepping. The SL s half line is the line where the bifurcation changes from SW to SS in the normal form (see Fig. 4).
we also show points obtained from time stepping where the LSS with an MW background exist as green stars. The SL s half line is the line where the branch ends on a heteroclinic bifurcation (see Fig. 4). Beyond this line, there are no SW, and time stepping evolves to a large-amplitude SS.

Localised TW
The second case we discuss is the one labelled as case I− with A < 0 in DK where D > 0 and M < 0. We consider two bifurcation scenarios to illustrate these cases, one above and one below the line µ = ν and plot them in Fig. 9 (b). The bifurcation below the diagonal in (ν, µ)-plane (see Fig. 9 (b) (ii)) has an unstable SS branch bifurcating from a pitchfork bifurcation at µ = 0 with ν > 0. It also has a subcritical TW branch and unstable SW branch bifurcating from the trivial state at a Hopf bifurcation where ν = 0, µ < 0. The SW branch undergoes saddle-node (SN) bifurcation at SN s2 and terminates at the SS branch at L M . From global stability, we expect the unstable branch of TW to regain stability at a saddle-node bifurcation, giving rise to a large amplitude branch of TW solutions (as shown by red lines in Fig. 9 (b). The bifurcation above the diagonal in (ν, µ)-plane (see Fig. 9 (b)(i)) has stable SS branch bifurcating from a pitchfork bifurcation at µ = 0 with ν > 0 which becomes unstable after passing the half lines L m . The subcritical TW branch bifurcates from the SS branch at L m . The trivial state is stable in the region where µ < 0 and ν < 0. This implies that we can expect coexistence between the large amplitude TW and the trivial state for values before the pitchfork bifurcation and we can expect coexistence between large amplitude TW and stable SS solutions in the range of values past the pitchfork bifurcation.
In order to explore the fully nonlinear behaviour of the system, we run timestepping from different initial conditions for a variety of parameters at three different radii r = 0.1, 0.7, 0.9 in a small (one wavelength) domain and plot the results in Fig. 9(c). In these calculations, we use 32 grid points per wavelength.
The trivial state is stable for the region ν < 0 and µ < 0 for all radii r = 0.1, 0.7, 0.9. The smallamplitude SS bifurcate at a pitchfork bifurcation θ = 180 • where µ = 0, ν < 0 and loses stability at L m , which has the slope µ ν ≈ −0.502. There is a large-amplitude TW solution around the whole circle in the (ν, µ)-plane at radius r = 0.1. At radii r = 0.7 and r = 0.9, the amplitude of the TW decreases in the . The blue and red curves refer to u and u t , and the waves travel to the right. A movie of the state in (a) is available at (Alrihieli et al., 2020).
region where the trivial state or the small-amplitude SS are stable and we lose the branch solutions (at potential saddle-node bifurcations). This confirms that the fully nonlinear behaviour in this case allows for bistability between two different states: a large-amplitude TW with trivial states in the region where µ < 0, ν < 0 and a large-amplitude TW with small-amplitude SS in the region between the pitchfork bifurcation µ = 0, ν < 0 to the half line L m . Therefore, the LTW-TS and LTW-SS can be sought. In order to obtain the localized state we increase the domain size to allow 64 wavelengths in the domain and perform time stepping to find an extended TW. To obtain the localized state, we use a sechenvelope with different widths and do time stepping to obtain the localised state. Using this method we are able to get LTW with two different backgrounds as shown in Figs. 10 and 11.
First, we find LTW-TS in the region where µ < 0 and ν < 0. Figure 10  . The blue and red curves refer to u and u t and the wave travels to the right. A movie of the state in (a) is available at (Alrihieli et al., 2020).
Second, we find LTW-SS in the region between the pitchfork bifurcation at µ = 0 where ν < 0 to the half line L m . Figure 11 shows two examples of LTW with an SS background for (a) radius r = 0.1 and θ = 170 • where µ = 0.017, ν = −0.098 and (b) for radius r = 0.4 and θ = 160 • where µ = 0.14, ν = −0.038. The LTW move from left to right. Again, we find only LTW-SS with one chosen width in this case.
In all these LTW examples, we started with initial conditions with a wide variety of widths, but always found a localised solution with the same width, unlike in the LSS cases.

Conclusions
In this paper, we have developed a simple nonlinear pattern-forming PDE model that has a Takens-Bogdanov primary bifurcation. The model is based on the Swift-Hohenberg equation, which was originally derived to describe the effects of thermal fluctuations and the evolution of roll patterns close to the onset of Rayleigh-Bénard convection and later used as a model of pattern formation in many physical problems. The new model can be reduced further using weakly nonlinear theory to the Takens-Bogdanov normal form where the pitchfork and Hopf bifurcations coincide. The advantage of the model lies in the relative ease of investigating the nonlinear behaviour numerically and analytically, as compared to the full PDEs of double-diffusive convection. Alongside the numerical results, the model is important for helping to understand the bifurcation structure and the solution behaviour close the Takens-Bogdanov point in an extended system and in particular for investigating localized solutions.
We identified two types of localized states which have previously been observed in various systems. The first of these is the localized steady state with the trivial state as a background (LSS-TS), which was observed numerically in binary convection (Batiste et al., 2006;Mercader et al., 2009) and in thermosolutal convection (Beaume et al., 2011). The second localised state is that of a localized travelling wave with the trivial state as a background (LTW-TS), which was observed in binary convection (Niemela et al., 1990;Barten et al., 1991Barten et al., , 1995bZhao & Tian, 2015;Jung & Lücke, 2005;Surko et al., 1991;Watanabe et al., 2012;Kolodner, 1991a,c,b). We have further identified two new spatially localised states: that of a localised steady state in a modulated wave background (LSS-MW) and that of a localised travelling wave in a steady state background (LTW-SS).
To find localized states, we looked for subcritical pitchfork and Hopf bifurcations, since we expected that subsequent saddle-node bifurcations would lead to stable large-amplitude solutions coexisting with the stable trivial solutions, and possibly then to localized solutions. We identified a subcritical pitchfork bifurcation for the case IV− with A > 0 (see Fig. 4) and a subcritical Hopf bifurcation for the case I− with A < 0 (see Fig. 9). From solving the model numerically, we obtained different types of localized states. In case IV− with A > 0, we obtained LSS in the region where there is bistability between the trivial state and a branch of periodic steady states, with µ < 0 and ν < 0. We used numerical continuation of the PDE model (3.12) to compute the branches of localized states. The continuation method we used can only find steady solutions, so the model is effectively the steady Swift-Hohenberg equation with solutions depending only on µ -though the stabilities depend on both µ and ν. The solutions are associated with homoclinic connections to the trivial state, in the same manner as spatially localized solutions in the Swift-Hohenberg equation (Burke & Knobloch, 2006Burke & Dawes, 2012). The two localized branches with odd and even numbers of peaks add an oscillation on each side as they snake back and forth until they reach the width of the domain, where they terminate on the steady state branch, at the saddle-node bifurcation (see Fig. 6 a and 7 i). In the absence of the saddle-node bifurcation in the periodic state where the left side of the snaking branch cannot be reached, we found isolated branches of LSS (see Fig. 6 (b) and 7 (ii)). The localized solutions we obtained still exist but are unstable in the region where the trivial state becomes unstable, where µ < 0 and ν > 0.
From time-stepping, we also found LSS with an MW background in the region where the largeamplitude SS branch and the small-amplitude SW branch are both stable (see Fig. 4, 5 (b)).
In case I− with A < 0 and from time stepping, we found LTW with the trivial state background in the region where the trivial state and a large-amplitude branch of TW are stable, with µ < 0 and ν < 0 (see Fig. 9). We also found LTW with SS background in the region where the small-amplitude SS and large-amplitude TW are stable (see Fig. 9 and 11). For the given parameter values, the LTW we found all have the same width, regardless of initial conditions. In contrast, LSS exist with a wide range of widths, with different numbers of peaks. In future work we will investigate why we get uniquely selected widths of LTW.
LTW with trivial state background and with uniquely selected widths have also observed in experimental (Kolodner, 1994(Kolodner, , 1991cNiemela et al., 1990;Kolodner, 1991a) and numerical (Barten et al., 1991;Taraut et al., 2012;Barten et al., 1995b) studies of binary convection. These studies were not carried out close to the Takens-Bogdanov point, so our model does not directly apply here. Using continuation to compute the LTW solutions would need more effort due to the time and space dependence. The numerical continuation would then require additional unknown variables: the group velocity and