Three-dimensional surface gravity waves of a broad bandwidth on deep water

Abstract A new nonlinear Schrödinger equation (NLSE) is presented for ocean surface waves. Earlier derivations of NLSEs that describe the evolution of deep-water waves have been limited to a narrow bandwidth, for which the bound waves at second order in wave steepness are described in leading-order approximations. This work generalizes these earlier works to allow for deep-water waves of a broad bandwidth with large directional spreading. The new NLSE permits simple numerical implementations and can be extended in a straightforward manner in order to account for waves on water of finite depth. For the description of second-order waves, this paper proposes a semianalytical approach that can provide accurate and computationally efficient predictions. With a leading-order approximation to the new NLSE, the instability region and energy growth rate of Stokes waves are investigated. Compared with the exact results based on McLean (J. Fluid Mech., vol. 511, 1982, p. 135), predictions by the new NLSE show better agreement than by Trulsen et al. (Phys. Fluids, vol. 12, 2000, pp. 2432–2437). With numerical implementations of the new NLSE, the effects of wave directionality are investigated by examining the evolution of a directionally spread focused wave group. A downward shift of the spectral peak is observed, owing to the asymmetry in the change rate of energy in a more complex manner than that for uniform Stokes waves. Rapid oblique energy transfers near the group at linear focus are observed, likely arising from the instability of uniform Stokes waves appearing in a narrow spectrum subject to oblique sideband disturbances.

