---
abstract: |
    The Hyperbolic Nonlinear Schrödinger equation () arises as a model for the dynamics of three–dimensional narrowband deep water gravity waves. In this study, the Petviashvili method is exploited to numerically compute bi-periodic time-harmonic solutions of the equation. In physical space they represent non-localized standing waves. Non-trivial spatial patterns are revealed and an attempt is made to describe them using symbolic dynamics and the language of substitutions. Finally, the dynamics of a slightly perturbed standing wave is numerically investigated by means a highly acccurate Fourier solver.

    **:** Hyperbolic equations; NLS equation; wave patterns; deep water waves
address:
- 'LAMA, UMR 5127 CNRS, Université de Savoie, Campus Scientifique, 73376 Le Bourget-du-Lac Cedex, France'
- 'University College Dublin, School of Mathematical Sciences, Belfield, Dublin 4, Ireland LAMA, UMR 5127 CNRS, Université de Savoie, Campus Scientifique, 73376 Le Bourget-du-Lac Cedex, France'
- 'School of Civil and Environmental Engineering & School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, USA'
author:
- Laurent Vuillon
- 'Denys Dutykh$^*$'
- Francesco Fedele
bibliography:
- 'biblio.bib'
title: Some special solutions to the Hyperbolic NLS equation
---

[^1]

Introduction
============

The celebrated cubic Nonlinear Schrödinger (NLS) equation is one of the most important equations in nonlinear science [@Sulem1999]. For example, it arises in plasma physics [@Sen1978] and in normally dispersive optical waveguide arrays modeling [@Conti2003; @Lahini2007]. In the context of water waves the equation is the leading order model of the wave envelope evolution. It was derived for the first time by V. <span style="font-variant:small-caps;">Zakharov</span> (1968) [@Zakharov1968] and rediscovered later by several other authors [@Chu1970; @Hasimoto1972]. In one dimension (1-D), the NLS equation for unidirectional water waves is integrable [@Zakharov1972] and of focusing type. As a result, a steady or periodic balance of the cubic nonlinearities and wave dispersion can be attained, and this yields the formation of localized traveling waves (solitons) or homoclinic orbits to a plane wave (breathers). Analytical solutions for solitons follow via the inverse scattering transform [@Zakharov1972; @John; @Ablowitz1979; @Pelinovsky2001a] and breathers can be easily obtained via the Darboux transformation [@Osborne2010]. The two-dimensional (2-D) propagation of deep water narrowband waves is instead governed by the 2-D Hyperbolic NLS equation [@Ablowitz1979; @Yuen1982; @Sulem1999; @Pelinovsky2003]. In this case, energy can spread along the transversal direction to the main propagation. A consequence of this defocusing is that localized traveling waves cannot occur. Indeed, their non-existence was proved in ([@Ghidaglia1996a], see also [@Ghidaglia1993]). However, this does not exclude the existence of nontrivial *nonlocalized* travelling wave patterns that may arise due to a balance between nonlinearities and wave dispersion in both directions under toric constraints.

An extended equation was derived by <span style="font-variant:small-caps;">Trulsen</span> & <span style="font-variant:small-caps;">Dysthe</span> (1996) [@Trulsen1996] and additionally to the classical hyperbolic (D’Alembert) operator, it contains also higher order dispersive and nonlinear terms. The generalization to the finite depth case leads to the Davey–Stewartson equations [@Benney1967; @Davey1974; @Ghidaglia1990]. Several important analytical solutions to the the Davey–Stewartson model were derived in [@Ebadi2011; @Ebadi2011a; @Yildirim2012; @Jafari2012].

For mathematical/numerical studies it is often convenient to restrict the attention to a particular class of solutions. For example, the very first mathematical description of plane permanent waves is known at least since G. <span style="font-variant:small-caps;">Stokes</span> (1847) [@Stokes1847]. Permanent waves can be periodic or localized in space. In this study we focus on bi-periodic time-harmonic solutions of the form $A = F(x - c_g t, y){\mathrm{e}}^{{\mathrm{i}}\omega t}$, where $A$ is the complex envelope and $F$ a complex function of its argument. These are stationary solutions in the frame of reference moving with the group speed $c_g$. In the physical domain, the associated wave surface displacements is that of standing waves, which have been the subject of many studies. In particular, radial standing solutions of the equation were investigated in [@Kevrekidis2011]. The existence of standing waves in deep waters were proved in [@Iooss2005]. In the shallow-water regime standing wave patterns have been found in the context of the Boussinesq equations by M. <span style="font-variant:small-caps;">Chen</span> & G. <span style="font-variant:small-caps;">Iooss</span> [@Chen2005; @Chen2006; @Chen2008]. 2-D bi-periodic travelling wave solutions to the Euler equations were studied by W. <span style="font-variant:small-caps;">Craig</span> & D. <span style="font-variant:small-caps;">Nicholls</span> (2002) [@Craig2002]. Recently, 2-D wave patterns of the free surface were investigated experimentally by D. <span style="font-variant:small-caps;">Henderson</span> *et al.* (2010) [@Henderson2010] along with a theoretical stability analysis.

