## Abstract

A Kuramoto–Sivashinsky equation in two space dimensions arising in thin film flows is considered on doubly periodic domains. In the absence of dispersive effects, this anisotropic equation admits chaotic solutions for sufficiently large length scales with fully two-dimensional profiles; the one-dimensional dynamics observed for thin domains are structurally unstable as the transverse length increases. We find that, independent of the domain size, the characteristic length scale of the profiles in the streamwise direction is about 10 space units, with that in the transverse direction being approximately three times larger. Numerical computations in the chaotic regime provide an estimate for the radius of the absorbing ball in

## 1. Introduction

The one-dimensional (1D) Kuramoto–Sivashinsky equation (KSE) is
*L*], and an *L*-periodic initial condition *u*(*x*,0)=*u*_{0}(*x*). Owing to the conservative nature of (1.1) and the presence of a Galilean invariance, attention may be restricted to zero-mean solutions. This equation, and variants in higher space dimensions, arise in the study of spatio-temporal organization in reaction–diffusion systems [1], the propagation of flame fronts [2–4] and thin film flows down a vertical plane [5]. Variants also arise in two-phase flows [6–8]. Moreover, it is found to emerge in numerous applications in physics, including plasma physics [9], ion sputtering [10] and chemical physics for the propagation of concentration waves [11,12]. A wide range of dynamical behaviours are observed depending on the length *L* of the periodic domain. Increasing *L* above 2*π* (below 2*π* all solutions converge uniformly to zero), steady states, travelling waves and time-periodic bursts are observed, with the onset of chaos for large enough *L* [13]. The scaling of the system energy with the length parameter can be quantified by considering the *L*-dependent radius of the absorbing ball in the space *L*-independent constant *C*, and some function *P*(*L*). An estimate of this form was first constructed using a Lyapunov function approach for odd solutions of (1.1), giving *P*(*L*)=*L*^{5/2} [14], a result which was later improved and generalized to non-parity solutions, implying the existence of a finite-dimensional global attractor [15]. After many intermediate developments [16–19], the most recent analytical improvement to this bound shows that (1.2) is satisfied for all solutions to (1.1) with *P*(*L*)=*L*^{q} for any *q*>5/6 [20,21]. Numerical work provides strong evidence that the optimal estimate for (1.2) is given by *P*(*L*)=*L*^{1/2} [22]; this was shown to be sharp for steady solutions of (1.1) using a dynamical systems approach by proving the stronger property of uniform boundedness of solutions independent of *L* [23] (this *L* in the chaotic regime, suggesting an invariant energy distribution in the thermodynamic limit as

In this paper, we present numerical results for the spatially periodic initial value problem for a KSE in two space dimensions over rectangles *Q*=[0,*L*_{1}]×[0,*L*_{2}], given by
*u*(*x*,*y*,0)=*u*_{0}(*x*,*y*) and dispersion parameter *δ*≥0. This was derived by Nepomnyashchy [27,28] with *δ*=0, and in general by Frenkel & Indireshkumar [29] and Topper & Kawahara [30] to describe the weakly nonlinear evolution of the interface of a thin film flow down a vertical plane (see [31] for a discussion of the derivations of this model for different fluid dynamical regimes). Without loss of generality, we can restrict our attention to zero-mean solutions because the spatial average of a solution to (1.3) is conserved and the equation is invariant under a Galilean transformation as in the 1D case. In the absence of dispersion, i.e. *δ*=0, equation (1.3) was studied analytically by Pinto [32,33] in the case of *L*_{1}=*L*_{2}=*L*. He proved global existence of solutions, the existence of a compact global attractor and analyticity of solutions. Using the Lyapunov function method, he obtained the estimate for the radius of the absorbing ball in *L*_{1} and *L*_{2}, we define the *u*. The non-dispersive problem was also considered numerically by Akrivis *et al.* [34] on a square domain. They found that
*δ*=0 in the space *L*_{1} and *L*_{2}. The result (1.7) implies that the solutions in the chaotic regime possess a finite energy density. We obtain a similar picture for the energy distribution of the Fourier modes as is found for the 1D KSE (1.1); a plateau of the energy for the low modes, rising to a peak and then decaying exponentially for the higher Fourier modes. The addition of the extra dimension in the dissipative fourth-order term of (1.3) produces an asymmetric energy distribution. By considering the decay of the Fourier spectra for large wavenumbers, we observe an increased spatial analyticity due to two-dimensionality of the solutions on domains that are not thin.

Next, we introduce dispersion to the problem, and consider how a small positive value of *δ* affects the dispersionless solution dynamics, energy distribution and the absorbing ball estimate (1.7). Dispersive effects are often included in the 1D KSE (1.1); Akrivis *et al.* [35] considered the addition of both third- and fifth-order dispersion by studying the Benney–Lin equation
*μ*=0 (which is of interest here), formal asymptotics show that for *δ* is increased (see fig. 4.2 in [35]). This laminarizing effect of dispersion in the 1D problem was additionally investigated by Chang *et al.* [39], where the authors showed that increasing dispersion diminishes the family of steady and travelling wave solutions—only KdV pulses remain for large enough *δ*, with a large basin of attraction. Gotoda *et al.* [40] studied the route of the full dynamics as dispersion is strengthened; they additionally estimated the critical value of *δ*≈0.2 (which appears to be independent of system length) where high-dimensional chaos crosses to low-dimensional chaos.

In this paper, we are interested in weak dispersive effects which do not regularize the chaotic dynamics; we study the effect of the fixed values of *δ*=0.01, 0.1, and 1 on the dynamics of the two-dimensional (2D) KSE (1.3). We provide numerical evidence that, given a fixed value of *δ*, the *δ*≫1) here, but briefly comment on known results. Travelling wave attractors of 2D solitary pulses are found, as observed by Toh *et al.* [41] and Indireshkumar & Frenkel [42], in analogy with the results in [35] for the 1D equation. Saprykin *et al.* [43,44] also studied this problem on an infinite domain, and provided a detailed analysis of the interaction between pulses. It can be shown with formal asymptotics (proved rigorously in [45]) that the solutions of (1.3) in this large dispersion limit converge to solutions of the Zakharov–Kuznetsov equation
*U* and *τ* are rescalings of *u* and *t*, respectively. This is a higher-dimensional KdV equation yielding 2D solitons whose stability has been studied analytically in [46,47].

