Inverse scattering transform for two-level systems with nonzero background

TL;DR

This paper develops an inverse scattering transform for the Maxwell-Bloch system with nonzero background, explicitly constructs soliton solutions, and analyzes their properties and stability under various initial conditions.

Contribution

It introduces a formalism for the inverse scattering transform with nonzero background and derives explicit soliton, periodic, and rational solutions for the Maxwell-Bloch system.

Findings

Pure background states generally do not exist with nonzero background.

Explicit one-soliton solutions are derived in determinant form.

Solutions are stable only when initial population inversion is negative.

Abstract

We formulate the inverse scattering transform for the scalar Maxwell-Bloch system of equations describing the resonant interaction of light and active optical media in the case when the light intensity does not vanish at infinity. We show that pure background states in general do not exist with a nonzero background field. We then use the formalism to compute explicitly the soliton solutions of this system. We discuss the initial population of atoms and show that the pure soliton solutions do not correspond to a pure state initially. We obtain a representation for the soliton solutions in determinant form, and explicitly write down the one-soliton solutions. We next derive periodic solutions and rational solutions from the one-soliton solutions. We then analyze the properties of these solutions, including discussion of the sharp-line and small-amplitude limits, and thereafter show that the two limits do not commute. Finally, we investigate the behavior of general solutions, showing that solutions are stable (i.e., the radiative parts of solutions decay) only when initially atoms in the ground state dominant, i.e., initial population inversion is negative.

Equations

$$\partial q/\partial z=-\int P(t,z,\lambda )g(\lambda )\mathrm{d}\lambda ,$$

$$\partial D/\partial t=2\mathrm{Re}(q^{\ast }P),\partial P/\partial t-2i\lambda P=-2qD.\text{vspace}\ast -1ex$$

$${\rho }_t=[i\lambda {\sigma }_3+Q,\rho ],$$

$$Q_z=-\frac{1}{2}\int [{\sigma }_3,\rho ]g(\lambda )\mathrm{d}\lambda ,$$

$$Q(t,z)=\left(\begin{array}{cc}0& q\\ -q^{\ast }& 0\end{array}\right),\rho (t,z,\lambda )=\left(\begin{array}{cc}D& P\\ P^{\ast }& -D\end{array}\right),$$

$${\sigma }_1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),{\sigma }_2=\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right),{\sigma }_3=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right).$$

$${\varphi }_t=(i\lambda {\sigma }_3+Q)\varphi ,$$

$${\varphi }_z=V\varphi ,$$

$$V(t,z,\lambda )=\frac{i\pi }{2}{\mathcal{H}}_{\lambda }[\rho (t,z,{\lambda }^{\prime })g({\lambda }^{\prime })],$$

$${\mathcal{H}}_{\lambda }[f({\lambda }^{\prime })]=\frac{1}{\pi }\int -\frac{f({\lambda }^{\prime })}{{\lambda }^{\prime }-\lambda }\mathrm{d}{\lambda }^{\prime }.$$

$$\rho (t,z,\lambda )={\mathrm{e}}^{X_{o}t}C{\mathrm{e}}^{-X_{o}t},$$

$$\gamma (\lambda )=\mathrm{sign}(\lambda )\sqrt{{\lambda }^2+q_o^2},\lambda \in \mathbb{R}.$$

$$Y_o=I+(i/\zeta ){\sigma }_3{Q}_o,$$

$$\zeta =\lambda +\gamma ,$$

$$\rho (t,z,\lambda )=Y_o{\mathrm{e}}^{i\gamma t{\sigma }_3}{\rho }_o{\mathrm{e}}^{-i\gamma t{\sigma }_3}{Y}_o^{-1},$$

$$\lambda =\frac{1}{2}(\zeta -q_o^2/\zeta ),\gamma =\frac{1}{2}(\zeta +q_o^2/\zeta ).$$

$$Y_o^{-1}=[I-(i/\zeta ){\sigma }_3{Q}_o]/\det Y_o,\det Y_o=2\gamma /\zeta =1+q_o^2/{\zeta }^2.$$

$$\rho (t,z,\lambda )=\mathbf{h}\cdot \mathbf{\rho },$$

$${\rho }_j(t,z,\lambda )=Y_o{\mathrm{e}}^{i\gamma t{\sigma }_3}{\sigma }_j{\mathrm{e}}^{-i\gamma t{\sigma }_3}{Y}_o^{-1},$$

$$\displaystyle\rho_{1}=\frac{1}{2\gamma\zeta}\big{(}\zeta^{2}\,\mathrm{e}^{i\gamma t\sigma_{3}}\sigma_{1}\mathrm{e}^{-i\gamma t\sigma_{3}}-2\zeta\mathop{\rm Im}\nolimits(\mathrm{e}^{-2i\gamma t}q_{o})\,\sigma_{3}+\mathrm{e}^{-i\gamma t\sigma_{3}}Q_{o}\sigma_{1}Q_{o}\mathrm{e}^{i\gamma t\sigma_{3}}\big{)}\,,$$

$$\displaystyle\rho_{2}=\frac{1}{2\gamma\zeta}\big{(}\zeta^{2}\,\mathrm{e}^{i\gamma t\sigma_{3}}\sigma_{2}\mathrm{e}^{-i\gamma t\sigma_{3}}-2\zeta\mathop{\rm Re}\nolimits(\mathrm{e}^{-2i\gamma t}q_{o})\,\sigma_{3}+\mathrm{e}^{-i\gamma t\sigma_{3}}Q_{o}\sigma_{2}Q_{o}\mathrm{e}^{i\gamma t\sigma_{3}}\big{)}\,,$$

$${\rho }_3=\frac{1}{\gamma }(\lambda {\sigma }_3-i{Q}_o).$$

$$\displaystyle[\sigma_{3},\rho_{1}]=(1/\gamma\zeta)\,\mathrm{e}^{-i\gamma t\sigma_{3}}\big{(}i\zeta^{2}\sigma_{2}+\sigma_{3}\,Q_{o}\sigma_{1}Q_{o}\big{)}\,\mathrm{e}^{i\gamma t\sigma_{3}}\,,$$

$$\displaystyle[\sigma_{3},\rho_{2}]=(1/\gamma\zeta)\,\mathrm{e}^{-i\gamma t\sigma_{3}}\big{(}-i\zeta^{2}\sigma_{1}+\sigma_{3}\,Q_{o}\sigma_{2}Q_{o}\big{)}\,\mathrm{e}^{i\gamma t\sigma_{3}}\,,$$

$$[{\sigma }_3,{\rho }_3]=-(i/\gamma )[{\sigma }_3,Q_o].$$

$$\int [{\sigma }_3,{\rho }_3]g(\lambda )\mathrm{d}\lambda =-i{w}_o[{\sigma }_3,Q_o],$$

$$w_o=\int g(\lambda )/\gamma \mathrm{d}\lambda .$$

$$q(t,z)=q_o(z)=q_o(0){\mathrm{e}}^{-i{h}_3{w}_{o}z},\rho (t,z,\lambda )=h_3(\lambda {\sigma }_3-i{Q}_o)/\gamma .$$

$$q(t,z)=q_o,\rho (t,z,\lambda )=h_3(\lambda {\sigma }_3-i{Q}_o)/\gamma .$$

$$\gamma (\lambda )=(q_o^2+{\lambda }^2{)}^{1/2}.$$

$$\gamma (\lambda )=\sqrt{r_1{r}_2}{\mathrm{e}}^{i\Theta },\Theta =({\theta }_1+{\theta }_2)/2+m\pi ,$$

$${\Gamma }^{+}=\{\zeta \in \mathbb{C}:(\mathrm{Im}\zeta )(|\zeta {|}^2-q_o^2)>0\},{\Gamma }^{-}=\{\zeta \in \mathbb{C}:(\mathrm{Im}\zeta )(|\zeta {|}^2-q_o^2)<0\}.$$

$$Y_{\pm }=I+(i/\zeta ){\sigma }_3{Q}_{\pm },$$

$$X_{\pm }{Y}_{\pm }=Y_{\pm }i\gamma {\sigma }_3,$$

Full text

\cftsecindent

0em \cftsubsecindent1.4em

Inverse scattering transform for two-level systems withnonzero background

Gino Biondini1, Ildar Gabitov2, Gregor Kovačič3 and Sitai Li4

Abstract

We formulate the inverse scattering transform for the scalar Maxwell-Bloch system of equations describing the resonant interaction of light and active optical media in the case when the light intensity does not vanish at infinity. We show that pure background states in general do not exist with a nonzero background field. We then use the formalism to compute explicitly the soliton solutions of this system. We discuss the initial population of atoms and show that the pure soliton solutions do not correspond to a pure state initially. We obtain a representation for the soliton solutions in determinant form, and explicitly write down the one-soliton solutions. We next derive periodic solutions and rational solutions from the one-soliton solutions. We then analyze the properties of these solutions, including discussion of the sharp-line and small-amplitude limits, and thereafter show that the two limits do not commute. Finally, we investigate the behavior of general solutions, showing that solutions are stable (i.e., the radiative parts of solutions decay) only when initially atoms in the ground state dominant, i.e., initial population inversion is negative.

1 Introduction

Resonant interaction between pulses of monochromatic light and an active optical medium is one of the most useful phenomena studied in applied optics, having lead to such indispensable devices as lasers and optical amplifiers [1, 2, 3, 4, 5]. The optical medium is often assumed to have two working atomic levels, with the transition frequency between these two levels roughly equal to the frequency of the impinging light, which is known as the two-level medium approximation [3]. In addition to lasing and amplification, the interaction between light and two-level media has given rise to several classical nonlinear-optics effects including self-induced transparency [2, 6, 7, 8] (undisturbed propagation of pulses with sufficiently high amplitudes through the medium and absorption of weaker pulses), superfluorescence (generation of short optical pulses from material polarizability fluctuations in an excited medium [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]), optical nutation (a phenomenon during which the material properties exhibit oscillations akin to nutation in rigid bodies [3]), and photon echo [20, 21, 9, 22, 23] (generation of a third pulse from the injection of a pair of pulses into the medium).

The wealth of physical phenomena exhibited by light interacting with two-level active optical media has provided an abundant source of fundamental mathematical problems used in their description. This description is afforded by various versions of the two-level Maxwell-Bloch equations (MBEs) [12, 24, 25, 26, 27, 3]. Frequently, such equations are derived by assuming unidirectional propagation and using the slowly-varying envelope approximation, i.e., separating the fast-oscillating plane wave representing the color of the light from the slowly-varying envelope representing its amplitude and phase, and the rotating wave approximation, i.e., averaging over the fast oscillations [27]. Mathematically, these approximations are equivalent to the more general method of multiple scales [28]. The properties of the material are typically described in terms of the quantum density matrix that yields a Bloch vector analogous to that used in spin dynamics [24]. Depending on the values of their parameters and inputs, the two-level MBEs exhibit a rich variety of dynamical regimes. While some regimes of laser operation are chaotic [29, 30, 31, 32, 33], and can be described by the Lorenz equations [34, 35], or are even turbulent [36, 37, 38, 39], and have a finite-dimensional attractor [40, 41, 42, 43] or a slow manifold [44], other regimes of light propagating through rarefied gases may be approximated by completely-integrable, soliton-type equations [9, 10, 45, 46, 47, 48, 49, 50, 51, 17, 52, 4].

The integrable MBEs are derived under the above-mentioned approximations, and have addressed phenomena including self-induced transparency [9, 10, 45], area theorem [46], photon echo [48], amplification [47, 49, 50], and superfluorescence [51, 17, 52]. Typically, these phenomena are addressed from the viewpoint of the initial-boundary-value or signaling problem for a finitely or infinitely long narrow tube containing the active medium, into which a narrow optical pulse is injected at one end, and whose initial state is given far in the past. The MBEs can be represented in terms of a Lax pair [45, 52], thus, the initial-value problem can be solved using the inverse-scattering transform (IST) [45, 53, 52, 54]. The temporal piece of the Lax pair is in the AKNS form [55, 56], so just as for the nonlinear Schrödinger equation appearing in [57, 58], which stretches from the infinite past to the infinite future, and pulses are evolved along the spatial variable, describing its propagation from the location, where the pulse is launched into the medium. The evolution operator in the Lax pair of the MBEs has an unusual form containing a Hilbert transform, which complicates the derivation and form of the equations describing the evolution of the spectral data, as well as reflects the fact that the phenomena described by the MBEs are irreversible [45].

All the above classic works on the IST for the integrable two-level MBEs addressed narrow optical pulses, i.e., pulses that decay both forward and backward in time. A few studies have discussed pulses riding on top of a continuous-wave (CW) background [59, 60, 61, 62, 63, 64, 65, 66, 67, 68]. However, these phenomena have so far attracted much less attention. In this paper, we make a systematic attempt to fill this gap by developing the IST for studying pulses that asymptotically approach plane waves with equal amplitudes and frequencies in both backward and forward times, i.e., for MBEs with nonzero background (NZBG).

The direct and inverse scattering components of our study closely parallel those in our work on the focusing nonlinear Schrödinger (NLS) equation with nonzero boundary conditions (NZBC) [69]. More generally, the study of nonlinear systems with NZBG has also received considerable interest thanks in part to its connections to the theory of modulational instability, rogue waves and integrable turbulence [70, 71, 72, 73, 74, 75, 76, 77, 78, 79]. The continuous spectra for the scattering problems of both the NLS and MBEs comprise the real axis plus a symmetric interval along the imaginary axis with length twice that the amplitude of the CW background. This symmetric interval is the branch cut for the double-sheeted Riemann surface on which both the scattering and inverse-scattering take place.

The main difficulty of the IST for MBEs lies in evaluating the evolution of the spectral data. Deriving this evolution is aided by the Liouville form of the equation for the material variables and the dual role of the spectral parameter as the frequency detuning [45, 52]. Especially intricate is finding the time-dynamics of the reflection coefficients, which requires careful evaluation of Fourier-type integrals using the Riemann-Lebesgue lemma and results in an inhomogeneous linear differential equation. While the result is somewhat similar to that derived in [52] for MBEs with vanishing background, the details of the derivation are significantly more involved. We find that already the single-soliton solutions of the MBEs with NZBG bring some surprises. The electric field envelope of their counterparts on zero background (ZBG) has a very simple, sech-like shape, possibly multiplied by an oscillating complex phase factor. On NZBG, this shape only remains for the solitons corresponding to the eigenvalue lying on the imaginary axis above the branch cut (and, of course, its complex conjugate, which we ignore from now on). An eigenvalue in the general position in the upper half-plane gives rise to a traveling breather-like structure, with an exponentially localized envelope containing sinusoidal ripples that move faster than the envelope in the same direction. The entire structure is again multiplied by an oscillating complex phase factor. Away from its center, the structure rapidly approaches the background plane wave. Moreover, as the corresponding eigenvalue approaches the branch cut, the width of the solution envelope spreads to infinity, and what remains is a periodic traveling wave. Moreover, the eigenvalue at the branch point produces a traveling wave with a rational form. If we consider the spectral parameter to be a bifurcation parameter, these families of solitons represent the unfolding of the solutions around the spectral singularity at the branch point. The four types of one-soliton solutions have been previously reported in [80]. Finally, we investigate soliton solutions in the case in which the frequency of the impinging light is detuned from the peak frequency of the atomic transition in the medium, and describe further unfolding properties of the solutions also in terms of the detuning as a parameter.

The outline of this work is the following: In Section 2 we present the MBEs and its Lax pair, and we discuss its background solutions and the treatment of the branch cut in the complex plane in the formulation of the IST. In Section 3 we formulate the direct scattering. We derive all scattering data including the transmission and reflection coefficients, discrete spectrum and norming constants. We then discuss the symmetries and asymptotics of eigenfunctions and scattering data, and the asymptotic behavior of the density matrix. In Section 4 we determine propagation equations for all the relevant quantities, including the scattering data and the continuous-wave background. In Section 5 we formulate the inverse problem in terms of a Riemann-Hilbert problem (RHP). We then reconstruct general solutions of the MBEs from the solutions of the RHP. Moreover, we derive a “trace” formula which recovers the analytic scattering coefficient, and derive a “theta” condition which yields the asymptotic phase difference of the solution. In particular, we write down an explicit formula for pure $N$-soliton solutions, for which the reflection coefficient is identically zero. In Section 6, we focus on one-soliton solutions. We show that there are four types of solitons depending on the location of the corresponding discrete eigenvalue. They are traveling-wave solitons, oscillatory solitons, periodic solutions and rational solutions. We also discuss their properties in detail and we analyze the stability of soliton solutions. In Section 7 we discuss soliton solutions resulting from a shifted Lorentzian as the spectral-line shape. Such a shift originates from a detuning between the resonance frequency of the atoms (corresponding to the transition between the ground state and the excited state) and the frequency of the incident light pulse. Importantly, we show that the presence of such a shift produces markedly different properties for the resulting soliton solutions. In particular, we will show that some of the solutions actually become rogue-wave in character, in the sense that they appear to arise from virtually nothing, attain a peak value which is several times larger than the average, and eventually revert to the background. We conclude this work with some final remarks in Section 8.

Auxiliary material is confined to the background. In Appendix A.1 we describe some frequently used notation. In Appendix A.2 we derive the symmetries of various quantities appearing in the direct problem of the IST. In Appendix A.3 we compute the asymptotics of the eigenfunctions. In Appendix A.4 we derive the explicit formula for an auxiliary matrix that appears in the propagation part of the IST. In Appendix A.5 we derive the propagation equations for the norming constants and for the reflection coefficient. In Appendix A.6 we show how to solve the RHP that appears in the inverse problem. In Appendix A.7 we derive the trace formula and the “theta” condition for the solutions of the MBEs with NZBG. In Appendix A.8 we calculatie the auxiliary matrix explicitly with a Lorentzian as the spectral-line shape. In Appendix A.9 we briefly discuss an alternative formulation of the IST, and prove the equivalence between the two versions of the IST.

2 Preliminaries

2.1 The Maxwell-Bloch system of equations and its Lax pair

Up to rescalings of the dependent and independent variables, the scalar MBEs are, in component form

[TABLE]

where $z=z_{\mathrm{lab}}$ is the propagation distance, $t=t_{\mathrm{lab}}-z_{\mathrm{lab}}/c$ is a retarded time ($c$ is the speed of light in vacuum), $q(t,z)$ is the optical field, i.e., the electric field envelope corresponding to the transitions between the ground state and the excited state, level inversion $D(t,z,\lambda)$ and medium polarization envelope $P(t,z,\lambda)$ are the entries of a $2\times 2$ Hermitian matrix representing the density matrix, and $g(\lambda)$ is the (known) shape of the spectral line. All integrals in this work run from $-\infty$ to $\infty$ unless indicated otherwise.

The MBEs (1) and (2) can be written compactly as

[TABLE]

where subscripts $t$ and $z$ denote partial differentiation, $[A,B]=AB-BA$ is the matrix commutator, the $2\times 2$ optical matrix and density matrix are

[TABLE]

and $\sigma_{1},\sigma_{2},\sigma_{3}$ are the Pauli matrices

[TABLE]

The Lax pair for the MBEs (1) and (2) is:

[TABLE]

where Eq. (6) and $q(t,z)$ are referred to as the scattering problem and as the scattering potential, respectively, and $\phi=\phi(t,z,\lambda)$ is the wave function, with

[TABLE]

where $\mathcal{H}_{\lambda}[f(\lambda^{\prime})]$ is the Hilbert transform of the function $f(\lambda)$ defined by

[TABLE]

Throughout this work, primes in integration variables do not denote differentiation.

The IST for the MBEs (3) and (4) with $q(t,z)\to 0$ as $t\to\pm\infty$ — hereafter referred to as ZBG — was carried out in [45, 52] (while some earlier results appeared in [9, 10]). Here we develop the IST for the scalar MBEs (3) and (4) with NZBG. Namely, we solve system (3) and (4) with $q\to q_{\pm}$ as $t\to\pm\infty$. We assume $|q_{\pm}|=q_{o}>0$. Similarly to what has been done in [69] for the focusing nonlinear Schrödinger (NLS) equation with NZBC, we will formulate the IST such that it allows one to take the reduction $q_{o}\to 0$ explicitly throughout.