In this work, symbolic dynamics and associated techniques of substitutions are exploited to investigate the structure of standing wave patterns (see [@Morse1940; @Lind1995; @Berthe2000]). In physics, their application in studies of dynamical systems led to unveiling the structure of quasi-periodic tilings (see [@Senechal1995]). Symbolic dynamics allows coding the nonlinear behavior of a complex system and pattern formation by means of 1-D or 2-D words of a finite alphabet. Such approach was applied to code the non-periodic trajectories on a unit circle with a particular partition on two intervals (see [@Morse1940]). Coding of finite, periodic and non-periodic infinite patterns using 2-D words and tilings was done in [@Berthe2000; @DeBruijn1981]. The numerical standing waves investigated in this work are periodic in space, thus the associated patterns are described up to toric constraints. This is the first step for developing a theory for the description of periodic or quasi-periodic patterns associated to trajectories $\{a_n(t)\}$ of the dynamics in the infinite phase space spanned by, for example, generalized Fourier basis $\phi_n(x,y)$ associated to the formal series for $A(x,y,t) = \sum_n a_n(t)\phi_n(x,y)$ on a periodic domain.

The present study is organized as follows. First, the equation is introduced in the context of deep water waves. Then, the Petviashvili method [@Petviashvili1976] used to compute a class of special standing wave solutions is presented. Symbolic dynamics is then applied to describe the associated spatial patterns using words and substitutions. Finally, a highly accurate Fourier-based solver is exploited to investigate the dynamics of a perturbed standing wave.