A related equation of interest is the multi-dimensional KSE,
*et al.* [54] provided a comprehensive numerical study to complement this analytical work. They give an exhaustive picture of the dynamics present for varying domain dimensions.

The structure of the paper is as follows. In §2, we briefly discuss the numerical schemes and data analysis tools employed for our simulations. The computations of (1.3) with *δ*=0 follow in §3, and the dispersive case follows in §4. In §5, we discuss our results and future work.

## 2. Numerical methods and data analysis tools

Equation (1.3) is solved numerically using implicit–explicit backwards differentiation formulae (BDF) for the time discretization, and spectral methods in space. The BDF belong to the family of linearly implicit methods constructed and analysed by Akrivis & Crouzeix [55] for a class of nonlinear parabolic equations. It was shown by Akrivis *et al.* [34] that such numerical schemes are convergent, and also that they are efficient and unconditionally stable under various conditions on the linear and nonlinear terms of the problem. We do not go into further details of these schemes here, because their applicability for our problem has been checked in [34]. As we are considering (1.3) on rectangular periodic domains, the solution may be written in the form of a Fourier series
*u*, and *u* is real-valued, *k*_{1}|≤*M* and |*k*_{2}|≤*N*, corresponding to a discretization of the spatial domain *Q* into (2*M*+1)×(2*N*+1) equidistant points.