Note that: (i) System (3) and (4) implies $\partial_{t}({\mathrm{tr}}\rho)=\partial_{t}(\det\rho)=0$, (ii) System (3) and (4) is invariant under the transformation $\rho(t,z,\lambda)\mapsto\tilde{\rho}(t,z,\lambda)=\rho(t,z,\lambda)+f(z,\lambda)\,I$, where $f$ is any scalar function independent of $t$. Hence, without loss of generality we can take $\rho$ to be traceless and with determinant equal to $-1$. Moreover, we assume that this normalization has been made for all $z\geqslant 0$.

2.2 Background solutions

Before we formulate the IST, we investigate whether there exist exact “constant” solutions of the MBEs (3) and (4) with NZBG. That is, we look for solutions $q(t,z)=q_{o}(z)$ $\forall t\in\mathbb{R}$, with $q_{o}\neq 0$. The right-hand side (RHS) of Eq. (3) contains the matrix $X_{o}=i\lambda\sigma_{3}+Q_{o}$, where $Q_{o}$ is the same as $Q$, but with $q_{o}$ instead of $q$. Since $Q_{o}$ is independent of $t$, the solution of Eq. (3) is simply

[TABLE]

for any $2\times 2$ matrix $C$ independent of $t$. The eigenvalues of $X_{o}$ are $\pm i\gamma$, where $\gamma^{2}=\lambda^{2}+q_{o}^{2}$. For now we limit our discussion to real values of $\lambda$. (The extension of $\gamma$ to complex values of $\lambda$ will be discussed in Section 2.3.) In order to be able to take the limit $q_{o}\to 0$ continuously we need to choose

[TABLE]

Of course, with this choice $\gamma$ has a sign discontinuity across $\lambda=0$. We can write an eigenvector matrix of $X_{o}$ compactly as

[TABLE]

where we introduced the shorthand notation

[TABLE]

which will be used throughout this work. Explicitly, $X_{o}Y_{o}=Y_{o}\,i\gamma\sigma_{3}$. Thus,

[TABLE]

where $\rho_{o}=Y_{o}^{-1}CY_{o}$. The inverse transformation to Eq. (12) is

[TABLE]

For future reference, note that

[TABLE]

Also note that $\det\rho=\det\rho_{o}$ and ${\mathrm{tr}}\rho={\mathrm{tr}}\rho_{o}$. Hence, equations $\rho^{\dagger}=\rho$, ${\mathrm{tr}}\rho=0$ and $\det\rho=-1$ imply that $\rho_{o}^{\dagger}=\rho_{o}$, ${\mathrm{tr}}\rho_{o}=0$ and $\det\rho_{o}=-1$. Therefore, we can write $\rho_{o}=\mathbf{h}\cdot{\boldsymbol{\sigma}}$, where ${\boldsymbol{\sigma}}=(\sigma_{1},\sigma_{2},\sigma_{3})^{T}$ and $\mathbf{h}\in\mathbb{R}^{3}$, with $\mathbf{h}\cdot\mathbf{h}=-\det\rho=1$. The superscript $T$ denote matrix transpose. So, we have

[TABLE]

where ${\boldsymbol{\rho}}=(\rho_{1},\rho_{2},\rho_{3})^{T}$ and

[TABLE]

with $\rho_{j}^{\dagger}=\rho_{j}$, ${\mathrm{tr}}\rho_{j}=0$, and $\det\rho_{j}=-1$. Explicitly,

[TABLE]

All these results reduce to the correct behavior in the ZBG limit (i.e., $q_{o}\to 0$). In particular, $\gamma\to\lambda$, $\zeta\to 2\lambda$ and $Y\to I$ when $q_{o}\to 0$, implying ${\boldsymbol{\rho}}\to\mathrm{e}^{i\lambda t\sigma_{3}}{\boldsymbol{\sigma}}\,\mathrm{e}^{i\lambda t\sigma_{3}}$ in the limit $q_{o}\to 0$. Remarkably, even when $q_{o}\neq 0$, $\rho_{3}$ is always independent of $t$ as in the case with ZBG, even though the ODEs (3) are more complicated when $Q_{o}\neq 0$.

Now we insert the behavior described by Eqs. (16), (17) and (18) into Eq. (4). It is trivial to see that $[\sigma_{3},\rho]=\mathbf{h}\cdot[\sigma_{3},{\boldsymbol{\rho}}]$. Moreover, direct calculation yields

[TABLE]

The third commutator above is $t$-independent, leaving the field invariant. Explicitly,

[TABLE]

where

[TABLE]

On the other hand, the two commutators in Eqs. (19) and (20) contain $t$-dependent oscillating exponentials, which in general do not cancel even upon integration. Hence, in general the material polarization does not preserve the $t$-independence of the optical field, i.e., they are not consistent with a constant solution $q(t,z)=q_{o}(z)$. Accordingly, the only self-consistent background solution is the one obtained with $h_{1}=h_{2}=0$ and $h_{3}=\pm 1$, resulting in

[TABLE]

With a generic spectral-line shape $g(\lambda)$, it is clear from Eq. (22) that $w_{o}\neq 0$. However, the simplest situation occurs when the distribution of atoms is even with respect to the normalized resonant frequency $\lambda=0$. In such a case, the evenness of the function $g(\lambda)$ and oddness of the function $\gamma(\lambda)$ imply $w_{o}=0$. The above background solution (23) reduces to

[TABLE]

Importantly, note that $\rho(t,z,\zeta)$ is discontinuous at $\lambda=0$. Also recall that in the case of the MBEs with ZBG, $h_{3}=-1$ implies the atoms are in the ground state, whereas, $h_{3}=1$ implies all the atom are in the excited state. However, in the case with NZBG, the situation is seemingly more complicated. Because of the presence of background radiation, a polarization-free state does not exist in the background solutions (24), i.e., the medium is always polarized. This polarization is reflected in the off-diagonal entries of the density matrix $\rho$, which are never vanish, i.e., $P=-ih_{3}q_{o}/\gamma$, which arises from the background radiation $q_{o}(z)$. In fact, later we will show that this statement holds for general solutions as well. Nonetheless, we will also show that a natural decomposition still exists for the density matrix. We will also show in Section 3.4 that, when $q(t,z)$ depends on $t$, as long as $q\to q_{\pm}$ as $t\to\pm\infty$, the asymptotic behavior of the density matrix [i.e., the solutions of Eq. (3)] is consistent with the above expressions obtained with a background value of the potential.

2.3 Riemann surface and uniformization variable

Before we can start to develop the IST for the MBEs (3) and (4), it is useful to discuss the branching of the eigenvalues of the scattering problem. The asymptotic scattering problem as $t\to\pm\infty$ is $\phi_{t}=X_{\pm}\,\phi$, where $X_{\pm}=i\lambda\sigma_{3}+Q_{\pm}$. The eigenvalues of $X_{\pm}$ are $\pm i\gamma$, with $\gamma^{2}=q_{o}^{2}+\lambda^{2}$, where $q_{o}=|q_{\pm}|$. As with the NLS equation [81], to address the branching, we introduce the two-sheeted Riemann surface defined by the complex square root

[TABLE]

Specifically, the Riemann surface is obtained by introducing two copies of the complex plane, called $\mathbb{C}_{I}$ and $\mathbb{C}_{\mathrm{II}}$, in which $\gamma(\lambda)$ takes either of the two possible signs of the complex square root. The branch points are the values of $\lambda$ for which $\gamma(\lambda)=0$, i.e., $\lambda=\pm iq_{o}$. Let $\lambda+iq_{o}=r_{1}\,\mathrm{e}^{i\theta_{1}}$ and $\lambda-iq_{o}=r_{2}\,\mathrm{e}^{i\theta_{2}}$. In this way we can write

[TABLE]

where $m=0,1$ respectively on sheet I and II. Now take $-\pi/2\leqslant\theta_{j}<3\pi/2$ for $j=1,2$. With these conventions, the discontinuity of $\gamma$ (which defines the location of the branch cut) occurs on the segment $i[-q_{o},q_{o}]$. The Riemann surface is then obtained by gluing the two copies of the complex plane along the cut. Formally, $\hat{\mathbb{C}}=\mathbb{C}_{\mathrm{I}}\cup\mathbb{C}_{\mathrm{II}}$. Along the real $\lambda$ axis we then have $\gamma(\lambda)=\pm\mathop{\rm sign}\nolimits(\lambda)\sqrt{q_{o}^{2}+\lambda^{2}}$, where the plus and minus signs apply respectively on sheet I and sheet II of the Riemann surface, and where the square root sign denotes the principal branch of the real-valued square root function. For later use, we use the subscripts $\mathrm{I}$ and $\mathrm{II}$ to denote that a quantity is evaluated on the first or second sheet, respectively.

The reason why we take the branch cut along $[-iq_{o},iq_{o}]$ is that, in this way, $\gamma(\lambda)$ reduces to (10) for $\lambda\in\mathbb{R}$ which in turn reduces to $\lambda$ as $q_{o}\to 0$ . Not only does this allow us to take the limit $q_{o}\to 0$ throughout, but it also allows us to choose the initial state for the MBEs motivated by the physical intuition gained from the case of ZBG. Of course, one could just as well formulate the IST with the branch cut on $i\mathbb{R}\backslash[-iq_{o},iq_{o}]$. In Appendix A.9, we show that doing so is equivalent to a redefinition of the initial states of the MBEs. Note that Eqs. (14) allows us to can express all $\lambda$ dependence in terms of $\zeta$, and we will do so throughout this work.

It will also be convenient to introduce a uniformization variable. As in [82], this is easily done by defining $\zeta=\lambda+\gamma$ as in Eq. (12), which is inverted by Eq. (14). Let $C_{o}$ be the circle of radius $q_{o}$ centered at the origin in the complex $\zeta$-plane. With these definitions: (i) the branch cut on either sheet is mapped onto $C_{o}$; (ii) in particular: $\zeta(\pm iq_{o})=\pm iq_{o}$ from either sheet, $\zeta(0_{\mathrm{I}}\pm\varepsilon)\to\pm q_{o}$ and $\zeta(0_{\mathrm{II}}\pm\varepsilon)\to\mp q_{o}$ as $\varepsilon\to 0^{+}$; (iii) $\mathbb{C}_{\mathrm{I}}$ is mapped onto the exterior of $C_{o}$; (iv) $\mathbb{C}_{\mathrm{II}}$ is mapped onto the interior of $C_{o}$; (v) in particular: $\zeta(\infty_{I})=\infty$ and $\zeta(\infty_{\mathrm{II}})=0$; (vi) the first/second quadrant of $\mathbb{C}_{\mathrm{I}}$ are respectively mapped into the first/second quadrant outside $C_{o}$; (vii) the first/second quadrant of $\mathbb{C}_{\mathrm{II}}$ are respectively mapped into the second/first quadrant inside $C_{o}$; (viii) note also $\zeta_{\mathrm{I}}\zeta_{\mathrm{II}}=-q_{o}^{2}$. The uniformization map can be seen clearly from Fig. 1.

Note that, in general, the density matrix $\rho(t,z,\lambda)$ is only defined for $\lambda\in\mathbb{R}$. Writing it as a function of the uniformization parameter $\zeta$, we can evaluate it for all $\zeta\in\mathbb{R}$, specifically for $|\zeta|>q_{o}$ on sheet I, and for $|\zeta|<q_{o}$ on sheet II, where obviously $\rho(t,z,-q_{o}^{2}/\zeta)=\rho(t,z,\zeta)$ [since $\rho(t,z,\lambda)$ is a single-valued function of $\lambda$]. But since a priori we do not have any guarantee that $\rho(t,z,\lambda)$ can be analytically extended off the real $\lambda$ axis, we will not be allowed to evaluate $\rho(t,z,\zeta)$ off the real $\zeta$ axis either.

We next discuss the difference between the two sheets. Recall in the case with ZBG, the background density-matrix solution is $\rho=h_{3}\,\sigma_{3}$ with $h_{3}=\pm 1$. The value $h_{3}=\mp 1$ indicates that all atoms are in the ground state or in the excited state, respectively. Therefore the background solution (24) with NZBG should reduce to $\rho=h_{3}\,\sigma$ as $q_{o}\to 0$. On sheet I, this is obviously true. However, on sheet II, the solution (24) reduces to $\rho=-h_{3}\,\sigma_{3}$ as $q_{o}\to 0$. Hence, to maintain consistency with the ZBG reduction, one must choose opposite values for $h_{3}$ on sheet II in the background solution (24). We will discuss this situation in detail in Section 3.4.

Hereafter, with some abuse of notation we will rewrite all the $\lambda$ dependence as dependence on $\zeta$. Note that Eq. (26) implies that \mathop{\rm Im}\nolimits\,\gamma\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\m@th\displaystyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\textstyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\scriptstyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\scriptscriptstyle\hfil#\hfil\cr>\crcr<\crcr} }}}0 respectively in $\Gamma^{+}$ and $\Gamma^{-}$, where

[TABLE]

These two regions can be seen in Fig. 1 with gray and yellow colors. As we show next, this property will determine the analyticity regions of the Jost eigenfunctions.

3 Direct problem

We are now ready to start formulating the IST for the MBEs (3) and (4) with NZBG.

3.1 Jost solutions, analyticity and scattering matrix

We first go back to the asymptotic scattering problem, i.e., the ODE (6) with $t\to\pm\infty$. We find it convenient to consider the wave function $\phi(t,z,\lambda)$ as a $2\times 2$ matrix (instead of a 2-component vector). On either sheet, we can write the asymptotic eigenvalue matrix and eigenvector matrices as $i\gamma\sigma_{3}$ and

[TABLE]

respectively, where $I$ denotes the $2\times 2$ identity matrix and $Q\to Q_{\pm}$ as $t\to\pm\infty$, so that

[TABLE]

with

[TABLE]

Similarly to what happened in Section 2.2, one can easily calculate that

[TABLE]

By definition [69], the continuous spectrum $\Sigma_{\lambda}$ consists of all values of $\lambda$ (on either sheet) such that $\gamma(\lambda)\in\mathbb{R}$. Based on the above choices (25) [equivalent to Eq. (10)], it is $\Sigma_{\lambda}=\mathbb{R}\cup i[-q_{o},q_{o}]$. The corresponding set in the complex $\zeta$-plane is $\Sigma_{\zeta}=\mathbb{R}\cup C_{o}$ from the definition (12). Hereafter we will omit the subscripts on $\Sigma$; the intended meaning should be clear from the context.

We expect that, for all $\zeta\in\Sigma$, the scattering problem (6) admits solutions $\phi_{\pm}$ which behave as

[TABLE]

As in Ref. [69], we introduce modified eigenfunctions by removing the asymptotic exponential oscillations:

[TABLE]

so that $\lim_{t\to\pm\infty}\mu_{\pm}(t,z,\zeta)=Y_{\pm}$. One can then formally integrate the resulting ODEs for $\mu_{\pm}$ found from the scattering problem (6), and obtains integral equations for $\mu_{\pm}$ as follows

[TABLE]

where $\Delta Q_{\pm}=Q-Q_{\pm}$. In particular, after some calculations, the integrand in Eqs. (34) and (35) can be written as

[TABLE]

and where the subscripts 1 and 2 on $\mu_{\pm}$ identify the matrix columns, i.e., $\mu_{\pm}=(\mu_{\pm,1},\mu_{\pm,2})$. Note that the limits of integration imply that $s=t-\tau$ is either always positive (for $\mu^{-}$) or always negative (for $\mu^{+}$). Recall that \mathop{\rm Im}\nolimits\,\gamma\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\m@th\displaystyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\textstyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\scriptstyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\scriptscriptstyle\hfil#\hfil\cr>\crcr<\crcr} }}}0 respectively in $\Gamma^{+}$ and $\Gamma^{-}$ (27) (and, in contrast, they are not sign-definite in $\mathbb{C}_{\mathrm{I}}$ and $\mathbb{C}_{\mathrm{II}}$).

Requiring boundedness as $\tau\to\pm\infty$, one can then show that the eigenfunctions can be analytically extended in the complex $\lambda$-plane into the following regions:

[TABLE]

The above arguments are made rigorous using Neumann series, as in [83]. The analyticity properties of the columns of the eigenfunction matrices $\phi_{\pm}$ follow trivially from those of $\mu_{\pm}$. Hereafter, we will consistently use the subscripts $\pm$ to denote limiting values as $t\to\pm\infty$, whereas the superscripts $\pm$ will denote analyticity (or more in general meromorphicity) in the regions $\Gamma^{\pm}$. The analyticity region of all eigenfunctions can be seen in Fig. 2(left).

Scattering matrix.

If the wave matrix $\phi(t,z,\zeta)$ solves the scattering problem (6), Abel’s theorem implies

[TABLE]

Also for any spectral-parameter value $\zeta\in\Sigma$, we have $\lim_{t\to\pm\infty}\phi_{\pm}(t,z,\zeta)\,\mathrm{e}^{-i\gamma t\sigma_{3}}=Y_{\pm}$. Hence for any $\zeta\in\Sigma$ we have

[TABLE]

where $\det Y_{\pm}$ was given in Eq. (29). So for any $\lambda\in\Sigma\backslash\{\pm iq_{o}\}$, $\phi_{-}$ and $\phi_{+}$ are two fundamental matrix solutions of the scattering problem. As in the case of NLS equation with NZBC, the behavior at the branch points can be singular. (See for example [82] for a discussion of the possible situations that can arise. A rigorous formulation of the IST for the NLS equation with NZBC considering the singularity at the brunch points is recently presented in [84]. For now we ignore any such situations.) Hence we can write, for any $\zeta\in\Sigma$,

[TABLE]

where $S(\zeta,z)=(s_{i,j})$ is the scattering matrix. The entries $s_{i,j}$ are called the scattering coefficients. It is convenient to introduce the notation  $\phi_{\pm}=(\phi_{\pm,1},\phi_{\pm,2})$, as before. With this notation, the scattering relation (38) becomes the columns system

[TABLE]

Moreover, Eqs. (29) and (37) also imply

[TABLE]

The reflection coefficients that will be needed in the inverse problem are

[TABLE]

Note that $\mathrm{Wr}\phi_{-}=\mathrm{Wr}Y_{\pm}=2\gamma(\gamma-\lambda)/q_{o}^{2}\geqslant 0$ from Eq. (37), where $\mathrm{Wr}$ denotes the Wronskian of a matrix. Using the scattering relation (38), we obtain

[TABLE]

So, from the analyticity properties (36) of the eigenfunctions, one can show that the diagonal entries of the scattering matrix are analytic in the following regions in the complex plane,

[TABLE]

As usual [69], the off-diagonal scattering coefficients are nowhere analytic in general.

We can also write an integral representation of the scattering matrix $S(\zeta,z)$ from Eq. (38), taking the limit of the relation (38) as $t\to-\infty$ and using the integral equations (34) and (35), to get

[TABLE]

Note however that care is necessary to evaluate Eq. (43), since (unlike the case of ZBG) one cannot separate the two terms in parentheses (because the individual terms do not admit separate limits).

3.2 Symmetries and discrete spectrum

Symmetries.

The symmetries of the eigenfunctions and scattering coefficients are complicated by: (i) the presence of the Riemann surface, because, while with ZBG one simply deals with $\lambda\mapsto\lambda^{*}$, here one must also keep track of the sheets of the Riemann surface; (ii) after removing the asymptotic oscillations, the Jost solutions do not tend to the identity matrix as $t\to\pm\infty$.

Recall from Sections 2.2 and 2.3 that  $\gamma_{\mathrm{II}}(\lambda)=-\gamma_{\mathrm{I}}(\lambda)$,  $\zeta=\lambda+\gamma$,  $q_{o}^{2}/\zeta=\gamma-\lambda$,  $\gamma=(\zeta+q_{o}^{2}/\zeta)/2$,  $\lambda=(\zeta-q_{o}^{2}/\zeta)/2$. One can consider several transformations compatible with Eq. (25). As we show next, three of them correspond to actual symmetries of the scattering problem. Of these, two are independent:

  1. $\zeta\mapsto\zeta^{*}$  (upper/lower half plane),  implying  $(\lambda,\gamma)\mapsto(\lambda^{*},\gamma^{*})$  (same sheet);

  2. $\zeta\mapsto-q_{o}^{2}/\zeta$  (outside/inside $C_{o}$),  implying  $(\lambda,\gamma)\mapsto(\lambda,-\gamma)$  (opposite sheets).