essential in a wide range of fields, such as in engineering practices and for research purposes. A theoretical study has obvious advantages in elucidating the underlying physics, thereby advancing the understanding of realistic wave problems, compared with other approaches such as experiments, field studies and direct numerical computations. Theoretical findings have contributed to providing possible explanations of the formation mechanism of extremely large waves that appear suddenly with much larger amplitude than their surroundings, known also as 'rogue' or 'freak' waves. A few examples of the possible mechanisms proposed are modulational or Benjamin-Feir instability in deep water (Onorato et al. 2009), refraction by ambient currents or bathymetry (Janssen & Herbers 2009;Onorato, Proment & Toffoli 2011) and effects on weakly nonlinear waves of a depth change in shallow or intermediate water depth (Trulsen et al. 2020;Li et al. 2021a,b).
The modulational instability of weakly nonlinear Stokes waves subject to sideband disturbances was discovered by Benjamin & Feir (1967). Theoretical advances made in the late 1960s in this perspective include contributions by Lighthill (1965), Benjamin (1967), Whitham (1967), Zakharov (1968) and Benney & Roskes (1969). Since the discovery of the instability of a train of Stokes waves, progress has been made both numerically (Janssen 1983;Lo & Mei 1985;Trulsen & Dysthe 1997;Dysthe et al. 2003) and experimentally (Lake et al. 1977;Melville 1982;Su 1982;Chabchoub, Hoffmann & Akhmediev 2011) in understanding the evolution properties of surface waves on deep water. For smaller values of wave steepness, symmetrical upper and lower sidebands tend to grow equally in the initial stage. This is due to the degenerate resonant interaction in the prestigious 'figure of eight' quartet resonance loop (Phillips 1960(Phillips , 1967Longuet-Higgins 1976;Lake et al. 1977;McLean 1982), which corresponds to the special cases where two out of the four wavenumbers obeying the 'figure of eight' loop are identical. As nonlinearity increases, the sidebands appear to grow unequally, with the lower sideband growing faster than the upper sideband (Lake et al. 1977;Lo & Mei 1985), reaching a maximum larger than the minimum to which the spectrum peak drops. This unequal growth in energy, as a result of a combination of dissipation, wave breaking and nonlinear wave evolution, likely causes the downward shift of the spectral peak of wind waves (Lake et al. 1977;Lo & Mei 1985;Trulsen & Dysthe 1997).
The so-called nonlinear Schrödinger equation (NLSE) has been widely known as a convenient approach for analysing the instability of Stokes waves. For deep-water waves, Zakharov (1968) found that the wave envelope satisfies the NLSE which is cubic in wave slope (or wave steepness). The NLSE is derived based on the assumption of a small wave slope and the so-called narrow-banded assumption that restricts the modulation of the wave envelope to be slow in both space and time relative to a rapidly varying wave phase.
Many attempts have been made to expand the applicability scope of the NLSE for deep-water waves through adding higher-order terms, with a primary focus on reducing the restrictions due to bandwidth, as shown in table 1. Table 1 indicates the lower-order wave fields considered in a NLSE for the nonlinear and dynamic evolution of the (potential or elevation) envelope of the first harmonic. Specifically, the second-order wave fields do not appear alone in a NLSE but in combination with the first-order wave fields. Let δ x,y be a small non-dimensional parameter as a measure of the bandwidth in the main propagation (δ x ) and the transverse (δ y ) directions, respectively. With one further step, Dysthe (1979) has obtained a NLSE that is correct to O( 3 δ x,y , δ n x δ 3−n y ) (n = 0, 1, 2, 3), known as Dysthe's equation, or the fourth-order NLSE, as O(δ x,y ) ∼ O( ) is implied. With the fourth-order NLSE, Dysthe (1979) found that second-order mean flows play a significant Order in wave steepness Overall order of accuracy First order Second order Third order Dysthe (1979) (n = 0, 1, 2, 3) Trulsen & Dysthe (1996) O( δ n x δ 5−n y ) ) (n = 0, 1, . . . , 5) Trulsen et al. (2000) O Table 1. A summary of NLSEs for waves on deep water and the regimes of applicability which refer to the accuracy of lower-order wave fields considered in an NLSE for the dynamic evolution of the (potential or elevation) envelope of the first harmonic, in addition to the NLSEs for an arbitrary depth. Here denotes the dimensionless wave steepness; δ x and δ y denote the dimensionless bandwidth in the main propagation direction and in the direction normal to the main, respectively; δ ∞ x,y denotes arbitrary order in bandwidth.
role in the instability of Stokes waves. It is shown by Stiassnie (1984) that Dysthe's equation can be derived from Zakharov's integral equation (Zakharov 1968) through a narrow-banded assumption. Higher-order extensions in bandwidth are derived by Trulsen & Dysthe (1996) with a modified NLSE correct to O( δ n x δ 5−n y , 3 δ x,y ) (n = 0, 1, . . . , 5) and by Trulsen et al. (2000) in which the linear evolution of the envelope of linear wave is described by the exact linear dispersion relation, giving an NLSE correct to O( δ ∞ x,y , 3 δ x,y ). Using some of the earlier versions of NLSEs, Trulsen & Dysthe (1997) and Dysthe et al. (2003) found that the NLSEs can recover a downward shift of the spectral peak in a three-and two-dimensional narrow spectrum, respectively. The latter occurs after the spectrum has reached a quasi-steady state and hence, is unlikely due to the modulational instability which results in symmetrical growth of sidebands. Strong frequency dependence developed in the temporal evolution of directional waves initially of no frequency dependence is observed in NLSE-based numerical simulations and field data (Simanesew et al. 2016). Simanesew et al. (2017) found the modified NLSE of Trulsen et al. (2000) gives less accurate predictions for short crest waves with larger directional spread.
Many other approaches have been developed as alternative tools to the NLSE for understanding the properties of surface deep-water waves, such as the Zakharov integral equation and its leading-order approximations (see Zakharov 1968; Crawford, Saffman & 926 A34-3 Yuen 1980; Crawford et al. 1981;Krasitskii 1994), numerical wave tanks (e.g. Bateman, Swan & Taylor 2001;Ducrozet et al. 2012) based on the high-order spectral (HOS) method (Dommermuth & Yue 1987;West et al. 1987) and direct numerical solutions of fully nonlinear potential-flow equations (Engsig-Karup, Bingham & Lindberg 2009;Bihs et al. 2020;Zheng et al. 2020). Through studying the temporal evolution of three-dimensional steep focused wave groups, rapid energy changes in a wave spectrum are reported by earlier works e.g. Gibbs & Taylor (2005) and Gibson & Swan (2007). They attribute this phenomenon to the third-order resonant interactions, which are considered as the cause for the formation of the so-called 'wing waves' that appear in the transverse direction after the group at nonlinear focus reported in Barratt, Bingham & Adcock (2020). Direct numerical solutions were used by Barratt et al. (2020) for the evolution of a steep focused wave group with directional spread. In three dimensions, oblique energy transfer in a wave spectrum and a downward shift of the spectral peak are observed in previous works, e.g. Trulsen & Dysthe (1997) and Barratt et al. (2021).
There is still room and the need for further development of the models based on an NLSE, with two aspects addressed here. Firstly, ocean surface waves are generally three-dimensional and broad banded. Nevertheless, it should be noted that the NLSEs mentioned above, e.g. Dysthe (1979) and Trulsen et al. (2000), are only correct to order ( 3 δ x,y ), as shown in table 1. This suggests the need for extending the bandwidth at third order in wave steepness. Secondly, surface waves interact with ambient environments such as subsurface currents, turbulence and wind actions in the atmosphere, which often require explicit and accurate vertical structures of flow fields below the surface to allow for coupling. In contrast, it is known that a NLSE is based on the description of the envelope of the surface displacement, accompanied by mostly lowest-order approximations to the vertical profiles of flow fields at second order in wave steepness, as shown in table 1.
This paper aims to make attempts to fill in gaps in the aforementioned two aspects, through the development of a new framework that would allow for the study of three-dimensional waves of a broad bandwidth and would take into account the structures of flow fields below the surface without compromising much the computational efficiency of an NLSE-based model. The objective of this paper is threefold. First, a new NLSE for the evolution of three-dimensional deep-water surface waves is derived in § § 2-4 that does not rely on the assumption of a narrow bandwidth as aforementioned NLSEs. Specially, the newly developed NLSE should allow for δ x,y > 1 since it has relaxed the narrowband assumption for wave fields at second order, as shown in table 1. Simple numerical implementations of the new NLSE as performed with earlier versions of NLSEs are explained. Second, a semianalytical approach for the description of wave fields for bound waves at second order in wave steepness is proposed in § 3 which allows for simpler and more efficient numerical implementations than by Dalzell (1999). Finally, a study of the instability region and energy growth rate of Stokes waves and the temporal evolution of a directionally spread focused wave group are presented in § 6, where comparisons with Longuet-Higgins (1978), Dysthe (1979), Crawford et al. (1981), McLean (1982) and Trulsen et al. (2000) are made.

Problem definition
We consider ocean surface waves propagating on deep water in the framework of potential-flow theory, thereby assuming incompressible inviscid flows, irrotational fluid motions and negligible effects of surface tension. A Cartesian coordinate system is chosen with the undisturbed water surface located at z = 0. The system can be described as a boundary value problem governed by the Laplace equation, where Φ(x, z, t) denotes the velocity potential, ζ(x, t) is the free surface elevation, x is the position vector in the horizontal plane, t is the time and ∇ 3 = (∇, ∂ z ), with ∇ = (∂ x , ∂ y ) denoting the gradient in the horizontal plane. Equation (2.1) should be solved subject to the nonlinear kinematic and dynamic boundary conditions (cf. Davey & Stewartson 1974) at the free water surface z = ζ(x, t), as follows: and where the operator Γ is defined as Γ = ∂ tt + g∂ z with g the gravitational acceleration; a deep-water boundary condition is

Stokes expansion and separation of harmonics
In order to solve the boundary value problem (2.1)-(2.3), we seek the solutions of unknown Φ and ζ in a form of power series in wave steepness defined as = k 0 A 0 (a so-called Stokes expansion), with k 0 and A 0 denoting the characteristic wavenumber and wave amplitude, respectively, where we consider the first three orders and the superscripts denote the order in . Substituting (2.4) into the boundary value problem (2.1)-(2.3) leads to the decomposition of the fully nonlinear system into different problems through a collection of the terms at the same order in . The decomposed problems can be solved successively from the lowest to higher orders, as presented in the following from the first ( § 2.3) to the third order ( § 5). Let linear surface elevation ζ (1) be expressed in two equivalent forms, as follows: in which c.c. denotes the complex conjugates; A denotes the complex wave envelope of the surface displacement of the carrier wave, with k 0 = (k 0 , 0) the wavenumber vector chosen in the x direction; ω 0 denotes the angular frequency of the carrier wave that satisfies the dispersion relation ω 0 = ω(k 0 ) given by ω(k) = √ gk (where k = |k|); k = (k x , k y ) denotes a wavevector in the horizontal plane. Equation (2.5a) denotes the first-order elevation being expressed in an envelope-type form and (2.5b) a form through linear superposition of monochromatic waves in the Fourier k plane. This paper focuses on the 926 A34-5

Y. Li
former which leads to a nonlinear evolution equation for A that allows for the relaxation of a narrow banded assumption (cf. Chu & Mei 1971;Davey & Stewartson 1974).
Equating (2.5a) and (2.5b) we obtain a relation between A andζ , as follows: where the integral variable was replaced with κ κ κ through k = κ κ κ + k 0 . Introducing a Fourier transform for A defined as follows: (2.7b) A linear evolution equation for A for waves of a broad bandwidth is derived by Trulsen et al. (2000) in two equivalent forms, as follows: and (2.8c) where (2.8a) is used for constructing the solutions in § § 3, 4 and 6.2; (2.8c) is in a form convenient for numerical implementations, as explained in § 4.4. It is clear that a narrow-banded assumption implies k 0 LA O( ), which, again, this paper aims to drop. For convenience and later reference, we introduce a dimensionless bandwidth parameter defined, as follows: where δ x ( k x ) and δ y ( k y ) denote the dimensionless (dimensional) bandwidth in the x and y direction, respectively. The bandwidth is not assumed small, in contrast to conventional, reduced-form evolution equations for A by earlier works, e.g. Dysthe (1979) and Stiassnie (1984). Specifically, δ x,y > 1 are permitted in this paper. Similar to (2.5a), ζ and Φ are expressed in an envelope-type form through the separation of wave harmonics, as follows: in which the superscripts (ij) denote O( i ) and jth harmonic, are the modulated amplitude and potential, respectively, and the potentials at second and third harmonics are given by, respectively, It is worth noting, that although (2.10a,b) are in a form similar to the so-called harmonic expansion, the wave fields at second and third order in such a form are introduced for convenience. They are due to the separation of wave harmonics arising mainly from the forcing equations (e.g. (3.1b) and (4.1)) at a still water surface, presented in § 3 at second order and in § 4 at third order.

Linear wave fields
Solutions to a linearized boundary value problem from (2.1)-(2.3) are obtained for linear wave fields. They are given here without detailed derivations (see § 13 in Mei, Stiassnie & Yue 2005), expressed in terms ofÂ(k), as follows: (2.12c) in which V denotes linear velocity vector andV its magnitude in the envelope-type form; u (1) (x, z, t) and w (1) (x, z, t) denote linear velocity vector in the horizontal plane and vertical component, respectively;ū andw are their corresponding magnitude in the envelope-type form; and c p (k) = √ g/|k| denotes the phase velocity of wave k. It is worth noting that (2.12c) is fundamental to the theory presented in this paper due to the following two aspects: it allows for a nonlinear evolution equation for A being expressed in terms of V ; and it also facilitates simple numerical implementations of the evolution equation for A presented in § 4.
Using the definition of linear velocity potential Φ (1) and following similar procedures that lead to the first-order evolution equation (2.8a), it is understood that the following identities hold (to the first order in wave steepness):

Quadratic property of linear waves: velocity head
Before presenting the wave solutions at second and third order in , we introduce a leading-order approximation to the wave-velocity head, H v (x, z, t), defined as which is correct to O( 2 ). The wave-velocity head plays an essential role in the forcing of bound waves at second order and, thus, the nonlinear evolution equation for A presented in § 4. Inserting the envelope-type expression for Φ (1) into (2.14) and using (2.12b), we obtain in which the asterisk '*' denotes the complex conjugates, and we note thatV ·V = u 2 +v 2 +w 2 withū andv denoting the velocity component ofū (ū = [ū,v]) in the x and y direction, respectively. It is noticeable that H where U SD and c g0 denote Stokes drift and the group velocity, respectively. We understand from (2.16) that the mean velocity head H (20) v maintains the energy for the propagation of the Stokes drift at group velocity c g0 for a quasi-monochromatic wave group. We show in § 3.3.2 that it is also responsible for the generation of the Eulerian return flow at second order.
In contrast, the rapidly varying velocity head, associated with the factor of second harmonic exp(2ik 0 · x − 2iω 0 t), has been rarely investigated as it is zero for unidirectional waves. In particular, we note 1 4V ·V e 2k 0 z e 2i(k 0 ·x−ω 0 t) + c.c. = 0, (2.17) which holds for unidirectional waves of a broad bandwidth, as can also be inferred from earlier works, e.g. equations (23) and (25) in Dalzell (1999). Due to this, we examine the effects of the rapidly varying velocity head in § § 3.3.1 and 6.3.1 for multidirectional waves.

Second-order solutions O( 2 )
In this section, the solutions for the second-order superharmonic and subharmonic waves forced by linear surface waves are presented. Compared with Dalzell (1999), the solutions in this section are as exact but expressed in an envelope-type form. In contrast to Chu & Mei (1971), Davey & Stewartson (1974) and Dysthe (1979), the derivations do not rely on a narrow-banded assumption except for § 3.3 where leading-order approximations are presented.

Governing equation and boundary conditions
The velocity potential at second order in wave steepness, Φ (2) = Φ (22) + c.c. + Φ (20) , are described by the following boundary value problem (Dalzell 1999;Li et al. 2021c): (3.1c) in which Γ z is defined as Γ z = ∂ z Γ . As highlighted in (3.1b), the first term on the right-hand side does not contribute to the forcing on a still water surface as it vanishes and, thus, only the second term needs to be considered. Using the second term and (2.15b) we define a forcing term S for z = 0 as follows: and The forcing term S (or the velocity head on a still water surface) is responsible for the forcing of the second-order bound waves. The negative sign on the right-hand side of (3.3a,b) implies a flow generated below a still water surface, propagating in the opposite direction to the change of the velocity head in time at z = 0 -a so-called return flow which is known for second-order subharmonic waves. In order to obtain the solution to boundary value problem (3.1a)-(3.1c) for Φ (2) , two different and novel approaches are proposed in this paper and presented in § § 3.2 and 3.3, respectively. In particular, the first approach is applicable to three-dimensional waves of a broad bandwidth, and the second is an approximate method based on the assumption of a narrow bandwidth with δ x,y 1.
3.2. Approach I: semianalytical approach for Φ (2) Approach I, referred to as a semianalytical approach, proposes solving boundary value problem (3.1a)-(3.1c) in the Fourier k plane by using a pseudospectral and a finite difference method. The boundary value problem for Φ (2j) from (3.1a)-(3.1c) is given by in which j = 0 and j = 2 denote the subharmonic and superharmonic wave, respectively. The solutions of (3.4a,b,c) can be expressed in a form, as follows: for j = 0 and j = 2. Substituting (3.5) into (3.4b) leads to the second-order differential equation forB (2j) , as follows: whereŜ (2j) denotes the Fourier transform for S (2j) for j = 0 and j = 2. It is understood that the second-order waves are bound, and which do not admit a homogeneous solution of (3.6) (Phillips 1960;Hasselmann 1962). As a result, (3.6) forB (2j) can be solved numerically by a semianalytical approach that proposes evaluating S (2j) (x, t) by using a pseudospectral method and obtaining the numerical solution of (3.6) by available time-marching methods. The numerical procedures are explained in § 4.4. For numerical implementations, prescribed conditions at an initial instant are required for Φ (2) (x, 0, t 0 ) and ∂ t Φ (2) (x, 0, t 0 ). In practice, multiple choices are available to this end, depending on the purpose of wave predictions. For instance, a second-order wavemaker theory based on Schäffer (1996), periodic boundary conditions as in Dommermuth & Yue (1987), the framework by Bonnefoy, Le Touzé & Ferrant (2006) for the waves generation in a numerical wave tank, stationary waves based on Dalzell (1999) and approximate initial conditions as presented in § 3.3. If the waves generation is based on linear wave theory, Φ (2) (x, 0, t 0 ) = 0 and ∂ t Φ (2) (x, 0, t 0 ) = 0 are implied, leading to inconsistency and additional generation of spurious waves (Schäffer 1996). The semianalytical approach leads to improved computational efficiency compared with the analytical method by Dalzell (1999) based on Fourier integrals. For N Fourier modes, the semianalytical approach requires computational operations at O(N In(N)), whereas Dalzell (1999) requires O(N 2 ). The validation of the semianalytical approach is presented in § 6.1.

Approach II: an approximate method for Φ (2)
With a narrow-banded assumption, this section presents an approximate method for Φ (22) ( § 3.3.1) and Φ (20) ( § 3.3.2). For identifying the correct order in bandwidth, we assume small and non-dimensional bandwidth scaling parameters δ x and δ y that denote the bandwidth in the x and y direction, respectively. The order of accuracy of the approximate solutions is indicated by the product of 2 and δ j x,y with j non-zero integers. With a narrow-banded assumption, we seek the solutions of (3.1a)-(3.1c) for unknown B (22) and Φ (20) in a form of power series in bandwidth, as follows: in which the subscripts '≈, j' denote an approximation at O( 2 δ j x,y ), and we consider leading-order approximations up to j = 2.

Forcing of second-order superharmonic waves by multidirectional waves
It is understood that the forcing term on the right-hand side of (3.1b) (see also (2.17)) for Φ (22) for unidirectional waves is zero, which suggests Φ (22) (x, z, t) = 0. Hence, the derivations for Φ (22) are only needed for multidirectional waves. Based on the boundary value problem described by (3.1a)-(3.1c), the forcing equation of the superharmonic bound waves for B (22) on a still water surface reads in which the first term on the right-hand side is of O( 2 δ y ) and the second term of O( 2 δ y δ x ). Inserting the expanded form (3.7a) for B (22) into (3.8) and the Laplace equation for Φ (22) and following the conventional procedures for perturbed solutions, we obtain sq,t denote the Fourier transform forV ·V andV · ∂ tV with respect to x for z = 0, respectively.

Potential for the second-order mean flows in two and three dimensions
Similar to the second-order superharmonic waves, the forcing equation of the subharmonic bound waves at second order for Φ (20) on a still water surface (cf. (3.1b)) reads With a narrow-banded assumption and following the same procedures as for B (22) , we obtain Φ (20) g|k| e ik·x+|k|z dk 2 and Φ (20) (20) and ∂ t S (20) with respect to x, respectively.

Wave elevation ζ (2) and velocity V (2) at second order
With second-order potential Φ (2) given by the semianalytical approach presented in § 3.2 or the approximate method presented in § 3.3, the surface elevation at second order is obtained from where Φ (2) = Φ (20) for unidirectional waves and Φ (2) = Φ (20) + Φ (22) + c.c. for multidirectional waves. Inserting the envelope-type expression for Φ (1) and Φ (2) , we arrive and With Φ (22) = 0 in mind for unidirectional waves, (3.13b) gives an exact expression for the elevation for second-order superharmonic waves that depend only on linear modulational parameters, e.g. A andV . Similarly, the velocity at second order, V (2) (x, z, t), is given by In this section, a new NLSE for linear wave envelope A is presented in § 4.1 for multidirectional waves of a broad bandwidth. Based on a narrow-banded assumption, leading-order approximations to the new NLSE are derived in § 4.2, and how this new evolution equation recovers existing NLSEs is explained in § 4.3. The numerical implementations of the new NLSE are presented in § 4.4.

An evolution equation correct to O( 3 ) for waves of a broad bandwidth
Collecting the terms at third order in wave steepness in (2.2), we obtain where the forcing term Q (3) is given by which are not necessarily so for a finite depth. The terms in the square bracket on the right-hand side of (4.2) can be simplified due to (3.1b). In particular, it is noticed that the following identity holds: Inserting the harmonic expansion for Φ (1) and Φ (2) into (4.2) leads to Q (3) being expressed in a form as follows: where Q (31) denotes the forcing term with a factor exp (ik 0 · x − iω 0 t) and Q (33) the term with a factor exp (3ik 0 · x − 3iω 0 t). Leaving the expression for Q (33) in § 5 and focusing on Q (31) in this section, we obtain (see details in Appendix A) in which the subscripts 'all' and 'dir' are used to distinguish two different cases; the former means the parameter expressed in terms that are non-constant for both unidirectional and multi-unidirectional waves, and 'dir' means the parameter composed of terms that are non-zero only for directional waves. Hence,Q (31) admits a much simpler expression for waves considered in two dimensions or for long-crested waves due toQ Q (31) would lead to secular solutions that should be removed (see the discussion on page 376 in Madsen & Fuhrman (2006)). In order to achieve this, a conventional approach is to introduce a nonlinear amplitude-dependent frequency ω 2 (∼O( 2 )) for the first-order wave fields through replacing factor exp Thereby, it leads to an additional contribution of the updated first-order potential Φ (1) to (4.1) at O( 3 ). Especially arising from ∂ tt Φ (1) , we obtain an additional term associated with first harmonic at Mathematically, the secular solutions can be removed if the sum of this additional term andQ (31) equal zero. Hence, we arrive at Noticing that ζ (1) = −∂ t Φ (1) /g for z = 0 which gives (∂ t − iω 0 )B = −gA, substituting this expression for B into (4.6) leads to the equation for A, as follows: The evolution equation for A with A = A exp(−iω 2 t), correct to third order in wave steepness, can now be expressed as Inserting (4.7) into (4.8) and omitting the prime leads to (4.9) where N = −Q (31) /g denotes the nonlinear term, given by Equation (4.9) denotes the new NLSE that describes the evolution of linear wave envelope A in time and space. In addition to the linear term associated with the time derivative ∂ t (. . .), (4.9) is composed of a term that denotes linear dispersion relation due to iω 0 L and a nonlinear term denoted by N that is correct to O( 3 ). For (4.9), no assumptions associated with the bandwidth have been introduced. Thereby, it has relaxed the limitation arising from a narrow-banded assumption on which most NLSEs are based, such as the NLSEs by Dysthe (1979), Trulsen & Dysthe (1996) and Trulsen et al. (2000); i.e. the new NLSE should permit the prediction of waves with bandwidth δ x,y > 1. The implementation of the new NLSE for the wave envelope, A, of the first-harmonic elevation requires the evaluation of linear term ∂ t A in the Fourier plane based on (2.8c) and the evaluation of the linear envelope (vector) of both the linear and second-order velocity for N , which will be explained in detail in § 4.4. Extensions of the new NLSE (i.e. (4.9)) to more general cases are simple and straightforward. For example, the new NLSE can be extended to consider a finite water depth following the derivations in this paper, except that one would expect that the nonlinear term N may contain a few more terms that introduce a negligible additional computational cost. Following Dias et al. (2008) and Kharif et al. (2010), the new NLSE would allow for weakly damped/forced free surface flows through adding the viscous effects in the dynamic and kinematic boundary condition (cf. equations (35, 37) by Dias et al. (2008)) at the water surface. Moreover, despite that the new NLSE is formulated in the framework of potential flow theory for simplicity, it is essentially an equation associated with linear envelope A, linear velocityV and second-order velocityV (2) . This suggests the weak effects due to viscosity and vorticity in a subsurface layer can be easily incorporated, allowing for waves interacting with ambient environments, e.g. turbulence and background rotational flows (cf. McWilliams, Restrepo & Lane 2004). These aspects will be considered in future work.

Leading-order approximations to the new NLSE
Approximations to (4.9) correct to O( 3 δ j x,y ), with j = 0, 1, 2, 3, are presented in this section with a narrow-banded assumption, i.e. O(δ) 1. In particular, we assume that the bandwidths both in the longitudinal (x) direction and the transverse direction to the propagation of the carrier wave are small and that O(δ y ) ∼ O(δ 2 x ). Small directional spread is, hence, implied. Based on previous papers (e.g. Chu & Mei 1971;Li et al. 2021c), the order of magnitude associated with the relevant terms in (4.5) is assumed, as follows: and Moreover, it is understood that the derivative operators (∂ t , ∇ 3 ) on the terms in (4.11) would lead to an order of magnitude higher in bandwidth δ x,y . With the right orders of magnitude in mind, collecting the terms up to a particular order in and δ based on (4.5) and (4.9b), and, thus, leading-order approximations to (4.9) are obtained. In particular, they are expressed in a form, as follows: x,y ), with j = 0, 1, 2 and 3, given by where (4.13e) denotes that N (3,4) has an identical expression to N . As noted earlier, the subscript 'dir' means the parameter being composed of nonlinear terms that are non-zero only for multidirectional waves, given by Equation (4.13b) denotes that a leading-order approximation to N , correct to O( 3 δ), does not require an evaluation of the second-order mean potential. Equation (4.13c) suggests that the approximations to Φ (20) and B (22) that are presented in § 3.3 are sufficient to give a leading-order approximation to (4.9), correct up to O( 3 δ 3 ). Any approximations correct to an order higher than O( 3 δ 2 ) would require the same computational cost as N because the approximations Φ ≈,2 and B ≈,2 do not lead to enhanced computational efficiency, compared with solving for Φ (2) by the semianalytical method.