The linear dispersion relation for (1.3) is
*s*, as a function of *L*_{1} and *L*_{2}, the stability of the (*k*_{1},*k*_{2})-mode (where *s*) that the purely transverse modes (*k*_{1}=0) are always linearly stable. If *L*_{1}≤2*π* (implying *u* and integrating over *Q*, it can be easily shown that the solution decays to zero exponentially. The purely imaginary component of

We write the domain lengths *L*_{1},*L*_{2} in a canonical form, taking *L*_{1}=*L* and *L*_{2}=*L*^{α}. We take *L* to present a view of the chaotic dynamics in the global attractor for many aspect ratios. For aspect ratios with *α*≤0, the domains are thin, as either *L*_{1}≤1 or *L*_{2}≤1. In the former case, we have trivial behaviour with solutions decaying to zero as mentioned before, and in the latter case only the purely streamwise modes may be linearly unstable, thus the dynamics of solutions are expected to be 1D (this is confirmed by numerical simulations). We use small-amplitude random initial conditions with unstable low wavenumber modes for our numerical simulations. For *L*,*α*). Owing to the existence of a global attractor [32], the large-time behaviour is independent of initial condition. Note that (2.5) does not contain contributions from the purely transverse modes (the summation excludes modes with *k*_{1}=0). For large values of *L*_{2}, the lower transverse modes have very small exponential decay rates and would affect the streamwise and mixed-mode dynamics at large times; taking transverse modes in the initial condition would only extend the transient phase of the dynamics.

We average desired quantities over one solution orbit in order to obtain an average of that quantity over the entire global attractor (orbits are assumed to be dense in the attractor). Instead of computing an estimate for
*T*_{1}≪*T*_{2} are two large times. We require *T*_{1} to be large enough that the solution has reached the global attractor, and *T*_{2} to be large enough that *T*_{1}=1×10^{4} and *T*_{2}=2×10^{4}, which proved to be suitable. To study the equipartition and analyticity properties, we consider the time-averaged power spectrum of solutions, given by
*T*_{1},*T*_{2}] as done for *E*_{L,α} is related to

## 3. Computations in the absence of dispersion

We first proceed with a numerical study of (1.3) with *δ*=0 on large periodic domains. As noted earlier, for *α*≤0, we either obtain trivial dynamics or 1D solutions corresponding to solutions of the 1D KSE (1.1), so we focus on domains with *α*>0 (not thin). Figure 2 shows instantaneous interfacial profiles of solutions in the chaotic regime at time *T*_{2}=2×10^{4}. A variety of aspect ratios are used: in figure 2*a* the domain is longer in the streamwise direction and has *L*_{1}=166.8, *L*_{2}=59.9 (i.e. *α*=0.8); figure 2*b* shows a solution on a square domain with *L*_{1}=*L*_{2}=122.5; and the domain in figure 2*c* is longer in the spanwise direction with *L*_{1}=46.4, *L*_{2}=215.4 (here *α*=1.4). In all cases shown, activity in the mixed modes promotes fully 2D solutions. These profiles highlight distinct features of solutions to (1.3) on sufficiently large domains; the behaviour is dominated by the streamwise dynamics, with solutions varying weakly in *y*, but maintaining the characteristic cellular behaviour in the *x*-direction associated with solutions to the 1D equation (1.1)—a streamwise slice of the solution profile is very similar to the typical profiles observed in the chaotic regime for the 1D equation. Movie S1 in the electronic supplementary material presents the time evolution of solutions to (1.3) for these aspect ratios; the solution profiles at the final time are those shown in figure 2. For all profiles shown in figure 2, the characteristic length of the nonlinear cellular structures in the streamwise direction is about 10 units; this corresponds to the most active streamwise Fourier mode which has a wavenumber slightly smaller than the most linearly unstable streamwise mode,

### (a) Computational estimation of the radius of the absorbing ball