In Appendix A.2, we show the following is true:

1. If $\phi(t,z,\zeta)$ is a solution of the scattering problem (6), we have

[TABLE]

where $\dagger$ denotes Hermitian adjoint. Explicitly, the columns satisfy,

[TABLE]

where $\sigma_{*}$ is defined by

[TABLE]

Similarly, for the scattering matrix, we have

[TABLE]

We therefore have the following relations between the scattering coefficients:

[TABLE]

2. If $\phi(t,z,\zeta)$ is a solution of the scattering problem, $\forall\zeta\in\Sigma$ it is

[TABLE]

Recalling Eq. (199), we have, for the columns:

[TABLE]

Again for the scattering matrix we have, $\forall\zeta\in\Sigma$,

[TABLE]

Recalling Eq. (199) we then have, elementwise,

[TABLE]

In Appendix A.2, we also consider the combination of the two symmetries, i.e., 1+2. The result for the scattering matrix is $S^{*}(\zeta^{*},z)=-\sigma_{*}(\sigma_{3}Q_{-})^{-1}\,S(-q_{o}^{2}/\zeta,z)\,\sigma_{3}Q_{+}\sigma_{*}$. Or, elementwise,

[TABLE]

Even though Eqs. (46) and (50) are only valid when $\zeta\in\Sigma$, the scattering coefficients $s_{1,1}$ and $s_{2,2}$ in Eqs. (47), (51), (52), (53) and (54) are analytic in appropriate regions [cf. Eq. (42)]. Hence those three equalities with $s_{1,1}$ and $s_{2,2}$ can be extended uniquely to the appropriate regions shown in Eq. (42) of the $\zeta$-plane via the Schwartz reflection principle.

The above symmetries yield immediately the symmetries for the reflection coefficients:

[TABLE]

Recall that the two reflection coefficients $b$ and $\tilde{b}$ are defined in Eq. (39).

Note that:

(i) Unlike the case of ZBG, where there is only one symmetry, here there are two.

One of them, $\zeta\mapsto\zeta^{*}$, is the same as for ZBG; the other, $\zeta\mapsto-q_{o}^{2}/\zeta$, relates to the definition of $\gamma$.

(ii) Unlike the case of ZBG, and unlike the defocusing NLS with NZBC, here even the symmetries of the non-analytic scattering coefficients involve the map $\zeta\mapsto\zeta^{*}$. This is because here the continuous spectrum is not just a subset of the real $\zeta$-axis.

(iii) The second involution, $\zeta\mapsto-q_{o}^{2}/\zeta$, simply expresses the switch from one sheet to the other. (The corresponding transformation for the defocusing NLS has no minus sign because of the different location of the branch points.) This transformation does not affect $\lambda$. That is, if $f(\lambda)$ is any single-valued function of $\lambda$, it is $f_{\mathrm{I}}(\lambda)=f_{\mathrm{II}}(\lambda)$. Thus, $f$, when expressed as a function of $\zeta$, satisfies the symmetry $f(\zeta)=f(-q_{o}^{2}/\zeta)$. That is because $f$ depends not on $\zeta$ directly, but only through the combination $\lambda=(\zeta-q_{o}^{2}/\zeta)/2$. So, we have $\rho(t,z,\zeta)=\rho(t,z,-q_{o}^{2}/\zeta)$ and $V(t,z,\zeta)=V(t,z,-q_{o}^{2}/\zeta)$. More generally, Eqs. (48), (51) and (52) relate the values of the Jost eigenfunctions and scattering coefficients on opposite sheets of the Riemann surface.

Discrete spectrum and residue conditions.

The discrete spectrum of the scattering problem is the set of all values $\zeta\in\mathbb{C}\setminus\Sigma$ such that decaying eigenfunctions exist. As usual, the discrete spectrum plays an important role in the inverse problem (see Section 5.1). We next show that these values are the zeros of $s_{1,1}(\zeta,z)$ in $\Gamma^{+}$ and those of $s_{2,2}(\zeta,z)$ in $\Gamma^{-}$. Note that, unlike defocusing NLS, we cannot exclude the possible presence of zeros along $\Sigma$. Such zeros would give rise to the so-called “embedded” eigenvalues. In this work we ignore this possibility.

Let $\zeta_{1},\dots,\zeta_{N}$ be the zeros of $s_{1,1}(\zeta,z)$ on the portion of $\Gamma^{+}$ lying in the upper-half plane (UHP). We assume these $N$ zeros are simple. That is, we assume $s_{1,1}(\zeta_{n},t)=0$ and $s^{\prime}_{1,1}(\zeta_{n},t)\neq 0$, with $|\zeta_{n}|>1$ and $\mathop{\rm Im}\nolimits\zeta_{n}>0$ for $n=1,\dots,N$, and where the prime denotes differentiation with respect to $\zeta$.

Recalling the asymptotic behavior of the individual columns of $\phi_{\pm}$ as $t\to\pm\infty$ and the fact that \mathop{\rm Im}\nolimits\gamma(\zeta)\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\m@th\displaystyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\textstyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\scriptstyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\scriptscriptstyle\hfil#\hfil\cr>\crcr<\crcr} }}}0 for $\zeta\in\Gamma^{\pm}$, we have that $\forall\zeta\in\Gamma^{+}$, $\phi_{+,1}(t,z,\zeta)\to 0$ as $t\to\infty$, and $\phi_{-,2}(t,z,\zeta)\to 0$ as $t\to-\infty$. Recalling Eq. (40), however, if $s_{1,1}(\zeta,z)=0$ at $\zeta=\zeta_{n}$ the eigenfunctions $\phi_{+,1}$ and $\phi_{-,2}$ at $\zeta=\zeta_{n}$ must be proportional, i.e.:

[TABLE]

with $b_{n}$ dependent on $z$ in general, but independent of $t$ and $\zeta$. We therefore have an eigenfunction that decays as $t\to\pm\infty$.

Owing to the symmetries (47), (51), (52), (53) and (54), we have that

[TABLE]

For each $n=1,\dots,N$ we therefore have a quartet of discrete eigenvalues. That is, the discrete spectrum is the set

[TABLE]

This is similar to what happens to the focusing and defocusing vector NLS with NZBC and scalar focusing NLS with NZBC. [Instead, for the scalar NLS, in the focusing case (with ZBC) and in the defocusing case (with NZBC) one has symmetric pairs, respectively in the $\lambda$ plane and in the $\zeta$ plane.] The symmetries among a quartet can be seen in Fig. 2(right).

Next we derive the residue conditions. We can write relation (56) equivalently as

[TABLE]

Thus,

[TABLE]

where $C_{n}=b_{n}/s_{1,1}^{\prime}(\zeta_{n},z)$ is called the norming constant. Again, prime denotes differentiation with respect to $\zeta$. Equation (57) is the first of the residue conditions that will be used in the inverse problem. Similarly, from Eq. (41) we have that, if $s_{2,2}(\zeta_{n}^{*},z)=0$ it is

[TABLE]

Or, equivalently, $\mu_{+,2}(t,z,\zeta_{n}^{*})=\tilde{b}_{n}\mathrm{e}^{2i\gamma(\zeta_{n}^{*})t}\mu_{-,1}(t,z,\zeta_{n}^{*})$, and as a result with $\tilde{C}_{n}=\tilde{c}_{n}/s_{2,2}^{\prime}(\zeta_{n}^{*},z)$,

[TABLE]

The above norming constants are related by the symmetries of the problem. Applying relation (44) to Eq. (56) and comparing with Eq. (58) one easily obtains $\tilde{b}_{n}=-b_{n}^{*}$. It is also easy to see that symmetry (47) implies $s_{1,1}^{\prime}(\zeta_{n},z)=(s_{2,2}^{\prime}(\zeta_{n}^{*},z))^{*}$. Hence we have, $\forall n=1,\dots,N$,

[TABLE]

Finally, we need to discuss the remaining two points of the eigenvalue quartet. Applying symmetry (49) to Eqs. (56) and (58) we have the relations

[TABLE]

Moreover, differentiating symmetries (51) and (52), using Eq. (47), and evaluating at $\zeta=\zeta_{n}$ or $\zeta=\zeta_{n}^{*}$, we have

[TABLE]

Combining these relations, we then have

[TABLE]

where for brevity we defined

[TABLE]

Note that $\tilde{C}_{N+n}=-C_{N+n}^{*}$, consistently with symmetry (59).

Hereafter it will be convenient to define $\zeta_{N+n}=-q_{o}^{2}/\zeta_{n}^{*}$ for $n=1,\dots,N$, so that $\zeta_{1},\dots,\zeta_{2N}$ are all the discrete eigenvalues in $\mathbb{C}^{+}$.

3.3 Asymptotics as $\zeta\to 0$ and $\zeta\to\infty$

As usual, the asymptotic properties of the eigenfunctions and the scattering matrix are instrumental in properly normalizing the RHP when we formulate the inverse problem later. Moreover, the next-to-leading-order behavior of the eigenfunctions will allow us to reconstruct the potential from the solution of the RHP. Again, the asymptotics with NZBG is more complicated than in the case of ZBG. The calculations are somewhat simpler and cleaner, however, in the uniformization variable. Note that the limit $\lambda\to\infty$ corresponds to $\zeta\to\infty$ in $\mathbb{C}_{\mathrm{I}}$ and to $\zeta\to 0$ in $\mathbb{C}_{\mathrm{II}}$. It is necessary to study both of these limits to normalize the RHP.

Recall from Eqs. (28) and (30) that the asymptotic eigenvector matrices are $Y_{\pm}=I+(i/\zeta)\sigma_{3}Q_{\pm}$ and $Y_{\pm}^{-1}=(I-(i/\zeta)\sigma_{3}Q_{\pm})/(1+(q_{o}/\zeta)^{2})$. Also recall that $\gamma(\zeta)={\textstyle\frac{1}{2}}(\zeta+q_{o}^{2}/\zeta)$. Consider now the following formal expansion for $\mu_{-}(t,z,\zeta)$:

[TABLE]

with

[TABLE]

Let $A_{d}$ and $A_{o}$ denote respectively the diagonal and off-diagonal parts of a matrix $A$.

In Appendix A.3, we prove that Eq. (63) provides an asymptotic expansion for the columns of $\mu_{-}(t,z,\zeta)$ as $\zeta\to\infty$ in the appropriate region of the $\zeta$-plane described by Eq. (36), with $\forall m\in\mathbb{N}$ and

[TABLE]

Explicitly, the result holds with $\mathop{\rm Im}\nolimits\zeta\leqslant 0$ for the first column and $\mathop{\rm Im}\nolimits\zeta\geqslant 0$ for the second column. A similar asymptotic behavior to Eqs. (66) and (67) holds for $\mu^{+}$, and can be derived in a similar way.

Next we consider the asymptotics as $\zeta\to 0$. We use the same formal expansion (63). In Appendix A.3, we prove that Eq. (63) provides an asymptotic expansion for the columns of $\mu_{-}(t,z,\zeta)$ as $\zeta\to 0$ in their analytic regions of the $\zeta$-plane, with $\forall m\in\mathbb{N}$,

[TABLE]

Again, similar asymptotic behavior to Eqs. (68) and (69) holds for $\mu_{+}$, and can be derived in a similar way.

In particular, by computing explicitly the first five terms in expansion (63) we have that, as $\zeta\to\infty$,

[TABLE]

Equation (70) will provide the means to reconstruct the scattering potential $Q(t,z)$ from the solution of the RHP in the inverse problem.

Finally, inserting the above asymptotic expansions for the Jost eigenfunctions into the Wronskian representations (40) and (41) one shows that, as $\zeta\to\infty$ in the appropriate regions of the complex $\zeta$-plane,

[TABLE]

Explicitly, the above estimate holds with $\mathop{\rm Im}\nolimits\zeta\geqslant 0$ and $\mathop{\rm Im}\nolimits\zeta\leqslant 0$ for $s_{1,1}$ and $s_{2,2}$, respectively, and with $\mathop{\rm Im}\nolimits\zeta=0$ for $s_{1,2}$ and $s_{2,1}$. Similarly, one shows that, as $\zeta\to 0$,

[TABLE]

again in regions of the $\zeta$-plane from Eq. (42).

3.4 Asymptotics of the density matrix

One can see (by direct substitution) that, as in the case of ZBG, if $\phi(t,z,\zeta)$ is any fundamental matrix solution of the scattering problem, it is

[TABLE]

Hence, the quantity in square brackets must be independent of $t$. So, $\forall\zeta\in\Sigma$, we can define

[TABLE]

Conversely,

[TABLE]

Then, taking the limit as $t\to\pm\infty$ and using Eqs. (34) and (35), we obtain

[TABLE]

However, note that the density matrix $\rho(t,z,\zeta)$ by itself does not have a limit. In other words, one should remember that despite the $\pm$ signs on the left-hand side of Eq. (72), $\rho_{\pm}(\zeta,z)$ are not simply the limits of $\rho(t,z,\zeta)$ as $t\to\pm\infty$. Note also that, even though $\rho$ is a single-valued function of $\lambda$, the Jost solutions $\phi_{-}$ and $\phi_{+}$ are not. Thus, like $\phi_{\pm}$, in general $\rho_{\pm}$ is only single-valued on $\mathbb{R}$, unlike in the case of ZBG. (In other words, $\rho_{\pm,{\mathrm{II}}}\neq\rho_{\pm,{\mathrm{I}}}$ in general.)

One can also see by direct calculation that  $S\,\rho_{+}=\rho_{-}S=\phi_{-}^{-1}\rho\,\phi_{+}$. Hence,

[TABLE]

Equation (75) relates the asymptotic values of the density matrix as $t\to\pm\infty$, and allows one to obtain $\rho_{+}$ from the knowledge of $\rho_{-}$ and $S$ [which in turn is completely determined by $q(t,z)$]. Thus, one can only pick one of $\rho_{\pm}$, after which the other one is fixed. Due to the causality, we choose $\rho_{-}$.

Note that properties of the density matrix ${\mathrm{tr}}\rho=0$ and $\det\rho=-1$ imply that ${\mathrm{tr}}\rho_{\pm}=0$ and $\det\rho_{\pm}=-1$. Note also that since $\rho^{\dagger}=\rho$ if $\zeta\in\mathbb{R}$, it is also $\rho_{\pm}^{\dagger}=\rho_{\pm}$ if $\zeta\in\mathbb{R}$. In the rest of this work, we assume that $\rho^{\dagger}=\rho(\zeta^{*})$ when $\zeta\in C_{o}$, which implies $\rho_{\pm}^{\dagger}=\rho_{\pm}(\zeta^{*})$ for all $\zeta\in\Sigma$. [Notice that this assumption will not affect our final solution $\rho$ since it is defined only on the real axis.] Thus, we denote the entries of $\rho_{\pm}$ as

[TABLE]

with $P_{\pm}(\zeta)=P_{\pm}(\zeta^{*})$. Note however that $D_{\pm}$ and $P_{\pm}$ are not the limits of $D$ and $P$ as $t\to\pm\infty$. (Indeed, as we show below, such limits do not exist in general.) Combining the symmetry (48) and definition (72), we have

[TABLE]

This implies that one does not have freedom to pick the asymptotic state $\rho_{\pm}$ for all values of $\lambda$ (or $\zeta$). One can only pick $\rho_{\pm}$ for $\lambda\in\mathbb{R}$ on the first sheet [or $\zeta\in(-\infty,-q_{o}]\cup[q_{o},\infty)$]. Thus, let us denote $\rho_{\pm,{\mathrm{I}}}$ and $\rho_{\pm,{\mathrm{II}}}$ to the values of the asymptotic state on the first or second sheet, respectively. Then, we define the asymptotic state on the two $\lambda$-sheets as

[TABLE]

or, in the $\zeta$-plane,

[TABLE]

Equation (75) can then be used to obtain $D_{+}$ and $P_{+}$ from $D_{-}$ and $P_{-}$. Note that, in principle, $D_{\pm}$ and $P_{\pm}$ depend on $\lambda$ on each sheet. Here however we make the assumption that the preparation of the medium on each $\lambda$-sheet is independent of $\lambda$, which is the simplest possible situation from a physical point of view.

Because the asymptotic state (77), (78), (79) and (80) satisfies the symmetry (76), we know that the calculations performed on both $\lambda$-sheets are equivalent. Hereafter, we will focus on the first $\lambda$-sheet.

Equation (74) gives the explicit relation between $\rho_{\pm}$ and $\rho$. In component form, it is

[TABLE]

Conversely, Eq (76) also implies that $\rho(t,z,\zeta)=Y_{\pm}\mathrm{e}^{i\gamma t\sigma_{3}}\rho_{\pm}(\zeta,z)\mathrm{e}^{-i\gamma t\sigma_{3}}Y_{\pm}^{-1}+o(1)$ as $t\to\pm\infty$, which in component form, as $t\to\pm\infty$, is,

[TABLE]

Let us now discuss the implications of the relations (81), (82), (83) and (84). First, just as in the case of ZBG, Eq. (84) implies that, in general, $P$ does not have a limit as $t\to\pm\infty$, but instead oscillates in time. Unlike in the case of ZBG, however, Eq. (83) implies that, in general, due to the nonzero background radiation $q_{\pm}$, the quantity $D$ also does not have a limit in time and instead oscillates. Second, unlike in the case of ZBG, $P\neq 0$ as $t\to\pm\infty$ even in the particular case in which $D$ and $P$ are time-independent and do tend to a limit as $t\to\infty$ or $t\to-\infty$, which is when $P_{+}=0$ or $P_{-}=0$, respectively. The nonzero contribution arises from the polarization induced by the limiting electric field. (If one takes the limit $q_{o}\to 0$, one of course recovers the relation $\lim_{t\to\pm\infty}D=D_{\pm}$ that holds in the case of ZBG.)

To further elaborate on the above, let us focus on the asymptotic state at $t\to-\infty$. Recall that the normalization $\det\rho_{-}=-1$ implies the constraint

[TABLE]

which implies that one does not have the freedom to assign $D_{-}$ and $P_{-}$ independently. In particular, in the special case when $\rho_{-}$ is diagonal, i.e., $D_{-}=\pm 1$ and $P_{-}=0$, we have

[TABLE]

which yields $P=-iq_{-}D/\gamma+o(1)$ as $t\to-\infty$. Therefore, the polarization in this case is entirely induced by the limiting values of the electric field. On the other hand, the time-dependent terms in Eq. (84) describe the intrinsic part of the polarization, which is due to the excitation of the medium as $t\to-\infty$, and has parts proportional to $P_{-}$. Note that, generally, the induced polarization is still present in the initial state for the density matrix $\rho$ even in the presence of nonzero intrinsic polarization, except that it is reduced in magnitude due to the constraint (85) together with $D_{-}$.

Finally, we remark that, even in the absence of intrinsic initial polarization, the medium is not in a pure state, because of the initial induced polarization, which must exist (and so $\lim_{t\to-\infty}|P|\neq 0$). Instead, initially the medium is necessarily in a superposition of two states. One can also obtain this result directly from $\lim_{t\to-\infty}|D|=|\lambda/\gamma|<1$ for any $q_{o}>0$.

The initial amount of population inversion.

The above discussion shows that although the quantity $D_{-}$ is not the direct limit of $D$ even in an initial state with no intrinsic polarization, it still provides an important indicator of the initial state of the medium. In fact, $D_{-}$* is directly proportional to the amount of the initial population inversion in the medium*, i.e., it describes the amount of the population inversion allowed initially by the background electric field. More generally, $D_{-}$* is proportional to the time-averaged amount of the initial population inversion in the medium*, i.e., it describes the amount of the time-averaged population inversion allowed initially by the background electric field and the intrinsic polarization. Recall that the population inversion is defined as the difference between the populations of the atoms in the excited state and in the ground state. Thus, if $-1\leqslant D_{-}<0$, initially there are more atoms in the ground state than in the excited state on average, while if $0<D_{-}\leqslant 1$, initially there are more atoms in the excited state than in the ground state on average. As we show in Section 6.6, the former condition leads to stable and the latter to unstable electric-field propagation through the medium. This is consistent with what happens to the case of ZBG ($D_{\pm}=\lim_{t\to\pm\infty}D$).

We further will see in Section 5.2 that pure soliton solutions can only exist when the initially intrinsic polarization of the medium vanishes. In the case when the ground state is more occupied initially than the excited state these solitons are subluminal and stable, and thus physical, while in the opposite case they are superluminal and unstable, and thus unphysical.