The Mathematical Model {#sec:math}
======================

Consider a three-dimensional fluid domain with a free surface. The water is assumed to be infinitely deep. The Cartesian coordinate system $Oxyz$ is chosen such that the undisturbed water level corresponds to $z = 0$, and the free surface elevation is $z = \eta(x,y,t)$.

The Euler equations that describe the irrotational flow of an ideal incompressible fluid of infinite depth with a free surface are of fundamental relevance in fluid mechanics, ocean sciences and both pure and applied mathematics (see for example, [@Stoker1957; @Zakharov1968; @Johnson1997]). The structure of the Euler equations is given in terms of the free-surface elevation $\eta(x, y, t)$ and the velocity potential $\varphi(x, y, t) = \phi(x, y, z = \eta(x,y,t), t)$ evaluated at the free surface of the fluid. In late 70s Dysthe used the method of multiple scales to derive from the Euler equations a modified Nonlinear Schrödinger (NLS) equation [@Dysthe1979] for the time evolution of the unidirectional narrowband envelope $A$ of the velocity potential $\varphi$ with carrier wave ${\mathrm{e}}^{{\mathrm{i}}(k_0 x - \omega_0 t)}$, $k_0$ and $\omega_0$ being, respectively, the wavenumber and frequency of the carrier wave. The equation for $A$ can be formulated in a frame moving with the group velocity $c_g = \omega_0/(2k_0)$ as follows. Define ${\varepsilon}$ as a small parameter, $a_0$ as a characteristic wave amplitude and rescale space, time and the envelope as $x \to k_0 x-c_g t$, $y \to k_0 y$, $t \to \omega_0 t$ and $A \to {\varepsilon}a_0 A$ respectively. Then, the 2-D Dysthe equation for $A$ is given by [@Trulsen1996]: $$\begin{gathered}
\label{eq:Dysthe}
  {\mathrm{i}}A_t = \frac{1}{8}A_{xx} - \frac{1}{4}A_{yy} + \frac{1}{2}|A|^2A \\ - {\mathrm{i}}{\varepsilon}\bigl(\frac{1}{16}A_{xxx} - \frac{3}{8}A_{xyy} - \frac{3}{2}|A|^2A_x + \frac{1}{4}A^2A^{\star}_x - {\mathrm{i}}A{\mathcal{H}}[|A|^2]\bigr),\end{gathered}$$ where ${\mathcal{H}}$ is the Hilbert transform and the subscripts $A_t = \partial_t A$, $A_x = \partial_x A$ denote partial derivatives with respect to $x$, $y$ and $t$ respectively, and $A^{\star}$ denotes complex conjugation. The envelope of the free surface $\eta$ relates to $A$ by a simple transformation that involves only $A$ and its derivatives (see [@Dysthe1979]). To ${\mathcal{O}}(1)$ in ${\varepsilon}$, the 2-D Dysthe equation reduces to the equation: $$\label{eq:hNLS}
  {\mathrm{i}}A_t = \frac{1}{8}A_{xx} - \frac{1}{4}A_{yy} + \frac{1}{2}|A|^2A.$$ By rescaling $A \to 2A$, $x \to 2\sqrt{2} x$ and the transverse coordinate $y \to 2y$, takes the form [@Sulem1999; @Pelinovsky2003]: $$\label{eq:adim}
  {\mathrm{i}}A_t = A_{xx} - A_{yy} + 2A|A|^2.$$ This equation admits the three invariants ${\mathbb{H}}$, ${\mathbb{A}}$ and ${\mathbb{M}}$, which have the meaning of energy, wave action and momentum respectively: $$\begin{aligned}
  {\mathbb{H}}&=& \iint_{{\mathds{R}}^2}\Bigl\{|A_y|^2 - |A_x|^2 + |A|^4\Bigr\}\,{\mathrm{d}}x\,{\mathrm{d}}y, \\
  {\mathbb{A}}&=& \iint_{{\mathds{R}}^2}|A|^2\,{\mathrm{d}}x\,{\mathrm{d}}y, \\
  {\mathbb{M}}&=& \frac{{\mathrm{i}}}{2}\iint_{{\mathds{R}}^2}\Bigl\{A{\boldsymbol{\nabla}}A^\star + A^\star{\boldsymbol{\nabla}}A\Bigr\}\,{\mathrm{d}}x\,{\mathrm{d}}y.\end{aligned}$$ Note that the total energy ${\mathbb{H}}$ is also the Hamiltonian for the equation. In the following, we will solve for a special class of standing wave solutions to .

Standing wave patterns {#sec:gs}
----------------------

Consider the ansatz for standing wave solutions that oscillate harmonically in time $$A(x,y,t) = {\mathrm{e}}^{-{\mathrm{i}}\omega t}B(x,y),$$ where the real function $B$ describes the spatial pattern of a standing wave. According to the unscaled equation , $B$ satisfies the following real nonlinear hyperbolic PDE $$\label{eq:Beq}
  \omega B + B_{xx} - B_{yy} = 2B^3,$$ which can be solved numerically using the classical Petviashvili method [@Petviashvili1976; @Pelinovsky2004; @Lakoba2007]. To do so, is rewritten in the operator form $${\mathcal{L}}\cdot B = {\mathcal{N}}(B), \qquad {\mathcal{L}}:= \omega + \partial_{xx} - \partial_{yy}, \qquad {\mathcal{N}}(B) := 2B^3.$$ and the iteration scheme is given by $$B_{n+1} = {\ensuremath{\mathcal{S}}}^\gamma {\mathcal{L}}^{-1}\cdot {\mathcal{N}}(B_n), \quad
  {\ensuremath{\mathcal{S}}}= \frac{\langle B_n, {\mathcal{L}}\cdot B_n\rangle}{\langle B_n, {\mathcal{N}}(B_n)\rangle},$$ where ${\ensuremath{\mathcal{S}}}$ is the so-called stabilyzing factor and the exponent $\gamma$ is usually defined as a function of the degree of nonlinearity $p$ ($p = 3$ for the equation). The rule of thumb prescribes the following formula $\gamma = \frac{p}{p-1}$. The scalar product is defined in the $L_2$ space. The inverse operator ${\mathcal{L}}^{-1}$ can be efficiently computed in Fourier space using the Fast Fourier Transform (FFT) (see, for example, [@Frigo2005]), and spatial periodicity across the 2D computational box boundaries is implicitly imposed. The iterative process converges for a large class of smooth initial guesses[^2]. Convergence is attained when the $L_\infty$ norm between two successive iterations is less than a prescribed tolerance $\varepsilon$ (usually of the order of machine precision). Additionally, the residual error $E_r$ in approximating the nonlinear equation is checked by substituting in the converged numerical solution. In the present work, the convergence of the algorithm is checked numerically in the extended floating point arithmetics using $30$ significant digits [@MATLAB2012].

For example, consider the domain $\Omega := [-\ell_x/2, \ell_x/2] \times [-\ell_y/2, \ell_y/2] \subset {\mathds{R}}^2$ to be a square with the side length equal to $210$ ($\Omega = [-105, 105]^2$). For $\omega = 0.012$, the Petviashvili method on a $1024 \times 1024$ Fourier grid yields the strictly periodic regular pattern shown in Figure \[fig:f0\](a). The convergence is attained in $N\sim 80$ iterations as clearly seen in Figure \[fig:f0\](b) and the associated error $L_\infty \sim {\mathcal{O}}(10^{-33})$, and the residual $E_r \sim {\mathcal{O}}(10^{-32})$. Note that $L_\infty \sim~{\mathrm{e}}^{-0.095 N}$ decays exponentially in agreement with the theoretical geometric convergence rate [@Alvarez2014]. For the sake of efficiency, the numerical solutions presented below are computed on the same Fourier grid using standard double-precision arithmetics. It is verified that residual errors are within the prescribed tolerance parameter $\varepsilon \sim 10^{-15}$ and never exceed $10\varepsilon$.

![*A bi-periodic wave pattern $B(x,y)$, $\Omega = [-30, 30]\times[-15, 15]$ and $\omega = 0.18$. $B$ can be described by the $1\times 4$ discrete pattern $E = {}^t(r, b, r, b)$ on the two-letter alphabet $\{r, b\}$ (red spots code negative values far from zero and blue spots refer to positive values far from zero).*[]{data-label="fig:f1"}](figs/gstates/60x30_0_18_select.eps){width="75.00000%"}

![*A bi-periodic wave pattern $B(x,y)$, $\Omega = [-30, 30]^2$ and $\omega = 0.18$. The box delimits the elementary $6\times 10$ discrete pattern of blue ($b$), red ($r$) and white ($w$, denoting the values around zero) spots identified to describe $B$ by substitutions on the three-letter alphabet $\{r, b, w\}$. For symbolic coding see caption of Figure \[fig:f0\].*[]{data-label="fig:f2"}](figs/gstates/60x60_0_18_select.eps){width="75.00000%"}

[figs/gstates/60x40\_0\_18\_select.eps]{} (24,54.5) [$\sigma(b)$]{} (43,54.5) [$\sigma(a)$]{} (24,47.5) [$\sigma(c)$]{} (43,47.5) [$\sigma(d)$]{} (24,39.5) [$\sigma(\bar{b})$]{} (43,39.5) [$\sigma(\bar{a})$]{} (24,30) [$\sigma(\bar{a})$]{} (43,30) [$\sigma(\bar{b})$]{} (24,22) [$\sigma(d)$]{} (43,22) [$\sigma(c)$]{} (24,14) [$\sigma(a)$]{} (43,14) [$\sigma(b)$]{}

![*A bi-periodic wave pattern $B(x,y)$, $\Omega = [-30, 30]\times[-60, 60]$ and $\omega = 0.18$. The box delimits the elementary $32\times 66$ discrete pattern of blue ($b$) red ($r$) and white ($w$) spots identified to describe $B$ by substitutions on the three-letter alphabet $\{r, b, w\}$. For symbolic coding see Figure \[fig:f0\].*[]{data-label="fig:f4"}](figs/gstates/60x120_0_18_select.eps){width="99.00000%"}

![*A bi-periodic wave pattern $B(x,y)$, $\Omega = [-\frac{55}{2}, \frac{55}{2}]\times[-\frac{89}{2}, \frac{89}{2}]$ and $\omega = 0.51$.*[]{data-label="fig:f5"}](figs/gstates/55x89_0_51.eps){width="99.00000%"}

![*(Right) a bi-periodic wave pattern $B(x,y)$, $\Omega = [-105, 105]^2$ and $\omega = 1.3295$; (Left) zoom on the sub-region $[-40, -10] \times [30, 80]$.*[]{data-label="fig:f6"}](figs/gstates/210x210_1_3295_select.eps){width="99.00000%"}

Symbolic coding {#sec:descr}
===============

The theory of bi-dimensional words and tilings (see [@Berthe2000; @DeBruijn1981]) can be exploited to describe standing wave patterns of the equation. The key tool in tiling theory is the substitution operation that consists in replacing letters by group of letters and constructing iteratively words on a finite alphabet. For example, the Fibonacci substitution $\sigma(a) = ab$ and $\sigma(b) = a$ constructs a fixed point (limit of infinite substitutions) with a non-periodic structure ([@Arnoux1991; @Lind1995; @Allouche2003]. This word constitutes a model for 1-D quasicrystals in physics (see [@Senechal1995]) and a coding of a discrete line with irrational slope in discrete geometry (see [@Vuillon2003]).

The method of substitutions empowers the construction of either non-periodic or periodic words $G$ as fixed points of the mapping $G=\sigma^n(G)$ as $n\to +\infty$ (see [@Lind1995]). On the one hand, non-periodic infinite words can be constructed iteratively by the use of morphism properties. In particular, a substitution $\sigma$ applied to a word $w=w_1w_2\cdots w_n$ is in fact the image of the substitution applied to each letter of $w$, that is $\sigma(w_1w_2\cdots w_n) = \sigma(w_1) \sigma(w_2) \cdots \sigma(w_n)$, where $w_i$ are letters of the alphabet. As an example, the first few iterations on the letter $a$ using the Fibonacci substitution are given by $\sigma^1(a) = \sigma(a) = ab$, $\sigma^2(a) = \sigma(\sigma(a)) = \sigma(ab) = \sigma(a) \sigma(b) = aba$, $\sigma^3(a) = \sigma(\sigma^2(a)) = \sigma(aba) = \sigma(a)\sigma(b)\sigma(a) = abaab$, and so on (see [@Arnoux1991; @Vuillon2003]). Note that each word $\sigma^{n}(a)$ is the beginning of the next word $\sigma^{n+1}(a)$. Thus, as $n\to +\infty$, a non-periodic fixed point is constructed, viz. $$\sigma(F) = F =  
  abaababaabaababaababaabaababaabaababaababaabaababaababaabaababaaba\cdots.$$

On the other hand, a periodic word can be constructed by means of the substitution $\sigma(a)=ab$ and $\sigma(b)=ab$ by repetition of the pattern $ab.$. Indeed $\sigma^1(a)=\sigma(a)=ab, $$\sigma^2(a) = \sigma(\sigma(a)) = \sigma(ab) = \sigma(a)\sigma(b) = abab$, $\sigma^3(a) = \sigma(\sigma^2(a)) = \sigma(abab) = \sigma(a)\sigma(b)\sigma(a)\sigma(b) = abababab,$ and so on. The fixed point is periodic and given by $$\sigma(P) = P = 
  abababababababababababababababababababababababababababa\cdots.$$

In order to code periodic spatial patterns of the equation, one can exploit bi-dimensional words (see [@Frougny2005; @Berthe2000; @Lind1995]) and proper substitutions (see [@Arnoux1991; @Berthe2000; @Frougny2005]). For example, the bi-dimensional Thue–Morse substitution (see [@Allouche2003]) is defined by $$\sigma(a) =
\begin{pmatrix}
  a & b \\ 
  b & a
\end{pmatrix}
\quad \mbox{and} \quad 
\sigma(b) = 
\begin{pmatrix}
  b & a \\ 
  a & b 
\end{pmatrix}.$$

One can just iterate a finite number of substitutions, viz. $\sigma^j(a)$ with $j$ finite, and obtain a finite world. For example, the first two iterations of the Thue-Morse substitution yield, respectively, the two words $$\sigma^1(a) = \sigma(a) = 
\begin{pmatrix}
  a & b \\ 
  b & a
\end{pmatrix}$$ and $$\sigma^2(a) = \sigma
\begin{pmatrix}
  a & b \\ 
  b & a
\end{pmatrix} = 
\begin{pmatrix}
  \sigma(a) & \sigma(b) \\ 
  \sigma(b) & \sigma(a)
\end{pmatrix} = 
\begin{pmatrix}
  a & b & b & a \\
  b & a & a & b \\
  b & a & a & b \\
  a & b & b & a
\end{pmatrix}$$

In the limit of infinite iterations, an infinite word is obtained in the form of the non-periodic fixed point $$\sigma(M) = M = 
\begin{pmatrix}
  a & b & b & a & b & a \cdots \\
  b & a & a & b & a & b \cdots \\
  b & a & a & b & a & b \cdots \\
  a & b & b & a & b & a \cdots \\
  b & a & a & b & a & b \cdots \\
  a & b & b & a & b & a \cdots \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots ~ \ddots
\end{pmatrix}$$

Hereafter, a discrete description of several spatial patterns $B(x,y)$ of the equation obtained via the Petviashvili method is presented. Each continuous pattern is coded on a finite alphabet using proper substitutions iterated on a discrete elementary pattern of letters.

Consider the bi-periodic pattern $B$ of the equation shown in Figure \[fig:f0\]. This can be easily described by repetition of the elementary 3x3 discrete pattern delimited by a box in the same Figure. To do so, coding with a three–letter alphabet is used since the pattern is characterized by a discrete structure with three different elementary spots. The letters $r$,$b$ and $w$ are used to code red spots (negative values far from zero), blue spots (positive values far from zero) and white spots (values around zero) respectively.

As clearly seen in Figure \[fig:f0\], the three spots are arranged on a rectangular sub-region of the ${\mathbb{Z}}^2$ grid. Thus, the continuous periodic pattern $B$ can be easily described by bi-dimensional words constructed using the two simple substitutions $$\sigma(r) = 
\begin{pmatrix}
  r & b \\ 
  b & r
\end{pmatrix} \quad \mbox{and}\quad \sigma(b) = 
\begin{pmatrix}
 r & b \\ 
 b & r 
\end{pmatrix}.$$

By iterating the above substitutions yields the strictly periodic word $$\sigma(R) = R = 
\begin{pmatrix}
  r & b & r & b \cdots\\
  b & r & b & r \cdots\\
  r & b & r & b \cdots \\
  b & r & b & r \cdots\\
  \vdots & \vdots & \vdots & \vdots ~ \ddots
\end{pmatrix}$$ which describes the continuous wave pattern of Figure \[fig:f0\].

Consider now the wave pattern solutions for $B(x,y)$ obtained for $\omega = 0.18$ on various rectangular grid sizes. For $\Omega = [-30, 30]\times [-15, 15]$, $B$ can be described by a $1 \times 4$ discrete pattern delimited by a box in Figure \[fig:f1\] and coded as $E = {}^t(r, b, r, b)$, where the superscript ${}^t$ denotes matrix transposition. $E$ can be decomposed as $\sigma(e) = {}^{t}(r, b)$ and the initial pattern $E$ is given by the substitution $\sigma$ applied to $E = {}^t(e, e)$ that is: $$\sigma(E) = {}^t\bigl(\sigma(e), \sigma(e)\bigr) = {}^t\bigl(r, b, r, b\bigr).$$

For $\Omega = [-30, 30]\times [-30, 30]$, the associated continuous pattern $B$ is shown in Figure \[fig:f2\]. A box delimits the elementary $6 \times 10$ discrete pattern $F$ that can describe $B$ by successive substitutions, viz. $$F = \begin{pmatrix}
  b & w & b & w & b & w \\
  r & w & r & w & r & w \\
  b & r & b & r & b & r \\
  r & b & r & b & r & b \\
  w & r & w & r & w & r \\
  w & b & w & b & w & b \\
  w & r & w & r & w & r \\
  r & b & r & b & r & b \\
  b & r & b & r & b & r \\
  r & w & r & w & r & w
\end{pmatrix}.$$

This discrete pattern can be decomposed as $$\sigma(f) = 
{}^t\begin{pmatrix}
  w & w & r & b & r & b & r & b & r & w \\
  b & r & b & r & w & w & w & r & b & r
\end{pmatrix}.$$ As a result, the initial pattern $F$ is given by the substitution $\sigma$ applied to $F = (f~f~f)$, that is $$\sigma(F) = \bigl(\sigma(f) ~ \sigma(f) ~ \sigma(f)\bigr) =
\begin{pmatrix}
  b & w & b & w & b & w \\
  r & w & r & w & r & w \\
  b & r & b & r & b & r \\
  r & b & r & b & r & b \\
  w & r & w & r & w & r \\
  w & b & w & b & w & b \\
  w & r & w & r & w & r \\
  r & b & r & b & r & b \\
  b & r & b & r & b & r \\
  r & w & r & w & r & w
\end{pmatrix}.$$

For the larger domain $\Omega = [-30, 30]\times [-20, 20]$, the associated standing wave pattern $B(x,y)$ is reported in Figure \[fig:f3\]. An elementary $32\times 22$ discrete pattern of blue ($b$), red ($r$) and white ($w$) spots is identified and delimited by a box to describe $B$. This can be obtained by four iterations of the following 6 substitutions on the three-letter alphabet $\{r, b, w\}$, that is $$\sigma(a) =
\begin{pmatrix}
  r & w & r & w & r & w & r & w \\
  w & r & w & r & w & r & w & w \\
  r & w & r & w & r & w & r & w \\
  w & w & w & w & w & w & w & w
\end{pmatrix},$$ $$\sigma(\bar{a}) =
\begin{pmatrix}
  w & r & w & r & w & r & w & w \\
  r & w & r & w & r & w & r & w \\
  w & r & w & r & w & r & w & w \\
  w & w & w & w & w & w & w & w
\end{pmatrix},$$ $$\sigma(b) =
\begin{pmatrix}
  w & b & w & b & w & b & w & w \\
  b & w & b & w & b & w & b & w \\
  w & b & w & b & w & b & w & w \\
  w & w & w & w & w & w & w & w
\end{pmatrix},$$ $$\sigma(\bar{b}) =
\begin{pmatrix}
  b & w & b & w & b & w & b & w \\
  w & b & w & b & w & b & w & w \\
  b & w & b & w & b & w & b & w \\
  w & w & w & w & w & w & w & w
\end{pmatrix},$$ $$\sigma(c) =
\begin{pmatrix}
  r & w & r & w & r & w & r & w \\
  w & r & w & r & w & r & w & w \\
  r & w & r & w & r & w & r & w
\end{pmatrix}$$ $$\sigma(d) =
\begin{pmatrix}
  w & b & w & b & w & b & w & w \\
  b & w & b & w & b & w & b & w \\
  w & b & w & b & w & b & w & w
\end{pmatrix}$$ As described above, this substitution is applied four times under toric constraints starting from the finite $4 \times 6$ pattern $G$ given by $$G =
\begin{pmatrix}
  b & a & b & a \\
  c & d & c & d \\
  \bar{b} & \bar{a} & \bar{b} & \bar{a} \\
  \bar{a} & \bar{b} & \bar{a} & \bar{b} \\
  d & c & d & c \\
  a & b & a & b
\end{pmatrix}.$$ Here, the difference between $\sigma(a)$ and $\sigma(\bar{a})$ resides in the interchange of the letters $r$ and $w$ in the first 3 rows of the pattern except the last column, and similarly for $\sigma(b)$ and $\sigma(\bar{b})$. This switching between components is exactly the typical dynamical property that we expect for symbolic coding of patterns. For example, the first iteration $\sigma(G)$ yields $$\sigma(G) = \sigma\left(\begin{array}{llllllllllllllll}
  b & a & b & a \\
  c & d & c & d \\
  \bar{b} & \bar{a} & \bar{b} & \bar{a} \\
  \bar{a} & \bar{b} & \bar{a} & \bar{b} \\
  d & c & d & c \\
  a & b & a & b
\end{array}\right) =
\begin{pmatrix}
\sigma(b) & \sigma(a) & \sigma(b) & \sigma(a) \\
\sigma(c) & \sigma(d) & \sigma(c) & \sigma(d) \\
\sigma(\bar{b}) & \sigma(\bar{a}) & \sigma(\bar{b}) & \sigma(\bar{a}) \\
\sigma(\bar{a}) & \sigma(\bar{b}) & \sigma(\bar{a}) & \sigma(\bar{b}) \\
\sigma(d) & \sigma(c) & \sigma(d) & \sigma(c) \\
\sigma(a) & \sigma(b) & \sigma(a) & \sigma(b)
\end{pmatrix}.$$ The whole discrete description of the continuous pattern $B(x,y)$ follows after four iterations as $$\sigma(G) =
\left(\begin{array}{lllllllllllllllllllllllllllllllllll}
  w & b & w & b & w & b & w & w & r & w & r & w & r & w & r & w \cdots \\
  b & w & b & w & b & w & b & w & w & r & w & r & w & r & w & w \cdots \\
  w & b & w & b & w & b & w & w & r & w & r & w & r & w & r & w \cdots \\
  w & w & w & w & w & w & w & w & w & w & w & w & w & w & w & w \cdots \\
  r & w & r & w & r & w & r & w & w & b & w & b & w & b & w & w \cdots \\
  w & r & w & r & w & r & w & w & b & w & b & w & b & w & b & w \cdots \\
  r & w & r & w & r & w & r & w & w & b & w & b & w & b & w & w \cdots \\
  b & w & b & w & b & w & b & w & w & r & w & r & w & r & w & w \cdots \\
  w & b & w & b & w & b & w & w & r & w & r & w & r & w & r & w \cdots \\
  b & w & b & w & b & w & b & w & w & r & w & r & w & r & w & w \cdots \\
  w & w & w & w & w & w & w & w & w & w & w & w & w & w & w & w \cdots \\
  w & r & w & r & w & r & w & w & b & w & b & w & b & w & b & w \cdots \\
  r & w & r & w & r & w & r & w & w & b & w & b & w & b & w & w \cdots \\
  w & r & w & r & w & r & w & w & b & w & b & w & b & w & b & w \cdots \\
  w & w & w & w & w & w & w & w & w & w & w & w & w & w & w & w \cdots \\
  \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots ~ \ddots
\end{array}\right)$$

For $\Omega = [-30, 30]\times [-60, 60]$, the associated water pattern $B$ can be described by the elementary $32 \times 66$ discrete pattern $K$ delimited by a box and shown in Figure \[fig:f4\]. In this case, the symbolic description of $K$ follows from the substitution $K=\sigma(H)$, where $\sigma$ is defined as $$\sigma(a) =
\begin{pmatrix}
  r & w & r & w & r & w & r & w \\
  w & r & w & r & w & r & w & w \\
  r & w & r & w & r & w & r & w \\
  w & w & w & w & w & w & w & w
\end{pmatrix},$$ $$\sigma(\bar{a}) =
\begin{pmatrix}
  w & r & w & r & w & r & w & w \\
  r & w & r & w & r & w & r & w \\
  w & r & w & r & w & r & w & w \\
  w & w & w & w & w & w & w & w
\end{pmatrix},$$ $$\sigma(b) =
\begin{pmatrix}
  w & b & w & b & w & b & w & w \\
  b & w & b & w & b & w & b & w \\
  w & b & w & b & w & b & w & w \\
  w & w & w & w & w & w & w & w
\end{pmatrix},$$ $$\sigma(\bar{b}) =
\begin{pmatrix}
  b & w & b & w & b & w & b & w \\
  w & b & w & b & w & b & w & w \\
  b & w & b & w & b & w & b & w \\
  w & w & w & w & w & w & w & w
\end{pmatrix},$$ $$\sigma(c) =
\begin{pmatrix}
  r & w & r & w & r & w & r & w \\
  w & r & w & r & w & r & w & w \\
  r & w & r & w & r & w & r & w
\end{pmatrix}$$ $$\sigma(d) =
\begin{pmatrix}
  w & b & w & b & w & b & w & w \\
  b & w & b & w & b & w & b & w \\
  w & b & w & b & w & b & w & w
\end{pmatrix}$$ and applied to the following $4 \times 12$ pattern $H$ under toric constraints, viz. $$H =
\begin{pmatrix}
  G \\
  G
\end{pmatrix} = {}^t\left(
\begin{array}{llllllllllll}
  a & d & \bar{a} & \bar{b} & c & b & a & d & \bar{a} & \bar{b} & c & b \\
  b & c & \bar{b} & \bar{a} & d & a & b & c & \bar{b} & \bar{a} & d & a \\
  a & d & \bar{a} & \bar{b} & c & d & a & d & \bar{a} & \bar{b} & c & b \\
  b & c & \bar{b} & \bar{a} & d & a & b & c & \bar{b} & \bar{a} & d & a
\end{array}\right)$$

Note that $H$ is made of two copies of the discrete structure $G$ identified for the pattern $B(x,y)$ relative to the domain $\Omega = [-30, 30]\times [-20, 20]$, see Figure \[fig:f3\]. Thus, $\sigma(H)$ is exactly the substitution applied to two copies of $G$, viz. $\sigma(H) = {}^t\bigl(\sigma(G), \sigma(G)\bigr)$. This suggests that for a given value of $\omega$ the patterns associated to various domains $\Omega$ have in common the same substitutions. These capture the structure of the discrete patterns and the nature of the dynamical system.

Future research aims at finding a sequence of increasing domain sizes whose associated discrete patterns share the same substitution. This leads to construct periodic orbits and fixed points of the dynamics. For example, Figure \[fig:f5\] shows a patterns that may be generated by substitutions of elementary cell patterns. Moreover, for given domain size, one could explore which frequencies $\omega$ gives the same kind of substitutions and try to explain this regularity by the arithmetic nature of $\omega$. More precisely, suppose that for a given $\omega$ one identifies two discrete patterns generated by repeated substitutions $\sigma^k(g)$ and $\sigma^\ell(g)$, respectively, with $k$ and $\ell$ as the number of repeations. Then, one may find a sequence of increasing domain sizes $(\Omega_i)_{i \in I}$ and try to associate substitutive patterns of the form $\sigma^i(g).$ This could yield characteristic scales with invariance of patterns and could suggest that the dynamics is given by a coding of iterated substitutions of the form $(\sigma^i(g)))_{i \in I}$. A fixed point $G$ such that $\sigma(G) = G$ could then exist (see [@Morse1940; @Lind1995]). Furthermore, note that there is a link between the complexity of the patterns and the decomposition in prime factors of the domain $\Omega$ size. For example, the domain $\Omega = [-\ell/2, \ell/2] \times [-\ell/2, \ell/2]$ with length $\ell = 2 \times 3 \times 5 \times 7 = 210$ is the product of the four first prime numbers and the associated complex pattern can be seen in Figure \[fig:f6\]. The apparent complexity of the solution may be explained by a combinations of well chosen substitutions. One expects that the larger domain $\Omega$ with length $\ell = 2 \times 3 \times 5 \times 7 \times 11 = 2310$ leads to even more complicated pattern solutions. Indeed, this decomposition in prime factors increases the number of divisors of the pattern size and the complexity of patterns. To summarize, this part presents the first step for a symbolic dynamic coding of patterns arising from non linear waves. We will use in next works these substitutions in order to understand the complexity of such patterns.

Dynamics of pertubed standing waves {#sec:descr}
===================================

Pseudo-spectral scheme
----------------------

A highly accurate Fourier-type pseudo-spectral method [@Trefethen2000; @Boyd2000] will be used to solve the unsteady equation in order to investigate the dynamics of a perturbed standing wave pattern. To do so, equation is recast in the following form: $$\label{eq:oper}
  A_t + {\mathrm{i}}{\mathcal{L}}\cdot A = {\mathcal{N}}(A),$$ where the operators ${\mathcal{L}}$ and ${\mathcal{N}}$ are defined such as: $${\mathcal{L}}:= \partial_{xx} - \partial_{yy}, \qquad {\mathcal{N}}(A) = -2{\mathrm{i}}|A|^2A.$$ With this setting, is discretized by applying the 2-D Fourier transform in the spatial variables $(x, y)$. The nonlinear terms are computed in physical space, while spatial derivatives are computed spectrally in Fourier space. The standard $3/2$ rule is applied for anti-aliasing [@Trefethen2000; @Clamond2001; @Fructus2005]. The transformed variables will be denoted by $\hat{A}(t, {\boldsymbol{k}}) = {\mathcal{F}}\{A(t, x, y)\}$, with ${\boldsymbol{k}}= (k_x, k_y)$ being the Fourier transform parameter.

In order to improve the stability of the time discretization procedure, the linear part of the operator is integrated exactly by a change of variables [@Milewski1999; @Fructus2005] that yields $$\hat{{\mathcal{A}}}_t = {\mathrm{e}}^{(t-t_0){\mathcal{L}}}\cdot{\mathcal{N}}\Bigl\{{\mathrm{e}}^{-(t-t_0){\mathcal{L}}}\cdot\hat{{\mathcal{A}}}\Bigr\}, \qquad \hat{{\mathcal{A}}}(t) := {\mathrm{e}}^{(t-t_0){\mathcal{L}}}\cdot\hat{A}(t), \qquad \hat{{\mathcal{A}}}(t_0) = \hat{A}(t_0).$$ The exponential matrix $\hat{{\mathcal{L}}}$ is explicitly computed in Fourier space as $${\mathrm{e}}^{(t-t_0)\hat{{\mathcal{L}}}} = {\mathrm{e}}^{-{\mathrm{i}}(k_x^2 - k_y^2)(t-t_0)}.$$ Finally, the resulting system of ODEs is discretized in space by the Verner’s embedded adaptive 9(8) Runge–Kutta scheme [@Verner1978]. The step size is chosen adaptively using the so-called H211b digital filter [@Soderlind2003; @Soderlind2006] to meet the prescribed error tolerance, set as of the order of machine precision.

Numerical results
-----------------

Consider the bi-periodic pattern computed for $\Omega=[-30, 30]^2$ and $\omega = 0.30$ and shown in Figure \[fig:f30\]. To simulate the dynamics the pseudo-spectral method described above is used with $256\times 256$ Fourier modes. The numerical solver independently confirmed that the time-harmonic standing wave $A(x,y,t) = B(x,y){\mathrm{e}}^{-{\mathrm{i}}\omega t}$ associated to the spatial pattern $B(x,y)$ of Figure \[fig:f30\] computed using the Petviashvili scheme is effectively a solution of the original equation . Then, another experiment was performed where the initial condition for the solver was set as $A(x,y,0) = B(x,y) + w(x,y)$, where $w(x,y)$ is approximatively a $7\%$ double-periodic perturbation with the wavelength four times smaller than that of the unperturbed pattern $B(x,y)$. The simulations were carried out up to the dimensionless time $T = 14.0$. The energy (Hamiltonian) ${\mathbb{H}}$ and action ${\mathbb{A}}$ were conserved with 12 digits accuracy during the whole simulation. The momentum ${\mathbb{M}}$ was preserved to machine precision. On short time scales slight oscillations around the unperturbed solution occur. However, on a much longer time scale a transition to another solution is observed, which is quite similar in shape to the $B(x,y)$, but slightly shifted in space. A few snapshots taken from the dynamical simulation are depicted in Figure \[fig:dyn\] (see also video [@Vuillon2014]).

To visualize the dynamics, $A$ is projected onto the subspace ${\ensuremath{\mathcal{S}}}= \operatorname{span}\bigl\{\phi_1, \phi_2, \phi_3\bigr\}$ spanned by the first three leading Karhunen–Loève (KL) eigenmodes $\phi_j$, $j=1,\ldots, 3$ (see, for example [@Ghanem2003]). These are estimated from the numerical simulations using the method of snapshots (see [@Sirovich1987; @Rowley2005]) after the time average is removed. The associated trajectory $\gamma(t)$ in ${\ensuremath{\mathcal{S}}}$ and the three KL modes are shown in Figure \[fig:dyn\]. The first two modes represent the most energetic structures of the imposed perturbation, whereas the 3[$\,^{\mathrm{\small\textsf{rd}}}$]{} mode arises due to the nonlinear interaction between the unperturbed standing wave and the perturbation. Note that $\gamma$ lies approximately on a cylindrical manifold. The motion is circular on the $x_1$–$x_2$ plane with oscillations in the vertical $x_3$ axis. In physical space the dynamical wave patterns smoothly vary between the 1[$\,^{\mathrm{\small\textsf{st}}}$]{} and 2[$\,^{\mathrm{\small\textsf{nd}}}$]{} KL mode in a periodic fashion while being modulated by the 3[$\,^{\mathrm{\small\textsf{rd}}}$]{} mode. A new dynamical state is reached, which is not a standing wave.

![*A bi-periodic wave pattern, $\Omega = [-30, 30]^2$ and $\omega = 0.30$.*[]{data-label="fig:f30"}](figs/gstates/60x60_0_30.eps){width="99.00000%"}

Conclusions {#sec:concl}
===========

In this study a class of special solutions to the hyperbolic NLS equation have been investigated. In particular, bi-periodic standing waves are obtained numerically using the iterative Petviashvili scheme [@Petviashvili1976; @Pelinovsky2004; @Lakoba2007]. Non-trivial wave patterns are revealed by varying the computational domain $\Omega$ and the frequency $\omega$ of the standing wave. These are described by means of symbolic dynamics and the language of substitutions. For given value of $\omega$, the patterns associated to different domains $\Omega$ have in common the same substitution rule $\sigma$. The dynamics of a perturbed standing wave is also numerically investigated by means a highly acccurate Fourier solver in the reduced state space $S$ spanned by the first three dominant KL eigenmodes. The trajectory in $S$ lies approximately on a cylindrical manifold and in physical space the wave pattern appears to vary both in space and time. The discrete symbolic construction is the first step for developing a whole theory of description of periodic patterns of the equations. The theoretical explanation of steady and unsteady wave patterns via coding remains a major challenge and also a perspective opened by this study.

Acknowledgements {#acknowledgements .unnumbered}
================

D. <span style="font-variant:small-caps;">Dutykh</span> acknowledges the support from ERC under the research project ERC-2011-AdG 290562-MULTIWAVE. The authors would like to thank Angel <span style="font-variant:small-caps;">Duran</span>[^3] for helpful discussions on the Petviashvili method and Pavel <span style="font-variant:small-caps;">Holodoborodko</span> for providing us with the MC Toolbox for Matlab.

[^1]: $^*$ Corresponding author

[^2]: In our computations we just took an initial localized bump in the center of the domain and the iterative method turned it into a pattern.

[^3]: University of Valladolid, Spain