In what follows, we present the results of extensive numerical experiments that were used to obtain an optimal numerical bound on the radius of the absorbing ball in *α*=1) in [34]. To obtain the results that follow, *α* was fixed to take values in the interval 0.6≤*α*≤2, and *L* increased to cover a sufficiently large range of domains that support complex dynamics (recall that the rectangular domain has dimensions *L*×*L*^{α}). For a given *α*, computations were carried out and the time-averaged quantity *L*. We do not consider *α*≤0 because the dynamics are 1D as noted earlier. The variation of *a* considers *α*≥1, i.e. domains that are longer in the spanwise direction, and figure 3*b* corresponds to *α*≤1, giving domains that are longer in the streamwise direction, with *α*=1 providing a reference between the two panels. We observe that a regime of direct proportionality between

A quantification of the linear behaviour apparent in figure 3 was carried out using a least-squares approximation to the slope of the different lines and their intercepts with the vertical axis. The slopes are found to be *α*+1 with an accuracy of 0.02 or less, and the vertical intercepts are all found to be zero, also with an accuracy of approximately 0.02. Hence, on sufficiently large domains with *α*>0, we observe that
*E*_{L,α}≈*L*^{1+α}. Surprisingly, the unit constant of proportionality in this expression for *E*_{L,α} is not found in the case of the 1D problem (1.1) where it is approximately 1.7. Recalling the definitions of *L*_{1} and *L*_{2}, and with the established proportionality of *E*_{L,α} with the quantity *C* is a constant which is independent of the length parameters. This result is also valid for *α*≤0, given the previous discussions on the dynamics in this regime. In fact, an even stronger result than (3.2) appears to hold; the numerical results provide evidence that the *u*| over *Q* at a fixed time) of solutions in the chaotic regime is bounded above by a constant, in direct analogy with the numerical results for the 1D equation (1.1)—this can be seen in figure 2, where the solution amplitude appears to be independent of aspect ratio and length parameters. We find that the mean of the *T*_{1} and *T*_{2} is approximately 2.4, with the maximum value often being as large as 3.5; this result appears to be independent of *Q* as long as the domain is sufficiently large. This computational evidence that *u* is bounded by an *O*(1) constant (for example, 4 would suffice) over all choices of *Q* trivially implies (3.2).

### (b) Equipartition of energy and analyticity of solutions

In the previous subsection, we presented numerical evidence that predicts how the time-averaged energy of the system scales with the underlying lengths *L*_{1} and *L*_{2} for large domains supporting chaotic solutions. It is also particularly interesting to understand how this energy is spread among the Fourier modes. Recall that for the 1D KSE (1.1), it was observed that the energy was equipartitioned, i.e. spread equally among the lower Fourier modes (see [22] for example). The energy distribution rises to a peak for the most active mode (this is slightly less than the most linearly unstable mode) and then decays exponentially after an inertial range where the power spectrum behaves like *u*=*v*_{x}. This power law behaviour has been attributed to the balance of the destabilizing and dissipative linear terms for *O*(1) wavenumbers [24]. By the Paley–Wiener–Schwartz theorem (see [56] for example), the exponential decay of the high-frequency modes informs us of the spatial analyticity properties of solutions. For the 1D equation on an *L*-periodic domain, it is observed numerically that
*β*(*L*) converges to approximately 3.5 as *u* analytically about the real axis into the complex plane in a strip with |*Im* *x*|<*β*(*L*). It is noted that *β*(*L*) converges to 3.5 from above (meaning that solutions lose spatial analyticity as *L* increases), and the limit value can thus be surmised to be the optimal lower bound of the width of the analytic extension.