4 Propagation

Recall that in the MBEs, the traditional role of evolution variable is played by the the physical propagation distance $z$. As a result, we will refer to the $z$-dependence as the propagation.

4.1 Evolution of the background

Let us first discuss the $z$ dependence of the asymptotic values of the optical field. The propagation of these values is given by the limit as $t\to\pm\infty$ of Eq. (4). That is,

[TABLE]

Using Eq. (72), we can express the solution of Eq. (3) as

[TABLE]

Note that the off-diagonal elements of $\rho$ can never be zero as $t\to\pm\infty$, even when the off-diagonal elements of $\rho_{\pm}$ are identically zero. This is because of the background-induced polarization from Eqs. (83) and (84). Denote the diagonal and off-diagonal part of a matrix the subscripts $d$ and $o$, respectively. In particular, let $\rho_{\pm}=\rho_{\pm,d}+\rho_{\pm,o}$. Note that $[\sigma_{3},\rho_{\pm,d}]=0$ and $[\sigma_{3},\rho_{\pm,o}]=2\sigma_{3}\rho_{\pm,o}$. (Also, $Q_{\pm}\sigma_{3}+\sigma_{3}Q_{\pm}=0$.) Using definitions (87) and (28) we then have

[TABLE]

All entries in the second term on the RHS, however, contain oscillating exponentials, which, when integrated, tend to zero in the limit $t\to\pm\infty$ thanks to the Riemann-Lebesgue lemma. Therefore, the effective part of $[\sigma_{3},\rho]$ as $t\to\pm\infty$ is

[TABLE]

Notice that the first and last commutators on the RHS vanish, because $\sigma_{3}$, $\rho_{\pm,d}$ and $\sigma_{3}Q_{\pm}\rho_{\pm,d}\sigma_{3}Q_{\pm}$ are all diagonal. Finally, note that $[\sigma_{3},[\sigma_{3}Q_{\pm},\rho_{\pm,d}]]=2[Q_{\pm},\rho_{\pm,d}]$, and recall $\det Y_{\pm}=2\gamma/\zeta$.

Evaluating Eq. (4) in the limit $t\to\pm\infty$ (assuming the derivative with respect to $z$ and the limits as $t\to\pm\infty$ commute) we then obtain the propagation of $Q_{\pm}$ with respect to $z$:

[TABLE]

where

[TABLE]

Recall that ${\mathrm{tr}}\,\rho_{\pm}={\mathrm{tr}}\,\rho=0$ and $D_{\pm}$ has opposite signs on sheet I and sheet II from symmetries (77), (78), (79) and (80). Since $\gamma$ also takes opposite signs on both sheets, it is easy to verify that $w_{\pm}$ is single-valued, as it should be. Equation (88) can be integrated to obtain

[TABLE]

or simply,

[TABLE]

Note that $|q_{\pm}(z)|^{2}=|q_{\pm}(0)|^{2}=q_{o}^{2}$, $\forall z\in\mathbb{R}^{+}$. In particular, if $w_{-}$ is independent of $z$, taking $q_{-}(0)=q_{o}$ without loss of generality we simply have

[TABLE]

Note also that Eq. (90) is similar to Eq. (23). The difference is that Eq. (90) describes the background as $t\to-\infty$, whereas Eq. (23) describe the entire background solution for $\forall t\in\mathbb{R}$. The corresponding propagation of the asymptotic eigenvector matrices is given by

[TABLE]

and the asymptotic behavior of the Jost solutions is

[TABLE]

Using relation (75), one can then express $w_{+}$ in terms of $\rho_{-,d}$. In general, $w_{+}\neq w_{-}$, and therefore $q_{+}$ evolves at a different rate than $q_{-}$. However,

[TABLE]

which shows $w_{+}=w_{-}$ for reflectionless solutions (i.e., when $S_{o}=0$) or when $\rho_{-,o}=0$. (We will see later that indeed $\rho_{-,o}=0$ is necessary for radiation to remain zero as a function of $z$.)

4.2 Evolution of reflection coefficients and norming constants

We now discuss the propagation of the reflection coefficients and norming constants. As in the case with ZBG, this is the place where the theory deviates the most from the IST for AKNS systems such as the focusing NLS equation with NZBC [69]. Here the situation is further complicated by the nontrivial behavior of the Jost solutions as $t\to\pm\infty$.

Simultaneous solutions of the Lax pair and auxiliary matrix.

Since the asymptotic behavior of $\phi_{\pm}(t,z,z)$ as $t\to\pm\infty$ is fixed, in general they will not be solutions of Eq. (7). Because both $\phi_{+}$ and $\phi_{-}$ are fundamental matrix solutions of the scattering problem, any other solution $\Phi(t,z,\zeta)$ can be written as

[TABLE]

where $C_{\pm}(\zeta,z)$ are $2\times 2$ matrices independent of $t$. Suppose that $\Phi(t,z,\zeta)$ is a simultaneous solution of both parts of the Lax pair (6) and (7), so $\Phi_{z}=V\,\Phi$. Then some algebra shows that,

[TABLE]

where the auxiliary matrices $R_{\pm}$ are given by

[TABLE]

It is not obvious a priori that the RHS of Eq. (95) is independent of $t$, but Eq. (94) shows it must be. Also, even though $g(\lambda)$ and $\rho(t,z,\lambda)$ are only defined for $\lambda\in\mathbb{R}$, definition (95) can be evaluated $\forall\zeta\in\Sigma$. But since $\phi_{\pm}$ are not single-valued functions of $\lambda$, neither are $C_{\pm}$ and $R_{\pm}$. That is, in general we have $C_{\pm,{\mathrm{II}}}\neq C_{\pm,{\mathrm{I}}}$ and $R_{\pm,{\mathrm{II}}}\neq R_{\pm,{\mathrm{I}}}$.

Moreover, recalling definition (72), we can write $V(t,z,\zeta)$ as

[TABLE]

with shorthand notation $\zeta^{\prime}=\zeta(\lambda^{\prime})$. As we show below, Eqs. (95) and (96) allow one to compute the propagation of the reflection coefficients and norming constants, as in the case with ZBG. Note that, even though the individual terms on the RHS of Eq. (96) are only single-valued on the real line, the whole RHS is a single-valued function for $\lambda\in\mathbb{C}$. In fact, one can use the symmetries of $\phi_{\pm}$ and $\rho_{\pm}$ to verify that $V(t,z,\zeta)=V(t,z,-q_{o}^{2}/\zeta)$.

We next express $R_{\pm}$ in terms of known quantities for all $z\geqslant 0$. Recalling the definition (95), and under the assumption that $z$-derivatives and the limits as $t\to\pm\infty$ commute, we have

[TABLE]

with $\gamma^{\prime}=\gamma(\lambda^{\prime})$. Above and throughout this subsection, the prime will denote the functional argument in the integrand of the Hilbert transform, not derivative with respect to $\zeta$.

Note that: (i) unlike in the case of ZBG, the last term in the curly bracket does not vanish; (ii) one cannot evaluate the terms in the curly brackets individually. The detailed calculation is presented in Appendix A.4, where we show that $R_{\pm}=R_{\pm,d}+R_{\pm,o}$, with

[TABLE]

where $\nu=1$ for $\zeta\in(-\infty,-q_{o}]\cup[q_{o},\infty)$ and $\nu=-1$ for $\zeta\in(-q_{o},q_{o})$, and where the subscripts $d$ and $o$ denote diagonal and off-diagonal parts, as before. Equations (98) and (99) in component form is

[TABLE]

where $R_{\pm,i,j}$ the $(i,j)$-elements of $R_{\pm}$. One should notice that in Eq. (98), both $\rho_{\pm,d}^{\prime}$ and $\gamma^{\prime}$ take opposite signs on sheet I and II, so the Hilbert transform achieves a unique result on the two $\lambda^{\prime}$-sheets, as it should. Therefore, we conclude that the matrix $R_{\pm,d}$ is determined independently of the choice of the integration variable.

Note that, even though the calculation is considerably more involved than that in the case of ZBG, the final result appears relatively simple. As in the case of ZBG, the matrices (98) and (99) enable one to calculate the propagation equations for the reflection coefficients and norming constants. Importantly, however, note that, since $R_{\pm}(\zeta,z)$ are defined in terms of a principal value integral, even when it admits an extension to the complex $\lambda$-plane, their values are discontinuous across the real $\lambda$-axis, both from above and from below.

In general, it is not possible to extend all the entries of $R_{\pm}$ off the continuous spectrum. Some of the entries, however, can be extended, as we show next. Recalling Eq. (97) and following the calculations in Appendix A.4, the off-diagonal entries of $R_{\pm}$ can be written as

[TABLE]

We point out that Eq. (102) is equivalent to Eq. (99) via Lemma A.1. For each matrix element, the Hilbert transform on the RHS of Eq. (102) is analytic and bounded in the complex $\zeta$-plane wherever the exponential inside tends to zero as $t\to\pm\infty$. Hence, looking at the regions where \mathop{\rm Im}\nolimits\gamma\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\m@th\displaystyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\textstyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\scriptstyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\scriptscriptstyle\hfil#\hfil\cr>\crcr<\crcr} }}}0, we have:

[TABLE]

Below we will show that, just as in the case of ZBG, these conditions are all that is necessary to determine the propagation equations for the reflection coefficient and norming constants.

Evolution equations for the reflection coefficient.

Since $\phi_{+}=\phi_{-}\,S$, Eq. (93) implies,

[TABLE]

Some algebra then shows that,

[TABLE]

Recall that the reflection coefficients are $b=s_{2,1}/s_{1,1}$ and $\tilde{b}=s_{1,2}/s_{2,2}$, with $\tilde{b}(\zeta,z)=-b^{*}(\zeta,z)$. Let us introduce the matrix

[TABLE]

where once more the subscripts $d$ and $o$ denote the diagonal and off-diagonal parts, respectively. In Appendix A.5 we show that

[TABLE]

where $\nu=1$ for $\zeta\in(-\infty,-q_{o}]\cup[q_{o},\infty)$ and $\nu=-1$ for $\zeta\in(-q_{o},q_{o})$, similarly to Eqs. (98) and (99). (Note that this result is formally identical to that of the case with ZBG on sheet I.) In particular,

[TABLE]

where $\nu=1$ if $\zeta\in(-\infty,-q_{o}]\cup[q_{o},\infty)$, $\nu=0$ if $\zeta\in C_{o}$, and $\nu=-1$ if $\zeta\in(-q_{o},q_{o})$, and

[TABLE]

Note that the propagation equation for the reflection coefficient is completely determined by the asymptotic values as $t\to-\infty$, as it should be. Notice also that because of the symmetry of the reflection coefficient (55), in general, $b(\zeta,z)$ with $\zeta\in(-q_{o},q_{o})$ evolves differently from $b(\zeta,z)$ with $\zeta\in(-\infty,-q_{o}]\cup[q_{o},\infty)$. This can also be seen from Eq. (107) with $\nu=\pm 1$, and from the fact that $D_{-}$ and $P_{-}$ take different values depending on $\zeta$. However, everything is still consistent, as $D_{-}$ and $P_{-}$ also satisfy the symmetries (77), (78), (79) and (80). Thus, again, we only focus on the first $\lambda$-sheet, i.e., $\zeta\in(-\infty,-q_{o}]\cup[q_{o},\infty)$.

Evolution equations for norming constants.

Recall that $C_{n}$ is defined by Eq. (57):

[TABLE]

In Appendix A.5 we show that for $n=1,\dots,N$,

[TABLE]

where $R_{-,1,1}$ is given in Eq. (100). Hereafter, one can use Eqs. (59) and (62) to obtain the propagation for $C_{n}$ with $n=N+1,\dots,2N$. Like the propagation of the reflection coefficient, the propagation equation for the norming constant is completely determined by the asymptotic values as $t\to-\infty$, as it should be. Also, in the ODE (110), the quantity $R_{-,1,1}$ is evaluated at $\zeta=\zeta_{n}\in\Gamma^{+}$, i.e., in the first $\lambda$-sheet. Therefore, norming constants $C_{n}$ are uniquely determined without any ambiguities.

5 Inverse problem

5.1 Riemann-Hilbert problem and reconstruction formula

The development of the inverse problem begins from the scattering relation (38). One can rewrite this scattering relation as $(\mu_{+,1},\mu_{+,2})=(\mu_{-,1},\mu_{-,2})\,\mathrm{e}^{i\gamma t\sigma_{3}}S\,\mathrm{e}^{-i\gamma t\sigma_{3}}$. Equivalently, one can write

[TABLE]

where the $(t,z,\zeta)$-dependence is omitted for brevity and the reflection coefficients $b$ and $\tilde{b}$ are defined in Eq. (39). One can view Eq. (111) as a relation between eigenfunctions analytic in $\Gamma^{+}$ and those in $\Gamma^{-}$. Thus, we introduce the meromorphic matrices

[TABLE]

As mentioned in Section 3.1, subscripts $\pm$ indicate normalization as $t\to\pm\infty$, while superscripts $\pm$ distinguish between analyticity in $\Gamma^{+}$ and $\Gamma^{-}$, respectively. We then have,

[TABLE]

where the jump matrix $G$ is

[TABLE]

Equation (113) is the jump condition for the desired RHP. This RHP is solved explicitly in Appendix A.6, where it is shown that the solution is

[TABLE]

with $\zeta_{N+n}=-q_{o}^{2}/\zeta_{n}^{*}$ for $n=1,\dots,N$. Note that the expressions for $M^{+}$ and $M^{-}$ are formally identical, except for the fact that $M$ is evaluated in different regions of the complex plane.

With the relabeling in definition (62), the residue relations (57), (60) and (61) imply that only the first column of $M^{+}$ has a pole at $\zeta=\zeta_{n}$ with $n=1,\dots,2N$, and its residue is proportional to the second column of $M^{+}$ at that point, explicitly:

[TABLE]

Similarly,

[TABLE]

We can therefore evaluate the second column of the solution (115) at $\zeta=\zeta_{n}$, obtaining:

[TABLE]

for $n=1,\dots,2N$, and where $G_{2}$ denotes the second column of $G$. Similarly, we can evaluate the first column of the solution (115) at $\zeta=\zeta_{n}^{*}$, obtaining:

[TABLE]

again for $n=1,\dots,2N$, and where $G_{1}$ denotes its first column. Finally, evaluating the solution (115) with $\zeta\in\Sigma$, we obtain together with equations (118) and (119), a closed linear system of algebraic-integral equations for the solution of the RHP.

The last remaining task in the inverse problem is to reconstruct the potential and the density matrix from the solution (115) of the RHP. It is easy to compute the asymptotic behavior of $M(t,z,\zeta)$

[TABLE]

where the residues are given by Eqs. (116) and (117). Taking $M=M^{+}$ and comparing the $(1,2)$-element of Eq. (120) to the one of Eq. (70), we then obtain the reconstruction formula for the potential:

[TABLE]

where $\mu_{-,2}(t,z,\zeta_{n}^{*})$ and $M^{+}(t,z,\zeta)$ are obtained from the simultaneous solution (115) for $\zeta\in\Sigma$ and Eqs. (118) and (119).

To reconstruct the density matrix, one also starts from the solution (115) and Eqs. (118) and (119), then computes $M^{-}(t,z,\zeta)$ for $\zeta\in\Sigma$. Recalling definition of $M^{\pm}$ (112), we then have

[TABLE]

In turn, from the definition (33) of $\mu$ and relation (73), we obtain the density matrix as

[TABLE]

Similarly to the NLS equation, it is also possible to derive a “trace” formula, which allows one to recover the analytic scattering coefficient from the knowledge of the reflection coefficient and discrete eigenvalues, and also a so-called “theta” condition, which yields the asymptotic phase difference of the potential in terms of the discrete eigenvalues:

[TABLE]

as well as

[TABLE]

See Appendix A.7 for details. In fact, the trace formulae (123) and (124) and the “theta” condition (125) are exactly the same as the ones obtained in focusing NLS with NZBC [69].

5.2 Reflectionless potentials

We now look at a special case of potentials $Q(t,z)$ for which the reflection coefficient $b(\zeta,z)$ vanishes identically. Note that in order for the solution to remain reflectionless for $z>0$ it is not sufficient that $b(\zeta,0)=0$. Indeed, the propagation equation (107) for the reflection coefficient $b(\zeta,z)$ shows that $b(\zeta,z)=0$ for $z>0$ only if $\rho_{-,o}(\zeta,z)=0$ $\forall z>0$. We therefore assume that this is the case for the rest of this section. Recalling the general asymptotic state (77), (78), (79) and (80) we then define the initial state for the medium as follows

[TABLE]

where $\nu=1$ for $\zeta\in(-\infty,-q_{o}]\cup[q_{o},\infty)$, $\nu=-1$ for $\zeta\in(-q_{o},q_{o})$, and where $h_{-}=\pm 1$ is related to the initial population inversion via Eq. (86). Notice that the sign discontinuities introduced by $\nu$ are necessary, because of the symmetry (76). Also recall that $h_{-}=\pm 1$ indicates whether, initially, more atoms are in the excited state than in the ground state or vice versa, respectively (cf. the discussion at the end of Section 3.4). Again, we emphasize the fact, shown in Section 3.4, that, due to the presence of background radiation, the medium is always polarized, and as a result, pure states do not exist in general, which differs from the case of ZBG. Moreover, without calculating the density matrix $\rho$ from the inverse problem explicitly, we know that $D\to\lambda h_{-}/\gamma$ and $P\to-ih_{-}q_{-}/\gamma$ as $t\to-\infty$ from Eqs. (83) and (84). In other words, the initial state of the medium is independent of discrete eigenvalues and soliton types (discussed later), but is dependent on $\lambda$.

As usual, for reflectionless potentials there is no jump from $M^{+}$ to $M^{-}$ across the continuous spectrum, and the RHP therefore reduces to an algebraic system, whose solution yields the soliton solutions of the nonlinear system.

It is convenient to introduce scalar functions

[TABLE]

Also recall the symmetry (62), and $\gamma(\zeta^{*})=\gamma^{*}(\zeta)$. The algebraic system obtained from the RHP can then be expressed as

[TABLE]

Substituting Eq. (127) into Eq. (128) yields, for all $n=1,\dots,2N$,

[TABLE]

Note that only the first component of the above eigenfunctions is needed in the reconstruction formula (121) for the potential. Let us write the resulting system in matrix form. Let $\mathbf{X}=(X_{1},\dots,X_{2N})^{T}$ and $\mathbf{B}=(B_{1},\dots,B_{2N})^{T}$, where

[TABLE]

Also, let us define the $2N\times 2N$ matrix $A=(A_{n,l})$, where

[TABLE]

The system (129) then becomes simply $K\,\mathbf{X}=\mathbf{B}$, where

[TABLE]

The solution of the system is $X_{n}=\det\,K_{n}/\det K$ for $n=1,\dots,2N$, where the $2N\times 2N$ matrices $K_{1},\dots,K_{2N}$ are $K_{n}=(\mathbf{K}_{1},\dots,\mathbf{K}_{n-1},\mathbf{B},\mathbf{K}_{n+1},\dots,\mathbf{K}_{2N})$. Finally, upon substituting $X_{1},\dots,X_{2N}$ into the reconstruction formula (121), one can write the result for the potential compactly as

[TABLE]

with $\mathbf{D}=(D_{1},\dots,D_{2N})^{T}$ and $D_{n}=\tilde{C}_{n}(z)\,\mathrm{e}^{2i\gamma(\zeta_{n}^{*})t}$ for $n=1,\dots,2N$.

The reconstruction formula for the density matrix also takes on a simpler form in the reflectionless case. In this case we also need the second component of the eigenfunctions in the linear system (127) and (128). But it is trivial to see that, from Eq. (129), one can obtain them from the solution of

[TABLE]

where $\mathbf{Y}=(Y_{1},\dots,Y_{2N})^{T}$ are the unknowns and $\mathbf{L}=(L_{1},\dots,L_{2N})^{T}$, with

[TABLE]

and where the coefficient matrix $K$ is still given by definition (130). Once $\mathbf{X}$ and $\mathbf{Y}$ have been obtained, from Eqs. (127) and (128) we have

[TABLE]

which are related by Eq. (127), i.e.,

