Provable bounds for the Korteweg-de Vries reduction in multi-component Nonlinear Schrodinger Equation

TL;DR

This paper rigorously derives bounds and conditions under which the multi-component VNLS can be approximated by a KdV system, providing insights into hydrodynamic effects in Bose gases and related fields.

Contribution

It offers a rigorous linear analysis and a perturbative reduction of VNLS to KdV, clarifying the connection and conditions for validity in multi-component systems.

Findings

Rigorous theorems on eigenvalues and eigenvectors of the linear problem.

Validation of KdV approximation through soliton solutions.

Relevance to Bose gases, nonlinear optics, and related fields.

Abstract

We study the dynamics of multi-component Bose gas described by the Vector Nonlinear Schr\"{o}dinger Equation (VNLS), aka the Vector Gross--Pitaevskii Equation (VGPE) . Through a Madelung transformation, the VNLS can be reduced to coupled hydrodynamic equations in terms of multiple density and velocity fields. Using a multi-scaling and a perturbation method along with the Fredholm alternative, we reduce the problem to a Korteweg de-Vries (KdV) system. This is of great importance to study more transparently, the obscure features hidden in VNLS. This ensures that hydrodynamic effects such as dispersion and nonlinearity are captured at an equal footing. Importantly, before studying the KdV connection, we provide a rigorous analysis of the linear problem. We write down a set of theorems along with proofs and associated corollaries that shine light on the conditions of existence and nature of eigenvalues and eigenvectors of the linear problem. This rigorous analysis is paramount for understanding the nonlinear problem and the KdV connection. We provide strong evidence of agreement between VNLS systems and KdV equations by using soliton solutions as a platform for comparison. Our results are expected to be relevant not only for cold atomic gases, but also for nonlinear optics and other branches where VNLS equations play a defining role.

Equations

$$i\hslash \frac{\partial {\psi }_k}{\partial t}=-\frac{{\hslash }^2}{2m}\frac{{\partial }^2{\psi }_k}{\partial x^2}+\sum _{j=1}^N{\alpha }_{kj}|{\psi }_j{|}^2{\psi }_k$$

$$\mathcal{H}=\int dx\sum_{k=1}^{N}\bigg{(}\frac{\hbar^{2}|\partial_{x}\psi_{k}|^{2}}{2m}+\sum_{j=1}^{N}\frac{\alpha_{kj}}{2}|\psi_{k}|^{2}|\psi_{j}|^{2}\bigg{)}$$

$$\frac{\partial {\psi }_k}{\partial t}=\{{\psi }_k,\mathcal{H}\}\forall k$$

$${\psi }_k(x,t)=\sqrt{{\rho }_k(x,t)}{e}^{i(m/\hslash ){\int }_0^x{v}_k(x^{\prime },t)d{x}^{\prime }}$$

$$\frac{\partial\rho_{k}}{\partial t}+\frac{\partial}{\partial x}\big{(}\rho_{k}v_{k}\big{)}=0,$$

$$\frac{\partial v_{k}}{\partial t}=-\frac{\partial}{\partial x}\bigg{[}\frac{v_{k}^{2}}{2}+\frac{1}{m}\sum_{j=1}^{N}\alpha_{kj}\rho_{j}-\bigg{(}\frac{\hbar^{2}}{2m^{2}}\bigg{)}\frac{\partial_{x}^{2}\sqrt{\rho_{k}}}{\sqrt{\rho_{k}}}\bigg{]}.$$

$$\frac{\partial }{\partial t}\delta {\rho }_k=-{\rho }_{0k}\frac{\partial }{\partial x}\delta v_k,$$

$$\frac{\partial }{\partial t}\delta v_k=-\frac{1}{m}\frac{\partial }{\partial x}\sum _{j=1}^N{\alpha }_{kj}\delta {\rho }_j+\frac{{\hslash }^2}{4m^2{\rho }_{0k}}\frac{{\partial }^3}{\partial x^3}\delta {\rho }_k.$$

$${\rho }_k$$

$$v_k$$

$${\partial }_t\left(\begin{array}{c}\delta \rho \\ \delta v\end{array}\right)=-{\partial }_x\mathcal{A}\left(\begin{array}{c}\delta \rho \\ \delta v\end{array}\right),$$

$$\mathcal{A}=\left(\begin{array}[]{cc}\mathbf{0}_{N\times N}&\rho\\ \alpha&\mathbf{0}_{N\times N}\end{array}\right),$$

$$\alpha_{ij}=\left\{\begin{array}[]{ll}g_{i},&i=j,\\ h,&i\neq j,\end{array}\right.$$

$${\mathcal{S}}^k\left(\left\{{\gamma }_i\right\}\right)=\sum _{\sigma \in C_{Nk}}\prod _{{\gamma }_j\in \sigma }{\gamma }_j$$

$${\mathcal{S}}^0\left(\left\{{\gamma }_i\right\}\right)=1,{\mathcal{S}}^m\left(\left\{{\gamma }_i\right\}\right)=0,m>N$$

$${\mathcal{S}}^N\left(\left\{\frac{{\rho }_{0i}{g}_i-{\lambda }^2}{{\rho }_{0i}h}\right\}\right)$$

$$+$$

$${\left(\frac{{\rho }_0^{\ast }{g}^{\ast }-{\lambda }^2}{{\rho }_0^{\ast }h}-1\right)}^{m-1}{\left(\frac{{\rho }_0^{\prime }{g}^{\prime }-{\lambda }^2}{{\rho }_0^{\prime }h}-1\right)}^{m^{\prime }-1}\psi ({\lambda }^2)=0$$

$$\psi ({\lambda }^2)$$

$$\displaystyle\quad\quad\quad\>\left.(-1)^{N-m-m^{\prime}-p}f(\gamma,\gamma^{\prime},p)\bigg{.}\right]$$

$$f(\gamma ,{\gamma }^{\prime },p))$$

$$+(\gamma -1)({\gamma }^{\prime }-1)(1+N-m-m^{\prime }-p)$$

$$\gamma =\left(\frac{{\rho }_0^{\ast }{g}^{\ast }-{\lambda }^2}{{\rho }_0^{\ast }h}\right),{\gamma }^{\prime }=\left(\frac{{\rho }_0^{\prime }{g}^{\prime }-{\lambda }^2}{{\rho }_0^{\prime }h}\right).$$

$$\alpha \rho q^i={\rho }^{-1/2}{\rho }^{1/2}\alpha {\rho }^{1/2}{u}^i={\lambda }_i^2{\rho }^{-1/2}{u}^i={\lambda }_i^2{q}^i,$$

$$v^{i}_{\pm}=\left[\begin{array}[]{c}\dfrac{\pm 1}{\lambda_{i}}\rho\>q^{i}\vspace{5pt}\\ q^{i}\end{array}\right].$$

$$V=\left[\begin{array}[]{cc}\rho Q\Lambda^{-1}&-\rho Q\Lambda^{-1}\\ Q&Q\end{array}\right],$$

$$\det (V)$$

$$=\det (2\rho Q{\Lambda }^{-1}Q)\ne 0,$$

$$v^{i}_{\pm}=\left[\begin{array}[]{c}\pm\dfrac{1}{\lambda_{i}}\rho q^{i}\vspace{5pt}\\ q^{i}\end{array}\right],$$

$$v^{(k)}=\left[\begin{array}[]{c}\dfrac{\pm 1}{\sqrt{\rho^{*}_{0}(g^{*}_{0}-h)}}\rho q^{(k)}\vspace{5pt}\\ q^{(k)}\end{array}\right],$$