For completeness and to check our numerical work, we have recovered the above results for (1.3) in the special limit that our periodic domain is thin in the transverse dimension, and we again concentrate on numerical results for aspect ratios with *α*>0. The key quantity is the time-averaged power spectrum *T*_{1},*T*_{2}] as done for the energy in (2.8). Figure 4 depicts the spectrum *Q*|=10^{4}. The values of *α* used in figure 4 are 0.8, 1 and 1.4 for parts (*a*)–(*c*), respectively, and the corresponding values of *L* are 10^{20/9}, 10^{2} and 10^{5/3} (recall that |*Q*|=*L*^{1+α}). The respective streamwise–spanwise aspect ratios are 10^{4/9}≈2.7826, 1 and 10^{−2/3}≈0.2154. The figure shows logarithmic (base 10) contour plots for the three cases in the positive wavenumber quadrant corresponding to *k*_{1},0)-modes (the spacing between the contours is the largest in this direction). Secondly, the power spectrum remains *O*(1) as *Q* used to produce the figure). In figure 5, we zoom in on the low mode region of the power spectrum shown in figure 4*c*; the other cases provide similar plots. The tongue of most active modes (depicted by the white region in the figure) is consistent with the characteristic length scales of the cellular structures of profiles shown in figure 2. The most active streamwise mode with wavenumber 0.6 gives a length of approximately 2*π*/0.6≈10 space units. The longer length scale in the transverse direction is compatible with the observation that the tongue only extends to mixed modes with transverse wavenumbers around 0.3.

Figure 6 provides a better description of the behaviour of the low modes, plotting the power spectrum against the size of the scaled wavenumber vector, *a* compares three different aspect ratios, with two sets of data for the square domain case—all simulations exhibit fully 2D chaotic dynamics. The purely streamwise modes are interpolated with a cubic spline which appears to bound the data points; these modes carry the most energy, which is unsurprising given the anisotropy of figure 4, and the fact that they are the most linearly unstable modes for a given value of *a*; we propose that the disappearance of the inertial range is due to the mixed mode activity when the transverse length is sufficiently large. This is investigated in figure 6*b*, where for a fixed value of *L*_{1}=250, we observe how increasing *L*_{2} affects the interpolant through the streamwise mode data points. For *L*_{2}=1, 10, the mixed modes are not active in the solutions and the resulting dynamics are 1D—the solutions are just elongations of solutions to the 1D KSE (1.1) in the transverse direction. The dotted line in figure 6*b* corresponding to *L*_{2}=1 matches the curve in [22] (for *L*_{2}=1, the definition of the power spectrum (2.9) reduces to the definition for the 1D case), and the curve for *L*_{2}=10 is simply a factor of 10 greater. Increasing *L*_{2} further, the spectrum begins to widen with increased activity in the mixed modes, and the streamwise component of the power spectrum tends towards the solid line shown in figure 6*b* for *L*_{2}=100. For this choice of *L*_{2}, the mixed modes are fully active in solutions, although we omit the data points lying on and below the interpolant in figure 6*b* because they are shown in figure 6*a*. Note also that no mixed modes are linearly unstable until approximately *L*_{2}=40, although activity is seen for much smaller *L*_{2} due to the energy transfer through the nonlinear term; equivalently, the mixed modes are linearly unstable about 1D chaotic solutions for much smaller *L*_{2}—this can be observed from a crude truncation of the set of ODEs for the Fourier modes (2.4). The inertial range (the linear behaviour for wavenumbers beyond the most active wavenumber) visible for *L*_{2}=1, 10, is no longer discernible for *L*_{2}=100, and the most active mode shifts even further towards the longer waves. This is consistent with the finding that the characteristic streamwise cell size of the profiles in figure 2 is larger than that found in simulations of the 1D equation, the respective values being 10 and 9 space units. Note that the equipartition observed in figure 6 is expected given its relation to the solution energy (2.10) which scales with |*Q*|—there is a constant energy density of solutions in the large domain limit.