[TABLE]

Recalling that in the reflectionless case there is no jump, we then have $\mu_{-}(t,z,\zeta)=\left(M_{-,1},M_{+,2}\right)$ $\forall\zeta\in\Sigma$, and therefore, from Eq. (115) we find

[TABLE]

Recall the relationship between the density matrix $\rho$ and the eigenfunction matrix $\mu$ is given by relation (122). In the reflectionless situation, it simplifies to

[TABLE]

6 Exact solutions and their behavior

In this section we use the IST formalism developed in previous sections to obtain a variety of exact solutions of the MBEs with NZBG. In particular, we will focus on variation of one-soliton solutions with a single discrete eigenvalue $\zeta_{1}$. According to the location of this discrete eigenvalue, there exist four types of solutions, shown in Fig. 3. We will discuss all types separately in detail. In addition, we will also discuss the stability of soliton solutions of the MBEs.

6.1 One-soliton solutions: representation formulae

Without loss of generality, we take the input background to be $q_{-}(0)=q_{o}$ (owing to the phase invariance of the MBEs), so that the background propagates as $q_{-}(z)=q_{o}\mathrm{e}^{2iW_{-}(z)}$ with $W_{-}(z)\in\mathbb{R}$ defined in Eq. (90). We parametrize the discrete eigenvalue $\zeta_{1}$ and the norming constant $C_{1}$ respectively as

[TABLE]

with $\eta>1$, $\alpha\in(0,\pi)$ and both $\xi(z)$ and $\varphi(z)$ real. For later use, we also defined quantities

[TABLE]

Recall that, for all the soliton solutions, the initial state for the medium is given by Eq. (126) with $h_{-}=\pm 1$. We further assume $h_{-}$ is independent of $z$ and $\zeta$. We take $\rho_{-}$ to be diagonal because only in this case a reflectionless solution remains reflectionless upon propagation. The formula for optical field will be given explicitly below. However, due to complexity of the expressions, we will only give explicit formulas for the modified eigenfunction matrix $\mu$ instead of the density matrix $\rho$. One can then simply use the relationship (133) to reconstruct $\rho(t,z,\zeta)$. Moreover, it turns out that $\mu$ satisfies the symmetries $\mu_{2,1}=-\mu_{1,2}^{*}$ and $\mu_{2,2}=\mu_{1,1}^{*}$, so it is sufficient to give only the $(1,1)$ and $(1,2)$ components of $\mu$.

Type 1. Traveling-wave soliton: solutions with purely imaginary eigenvalues.

Let us start from the simplest case $\alpha=\pi/2$ in Eq. (134). The theta condition (125) implies that the corresponding asymptotic phase difference is $2\pi$ (i.e., no phase difference). The general soliton solution (131) and (132) reduces to

[TABLE]

with

[TABLE]

Recall $W_{-}(z)$ was defined in Eq. (90) and all quantities in Eq. (135). Note that $\Delta_{+}>2$, so the above soliton solution is always non-singular. Equation (136) is the analogue of the Kuznetsov-Ma solution of the focusing NLS equation with NZBC [85, 86]. Indeed, the calculations to obtain the soliton solution (136) are similar to those for the focusing NLS equation (e.g., see [69] for details), but the resulting solutions exhibit quite different behavior upon propagation.

Importantly, the $z$ dependence of the solution is entirely determined by the norming constant via $\xi(z)$, $\varphi(z)$ and the background $q_{-}(z)$. We will discuss the $z$ dependence later.

Type 2. Oscillatory soliton: solutions with eigenvalues in general position.

A more general expression of the one-soliton solution for a general discrete eigenvalue can also be found, similarly to that for the focusing NLS equation. (Again, see [69] for details.) Recalling the parametrization (134) for the discrete eigenvalue and norming constant, after tedious but straightforward calculations, we obtain

[TABLE]

where

[TABLE]

and $\Delta_{\pm}$, $d_{\pm}$ are defined in Eq. (136) and $W_{-}(z)$ is defined in Eq. (90). Note that the term $\cosh(\chi-2i\alpha)$ can be easily expressed in terms of real-valued functions using addition formulae. Equation (141) is the analogue of the Tajiri-Watanabe solutions of the focusing NLS equation with NZBC [87], but exhibits different spatial-temporal behavior.

The same as for the traveling-wave soliton solutions, the propagation of the above solution is also determined by the three quantities $\xi(z)$, $\varphi(z)$ and $q_{-}(z)$ from Section 4, which will be discussed later. The names “traveling-wave soliton” and “oscillatory soliton” will become clear after we compute the $z$ propagation explicitly.

Soliton amplitude.

We define the instantaneous amplitude compared to the background as

[TABLE]

It is easy to see that the soliton achieves its maximum when $\chi=0$ in both cases (136) and (141). For the oscillatory soliton solution (141), the explicit formula (148) is complicated and does not simplify in general, so it is omitted for brevity. However, for the traveling-wave soliton (136), simple calculations show that the amplitude (148) becomes

[TABLE]

The maximal possible amplitude is attained when $s(z)=(2n+3/2)\pi$ with $n\in\mathbb{Z}$, and is

[TABLE]

Notice that $A_{\mathrm{max}}$ is an increasing function of $\eta$. That is, traveling-wave soliton solutions have a larger maximal possible amplitude when the discrete eigenvalues move farther away from the continuous spectrum.

Importantly, however, note that, for some solutions, the above theoretical maximum may not be achieved in practice. This is because, as we will show in Section 6.2, for discrete eigenvalues on the imaginary axis, $s(z)$ will turn out to be constant, and therefore the condition $s(z)=(2n+3/2)\pi$ may not be achieved for any value of $z$ or $n$. In such cases, the actual maximum of the solution depends on the phase of the norming constant. In fact, the entire shape of the traveling-wave soliton solution depends on the phase of the norming constant. This is unlike the case of ZBG and also unlike the solutions of the focusing NLS equation.

Soliton velocity.

By inspecting soliton solutions (136) or (141), it is evident that solitons are localized along the curve $\chi(t,z)=y$. Our assumption that the initial state for the medium (126) is independent of $z$ implies that $\xi(z)$ from Eqs. (136) and (141) has the form $\xi(z)=\xi_{1}\,z+\xi(0)$, where $\xi_{1}$, $\xi(0)\in\mathbb{R}$ and are independent of $z$. Then, provided that $\xi_{1}\neq 0$, the soliton is located along the line $z-Vt=z_{o}$, and travels with velocity $V$ in the $(t,z)$ frame, where

[TABLE]

Now recall that the MBEs (1) and (2) are written in a comoving frame of reference. Straightforward algebra shows that the physical soliton velocity in the laboratory frame is,

[TABLE]

where $c$ is the speed of light in vacuum. That is, in laboratory coordinates, the soliton is located along the line $z-V_{\mathrm{lab}}t=z_{o}$. Conversely, $V=V_{\mathrm{lab}}/(1-V_{\mathrm{lab}}/c)$. Therefore, the requirement $V_{\mathrm{lab}}\leqslant c$ of physical signals translates into the requirement $V\geqslant 0$. Moreover, $V_{\mathrm{lab}}\to 0$ implies $V\to 0$ and $V\to\infty$ as $V_{\mathrm{lab}}\to c$.

6.2 Soliton solutions with inhomogeneous broadening: propagation

As usual, when studying the case of inhomogeneous broadening we take the spectral-line shape to be a Lorentzian, namely:

[TABLE]

The propagation of the optical background $q_{-}(z)$ is governed by Eq. (90). But because $g(\lambda)$ and $\gamma(\lambda)$ are respectively an even and an odd function of $\lambda$, the integral in the definition (89) is zero, i.e., $w_{-}(z)=0$. Thus, the initial state for $q(t,z)$ is independent of $z$, namely $q_{-}(z)=q_{o}$. We point out that a more general scenario, i.e., solutions with a shifted Lorentzian, is discussed in Appendix 7.

The propagation of the norming constant $C_{1}$ is governed by ODE (110). Recall that we parameterized the norming constant in Eq. (134). In order to obtain $\xi(z)$ and $\varphi(z)$ one must compute $R_{\pm,d}$ explicitly by using Eq. (98). In this case the calculations are rather involved, and are presented in Appendix A.8. Since we will substitute the discrete eigenvalue into $R_{-,d}(\zeta,z)$, we know $\zeta\in\Gamma^{+}$, and therefore we consider the case with $\zeta\in\mathbb{C}/\mathbb{R}$, i.e., $\lambda\in\mathbb{C}/\mathbb{R}$. At the end, we write the matrix in the following form

[TABLE]

with $\Theta(\lambda)$ defined as

[TABLE]

Importantly, in Appendix A.8, we also show that $R_{-,d}$ has a different form on the continuous spectrum $\Sigma$ as follows,

[TABLE]

Therefore, the matrix $R_{-,d}$ is discontinuous across the real $\zeta$-axis, i.e., the real $\lambda$-axis. (This is the same as in the case of ZBG, and follows immediately from its definition as a principal value integral.) Clearly, the value of $R_{-,d}$ (152) depends on whether $\varepsilon<q_{o}$ or $\varepsilon>q_{o}$.

Using the expression above, we can now compute the propagation for the norming constants. Recalling the parametrization (134), we have that

[TABLE]

Explicitly, from Eq. (110) and (154) we have

[TABLE]

It is easy to show that the one-soliton solution (141) oscillates with temporal frequency

[TABLE]

In particular, if $\alpha=\pi/2$ (traveling-wave soliton), simple calculations show that $\omega=0$, i.e., there are no oscillations. This is why we name the solitons as traveling-wave or oscillatory solitons.

Two traveling-wave soliton solutions are shown in Figs. 5 and 5, whereas an oscillatory soliton solution is shown in Fig. 7. Notice that, since $\gamma(\lambda)$ is discontinuous at $\lambda=0$, the density matrix $\rho(t,z,\lambda)$ is also discontinuous there. Thus, in the above figures as well as the subsequent ones, we plot $\rho(t,z,0^{+})$. Moreover, we take $\lambda$ on the first sheet, i.e., we take $\nu=1$ in the initial state (126). Notice also that oscillatory behavior can be easily observed from Fig. 7 just like for the focusing NLS equation with NZBC. Moreover, one can easily observe from all the figures that $\lim_{t\to-\infty}D$ does not equal $h_{-}$, and is a function of $\lambda$ even though $h_{-}$ is not. This is consistent with our discussion in Section 3.4 and with Eqs. (83) and (84).

Now that we have derived the propagation of the norming constant, we can explicitly compute the amplitude and velocity of the soliton. In Fig. 7(left) we show the maximal possible amplitude of the soliton, $A_{\mathrm{max}}(z)$, as a function of $\eta$ and $\alpha$, which is computed using Eq. (148). Also, taking $\xi(z)$ from Eq. (156), if \mathop{\rm Im}\nolimits\{g(\lambda)\big{[}\Theta(\lambda)-\gamma\Theta(i\varepsilon)/(q_{o}^{2}-\varepsilon^{2})^{1/2}\big{]}\}\neq 0 we can obtain the soliton velocity (150) as

[TABLE]

where we use the notation from Section 6.1. In Fig. 7 (center and right), we show the dependence of the soliton velocity on the parameters $\eta$ and $\alpha$, in the cases $\varepsilon<q_{o}$ and $\varepsilon>q_{o}$, respectively. It is apparent that the velocity is an increasing function as $\eta\to\infty$. Thus, as the discrete eigenvalue moves closer to the real line, the corresponding soliton travels slower.

Note that, just like in the case of ZBG, the soliton velocity is symmetric with respect to the imaginary $\lambda$-axis. That is, solitons corresponding to discrete eigenvalues $\lambda_{n}$ and $-\lambda_{n}^{*}$ (i.e., with the same imaginary part and opposite real part) have the same velocity. (This is unlike what happens to the NLS equation, where such solitons would have opposite velocities.)