$$(q^{(k)})^{i}=\left\{\begin{array}[]{ll}1,&i=i_{1},\\ -1,&i=i_{k+1},\\ 0,&\mbox{else,}\end{array}\right.$$

$$\alpha \rho q={\lambda }^{2}q$$

$$\alpha \rho ={\alpha }_0\rho +h{\alpha }_1\rho$$

Full text

I Introduction

Multi-component coupled systems are ubiquitous in physics ranging from cold atomic systems Smerzi et al. (2003); Burchianti et al. (2018); Roati et al. (2007); Thalhammer et al. (2008); Ejnisman et al. (1998); Wacker et al. (2015); McCarron et al. (2011); Papp et al. (2008); Wang et al. (2015); Matthews et al. (1999) to nonlinear optics Chen et al. (1997); Ostrovskaya et al. (1999); Mitschke and Mollenauer (1987); Hasegawa (1980); Andrekson et al. (1991); Mitchell et al. (1998); Mitchell and Segev (1997). Such systems are typically nonlinear, i.e., with considerable interactions and often have an intricate interplay between the various species. Given these rich interactions, both intra-species and inter-species, the cutting-edge technologies Andrews et al. (1997, 1996); Anderson et al. (1995); Mewes et al. (1996) to image their collective behaviour and the ability to engineer these systems makes them a rich platform to study far-from-equilibrium physics in multi-component systems.

Often arriving at a Hamiltonian or a set of differential equations to describe the collective behaviour of particles is in itself a difficult task. However, a substantial work is done in this direction and there is a reasonable understanding of an effective Hamiltonian or differential equations that could describe multi-species systems in certain parameter regimes and conditions Agrawal (2000); Kevrekidis et al. (2008). However, the complex nature of these systems results in dealing with equations which are often cryptic and the consequences of which are difficult to understand. For example, if the collective behaviour of multi-species systems have nonlinearities and higher derivatives, one would expect to see nonlinear and dispersive effects. Such fingerprints of hydrodynamics are often completely elusive. Even the linearized version of the problem, existence of stable modes are unclear. Hence, it is of great importance to develop a systematic theory that will lead to a universal framework Kulkarni and Abanov (2012); Erdős et al. (2007); Dalfovo et al. (1999); Zakharov and Faddeev (1971); Ablowitz (2011); Pethick and Smith (2008); Horikis and Frantzeskakis (2014); Leblond (2008); Zakharov and Kuznetsov (1986); Spiegel (1980); Gardner and Morikawa (1960), that captures the various hallmarks of collective field theory or hydrodynamics.

We will describe here several systems where multi-component physics comes into play. Nonetheless, we will keep our main motivation as non-equilibrium dynamics in Bose mixtures Burchianti et al. (2018). In the context of cold atoms, VNLS (aka multi-component Gross$-$Pitaevskii equation) appears in multiple coupled species of bosons Liu et al. ; Zhou et al. (2008); Liu et al. (2018a); Kasamatsu and Tsubota (2005, 2006); Mareeswaran and Kanna (2016); Kasamatsu and Tsubota (2004); pat (2014); Feng (2014); Manikandan et al. (2016); Kuo and Shieh (2008); Caliari and Squassina (2008) or bosons with hyperfine degrees of freedom (spinor BECs Wen and Yan (2017); Sun and Wang (2018); Li and Yu (2017); Liu et al. (2018b); Belobo and Meier (2018); Massignan et al. (2015); Oztas (2019); Yakimenko et al. (2009); Hong-Qiang et al. (2011)). The nonlinearities or interactions one sees are of inter-species and intra-species type. Such systems can be placed out of equilibrium and their collective density dynamics can be imaged in-situ using cutting-edge technologies in absorption imaging techniques. On the other hand, non-linear optical setups also provide a great platform for studying multi-component NLS systems Afanasyev et al. (1989); Trillo et al. (1988). Here, typically, the role of time is played by an additional spatial axis Agrawal (2000). The intensity of light can be directly measured Mitchell et al. (1996); Martienssen and Spiller (1964) which is captured by a set of NLS equations. Apart from these two main avenues, NLS-type equations also appears in a variety of other contexts like Quantum Mechanics Rosales and Sánchez-Gómez (1992), accelerator dynamics Fedele et al. (1993), biomolecular dynamics Davydov et al. (1985); Daniel and Latha (2002); Qin et al. (2010), plasma and water waves Kourakis and Shukla (2006); Infeld and Rowlands (2000).

Keeping in mind the above motivation, we start with the description of $N$-component coupled NLS equations. Whenever appropriate, we will discuss the physical relevance mainly keeping collective description of cold atomic systems in mind. The $N$-component Nonlinear Schrödinger Equation (NLS) in 1D is given by

[TABLE]

where $\psi_{k}$ is the macroscopic wavefunction and $\alpha$ is the matrix of coupling constants. It is to be noted that the diagonal elements of the $\alpha$ matrix correspond to intra-species interaction and the off-diagonal elements correspond to inter-species interaction. We assume a symmetric coupling and hence $\alpha$ is a symmetric matrix. It turns out that in cold atomic systems, both intra-species and inter-species coupling are tunable via sweeping across a Feshbach resonance Inouye et al. (1998); Roati et al. (2007); Wacker et al. (2015); Thalhammer et al. (2008). The quantity $|\psi_{k}|^{2}$ gives the density of particles of species type $k$ and the angle associated with the complex number $\psi$ gives the phase both of which are measurable in experiments. The equations have a Hamiltonian structure given by the Hamiltonian

[TABLE]

equipped with Poisson brackets $\{\psi_{j}^{*}(x),\psi_{k}(y)\}=\frac{i}{\hbar}\delta_{jk}\delta(x-y)$. The NLS equations Ablowitz et al. (2004) of motion can be obtained from

[TABLE]

It is to be pointed out that single component version has been studied intensively in literature both from a point of mathematical interest Chiron (2012); Spiegel (1980); Huang et al. (2001); Kamchatnov and Pavloff (2012); Yan and Konotop (2009); Huang (2001) and experimentally in physical systems Anderson et al. (1995); Andrews et al. (1997); Davis et al. (1995); Andrews et al. (1996) and finds applications in a variety of fields Rosales and Sánchez-Gómez (1992); Fedele et al. (1993); Davydov et al. (1985); Mollenauer et al. (1980); Chang (2012); Kivshar (1990); Kuwamoto et al. (2004); Kevrekidis et al. (2007). Typical aspects studied both theoretically and experimentally include non-equilibrium evolution of density profiles, solitons Burger et al. (1999); Khaykovich et al. (2002); Kevrekidis and Frantzeskakis (2016), quenches Franchini et al. (2016), problems in presence of defects and disorder Paiva et al. (2015); Hulet et al. (2009); Chen et al. (2009, 2008). In addition to systems which have a Hamiltonian structure, there has been a lot of work on Driven-dissipative (gain-loss Smirnov et al. (2014); Kivshar and Luther-Davies (1998); Sich et al. (2012); Amo et al. (2011); Menon et al. (2010); Hivet et al. (2012); Lagoudakis et al. (2008); Amo et al. (2009)) and PT-symmetric systems El-Ganainy et al. (2007); Makris et al. (2008); Rüter et al. (2010). These are interesting open system or non-Hermitian generalizations of the NLS family of equations and are not a subject of our current work.

From (1) it is difficult to understand non-equilibrium phenomena. More precisely, given the initial condition $\psi_{i}(x,0)$, one is interested in the time evolution $\psi_{i}(x,t)$ from which various experimentally relevant quantities can be extracted. Experimentally, one can prepare an initial density profile $|\psi_{i}(x,0)|^{2}$ and associated phase (or its derivative which is akin to the velocity field). Then the system can be made to evolve and the time evolution of these quantities can be obtained. The central goal of our work is to provide a universal framework to understand the linear and non-linear properties of these evolving density and velocity fields of all species. A natural well-known conservative partial differential equation that is expected to capture fingerprints of hydrodynamics (namely, dispersion and nonlinearity) is the KdV equation. However, a systematic understanding of the linear problem and then the non-linear problem is far from obvious. Before going into the contents and details of our paper, below we provide a summary of our main findings.

The main contribution of the present manuscript is the systematic derivation of the qualitative long time dynamics of equation (1). In particular, we show that the evolution of small perturbations to the trivial state are governed by the KdV equation. For the single component case ($N=1$), the derivation of KdV from NLS is well known Gardner and Morikawa (1965); Su and Gardner (1969); Jeffrey and Kawahara (1982); Johnson (1997); Newell (1985). The present work provides the systematic and explicit derivation of KdV for the multi-component case. The coefficients for the KdV and the speed of sound (i.e. in the frame of reference for the evolution) are given in terms of the background trivial state, mass and the coupling coefficients $\alpha_{jk}$. We achieve this by a systematic and complete analysis of the spectral problem for the equation linearised about the trivial state. We derive necessary and sufficient conditions for real sound speeds dependent solely on the coupling matrix. For a specific case of the coupling matrix, we obtain necessary and sufficient conditions for sound speeds to be distinct. The coefficients of the KdV dynamics are given in terms of the eigenvectors of the linearised problem. For the case of repeated eigenvalues, we provide the eigenvectors explicitly. For the simple eigenvalue case, we construct the eigenvectors of the linearised problem in terms of the eigenvectors of the coupling matrix. We also provide an accurate efficient and stable numerical algorithm to compute the eigenvectors of the coupling matrix individually that also provides us with information on how eigenvalues (sound speeds) and eigenvectors (KdV coefficients) change as the cross-component coupling coefficient varies.

Though the bulk of the present work analyses the multi-component system (1) under Assumption 1, our analysis is general enough to be readily extended to other coupling matrices. We have attempted to provide the reader with sufficient details of the proofs to make this generalisation obvious. We will consider alternative couplings in future works. We have also attempted to highlight the role of a key mathematical idea, namely the Fredholm alternative, in reductive perturbative theory. To this end, we present perhaps more detail than is typical so that the interested reader can appreciate the systematic nature of the derivation.

The contents of the paper are organized as follows. In subsection 1 the coupled NLS system (1) is transformed into hydrodynamic form where the relevant physical quantities are the density and velocity fields. The linear regime of the thus obtained hydrodynamic equations is studied in subsection 2. Rescaling the independent variables, $x$ and $t$, so that we focus our attention on the dynamics on long space and time scales, we are lead to, at linear order, an equation that governs the linear stability of perturbations to a trivial background state. In subsection 3 we determine necessary and sufficient conditions for these perturbations to propagate stably in terms of the coupling matrix of the original system (1). Mathematically, this involves a complete analysis of the eigen-system for a particular matrix $\mathcal{A}$. This analysis includes an explicit formula for the characteristic polynomial; necessary and sufficient conditions for repeated roots; and expressions for the associated eigenvectors. The eigenvalues and eigenvectors of $\mathcal{A}$ play an important role in the qualitative dynamics of perturbations to the trivial state: (i) the eigenvalues determine the sound speed and (ii) eigenvectors determine coefficients of the effective KdV equation governing these perturbations. After introducing the relevant mathematical ideas in subsection 4, we proceed with the reductive perturbation method applied on the density and velocity fields to derive the effective KdV equation for $N$-components in § III. We present explicit results for case of few components ($N=2,3$) in § IV. In § V, we discuss the results of numerical comparison between the coupled ($N=2$) NLS and corresponding KdV by simulating a solitary wave profile which further explicates nontrivial features of the NLS in the reduced KdV. We finally conclude in § VI along with an outlook.

II Hydrodynamic Model and Linearization

1 Modeling

As mentioned in the introduction, one of the primary contributions of the present manuscript is a characterization of a multicomponent system, particularly in the small-amplitude long-wavelength regime. This regime naturally leads to a coupled multi-species KdV-like model. In order to derive the associated KdV model, we first perform the usual Madelung transform Madelung (1927) to obtain a set of hydrodynamic equations for the density and velocity,

[TABLE]

where $\psi_{k}$ is the macroscopic wavefunction of the $k^{th}$ condensate, $\rho_{k}(x,t)$ is the corresponding density field and $v_{k}(x,t)$, the velocity field where $k=1,\ldots N$. The resultant equations of motion are an equation of continuity (for the density)

[TABLE]

and the Euler equation (for the velocity)

[TABLE]

We remark that the density equation of $k-$th component is uncoupled; the coupling matrix only appearing in the velocity equations. It is worth emphasizing that the there is no approximation in obtaining the above equations. They are fully equivalent to the original multicomponent system.

2 Linearized dynamics

A trivial solution to the hydrodynamic equations is given by setting the densities to non-negative constants $\rho_{0k}$ and the velocities to zero. A natural question is how small perturbations to this background state evolve. The standard approach to this question is given by linearising the equations about the background state, namely setting $\rho_{k}=\rho_{0k}+\delta\rho_{k}$ and $v_{k}=\delta v_{k}$, to obtain the following linear evolution equations for the perturbations $\delta\rho_{k},\delta v_{k}$

[TABLE]

These constant coefficient equations are readily solved using Fourier methods. Furthermore, for the simpler case of $\alpha$ equal to a diagonal matrix (the uncoupled case), the above equations have solutions of the form $\mathrm{exp}(ipx+i\omega t)$ with $\omega\approx c_{1}p+c_{2}p^{3}$ for some real constants $c_{1},c_{2}$. This indicates waves travel at speed $c_{1}$ but also disperse due to the presence of the cubic term. This suggests that if we consider long-wavelength perturbations such that the dispersive term balances with the nonlinear corrections, we may arrive at a KdV-like model. This same argument applies also to the coupled case, i.e. for a generic matrix $\alpha$.

To obtain a balance between nonlinearity and dispersion, we assume the following form for the original full density and velocity:

[TABLE]

where $\epsilon$ is a small formal parameter. From here on, we limit ourselves to the case where $m$ is a scalar and common to all species. All of our analysis on the linear system (and the resultant perturbation scheme) extends to the case when each species has a corresponding distinct value for $m$, without any change in our conclusions. However, for the sake of simplicity of presentation we limit ourselves to a single common value for $m$. Moreover, we set that value to $m=1$ without loss of generality.

Substituting the form of the perturbation (8-9) into the hydrodynamic equations and dropping terms of $\mathcal{O}(\epsilon^{2})$ we obtain

[TABLE]

where

[TABLE]

with $\rho$ an $N\times N$ diagonal matrix with elements $\rho_{0k}>0$ and $\delta\rho\,,\delta v$ are the $N\times 1$ vectors for the perturbations in density and velocity. Real eigenvalues of the matrix $\mathcal{A}$ correspond to traveling wave solutions, either to right or left, depending on the sign of the eigenvalue. Complex eigenvalues, on the other hand, correspond to unstable exponentially growing modes. Evidently, for a stable background state we must choose the coupling matrix $\alpha$ such that all eigenvalues of $\mathcal{A}$ are real.

If all eigenvalues are indeed real, the solution consists of pulses propagating at speeds given by each of the eigenvalues. These speeds, referred to as sound speeds for the system, depend on the values of the coupling constants, background densities and mass of the condensate atoms. These may be readily measured in an experiment and compared to the theoretically predicted values. As we show in the subsequent section, for a specific model of the coupling matrix, the sound speeds are easily computed.

3 Spectral analysis of $\mathcal{A}$

As emphasized in the previous section, the eigenvalues of $\mathcal{A}$ play a crucial role in designing the multicomponent system. As we will see in the following section, both eigenvalues and eigenvectors play a key role in deriving equations governing the nonlinear dynamics of small amplitude pulses. The purpose of the present section is to obtain a full characterisation of the spectrum of $\mathcal{A}$. Indeed, we present arguments to compute eigenvalues, assure their reality, determine their multiplicity and calculate eigenvectors. As a consequence, we are able to determine the spectral decomposition of the $2N\times 2N$ matrix $\mathcal{A}$ as a function of the coupling coefficients $\alpha$. Our main assumption is that the coupling matrix is a real symmetric positive definite matrix. Under this assumption alone, Theorems 1 and 5 characterize the spectrum (eigenvalues and eigenvectors) of $\mathcal{A}$ in terms of the spectrum of $\alpha$ and the background densities. It is this information that is needed to perform the reductive perturbation theory of the following section.

To obtain more information on the eigenvalues/eigenvectors of $\mathcal{A}$ (mulitiplicity of eigenvalues, the characteristic polynomial, assure reality and positivity of eigenvalues of $\alpha$, etc.) we make the following assumption for $\alpha$

Assumption 1

[TABLE]

where $g_{i}\,,h$ are all positive constants.

Theorems 2, 3, 4 and 6 invoke the above assumption. On the other hand, the arguments presented in the proof of these theorems can be adapted to alternate forms of the coupling matrix. We present more detail than typical for the interested reader to adapt these arguments for alternate coupling matrices.

We begin with a sufficient condition to ensure real eigenvalues for $\mathcal{A}$. Note, that since $\mathcal{A}$ is not symmetric, we are not readily guaranteed real eigenvalues.

Theorem 1

If $\alpha$ is symmetric positive definite, then the eigenvalues of $\mathcal{A}$ are non-zero, real and come in pairs with opposite sign.

Proof. Any eigenvalue $\lambda$ of $\mathcal{A}$ satisfies the characteristic polynomial $\mathrm{det}(\mathcal{A}-\lambda)=\mathrm{det}(\lambda^{2}-\alpha\rho)$ where the equality follows from the block nature of $\mathcal{A}$. Hence, all $\lambda$ are real if and only if $\alpha\rho$ has positive eigenvalues. Since $\rho$ is a diagonal matrix with positive elements, we have $\alpha\rho=\rho^{-1/2}(\rho^{1/2}\alpha\rho^{1/2})\rho^{1/2}$. Hence $\alpha\rho$ is similar to $\rho^{1/2}\alpha\rho^{1/2}$ and therefore these matrices have the same eigenvalues. On the other hand, $\rho^{1/2}\alpha\rho^{1/2}$ is congruent to $\alpha$ and so these matrices have the same number of positive eigenvalues. Since $\alpha$ by assumption is positive definite, $\rho^{1/2}\alpha\rho^{1/2}$ has $n$ positive eigenvalues and consequently, $\alpha\rho$ has $N$ positive eigenvalues. Thus $\mathcal{A}$ has $N$ positive eigenvalues and $N$ negative eigenvalues.

We now present a necessary condition for the coupling matrix to be positive definite. Thus if the condition is violated, $\alpha$ cannot be positive definite and we expect perturbations to the background to be unstable. To be specific, we have the following theorem.

Theorem 2

Let $\alpha$ satisfy Assumption 1. Let the components of the system be ordered such that $g_{i}\leq g_{i+1}$, $i=1,2,\ldots$. Then $\alpha$ is a positive definite matrix implies $h<\sqrt{g_{1}g_{2}}$.

Proof. See Appendix A.

It is to be noted that the above is a necessary condition but not sufficient. A simple sufficient condition is that $h<g_{1}$. To see this, let $x$ be any vector in $\mathbb{R}^{N}$. A straightforward computation shows $x^{T}\alpha x=\sum_{i}x_{i}^{2}(g_{i}-h)+h\left(\sum_{i}x_{i}\right)^{2}$. If $h<g_{1}\Rightarrow h<\min g_{i}$, then $x^{T}\alpha x$ is automatically positive for any vector $x$ and so $\alpha$ is in fact positive definite.

Corollary 2.1

The theorem above assumed a particular ordering. However, the eigenvalues of the system do not depend on the ordering; different orderings being obtained as mere permutations of the same system of equations. Hence the conclusion of the theorem holds when $g_{1},g_{2}$ are interpreted as the two smallest diagonal elements of $\alpha$.

Assumption 1 allows us to compute the characteristic polynomial of $\mathcal{A}$ in closed form. We first, however, introduce some notation.

Definition 1

Let $\{\gamma_{i}\}_{i=1}^{N}$ be a list of real numbers and $C_{Nk}$ be the set of all possible ways to choose any $k$ of these $N$ real numbers. We denote by $\mathcal{S}^{k}\left(\left\{{\gamma_{i}}\right\}\right)$ the symmetric product

[TABLE]

The definition above gives symmetric products of $\gamma_{i}$. For example, $\mathcal{S}^{n}\left(\left\{{\gamma_{i}}\right\}\right)=\prod_{j}\gamma_{j}$ and $\mathcal{S}^{1}\left(\left\{{\gamma_{i}}\right\}\right)=\sum_{j}\gamma_{j}$. In the case of only three elements $\gamma_{i},\ i=1,2,3$, then $\mathcal{S}^{2}\left(\left\{{\gamma_{i}}\right\}\right)=\gamma_{1}\gamma_{2}+\gamma_{2}\gamma_{3}+\gamma_{3}\gamma_{1}$.

Definition 2

We define

[TABLE]

Using the above notation, we state the following theorem.

Theorem 3

The characteristic polynomial of $\mathcal{A}$ is

[TABLE]

when $\alpha$ satisfies Assumption 1.

The above theorem may be proved using the principle of induction; the base case $N=2$ is easily checked by hand. Though straightforward, the proof by induction is involved and employs some particular properties of the symmetric products we have defined. The details are presented in Appendix B.

The characteristic polynomial of $\mathcal{A}$ is evidently a polynomial in $\lambda^{2}$ (as expected from Theorem 1). However, it is also a polynomial in $h$ as seen by multiplying the entire expression (14) by $h^{N}$. Indeed, the characteristic polynomial is a polynomial in $\lambda^{2}$ with coefficients that are polynomials in $h$. Such an expression is called an algebraic curve. These objects are well studied amongst mathematicians, though we only require some very basic properties of such polynomials.

Note for each $h$, there are $2N$ values of $\lambda$ that are roots of the characteristic polynomial (choosing $h$ suitably so that the values of $\lambda$ are all real). Since the eigenvalues of $\mathcal{A}$ represent the physical speeds of the small amplitude pulse-like initial perturbations to the base state, we are especially interested in knowing whether all sound speeds of the system are distinct. Distinct speeds correspond to pulses propagating such that eventually they do not interact. Consequently, determining whether the sound speeds are distinct is equivalent to determining whether the eigenvalues of $\mathcal{A}$ are distinct.

An important property of algebraic curves is that either (a) for most values of $h$ (in particular, except for a finite number of isolated values of $h$), the roots $\lambda^{2}$ of (14) are distinct, or (b) the number of distinct roots of (14) is less than $2N$ for all values of $h$.

A polynomial under case (b) is said to be permanently degenerate. We will state necessary and sufficient conditions on the coupling matrix $\alpha$ such that case (b) holds. Moreover, we will also determine an analytic expression for the repeated eigenvalues. A multicomponent system prepared such that case (a) is true, namely, when (b) does not hold will almost surely have distinct sound speeds since to have repeated eigenvalues, very precise values of $h$ must be chosen. Any perturbation of these particular values for $h$ will immediately lead to distinct sound speeds. If the conditions for case (b) are not satisfied, then we will assume the eigenvalues are simple for a given $h$. In other words, the sound speeds are distinct. The following theorem states conditions for case (b).

Theorem 4

Suppose $\alpha$ satisfies Assumption 1. The characteristic polynomial of $\mathcal{A}$ is permanently degenerate with a root of multiplicity $m-1$ if and only if $m$ pairs of $(\rho_{0i}g_{i},\rho_{0i})$ are equal. Furthermore, if the common value of the $(\rho_{0i}g_{i},g_{i})$ pairs is denoted by $(\rho_{0}^{*}g^{*},\rho_{0}^{*})$, then $\pm\sqrt{\rho_{0}^{*}(g^{*}-h)}$ are the associated repeated eigenvalue of $\mathcal{A}$.

Proof. See Appendix C.

Corollary 4.1

Some implications of the above theorem are,

• •

if two $(\rho_{0i}g_{i},\rho_{0i})$ pairs are equal (to say $\rho_{0}^{*}g^{*},\rho_{0}^{*}$), then $\pm\sqrt{\rho_{0}^{*}(g^{*}-h)}$ must be eigenvalues of $\mathcal{A}$.

• •

the above theorem is true for any repeated pair $(\rho_{0i}g_{i},\rho_{0i})$ and thus holds for each repeated pair.

• •

suppose another $m^{\prime}$ pairs of $(\rho_{0i}g_{i},\rho_{0i})$ were equal (and distinct from the first $m$ pairs). Then the characteristic polynomial takes the following form

[TABLE]

with

[TABLE]

where

[TABLE]

and

[TABLE]

Having fully characterized the eigenvalues of $\mathcal{A}$, we now proceed to investigate the eigenvectors. Our first result states that $\mathcal{A}$ is diagonalisable. Note this result does not require Assumption 1.

Theorem 5

Assuming the matrix $\alpha$ is positive definite, the matrix $\mathcal{A}$ has $2N$ independent eigenvectors. In particular, the algebraic and geometric multiplicities of any eigenvalue of $\mathcal{A}$ are equal for all permissible values of $h$. In other words, the matrix $\mathcal{A}$ is diagonalisable.

Proof. The $N\times N$ matrix $\rho^{1/2}\alpha\rho^{1/2}$ is real, symmetric and positive definite whenever $\alpha$ is real, symmetric and positive definite. Hence there exist $N$ mutually orthogonal eigenvectors $u^{i},\ i=1,2,\ldots N$ that span $\mathbb{R}^{N}$. Denote the eigenvalue associated with $u^{i}$ by $\lambda^{2}_{i}$ and let $q^{i}=\rho^{-1/2}u^{i}$. Then we have

[TABLE]

and hence $q^{i}$ is an eigenvector of $\alpha\rho$ with eigenvalue $\lambda^{2}_{i}$. This is true for each $i$ and hence we have determined $N$ eigenvectors for $\alpha\rho$. Consequently, the matrix $\alpha\rho$ is diagonalisable: the algebraic and geometric multiplicities are equal. We claim each $q^{i}$ induces two $2N-$dimensional eigenvectors $v^{i}_{\pm}$ for the matrix $\mathcal{A}$ corresponding to eigenvalues $\pm\lambda_{i}$. Indeed, define

[TABLE]

It then follows that $\mathcal{A}v^{i}_{\pm}=\pm\lambda_{i}v^{i}_{\pm}$ where $\lambda_{i}$ is the positive root of $\lambda^{2}_{i}$, the associated eigenvalue for $\alpha\rho$. Let $Q$ represent the matrix with columns $q^{i}$ and $\Lambda$ be the diagonal matrix with diagonal elements $\lambda_{i}$. Then the matrix whose columns are eigenvectors of $\mathcal{A}$ is given by

[TABLE]

with determinant

[TABLE]

since $\rho,Q,\Lambda^{-1}$ are all invertible matrices. Thus the columns of $V$ are linearly independent and hence there are $2N$ independent eigenvectors for $\mathcal{A}$.

Corollary 5.1

If $q^{i}$ is an eigenvector of $\alpha\rho$ with eigenvalue $\lambda^{2}_{i}$, then

[TABLE]

is an eigenvector of $\mathcal{A}$ with eigenvalue $\pm\lambda_{i}$. Hence computing eigenvectors of $\mathcal{A}$ is equivalent to computing those of $\alpha\rho$.

When $\alpha$ satisfies Assumption 1 and the characteristic polynomial is permanently degenerate (has repeated eigenvalues for all values of $h$), Theorem 4 states the exact expression of these repeated eigenvalues and their multiplicity. The following theorem provides an exact form for the associated eigenvector in this case.

Theorem 6

Suppose $\alpha$ satisfies Assumption 1 and the matrix $\mathcal{A}$ is permanently degenerate, i.e. there is an eigenvalue of multiplicity $m>1$ for all suitable values of $h$. Let $(\rho^{*}_{0}g^{*},\rho^{*}_{0})$ represent the repeated pair resulting in the degenerate eigenvalue. A set of $m$ independent eigenvectors for these permanently degenerate eigenvalues $\lambda=\pm\sqrt{\rho^{*}_{0}(g^{*}_{0}-h)}$ of multiplicity $m$ are given by

[TABLE]

where the $i-$th component of the $k-$th eigenvector is given by

[TABLE]

for $k=1,2,\ldots m$. The indices $i_{k},\>k=1,2,\ldots,m+1$ are such that the diagonal elements of $\alpha\rho$ at these locations $(\alpha\rho)_{i_{k}i_{k}}=\rho^{*}_{0}g^{*}$.

Proof. See Appendix  D.

3.a Numerical method to compute eigenvalues and eigenvectors of $\mathcal{A}$

The previous theorems establish the reality of the eigenvalues of $\mathcal{A}$. Moreover, we have explicit formulae for the permanently repeated eigenvalues (if any) and their associated eigenvectors. It remains to investigate the eigenvalues and eigenvectors which are typically simple, i.e. simple for most values of $h$. Although one could simply use a standard numerical solver to compute roots of the characteristic polynomial for various $h$, we present some further analytic results and a simple iterative procedure that determines both eigenvalues and eigenvectors.

From Theorem 4, non-degenerate eigenvalues correspond to those $(\rho_{0i}g_{i},\rho_{0i})$ pairs which do not repeat. Notice that when $h=0$, the characteristic polynomial (3) has roots $\lambda=\pm\sqrt{g_{i}\rho_{0i}}$. We conclude that at $h=0$, the non-repeating $\rho_{0i}g_{i}$ are simple eigenvalues. It is known that this behaviour, namely the simple nature of the eigenvalue, must persist at least for small $h$ Kato (2013). The claim essentially follows from the implicit function theorem. We now present a scheme to compute both eigenvalues and eigenvectors as a function of $h$ that limit to a non-repeated $\rho_{0i},g_{i}$ when $h=0$. This presentation employs the particular structure of the coupling matrix dictated by Assumption 1.

We work directly with the matrix $\alpha\rho$ instead of $\mathcal{A}$. Theorems 1 and 5 readily allow us to translate spectra between the two matrices. An eigenvalue-eigenvector pair for $\alpha\rho$ satisfies the following equation

[TABLE]

where $q$ is the eigenvector. When $h=0$, eigenvalues and eigenvectors are readily available: $\rho_{0i}g_{i}$ are eigenvalues with the canonical basis in $\mathbb{R}^{N}$ as the associated eigenvectors. We proceed to compute eigenvalues and eigenvectors for non-zero $h$ as follows.

Let $q^{(k)}=e_{k}+hq_{h}^{(k)}$ be the eigenvector for $\alpha\rho$ associated with eigenvalue $\lambda^{2}=\rho_{0k}g_{k}+h\mu_{h}^{(k)}$ for $h\neq 0$. Here we consider only those cases when $\rho_{0k},g_{k}$ is not repeated and $e_{k}$ represents the $k-$th canonical unit vector in $\mathbb{R}^{N}$. Under Assumption 1 the matrix $\alpha\rho$ may be written as

[TABLE]

where

[TABLE]

Substituting for $\lambda^{2}$, $q$ and $\alpha\rho$, we have after some rearrangement

[TABLE]

Note that $\alpha_{0}\rho-\rho_{0k}g_{k}$ is diagonal and has a null-space: $e_{k}$. Thus we stipulate $q_{h}^{(k)}$ be orthogonal to $e_{k}$ and further require that the right-hand side of the above equation also be orthogonal to $e_{k}$. In more mathematical parlance, we are invoking the Fredholm alternative. The orthogonality condition for the right-hand side leads to

[TABLE]

In other words, the correction to the eigenvalue is given by $h$ times the $k-$th component of $\alpha_{1}\rho\>q_{h}^{(k)}$. Substituting this expression into the equation for $q_{h}^{(k)}$ we obtain

[TABLE]

This equation is iteratively solved for $q_{h}^{(k)}$ since the left-hand side matrix is invertible (when $q_{h}^{(k)}$ is orthogonal to $e_{k}$). Notice that $\alpha_{0}\rho-\rho_{0k}g_{k}$ is diagonal and hence readily inverted. Once we converge to a $q_{h}^{(k)}$ that satisfies the above equation for some $h>0$, we evaluate the eigenvalue as

[TABLE]

We notice that the correction to the eigenvalue for $h>0$ is quadratic in $h$. If one expands $q_{h}^{(k)}$ in a power series of $h$, we find

[TABLE]

where the prime indicates the $k-$th term is skipped.

The iteration procedure described above may be justified by appealing to the implicit function theorem. Here we consider the correction to the eigenvector $q^{(k)}_{h}$ as a function of $h$. The requirements of the implicit function theorem hold at the point $h=0,q^{(k)}_{h}=0$, i.e. the linearisation of expression (32) at $h=0,q^{(k)}_{h}=0$ leads to an invertible matrix.

At the outset, we do not know how large the radius of convergence (in $h$) of the resultant series is. However, given an eigenvalue-eigenvector pair for $\alpha\rho$ for $h\neq 0$, say $(e_{k}+hq^{(k)}_{h},g_{k}\rho_{k}+h\mu_{h}^{(k)})$, we may repeat the perturbation argument and restart the series around a non-zero value of $h$. Hence setting $q^{(k)}=e_{k}+hq^{(k)}_{k}+\delta hp^{(k)}_{\delta h}$ and $\lambda^{2}=g_{k}\rho_{0k}+h\mu^{(k)}_{h}+\delta h\nu^{(k)}_{\delta h}$ we obtain the following equation for $p^{(k)}_{\delta h}$

[TABLE]

where

[TABLE]

Once again by appealing to the implicit function theorem, one can establish that equation (35) can be solved for $p^{(k)}_{\delta h}$ for sufficiently small $\delta h$. By repeatedly using the above argument, we may obtain the eigenvector-eigenvalue for $\alpha\rho$ for all suitable $h$.

4 Inhomogeneous linear dynamics

With an eye towards the calculations in the next section, we now discuss the solution procedure for inhomogeneous equations of the form

[TABLE]

where we assume $f$ is a known $N\times 1$ vector valued function and we wish to determine the $N\times 1$ vector $s$. The main tool we employ is the Fredholm alternative.

The Fredholm alternative is a statement on the solvability of linear equations. Consider a matrix equation $Lx=b$, where $L$ is a square matrix and $b$ is known. If $L$ is invertible, the solution is readily available: $x=L^{-1}b$. If however $L$ is not invertible, a necessary condition for a solution is that $b$ must be orthogonal to all $y$ such that $L^{T}y=0$. This is the alternative. Note if $\langle y,b\rangle=0$ for all such $y$, we may have an infinite number of solutions. A unique solution can be obtained from an infinite possible set, if we also suppose $x$ is orthogonal to the null space of $L$.

The above considerations for a matrix apply also to differential operators. Consider the equation $\partial_{X}\psi(X)=\zeta(X,T)$, where $\zeta$ is known. Clearly any constant is in the null space of the operator $\partial_{X}$. The adjoint of $\partial_{X}$ is $-\partial_{X}$ which also has constants as its null space. Hence we require the function $\zeta$ to be orthogonal to constants. To make these statements rigorous we need to state appropriate Hilbert spaces and inner products. We will avoid such technicalities presently.

Coming back to equation (39), we recall that $\mathcal{A}$ has a spectral decomposition

[TABLE]

where $\Lambda$ is a diagonal matrix with positive elements such that $\alpha\rho Q=Q\Lambda^{2}$. Substituting this into the inhomogeneous equation (39) we obtain

[TABLE]

The null space of the differential operator on the left-hand side above consists of vectors of the form $\psi(x-\tilde{\lambda}_{j}t)e_{j}$, where $e_{j}$ is the $j-$th canonical unit vector in $\mathbb{R}^{2N}$, $\tilde{\lambda}_{j}$ is the $j-$th element along the diagonal of $\tilde{\Lambda}$ and $\psi$ is any function. The null space of the adjoint is the same. Thus the condition to solve the above equation for $V^{-1}s$ is that $\langle e_{j},V^{-1}f\rangle$ should not be a function of $(X-\tilde{\lambda}_{j}T)$. In other words $\langle e_{j},V^{-1}f\rangle$ should not be proportional to $\phi(X-\tilde{\lambda}_{j}T),$ for any function $\phi$. In the next section we will see how this analysis of the linear inhomogeneous equation serves us in deriving equations governing the slow evolution of perturbations to the multi-component system (1).

III Reductive Perturbation method for N-component NLS

In the previous section we analysed the linearised equations for the perturbations $\delta\rho_{k},\delta v_{k}$ about the trivial state (7). The fully nonlinear equations for the perturbations (in hydrodynamic form) without any additional scaling are

[TABLE]

where

[TABLE]

Let us rescale the variables so that $\partial_{t}\to\epsilon\partial_{t}$, $\partial_{x}\to\epsilon\partial_{x}$ and $\delta\rho,\delta v\to\epsilon^{2}\delta\rho,\epsilon^{2}\delta v$. This scaling is equivalent to assuming the following for the original physical variables

[TABLE]

This leads to

[TABLE]

where $X=\epsilon x,T=\epsilon t$ and

[TABLE]

and $\delta\rho_{k}^{\prime\prime}=\partial_{X}^{2}\delta\rho_{k}$.

We now solve the above equation perturbatively. Assuming an expansion in $\epsilon^{2}$ for the unknowns

[TABLE]

and substituting into the above equation we obtain to lowest order

[TABLE]

Since $\mathcal{A}=V\tilde{\Lambda}V^{-1}$ is diagonalisable, this equation is equivalent to

[TABLE]

which has a solution

[TABLE]

or

[TABLE]

where $f^{(0)}_{j}(\xi)$ is any function, $e_{j}$ is the unit vector in $\mathbb{R}^{2N}$ and $\tilde{\lambda}_{j}$ is any of the eigenvalues of $\mathcal{A}$. As common in the method of multiple scales, we will assume $f^{(0)}_{j}$ depends on $X-\tilde{\lambda}T$ as well as a new slow time scale $\tau=\epsilon^{2}T=\epsilon^{3}t$. Hence

[TABLE]

The equations at order $\epsilon^{2}$ are then given by

[TABLE]

which is equivalent to

[TABLE]

Notice this is a linear inhomogeneous equation for the order $\epsilon^{2}$ correction to $\delta\rho,\delta v$. The right-hand side is essentially a known function since every term on the right hand side can be written in terms of $f^{(0)}_{j}(X-\tilde{\lambda}_{j}T,\tau)$. Moreover, the adjoint of the linear operator on the left-hand side has a null space: precisely those functions of the form $\psi(X-\tilde{\lambda}_{j}T)e_{j}$. From the Fredholm alternative, the right-hand side should be orthogonal to this null space. Notice all terms on the right are of the form $\psi(X-\tilde{\lambda}_{j}T)$. Hence we have the solvability condition

[TABLE]

where $e_{j}$ being the $j-$th unit vector of the $2N\times 2N$ identity matrix. The above equation can be simplified using the expression for the zeroth order solution to

[TABLE]

or in other words

[TABLE]

Rewriting $\mathcal{N}_{1},\mathcal{N}_{2}$ entirely in terms of $f^{(0)}_{j}$ we have the required KdV equation. This is true for each $j$ and hence we have $2N$ KdV equations, $N$ of which correspond to perturbations traveling to the right and $N$ of which correspond to perturbations traveling to the left. If all sound speeds (i.e. eigenvalues $\lambda_{j}$ of $\mathcal{A}$) are distinct, these equations are uncoupled since the terms $\mathcal{N}_{1},\mathcal{N}_{2}$ are given entirely in terms of the profile $f^{(0)}_{j}$. Physically this corresponds to moving into different traveling frames centered around each pulse.

To summarize, if $V$ is the eigenvector matrix of $\mathcal{A}$, namely $\mathcal{A}=V\tilde{\Lambda}V^{-1}$ then the density and velocity vectors are given by

[TABLE]

where $\lambda_{j}$ is eigenvalue associated with the eigenvector $Ve_{j}$ ($e_{j}$ being the $j-$th unit vector of the $2N\times 2N$ identity matrix).

The upshot of the above analysis is a reduction of the dynamics of the coupled GPE (1) in terms of $2N$ KdV equations. Physically speaking, a generic perturbation to the background densities, on the shortest timescale evolves according to (1) in such a way so as to give rise to $2N$ small-amplitude waves traveling at the sound speeds (given by the eigenvalues of $\mathcal{A}$). On a longer timescale, the waves evolve according to the corresponding KdV equation given in (63). Hence a generic perturbation to the background state resolves into $2N$ waves evolving according to KdV ($N$ going to the right; $N$ going to the left).

1 Coupled KdV equations (non-distinct speeds)

In the derivation presented in the previous section we assumed a solution to the homogeneous problem that depended only on one profile $f^{(0)}_{j}$. This is however not the most general solution. Indeed one may have well assumed

[TABLE]

Evidently now

[TABLE]

contains functions of all $X-\lambda_{j}T$ leading to what one may consider to be a coupled system of KdV. However, it must be noted, that when we project onto $e_{j}$ to obtain the equation of evolution for $f^{(0)}_{j}(X-\lambda_{j}T)$ we only retain those terms for the right-hand side which are functions of $X-\lambda_{j}T$ alone (and not products of functions of multiple $X-\lambda_{j}T)$. Hence once again we end up with uncoupled equations.

Note the above argument fails when $\mathcal{A}$ has repeated eigenvalues. This is precisely why we determined necessary and sufficient conditions for simple non-repeating eigenvalues. In the case of repeated eigenvalues however, the zeroth order solution is given as

[TABLE]

where now the sum only extends over those vectors $Ve_{k}$ which correspond to the same eigenvalue $\lambda$. All the functions $f_{j}^{(0)}$ depend on the same spatial variable $\xi=X-\lambda T$ and slow time scale $\tau=\epsilon^{2}T$. In the case of repeated eigenvalues, we necessarily obtain a coupled system of KdV equations; the number of equations is equal to the multiplicity of the eigenvalue. The dynamics of the coupled KdV equations arising out of repeated eigenvalues for $\mathcal{A}$ will be discussed in a future paper.

A final scenario that may also lead to coupled equations is when the eigenvalues are close together, indeed when $|\lambda_{j+1}-\lambda_{j}|<\epsilon^{2}$. In such a case, although the asymptotic behaviour of such a system is described by two uncoupled KdV equations (since they correspond to two different traveling frames of reference), due to the small difference in sound speeds, the dynamics may appear to be coupled even on the longer time scale for KdV-type equations. Note however, the resultant coupled system will typically have different coefficients than the one corresponding to repeated eigenvalues (when the sound speeds are exactly the same) since the associated eigenvectors are different in either case.

2 Some useful relations

Suppose $Q$ is the matrix of eigenvectors of $\alpha\rho$ with eigenvalues given by the diagonal matrix $\Lambda^{2}$. In other words $\alpha\rho Q=Q\Lambda^{2}$. All diagonal entries of $\Lambda^{2}$ are positive (see proof of Theorem 5). Then the matrix of eigenvectors for $\mathcal{A}$ is given by

[TABLE]

with inverse

[TABLE]

For the purposes of deriving KdV, the relevant matrix is $V$ and $(V^{-1})^{T}$. It turns out, one may express $(V^{-1})^{T}$ explicitly in terms of $\rho,Q,\Lambda$. Indeed one has

[TABLE]

where $L=Q^{T}\rho Q$. Moreover $L$ is a diagonal matrix with positive elements and hence $L^{-1}$ is readily computed. Note that $L=Q^{T}\rho Q$ implies $I=L^{-1}Q^{T}\rho Q$ and hence $Q^{-1}=L^{-1}Q^{T}\rho$.

The above statements expressing $Q^{-1}$ in terms of $Q$ and $\rho$ are explained as follows. The matrix $\rho^{1/2}\alpha\rho^{1/2}$ is real symmetric and positive definite (assuming $\alpha$ is symmetric positive definite). As a result there exists a orthonormal matrix $U$ which is the eigenvector matrix of $\rho^{1/2}\alpha\rho^{1/2}$. Then

[TABLE]

But then we have $\rho^{-1/2}U$ is also an eigenvector matrix of $\alpha\rho$ and $\Lambda^{2}=\mu^{2}$. This means $\rho^{-1/2}UM=Q$ where $M$ is a real diagonal matrix. In other words, the columns of $Q$ are parallel to columns of $\rho^{-1/2}U$. Expressing $U$ in terms of $Q,M,\rho^{1/2}$ and substituting in $UU^{T}=I$, $I$ being the $N\times N$ identity matrix, leads to $M^{2}=Q^{T}\rho Q$. We then define $L=M^{2}$, which is the matrix that appears in (76).

We also note that if $D_{1},D_{2}$ are diagonal matrices, then elements of $D_{1}\alpha D_{2}$ are given by $(D_{1}\alpha D_{2})_{ij}=(D_{1})_{i}\alpha_{ij}(D_{2})_{j}$ where $(D_{1})_{i},(D_{2})_{j}$ denote the $i-$th and $j-$th diagonal entry of $D_{1},D_{2}$ respectively. Hence any column of $V$ and $(V^{-1})^{T}$ is readily computed once the relevant column of $Q$ is determined. We recall Theorem 6 and the procedure detailed in subsection subsubsection 3.a allow us to compute a column of $Q$ independently of other columns. Of course, standard libraries provide all eigenvalues and eigenvectors simultaneously.

With these definitions the coefficients of the respective KdV equations are obtained in a straightforward manner by considering the relevant column of $(V^{-1})^{T}$ for

[TABLE]

is equivalent to

[TABLE]

where

[TABLE]

We re-emphasize there are in total $2N$ KdV equations in (79).

3 A special case: KdV with zero nonlinearity

The matrix $\alpha$ represents the coupling between the different species and the matrix $\rho$ represents the trivial background states for the different species. Let us consider a case when two self-couplings (diagonal elements of $\alpha$) and their corresponding background states (the respective diagonal elements of $\rho$) are equal. In other words we assume $g_{i_{1}}=g_{i_{2}}=g^{*}$ and $\rho_{0i_{1}}=\rho_{0i_{2}}=\rho^{*}_{0}$ for some indices $i_{1},i_{2}$. From Theorem 4 we are guaranteed that $\pm\sqrt{\rho^{*}_{0}(g^{*}-h)}$ are eigenvalues of $\mathcal{A}$. If more than two self-coupling–density pairs are equal, then the eigenvalues will have a multiplicity greater than one. Higher order multiplicities will be the focus of a future work and here we consider only the case of simple eigenvalues $\pm\sqrt{\rho^{*}_{0}(g^{*}-h)}$. We also limit the present discussion to the eigenvalue corresponding to waves traveling to the right. The analysis in this section extends similarly to the one traveling to the left.

Using Theorem 6 we also know the exact form of the eigenvector associated with the eigenvalue $\sqrt{\rho^{*}_{0}(g^{*}-h)}$. Indeed it is

[TABLE]

where $\rho$ is the matrix with diagonal elements $\rho_{0i}$ and the $i-$th component of $q$ is given by

[TABLE]

To obtain the relevant coefficients of KdV for this case, we need the relevant column of $(V^{-1})^{T}$ which is

[TABLE]

where $l$ is the element of $L=Q^{T}\rho Q$ corresponding to the vector $q$. Since $L$ is diagonal and $L=(\rho^{1/2}Q)^{T}\rho^{1/2}Q$ then $l=\sum_{i}\rho_{0i}q^{2}_{i}=2\rho^{*}_{0}$. Thus

[TABLE]

Then the KdV equation is given by

[TABLE]

which upon simplifying is

[TABLE]

where $f^{\prime\prime\prime}$ denotes the third-derivative with respect to the profile variable $X-\sqrt{(\rho^{*}_{0}(g^{*}-h)}\,T$. Remarkably this equation is linear and thus readily solved using Fourier transform techniques.

IV Few component GPE/NLS and associated KdV and other findings

In this section we present explicit results for two specific cases, that of two and three component coupled systems. The reality of the eigenvalues (positivity of elements of $\Lambda^{2}$) ensures that all coefficients in the final KdV equation are well-defined real numbers.

1 $N=2$ case

The coupled NLS are as follows

[TABLE]

In other words, the $\alpha$ matrix (under the assumption $m=1$; see text below Assumption 1) is given by

[TABLE]

The above equations, (87) can be written in a hydrodynamic form for a perturbation of a trivial state as given in Eqs. (47) to (49) where $k={1,2}$. We perform the perturbation series (eqs. 50 and 51) to arrive at Eq. (55) which is the equation for the lowest order where the $\mathcal{A}$ is given by,

[TABLE]

where the subscript in $\mathcal{A}_{2}$ means that we are dealing with the two component case. The eigenvalues of the matrix $\alpha\rho$ are the diagonal elements of

[TABLE]

and eigenvectors are respectively, the columns of,

[TABLE]

where, $A,B,C$ are given by,

[TABLE]

Note that the eigenvalues (diagonal elements of $\Lambda^{2}$) are both positive under the assumption $h<\sqrt{g_{1}g_{2}}$. This condition on $h$ is indeed not only necessary (Theorem 2) but also sufficient for the two component case. Next we employ the relations of Section 2, using the definitions of $Q$ and $\Lambda$ given above, to determine the coefficients of the respective KdV equations. Specifically, we determine the eigenvector matrix of $\mathcal{A}_{2}$ using (70). Similarly, we may compute $(V^{-1})^{T}$ using (76). In the following, we only present the results for the right chiral sector i.e. two positive eigenvalues (sound speeds). Eq. 79 for the two component species case (i.e., $N=2$) explicitly reads as follows for the two eigenvalues $\lambda=\sqrt{\dfrac{A\pm B}{2}}$ respectively,

[TABLE]

where,

[TABLE]

Here, $A_{1},B_{1}$ are the KdV coefficients for the largest positive eigenvalue, $\lambda_{1}=\sqrt{\frac{A+B}{2}}$ and $A_{2},B_{2}$ are the KdV coefficients for the second largest positive eigenvalue, $\lambda_{2}=\sqrt{\frac{A-B}{2}}$. Here $f^{\prime\prime\prime}_{j}$ denotes the third-derivative with respect to the profile variable $X-\lambda_{j}T$.

Equations (93) give the dynamics of perturbations to the background state in reference frames moving to the right with with speeds $\sqrt{(A\pm B)/2}$. Needless to say, there are two KdV equations for the other chiral sector namely for perturbations moving to the left with speeds $-\sqrt{(A\pm B)/2}$. The KdV equations for perturbations traveling to the left are obtained by setting $A_{j}\to-A_{j}$ in (93). Hence there are in total four KdV equations. § V contains the brute force numerical comparison between above KdV equation and the $N=2$ NLS case.

2 $N=3$ case

The $N=3$ case poses an interesting scenario. In general, for arbitrary $g_{1},g_{2},g_{3}$ and $\rho_{01},\rho_{02},\rho_{03}$ (all different) the eigenvalues are very cumbersome but with our prescription outlined in the previous sections one can explicitly write it down. Below, we describe the situation when $g_{1}=g_{2}$ and $\rho_{01}=\rho_{02}$. Here the eigenvalues and eigenvectors are still different. However, one of the eigenvalues and its corresponding eigenvector takes a particularly simple form.

The three-coupled NLS are as follows

[TABLE]

In other words, the $\alpha$ matrix (with $m=1$) is given by

[TABLE]

Correspondingly, we get

[TABLE]

The eigenvalues of the matrix $\alpha\rho$ are the diagonal elements of

[TABLE]

where we define,

[TABLE]

Corresponding to each of the three positive eigenvalues $\lambda_{j}$ of $\mathcal{A}$ (diagonal elements of $\Lambda$ in (101)), we obtain a KdV equation of the form

[TABLE]

Remarkably the corresponding KdV equation for the first eigenvalue $\lambda_{1}=\sqrt{(g_{1}-h)\rho_{01}}$ is linear. Indeed

[TABLE]

For the other two eigenvalues, i.e., $\lambda_{2}=\sqrt{\frac{X-Y}{2}}$ and $\lambda_{3}=\sqrt{\frac{X+Y}{2}}$, we get $A_{2},B_{2},A_{3},B_{3}$ as follows,

[TABLE]

where $Z\,,W$ are given by,

[TABLE]

In other words, the two eigenvalues other than $\lambda_{1}=\sqrt{(g_{1}-h)\rho_{01}}$ have corresponding KdV equations for the nonlinear problem.

We emphasize once more that equations (104) give the dynamics of perturbations that travel in the positive direction with speeds $\lambda_{j}$. There are three KdV equations for the other chiral sector namely, perturbations traveling to the left with speeds $-\lambda_{j}$. These KdV equations are obtained by setting $A_{j}\to-A_{j}$ in (104). Hence there are in total six KdV equations.

V Numerical Results for two component case

Following the results obtained in Sec. IV for $N=2$, in this section, we explore through numerical simulations the comparison between the obtained KdV and the coupled NLS. From our asymptotic analysis in that section, recall that a generic perturbation to the constant background density resolves into $2N$ weakly-nonlinear waves that each evolve according to KdV dynamics. To facilitate the comparison between NLS and KdV dynamics we choose specific initial conditions that generate perturbations which evolve according to only one of the KdV equations.

Because of the integrability of the KdV (93), a soliton (or solitary wave profile) solution is chosen as a platform for comparison. Due to the special nature of these profiles, i.e. their ability to propagate and retain their structure without breaking by carefully balancing the effects of nonlinearity and dispersion, they also serve as a check on the numerics.

Additionally, the validity of the simulations is also ensured by checking the conservation of wavefunction density and the Hamiltonian. The solitary wave solution is obtained for the KdV (93) which is moving to the right with a velocity of $\mathcal{V}$

[TABLE]

where the index $j$ labels the four eigenvalues corresponding to two left and right movers each. The subscript $j$ on $f$ which is the eigenvalue index should not be confused with the subscript $k$ on quantities like $\rho,\;v,\;\psi$ which is the species index. We choose $j=1,2$ denoting both right movers with speed $+\sqrt{\frac{A\pm B}{2}}$. In order to compare with the coupled NLS, $\delta\rho$ and $\delta v$ need to be calculated from (80). For the fastest mover, $j=1$ eigenvalue $\lambda_{1}$ is chosen.

[TABLE]

We can now obtain $\rho(x,t)$ and $v(x,t)$ from (64), given, $\xi=\epsilon(x-\lambda_{j}t)$ and $\tau=\epsilon^{3}t$.

[TABLE]

Subsequently, $\psi(x,t)$ can be obtained by using the transformation (4). This requires the calculating the integral $\int_{0}^{x}\text{sech}^{2}(ax^{\prime},t)dx^{\prime}=\frac{1}{a}\tanh(ax,t)$. By repeating this exercise for the next larger eigenvalue $\lambda_{j=2}$, we thus obtain the initial profiles for both the components as in (116),(117).

We simulate the dynamics of both the right chiral sectors having speeds $\lambda_{1}$ and $\lambda_{2}$. The initial conditions in the NLS language corresponding to the two eigenvalues $\lambda_{1},\lambda_{2}$ are

[TABLE]

[TABLE]

For simulating the coupled NLS, we use Classical Explicit Method as the time stepping method as outlined in Ref. Taha and Ablowitz, 1984. Keeping in mind the conditions on coupling constants, we choose the following values for the parameters: $g_{1}=1,\;g_{2}=1,\;h=0.5,\;\epsilon=0.2,\;\rho_{01}=1,\;\rho_{02}=0.1$. Thus, $A=1.1$, $B=0.954$ and $C=0.9$. The system size is $L=300$ and spacial axis runs from $x=-150\text{ to }150$ which has been discretized into $n=12000$ steps. Thus, $dx=0.025$ and $dt=7.8125\times 10^{-5}$. The speed of the solitary wave is chosen as $\mathcal{V}=2.5$.

Now that we have the correct initial conditions for both NLS (116) and KdV (113), we let them evolve in time. We study the time evolution of $|\psi_{k}(x,t)|^{2}$. At this stage, the coupled NLS and KdV profiles cannot be compared because they are not the same physical quantities. Since, the KdV problem only involves $f_{j}$’s, while the NLS profile $\psi_{k}^{\lambda_{j}}$ is a combination of $f_{j}$’s multiplied by appropriate coefficients and added to a background, for a meaningful comparison, the KdV variables need to be rescaled via a suitable transformation. The appropriate physical quantities are the density fields and can be obtained by using the transformation (50). We also take care that the two profiles are in the same frame of reference. We choose the lab frame $(x,t)$ for both the profiles. Since, $f_{j}(\xi,\tau)$ is a function of $(\xi-A_{j}\mathcal{V}\tau)$, using the inverse scaling relations we obtain the speeds in the lab frame.

[TABLE]

Thus, the speed of sound in the lab frame is $\varLambda_{j}=\lambda_{j}+A_{j}\mathcal{V}\epsilon^{2}$. For the chosen coupling constants, the parameters are $A_{1}=-0.123$, $B_{1}=1.372$, $A_{2}=-0.462$, $B_{2}=0.728$, $\lambda_{1}=1.013$, $\lambda_{2}=0.270$, $\varLambda_{1}=1.001$, and $\varLambda_{2}=0.224$. Note that opposite signs of $A_{j}$ and $B_{j}$ for both $j=1,2$ ensure an overall negative sign on $f_{j}(x,t)$ (113). Additionally, since the coefficient $C-B=-0.054<0$, $|\psi_{1}^{\lambda_{2}}|^{2}$ profile is expected to be a bump as opposed to the other profiles being a dip. For completeness, we explicitly write out the density fields for both eigenvalues $j=1,2$ for both components $k=1,2$, given, $\xi=\epsilon(x-\lambda_{j}t)$ and $\tau=\epsilon^{3}t$,

[TABLE]

We plot four different time snapshots of the evolution of $|\psi_{k}(x)|^{2}$ through the coupled NLS equation (87) along with the evolution of $\rho_{k}(x)$ through the KdV equation (93) in Fig. (1) for the largest eigenvalue ($\lambda_{1}$) and in Fig. (2) for the second-largest eigenvalue ($\lambda_{2}$) . It is surprising to note that they have quantitative agreement for significant times. Thus, it turns out that if we just evolve the NLS problem with the initial conditions (116),(117) without the knowledge of any scaling or transformations applied, the evolution has a significant match with the independent KdV evolution of (119). Consequently, we see a strong correspondence between the two equations, namely, coupled NLS and the KdV equation.

VI Conclusions

In this paper, we have analyzed the linear and nonlinear problem for the multi-component NLS which is a physically relevant system spanning a broad range of fields. We have systematically studied the qualitative long time dynamics of non-equilibrium profiles. We started with writing a hydrodynamic form. In the linearized regime, we stated and proved a set of theorems. We obtained necessary and sufficient conditions for real speeds of sound that depend solely on the coupling matrix. For the nonlinear problem, using the key mathematical concept of the Fredholm alternative, we show that the coefficients of the KdV dynamics are given in terms of the eigenvectors of the linearised problem. We also discuss numerical protocols to compute the eigenvectors of the coupling matrix individually that also provides us with information on how eigenvalues (sound speeds) and eigenvectors (KdV coefficients) change as the cross-component coupling coefficient varies. This is of high experimental relevance given the tunability of coupling constants. We show compelling evidence of agreement between KdV and multi-component NLS in the nonlinear dynamics using soliton profiles as a platform for comparison This kind of effective mapping shines light on the complex non-equilibrium dynamics of interacting multi-component coupled systems.

The present manuscript investigated the qualitative dynamics of small amplitude perturbations of a trivial state when the speeds of sound (namely the eigenvalues $\mathcal{A}$) are distinct The case of repeated eigenvalues leading to coupled KdV will be investigated in future works. The future outlook also includes generalizations of the coupling matrix $\alpha$. This is important given various physical systems where coupling can vary spatially. Understanding the role of external potential from a rigorous perspective remains unexplored. Although this work is restricted to Hamiltonian systems, it can be extended to open systems which are connected to reservoirs Satpathi et al. (2019) and much remains unexplored in that avenue of driven-dissipative systems.

Acknowledgments

We thank Bernard Deconinck, Chiara D’Errico, Nicolas Pavloff, Urbashi Satpathi and Raghavendra Nimiwal for useful discussions. MK gratefully acknowledges the Ramanujan Fellowship SB/S2/RJN-114/2016 from the Science and Engineering Research Board (SERB), Department of Science and Technology, Government of India and support the Early Career Research Award, ECR/2018/002085 from the Science and Engineering Research Board (SERB), Department of Science and Technology, Government of India. MK would like to acknowledge support from the project 6004-1 of the Indo-French Centre for the Promotion of Advanced Research (IFCPAR). SS acknowledges ICTS for support and hospitality during the S. N. Bhatt Memorial Excellence Fellowship Program 2017.

Appendix A Proof of Theorem 2

Suppose $\alpha$ is positive definite. From Sylvester’s criterion, determinants of all leading principal minors of a positive definite matrix are positive. Then under the supposed ordering, the determinant of the $2\times 2$ leading principal minor, i.e. the $2\times 2$ matrix in the top left corner, $g_{1}g_{2}-h^{2}$ is positive. Hence $h<\sqrt{g_{1}g_{2}}$.

Appendix B Proof of Theorem 3

We first introduce a definition for the symmetric products of $\gamma_{i}$.

Definition 3

We denote the set of $\{\gamma_{i}\}$, obtained by eliminating $\gamma_{k}$ for some $k$, by $\{\gamma_{i}\}-\gamma_{k}$. Similarly, the set obtained adding an element $\gamma_{k}$ is denoted by $\{\gamma_{i}\}+\gamma_{k}$. Suppose $\gamma_{j}=\gamma_{k}$ for some $j,k$ for a list of $\{\gamma_{i}\}$. Let $\beta_{i}$ be the distinct $\gamma_{i}$. Then $\{\gamma_{i}\}=\{\beta_{i}\}+\gamma_{k}$, where $\beta_{i}$ are all distinct. Consequently we describe the replacement of $\gamma_{k}$ with $\gamma_{j}$ in a list by $\{\gamma_{i}\}-\gamma_{j}+\gamma_{k}$.

Keeping in mind that

[TABLE]

we have the following consequence of these definitions.

Theorem 7

[TABLE]

In particular, suppose one of $\gamma_{i}$ equals $1$. Then setting $\{\gamma_{i}\}=\{\beta_{i}\}+1$ we have

[TABLE]

Theorem 3 is proved using the principle of induction. Consider the case when $\alpha,\rho$ are $2\times 2$ matrices. Then

[TABLE]

and so the theorem is true for $N=2$. Furthermore, note that since $\alpha\rho=\rho^{-1}\rho\alpha\rho$, then $\alpha\rho$ and $\rho\alpha$ have the same eigenvalues. Let us now suppose the theorem holds for matrices of some size $n$. We denote the relevant matrices by $\rho^{(n)}$, $\alpha^{(n)}$ so that $\rho^{(n)}\alpha^{(n)}-\mu$ is given by

[TABLE]

where $\mu$ is an eigenvalue of $\rho^{(n)}\alpha^{(n)}$. The characteristic polynomial is

[TABLE]

where

[TABLE]

Since the theorem holds for the $n-$th order matrix, we have

[TABLE]

The matrix $X_{n+1}$ is obtained from $X_{n}$ by the following relation

[TABLE]

where $\vec{1}_{n}$ is a column vector of length $n$ consisting of only ones. The characteristic polynomial then is obtained by computing $\det(X_{n+1})$. This determinant is obtained by a linear combination of the determinants of its co-factor matrices. The co-factor matrices of $X_{n+1}$ are either $X_{n}$ or the matrix obtained by replacing the relevant column of $X_{n}$ by $\vec{1}$. The determinants of those co-factor matrices, obtained by replacing a column of $X_{n}$ by $\vec{1}$, are equal to (up to sign) the determinants of matrices $X_{n}^{(j)}$ where

[TABLE]

Taking into account the signs, we have

[TABLE]

We consider each term on the right-hand side of the above equation individually. By definition $\det(X_{n}^{(j)})$

[TABLE]

where we have used theorem 7. It is straightforward to show that if a set $\{\beta_{i}\}$ contains $n$ elements, then for $m<n$

[TABLE]

Using this relation we have

[TABLE]

On the other hand,$\gamma_{n+1}\det{(X_{n})}$

[TABLE]

Combining the expressions for either of the right-hand side terms of (132) we have

[TABLE]

which proves the theorem for $n+1$ and the statement of the Theorem 3 follows.

Appendix C Proof of Theorem 4

Let $\mathcal{P}(\lambda^{2})$ denote the characteristic polynomial of $\alpha\rho$. We first claim that repeated eigenvalues can only occur when $N>2$. This is easy to see since the characteristic polynomial for $N=2$ is readily computed as

[TABLE]

which has roots

[TABLE]

Since the discriminant is positive for all real values of $\rho_{0i},g_{i},h$, there are no repeated roots when $N=2$.

1 Proof of sufficiency

We now proceed to the general case. We consider first $m$ repeated pairs of $(\rho_{0i}g_{i},\rho_{0i})$ and prove this is a sufficient condition to guarantee a repeated eigenvalue. Let $(\rho^{*}_{0},g^{*})$ denote the common value of $m$ repeated pairs of $(\rho_{0i},g_{i})$. Since the total number of eigenvalue pairs is $N$, there are $N-m$ not-necessarily-repeated pairs. Define

[TABLE]

where the $\gamma_{i}$ are shorthand for the not necessarily repeated pairs. The characteristic polynomial is

[TABLE]

Here $\{\gamma_{i}\}+m\gamma$ in the argument to the symmetric polynomial indicates $m$ repetitions of $\gamma$ in addition to the list of $\gamma_{i}$. We note that

[TABLE]

Also notice that $\mathcal{S}^{N}\left(\left\{\gamma_{i}\right\}+m\gamma\right)=\gamma^{m}\mathcal{S}^{N-m}\left(\left\{{\gamma_{i}}\right\}\right)$ which also follows from the identity above when $l=N$ and recalling $\mathcal{S}^{l}\left(\left\{{\beta_{i}}\right\}\right)=0$ when $l$ is larger than the number of elements in the list $\{\beta_{i}\}$. Substituting this identity into the expression for the characteristic polynomial we have

[TABLE]

Note that $\mathcal{S}^{m-k}\left(\left\{{m\gamma}\right\}\right)$ represents products of $m-k$ $\gamma$’s. Of course there are ${}^{m}C_{m-k}$ ways to choose these products and thus we have

[TABLE]

Substituting the above in to the expression for the characteristic polynomial we have

[TABLE]

Evidently, the above expression is true for $m<N$. Indeed if $m=N$ we have

[TABLE]

or in other words

[TABLE]

Moving ahead with the case $m<N$, we have

[TABLE]

The second term on the left-hand side may be simplified as follows

[TABLE]

The first equality is true since $\mathcal{S}^{N-k-p}\left(\left\{{m\gamma}\right\}\right)=0$ unless $N-k-p\leq m$. To replace the lower limit of the $k-$sum we also need to assure $N-p-m\geq 2$ which implies $p\leq N-m-2$. Thus the only possible exception is when $p=N-m-1$, i.e. the upper limit of the $p-$sum. However, it is easy to see that this term has no contribution for $k=0,1$ since $\mathcal{S}^{m+1-k}\left(\left\{{m\gamma}\right\}\right)=0$ and $k-1=0$ when $k=0,1$ respectively. Finally we obtain the following expression for the characteristic polynomial

[TABLE]

2 Proof of necessity

It is also necessary that at least three pairs of $(\rho_{0i}g_{i},\rho_{0i})$ be equal for the characteristic polynomial to be permanently degenerate. To show this we prove the contrapositive, i.e. we show that if only $m=1,2$ of the $(\rho_{0i}g_{i},\rho_{0i})$ pairs are equal, then the polynomial is not permanently degenerate. Consider first the case $m=1$ i.e. when none of the $(\rho_{0i}g_{i},\rho_{0i})$ are equal. Then a standard implicit function theorem argument applied to

[TABLE]

using the fact that (i) $\mathcal{Q}(\mu)$ has distinct zeros $\rho_{0i}g_{i}$ when $h=0$ and, (ii) $\partial\mathcal{Q}/\partial\mu$ evaluated at $h=0$ is non-zero (due to the distinct values of $\rho_{0i},g_{i}$), we have an open neighborhood of $h=0$ where there are $N$ distinct roots to the polynomial $\mathcal{Q}$ and so the polynomial cannot be permanently degenerate. We now consider the case $m=2$. Suppose only two of the $(\rho_{0i}g_{i},\rho_{0i})$ are equal. Let the common value be $(\rho^{*}_{0}g^{*},\rho^{*}_{0})$. By assumption this common value is distinct from all remaining $N-2$ $(\rho_{0i}g_{i},\rho_{0i})$ pairs. Factoring as we did previously, but for the case $m=2$, we have

[TABLE]

This polynomial has a root $\mu=\rho^{*}_{0}(g^{*}-h)$ and the roots of

[TABLE]

Applying the implicit function theorem to the above polynomial we have that in an open neighborhood of $h=0$, the above polynomial has distinct roots. The only remaining possibility is that $\rho^{*}_{0}(g^{*}-h)$ itself is a root of this polynomial. But this leads to the following expression valid for all suitable $h$.

[TABLE]

which for $h=0$ is

[TABLE]

which is not possible since $\rho^{*}_{0}g^{*}\neq\rho_{0i}g_{i}$. Consequently when $m=1,2$ the characteristic polynomial cannot be permanently degenerate.

Appendix D Proof of Theorem 6

From Theorem 5, it suffices to find eigenvectors for the matrix $\alpha\rho$ with eigenvalues $\lambda^{2}=\rho^{*}_{0}(g^{*}_{0}-h)$. A straightforward computation gives the components of $\alpha\rho-\lambda^{2}$

[TABLE]

Since the matrix $\mathcal{A}$ is permanently degenerate, $m+1$ pairs of $(\rho_{0i}g_{i},\rho_{0i})$ are equal to $(\rho^{*}_{0}g^{*},\rho^{*}_{0})$. From this, it follows that $m+1$ columns of the matrix $(\alpha\rho-\lambda^{2})$ are parallel. Indeed all elements of such columns are $\rho^{*}_{0}h$. We identify these columns by $i_{k},\ k=1,2,\ldots m+1$. Define the $i-$th element of the vector $q^{(k)}$ by

[TABLE]

Then $\alpha\rho q^{(k)}=\rho^{*}_{0}(g^{*}-h)q^{(k)}$. Using the construction of the previous theorem, we obtain the associated eigenvector for $\mathcal{A}$.