The effect of the mixed mode activity can be seen more drastically in the analyticity of solutions. In [57], the authors give the generalization of the connection between the decay rate of the Fourier spectrum and the analyticity of solutions to the 2D case. Informally, the observation that
*u* may be extended holomorphically into *β*(*L*_{1},*L*_{2}). They also provide the analytical estimate for the problem (1.3) with *δ*=0 and *α*=1 that
*β* may be computed numerically using a least-squares approximation from the slope of *L* as mentioned earlier, and we find the same result in the 2D case, contrasting the analytical result (3.5). We recover the convergence of *β*(*L*_{1},*L*_{2}) to approximately 3.5 for thin domains in the transverse dimension—the first two rows of table 1 take lengths *L*_{1} and *L*_{2} which result in no mixed mode activity, hence the dynamics are that of the 1D equation. We observe that increasing *L*_{2} so that mixed modes are active in the chaotic solutions actually improves the radius of analyticity, as observed in rows 3 to 7 of the table. For large domains with solutions exhibiting fully 2D spatio-temporal chaos, we are able to estimate that solutions can be extended holomorphically into *β* different from those obtained in [34] for the case of *α*=1, where the analyticity of solutions is found to be less than that observed for the 1D equation. Indeed, we find an increase in *β* from the 1D value, something which would be expected given the additional dissipation. Obviously, this does not improve the optimal lower bound on the analytic extensibility of solutions in the attractor, but it tells us that increasing two-dimensionality improves analyticity of the solutions (we assume that this is due to the activity of the mixed modes promoting energy in the dissipative range to move away from the purely streamwise modes).

The 2D KSE (1.11) studied in [54] is symmetric, and the resulting power spectrum is thus a function of