Importantly, note also that \mathop{\rm Im}\nolimits\{g(\lambda)\big{[}\Theta(\lambda)-\gamma\Theta(i\varepsilon)/(q_{o}^{2}-\varepsilon^{2})^{1/2}\big{]}\}>0 for all $|\zeta|>q_{o}$. Thus, V\mathrel{\mathchoice{\vbox{\offinterlineskip\halign{\m@th\displaystyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\textstyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\scriptstyle\hfil#\hfil\cr>\crcr<\crcr} }}{\vbox{\offinterlineskip\halign{\m@th\scriptscriptstyle\hfil#\hfil\cr>\crcr<\crcr} }}}0 (solitons are subluminal or superluminal) when $h_{-}=\mp 1$ (i.e., more atoms initially are in the ground than in the excited state or vise versa), respectively, similarly to the case with ZBG.

6.3 Soliton solutions in the sharp-line and zero-background limits

We next discuss the two limits $\varepsilon\to 0$ and $q_{o}\to 0$. The former is called sharp-line limit, and the latter is called zero-background limit. Importantly, in the development of the IST, one could take the limit $q_{o}\to 0$ throughout to recover the formalism in the case of ZBG. However, we will see shortly that the behavior of solutions differs depending on whether one lets $q_{o}\to 0$ before or after taking the limit $\varepsilon\to 0$. In other words, the limits $q_{o}\to 0$ and $\varepsilon\to 0$ do not commute. In practice, what this means is that one must know a priori whether the physical scenario involves $\varepsilon<q_{o}$ or $q_{o}<\varepsilon$.

Case 1: Taking the limit $q_{o}\to 0$ first.

Taking $q_{o}\to 0$ in Eq. (90) we immediately obtain $q_{-}=0$ for all $z>0$. Also, as it was shown in [69] that, in the limit $q_{o}\to 0$, the solution profile (136) reduces to a hyperbolic secant. Moreover, taking the limit $q_{o}\to 0$ in Eq. (152), we obtain

[TABLE]

Note that this expression is identical to the one in the case of ZBG [45]. Also, the expression (152) was derived on sheet I. The same applies to $\rho(t,z,\lambda)$. Thus, the entire solution reduces to that in the case of ZBG. One obtains a similar result while computing everything on sheet II.

If we now take the limit $\varepsilon\to 0$, we recover the sharp-line limit in the case of ZBG (in particular, $R_{-}=(h_{-}/\lambda)\,\sigma_{3}$).

Case 2: Taking the limit $\varepsilon\to 0$ first.

We now take the limit $\varepsilon\to 0$ first. Using Eq. (152), it is easy to show that the matrix $R_{-,d}(\zeta_{1})$ is identically zero in this limit. Therefore, the soliton solution is constant with respect to the spatial variable $z$.

If we now take the limit $q_{o}\to 0$, we preserve a soliton solution that is constant with respect to $z$, which is not the same as the soliton solution of the case with ZBG in the sharp-line limit. Thus, as anticipated, the two limits $\varepsilon\to 0$ and $q_{o}\to 0$ do not commute. This is consistent with the fact that the limit $\varepsilon\to 0$ is subtle, and must be carried out carefully (e.g., see [88]).

Direct sharp-line limit.

We can also compute the sharp-line limit directly by using the Dirac delta with the initial state for the medium (126), meaning that we take the spectral-line shape as follows

[TABLE]

Again, the calculation requires evaluating the integral in Eq. (98) on either $\lambda$-sheet. Importantly, with our definition (10) of $\gamma(\lambda)$, both $\zeta(\lambda)$ and $\gamma(\lambda)$ are discontinuous at $\lambda=0$. Therefore, we need to use identity (197) in the calculation. We then obtain $q_{-}(z)=q_{o}$ from the propagation equations (89) and (90). Calculating the integral in Eq. (98) we obtain $R_{-,d}=0$. Thus, the norming constant $C_{1}$ satisfies ${\partial C_{1}}/{\partial z}=0$ from the ODE (110), implying

[TABLE]

In other words, the soliton solution (136) is constant in $z$, and therefore the soliton travels at the speed of light. From these results, one can also reconstruct the density matrix $\rho(t,z,\lambda)$. This result is precisely the same as that obtained in case 2 above, namely by taking the limit $\varepsilon\to 0$ first.

6.4 Type 3. Periodic solutions

Here, we discuss the nontrivial limiting case of oscillatory soliton solutions (141), when the discrete eigenvalues tend to the branch cut of $\gamma$. More precisely, we take the limit as $\eta\to 1^{+}$ of the discrete eigenvalue $\zeta_{1}=q_{o}\eta\mathrm{e}^{i\alpha}$, as in the parameterization (134), with the phase $\alpha\in(0,\pi/2)\cup(\pi/2,\pi)$, and define $q_{\mathrm{P}}(t,z)=\lim_{\eta\to 1}q(t,z)$. [Several traces of the discrete eigenvalue as $\eta\to 1^{+}$ with fixed values of $\alpha$ are shown in Fig. 3(right).] After simple calculations, this solution is found to be

[TABLE]

where

[TABLE]

and $\xi(z)$ and $\varphi(z)$ will be discussed next. Notice that the $t$-dependence [in $s(t,z)$] of the solution (160) only appears in trigonometric functions. As a result, this solution is periodic in both $t$ and localized along the curve $\chi(z)=y$. We should point out that this periodic solution does not satisfy our initial data: $\lim\nolimits_{t\to\pm\infty}q(t,z)=q_{-}(z)$. This, however, is similar to what happens when the same limiting process is applied to the focusing NLS equation [69]. Indeed, solution (160) is the analogue of the so-called Akhmediev breathers [89] of the focusing NLS equation, which are spatially periodic. However, the Akhmediev breathers are not temporally periodic, but the solutions of the MBEs are. Of course, in our case the propagation depends on the choice of $\rho_{-,d}$ and the spectral-line shape $g(\lambda)$. We next compute the $z$ dependence explicitly with inhomogeneous broadening.

Periodic solution with inhomogeneous broadening.

Following the discussion in Section 6.2, taking the limit of Eq. (152) as $\eta\to 1$, we obtain the following quantities appearing in the soliton solution

[TABLE]

with

[TABLE]

Therefore, we have the following two quantities from Eqs. (156) and (157), which determine the norming constant,

[TABLE]

In this particular case, we can rewrite the periodic solution (160) explicitly as

[TABLE]

where

[TABLE]

with $\xi(0)$ and $\varphi(0)$ arbitrary real constants and $R_{\varepsilon}(\alpha)$ given by Eq. (166). The modified eigenfunction matrix $\mu$ can be computed in a similar way as in the previous cases and is omitted for brevity.

Due to the fact that both variables $t$ and $z$ appear in trigonometric functions, this solution $q_{\varepsilon}(t,z)$ is periodic in both $t$ and $z$ with frequencies $\omega_{t}=q_{o}\cos\alpha/\pi$ and $\omega_{z}=R_{\varepsilon}/(2\pi)$. Moreover, each peak travels with velocity $V=-2q_{o}h_{-}\cos\alpha/R_{\varepsilon}$. One such solution is shown in Fig. 8.

It is easy to compute the sharp-line limit, i.e., $\varepsilon\to 0$, of the periodic solution. By Eq. (165), we find that $R_{\varepsilon}\to 0$ as $\varepsilon\to 0$. Thus, the solution becomes

[TABLE]

with

[TABLE]

Notice that this sharp-line limit solution is independent of $z$. Again, notice that all the solutions discussed in this section are the Maxwell-Bloch analogue of the so-called “Akhmediev breathers” of the focusing NLS equation.

6.5 Type 4. Rational solutions

Recall that, for the focusing NLS with NZBC, one can obtain the Peregrine solution as a suitable limit of its one-soliton solutions [90]. We next show that it is possible to obtain rational solutions of the MBEs with NZBG in a similar way. We compute such a solution by taking the limit of the traveling-wave soliton solution (136) as the discrete eigenvalue tends to the branch point $iq_{o}$, i.e., $\eta\to 1$. (Recall that the discrete eigenvalue is $\zeta_{1}=iq_{o}\eta$.) Notice that this limit is different from the one in the previous section due to different phases of the discrete eigenvalue. To continue our calculation, we first expand the norming constant $C_{1}$ around $\eta=1$. Recall that we parameterized the norming constant $C_{1}$ in the parameterization (134). Therefore, we have the following expansions

[TABLE]

where

[TABLE]

and similarly for $s_{o}(z)$ and $s_{1}(z)$. In order to write down these two Taylor expansions explicitly, we need to compute the expansion of the matrix $R_{-,d}(\zeta_{1},z)$ as $\eta\to 1$, namely:

[TABLE]

Using Eq. (98), straightforward calculations show that the coefficients in the above expression are

[TABLE]

where $w_{-}(z)$ is defined in Eq. (89). Notice that $R^{(1)}(z)$ is nonzero in general, and can be computed on either $\lambda$-sheet. Let us recall the propagation equation for the norming constant (110) and the definition of $W_{-}(z)$ (90). The norming constant in this limit is given by $\lim_{\eta\to 1}C_{1}(iq_{o}\eta,z)=\mathrm{e}^{-2iW_{-}(z)}\lim_{\eta\to 1}C_{1}(iq_{o}\eta,0)$. Comparing this expression with Eq. (134), we find

[TABLE]

Therefore, both $\xi_{o}(z)$ and $s_{o}(z)$ are independent of $z$. Correspondingly, hereafter we will simply write them as $\xi_{o}$ and $s_{o}$.

Importantly, if one takes the limit $\eta\to 1$ with an arbitrary choice of ICs for $\xi_{o}$ and $s_{o}$, the soliton solution (136) reduces to the background solution (i.e., a trivial solution). On the other hand, choosing the ICs to be $\xi_{o}=-\Delta_{o}$ [recall $\Delta_{o}$ is defined in Eq. (140)] and $s_{o}=-\pi/2$, the limit of the solution (136) as $\eta\to 1$ yields a solution of the MBEs with rational dependence on $t$, namely:

[TABLE]

The center of the soliton and its velocity are given by the constraint $2q_{o}t+\xi_{1}(z)=0$. Furthermore, the density matrix is also a rational function of $t$. Indeed, computing the limit of the eigenfunction $\mu_{-}(t,z,\zeta)$ in Eq. (132) using the same methods as above, we find

[TABLE]

with $I$ the $2\times 2$ identity matrix and

[TABLE]

The $z$-dependence of all the entries was omitted for brevity. One can now use (133) to obtain the density matrix $\rho(t,z,\lambda)$. The final result is omitted for brevity due to its complexity. We refer to the solution (173) and the density matrix generated by Eqs. (174) and (133) as the rational solution of the MBEs with NZBG.

In the limit $q_{o}\to 0$ of the above solution, $q_{\mathrm{R}}(t,z)\to 0$ and $\mu_{\mathrm{R}}(t,z,\zeta)\to I$, implying $\rho_{\mathrm{R}}(t,z,\zeta)\to\rho_{-}(z,\zeta)$. Therefore, the rational solution reduces to the trivial solution of the MBEs with ZBG.

Rational solution with inhomogeneous broadening.

Recall that the propagation of the above solution (173) is determined by three quantities $W_{-}(z)$, $\xi_{1}(z)$ and $s_{1}(z)$. The latter two are obtained from the second term $R_{1}(z)$ in Eq. (172). We turn to calculating them next.

We again use the Lorentzian spectral-line shape (151). Recall from the discussion in Section 6.2 that in this case $w_{-}(z)=0$. Then from Eq. (172), we obtain the explicit expansion (171) for the matrix $R_{-,d}$ as $\eta\to 1$, with

[TABLE]

which again, can be computed on either sheet. Notice that $h_{-}=\pm 1$ denotes the initial state of the medium as in Eq. (126), and

[TABLE]

where $\Theta(\lambda)$ is defined in Eq. (153). Notice that $g_{o}$ and $\tilde{R}_{1}$ are both real, and both $R^{(0)}(z)$ and $R^{(1)}(z)$ are independent of $z$. By using Eqs. (110), (169) and (170), we then find

[TABLE]

Substituting all the components into Eq. (173), we then obtain the fully explicit rational solution

[TABLE]

where $\tilde{\xi}_{1}$ and $\tilde{s}_{1}$ are arbitrary real constants, and both $g_{o}$ and $\tilde{R}_{1}$ are given in Eq. (176). Notice that this solution is a traveling wave. In particular, if $\tilde{s}_{1}=0$, solution (177) is purely real. The parameter $\tilde{\xi}_{1}$ determines a spatial/temporal displacement, while $\tilde{s}_{1}$ determines a “phase shift”. The density matrix $\rho(t,z,\zeta)$ can also be calculated from Eq. (174), but is omitted for brevity due to its complexity. One such solution is shown in Fig. 9.

The amplitude compared to the background and the velocity of such a solution are, respectively,

[TABLE]

The maximal possible value for $A$ is $A_{\mathrm{max}}=2q_{o}$ when $\tilde{s}_{1}=0$.

It is easy to compute the sharp-line limit of the above rational solution. In the limit $\varepsilon\to 0$, we have $g_{o}\tilde{R}_{1}\to 0$, while all other parts of the solution remain the same as the ones in the inhomogeneous broadening case, implying

[TABLE]

Therefore, the solution travels at the speed of light in vacuum. Also, we can use relation (133) to write down the density matrix with $\tilde{\xi}_{1}=0$ explicitly

[TABLE]

where $\tilde{y}=4q_{o}^{2}t^{2}+\tilde{s}_{1}^{2}+1$ and $\sigma_{j}$ are defined in Eq. (5). Two such solutions are shown in Figs. 11 and 11.

6.6 On the stability of the soliton solutions

We now discuss the stability of general soliton solutions. To do so, we use a general form for the initial state $\rho_{-}$ consistent with Eqs. (77), (78), (79) and (80), which we rewrite in a convenient form as

[TABLE]

where $\nu=1$ for $\zeta\in(-\infty,-q_{o}]\cup[q_{o},\infty)$, $\nu=-1$ for $\zeta\in(-q_{o},q_{o})$. Functions $\beta(z)$ and $\phi(\zeta,z)$ are real-valued, and ensure that $\det\rho_{-}(z)=-1$. Moreover, $\phi(-q_{o}^{2}/\zeta,z)=-\phi(\zeta,z)+4W_{-}(z)$ with $W_{-}(z)$ defined in Eq. (90) and $\zeta\in(-\infty,-q_{o}]\cup[q_{o},\infty)$, because of the enforced symmetries (77), (78), (79) and (80). Recall the discussion in Section 3.4, $\beta(z)$ relates to the initial amount of population inversion of this system. If $0\leqslant\beta(z)<\pi/2$, initially, there are more atoms in the excited state than in the ground state, whereas if $\pi/2<\beta(z)\leqslant\pi$, initially, there are more atoms in the ground state than in the excited state. The special values $\beta(z)=0$ and $\beta(z)=\pi$ identify the situation where the largest proportion of the atoms is initially in the excited and ground states, respectively. Moreover, we consider the inhomogeneous broadening case, i.e., we take $g(\lambda)$ as in Eq. (151).

We first point out that the stability is characterized by the reflection coefficient $b(\zeta,z)$ appearing in the jump condition of the RHP (113). Therefore, we need to analyze the quantity $b(\zeta,z)$ whose dependence on $z$ is governed by the ODE (107). Since the coefficients of this ODE contain the entries of $R_{-,d}$, we consider the matrix $R_{-,d}$, as defined in Eq. (98). Because $\rho_{-,d}(\zeta,z)=\nu\,\cos\beta(z)\sigma_{3}$ from Eq. (179), $R_{-,d}$ can be calculated in the same way as in Appendix A.8. We thus obtain Eqs. (152) and (154) for $\lambda\in\mathbb{C}\backslash\mathbb{R}$ and $\lambda\in\mathbb{R}$, respectively, on either $\lambda$-sheet, with $h_{-}=\cos\beta(z)$. It is then convenient to write the resulting expression as a function of $\zeta$ in the form of

[TABLE]

with $\bar{R}_{\mathrm{re}}$ and $\bar{R}_{\mathrm{im}}$ both real. Comparing Eqs. (100), (101) and (108), we then obtain

[TABLE]

where $\nu_{o}=1$ if $\zeta\in\mathbb{R}$ and $\nu_{o}=0$ if $\zeta\in C_{o}$. Using the propagation equation (107) for the reflection coefficient, we find the following explicit expression

[TABLE]

where the primes denote the integration variable (not differentiation) as before, and where $\nu_{o}=1$ if $\zeta\in\mathbb{R}$, $\nu_{o}=0$ if $\zeta\in C_{o}$, and

[TABLE]

Note that Eq. (181) on the real axis is formally the same as the one in the case of ZBG. If $b(\zeta,0)$ (i.e., the reflection coefficient at $z=0$) and $\sin\beta(z)$ are both zero, then $b(\zeta,z)$ is identically zero (i.e., there is no radiation). We assume this is not the case.

Notice that the reflection coefficient (181) derived from the inhomogeneous ODE contains a homogeneous and an inhomogeneous portion, and $\Omega(z,z^{\prime},\zeta)$ appears in both of them. As a result, the long-distance behavior of the reflection coefficient $b(\zeta,z)$ depends on the absolute value of $\Omega(z,z^{\prime},\zeta)$. The situation is further complicated by the fact that both $\beta(z)$ and $\phi(\zeta,z)$ are in principle arbitrary. For simplicity, and by analogy with the case with ZBG, here we will study the simpler situation in which $\beta(z)$ and $\phi(\zeta,z)$ are constants with respect to $z$ and in which $g(\lambda)$ is the Lorentzian function (151). In this case, the expression (181) simplifies to:

[TABLE]

where again, $\nu_{o}=1$ if $\zeta\in\mathbb{R}$ and $\nu_{o}=0$ if $\zeta\in C_{o}$, $\bar{R}_{\mathrm{re}}$ and $\bar{R}_{\mathrm{im}}$ given by Eq. (180) as before. In this situation, the growth/decay of $b(\zeta,z)$ depends on the modulus of \exp\big{[}(\bar{R}_{\mathrm{im}}+\nu_{o}\pi g(\lambda)-i\bar{R}_{\mathrm{re}})\cos\beta\,z\big{]}, so the most important part by now is the unknown quantity $\bar{R}_{\mathrm{im}}$. Since any complex phase of the exponential produces inessential oscillations, the quantity $\bar{R}_{\mathrm{re}}$ is therefore not crucial for determining the long-distance behavior of the reflection coefficient. By this discussion, we next turn to characterize the imaginary part of the matrix $R_{-,d}$ given by Eqs. (152) and (154), depending on the location of $\zeta$. It is necessary to discuss separately two cases depending on whether $q_{o}>\varepsilon$ or $q_{o}<\varepsilon$.

Case I: $0\leqslant q_{o}<\varepsilon$.

Firstly, we discuss the case where the amplitude of the background field is small relative to the width of the spectral-line shape. In this case $\Theta(i\varepsilon)=2\mathrm{arcsech}(\varepsilon/q_{o})$ is pure imaginary. We can then rewrite $R_{-,d}$ from Eqs. (152) and (154) as

[TABLE]

where

[TABLE]

$C_{o}$ is the complex portion of the continuous spectrum and $g(\lambda)$ is given in Eq. (151) as before. Notice that the function $C(\lambda)$ is always real on the continuous spectrum, because $\zeta\in\mathbb{R}$ implies $\lambda\in\mathbb{R}$ and $\zeta\in C_{o}$ implies $\lambda\in i[-q_{o},q_{o}]$. Also $\gamma$ is real for all $\zeta\in\Sigma$. We then conclude that the quantity $R_{-,d}(\zeta,z)$ is always real. In other words, $\bar{R}_{\mathrm{im}}=0$, for all $\zeta\in\Sigma$.

Case II: $0<\varepsilon<q_{o}$.

Next, we consider the case in which the width of the spectral-line shape is small compared to the amplitude of the background optical field. In this case, $\Theta(i\varepsilon)=2\mathrm{arcsech}(\varepsilon/q_{o})$ is real. Similarly to the previous case, we rewrite the matrix from Eqs. (152) and (154) into the form of Eq. (180). We obtain

[TABLE]

where $C(\lambda)$ is defined in Eq. (185). Following a similar discussion as in case I, we know that the above expression is purely real, implying $\bar{R}_{\mathrm{im}}=0$ for all $\zeta\in\Sigma$.

Implications for stability.

In summary, we have that $\bar{R}_{\mathrm{im}}=0$ on the continuous spectrum, no matter what the values of $\varepsilon$ and $q_{o}$ are. Therefore, the solution (183) further reduces to

[TABLE]

where $\bar{R}_{\mathrm{re}}$ can be computed from Eq. (184) or (186) and is omitted here due to its lack of importance. Obviously, the reflection coefficient $b(\zeta,z)$ is always bounded on the portion $\zeta\in C_{o}$ for $z>0$, so we next turn to discuss the case $\zeta\in\mathbb{R}$.

It is easy to see that $b(\zeta,z)$ (and therefore the radiation) grows exponentially as $z\to\infty$ if $\cos\beta>0$, whereas $b(\zeta,z)$ (and therefore the radiation) decays exponentially as $z\to\infty$ if $\cos\beta<0$. Therefore, the stability of general solutions with NZBG is identical to the case of ZBG. Namely, radiation decays (implying that the background and the solitons are stable) when initially the atoms in the ground state dominate (population inversion is negative); radiation grows (implying that the background and the solitons are unstable) when initially the atoms in the excited state dominate (population inversion is positive).

In particular, for a purely diagonal initial state (i.e., $\sin\beta=0$) without any discrete eigenvalues, the solution of MBE with NZBG tends to the background value when the medium initially has the lowest population inversion. That is, $q(t,z)=q_{-}(z)+o(1)$ as $z\to\infty$ if $\rho_{-}=-\nu\sigma_{3}$ and if no solitons are present.

7 Soliton solutions with a shifted Lorentzian

In this section we present one-soliton solutions (141) obtained when the shape of spectral line is a shifted Lorentzian, namely

[TABLE]

where $\lambda_{o}$ is the center of the Lorentzian. In the sharp-line limit ($\varepsilon\to 0$), Eq. (187) becomes a shifted Dirac delta $g(\lambda)=\delta(\lambda-\lambda_{o})$.

Note that changing the shape of spectral line only affects the $z$-dependence of solutions. Therefore, all the calculations for the pure soliton solutions in Sections 5.2 and 6.1 remain valid. On the other hand, since the function (187) is not even, the calculations are more complex than with the standard Lorentzian (151). As a result, here we present the calculations in more detail than in Section 6.2. More importantly, we show that the modified $z$-dependence results in drastically different dynamical behavior of the soliton solutions.

Asymptotic background.

In the case of a standard Lorentzian, $q_{-}(z)$ is constant with respect to $z$ (cf. Section 6.2). This will not be true anymore in our case, however. Recall the relevant equations are (88), (89) and (90). Moreover, recall that the choice of $\rho_{-}$ resulting in pure soliton solutions is (126). In other words, on the first $\lambda$-sheet, $D_{-}=h_{-}$ with $h_{-}=\pm 1$ and $\lambda\in\mathbb{R}$. For this value of $D_{-}$, Eq. (89) becomes

[TABLE]

We calculate this $w_{-}$ as follows [where we recall $\gamma=\mathop{\rm sign}\nolimits(\lambda)\sqrt{\lambda^{2}+q_{o}^{2}}$ for $\lambda\in\mathbb{R}$]

[TABLE]

Performing the change of variable $u=\sqrt{\lambda^{2}+q_{o}^{2}}$, we then obtain

[TABLE]

Importantly, we point out that Eq (188) implies $w_{-}\in\mathbb{R}$ for all values of $q_{o}$, $h_{1}$, $\varepsilon$ and $\lambda_{o}$, even though this is not trivially seen from Eq. (189). If we consider the special case $\lambda_{o}=0$ (an unshifted Lorentzian), the two terms in Eq. (189) cancel out exactly, and we recover the result $w_{-}=0$ of Section 6.2. We next consider the limit of Eq. (189) as $\varepsilon\to 0$. After a careful calculation which takes into account the branch cut of the complex logarithmic function, we obtain the limit of Eq. (189) as

[TABLE]

This result coincides with direct calculations of the sharp-line limit where $g(\lambda)=\delta(\lambda-\lambda_{o})$. Thus, the $z$-dependence of the background is given by Eq. (91) with $w_{-}$ given by Eq. (189). Thus, with a shifted Lorentzian, the background value is not independent of the propagation variable $z$.

The auxiliary matrix and propagation of the norming constant.

We next discuss the $z$-dependence of the norming constant, which is governed by Eq. (110). As a result, we need to compute the auxiliary matrix $R$ first. In particular, we next calculate $R_{-,1,1}$ from Eq. (100). Notice that in Eq. (100), the quantity $R_{-,1,1}$ is expressed as a principal value integral, but later on is evaluated at $\zeta=\zeta_{n}$, i.e., at that discrete spectrum. In other words, eventually, the quantity $R_{-,1,1}$ is evaluated off the real line. Thus, the principal value is not needed for soliton solutions. Equation (100) becomes

[TABLE]

where $D_{-}=h_{-}=\pm 1$. We have already calculated the quantity $w_{-}$, so we only need to consider the integral in the first term. We name this integral $I_{1}$ and we can rewrite it as

[TABLE]

where we rename the integration variable $\lambda^{\prime}$ to $u$ for clarity. With a change of variable $x=\sqrt{u^{2}+q_{o}^{2}}$, the above integral becomes

[TABLE]

Decomposing this integral into partial fractions, using standard integration formulae, and combining all pieces, we then find

[TABLE]

Note that, in the limit $\lambda_{o}\to 0$ (an unshifted Lorentzian), we have

[TABLE]

This expression is equivalent to Eq. (152). We now have all the components needed to determine the $z$-dependence of the soliton solutions. Recall that there are four types of one-soliton solutions depending on the location of the discrete eigenvalue (cf. Fig. 3 left). Next we discuss all of them, and we show how they differ from the ones in the main text.

Soliton solutions of Types 1–3.

Recall that the norming constant is parameterized by Eq. (155), where $\zeta_{1}$ is the discrete eigenvalue and $R_{-,1,1}(\zeta)$ is given by Eq. (191). The formulas for the first three types of soliton solutions all remain the same except for the $z$ dependence. Expliticly, Eq. (136) for Type I, Eq. (141) for Type II and Eq. (160) for Type III all remain valid except that one simply substitutes $w_{-}$ from Eq. (189), and $\xi(z)$ and $\varphi(z)$ from Eq. (155) with $R_{-,1,1}$ from Eq. (191). The three types of solutions with a shifted Lorentzian are shown in Figs. 13, 13 and 15, respectively.

Figure 13 shows Type-1 soliton solutions obtained with a few different choices of $\lambda_{o}$. It is evident from the figure that, as the peak of detuning function is shifted from the origin to infinity, the solution becomes oscillatory. As a result, this is not a simple traveling wave anymore, but becomes a breather-like solution. Eventually, as $\lambda_{o}\to\infty$ the spatial period of the oscillation increases without bound. Moreover, it appears that, as $\lambda_{o}$ gets large, the soliton velocity tends to the speed of light. This is to be expected: as the atoms are detuned from the resonant frequency, the nonlinear interactions between atoms and light become weaker, so the medium becomes essentially transparent to the light pulse, which then becomes regular light traveling with speed $c$.

As shown in Fig. 13, Type-2 solitons retain the main characteristics of the corresponding solutions with an unshifted Lorentzian, namely a localized, traveling breather-like solution with an internal oscillatory structure. However, as the center of the Lorentzian approaches infinity, the soliton velocity, i.e., the group velocity, and the phase velocity both tend to the speed of light. Eventually, as $\lambda_{o}\to\infty$ this type of soliton solution becomes a localized traveling wave.

On the other hand, as shown in Fig. 15, there is a dramatic change in the characteristics of the solution for Type-3 solitons. Recall that, if $\lambda_{o}=0$ (i.e., for an unshifted Lorentzian), these solutions are periodic with respect to both variables $t$ and $z$ (cf. Section 6.4). If $\lambda_{o}\neq 0$, however, Type-3 soliton solutions are only periodic in $t$, and are not periodic in $z$ in general. This is because when $\lambda_{o}\neq 0$ and $\zeta_{1}=\mathrm{e}^{i\alpha}$, the quantity $R_{-,1,1}(\zeta_{1})$ is not purely real anymore. So both $\xi(z)$ and $\varphi(z)$ acquire a nontrivial dependence on $z$. Moreover, when $\lambda_{o}\neq 0$, the solution decays to the background as $z\to\infty$, i.e., the solution becomes localized in $z$. As a result, the character of this solution becomes similar to that of the Akhmediev breathers of the NLS equation.

Type 4. rational solutions.

The general expression (173) for the rational solutions still applies to the current case. The difference from before is in the spatial dependence of the three quantities $q_{-}(z)$, $\xi_{1}(z)$ and $s_{1}(z)$. As before, the background $q_{-}(z)$ is obtained from Eq. (90) with $w_{-}$ given by Eq. (189). The quantities $\xi_{1}(z)$ and $s_{1}(z)$ are obtained from the Taylor expansions (169) and (170). To calculate these quantities explicitly, we first expand $R_{-,1,1}(iq_{o}\eta,z)$ similarly to Eq. (171) as

[TABLE]

where $R_{-,1,1}^{(0)}=2w_{-}\in\mathbb{R}$ and

[TABLE]

Then, from Eq. (155), we can expand $\xi$ and $\varphi$ as

[TABLE]

as $\eta\to 1$, with

[TABLE]

and where $\xi_{o}^{(j)}$ and $\varphi_{o}^{(j)}$ for $j=0,1$ are arbitrary real constants. Following a similar procedure from Section 6.5, we can write down the rational solution as

[TABLE]

where $w_{-}$ is given by Eq (189) and $\xi^{(1)}(z)$ and $\varphi^{(1)}(z)$ are given by Eqs. (193) and (194), with $R_{-,1,1}^{(1)}$ given by Eq. (192). As before, a few examples of the resulting solutions are shown in Fig. 15 for a few different choices of $\lambda_{o}$. Importantly, it should be clear from Eq. (195) and Fig. 15 that, when $\lambda_{o}\neq 0$, Type-4 soliton solutions are not periodic traveling waves, but rather a localized structure in both $t$ and $z$, and exhibits rogue-wave characteristics. This is also evident from Fig. 15. In other words, rational solutions with $\lambda_{o}\neq 0$ are similar to the Peregrine solutions of the focusing NLS equation. Indeed, the very expression of the rational solution (195) is essentially identical to that of the Peregrine solution of the focusing NLS equation.

8 Final remarks

In conclusion, we formulated the inverse scattering transform for the Maxwell-Bloch system of equations describing resonant interaction between light and active optical media in the case when the light intensity does not vanish in the distant past and future. We characterized the asymptotics of the density matrix. We then used the formalism to compute explicitly the soliton solutions of the system. We obtained a representation for the soliton solutions in determinant form, and we wrote down explicitly the one-soliton solutions. Afterwards, we also derived periodic solutions and rational solutions. We then analyzed the properties of these solutions, we discussed the sharp-line and small-amplitude limits, and we show that the two limits do not commute. Finally, we investigated the behavior of radiation, showing that the background and solitons are stable (i.e., radiative part of the solution is decaying) when the initial population inversion is negative, which is similar to the case with ZBG.

An important future step will be to generalize the IST theory of this work in order to study coupled Maxwell-Bloch systems (and in particular system describing the so-called “lambda” configuration) with NZBG, using the recent results of [91]. The IST for such systems in the case with ZBG was done in [52, 88], and was further studied in [92]. The study of the coupled MBEs with NZBG may be instrumental to develop a mathematical theory of the phenomenon of slow light [5]. In this respect, we note that a few solutions, obtained using direct methods were reported in [63, 64], but a comprehensive description of the phenomenon is still missing.

Acknowledgments

We thank Barbara Prinari for many interesting discussions on topics related to this work. This work was partially supported by the National Science Foundation under grant numbers DMS-1615524 and DMS-1615859.

Appendix

A.1 Notation

Let $\delta(\lambda)$ be a suitable representation of the Dirac delta. If $f(\lambda)$ is continuous at $\lambda=0$,

[TABLE]

On the other hand, if $f(\lambda)$ has a jump discontinuity at $\lambda=0$,

[TABLE]

For any 2$\times$2 matrix $M$, one has

[TABLE]

where $\sigma_{*}$ is defined in Eq. (45).

For reference, note that

[TABLE]

A.2 Symmetries

We consider the two symmetries separately.

1. $\zeta\mapsto\zeta^{*}$  (upper/lower half plane),  implying  $(\lambda,\gamma)\mapsto(\lambda^{*},\gamma^{*})$  (same sheet). It can be shown that, if $\phi(t,z,\zeta)$ is a solution of the scattering problem (6), so is

[TABLE]

where $\sigma_{*}$ is defined in Eq. (45). [That is because  $\sigma_{*}^{2}=-I$,  $\sigma_{*}\sigma_{3}\sigma_{*}=-\sigma_{3}$  and  $\sigma_{*}Q^{*}\sigma_{*}=Q^{\dagger}=-Q$.] By the same token, so is $\tilde{\phi}\,C$, where $C$ is any constant $2\times 2$ matrix. Now restrict our attention to $\zeta\in\Sigma$. Let $\phi\equiv\phi_{\pm}$, and look at the asymptotic behavior of $\tilde{\phi}_{\pm}$ as $t\to\pm\infty$: Recall that $\phi(t,z,\zeta)\sim Y_{\pm}(\zeta)\,\mathrm{e}^{i\gamma(\zeta)t\sigma_{3}}$ as $t\to\pm\infty$. Note that $\gamma^{*}(\zeta^{*})=\gamma(\zeta)$ and $\sigma_{*}Y_{\pm}^{*}(\zeta^{*})=\sigma_{*}(I-i\sigma_{3}Q_{\pm}^{*}/\zeta)=Y(\zeta)\sigma_{*}$. Also, a little algebra shows that $\sigma_{*}\mathrm{e}^{ia\sigma_{3}}\sigma_{*}=-\mathrm{e}^{-ia\sigma_{3}}$. Therefore, $\tilde{\phi}_{\pm}(t,z,\zeta)\,\sigma_{*}\sim Y_{\pm}(\zeta)\,\mathrm{e}^{i\gamma(\zeta)t\sigma_{3}}$ as $t\to\pm\infty$. But since the solution of the scattering problem with given BCs is unique, we must have $\tilde{\phi}_{\pm}\,\sigma_{*}=\phi_{\pm}$. Hence $\forall\zeta\in\Sigma$ we have

[TABLE]

Recalling Eq. (198), the desired results for the eigenfunctions are obtained by rewriting the symmetry (200). Recalling the scattering relation (38) and using symmetry (200), we have, $\forall\zeta\in\Sigma$,

[TABLE]

Using Eq. (198), one can compute the corresponding symmetries for the scattering matrix.

2. $\zeta\mapsto-q_{o}^{2}/\zeta$  (outside/inside $C_{o}$),  implying  $(\lambda,\gamma)\mapsto(\lambda,-\gamma)$  (opposite sheets). Since $\lambda(-q_{o}^{2}/\zeta)=\lambda(\zeta)$, it is easy to show that if $\phi(t,z,\zeta)$ is a solution of the scattering problem, so is

[TABLE]

As before, this implies that $w\,C$ is also a solution, for any matrix $C$ independent of $t$. With $\phi=\phi_{\pm}$, it is $\tilde{\phi}_{\pm}(t,z,\zeta)C\sim Y_{\pm}(-q_{o}^{2}/\zeta)\,\mathrm{e}^{-i\gamma(\zeta)\,t\sigma_{3}}C$ as $t\to\pm\infty$, because $\gamma(-q_{o}^{2}/\zeta)=-\gamma(\zeta)$. Now note that  $Y_{\pm}(-q_{o}^{2}/\zeta)\,\mathrm{e}^{-i\gamma t\sigma_{3}}\,\sigma_{3}Q_{\pm}=-i\zeta Y_{\pm}(\zeta)\,\mathrm{e}^{i\gamma t\sigma_{3}}$ [because $\mathrm{e}^{-ia\sigma_{3}}Q=Q\mathrm{e}^{ia\sigma_{3}}$ and $(\sigma_{3}Q)^{2}=q_{o}^{2}I$]. Thus, taking  $C=(i/\zeta)\,\sigma_{3}Q$, we obtain the desired results for $\phi_{\pm}$ in this case. For the scattering matrix, from the relation (38) and by using symmetry (48) and Eq. (199), we finally obtain symmetry (50).

1+2. We can also combine the results of the first two symmetries to get

[TABLE]

The corresponding symmetry for the columns follows as before.

A.3 Asymptotics of the eigenfunctions as $\zeta\to\infty$ and $\zeta\to 0$

Here we will prove the asymptotic expansions (66), (67), (68) and (69).

We start deriving the asymptotic behavior (66) and (67). The derivation proceeds by induction. The result is trivially true for $\mu_{d}^{(0)}$ and $\mu_{o}^{(0)}$. Moreover, using definition (28) and separating the diagonal and off-diagonal parts of the expansion (63), we have

[TABLE]

[TABLE]

As $\zeta\to\infty$, the four terms on the RHS of Eq. (202) are, respectively: $O(\mu_{o}^{(n)})$, $O(\mu_{d}^{(n)}/\zeta)$, $O(\mu_{d}^{(n)}/\zeta^{2})$ and $O(\mu_{o}^{(n)}/\zeta^{3})$. The last two estimates are obtained using integration by parts, taking advantage of the fact that $\gamma(\zeta)=\zeta/2+O(1/\zeta)$ as $\zeta\to\infty$ from Eq. (14). Similarly, as $\zeta\to\infty$, the four terms on the RHS of Eq. (203) are, respectively: $O(\mu_{o}^{(n)}/\zeta)$, $O(\mu_{d}^{(n)}/\zeta^{2})$, $O(\mu_{d}^{(n)}/\zeta)$ and $O(\mu_{o}^{(n)}/\zeta^{2})$ (where again the last two estimates are obtained using integration by parts). $\square$

Next we describe the asymptotic behavior (68) and (69). Again, the result is derived by induction. The result is trivially true for $\mu_{d}^{(0)}$ and $\mu_{o}^{(0)}$. Decomposing expansion (63) into its diagonal and off-diagonal parts yields Eqs. (202) and (203) as before. Finally, one shows that the four terms on the RHS of Eqs. (202) and (203) are, respectively: for Eq. (202): $O(\zeta^{2}\mu_{o}^{(n)})$, $O(\zeta\,\mu_{d}^{(n)})$, $O(\zeta^{2}\mu_{d}^{(n)})$ and $O(\zeta\,\mu_{o}^{(n)})$; for Eq. (203): $O(\zeta\mu_{o}^{(n)})$, $O(\mu_{d}^{(n)})$, $O(\zeta^{3}\mu_{d}^{(n)})$ and $O(\zeta\,\mu_{o}^{(n)})$. Again, some estimates are obtained using integration by parts and $\gamma(\zeta)=O(1/\zeta)$ as $\zeta\to 0$ from Eq. (14). $\square$

A.4 Calculation of $R_{\pm}$

We first establish the following result, which will be needed in the calculation:

Lemma A.1

[TABLE]

where $\nu=\pm 1$ when $\lambda$ on sheet I or II, respectively, and where we used the shorthand notation $\gamma=\gamma(\lambda)$ and $\gamma^{\prime}=\gamma(\lambda^{\prime})$.

Proof: Notice that for both $\lambda^{\prime}$ and $\lambda$, we have two choices: sheet I and sheet II. Here, we compute the integral in the way that $\lambda^{\prime}$ and $\lambda$ are always on the same sheet. In this case, we can eliminate the mixed case: $\lambda^{\prime}$ is on sheet I and $\lambda$ is on sheet II or vise versa. We choose such a calculation because, later we will show that this version of IST naturally reduces to the ZBCs as $q_{o}\to 0$.

Denote the LHS of the definition (204) by $I_{\pm}(\zeta)$. Note that, since $\gamma(\lambda)$ takes on different values on each sheet of the Riemann surface, so do $I_{\pm}$. Performing the change of variable $\lambda^{\prime}\mapsto\gamma^{\prime}$ we have

[TABLE]

where again $\nu=\pm 1$ when both $\lambda$ and $\lambda^{\prime}$ are on sheet I or sheet II, respectively, $L=(-\infty,-q_{o})\cup(q_{o},\infty)$ and

[TABLE]

Further letting $y=\gamma^{\prime}-\gamma$ in Eq. (205), we have

[TABLE]

where $L^{\prime}=(-\infty,-q_{o}-\gamma)\cup(q_{o}-\gamma,\infty)$. From here, we need to discuss two cases depending on whether $0\in L^{\prime}$ or $0\not\in L^{\prime}$.

• •

If $\lambda\in\mathbb{R}$, i.e., $\gamma\in(-\infty,-q_{o}]\cup[q_{o},\infty)$, then we know $0\in L^{\prime}$ and we obtain

[TABLE]

from the formula $\lim_{t\to\pm\infty}\mathop{\mathop{\textstyle\mathop{\textstyle\int}\limits}\limits\kern-9.40002pt-\kern 1.99997pt}\limits\nolimits_{-\infty}^{\infty}\mathrm{e}^{\pm iyt}F(y)\,\mathrm{d}y/y=\pm i\pi\,F(0)$. Equation (204) then follows by noting that $\lim_{\gamma^{\prime}\to\gamma}(\lambda^{\prime}-\lambda)/(\gamma^{\prime}-\gamma)=d\lambda/d\gamma$, and therefore $\lim_{\gamma^{\prime}\to\gamma}j(\zeta,\zeta^{\prime})=1$.

• •

If $\lambda\in i[-q_{o},0)\cup i(0,q_{o}]$, i.e, $\zeta\in C_{o}$ and $-q_{o}<\gamma<q_{o}$, then $0\not\in L^{\prime}$. By using the Riemann-Lebesgue lemma, we know $I_{\pm}(\zeta)=0$. $\square$

Moreover, note that if the term $(\gamma^{\prime}-\gamma)t$ in the exponent of Eq. (204) were replaced by $\gamma^{\prime}+\gamma$ or $\gamma^{\prime}t$, one can show that the integral is zero using the Riemann-Lebesgue lemma.

We now proceed to compute $R_{\pm}$. We rewrite Eq. (97) for convenience as

[TABLE]

We will first evaluate the two terms on the RHS separately. Note first that, from Eq. (92) we get, as $t\to\pm\infty$,

[TABLE]

where $w_{\pm}$ was defined in Eq. (89), implying $\phi_{\pm}^{-1}(\phi_{\pm})_{z}=(\phi_{\pm}^{-1}(\phi_{\pm})_{z})_{\mathrm{eff}}+o(1)$, where $(\phi_{\pm}^{-1}(\phi_{\pm})_{z})_{\mathrm{eff}}$ denotes the leading order term in the limit given by

[TABLE]

Moreover, recalling definitions (28) and (30), we have

[TABLE]

Combining the first two terms on the RHS and recalling that $Q_{\pm}^{2}=-q_{o}^{2}I$ we then have

[TABLE]

where $\sigma_{1}$ was defined in Eq. (5). Note that $(\phi_{\pm}^{-1}(\phi_{\pm})_{z})_{\mathrm{eff}}\to 0$ as $q_{o}\to 0$, as it should be, since then the BCs for $\phi_{\pm}$ are independent of $z$. On the other hand, the second term on the RHS of Eq. (207) does not have a finite limit as $t\to\pm\infty$. One would therefore expect that this term will be canceled by the first term in Eq. (97). We will see that this is indeed the case. The calculation of the first term of $R_{\pm}$ is considerably more involved, however.

To calculate the first term in $R_{\pm}$, we first note that, as $t\to\pm\infty$,

[TABLE]

where $d$ and $o$ denote the diagonal and off-diagonal parts of a matrix, as before, and

[TABLE]

Now note that, as $t\to\pm\infty$,

[TABLE]

where as before we used primes to indicate functional dependence on $\lambda^{\prime}$. We therefore have

[TABLE]

Substituting these expressions in Eq. (208), and using the results we obtained earlier, it is straightforward to see that all the terms in Eqs. (211) and (212) containing oscillating exponentials drop out of $(\phi_{\pm}^{-1}V\phi_{\pm})_{\mathrm{eff},d}$ in the limit $t\to\pm\infty$, leaving

[TABLE]

We can simplify this expression, noting that $[\sigma_{3}Q,[\sigma_{3}Q,\rho_{d}]]=2(q_{o}^{2}\rho_{d}-2\sigma_{3}Q\rho_{d}\sigma_{3}Q)$. After some algebra, we then obtain simply

[TABLE]

Next, since $\rho$ and $\rho_{\pm}$ are traceless, it is $\sigma_{1}\rho_{\pm,d}\sigma_{1}=-\rho_{\pm,d}$. We can therefore combine the two terms in the curly brackets to obtain a simplified version (using the identity $\zeta^{2}-q_{o}^{2}=2\lambda\zeta$). So, we obtain

[TABLE]

Recall that the definition of $w_{\pm}$ is given in Eq. (89). The last step is to simplify Eq. (214). The key idea is to rewrite $w_{\pm}$ as a Hilbert transform, namely

[TABLE]

By using the identity $D_{\pm}\sigma_{3}=\rho_{\pm,d}$, we can combine the two Hilbert transforms in Eq. (214) and notice that $1/\zeta=(\gamma-\lambda)/q_{o}^{2}$. Finally, we obtain the desired result (98). Note that, as $q_{o}\to 0$, Eq. (213) reduces to the result with ZBC. Also, one should notice that this formula holds for all $\zeta\in\mathbb{C}$.

Our next task is to compute Eq. (210). We first look at the terms coming from the $t$-independent parts of $V_{d}$ and $V_{o}$, of which there are four. They yield:

[TABLE]

Now note that $\sigma_{3}Q[\sigma_{3}Q,\rho_{d}]\sigma_{3}Q=[\sigma_{3}Q,\sigma_{3}Q\rho_{d}\sigma_{3}Q]=-q_{o}^{2}[\sigma_{3}Q,\rho_{d}]$. The terms in the curly brackets are then simply $(2i/\zeta\zeta^{\prime})(\lambda^{\prime}-\lambda)[\rho_{\pm,d},\sigma_{3}Q_{\pm}]$. We thus see that $(\phi_{\pm}^{-1}V\phi_{\pm})_{\mathrm{eff},o.1}$ cancels exactly the second term on the RHS of Eq. (207), as expected.

Our last task is to compute the terms in Eq. (210) coming from the $t$-dependent parts of $V_{d}$ and $V_{o}$. There are five terms, which yield:

[TABLE]

where we used that $(\sigma_{3}Q_{\pm})^{2}=q_{o}^{2}$, and where for brevity we denoted $\rho_{\pm,o\sigma_{3}}^{\prime}=\mathrm{e}^{i\gamma^{\prime}\sigma_{3}t}\rho_{\pm,o}^{\prime}\mathrm{e}^{-i\gamma^{\prime}\sigma_{3}t}$. Now note that $[\sigma_{3}Q_{\pm},[\sigma_{3}Q,\rho_{\pm,o\sigma_{3}}^{\prime}]]=2(q_{o}^{2}\rho_{\pm,o\sigma_{3}}^{\prime}-\sigma_{3}Q_{\pm}\rho_{\pm,o\sigma_{3}}^{\prime}\sigma_{3}Q_{\pm})$, and recall $\sigma_{3}Q\,\mathrm{e}^{i\theta\sigma_{3}}=\mathrm{e}^{-i\theta\sigma_{3}}\sigma_{3}Q$. Thus we can rewrite $(\phi_{\pm}^{-1}V\phi_{\pm})_{\mathrm{eff},o.2}$ as

[TABLE]

Recalling Lemma A.1, we see that this integral yields

[TABLE]

Finally, inserting Eq. (215) in Eq. (97) we obtain Eq. (99). Note that the expression for $R_{\pm,o}$ is exactly the same as that for ZBC.

A.5 Propagation of reflection coefficients and norming constants

Here we derive the propagation equations for the reflection coefficients and norming constants.

We start to derive the ODE (106). It is obviously

[TABLE]

Decomposing Eq. (105) into its diagonal and off-diagonal parts, we then have

[TABLE]

We need to: (i) express the RHS in terms of the limiting values as $t\to-\infty$, and (ii) linearize the resulting equations. To accomplish these tasks we need to look at the last three terms on the RHS.

Recall that $\rho^{+}$ is expressed in terms of $\rho^{-}$ via relation (75). Decomposing Eq. (75) into its diagonal and off-diagonal parts (and recalling that $S^{-1}=S^{\dagger}$) we have

[TABLE]

Finally, since $S^{\dagger}S=SS^{\dagger}=I$, it is

[TABLE]

Substituting all of these expressions into ODE (216), and after some tedious but straightforward algebra, and using that $R_{+,o}$ is given by Eq. (99), we finally get the desired ODE (106).

We now derive the propagation equation for the norming constants (110). Recalling definition (109) and taking the derivative (and assuming the limit and the derivative commute), we have

[TABLE]

We compute each of the two derivatives on the RHS separately.

Recall that, for $z\in\Sigma$, the propagation of the scattering matrix is governed by ODE (105). Thanks to Eqs. (103) and (104) some elements of ODE (105) can be extended into the UHP. In particular:

[TABLE]

Also recall from definition (95) that, $\forall z\in\Sigma$,

[TABLE]

Again, thanks to Eqs. (103) and (104) some columns of the above equation can be extended into the UHP. Explicitly,

[TABLE]

Now recall that at $\zeta=\zeta_{n}$ it is $\phi_{+,1}(t,z,\zeta_{n})=b_{n}\phi_{-,2}(t,z,\zeta_{n})$. Thus

[TABLE]

Evaluating the above system in the limit $t\to-\infty$ we see that the terms in square bracket must vanish, yielding

[TABLE]

Substituting Eqs. (218) and (221) into Eq. (217) we then finally obtain the desired ODE (110).

A.6 Solution of the Riemann-Hilbert problem

As usual, in order to solve the RHP we need to take into account its normalization, namely, the asymptotic behavior of $M^{\pm}$ as $\zeta\to\infty$. Recalling the asymptotic behavior of the Jost eigenfunctions and scattering coefficients obtained in Section 3.3, it is easy to check that

[TABLE]

Moreover, from the asymptotic behavior of the Jost eigenfunctions and scattering data as $\zeta\to 0$ we see that

[TABLE]

Thus, as in the case of the defocusing NLS with NZBC, in addition to the behavior at $\zeta=\infty$ and the poles from the discrete spectrum one also needs to subtract the pole at $\zeta=0$ in order to obtain a regular RHP. On the other hand, the asymptotic behavior of the off-diagonal scattering coefficients implies that the jump matrix $G(t,z,\zeta)$ from Eq. (114) is $O(1/\zeta)$ as $\zeta\to\pm\infty$ and $O(\zeta)$ as $\zeta\to 0$ along the real axis.

To solve, we subtract out the asymptotic behavior and the pole contributions. Recall that discrete eigenvalues come in symmetric quartets: $\{\zeta_{n},\zeta_{n}^{*},-q_{o}^{2}/\zeta_{n},-q_{o}^{2}/\zeta_{n}^{*}\}_{n=1}^{N}$. More precisely, $s_{2,2}(\zeta_{n},z)=s_{2,2}(-q_{o}^{2}/\zeta_{n}^{*},z)=s_{1,1}(\zeta_{n}^{*},z)=s_{1,1}(-q_{o}^{2}/\zeta_{n},z)=0$ for $n=1,\dots,N$. We therefore rewrite the jump condition (113) as

[TABLE]

Recall we relabeled all the discrete eigenvalues at the end of Section 3.2 as $\zeta_{N+n}=-q_{o}^{2}/\zeta_{n}^{*}$. Therefore, the above equation reduces to

[TABLE]

Now note that the LHS is analytic in $\Gamma^{-}$ and is $O(1/\zeta)$ as $\zeta\to\infty$ there, while the sum of the first four terms of the RHS is analytic in $\Gamma^{+}$ and is $O(1/\zeta)$ as $\zeta\to\infty$ there. We then introduce modified Cauchy projectors $\mathscr{P}_{\pm}$ by

[TABLE]

where $\mathop{\textstyle\mathop{\textstyle\int}\limits}\limits\nolimits_{\Sigma}$ denotes the integral along the oriented contour shown in Fig. 1, and the notation $z\pm i0$ indicates that when $z\in\Sigma$, the limit is taken from the left/right of it. That is, $F(\zeta\pm i0)=\lim_{\varepsilon\to 0^{+}}F(\zeta\pm i\varepsilon)$. Recall Plemelj’s formulae: if $f_{\pm}$ are analytic in the $\mathbb{C}^{\pm}$ and are $O(1/\zeta)$ as $\zeta$ tends to $\infty$ there, it is

[TABLE]

Applying $\mathscr{P}_{+}$ and $\mathscr{P}_{-}$ to the jump condition (222) we then obtain the solution (115) of the RHP.

A.7 Trace formulae and “theta” condition

Recall that the components of the scattering matrix $s_{1,1}$ and $s_{2,2}$ are analytic respectively in $\Gamma^{+}$ and $\Gamma^{-}$ from Eq. (42). Also recall that the discrete spectrum is composed of quartets: $\zeta_{n},\zeta_{n}^{*},-q_{o}^{2}/\zeta_{n},-q_{o}^{2}/\zeta_{n}^{*}$ $\forall n=1,\dots,N$. Then define the quantities

[TABLE]

The functions $\beta^{\pm}$ are analytic respectively in $\Gamma^{\pm}$, like $s_{1,1}$ and $s_{2,2}$. But, unlike $s_{1,1}$ and $s_{2,2}$, they have no zeros. Moreover, $\beta^{\pm}(\zeta,z)\to 1$ as $\zeta\to\infty$ in the proper half planes. Finally, note that

[TABLE]

We can then take the logarithm of Eq. (225) and apply the Cauchy projectors (see Section A.6), obtaining

[TABLE]

where we used that $\det S=1$, implying $1/(s_{1,1}s_{2,2})=1-b\tilde{b}=1+|b|^{2}$ $\forall\zeta\in\Sigma$ by the definition (39) and symmetries (55) of the reflection coefficients. Now, using definition (224) to replace $\beta^{\pm}$ with $s_{1,1}$ and $s_{2,2}$, we finally obtain the so-called “trace” formulae (123) and (124), which express the analytic scattering coefficients in terms of the discrete eigenvalues and the reflection coefficient. For reflectionless solutions $s_{1,2}(\zeta,z)=0$, $\forall\zeta\in\Sigma$, so the integrals in Eqs. (123) and (124) are identically zero.

Now note that

[TABLE]

because $|b(\zeta^{\prime},z)|=|b(-q_{o}^{2}/\zeta^{\prime},z)|$ thanks to the symmetry (55). A similar relation holds between the integral from $-\infty$ to $-q_{o}$ and that from 0 to $q_{o}$ and for those in the upper/lower semicircle of radius $C_{o}$. However, due to the orientation of the continuous spectrum, these individual integrals do not cancel with each other, but they add together instead. Recalling the asymptotics (71), we then obtain the so-called “theta” condition (125).

A.8 Explicit evaluation of $R_{\pm,d}$

In this appendix, we evaluate the integral (98) that defines $R_{\pm,d}$. Recall that the auxiliary matrix $R_{\pm}$ has the same value when the Hilbert transform is computed on either $\lambda^{\prime}$-sheet, so we will compute this matrix on sheet I, i.e., we choose $\nu=1$ and $\rho_{-}(\lambda,z)=h_{-}\,\sigma_{3}$ from Eqs. (77), (78), (79) and (80). We start by rewriting the integral (98) for convenience. Also recall from the main text, we know in inhomogeneous broadening the auxiliary quantity $w_{-}=0$ in Eq. (98).

Firstly, we consider the case where $\lambda\in\mathbb{C}\backslash\mathbb{R}$.

Recall that the Hilbert transform is given by Eq. (8). There is no singularity on the integration contour, and the principal value is not needed. To compute the integral, it is convenient to write things in terms of the uniformization variable $\zeta$. Recall that $\lambda=(\zeta-q_{o}^{2}/\zeta)/2$ and $\gamma=(\zeta+q_{o}^{2}/\zeta)/2$ from Eq. (14). Letting $s=\zeta^{\prime}$ it is $d\lambda^{\prime}={\textstyle\frac{1}{2}}(1+q_{o}^{2}/s^{2})\,\mathrm{d}s$ and, taking the principal branch of the square root, the integration contour is $L_{o}=(-\infty,-q_{o})\cup(q_{o},\infty)$. Moreover,

[TABLE]

where

[TABLE]

Collecting all the pieces,

[TABLE]

The two roots of $\Delta_{2}(s)$ are

[TABLE]

while the four roots of $\Delta_{4}(s)$ are

[TABLE]

Let $s_{1},\dots,s_{6}$ denote the zeros of $\Delta_{2}(s)$ and $\Delta_{4}(s)$, and let $f(s)$ denote the numerator of the integrand in Eq. (226). Expanding the integrand in partial fractions we have:

[TABLE]

where the prime indicates that the term $m=j$ is omitted from the product. Hence

[TABLE]

After some tedious but straightforward algebra, we then get Eq. (152).

Secondly, we consider the case where $\lambda\in\mathbb{R}$.

In this case, there is a singularity in the integral and the principal value becomes necessary. Then we use the following relationship to compute the integral in Eq. (98).

[TABLE]

where the contour $L$ is: $L=(-\infty,\lambda-r)\cup\{z=r\,\mathrm{e}^{i\theta}|\,\theta\in[0,\pi]\}\cup(\lambda+r,\infty)$, provided $a$ is sufficiently small and there are no singularities on the contour. We will also use $\lambda^{\prime}$ on sheet I. The integral on the RHS of Eq. (227) can be computed normally, because $\lambda$ is not on the contour. So this integral is exactly the same as the one we computed in the first case where $\lambda\in\mathbb{C}\backslash\mathbb{R}$. Thus,

[TABLE]

The second term on the RHS of Eq. (227) can be computed easily,

[TABLE]

Therefore, by using the relationship (227), we obtain Eq. (154) for $\lambda\in\mathbb{R}$.

A.9 IST with branch cut outside

In this appendix, we will formulate an alternative version of IST to solve the same problem (3) and (4). The idea of the new version is to use a different branch cut in the complex $\lambda$-plane. As shown later, this version of IST is equivalent to the one presented in the main text.

Recall that the eigenvalue $\gamma(\lambda)$ was defined in the main text with $\gamma(\lambda)=\mathop{\rm sign}\nolimits(\lambda)\sqrt{q_{o}^{2}+\lambda^{2}}$ for $\lambda\in\mathbb{R}$. Alternatively, one could define it as $\gamma(\lambda)=\sqrt{q_{o}^{2}+\lambda^{2}}$ for $\lambda\in\mathbb{R}$ (i.e., without the signum function), and then formulate the IST for MBEs in a similar way.

Background states.

As before, we first investigate whether there exist exact “constant” solutions of the MBEs (3) and (4) with NZBC. So suppose for now that $q(t,z)=q_{o}(z)$ $\forall t\in\mathbb{R}$, with $q_{o}\neq 0$. Equation (9) is still the solution of Eq. (3), i.e., $\rho(t,z,\lambda)=\mathrm{e}^{X_{o}t}\,C\,\mathrm{e}^{-X_{o}t}$. The eigenvalues of $X_{o}$ are still $\pm i\gamma$, where $\gamma^{2}=\lambda^{2}+q_{o}^{2}$. Consistently with the above discussion, we now take

[TABLE]

With this choice, $\gamma$ is always continuous along the continuous spectrum. The branching structure of $\gamma(\lambda)$ will be discussed later. We can write an eigenvector matrix of $X_{o}$ as Eq. (11), i.e., $Y_{o}=I+(i/\zeta)\sigma_{3}Q_{0}$, where we use the short notation (12) again. Thus, we again have Eq. (13), i.e., $\rho(t,z,\lambda)=Y_{o}\,\mathrm{e}^{i\gamma t\sigma_{3}}\rho_{o}\,\mathrm{e}^{-i\gamma t\sigma_{3}}Y_{o}^{-1}$, where $\rho_{o}=Y_{o}^{-1}CY_{o}$. We again express all $\lambda$ dependence in terms of $\zeta$, noting that the inverse transformation to Eq. (12) is the same as Eq. (14).

Using similar arguments as in Section 2.2, we can again write $\rho(t,z,\zeta)=\mathbf{h}\cdot{\boldsymbol{\rho}}$, where ${\boldsymbol{\rho}}=(\rho_{1},\rho_{2},\rho_{3})^{T}$, $\mathbf{h}\cdot\mathbf{h}=1$ and $\rho_{j}(t,z,\zeta)$ is still given by Eq. (15), i.e., $\rho_{j}(t,z,\zeta)=Y_{o}\,\mathrm{e}^{i\gamma t\sigma_{3}}\sigma_{j}\,\mathrm{e}^{-i\gamma t\sigma_{3}}Y_{o}^{-1}$, with $\rho_{j}^{\dagger}=\rho_{j}$, ${\mathrm{tr}}\rho_{j}=0$, and $\det\rho_{j}=-1$. Explicitly, we still have formulas (16), (17) and (18).

Now we insert this behavior into Eq. (4) as before. Direct calculations yield Eqs. (19), (20) and (21) again. By a similar discussion as in the main text, we know that only the commutator in Eq. (21) is $t$-independent. Moreover, we have $\mathop{\textstyle\mathop{\textstyle\int}\limits}\limits[\sigma_{3},\rho_{3}]g(\lambda)\mathrm{d}\lambda=-iw_{o}[\sigma_{3},Q_{0}]$, where $w_{o}=\mathop{\textstyle\mathop{\textstyle\int}\limits}\limits g(\lambda)/\gamma\,\mathrm{d}\lambda$. Then, it is easy to compute the corresponding integral

[TABLE]

with the choice of the detuning function $g(\lambda)$ as in Eq. (151). Notice that due to different definitions of $\gamma$, this integral differs from the one in the main text. It follows that the only self-consistent solutions are

[TABLE]

These expressions differ from those in Eq. (24) for all $z>0$.

Formulation of the IST.

As before, we define the quantity

[TABLE]

Now however we take the branch cut on $i(-\infty,-q_{o}]\cup i[q_{o},\infty)$. This is done by taking $\lambda+iq_{o}=r_{1}\mathrm{e}^{i\theta_{1}}$ and $\lambda-iq_{o}=r_{2}\mathrm{e}^{i\theta_{2}}$, so that

[TABLE]

where $m=0,1$ indicates the first or the second sheet, respectively, and $-\pi/2\leqslant\theta_{1}<3\pi/2$, $-3\pi/2<\theta_{2}\leqslant\pi/2$. It is easily verified that with this choice the discontinuity of $\gamma$ occurs on the segment $i(-\infty,-q_{o}]\cup[q_{o},\infty)$. As before, the Riemann surface is obtained by gluing the two sheets of the complex $\lambda$-plane along the cut by defining the uniformization variable $\zeta=\lambda+\gamma$, which maps the first sheet into the right half plane and the second one into the left half plane. The continuous spectrum is the same as before, i.e., $\lambda\in\mathbb{R}\cup i[-q_{o},q_{o}]$. On the $\zeta$-plane, the continuous spectrum is also the same, $\zeta\in\mathbb{R}\cup C_{o}$. The left half of the first sheet and the right half of the second sheet are mapped into the interior of $C_{o}$. The branch cut becomes the imaginary axis on the $\zeta$-plane. Moreover, the inverse transformations are the same, i.e., Eq. (14).

Differently from before, however, for $\lambda\in\mathbb{R}$ we have $\gamma=\sqrt{q_{o}^{2}+\lambda^{2}}$, which has no discontinuity at $\lambda=0$. If one now takes the limit $q_{o}\to 0$, it is obvious that $\gamma=\pm|\lambda|$ for $\lambda\in\mathbb{R}$ (where the $\pm$ values apply on sheet I or II respectively). Thus, this formulation of the IST does not reduce to the one with ZBCs in the limit $q_{o}\to 0$ directly. In fact, we will show below that this formulation of the IST yields different solutions of the MBEs.

Even though the new $\gamma$ is different from the one in the main text, for $\lambda\in\mathbb{C}$, $\gamma(\lambda)$ is still given by definition (10). This means that the calculations in the two formulations of the IST are very similar. In fact, the direct and inverse problem are identical in both versions. This is because all of the IST was formulated using the uniformization variable $\zeta=\lambda+\gamma(\lambda)$, and the values of $\gamma(\lambda)$ for $\zeta\in\mathbb{C}$ are exactly the same in the two formulations. Essentially, the definition of $\gamma(\lambda)$ amounts to just switching the left half planes between sheet I and sheet II.

In other words, we can use the same definition of the regions $\Gamma^{\pm}$ of the complex $\zeta$-plane from definition (27). It is then easy to show that the analyticity properties of the columns of the eigenfunctions $\mu_{\pm}$ in the complex $\zeta$-plane are still the same, namely Eq. (36). The analyticity properties of the scattering data and the symmetries for eigenfunctions and scattering data are also the same as before.

The only part of the IST that is affected by the change is the propagation, because it involves integrals over $\lambda\in\mathbb{R}$. It is easy to repeat the calculation in Section 3.4 and 4.2 and obtain the integrals for the boundary data $w_{\pm}$ and $R_{\pm,d}$. The resulting integrals are formally the same as Eqs. (89) and (98). In other words, as long as the explicit expression of $\gamma$ is not used, the formulas are the same.

Equivalence between two versions of IST.

In this part, we will prove that the two versions of IST are equivalent, in the sense that one could obtain the same solutions simultaneously by using different boundary conditions in two versions of IST.

To begin with, it is convenient to introduce a new notation in this proof. The superscript ${\mathrm{I}}$ or ${\mathrm{II}}$ denotes the variable or function in the first or second version of IST, respectively. For example, $\gamma^{\mathrm{I}}$ denotes $\gamma$ in the first IST (i.e., the IST presented in the main text) and $\rho_{-}^{{\mathrm{II}}}$ denotes the boundary condition $\rho_{-}$ in the second version (i.e., the IST discussed in this appendix).

Firstly, one should notice that the solutions of MBEs are unique. More precisely, the solutions are uniquely determined by the input pulse $Q(t,0)$, together with the initial state $\rho_{-}(\zeta,z)$ which is defined by Eq. (74). Thus, to prove the equivalence between the two versions of IST, it is sufficient to find both sets of data in the two versions of IST, respectively, that produce the same solution.

Secondly, one should also notice that, both quantities $Q(t,0)$ and $\rho_{-}(\zeta,z)$ are defined with $\lambda\in\mathbb{R}$. So we can focus on real values of $\lambda$ (or $\gamma$, $\zeta$) instead of complex values. Moreover, because the fundamental difference between the two versions of IST comes from the definition of $\gamma$ with negative values of $\lambda$ (or equivalently, how to take the branch cut), it is sufficient to only consider the case where $\lambda<0$.

Now, we are ready to begin our proof. For any given solution $Q(t,z)$ and $\rho(t,z,\lambda)$, it is obvious that the corresponding inputs in the two versions of IST must be the same. [Recall that the input pulse is $Q(t,0)$.] It is then sufficient to consider the initial data $\rho_{-}^{\mathrm{I}}(\zeta^{\mathrm{I}},z)$ and $\rho_{-}^{{\mathrm{II}}}(\zeta^{{\mathrm{II}}},z)$. Let $q_{-}=\lim_{t\to-\infty}q(t,z)$ as before, and let the density matrix be

[TABLE]

where $\sigma_{j}$ are Pauli matrices from Eq. (5) and $h_{j}\in\mathbb{R}$ for $j=1,2,3$. By using relation (74), one can compute the corresponding initial states $\rho_{-}^{\mathrm{I}}$ and $\rho_{-}^{\mathrm{II}}$. Moreover, since both $\rho_{-}^{\mathrm{I}}$ and $\rho_{-}^{\mathrm{II}}$ are Hermitian matrices, it is sufficient to compute the $(1,1)$ and $(1,2)$ components. Consequently, the elements are (where $\gamma_{o}=\sqrt{q_{o}^{2}+\lambda^{2}}$ for $\lambda<0$)

[TABLE]

Thus, the two sets of initial states for the medium are found.

By comparing the above expressions (229), (230), (231) and (232), we obtain the following relationship between the two sets of data that produce the same solution with $\lambda<0$

[TABLE]

or the reversion

[TABLE]

where subscript $d$ or $o$ denotes the diagonal or off-diagonal part of the matrix, respectively, as before. Recall that for positive values of $\lambda$, the two sets of data are exactly the same. Therefore, if one solution is produced by either one of the versions of IST, the relationships (233) and (234) ensure that the same solution can also be produced by using the other version. Thus the two versions of IST have the same solutions set and the equivalence is proved.