Comparisons with two reduced-form equations for A
Two limiting cases are considered in this section in order to demonstrate that the new NLSE can recover two reduced-form equations for A, including the NLSE for the evolution of a train of Stokes wave ( § 4.3.1) and the NLSEs by Dysthe (1979), Trulsen & Dysthe (1996) and Trulsen et al. (2000) with a narrow-banded assumption ( § 4.3.2).

Uniform Stokes waves
For uniform Stokes waves on an infinite depth, that denote A andV having no slow modulation in x and z, the new NLSE leads to N = N (3,0) and hence, N = 1 2 ik 2 0 ω 0 |A| 2 A. (4.14) Introducing = k 2 0 |A| 2 and α 1 = 1 2 for Stokes waves, we obtain the evolution equation for A as follows: which agrees with the classic NLSE for a train of Stokes waves as presented in many papers, e.g. Benjamin & Feir (1967) and Zakharov (1968). Equation (4.15) will be used for the analysis of sideband instability of Stokes waves presented in § 6.2.

Waves of a narrow bandwidth (O( 3 δ x,y ))
In order to recover the previous versions of NLSEs by Dysthe (1979), Trulsen & Dysthe (1996) and Trulsen et al. (2000) from the NLSE (4.9), two aspects are addressed. First, the term associated with linear dispersion relation iω 0 L is identical to equations (12) and (13) in Trulsen et al. (2000) using the dimensional versions. It has been demonstrated in Trulsen et al. (2000) how this term can recover the terms associated with linear dispersion relation derived by Trulsen & Dysthe (1996) and Dysthe (1979). Second, the nonlinear term in the NLSEs by Trulsen et al. (2000) and Trulsen & Dysthe (1996) are directly obtained from the Dysthe equation by Dysthe (1979). Hence, we only have to demonstrate that N can recover the nonlinear terms in the Dysthe equation. To this end, we start with (4.9b) for unidirectional waves as the Dysthe equation has neglected the contribution from Φ (22) (which is non-zero only for multidirectional waves). Moreover, a narrow-banded assumption as in Dysthe (1979) With a narrow-banded assumption, the velocity (magnitude)V correct to O( δ x,y ) reads (cf. equations (13.2.21) and (13.2.30) by Mei et al. (2005) for deep-water waves) in which e x = [1, 0] denotes a unit vector along the x direction in the horizontal plane. Inserting (4.16) into (4.9b) leads to a leading-order approximation to gN , as follows (see Appendix A.2 for details): where C is defined as C ≡ B/2 = −iω 0 A/(2k 0 ). The subscript 'Dysthe' means that the parameter corresponds to the 'Γ ' in Dysthe (1979); i.e. (4.17) recovers 'Γ ' on the right-hand side of equation (2.17) of Dysthe (1979), while observing the following two aspects. Firstly, the notation C defined here is denoted by 'A' in Dysthe (1979 Figure 1 shows the detailed procedures of the numerical implementation of the new NLSE (4.9) for A, including the solution of the second-order equation, (3.6), for B (2j) and the evaluation of N by a pseudospectral method and a finite-difference method. We remark here that it is conventional to solve an NLSE numerically for A in a moving coordinate system through the transforms ξ = x − c g0 t and η = t, which would also follow similar numerical procedures as follows, but for A (ξ, y, η) with A (ξ, y, η) = A(x, t) instead.

Numerical implementation
At any temporal step t n = t 0 + n t with t 0 denoting the initial time, it is understood that the value of A(x, t n ) for all x s < x < x e , with x s and x e denoting the start and end position vector of a chosen spatial domain, respectively, is given from earlier computations. The spatial domain is discretized by 2N x × 2N y points at which A(x, t n ) is known, with 2N x and 2N y denoting the total number of points in x and y direction, respectively. Transforming A to the discretized Fourier space using a fast Fourier transform (FFT) (Frigo & Johnson 1998), we obtainÂ where dk x and dk y denote the interval between two adjacent points in the Fourier k x and k y direction, respectively and F denotes a FFT with respect to x. As highlighted in grey in figure 1 with the known value ofÂ(k, t n ), the new NLSE is solved numerically by using a split-step Fourier method (Tappert 1974). The main difference between the numerical procedures presented here and these presented in Lo & Mei (1985), lie in the fact that the evaluation of N per time step is obtained differently as N has different expressions. The numerical procedures for N at t = t n are decomposed into three consecutive steps, as explained in the following. ( Step I) The first step is to use the known value ofÂ(k, t n ) to evaluateV (x, 0, t n ) and its derivative with respect to t, x, y and z using an inverse FFT, as highlighted in figure 1 in blue. Let V be the transformed linear velocity,V , in the Fourier k space. We understand from (2.12) that (4.19a) and thus Hence,V and also its derivatives [∂ t , ∂ x , ∂ y , ∂ z ]V are evaluated for all spatial points through an inverse FFT using (4.19) for z = 0.
(Step II) The second step (highlighted in green in figure 1) is to solve the second-order equations (3.6) for B (2j) with j = 0 and j = 2 by the semianalytical approach presented in §3.2. The forcing terms S (2j) for j = 0 and j = 2 are computed based on (3.3) in the physical plane using ∂ tV ,V andV * evaluated in step I. The second-order potentialsB (2j) and their derivatives with respect to time, ∂ tB (2j) for j = 0 and j = 2, with prescribed initial conditions, can be evaluated based on (3.6) by using a time-marching method, e.g. the midpoint method as indicated in figure 1. With the potentials computed per time step, second-order velocity at different harmonics,V (2j) for j = 0 and j = 2, and their derivatives with respect to t, x, y and z can be readily obtained by an inverse FFT. (Step III) The third step, as highlighted in purple in figure 1, is to evaluate N in the physical space based on (4.9b) using the terms evaluated in step I and step II.
Some remarks on the numerical efficiency of the new NLSE are explained here for two-dimensional waves, which allow for a fair comparison with the efficiency of existing NLSEs, e.g. the Dysthe equation as explained in Lo & Mei (1985). In the pseudospectral method, the nonlinear terms are evaluated through FFT at each x and t n before the integration in time is performed. The evaluation requires 20 FFTs (cf. the areas highlighted in blue and green in figure 1). The third-order nonlinear terms are evaluated in the physical domain, and an additional eight FFT computations arise from the terms associated with ∂ t (∂ x , ∂ z )B (2j) and (∂ x , ∂ z )B (2j) due to an inverse FFT (cf. the areas highlighted in green and purple in figure 1). Two FFTs are needed in the linear equation (i.e. the equation in the last line highlighted in grey in figure 1). If N Fourier modes are used, a total of 30N In(N) operations are needed for advancing the solution by one step in time whereas 10N In(N) operations are required with NLSEs, e.g. the Dysthe equation and modified NLSE (Dysthe 1979;Lo & Mei 1985;Trulsen et al. 2000). Although a slightly increased computational cost relative to the previous versions is inevitable, as expected, the new NLSE offers much better computational efficiency compared with other methods, e.g. solving the Zakharov integral equation that requires roughly O(N 3 ) operations when an explicit time-differencing scheme is used. Moreover, it is worth noting that the kinematics of the system, i.e. the elevation and velocities at the first and second orders, are obtained at almost no additional cost throughout the numerical implementations of the NLSE.
In contract to the HOS method which is also based on efficient FFTs, a pseudospectral method and numerical methods for time marching (cf. Dommermuth & Yue 1987;West et al. 1987), the new NLSE has the following main features. First, as aforementioned, the new NLSE allows for a computational domain chosen based on the scaling of the wave envelope (or the maximum of the side bandwidth of a spectrum), in a manner similar to other NLSE-based models. This suggests a smaller number of discrete points than the HOS method can be chosen in a computational domain. It is mainly achieved by the semianalytical approach for second-order superharmonic bound waves, which are expressed in an envelope-type form. It is equivalent to shifting the spectrum of the second-order superharmonic bound waves by 2k 0 towards the origin of the Fourier plane. In practice, the choice of the carrier wave (i.e. k 0 ) allows for flexibility and, therefore, can contribute to reducing the computational cost, in particular for the ocean spectra of a long tail. With a reasonable selection of k 0 , the side bandwidth of an asymmetrical spectrum can be much reduced and, therefore, a less fine mesh would be allowed for the implementation of the new NLSE. For example, a numerical test for a JONSWAP spectrum (see Appendix C) was carried out with k 0 = 4k p chosen for computations, where k p denotes the spectrum peak. Secondly, the new NLSE separately obtains the wave fields due to free waves and second-order bound waves at a still water surface. Last, with proper separation of the perturbed solutions and multiple scales in time (in terms of wave steepness), the author conjectures that one would show that the solution of the third-order HOS method following Onorato, Osborne & Serio (2007) (i.e. equations (13) and (14) in the paper) can recover the new NLSE for the envelope of the first-harmonic wave elevation. This will be addressed in future work.

Solutions at third order O( 3 )
For completeness, this section presents the solutions at third order for both potential Φ (33) given in §5.1 and elevation ζ (3) in § 5.2.

Potential Φ (33) at third harmonic
Collecting the terms at third order in wave steepness in (2.1)-(2.3) and keeping the terms that have a factor exp(3ik 0 · x − 3iω 0 t), we obtain the following equations for Φ (33) : where the third-order term Q (33) (x, t) denotes the forcing of bound waves at third harmonic on a still water surface, as noted in § 4, expressed as and For unidirectional waves, it is understood that both B (22) and (V ·V ) have no contribution to the third-order nonlinear terms as they vanish, meaningQ (33) = 0 and thus, Φ (33) = 0. The solution of equations (5.1a,b,c) for Φ (33) can be obtained in a way similar to the semianalytical approach for Φ (2) , i.e. the solution can be obtained by a finite-difference method in the Fourier k space based on the forcing equation on a still water surface, 926 A34-19 Y. Li as follows: where, as defined in previous sections, '(. . .)' denotes a Fourier transform to the k plane with respect to x. Thereby, Φ (33) is expressed as Based on (5.5), it is understood ζ (3) can be expressed in a form, as follows: where amplitudes A (31) (x, t) and A (33) (x, t) are in a form given by, respectively, in which the subscript 'uni' refers to the cases for unidirectional waves. With a narrow-banded assumption O(δ) 1, we readily obtain from (5.9) (through omitting terms associated with the derivatives with respect to t and z and the second-order mean fields and notingw = −iω 0 A + O( δ)), as follows: in which the subscript '≈' denotes an approximate expression due to a narrow-banded assumption. Inserting (5.10a,b) into (5.6) leads to a leading-order approximation for ζ (3) that agrees with (2.12) in Lo & Mei (1985).

Results
This section focuses on three aspects in order to explore and validate the theory presented, including the second-order solutions presented in § 3, sideband instability of Stokes waves and the roles of wave directionality. They are examined in § § 6.1, 6.2 and 6.3, respectively. For numerical implementations, a Gaussian amplitude spectrum (Ξ(k)) in wavenumber and a (unnormalized) directional distribution (D(θ )) were chosen to generate a short (focused) wave group, given by, respectively, where k p denotes the peak wavenumber and θ p denotes the propagation direction of the peak wave at linear focus, k w and θ w denote the standard deviation of the spectrum in wavenumber and direction, respectively. An asymmetrical Gaussian spectrum is obtained from (6.1a) through using two different values of k w for the upper sideband k > k p and lower sideband k ≤ k p , respectively. For the relevant cases computed in this section, we used θ p = 0 and k p = 0.02769 m −1 which is typical of the North Sea (Barratt et al. 2021). Using the spectra described by (6.1a,b), linear elevation ζ (1) at t = t 0 is given as input by for unidirectional waves, where 0 denotes the dimensionless steepness of a wave group at linear focus; k n and N k denote the evenly spaced discrete wavenumbers and the total number chosen for numerical implementations, respectively; the subscript 'f ' denotes the position or time for the wave group at linear focus; ψ f denotes the wave phase at linear focus; for multidirectional waves, where θ m and M θ denote the evenly spaced discrete wave directions in the range of (−π, π) and the total number chosen for numerical implementations, respectively. Despite that the directional distribution described by (6.1b) is in an unnormalized form, (6.3a,b) suggest that the integral unity of the directions of the linear multidirectional waves generated has implicitly been satisfied due to the additional scaling factor, Ξ k,θ , defined in (6.3b).

Second-order solutions
We validate the second-order solutions presented in § 3 in this section through comparisons with the analytical method by Dalzell (1999) and with existing approximate methods. We focus on the temporal-spatial evolution of a wave group. For simplicity, the effects at third order are neglected in this section as they do not affect, qualitatively, the discussion presented here. Let a wave group start to propagate at t = t 0 at which the linear elevation ζ (1) (x, t 0 ) of a group is given in space by (6.2) for a unidirectional group and by (6.3a) for a directionally spread focused wave group. Inserting ζ (1) (x, t 0 ) into (2.6) leads to an expression for A(x, t). For a fair comparison, the input at t = t 0 for second-order wave fields is based on the results by Dalzell (1999). Following the numerical procedures explained in § 4.4, ζ (22) and ζ (20) are readily obtained from (3.13b) and (3.13c), respectively. Figure 2 shows a comparison of the second-order superharmonic wave elevation for a unidirectional wave group at two different times (noting that Φ (22) (x, t) = 0), between the results predicted by (3.13b), Dalzell (1999) and a leading-order approximation based on a second-order Stokes theory given by

Unidirectional and multidirectional wave group
Due to an asymmetrical Gaussian amplitude spectrum chosen in figure 2, a good estimate of the side bandwidth can be given by δ = 3k w /k p for k > k p , which measures from the wavenumber (i.e. k = k p + 3k w ) dropping by 99 % in magnitude relative to the spectrum's peak. This gives δ = 4.86 for the case examined in figures 2 and 3. The good agreement between (3.13b) and Dalzell (1999) is evident in figure 2, whereas the less satisfactory agreement between the second-order superharmonic elevation (6.4) and the other two methods is also shown. For completeness, the expressions for both ζ (22) and ζ (20) based on The present theory Exact (Dalzell 1999 -12 -10 -8 6 -9 -9.5 -10 4 2 -6 -4 -2  Dalzell (1999) are presented in Appendix B. It is understood Dalzell (1999) is capable of providing exact predictions of second-order wave fields based on Fourier integrals. The agreement with Dalzell (1999) in figure 2 clearly demonstrates that the semianalytical approach presented in § 3.2 can also work for this purpose. Similar observations are also shown in figure 4 for ζ (22) predicted for the evolution of a directionally spread focused wave group at different times. Figure 3 shows a comparison of the elevation for second-order subharmonic waves forced by a unidirectional linear wave group on a still water surface between the results predicted by (3.13c), Dalzell (1999) (cf. Appendix B) and Trulsen et al. (2000) (cf. Appendix D). It is clearly seen in figure 3 that (3.13c) can predict ζ (20) as well as Dalzell (1999), whereas Trulsen et al. (2000) leads to less satisfactory agreement with Dalzell (1999), owing to a narrow banded assumption.

Sideband instability of Stokes waves
In this section, we examine the capabilities of the new NLSE by investigating the sideband instability of Stokes waves through comparisons with earlier works, e.g. The present theory Exact (Dalzell 1999) Trulsen et al. (2000 -25 -20 -15 -10 -15 -14 -13  Dysthe (1979), Crawford et al. (1981), McLean (1982) and Trulsen et al. (2000). The approximations to N presented in § 4.2 are used in this section as they are sufficient to this end. The contribution due to N ≈,dir is neglected for simplicity. Based on (4.15), the stability of Stokes waves to sideband perturbations can be investigated by assuming small disturbances in amplitude and phase for A in a form given by (cf. Trulsen et al. 2000) in which a 0 denotes the amplitude of a train of Stokes wave, 0 = k 0 a 0 denotes the dimensionless steepness of the Stokes wave,â andθ are infinitesimal real parameters (i.e. O(1)) that denote small perturbations in amplitude and phase, respectively, ψ(x, t) = k 0 [(1 + δ x )x + δ y y] − ω 0 δ Ω t denotes the modulated phase, with k 0 and ω 0 the wavenumber and angular frequency of the Stokes wave, respectively, δ x and δ y the dimensionless sideband in wavenumber in the longitudinal (k x ) and transverse/spanwise (k y ) direction, respectively, and δ Ω the dimensionless sideband in frequency; as defined in § 4.3.1, α 1 is defined as α 1 = 1/2. Treating the expression (6.5) for A as the modulated amplitude of a plane wave k 0 and inserting (6.5) into (2.12) leads toV (x, z, t) given by The full solution forV in a form as (6.6) allows for min(|δ x |, |δ y |) > 1. It should be noted that, time dependant disturbances (small) forV except for these due to linear dispersion relation are neglected for simplicity. Otherwise, the derivations following Benjamin (1967) (e.g. starting from equations (11) and (21) in the paper) are needed that would lead to numerical analysis for the sideband instability. The insertion of (6.6) into (4.12) for j = 3 leads to an eigenvalues equation for δ Ω , with the terms at O( 2 0 δ 2 y ) and higher neglected. In particular, the eigenvalues equation for δ Ω reads denote the contribution of the linear dispersion, and coefficients α ij are given by in which the z dependent variables associated with V i,± (z) are evaluated for z = 0, δ 20 = δ 2 x + δ 2 y denotes the dimensionless magnitude of the sidebands away from the wavenumber of the Stokes wave, and the approximation L 2 ≈ 1 2 (L I + L r ) (cf. equations (11) and (15) by Crawford et al. (1981) for similar discussions) was used due to the terms associated with ∂ t (V ·V * ) and, therefore, V (20) , arising from the symmetry of the disturbances to the Stokes wave.
It is understood that the instability occurs if (6.8), a quadratic equation in δ Ω , has a positive imaginary root for δ Ω . Thereby, the instability occurs if the following inequality holds: (6.10) in which R Ω is defined for later reference. Figure 5 shows a comparison of the instability region predicted by (6.10) and its lower-order expression, equations (18-22) by Trulsen et al. (2000), and the exact computations by McLean (1982) for three different values of wave steepness. The equations of Trulsen et al. (2000) are presented in Appendix D for reference. In the cases chosen by McLean (1982) the value of a different steepness, defined as half the crest-to-trough height of Stokes waves and referred to as M here, is given. For a direct comparison, the identity derived by Trulsen & Dysthe (1996) (i.e. equation (30) therein) was used to evaluate 0 , as follows: (6.11) It is seen from figure 5(a-c) that, compared with McLean (1982), both approximations behave similarly; good agreement is shown for small δ x and small steepness 0 , whereas the agreement becomes less satisfactory for larger 0 . As approximations of higher-order accuracy are used, better agreement with McLean (1982) is clearly seen in figure 5(d-f ) compared with figure 5(a-c). In particular, as seen in figure 5(d-f ), equation (6.10) agrees with McLean (1982) in a nearly perfect manner for 0 ≈ 0.1 and 0 ≈ 0.2 in the range δ x 1, while the predictions by (6.10) are also satisfactory for 0 ≈ 0.33. Comparing the results by Trulsen et al. (2000) and (6.10) with McLean (1982) shown in figure 5(d-f ), it is clear that (6.10) provides a better prediction for all three cases. In particular, (6.10) is capable of proving more accurate predictions for the instability of Stokes waves subject to sidebands in the transverse direction, as clearly shown in figure 5(e, f ). Focusing on the instability due to transverse sidebands (δ y ), figure 5 shows that equations (19-22) in Trulsen et al. (2000) fail to improve the accuracy of the lowest-order approximation (i.e. equation (18) by Trulsen et al. (2000)) for small δ x for all three cases. The reason to this is explained as follows. Equations (19-22) of Trulsen et al. (2000), with the notation in this paper, lead to the instability region described by the inequality, as follows: where R T is introduced for later reference. Equation (6.12) gives the instability region based on equation (18) by Trulsen et al. (2000) by setting δ x to zero. For small 0 ( 0 0.5), the last term in R T is negligible, which gives two conditions for instability to occur, Noting that δ 2 x /δ 20 < 1 for δ x 1, the first condition in (6.13) dominates for small δ x . Using a leading-order approximation to L r , i.e. L r ≈ (2δ 2 y − δ 2 x )/8, the first condition is 926 A34-27 approximately 2δ 2 y − δ 2 x < 0, (6.14) which corresponds to a (neutral) instability slope near the origin defined by δ y /δ x = ±1/ √ 2 (corresponding to ∼ ±35.26 • to the carrier wave k 0 ), as pointed out by Longuet-Higgins (1976) for the analysis of resonant energy transfer in a narrow spectrum and clearly shown by the predictions based on Trulsen et al. (2000) in figure 5 (see also § 6.3.2).
In contrast, (6.10) is capable of providing an accurate prediction of the instability region of Stokes waves subject to spanwise (y direction) sideband perturbations for moderate values of wave steepness, 0 0.2, as can be clearly seen in figure 5(d,e). This is due to that (6.6) has taken into account the transverse sidebands δ y in δ z,± which denotes the effects on wave kinematics (i.e.V ) in addition to the wave envelope A. For 0 ≈ 0.33, the difference between (6.10) and McLean (1982) is likely due to that the terms described by N ≈,dir presented in § 4.2 at O( 2 0 δ 2 y ) were neglected in the computations. 926 A34-28 Equations (19-22) by Trulsen et al. (2000) Implementation of N ª (33) Crawford et al. (1981) Lake et al. (1977 for δ x = 0.2 Benjamin (1967) for δ x = 0.2 Lake et al. (1977) for δ x = 0.4 Equations (19-22) by Trulsen et al. (2000) Implementation of N ª Figure 6. The growth rate Im(δ Ω ) varying with the longitudinal bandwidth δ x (a) and 0 = k 0 a 0 (b) when the instability of Stokes waves occurs. In the figure, the predictions are computed based on Crawford et al. (1981), (6.15), and equations (19-22) by Trulsen et al. (2000) for three different values of steepness (a) and two different values of δ x (b). Panel (b) shows, in addition, experimental measurements by Benjamin & Feir (1967) and Lake et al. (1977).
We now proceed to examining the energy growth rate of Stokes waves, which can be defined as the positive imaginary component of the roots of (6.8), as follows: where Im denotes the imaginary component. When instability occurs, the rate of amplitude growth/decline behaves like exp[ 1 2 Im(δ Ω )ω 0 t] based on (6.5). This will be further examined in § 6.3.2. Figure 6 shows the energy growth rate, Im(δ Ω ), varying with the sideband δ x for five values of 0 and varying with 0 for two values of δ x . In addition to the experimental measurements and the predictions of (6.15), figure 6 also presents the results due to Crawford et al. (1981), for which small corrections following Krasitskii (1994) and Janssen (2004) (cf. § 4.11) were used. Crawford et al. (1981) is based on the Zakharov integral equation for the study of a uniform wave train and the results are correct to third order in wave steepness with all higher-order dispersion effects included. It is seen in figure 6 that (6.15) shows better agreement with Crawford et al. (1981) than Trulsen et al. (2000) (i.e. predicted by Im √ R T ), which is more so for 0.15 0 and 0.4 δ x . Figure 7 shows the instability boundary for growth of unstable perturbations for two-dimensional wave trains, predicted due to Benjamin & Feir (1967), Longuet-Higgins (1978), Crawford et al. (1981) and by this paper based on (6.15). Compared with the exact numerical results due to Longuet-Higgins (1976), (6.15) is in qualitative good agreement whereas the quantitative discrepancy is ∼23 % for the restabilization of steeper waves perturbed by long waves with δ x → 0. Although the results based on (6.15) in figure 7 show slightly worse performance than the Zakharov equation, it can obviously yield restabilization for larger wave steepness, being an improvement over earlier versions of NLSEs.
For large values of wave steepness, i.e. 0.2 0 , as shown in figures 6(a,b) and 7, the disagreement between (6.15) and Crawford et al. (1981) can be observed, which likely arises from that the additional time-dependant (small) disturbances are neglected in (6.15) 0.5 Benjamin & Feir (1967) Longuet-Higgins ( Figure 7. Stability boundary for growth of unstable perturbations for two-dimensional wave trains predicted based on (6.15), compared with the results of Benjamin & Feir (1967), Longuet-Higgins (1978) and the Zakharov equation due to Crawford et al. (1981).
to the linear velocity envelope (vector),V . The analysis in this section shows the results due to (6.15) are satisfactory. Figures 5-7 suggest that the new NLSE described by (4.9) should work for 0 0.25 and δ x,y > 1 for the evolution of surface waves with directional spread.

Effects of wave directionality
The presented theory permits us to investigate the effects of wave directionality based on an order in wave steepness. Thereby, we focus the wave-directionality effects on waves at second order ( § 6.3.1) and on nonlinear energy transfer in a narrow wave spectrum by implementing the new NLSE numerically in § 6.3.2.
6.3.1. Waves at second order As explained in § 2.4, one of the major effects of wave directionality arises from the fact that the superharmonic velocity head H (22) v throughout the water columns becomes non-zero for directionally spread waves. This leads to both non-zero Φ (22) , additional contributions to elevation ζ (22) due to both Φ (22) and H (22) v . These second-order effects on the spatial-temporal evolution of potential or elevation envelope of the first harmonic are considered negligible in most models based on an earlier version of NLSE, e.g. the NLSE by Dysthe (1979), Stiassnie (1984) and Trulsen et al. (2000). We focus on examining the velocity head H (22) v and H (20) v by studying the evolution of a directionally spread focused wave group at different times. For simplicity, third-order effects were neglected in the computations as they do not affect the discussions here in a qualitative manner. Figures 8 and 9 show the dimensionless superharmonic and subharmonic velocity head defined in § 2.4, respectively. It is seen from figure 8 that, except for the group at linear focus (figure 8b,e,h), larger values of directional spread lead to smaller magnitudes in superharmonic velocity head, owing to the fact that, as the total wave energy is kept the same, the wave energy distributes in larger areas for larger values of directional spread. The magnitude of the superharmonic velocity head is shown to be 10 ∼ 30 % × 0 in figure 8. Similar observations can be found in figure 9 for the subharmonic velocity head which is of a similar order of magnitude to the superharmonic velocity head shown in figure 8. This indicates both components can be equally important for waves with directional spread.

Nonlinear energy transfer in a narrow wave spectrum
In order to examine the nonlinear energy transfer in a narrow wave spectrum, the temporal evolution of a 'steep' wave group is investigated numerically with the implementation of the new NLSE and equations (19-22) by Trulsen et al. (2000). Figure 10 shows the evolution of the discrete amplitude spectra,Â(k, t), in the first quadrant of the Fourier space for a wave group at different times with 0 = 0.2. It is seen in figure 10 that the results based on (4.9) and Trulsen et al. (2000) do not show differences in a qualitative manner and the energy conservation described by  were verified (not shown) to be within 0.01 % for both models. The dimensionless sidebands in the longitudinal and transverse direction are defined as, respectively, (6.17a,b) in which k 0 was chosen to be the peak wavenumber of the initial wave spectrum, δ x > 0 and δ x < 0 denote the region of upper and lower sidebands in the longitudinal direction, respectively. Compared with the initial time shown in figure 10(a,d), figure 10(b,e) shows broadened spectral amplitude in both the lower and higher sidebands in the longitudinal direction but narrowed in the region of larger sidebands in the transverse (k y ) direction, suggesting energy transfers from the latter to the former.
From t = 0 to t = 24 × T 0 with T 0 denoting the wave period of the spectral peak, two main features can be observed in figure 10. Firstly, oblique energy transfers are clear and the change of the spectral amplitude in the regions of the lower (k x /k 0 < 1 or δ x < 0) and upper (k x /k 0 > 1 or δ x < 0) sidebands behaves differently. This indicates asymmetrical change rate of the spectral amplitude in sidebands. Particularly, in the region of the lower sidebands described by 0 < δ  is indicated whereas an energy expansion is obvious in the region of upper sidebands, especially for δ y > δ x / √ 2. Secondly, relative to the original spectral peak at k x /k 0 = 1, a small downward shift of the spectral peak can be observed in figure 10(c, f ), as is clearly shown by the nearest contour in the vicinity of (1, 0) that is no longer symmetrical relative to k x /k 0 = 1.
In order to examine the change rate of the spectral amplitude in time, |Â(k, t)| were calculated based on (4.9) at different times, and |Â(k, t)| at 35 different wavenumbers were chosen and shown in figure 11. Figure 11 shows the non-dimensional spectral amplitude as a function of the dimensionless time 1 2 2 0 ω 0 t. As is seen in figure 11, the change rate of |Â(k, t)| behaves differently at different wavenumbers. The most interesting range is obviously for 1 2 2 0 ω 0 (t − t 0 ) in ∼ [1, 3] in which the amplitude changes rapidly. This range is in the vicinity of t = t f at which the wave group focuses linearly. We refer to t t f and t f t in this time range as the prefocus and postfocus phase, respectively. Focusing on the prefocus phase first, figure 11 shows the amplitude grows for the wavenumbers in the range δ y 0.1 and wavenumbers in the upper sidebands with δ x = 0.36 for δ y 0.35. In contrast, the amplitudes for larger sidebands in the transverse direction, i.e. for 0.18 δ y , experience a decline except for δ x = 0.36. In the postfocus phase, the amplitudes that have experienced a growth (decline) in the prefocus phase tend to decrease (grow) exponentially. Especially, for the amplitudes at the wavenumbers in  Figure 11. Amplitude spectra k 0 |Â(k, t)| at 35 wavenumbers in the first quadrant of the Fourier k space as a function of the dimensionless time 1 2 2 0 ω 0 (t − t 0 ). The results are based on the new NLSE and the parameters were chosen the same as figure 10. In the figure, δ x = (k x − k 0 )/k 0 and δ y = k y /k 0 denote the dimensionless sidebands in the longitudinal and transverse/spanwise direction, respectively. Short vertical lines near 1 2 2 0 ω 0 (t − t 0 ) = 2 indicates for t = t f at which the wave group focuses linearly. small transverse sidebands, i.e. δ y 0.09, the rate of decline is larger than that of increase. The asymmetry in the rate of amplitude change in different regions is likely the cause for the downward shift of the spectral peak shown in figure 10(c, f ).
The change rate of In(k 0 |Â|) in the range for 1 figure 11 shows that In(k 0 |Â|) decreases/grows with 1 2 2 0 ω 0 (t − t 0 ) (approximately) linearly with absolute slopes 1. This suggests that |Â(k, t)| changes exponentially with the dimensionless time 1 2 2 0 ω 0 (t − t 0 ), unlikely arising from resonant energy transfers as resonance would lead to the change rate to be linear in time. The rapid change of the spectral amplitudes shown in figure 11 may be due to instability. This paper argues that it is likely due to a mechanism similar to oblique sideband instability of Stokes waves. To illustrate this, the amplitude change rate in the instability region of Stokes waves is shown in figure 12 for five different values of wave steepness. The results in figure 12 were computed based on (6.15) and equations (19)(20)(21)(22) in Trulsen et al. (2000) (i.e. Im √ R T ). As noted in § 6.2, it is understood that, when instability occurs, the amplitude behaves like |A| ∼ a 0 e Im(δ Ω )ω 0 (t−t 0 ) and thus, R c ≡ In(k 0 |A|) where R c denotes the change rate of amplitude in the logarithmic scale, introduced for convenience. It also corresponds to the gradient of In(k 0 |A|) with respect to 1 2 ω 0 2 0 (t − t 0 ) shown in figure 11.
Based on (6.18b), it is clear that R c depends on both 0 and the magnitude of Im(δ Ω ). Due to 0 < 0.5 in general, which means In( 0 ) < 0, whether or not the amplitude grows in time depends on the sign of Im(δ Ω ); a positive sign leads to a decline whereas a negative sign a growth. In practice 0.1 0 0.3 would give a good estimate of the steepness range, yielding −2.3 In( 0 ) −1.2. Figure 12 Figure 12. Instability growth rate of Stokes waves, defined as Im(δ Ω )/( 2 0 /2), subject to oblique sideband perturbations for five values of wave steepness 0 = k 0 a 0 , obtained using the prediction by (6.8) (a-e) and equations (19-22), ( f -j) by Trulsen et al. (2000) for five different values of steepness. transverse and longitudinal sidebands at five different values of 0 varying from 0.1 to 0.3. The contours of large growth rates shown in figure 12 suggest that oblique energy transfers would be favoured when instability occurs, although the maximum δ Ω always lies on the δ x axis for deep-water Stokes waves. The magnitude of δ Ω /( 1 2 2 0 ) predicted by both the present theory and Trulsen et al. (2000) is 1. Obvious 'coincidences' between figures 11 and 12 are that both indicate the same magnitude for R c , i.e. 1, and that the wave amplitude in both cases changes exponentially. Due to the two points, this paper argues the instability of Stokes waves subject to sidebands disturbances in oblique directions would be a possible cause for the rapid change of wave amplitude shown in figure 11 in the neighbourhood of the linear focused time t = t f .
Examining figures 10, 11 and 12 together, the following features are highlighted.
(i) The rapid energy transfers between upper and lower sidebands shown in figure 10 are likely due to the instability of Stokes waves in addition to the degenerate resonant interaction in the 'figure of eight' quartet resonant loop (Phillips 1967), as a result of that, the change rate of the amplitude shown in figure 11 is of the same order of magnitude as those in figure 12 when instability occurs. It is well understood that the instability of Stokes waves is mostly accompanied by the degenerate resonant interaction, as pointed out by Phillips (1967) where it is shown that the neutral instability modes (i.e. Im(δ Ω ) = 0) of Stokes waves can be directly derived from the degenerate resonant interaction (Hammack & Henderson 1993) and, most importantly, the instability of Stokes waves subject to oblique wave-train disturbances is likely to occur in a narrow spectrum. (ii) A downward shift of the spectral peak arises probably from the asymmetrical change rates of energy between the different regions of sidebands. This shift is consistent with the finding of Trulsen & Dysthe (1997). When a downward shift of spectral peak is observed, rapid energy transfers are shown between lower and upper sidebands in the longitudinal direction, and between the transverse and longitudinal sidebands. (iii) Oblique energy transfers towards large upper sidebands (0.3 δ x or 0.27 δ y ) are obviously seen in figure 10(c, f ), probably due to the different change rates within the instability region of Stokes waves; the contour levels of large change rates of amplitude have demonstrated that the change of energy favours more in oblique directions. This observation is consistent with Trulsen & Dysthe (1997) and the rapid energy transfer reported in Barratt et al. (2021). Barratt et al. (2021) focus on a numerical study of the propagation of a steep wave group based on direct numerical solutions of fully nonlinear potential flow equations. Barratt et al. (2021) reports an angle along ang(δ y /δ x ) = ±55 • to the spectral peak in the region of upper sidebands. This paper argues that the angle as such, ±55 • , may not exist. It is likely a coincidence of the effect shown in figure 11(d,e) that larger upper sidebands (0.3 δ x or 0.27 δ y ) tend to grow exponentially at a rate among the largest. Moreover, Barratt et al. (2021) attribute the rapid oblique energy transfers, similar to those shown in the region (i.e. δ y > δ x / √ 2 for 0.3 δ x ) of the upper sidebands in figure 10(c, f ), to non-degenerate resonant interaction. This argument is different from what is shown in figures 10, 11 and 12, i.e. that the amplitude changes exponentially is unlikely a result of resonant interactions.

Concluding remarks
Using a perturbation expansion, this paper has derived a new NLSE for the linear envelope evolution of three-dimensional surface gravity waves without the assumption of a narrow bandwidth. The new NLSE can be extended for more general cases, e.g. waves on water of finite uniform or slowly varying depth. Simple numerical implementations of the NLSE for the wave envelope using a pseudospectral, a split-step and a finite-difference method are demonstrated. The computational cost is of the same order as the numeral implementations of an existing NLSE, e.g. the Dysthe equation of Dysthe (1979) and the modified NLSE of Trulsen et al. (2000). With leading-order approximations to the NLSE, its capabilities for studying the instability region and energy growth rate of Stokes waves are demonstrated. The approximation is shown to be capable of providing satisfactory predictions of the instability region, compared with the exact computations by McLean (1982) for dimensionless wave steepness equal to as large as 0.33. It can also provide good estimates of the energy growth rate and stability boundary of unstable perturbations for two-dimensional wave trains, compared with the results due to Crawford et al. (1981) and the exact numerical results based on Longuet-Higgins (1978). In these comparisons, the approximation has been demonstrated to provide a better performance than the modified NLSE by Trulsen et al. (2000), particularly for sideband disturbances of waves in an oblique direction to a train of Stokes waves. This paper has also proposed a semianalytical approach for the description of wave fields at second order in wave steepness, based on a pseudospectral and a finite-difference method. It is demonstrated that the semianalytical approach is capable of providing exact predictions as that of Dalzell (1999). Compared with Dalzell (1999), the semianalytical approach can significantly improve the computational efficiency; the numerical implementations of Dalzell (1999) require computational operations at O(N 2 ) and the semianalytical approach requires O (N In(N)) if N Fourier modes are used. In the analysis of a velocity head at second order in wave steepness, it is suggested that the superharmonic and subharmonic waves may be equally important in directional sea states due to the same order of magnitude.
Through a study of the temporal evolution of a steep focused wave group, energy transfers in a narrow spectrum have been investigated. For the time range in the vicinity of the group at linear focus, rapid oblique energy transfers between different regions of sidebands are observed, likely arising from degenerate resonant interactions in the so-called 'figure of eight' quartet resonant loop (Phillips 1967) and the instability of Stokes waves subject to oblique sideband perturbations. The change rate of wave amplitude is found to be of the same order of magnitude as that of oblique instability of Stokes waves. A downward shift of the spectral peak is also observed in this paper, consistent with the finding in Trulsen & Dysthe (1997). This shift is observed in the course of rapid energy transfers between different sideband regions when oblique modulational instability occurs, in contrast to the temporal evolution of a unidirectional spectrum which is found to lead to a permanent shift of the spectrum peak only in a quasi-steady state in Dysthe et al. (2003). Rapid energy transfers towards the angle, ∼ ± 55 • , to the spectral peak in the region of upper sideband, that are reported in Barratt et al. (2021), are not supported by this study. Instead, the upper sideband range of large oblique directions (>45 • ) to the spectral peak is found to be favoured the most for the energy growth, coinciding with the sideband region where the instability growth rate of three-dimensional Stokes waves is among the largest. The present theory Exact (Dalzell 1999) Second-order Stokes where the relations ∂ t A = −ω 0 ∂ x A/(2k 0 ) and ∂ t ∂ x A * = −ω 0 A * A/(2k 0 ) were used in (A11c). Thus, we obtain for unidirectional waves Moreover, the second term in (4.5b) reads Combing (A12)