It is important to consider the possibility of a regime of dynamics beyond the length scales studied in this paper as discussed for the 1D problem (1.1) in [22]. It has been observed through extensive numerics and analysis that the large wavelength fluctuations of the 1D form of (1.11) (*L*, and is related to the *L*^{1/2} and the *L*. This scaling is observed for relatively small system sizes, although the other two critical exponents^{1} characterizing the KPZ universality class are not observed until *L* is much larger when full crossover to the KPZ scaling occurs. It is also worth noting that this asymptotic description is consistent with the observed energy spectrum. With this knowledge of the dynamics for the 1D problem, we conjecture that the energy behaviour (3.2) will not exhibit a crossover to a different scaling for even larger periodic domains. We do not attempt to compute the growth and dynamic exponents in this study, nor do we believe that the domain lengths used here are large enough to estimate these successfully; a number of studies have attempted to calculate these exponents for similar KS-type problems, but do so by resorting to very crude numerical discretizations in order to compute at large system sizes for a large number of time units. We are not certain that the form of the spectra observed for solutions in which mixed modes are active (figure 6) does not enter a different scaling regime which is computationally out of reach. In one of the less extreme cases used to compute the solution on a square domain with side *L*=100, there are 386 linearly unstable modes in total. This requires a numerical truncation with at least *M*=400, *N*=200 (approximately 320 000 Fourier modes) to obtain good accuracy (the spectrum is resolved to machine accuracy). Combining this with small time step requirements and large times of integration requires a large computing time.

## 4. Computations when dispersion is present

For the 1D KSE (1.1), it was observed in [35] that the strengthening of a physically derived third-order dispersion term can lead to a reverse period-doubling cascade. It is suggested that sufficiently large *δ* (i.e. a large amount of dispersion) can regularize chaotic dynamics for any system length *L*, as solutions are observed to be attracted to nonlinear travelling waves. Surprisingly, third-order dispersion acts as a destabilizing mechanism for this equation, competing with the stabilizing nonlinear term—it hinders the transfer of energy from low to high wavenumbers, and consequently analyticity of solutions reduces as *δ* is increased (see [61] for a discussion of this for the 1D case). Turning to the 2D problem, Toh *et al.* [41] and Indireshkumar & Frenkel [42] observed pulse solutions of (1.3) for large values of dispersion on large periodic domains—the usual streamwise cellular structures are found to be unstable and give way to the 2D pulses. An example of such a multi-pulse solution is given in Movie S2 in the electronic supplementary material, where the *O*(*δ*) pulses are seen to form an arrow head arrangement (the parameters taken are *L*_{1}=*L*_{2}=100 and *δ*=25). The arrow head of solitons is approximately time-periodic, with a period of about 10 time units; the pulses travel in the positive *x*-direction (streamwise) above a chaotic sea state of waves travelling upstream. Chaotic fluctuations of *O*(1) still exist in this case, but temporally periodic solutions are observed when *δ* is larger, where the *O*(1) component of solutions are time-periodic interactions at the bases of the pulses (see fig. 5 in [41]). We do not investigate the large *δ* limit here, nor questions concerning the regularization of chaotic dynamics. We are concerned with weak dispersive effects which do not fully regularize the chaotic behaviour, and observe how this affects the absorbing ball estimate and equipartition in the previous section.

The addition of dispersion qualitatively changes the profiles of solutions observed in the chaotic regime (see figure 2 for *δ*=0), yet they remain dominated by the streamwise dynamics as long as *δ* is not too large. The profile of a numerical solution in the chaotic attractor for *δ*=1 is shown in figure 7 for a square periodic domain with *L*=100. For *δ*=1, wavefronts are apparent, with higher peaks than the dispersionless case and flat trough regions in between. Streamwise slices of these profiles are similar to the solutions of 1D dispersive KS-type problems, for example the Benney–Lin equation (1.8)—the solutions observed are chaotic interactions of KdV pulses. These wavefronts cross and interact nonlinearly; this can be seen in Movie S3 in the electronic supplementary material, where the evolution of a solution to (1.3) with *δ*=1 is shown, and the profile at the final time is the same as that in figure 7. Our numerical simulations agree with the conjecture that chaotic dynamics may be regularized with sufficiently strong dispersion. We also note the path along which the regularization appears to occur as *δ* is increased: from the streamwise-dominated dynamics observed in the absence of dispersion, the cellular structures in the transverse direction begin to become more peaked in places forming wavefronts perpendicular to the streamwise direction. These fronts break up yielding pulse structures (for a square domain of side 100, this occurs around *δ*=5). Then, the arrangements of these 2D solitons are regularized fully to travelling waves for much larger *δ*.

In the dispersionless case, recall we observed that *E*_{L,α}≈*L*^{1+α}. To extend this result to the case of non-zero dispersion, we performed numerical simulations for *α*=0.8, 1 and 1.4, taking *δ*=0.01, 0.1 and 1. We obtained the same result as shown in figure 3, with a modification in the intercept of the straight lines with the vertical axis; this corresponds to the introduction of a constant *α*=1 with *δ*=0.01,0.1,1 is shown in figure 8*a*, and it is clear in this case, as in the other cases, that *δ* (the line for *δ*=0 is not included in figure 8*a* because it is unmistakeable from the *δ*=0.01 line at this scale). From our computations, we find roughly that *δ*∼*o*(1) to *δ*; the chaotic profiles for larger values of dispersion have larger amplitude solutions, and the dynamics appears to consist of the creation, interaction and annihilation of many 2D pulses. The scaling of the *δ* becomes linear in the large dispersion regime where the solution converges to travelling wave solutions of the ZKE (1.10), scaled by *δ*. Figure 8*b* shows how the increase of *δ* affects the energy distribution among the Fourier modes. The plot uses data from numerical simulations with *δ*=0.01, 0.1 and 1, for square domains of sides *L*=120, 130, 140, and shows the interpolants of the streamwise data points (these are found to be the most active modes as in the *δ*=0 case). Increasing *δ* results in a larger value of the small wavenumber asymptote and an increase in the energy in the low modes—this is consistent with the fact that *C*(*δ*) is an increasing function of its argument. The interpolant of the data points for *δ*=0.01 is almost identical to the dispersionless case shown in figure 6*a*, thus we do not plot the latter for comparison. For the moderate value of *δ*=0.1, the energy equipartition is skewed, as the active mode hump widens towards the low wavenumbers. The hump of active modes appears to cover the entire low wavenumber range for *δ*=1, and thus we recover the equipartition of energy among the low modes. For larger values of *δ* (for example, *δ*=25 as in Movie S2 in the electronic supplementary material), we recover peaks in the spectrum as observed for the 1D problem by Gotoda *et al.* [40]; however the mixed modes are much more active—further investigation of the dynamics of moderate to large *δ* is warranted.

In analogy with the 1D case, we see that the addition of dispersion decreases the radius of analyticity of solutions; for example, in the case of *δ*=1 and a domain which yields fully 2D solutions, it is observed that the Fourier coefficients decay as (3.4) with *β*≈3.5. As before, we find that the optimal numerical lower bound on the strip of analyticity occurs for thin domains (the smoothening of solutions due to two-dimensionality is independent of *δ*), and the corresponding 1D results are investigated in [35].

## 5. Conclusion

In this work, we have studied the dynamics of a physically derived dispersive KSE (1.3) in two spatial dimensions exhibiting extensive behaviour. Without dispersion, we observed that for sufficiently large domains, the system enters a regime of full spatio-temporal chaos, which is dominated by the streamwise dynamics (see Movie S1 in the electronic supplementary material). Furthermore, the system possesses a constant energy density since the *L*_{1} and *L*_{2} become large (figure 4) and the *Q*. These features are seen for the 1D KSE (1.1); however, the anisotropic 2D KSE (1.3) of interest in this paper does not present an inertial range in the simulations we have performed with mixed mode activity. In addition to this, we saw that the increase in two-dimensionality of solutions, through increasing the transverse length *L*_{2} until mixed modes become active, results in increased spatial analyticity. The optimal lower bound on the strip of analyticity is found when the domain is thin in the transverse direction, where the dynamics are governed by the 1D equation (1.1).

The addition of strong dispersion results in regularization of the spatio-temporal chaos, but moderate values of *δ* (dispersion parameter) change the nature of the chaotic dynamics, with interacting wavefronts that resemble KdV-type pulses emerging (see Movie S3 in the electronic supplementary material). The energy density is an increasing function of *δ*, and the constant *δ* fixes the pulse height, and the number of pulses scales with the size of the periodic rectangular domain. Much larger values of dispersion require smaller time steps for good accuracy, and a comprehensive study of the moderate to large *δ* regime for very large domains is numerically challenging.

It appears that finite energy density, corresponding to systems exhibiting equipartition, is a hallmark of the dynamics of KS-type systems with a *uu*_{x} nonlinearity. This property has been shown for multi-dimensional equations even with the addition of dispersion and variation in the linear and nonlinear terms. A non-local KSE in 2D was derived by Tomlin *et al.* [62] for the problem of a gravity-driven thin liquid film under the action of a normal electric field. Preliminary results appear to indicate a finite energy density for this problem also. Current work by the authors is the investigation of the extent of the class of PDE with quadratic nonlinearities exhibiting a finite energy density (corresponding to a roughness exponent of 0) by considering non-local variants of (1.1). Correspondingly, there is the related problem of finding the extent of the KPZ universality class by considering equations with the

## Data accessibility

An executable file, datafile and Matlab script required to run and analyse numerical simulations of (1.3) are available at: https://github.com/RubenJTomlin/Anisotropic-dispersive-2D-Kuramoto-Sivashinsky-Equation.

## Authors' contributions

All the authors contributed equally to this work.

## Competing interests

There are no competing interests.

## Funding

The work of D.T.P. was supported by EPSRC grant nos EP/L020564/1 and EP/K041134/1.

## Acknowledgements

R.J.T. acknowledges the support of a PhD scholarship from EPSRC and A.K. acknowledges funding by the Leverhulme Trust through an Early Career Fellowship.

## Footnotes

Electronic supplementary material is available online at http://dx.doi.org/10.6084/m9.figshare.c.4030363.

↵1 These are the growth and dynamic exponents which characterize the transient dynamics, i.e. before the solution orbit enters the chaotic attractor. These exponents are defined by how the surface roughness grows with time before saturation, and how the critical saturation time scales with the system lengths, respectively. Such exponents are not of current interest to us because we study large-time properties of solutions.

- Received September 27, 2017.
- Accepted February 27, 2018.

- © 2018 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.