Topologically protected strongly correlated states of photons

Light–matter interaction surrounds us everywhere in nature, and it has played a tremendous role in the development of current technology. Until recent decades it was sufficient to deal with this interaction on average, with many photons and many atoms involved. However, the increasing miniaturization of basic constituents toward nanoscale is a common trend in modern technology. Downscaling to single-atom and/or single-photon levels promotes some traditionally classical research areas into the quantum realm [1–4]. A quantum control over the light–matter interaction will eventually become a vital ingredient of emerging quantum devices, and it is equally important for the development of several related fields, including communication, signal processing, ultrafast optics, optomechanical cooling, imaging and spectroscopy, and quantum information [5]. However, efficient manipulation and control requires a relatively strong interaction at the level of a single atom or a single photon [6–12]. This presents a significant challenge, since the typical interaction scale for individual particles is given by the smallness of the QED coupling constant $\alpha \approx 1/137$. Two possible ways to overcome this natural limitation, and to increase the effects of correlations, are either to use artificial materials and devices, or to resort to many-body effects to produce nonlinearities.

Recent experimental progress in fabricating few-photon sources coupled to one-dimensional (1D) transmission lines [13–20] opens an avenue for creating and manipulating strongly correlated states of photons. It has also triggered a significant number of theoretical studies [21–27]. Our experience from condensed matter and atomic physics teaches us that combined effects of reduced dimensionality and interactions can often effectively enhance correlations, and eventually lead to new collective states of matter with properties which are very different from those of the individual particles (e.g., Luttinger liquids in 1D). This insight is one of the driving forces behind the quest for novel correlated states of photons [28–30].

A different class of collective phenomena is the topological insulating/superconducting state of matter. This state has been recently experimentally realized for electrons [31] and extensively studied theoretically [32]. The main signature of topological properties in the band structure of bulk materials is the existence of edge states, which are insensitive to local perturbations, impurities, or geometrical imperfections. In some—chiral—instances edge states propagate unidirectionally. Recently, the existence of chiral states in photonic systems has been theoretically predicted [33] and experimentally observed [34–39]. These edge states represent an optical analogue of the quantum Hall edge states, and they are characterized by a nontrivial Chern number [40].

In this paper, we suggest the use of 1D unidirectional edge states as a robust platform for the controllable generation of strongly correlated states of photons. To achieve this goal, we couple an array of emitters to a chiral edge state, see figure 1. Experimentally available novel hybrid photonic systems [41] can be used to realize this setup. Photons in the edge channel, populated by an external few-photon source, interact with an ensemble of emitters and produce outgoing states with robust, controllable, and universal properties originating from the topological nature of the edge state. In particular, we find that the outgoing photonic wavefunction does not contain any information about the positions of emitters; its nodes are rather determined by universal numbers—the zeroes of Laguerre polynomials. An initial single-photon wavepacket is fragmented into pieces between these nodes. In the case of two-photon scattering, we observe a clearly pronounced parity effect with respect to the number of emitters, which manifests in a transition from photons' bunching to antibunching as one changes the parity of the number of emitters from even to odd. We also show that the observed properties are robust with respect to fluctuations in the emitters' transition frequencies. The proposed setup can be experimentally realized in the GHz and optical domains with existing nanophotonic elements.

Zoom In Zoom Out Reset image size

We consider topologically protected states propagating unidirectionally at the edge of the topological photonic crystal. We propose to embed an array of two-level (in general—multi-level) systems—called emitters—into the chiral edge channel. We assume that the transition frequencies of emitters are commensurate with the frequency of the propagating chiral mode of light. Assuming also that the excited emitter states mainly decay into the 1D mode (with a decay rate ${{\Gamma }_{1{\rm d}}}$), we note that there are several sources of inevitable losses. Namely, the excited emitter can decay into the continuum of three-dimensional modes in ambient space (with a decay rate Γ₀), into the bulk modes of a 2D photonic crystal (decay rate ${{\Gamma }_{2{\rm d}}}$), or to its impurity (bound) states (decay rate ${{\Gamma }_{{\rm bs}}}$). Here we note that the decay into the two-dimensional bulk modes of a photonic crystal is suppressed, ${{\Gamma }_{2{\rm d}}}=0$, since their density of states is zero in the bulk bandgap. We also assume that a 2D photonic crystal is clean enough to neglect the coupling of edge states to eventual mid-gap impurity states.

To generate well-defined chiral modes, it is necessary to minimize losses occurring at the rate Γ₀. The waveguide Purcell factor [42]

$$\begin{eqnarray}&&\frac{{{\Gamma }_{1{\rm d}}}}{{{\Gamma }_{0}}}=\frac{3}{\pi \epsilon _{{\rm d}}^{3/2}}\frac{c}{{{v}_{{\rm g}}}}\frac{{{A}_{0}}}{{{A}_{{\rm eff}}}}\end{eqnarray} \tag{ 1 }$$

must then considerably exceed unity, which defines the strong coupling regime. In this expression ${{A}_{0}}=\lambda _{0}^{2}/4$ is the minimal cross-sectional area to confine the light in vacuum (for the lightʼs wavelength λ₀); _d is the dielectric constant of the host medium where the emitter is embedded; v_g is the group velocity of the propagating edge mode. The cross-sectional area ${{A}_{{\rm eff}}}$ of the effective confinement is estimated by $\sim b\xi$, where b is the thickness of the slab, and $\xi \sim c\hbar /B$ is the localization depth of the edge mode, B being the bulk bandgap. Thus, the waveguide Purcell factor (1) can be enhanced in three different ways: (1) by reducing the thickness b of the slab; (2) by increasing the bandgap B; and (3) by reducing the group velocity v_g. All these methods to achieve the strong coupling regime have been successfully applied in plasmonic nanowires [24, 43, 44] and in the Floquet photonic topological insulators made of helical waveguides [39].

Deep inside the bulk bandgap we can model the dispersion of the chiral edge mode by a linear dependence. In the following we measure all gauge-dependent energy scales from the spectrum linearization point. In the strong coupling regime we neglect losses into the three-dimensional continuum, and consider only the interaction of the edge mode with an array of M emitters. Emitters are modeled by two-level systems with transition frequencies Δ_a and couplings $\sqrt{{{\kappa }_{a}}}\sim \sqrt{{{\Gamma }_{1{\rm d}}}}$, and are placed at different positions x_a, which are separated from each by distances larger than the lightʼs wavelength (an arrangement opposite to Dickeʼs [45]). Thus, the Hamiltonian of our model reads

$$\begin{eqnarray}\begin{array}{rcl} H & = & -i\int {\rm d}x{{a}^{\dagger }}(x){{\partial }_{x}}a(x)+\mathop{\sum }\limits_{a=1}^{M}{{\Delta }_{a}}\left( S_{a}^{z}+\frac{1}{2} \right) \\ {} & {} & -\mathop{\sum }\limits_{a=1}^{M}\sqrt{{{\kappa }_{a}}}[S_{a}^{+}a({{x}_{a}})+{{a}^{\dagger }}({{x}_{a}})S_{a}^{-}], \\ \end{array}\end{eqnarray} \tag{ 2 }$$

where we use units such that ${{v}_{{\rm g}}}=\hbar =1$. The chiral field operators ${{a}^{\dagger }}(x),a(x)$ satisfy the standard commutation relations $[a(x),{{a}^{\dagger }}(x^{\prime} )]=\delta (x-x^{\prime} )$. Transitions between emitters' states are described by the operators $S_{a}^{\pm }$, which satisfy the standard spin algebra $[S_{a}^{z},S_{a^{\prime} }^{\pm }]=\pm S_{a}^{\pm }{{\delta }_{aa^{\prime} }}$, $[S_{a}^{+},S_{a^{\prime} }^{-}]=2S_{a}^{z}{{\delta }_{aa^{\prime} }}$. Spontaneous emission to other modes out of the 1D waveguide is modeled by attributing an imaginary part $-{\rm i}\kappa ^{\prime} /2\equiv -i({{\Gamma }_{0}}+{{\Gamma }_{{\rm bs}}})/2$ to the transition frequencies ${{\Delta }_{a}}\to {{\Delta }_{a}}-{\rm i}\kappa ^{\prime} /2$, in the spirit of the quantum jump picture [46].

The problem is further specified by an initial state. Inspired by experimental realizations of waveguides coupled to emitters [18, 19], we assume that the edge states are populated by an external few-photon source, while emitters are initially prepared in the ground state. Injected photons propagate unidirectionally and interact with an ensemble of emitters, and after a sufficiently long time a stationary state is eventually established: emitters typically relax back to the ground state, while the photonic wavefunction is modified. The evolution of the N-photon wavefunction is described in terms of the scattering matrix $S_{M;\{{{k}_{i}}\},\{k_{i}^{\prime }\}}^{(N)}\equiv S_{M;\{{{k}_{i}}\},\{k_{i}^{\prime }\}}^{(N)}(\{{{\Delta }_{a}}\},\{{{\kappa }_{a}}\})$ depending on the sets on all transition frequencies Δ_a and couplings κ_a

$$\begin{eqnarray}&&{{\phi }_{N,{\rm out}}}(\{{{k}_{i}}\})=\int {\rm d}\{k_{i}^{\prime }\}S_{M;\{{{k}_{i}}\},\{k_{i}^{\prime }\}}^{(N)}{{\phi }_{N,{\rm in}}}(\{k_{i}^{\prime }\}),\end{eqnarray} \tag{ 3 }$$

where ${{\phi }_{N,{\rm in}}}(\{k_{i}^{\prime }\})$ and ${{\phi }_{N,{\rm out}}}(\{{{k}_{i}}\})$ are envelope functions of outgoing and incoming states depending on sets of outgoing $\{{{k}_{i}}\}$ and incoming $\{k_{i}^{\prime }\}$ photonic energies, respectively, obeying the energy conservation $\sum _{i=1}^{N}{{k}_{i}}=\sum _{i=1}^{N}k_{i}^{\prime }$. The convolution in (3) is performed over the incoming set of momenta ${\rm d}\{k{{^{\prime} }_{i}}\}=\frac{1}{N!}\prod _{i=1}^{N}{\rm d}k_{i}^{\prime }.$

A general diagrammatic approach to calculate the scattering matrix of a local emitter with an arbitrary level structure and transition amplitudes has been developed in [47]. It gives results coinciding with a direct solution of the Lippmann–Schwinger equation [22].

A theoretical study of the scattering off an array of distributed scatterers is more involved, since one has to take into account interference effects. For a single photon scattering, an evaluation of the scattering matrix is facilitated by an application of the transfer matrix method [25, 48]. For the scattering of two or more photons, there is no general prescription on how to compute the exact scattering matrix, and the complexity of this task is determined by the interplay of interference and correlation effects. There are few numerical results in the literature [49, 50] on this problem.

The model (2), however, affords a considerable simplification based on the unidirectional propagation of light: no backscattering can happen during each scattering event. Moreover, all photons travel with the same group velocity. For this reason the interference does not occur, and the result of scattering does not depend on the travel time between emitters. Therefore, the scattering from one emitter happens independently of any other emitter, and the net result of the scattering off an array of emitters is represented by the convolution property [47]

$$\begin{eqnarray}\begin{array}{rcl} S_{M;\{{{k}_{{{i}_{M}}}}\},\{{{k}_{{{i}_{0}}}}\}}^{(N)}(\{{{\Delta }_{a}}\},\{{{\kappa }_{a}}\}) & = & \int \left( \mathop{\prod }\limits_{b=1}^{M-1}{\rm d}\{k_{{{i}_{b}}}^{{}}\} \right)S_{\{{{k}_{{{i}_{M}}}}\},\{{{k}_{{{i}_{M-1}}}}\}}^{(N)}({{\Delta }_{M}},{{\kappa }_{M}}) \\ {} & {} & \times S_{\{{{k}_{{{i}_{M-1}}}}\},\{{{k}_{{{i}_{M-2}}}}\}}^{(N)}({{\Delta }_{M-1}},{{\kappa }_{M-1}})\ldots \\ {} & {} & \times S_{\{{{k}_{{{i}_{2}}}}\},\{{{k}_{{{i}_{1}}}}\}}^{(N)}({{\Delta }_{2}},{{\kappa }_{2}})S_{\{{{k}_{{{i}_{1}}}}\},\{{{k}_{{{i}_{0}}}}\}}^{(N)}({{\Delta }_{1}},{{\kappa }_{1}}), \\ \end{array}\end{eqnarray} \tag{ 4 }$$

where ${{S}^{(N)}}({{\Delta }_{a}},{{\kappa }_{a}})$ is the N-photon scattering matrix of the ath emitter.

The property (4) is very basic. We summarize the condition under which it holds: (1) unidirectional nature of the spectrum; (2) a constant group velocity of the incoming wavepacket; and (3) linear and energy-independent coupling between the photons and the emitters. The independence of (4) of emitters' positions x_a lies at the origin of many universal properties of the outgoing photonic states which we discuss in the following.

A combination of methods to evaluate $S_{M=1}^{(N)}$ of a single emitter [22, 47] with the convolution property (4) provides a general platform for calculating scattering outcomes off arrays composed of emitters with an arbitrary complex structure of levels and transitions between them. For example, one can use three-level emitters with ${{S}_{+}}={{g}_{31}}\mid 3\rangle \langle 1\mid +{{g}_{32}}\mid 3\rangle \langle 2\mid$ (Λ-scheme), ${{S}_{+}}={{g}_{31}}|3\rangle \langle 1|+{{g}_{21}}|2\rangle \langle 1|$ (V-scheme), or ${{S}_{+}}={{g}_{32}}|3\rangle \langle 2|+{{g}_{21}}|2\rangle \langle 1|$ (Σ-scheme). One can even combine emitters of different types along the line of lightʼs propagation. In all such cases the scattering matrix $S_{M}^{(N)}$ can be explicitly determined. Once we understand how its properties depend on the emitters' parameters, we obtain a powerful tool to engineer correlated multiphoton states.

In this paper, we concentrate on the description of the model (2). We also remark that its alternative solution in the case of identical couplings ${{\kappa }_{a}}=\kappa$ was obtained by the means of the Bethe ansatz [51–54], which we use as a benchmark to verify our approach based on the usage of the convolution property (4).

2.1. Single-photon scattering

Let us now discuss an application of the general theory formulated above to the scattering of few-photon wavepackets off an array of M emitters.

We start from the most basic case of single-photon scattering. In the following we will neglect losses setting $\kappa ^{\prime} =0$. This scattering is purely elastic: a photon scattered off a single emitter with parameters Δ_a and κ_a just acquires an additional phase, which defines the scattering matrix

$$\begin{eqnarray}&&S_{kk^{\prime} }^{(1)}={{\delta }_{kk^{\prime} }}\frac{k-{{\Delta }_{a}}-{\rm i}{{\kappa }_{a}}/2}{k-{{\Delta }_{a}}+{\rm i}{{\kappa }_{a}}/2}\equiv {{\delta }_{kk^{\prime} }}{{{\rm e}}^{{\rm i}\varphi _{k}^{(a)}}}.\end{eqnarray} \tag{ 5 }$$

Furthermore, phases acquired on every individual scatterer are additive, in accordance with the convolution property (4). This gives the single-photon scattering matrix of an array of M emitters

$$\begin{eqnarray}&&S_{M;kk^{\prime} }^{(1)}={{\delta }_{kk^{\prime} }}\mathop{\prod }\limits_{a=1}^{M}{{\frac{k-{{\Delta }_{a}}-{\rm i}{{\kappa }_{a}}/2}{k-{{\Delta }_{a}}+{\rm i}{{\kappa }_{a}}/2}}^{{}}}\equiv \delta _{kk^{\prime} }^{{}}{{{\rm e}}^{{\rm i}\mathop{\sum }\limits_{a=1}^{M}\varphi _{k}^{(a)}}}.\end{eqnarray} \tag{ 6 }$$

The corresponding outgoing single-photon wavefunction at point x and time t can be decomposed as

$$\begin{eqnarray}\begin{array}{rcl} {{\phi }_{1,{\rm out}}}(x-t) & \equiv & \int {\rm d}x^{\prime} S_{M}^{(1)}(x-t,x^{\prime} ){{\phi }_{1,{\rm in}}}(x^{\prime} ) \\ {} & = & {{\phi }_{1,{\rm in}}}(x-t)+{{\phi }_{1,{\rm scatt}}}(x-t), \\ \end{array}\end{eqnarray} \tag{ 7 }$$

where $S_{M}^{(1)}(x,x^{\prime} )$ is the coordinate representation of (6), and ${{\phi }_{1,{\rm scatt}}}$ is the scattered part of the outgoing wavefunction.

To understand the structure of ${{\phi }_{1,{\rm scatt}}}$ we consider the limiting case of identical emitters ${{\Delta }_{a}}\to \Delta$, ${{\kappa }_{a}}\to \kappa$, and obtain the expression (see the appendix for details)

$$\begin{eqnarray}&&{{\phi }_{1,{\rm scatt}}}(x-t)=-\int {\rm d}x^{\prime} \Theta (x^{\prime} )\kappa L_{M-1}^{(1)}(\kappa x^{\prime} ){{{\rm e}}^{-({\rm i}\Delta +\kappa /2)x^{\prime} }}{{\phi }_{1,{\rm in}}}(x^{\prime} +x-t),\end{eqnarray} \tag{ 8 }$$

where $L_{M-1}^{(1)}(x)$ is the associated Laguerre polynomial. An appearance of the polynomial behavior is remarkable, and we next discuss its implications.

For a realistic scattering experiment we specify the initial wavepacket ${{\phi }_{1,{\rm in}}}(x)=\frac{1}{\sqrt{\sigma \sqrt{\pi }}}{{{\rm e}}^{{\rm i}(\Delta +\delta )x-\frac{{{x}^{2}}}{2{{\sigma }^{2}}}}}$. In the momentum space it corresponds to a Gaussian distribution around $k^{\prime} =\Delta +\delta$ with the variance $1/\sigma$, where $\delta$ is the detuning. Assuming that ${{\delta }^{-1}},{{\kappa }^{-1}}\gg \sigma$, we obtain

$$\begin{eqnarray}&&{{\phi }_{1,{\rm scatt}}}(x-t)\approx -\sqrt{2\sigma \sqrt{\pi }}\Theta (t-x)\kappa L_{M-1}^{(1)}(\kappa (t-x)){{{\rm e}}^{-\kappa (t-x)/2}}{{{\rm e}}^{{\rm i}\Delta (x-t)}}.\end{eqnarray} \tag{ 9 }$$

In figure 2 we plot $|{{\phi }_{1,{\rm scatt}}}{{|}^{2}}$ for various numbers of emitters M and detunings $\delta$.⁵ These results clearly manifest the universal character of scattering in the system under consideration: the dependence of the outgoing wavepacket on the positions of emitters is absent, while the minima are determined by universal numbers—the zeroes of the Laguerre polynomials. Thus, if the position of the first node is known, the subsequent nodes can be determined from $L_{M-1}^{(1)}(x)$. An emergence of nodes is accounted by time delays on each emitter, which eventually leads to the fragmentation of the incoming wavepacket into $M-1$ pieces.

Zoom In Zoom Out Reset image size

2.2. Two-photon scattering

Let us next consider two-photon scattering. It can happen in two different ways: (1) via the elastic scattering of two individual photons; (2) via the inelastic scattering of two photons exchanging energy with each other. The existence of the second possibility is characteristic for interacting systems, which leads to an emergence of correlated states of photons. The corresponding two-photon scattering matrix is represented by a sum of reducible (elastic) and irreducible (inelastic) terms

$$\begin{eqnarray}&&S_{{{k}_{1}}{{k}_{2}},k_{1}^{\prime }k_{2}^{\prime }}^{(2)}=S_{{{k}_{1}}k_{1}^{\prime }}^{(1)}S_{{{k}_{2}}k_{2}^{\prime }}^{(1)}+S_{{{k}_{1}}k_{2}^{\prime }}^{(1)}S_{{{k}_{2}}k_{1}^{\prime }}^{(1)}+i\mathcal{T}_{{{k}_{1}}{{k}_{2}},k_{1}^{\prime }k_{2}^{\prime }}^{(2)}.\end{eqnarray} \tag{ 10 }$$

This picture of scattering mechanisms is important for the interpretation of scattering results below.

The amplitude of inelastic scattering off a single emitter M = 1 is known [22, 47]

$$\begin{eqnarray}&&\mathcal{T}_{{{k}_{1}}{{k}_{2}},k_{1}^{\prime }k_{2}^{\prime }}^{(2)}=\frac{{{\kappa }^{2}}}{\pi }\frac{K-2\Delta +{\rm i}\kappa }{\mathop{\prod }\limits_{p}{{(p-\Delta +{\rm i}\kappa /2)}^{{}}}}{{\delta }_{{{k}_{1}}+{{k}_{2}},k_{1}^{\prime }+k_{2}^{\prime }}},\end{eqnarray} \tag{ 11 }$$

where $K={{k}_{1}}+{{k}_{2}}=k_{1}^{\prime }+k_{2}^{\prime }$ and where $p\equiv ({{k}_{1}},{{k}_{2}},k_{1}^{\prime },k_{2}^{\prime })$. On the basis of this expression one can explain the bunching property of photons in 1D waveguides. Most clearly this can be viewed in the coordinate representation, see below.

In order to find an explicit expression of $S_{M;{{k}_{1}}{{k}_{2}},k_{1}^{\prime }k_{2}^{\prime }}^{(2)}$ for arbitrary M it is necessary to compute $M-1$ convolutions in (4), each being represented by a twofold momentum integration. As both $S_{kk^{\prime} }^{(1)}$ in (5) and $\mathcal{T}_{{{k}_{1}}{{k}_{2}},k_{1}^{\prime }k_{2}^{\prime }}^{(2)}$ in (11) have a simple pole structure, this becomes a routine task.

A direct pathway to $S_{M}^{(2)}$ is available in the case when the couplings of the field to all emitters are identical, ${{\kappa }_{a}}=\kappa$. The two-photon (and also N-photon) scattering matrix in the coordinate representation can be then determined by the Bethe ansatz method, and it reads [53]

$$\begin{eqnarray}\begin{array}{rcl} S_{M}^{(2)}({{y}_{1}},{{y}_{2}},{{z}_{1}},{{z}_{2}}) & = & \frac{1}{{{(2\pi )}^{2}}}\int {\rm d}{{k}_{1}}{\rm d}{{k}_{2}}{\rm d}k_{1}^{\prime }{\rm d}k_{2}^{\prime }{{{\rm e}}^{{\rm i}{{k}_{1}}{{y}_{1}}+{\rm i}{{k}_{2}}{{y}_{2}}-{\rm i}k_{1}^{\prime }{{z}_{1}}-{\rm i}k_{2}^{\prime }{{z}_{2}}}}S_{M;{{k}_{1}}{{k}_{2}},k_{1}^{\prime }k_{2}^{\prime }}^{(2)} \\ {} & = & \Theta ({{z}_{1}}-{{z}_{2}}){{\int }_{{{\gamma }_{1}}}}\frac{{\rm d}{{\lambda }_{1}}}{2\pi }{{\int }_{{{\gamma }_{2}}}}\frac{{\rm d}{{\lambda }_{2}}}{2\pi }{{{\rm e}}^{{\rm i}{{\lambda }_{1}}({{y}_{1}}-{{z}_{1}})+{\rm i}{{\lambda }_{2}}({{y}_{2}}-{{z}_{2}})}} \\ {} & {} & \times \left[ 1-\frac{2{\rm i}\kappa \Theta ({{y}_{2}}-{{y}_{1}})}{{{\lambda }_{1}}-{{\lambda }_{2}}+{\rm i}\kappa } \right] \\ {} & {} & \times \left( \mathop{\prod }\limits_{a=1}^{M}\frac{{{\lambda }_{1}}-{{\Delta }_{a}}-{\rm i}\kappa /2}{{{\lambda }_{1}}-{{\Delta }_{a}}+{\rm i}\kappa /2} \right)\left( \mathop{\prod }\limits_{b=1}^{M}\frac{{{\lambda }_{2}}-{{\Delta }_{b}}-{\rm i}\kappa /2}{{{\lambda }_{2}}-{{\Delta }_{b}}+{\rm i}\kappa /2} \right) \\ {} & {} & +({{y}_{1}}\leftrightarrow {{y}_{2}})\cdot ({{z}_{1}}\leftrightarrow {{z}_{2}}). \\ \end{array}\end{eqnarray} \tag{ 12 }$$

Here the contours of integration γ₁ and γ₂ are chosen in such a way that ${{\lambda }_{1}}={{r}_{1}}-{\rm i}\kappa /2+i\epsilon$, ${{\lambda }_{2}}={{r}_{2}}+{\rm i}\kappa /2+2i\epsilon$, ${{r}_{1,2}}\in \mathbb{R}$, and $\epsilon \gt 0$ is an infinitesimal parameter. The expression (12) is much more compact than (4), since it contains a twofold integration instead of a ${{2}^{M-1}}$-fold one. Nevertheless, it is still useful to use both approaches of [47] and [53] to unravel all properties of the scattering matrix $S_{M}^{(2)}$.

One can immediately note the fundamental properties of the two-photon scattering matrix (12). First, it does not contain any dependence on the positions of emitters, therefore the scattering is robust with respect to variations of the latter. Second, the expression (12) is invariant under a permutation of emitters—the products over the set of emitter do not change under it. Therefore, the scattering results do not depend on the ordering of emitters. They rather appear to be characteristic of sets of emitters than of individual emitters. This feature of (12) allows us to re-order emitters for the purposes of computational efficiency, in particular in the expression (4).

Modeling the incoming state by the Gaussian two-photon wavepacket

$$\begin{eqnarray}\begin{array}{rcl} {{\phi }_{2,{\rm in}}}({{x}_{1}},{{x}_{2}}) & = & \frac{1}{\sqrt{2\sigma \mu \pi }}{\rm exp} [{\rm i}(\Delta +\delta )({{x}_{1}}+{{x}_{2}})] \\ {} & {} & \times {\rm exp} [-{{({{x}_{1}}+{{x}_{2}})}^{2}}/8{{\mu }^{2}}-{{({{x}_{1}}-{{x}_{2}})}^{2}}/2{{\sigma }^{2}}], \\ \end{array}\end{eqnarray} \tag{ 13 }$$

where μ and σ are the variances of the center-of-mass coordinate $\frac{{{x}_{1}}+{{x}_{2}}}{2}$ and the relative coordinate $d={{x}_{1}}-{{x}_{2}}$ distributions, respectively, we are mainly interested in the regime $\mu \gg \sigma$. In the limit $\mu \to \infty$, the total energy $K={{k}_{1}}+{{k}_{2}}=k_{1}^{\prime }+k_{2}^{\prime }$ of the two photons is approximately conserved at the value $K=2\Delta +2\delta$, and the incoming wavepacket (14) acquires the factorized form

$$\begin{eqnarray}\begin{array}{rcl} {{\phi }_{2,{\rm in}}}(K,d) & = & \frac{1}{\sqrt{2\pi }}\int {\rm d}X{{{\rm e}}^{-{\rm i}KX}}{{\phi }_{2,{\rm in}}}\left( X+\frac{d}{2},X-\frac{d}{2} \right) \\ {} & \approx & \frac{\sqrt{\mu }}{\sqrt{2\sqrt{\pi }}}{{{\rm e}}^{-\frac{{{\mu }^{2}}}{2}{{(K-2\Delta -2\delta )}^{2}}}}{{\phi }_{2,{\rm in}}}(d), \\ \end{array}\end{eqnarray} \tag{ 14 }$$

where ${{\phi }_{2,{\rm in}}}(d)=\frac{1}{\sqrt{\sigma \sqrt{\pi }}}{{{\rm e}}^{-\frac{{{d}^{2}}}{2{{\sigma }^{2}}}}}$.

In this setting, the entire effect of scattering is visible in the relative part ${{\phi }_{2}}(d)$ of the two-photon wavefunction: if the energy is conserved, then ${{\phi }_{2,{\rm out}}}(K,d)$ can be also factorized like (14), and we can express the relative part ${{\phi }_{2,{\rm out}}}(d)$ of ${{\phi }_{2,{\rm out}}}(K,d)$ through ${{\phi }_{2,{\rm out}}}(d)$

$$\begin{eqnarray}&&{{\phi }_{2,{\rm out}}}(d)=\int {\rm d}(d^{\prime} )S_{M}^{(2)}(d,d^{\prime} ){{\phi }_{2,{\rm out}}}(d^{\prime} )\end{eqnarray} \tag{ 15 }$$

by means of the relative scattering matrix $S_{M}^{(2)}(d,d^{\prime} )$ depending on relative coordinates of photons d and $d^{\prime}$ in final and initial states, respectively. We also note that $S_{M}^{(2)}(d,d^{\prime} )$ also depend parametrically on $\delta$ which measures a detuning of the total energy from the two-photon resonance $2\Delta$.

We note that the representation of a scattering wavefunction in relative coodinates of photons is very advantageous, because of its direct relation to the correlation function ${{G}^{(2)}}(\tau )=\langle {\rm out}\mid {{a}^{\dagger }}(x){{a}^{\dagger }}(x+\tau )a(x+\tau )a(x)\mid {\rm out}\rangle =4{{\left| {{\phi }_{2,{\rm out}}}(d=\tau ) \right|}^{2}}$, which is an important measure of correlation effects between photons.

In the case of identical emitters ${{\Delta }_{a}}=\Delta$, ${{\kappa }_{a}}=\kappa$, the relative scattering matrix $S_{M}^{(2)}(d,d^{\prime} )$ is explicitly given by

$$\begin{eqnarray}&&S_{M}^{(2)}(d,d^{\prime} )=\delta (|d|-|{{d}^{\prime }}|)\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}\begin{array}{cll} {} & - & \frac{i}{(M-1)!}\frac{{{\partial }^{M-1}}}{\partial {{s}^{M-1}}}\left\{ \left[ {{{\rm e}}^{{\rm i}\left| |d|-|{{d}^{\prime }}| \right|(\delta +{\rm i}\kappa /2-s)}}+{{{\rm e}}^{{\rm i}(|d|+|{{d}^{\prime }}|)(\delta +{\rm i}\kappa /2-s)}} \right] \right. \\ {} & {} & \times {{\left. \frac{{{(s-{\rm i}\kappa )}^{M}}{{(s-2\delta )}^{M}}}{{{(s-2\delta -{\rm i}\kappa )}^{M}}} \right\}}_{s=0}} \\ \end{array}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}\begin{array}{cll} {} & - & \frac{\kappa }{(M-1)!}\frac{{{\partial }^{M-1}}}{\partial {{s}^{M-1}}}\left\{ {{{\rm e}}^{{\rm i}(|d|+|{{d}^{\prime }}|)(\delta +{\rm i}\kappa /2-s)}}\frac{{{(s-{\rm i}\kappa )}^{M}}}{(s-\delta )} \right. \\ {} & {} & \times {{\left. \left[ \frac{{{(s-2\delta )}^{M}}}{{{(s-2\delta -{\rm i}\kappa )}^{M}}}-\frac{{{s}^{M}}}{{{(s+{\rm i}\kappa )}^{M}}} \right] \right\}}_{s=0}}, \\ \end{array}\end{eqnarray} \tag{ 18 }$$

see the appendix for details of the derivation. This expression contains interesting effects which we discuss below.

Let us first remark the basic properties of $S_{M}^{(2)}(d,d^{\prime} )$. First, we note that it is symmetric

$$\begin{eqnarray}&&S_{M}^{(2)}(d,d^{\prime} )=S_{M}^{(2)}(d^{\prime} ,d).\end{eqnarray} \tag{ 19 }$$

Second, it obeys the reality condition

$$\begin{eqnarray}&&S_{M}^{(2)*}(d,d^{\prime} ;\delta )=S_{M}^{(2)}(d,d^{\prime} ;-\delta ).\end{eqnarray} \tag{ 20 }$$

Third, the convolution property is also fulfilled in the relative coordinates

$$\begin{eqnarray}\begin{array}{rcl} S_{M}^{(2)}(d,d^{\prime} ) & = & \int _{0}^{\infty }\mathop{\prod }\limits_{b=1}^{M-1}{\rm d}({{d}_{b}}){{S}^{(2)}}(d,{{d}_{M-1}}){{S}^{(2)}}({{d}_{M-1}},{{d}_{M-2}})\times \\ {} & \ldots & \times {{S}^{(2)}}({{d}_{2}},{{d}_{1}}){{S}^{(2)}}({{d}_{1}},d^{\prime} ). \\ \end{array}\end{eqnarray} \tag{ 21 }$$

Fourth, the unitarity condition implies that

$$\begin{eqnarray}&&\delta (|d|-|{{d}^{\prime }}|)=\int _{0}^{\infty }{\rm d}({{d}_{1}})S_{M}^{(2)*}(d,{{d}_{1}})S_{M}^{(2)}({{d}_{1}},d^{\prime} ).\end{eqnarray} \tag{ 22 }$$

The listed properties impose rigid constraints on possible scattering outcomes. In particular, in the resonance case $\delta =0$ (which in the present context means only that the total energy K matches with $2\Delta$, while the energies of both photons can differ from each other) the scattering result does not depend on the number M, but rather on its parity

$$\begin{eqnarray}&&S_{M\,\,{\rm odd}}^{(2)}(d,d^{\prime} )=S_{M=1}^{(2)}(d,d^{\prime} ),\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&S_{M\,\,{\rm even}}^{(2)}(d,d^{\prime} )=S_{M=2}^{(2)}(d,d^{\prime} )=\delta (|d|-|{{d}^{\prime }}|).\end{eqnarray} \tag{ 24 }$$

Moreover, the latter case of an array with an even number of scatterers is transparent for incident light

$$\begin{eqnarray}&&\phi _{2,{\rm out}}^{M\,\,{\rm even}}(d)={{\phi }_{2,{\rm in}}}(d).\end{eqnarray} \tag{ 25 }$$

To substantiate these conclusions, it is sufficient to notice that the unitarity condition (22) for M = 1 in combination with the symmetry (19) and reality (20) conditions yields $S_{M=2}^{(2)}=\delta (|d|-|{{d}^{\prime }}|)$, which holds by virtue of (21). The generalization for arbitrary M follows from the further application of the convolution property (21).

The same result (23), (24) follows from the explicit expression for $S_{M=1}^{(2)}(d,d^{\prime} )$ given by (16)–(18). This consideration elucidates the crucial role of the inelastic processes for the emergence of the parity effect (23), (24): the elastic contribution (17) vanishes at $\delta =0$, and only the inelastic contribution (18) makes the scattering off an odd number of emitters and off an even number of emitters distinguishable. In addition, we note that the parity effects in a setup similar to ours and driven by a classical field have been discussed in [55].

The transparency effect expressed in (24) can be straightforwardly generalized to arrays of emitters with a set of different detunings $\{{{\delta }_{a}}\}$ obeying a constraint such that for each ${{\delta }_{a}}$ there exists ${{\delta }_{{\bar{a}}}}=-{{\delta }_{a}}$ (for even M). In fact, reshuffling the scattering matrices of individual scatterers in (21) (which is allowed by the permutation symmetry discussed above), we can arrange them pairwise. The convolution within each pair $({{\delta }_{a}},{{\delta }_{{\bar{a}}}})$ yields the identity ($\delta$-function) by virtue of the symmetry, the reality, and the unitarity conditions; a convolution of the identities is again the identity.

Since the M = 1 case is important for the understanding of the scattering off an odd number of emitters, we briefly revisit it. For $\delta =0$ we have

$$\begin{eqnarray}&&{{S}^{(2)}}(d,d^{\prime} )=\delta (|d|-|{{d}^{\prime }}|)-2\kappa {{{\rm e}}^{-\kappa (|d|+|{{d}^{\prime }}|)/2}}.\end{eqnarray} \tag{ 26 }$$

For the Gaussian initial condition we obtain

$$\begin{eqnarray}\begin{array}{rcl} {{\phi }_{2,{\rm out}}}(d) & = & \frac{1}{\sqrt{\sigma \sqrt{\pi }}}\left[ {{{\rm e}}^{-\frac{{{d}^{2}}}{2{{\sigma }^{2}}}}}-\sqrt{2\pi }\kappa \sigma {{{\rm e}}^{-\frac{\kappa |d|}{2}+\frac{{{\sigma }^{2}}{{\kappa }^{2}}}{8}}}{\rm erfc}\left( \frac{\sigma \kappa }{2\sqrt{2}} \right) \right] \\ {} & {} & \mathop{\approx }\limits^{\sigma \gg {{\kappa }^{-1}},|d|}\frac{1}{\sqrt{\sigma \sqrt{\pi }}}\left[ 1-4{{{\rm e}}^{-\frac{\kappa |d|}{2}}} \right]. \\ \end{array}\end{eqnarray} \tag{ 27 }$$

This wavefunction describes a bound state of two photons. If the initial variance σ of relative distances between photons is sufficiently large, $\sigma \gg {{\kappa }^{-1}}$, then the final distance is distributed on the smaller length scale κ⁻¹, which indicates an emerging effective attraction between the photons.

Another important consequence for the scattering off an odd number of emitters which can be derived from (27) is that for a special choice of the parameter $\sigma \kappa$ one can observe the antibunching behavior of photons. It is characterized by the vanishing value of ${{G}^{(2)}}(0)\propto {{\left| {{\phi }_{2,{\rm out}}}(0) \right|}^{2}}$, which indeed happens at $\sigma \kappa \approx 0.5$. In contrast, for an even number of emitters ${{G}^{(2)}}(0)$ does not vanish for any $\sigma \kappa$, which indicates a tendency toward bunching. Therefore we use the value $\sigma \kappa =0.5$ in figure 3(b) depicting $|{{\phi }_{2,{\rm out}}}(d){{|}^{2}}$ in order to emphasize the qualitative difference between scattering off even and odd numbers of emitters.

Zoom In Zoom Out Reset image size

A much richer structure of an outgoing wavefunction appears at a finite detuning $\delta$. Assuming now ${{\delta }^{-1}},{{\kappa }^{-1}}\gg \sigma$, we observe that it develops a polynomial dependence on d

$$\begin{eqnarray}&&{{\phi }_{2,{\rm out}}}(d)={{\phi }_{2,{\rm in}}}(d)-\sqrt{2\sigma \sqrt{\pi }}\kappa {{{\rm e}}^{-|d|(\kappa /2-{\rm i}\delta )}}{{P}_{M-1}}(\kappa |d|;\delta /\kappa ),\end{eqnarray} \tag{ 28 }$$

where the polynomial of the order $M-1$

$$\begin{eqnarray}&&{{P}_{M-1}}(x;\hat{\delta })=\frac{{{{\rm e}}^{x}}}{(M-1)!}\frac{{{\partial }^{M-1}}}{\partial {{s}^{M-1}}}{{\left\{ {{{\rm e}}^{-sx}}{{s}^{M}}\frac{(s-{\rm i}\hat{\delta }-\frac{1}{2}){{(s-2{\rm i}\hat{\delta }-1)}^{M}}}{(s-{\rm i}\hat{\delta }-1){{(s-2{\rm i}\hat{\delta })}^{M}}} \right\}}_{s=1}}\quad \end{eqnarray} \tag{ 29 }$$

has complex-valued coefficients.

At a large detuning $\delta \gg \kappa$ the elastic scattering dominates, ${{P}_{M-1}}(x;\hat{\delta })$ turns into the Laguerre polynomial $L_{M-1}^{(1)}(x)$ with the real-valued coefficients, and ${{\phi }_{2,{\rm out}}}(d)$ exhibits the same behavior as the one-photon scattering function (8): the probability $|{{\phi }_{2,{\rm out}}}(d){{|}^{2}}$ oscillates on the scale κ⁻¹ featuring precisely $M-1$ nodes. We note that the real-valuedness of the coefficients of ${{P}_{M-1}}(x;\hat{\delta })$ is important for the presence of the nodes in $|{{\phi }_{2,{\rm out}}}(d){{|}^{2}}$, which are destroyed at a small detuning $\delta \ll \kappa$ by the emerging effective interaction between photons.

Thus, tuning $\delta$ from $\delta \gg \kappa$ to $\delta \ll \kappa$ we observe a qualitative change in the scattering properties (see figure 3(a)), characterized by the increasing role of the inelastic processes. The formation of bound states, transitions from bunching to antibunching, and the clearly pronounced parity effect give full evidence of strong correlations between photons resulting from their scattering off arrays of emitters in a 1D chiral channel.

2.3. Robustness of correlated states

It has been already discussed above that scattering results are independent of the positions of emitters, which ensures their robustness against fluctuations of the latter. Next, we study how robust scattering results are with respect to randomness in the emitters' transition frequencies. To this end, we consider detunings ${{\delta }_{a}}$ as random variables which are normally distributed around the mean value $\delta$ with the variance Σ, and evaluate the probability distribution

$$\begin{eqnarray}&&\begin{array}{l} P[|{{\phi }_{2,{\rm out}}}(d){{|}^{2}}] \\ =\int \mathop{\prod }\limits_{a^{\prime} =1}^{M}\left( \frac{{{{\rm e}}^{-\frac{{{({{\delta }_{a^{\prime} }}-\delta )}^{2}}}{2{{\Sigma }^{2}}}}}{\rm d}{{\delta }_{a^{\prime} }}}{\Sigma \sqrt{2\pi }} \right)\delta \left( |{{\phi }_{2,{\rm out}}}(d){{|}^{2}}-|{{\phi }_{2,{\rm out}}}(d;\{{{\delta }_{a}}\}){{|}^{2}} \right) \\ \end{array}\end{eqnarray} \tag{ 30 }$$

with the help of the exact expression for $|{{\phi }_{2,{\rm out}}}(d;\{{{\delta }_{a}}\}){{|}^{2}}$ based on equation (12). In figure 4 we plot typical results of this averaging showing mean values as well as median and mean absolute deviations of the distribution (30). We observe that for not so large $\Sigma \lesssim \kappa$, this distribution is sufficiently narrow so that it does not mask the qualitative effects (parity, antibunching, etc) discussed above. Thus, we conclude that our results are robust with respect to the fluctuations in transition frequencies as well.

Zoom In Zoom Out Reset image size

2.4. Effect of losses

We have suggested a novel approach to producing strongly correlated states of two and more photons in a highly controllable way. Our setup uses the edge states of photonic topological insulators with multi-level emitters coupled to it. The scattering in this setup possesses the following universal features.

First, multi-particle scattering does not depend on the positions of emitters. This is a general consequence of the absence of backscattering, and it makes the scattering outcomes robust against the fluctuations in the emitters' coordinates. This property also provides us with a powerful theoretical tool for calculating the scattering results based on the convolution property (4).

Second, the scattered wavepacket has a polynomial structure and the minima of the outgoing pulse are given by the zeroes of this polynomial. In particular, single-particle scattering (and the reducible part of multi-particle scattering in general) is described by Laguerre polynomials. In single-photon scattering, the outcome looks like a wavepacket fragmented between the nodes—the zeroes of this polynomial. The emergence of the nodes is accounted by time delays on each emitter. This leads to the non-monotonic behavior of the G² correlation function thus signifying correlation between photons. In multi-photon scattering, this picture takes place at large detunings $\delta$, i.e., in the regime where the elastic processes dominate over the inelastic ones.

Third, the role of the inelastic processes enhances as we tune to small values of $\delta$. These processes lead to the emergence of the effective interaction between photons, which lies at the origin of the exciting effects: the parity-dependent scattering and the antibunching of photons. In particular, we have observed for the first time that at a vanishing detuning the two-photon scattering matrix in the chiral edge channel only depends on the parity of the number of emitters, and does not depend on the number itself. We substantiated our observation by the fundamental symmetry arguments. For an even number of emitters, we proposed certain array configurations which are transparent to arbitrary two-photon wavepackets. For an odd number of emitters, we discussed how to adjust the parameters of the initial two-phonon wavepackets in order to observe the antibunching behavior—the effective repulsion of photons—in the 1D geometry.

The physical picture behind the parity effect is the following. After the first incident photon is absorbed by the first emitter, the latter becomes transparent for the second photon. Thus, the first photon is delayed by the time $\sim {{\kappa }^{-1}}$ with respect to the second photon. This leads to antibunching of photons via scattering on a single emitter. If the second emitter is present, it absorbs the second photon and captures it during the time $\sim {{\kappa }^{-1}}$. This gives the first photon a possibility to compensate the lag, and both photons end up with the same relative distance as it was in the incoming state (note that the compensation is exact, only if the resonance condition is fulfilled, and if the emitters are identical). Repeating these arguments for odd and even numbers of identical emitters, we observe that the net result is the same as for one and two emitters, respectively. Deviations from the perfect conditions would lead to a difference between scattering results on, say, M = 1 and M = 3 emitters, but they still remain closer and qualitatively similar to each other, than to a scattering result on M = 2 emitters.

These exciting effects result from the combination of the reduced dimensionality, which increases the role of the photonic correlations, and the topological origin of the chiral edge states propagating without the backscattering. For their experimental observation, we specify below the necessary physical parameters. Photonic correlations are a vital ingredient for building up working schemes of future photonic quantum technology [5]. We believe that the effects described here will be used in future functional quantum photonic devices or for quantum simulation of condensed matter systems using photonic setups.

Physical realization. To observe the quantum many-body effects of photons described above one needs a photonic topological insulator analogous to the one studied in [34, 35]. The role of emitters can be played by either quantum dots or superconducting qubits, like, e.g., in [11] and [17]. Available single-photon emitters made of quantum dots and coupled to photonic crystal waveguides [9, 12, 20] can also serve for a realization of our model. Typical bandgap frequencies of existing photonic topological insulators lying in the GHz range match with the transition frequencies of single-photon emitters, which are necessary to inject photons into the 1D channel. The parameters provided in [34] allow us to estimate the topological photonic band gap $B\sim$ 0.27 GHz as well as the group velocity $c/{{v}_{{\rm g}}}\approx 175$. The electromagnetic field is bounded within an area of $\sim$ 1600 mm², which gives ${{A}_{0}}/{{A}_{{\rm eff}}}\sim 1$. The dielectric host material has a relatively large dielectric constant (see the Method summary section of [34]). On the basis of (1) we conclude that the Purcell factor of the setup of [34] is of the order of a few tens. The other setup described in [39] is made of an array of helix-like waveguides with the radius R. The spectral gap and the group velocity of the edge state (see figure 2(c) of [39]) depend nonmonotonically on R; the group velocity can even go down close to zero at $R\approx 17\;\mu$m, rendering $c/{{v}_{{\rm g}}}$ very large. This leads us to much more optimistic values of $c/{{v}_{{\rm g}}}\sim {{10}^{3}}$ and a Purcell factor of the same order. A similar analysis of the third setup presented in [38] yields a Purcell factor $\sim$ 20. A further enhancement of the Purcell factor in the near future can facilitate the experimental realization of our proposal.

In this paper we discussed the scattering results in arrays of two-level emitters. The present consideration, however, can be generalized to arrays of three- and four-level emitters, which are interesting due to the experimentally observed effects [8, 56, 57] of the single-emitter electromagnetically induced transparency and the transistor-like behavior, and the recent proposals to implement four-level emitters in the engineering of quantum gates, e.g. CNOT [58]. By combining emitters of different level structures one gains flexibility and the ability to build optical schemes with desired correlations of the outgoing photonic wavefunctions. For an explicit evaluation of the outgoing state one can use the scattering formalism of [47] in combination with the convolution property (4).

Currently available topological photonic crystals operate in the GHz frequency range. This is related to the difficulties in creating a large magneto-optical response with existing materials, which is needed for a nontrivial topology of the band structure. However, some recent developments [59] in material science might lead to a significant enhancement of the magneto-optical response in the optical domain, which can be used to engineer topological photonic crystals in the optical range of frequencies.

Finally, we note that interacting photons in two dimensions can lead to anyonic statistics and fractional Hall states [60]. The recent suggestion of realizing Dirac cone structure [62] may pave the way to exotic physics of interacting two-dimensional photonic fluids. Moreover, the recent classification of interacting bosonic topological insulators [61] suggests a possible route to engineer other exotic states of light.

We expect that the novel correlated states of photons and other propagating bosons will find applications in quantum information technology, optomechanics, and precision measurements.

We would like to thank D Baeriswyl, D Fioretto, and V Yudson for useful discussions. MR and VG are supported by MaNEP and Swiss NSF. MP is supported by DFG-FG 723 and 912. VG thanks KITP for hospitality.

In the case of two-level emitters and identical couplings the model (2) is exactly solvable by the Bethe ansatz [53]. For a more general problem involving arrays of emitters with an arbitrary level structure, a complementary diagrammatic approach to scattering is available [47].

Below we provide the derivation of formulas appearing in the main text based on these approaches.

Appendix. Single-photon scattering

First, we rewrite the single-photon scattering matrix (6) in the coordinate representation

$$\begin{eqnarray}\begin{array}{rcl} S_{M}^{(1)}(y,z) & = & \frac{1}{2\pi }\int {\rm d}k{\rm d}k^{\prime} {{{\rm e}}^{{\rm i}ky-{\rm i}k^{\prime} z}}S_{M;kk^{\prime} }^{(1)} \\ {} & = & \delta (y-z)-\Theta (z-y)\mathop{\sum }\limits_{a=1}^{M}{{\kappa }_{a}}\mathop{\prod }\limits_{a^{\prime} =1;a\ne a}^{M}\frac{{{\Delta }_{a}}-{\rm i}{{\kappa }_{a}}/2-{{\Delta }_{a^{\prime} }}-{\rm i}{{\kappa }_{a^{\prime} }}/2}{{{\Delta }_{a}}-{\rm i}{{\kappa }_{a/2}}-{{\Delta }_{a^{\prime} }}+{\rm i}{{\kappa }_{a^{\prime} }}/2} \\ {} & {} & \times {{{\rm e}}^{-({\rm i}{{\Delta }_{a}}+{{\kappa }_{a}}/2)(z-y)}}. \\ \end{array}\end{eqnarray} \tag{ A.1 }$$

The corresponding scattered part ${{\phi }_{1,{\rm scatt}}}$ is given by

$$\begin{eqnarray}\begin{array}{rcl} {{\phi }_{1,{\rm scatt}}}(x-t) & = & -\int {\rm d}x^{\prime} \Theta (x^{\prime} )\mathop{\sum }\limits_{a=1}^{M}{{\kappa }_{a}}\mathop{\prod }\limits_{a^{\prime} =1;a^{\prime} \ne a}^{M}\frac{{{\Delta }_{a}}-{\rm i}{{\kappa }_{a}}/2-{{\Delta }_{a^{\prime} }}-{\rm i}{{\kappa }_{a^{\prime} }}/2}{{{\Delta }_{a}}-{\rm i}{{\kappa }_{a}}/2-\Delta _{a^{\prime} }^{{}}+{\rm i}{{\kappa }_{a^{\prime} }}/2} \\ {} & {} & \times {{{\rm e}}^{-({\rm i}{{\Delta }_{a}}+{{\kappa }_{a}}/2)x^{\prime} }}{{\phi }_{1,{\rm in}}}(x^{\prime} +x-t). \\ \end{array}\end{eqnarray} \tag{ A.2 }$$

In the limit ${{\Delta }_{a}}\to \Delta$, $\kappa \to {{\kappa }_{a}}$ we obtain

$$\begin{eqnarray}\begin{array}{rcl} S_{M}^{(1)}(y,z) & = & \frac{1}{2\pi }\int {\rm d}k{{{\rm e}}^{{\rm i}k(y-z)}}{{\left( 1-\frac{{\rm i}\kappa }{k-\Delta +{\rm i}\kappa /2} \right)}^{M}} \\ {} & = & \delta (y-z)+\mathop{\sum }\limits_{m=1}^{M}C_{m}^{M}\frac{{{(-{\rm i}\kappa )}^{m}}}{2\pi }\int {\rm d}k\frac{{{{\rm e}}^{{\rm i}k(y-z)}}}{{{(k-\Delta +{\rm i}\kappa /2)}^{m}}} \\ {} & = & \delta (y-z)-i\Theta (z-y)\mathop{\sum }\limits_{m=1}^{M}C_{m}^{M}\frac{{{(-{\rm i}\kappa )}^{m}}{{[i(y-z)]}^{m-1}}}{(m-1)!}{{{\rm e}}^{({\rm i}\Delta +\kappa /2)(y-z)}} \\ \end{array}\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&=\delta (y-z)-\kappa \Theta (z-y)L_{M-1}^{(1)}(\kappa (z-y)){{{\rm e}}^{({\rm i}\Delta +\kappa /2)(y-z)}},\end{eqnarray} \tag{ A.4 }$$

where $L_{M-1}^{(1)}(x)$ are the associated Laguerre polynomials and C_m^M are the binomial coefficients.

Appendix. Two-photon scattering

The two-photon scattering matrix (12) can be parameterized by the detuning $\delta$ (recall that the total energy ${{k}_{1}}+{{k}_{2}}=k_{1}^{\prime }+k_{2}^{\prime }\equiv K=2\Delta +2\delta$ is conserved) and the two relative coordinates $d={{y}_{1}}-{{y}_{2}}$ and $d^{\prime} ={{z}_{1}}-{{z}_{2}}$.

First, we separate elastic and inelastic contributions to (12). Picking the unity in the square brackets in (12), we see that the two integrations disentangle, and we recognize the elastic contribution $S_{M}^{(1)}({{y}_{1}},{{z}_{1}})S_{M}^{(1)}({{y}_{2}},{{z}_{2}})+S_{M}^{(1)}({{y}_{1}},{{z}_{2}})S_{M}^{(1)}({{y}_{2}},{{z}_{1}})$ to the two-photon scattering matrix (10). The remaining term in the square brackets in (12) generates the inelastic contribution

$$\begin{eqnarray}\begin{array}{rcl} {\rm i}\mathcal{T}_{M}^{(2)}({{y}_{1}},{{y}_{2}},{{z}_{1}},{{z}_{2}}) & = & -2{\rm i}{{\kappa }^{3}}\Theta ({{z}_{1}}\gt {{z}_{2}}\gt {{y}_{2}}\gt {{y}_{1}}) \\ {} & {} & \times \mathop{\sum }\limits_{a,b}\frac{{{C}_{a}}{{C}_{b}}}{{{\Delta }_{a}}-{{\Delta }_{b}}+{\rm i}\kappa }{{{\rm e}}^{{\rm i}({{\Delta }_{a}}-{\rm i}\kappa /2)({{y}_{1}}-{{z}_{1}})+{\rm i}({{\Delta }_{b}}-{\rm i}\kappa /2)({{y}_{2}}-{{z}_{2}})}} \\ {} & {} & +({{y}_{1}}\leftrightarrow {{y}_{2}})\cdot ({{z}_{1}}\leftrightarrow {{z}_{2}}), \\ \end{array}\end{eqnarray} \tag{ A.5 }$$

where

$$\begin{eqnarray}&&{{C}_{a}}=\mathop{\prod }\limits_{a^{\prime} =1;a^{\prime} \ne a}^{M}\frac{{{\Delta }_{a}}-{{\Delta }_{a^{\prime} }}-{\rm i}\kappa }{{{\Delta }_{a}}-{{\Delta }_{a^{\prime} }}}.\end{eqnarray} \tag{ A.6 }$$

The elastic (reducible) part can be cast to

$$\begin{eqnarray}\begin{array}{rcl} S_{M;{\rm red}}^{(2)}(d,{{d}^{\prime }}) & = & \frac{1}{2\pi }\int {\rm d}k\left[ {{{\rm e}}^{{\rm i}k(d-{{d}^{\prime }})}}+{{{\rm e}}^{-{\rm i}k(d+{{d}^{\prime }})}} \right] \\ {} & {} & \times {{\left( \frac{\delta +k-{\rm i}\kappa /2}{\delta +k+{\rm i}\kappa /2} \right)}^{M}}{{\left( \frac{\delta -k-{\rm i}\kappa /2}{\delta -k+{\rm i}\kappa /2} \right)}^{M}} \\ {} & = & \frac{1}{2\pi }\int {\rm d}k\left[ {{{\rm e}}^{{\rm i}k(|d|-|{{d}^{\prime }}|)}}+{{{\rm e}}^{-{\rm i}k(|d|+|{{d}^{\prime }}|)}} \right] \\ {} & {} & \times {{\left( \frac{\delta +k-{\rm i}\kappa /2}{\delta +k+{\rm i}\kappa /2} \right)}^{M}}{{\left( \frac{\delta -k-{\rm i}\kappa /2}{\delta -k+{\rm i}\kappa /2} \right)}^{M}}. \\ \end{array}\end{eqnarray} \tag{ A.7 }$$

Evaluating this integral in terms of the Mth order residues, we obtain the contributions (16) and (17).

The inelastic (irreducible) part (A.5) admits the integral representation

$$\begin{eqnarray}\begin{array}{rcl} i\mathcal{T}_{M}^{(2)}({{y}_{1}},{{y}_{2}},{{z}_{1}},{{z}_{2}}) & = & -2{{\kappa }^{3}}\Theta ({{z}_{1}}\gt {{z}_{2}}\gt {{y}_{2}}\gt {{y}_{1}}) \\ {} & {} & \times \int _{0}^{\infty }{\rm d}\tau {{{\rm e}}^{-\kappa \tau }}\left( \mathop{\sum }\limits_{a}{{C}_{a}}{{{\rm e}}^{{\rm i}({{\Delta }_{a}}-{\rm i}\kappa /2)(\tau +{{y}_{1}}-{{z}_{1}})}} \right) \\ {} & {} & \times \left( \mathop{\sum }\limits_{b}{{C}_{b}}{{{\rm e}}^{{\rm i}({{\Delta }_{b}}-{\rm i}\kappa /2)(-\tau +{{y}_{2}}-{{z}_{2}})}} \right) \\ {} & {} & +({{y}_{1}}\leftrightarrow {{y}_{2}})\cdot ({{z}_{1}}\leftrightarrow {{z}_{2}}), \\ \end{array}\end{eqnarray} \tag{ A.8 }$$

which helps to decouple the sums over indices a and b. Assuming in the following ${{z}_{1}}\gt {{z}_{2}}\gt {{y}_{2}}\gt {{y}_{1}}$ and employing the identities (see equation (A.4))

$$\begin{eqnarray}&&\begin{array}{l} \mathop{{\rm lim} }\limits_{{{\Delta }_{a}}\to \Delta }\Theta (x)\mathop{\sum }\limits_{a}{{C}_{a}}{{{\rm e}}^{-{\rm i}({{\Delta }_{a}}-{\rm i}\kappa /2)x}} \\ =-\frac{\Theta (x)}{2\pi }\int {\rm d}s{{{\rm e}}^{-{\rm i}sx}}\mathop{\sum }\limits_{m=1}^{M}C_{m}^{M}\frac{{{(-{\rm i}\kappa )}^{m}}}{{{(s-\Delta +{\rm i}\kappa /2)}^{m}}}, \\ \end{array}\end{eqnarray} \tag{ A.9 }$$

$$\begin{eqnarray}&&\begin{array}{l} \mathop{{\rm lim} }\limits_{{{\Delta }_{a}}\to \Delta }\Theta (-x)\mathop{\sum }\limits_{a}{{C}_{a}}{{{\rm e}}^{-{\rm i}({{\Delta }_{a}}+{\rm i}\kappa /2)x}} \\ =-\frac{\Theta (-x)}{2\pi }\int {\rm d}s{{{\rm e}}^{-{\rm i}sx}}\mathop{\sum }\limits_{m=1}^{M}C_{m}^{M}\frac{{{({\rm i}\kappa )}^{m}}}{{{(s-\Delta -{\rm i}\kappa /2)}^{m}}}, \\ \end{array}\end{eqnarray} \tag{ A.10 }$$

we transform (A.8) into

$$\begin{eqnarray}\begin{array}{rcl} i\mathcal{T}_{M}^{(2)}({{y}_{1}},{{y}_{2}},{{z}_{1}},{{z}_{2}}) & = & -2\kappa \int \frac{{\rm d}s}{2\pi }\int \frac{{\rm d}r}{2\pi }\int _{0}^{{{z}_{1}}-{{y}_{1}}-{{0}^{+}}}{\rm d}\tau {{{\rm e}}^{-{\rm i}(r-s-{\rm i}\kappa )\tau -{\rm i}s({{z}_{1}}-{{y}_{1}})-{\rm i}r({{z}_{2}}-{{y}_{2}})}} \\ {} & {} & \times \mathop{\sum }\limits_{M}^{{{m}_{1}}=1}\mathop{\sum }\limits_{M}^{{{m}_{2}}=1}C_{{{m}_{1}}}^{M}C_{{{m}_{2}}}^{M}\frac{{{(-{\rm i}\kappa )}^{{{m}_{1}}}}}{{{(s-\Delta +{\rm i}\kappa /2)}^{{{m}_{1}}}}}\frac{{{(-{\rm i}\kappa )}^{{{m}_{2}}}}}{{{(r-\Delta +{\rm i}\kappa /2)}^{{{m}_{2}}}}} \\ {} & {} & +2\kappa \int \frac{{\rm d}s}{2\pi }\int \frac{{\rm d}r}{2\pi }\int _{{{z}_{1}}-{{y}_{1}}+{{0}^{+}}}^{\infty }{\rm d}\tau {{{\rm e}}^{-{\rm i}(s+r-{\rm i}{{0}^{+}})\tau -{\rm i}(s+{\rm i}\kappa )({{y}_{1}}-{{z}_{1}})-{\rm i}r({{z}_{2}}-{{y}_{2}})}} \\ {} & {} & \times \mathop{\sum }\limits_{M}^{{{m}_{1}}=1}\mathop{\sum }\limits_{M}^{{{m}_{2}}=1}C_{{{m}_{1}}}^{M}C_{{{m}_{2}}}^{M}\frac{{{({\rm i}\kappa )}^{{{m}_{1}}}}}{{{(s+\Delta +{\rm i}\kappa /2)}^{{{m}_{1}}}}}\frac{{{(-{\rm i}\kappa )}^{{{m}_{2}}}}}{{{(r-\Delta +{\rm i}\kappa /2)}^{{{m}_{2}}}}}. \\ \end{array}\end{eqnarray} \tag{ A.11 }$$

Integrating over τ we obtain

$$\begin{eqnarray*}\begin{array}{rcl} {\rm i}\mathcal{T}_{M}^{(2)}({{y}_{1}},{{y}_{2}},{{z}_{1}},{{z}_{2}}) & = & 2{\rm i}\kappa \int \frac{{\rm d}s}{2\pi }\int \frac{{\rm d}r}{2\pi }\frac{{{{\rm e}}^{-{\rm i}s({{z}_{1}}-{{y}_{1}})}}{{{\rm e}}^{-{\rm i}r({{z}_{2}}-{{y}_{2}})}}}{r-s-{\rm i}\kappa } \\ {} & {} & \times \mathop{\sum }\limits_{M}^{{{m}_{1}}=1}\mathop{\sum }\limits_{M}^{{{m}_{2}}=1}C_{{{m}_{1}}}^{M}C_{{{m}_{2}}}^{M}\frac{{{(-{\rm i}\kappa )}^{{{m}_{1}}}}}{{{(s-\Delta +{\rm i}\kappa /2)}^{{{m}_{1}}}}}\frac{{{(-{\rm i}\kappa )}^{{{m}_{2}}}}}{{{(r-\Delta +{\rm i}\kappa /2)}^{{{m}_{2}}}}} \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray}\begin{array}{l} {} & {} & -2{\rm i}\kappa \int \frac{{\rm d}s}{2\pi }\int \frac{{\rm d}r}{2\pi }\frac{{{{\rm e}}^{-\kappa ({{z}_{1}}-{{y}_{1}})}}{{{\rm e}}^{-{\rm i}r({{z}_{1}}-{{y}_{1}}+{{z}_{2}}-{{y}_{2}})}}}{s+r-{\rm i}{{0}^{+}}} \\ {} & {} & \times \mathop{\sum }\limits_{M}^{{{m}_{1}}=1}\mathop{\sum }\limits_{M}^{{{m}_{2}}=1}C_{{{m}_{1}}}^{M}C_{{{m}_{2}}}^{M}\frac{{{({\rm i}\kappa )}^{{{m}_{1}}}}}{{{(s+\Delta +{\rm i}\kappa /2)}^{{{m}_{1}}}}}\frac{{{(-{\rm i}\kappa )}^{{{m}_{2}}}}}{{{(r-\Delta +{\rm i}\kappa /2)}^{{{m}_{2}}}}}. \\ \end{array}\end{eqnarray} \tag{ A.12 }$$

The second term exactly compensates the contribution from the pole $s=r-{\rm i}\kappa$ in the first term. So we can rewrite (A.13) as

$$\begin{eqnarray}\begin{array}{rcl} {\rm i}\mathcal{T}_{M}^{(2)}({{y}_{1}},{{y}_{2}},{{z}_{1}},{{z}_{2}}) & = & -2{\rm i}\kappa {{\int }_{\gamma _{1}^{\prime }}}\frac{{\rm d}s}{2\pi }{{\int }_{\gamma _{2}^{\prime }}}\frac{{\rm d}r}{2\pi }\frac{{{{\rm e}}^{-{\rm i}s({{z}_{1}}-{{y}_{1}})}}{{{\rm e}}^{-{\rm i}r({{z}_{2}}-{{y}_{2}})}}}{s-r+{\rm i}\kappa } \\ {} & {} & \times \mathop{\sum }\limits_{M}^{{{m}_{1}}=1}\mathop{\sum }\limits_{M}^{{{m}_{2}}=1}C_{{{m}_{1}}}^{M}C_{{{m}_{2}}}^{M}\frac{{{(-{\rm i}\kappa )}^{{{m}_{1}}}}}{{{(s-\Delta +{\rm i}\kappa /2)}^{{{m}_{1}}}}}\frac{{{(-{\rm i}\kappa )}^{{{m}_{2}}}}}{{{(r-\Delta +{\rm i}\kappa /2)}^{{{m}_{2}}}}}, \\ \end{array}\end{eqnarray} \tag{ A.13 }$$

where the contours of integration are deformed to the small circles $\gamma _{1}^{\prime }$ and $\gamma _{2}^{\prime }$ embracing clockwise the poles $s=\Delta -{\rm i}\kappa /2$ and $r=\Delta -{\rm i}\kappa /2$, respectively.

Transforming (A.13) plus its complements $({{y}_{1}}\leftrightarrow {{y}_{2}})\cdot ({{z}_{1}}\leftrightarrow {{z}_{2}})$ to the mixed representation of the total energy $K=2\Delta +2\delta$ and the relative coordinates $d,d^{\prime}$, we obtain

$$\begin{eqnarray}\begin{array}{rcl} i\mathcal{T}_{M}^{(2)}(d,d^{\prime} ) & = & 2\kappa {{\int }_{\gamma _{1}^{\prime }}}\frac{{\rm d}s}{2\pi }{{\int }_{\gamma _{2}^{\prime }}}\frac{{\rm d}r}{2\pi }\frac{{{{\rm e}}^{{\rm i}(\Delta +\delta -s)(|d|+|{{d}^{\prime }}|)}}}{(s-r+{\rm i}\kappa )(2\Delta +2\delta -s-r)} \\ {} & {} & \times \mathop{\sum }\limits_{M}^{{{m}_{1}}=1}\mathop{\sum }\limits_{M}^{{{m}_{2}}=1}C_{{{m}_{1}}}^{M}C_{{{m}_{2}}}^{M}\frac{{{(-{\rm i}\kappa )}^{{{m}_{1}}}}}{{{(s-\Delta +{\rm i}\kappa /2)}^{{{m}_{1}}}}}\frac{{{(-{\rm i}\kappa )}^{{{m}_{2}}}}}{{{(r-\Delta +{\rm i}\kappa /2)}^{{{m}_{2}}}}} \\ \end{array}\end{eqnarray} \tag{ A.14 }$$

$$\begin{eqnarray}\begin{array}{cl} {} & = & 2\kappa {{\int }_{\gamma _{1}^{\prime\prime }}}\frac{{\rm d}s}{2\pi }{{\int }_{\gamma _{2}^{\prime\prime }}}\frac{{\rm d}r}{2\pi }\frac{{{{\rm e}}^{-{\rm i}(s-\delta -{\rm i}\kappa /2)D}}}{(s-r+{\rm i}\kappa )(2\delta -s-r+{\rm i}\kappa )} \\ {} & {} & \times \mathop{\sum }\limits_{M}^{{{m}_{1}}=1}\mathop{\sum }\limits_{M}^{{{m}_{2}}=1}C_{{{m}_{1}}}^{M}C_{{{m}_{2}}}^{M}\frac{{{(-{\rm i}\kappa )}^{{{m}_{1}}}}}{{{s}^{{{m}_{1}}}}}\frac{{{(-{\rm i}\kappa )}^{{{m}_{2}}}}}{{{r}^{{{m}_{2}}}}}, \\ \end{array}\end{eqnarray} \tag{ A.15 }$$

where $D=|d|+|{{d}^{\prime }}|$, and the contours $\gamma _{1}^{\prime\prime }$ and $\gamma _{2}^{\prime\prime }$ embraces the origin clockwise. Noting that

$$\begin{eqnarray}\begin{array}{l} {} & {} & \mathop{\sum }\limits_{M}^{{{m}_{1}}=1}\mathop{\sum }\limits_{M}^{{{m}_{2}}=1}C_{{{m}_{1}}}^{M}C_{{{m}_{2}}}^{M}\frac{{{(-{\rm i}\kappa )}^{{{m}_{1}}}}}{{{s}^{{{m}_{1}}}}}\frac{{{(-{\rm i}\kappa )}^{{{m}_{2}}}}}{{{r}^{{{m}_{2}}}}} \\ {} & = & \left[ \frac{{{(s-{\rm i}\kappa )}^{M}}}{{{s}^{M}}}-1 \right]\left[ \frac{{{(r-{\rm i}\kappa )}^{M}}}{{{r}^{M}}}-1 \right], \\ \end{array}\end{eqnarray} \tag{ A.16 }$$

we cast (A.15) to

$$\begin{eqnarray}&&\begin{array}{r} {\rm i}\mathcal{T}_{M}^{(2)}(d,d^{\prime} ) \\ =-\frac{2\kappa {{{\rm e}}^{({\rm i}\delta -\kappa /2)D}}}{{{[(M-1)!]}^{2}}}\frac{{{\partial }^{2\,M-2}}}{\partial {{s}^{M-1}}\partial {{r}^{M-1}}}{{\left[ \frac{{{{\rm e}}^{-{\rm i}sD}}{{(s-{\rm i}\kappa )}^{M}}{{(r-{\rm i}\kappa )}^{M}}}{(s-r+{\rm i}\kappa )(2\delta -s-r+{\rm i}\kappa )} \right]}_{s=r=0}}. \\ \end{array}\end{eqnarray} \tag{ A.17 }$$

Representing

$$\begin{eqnarray}\begin{array}{rcl} \frac{{{(r-{\rm i}\kappa )}^{M}}}{s-r+{\rm i}\kappa } & = & -{{(r-{\rm i}\kappa )}^{M-1}}\frac{1}{1-\frac{s}{r-{\rm i}\kappa }}=-\mathop{\sum }\limits_{l=0}^{\infty }{{s}^{l}}{{(r-{\rm i}\kappa )}^{M-1-l}} \\ {} & = & -\mathop{\sum }\limits_{l=0}^{M-1}{{s}^{l}}{{(r-{\rm i}\kappa )}^{M-1-l}}-{{s}^{M}}\mathop{\sum }\limits_{l=0}^{\infty }{{s}^{l}}{{(r-{\rm i}\kappa )}^{-1-l}} \\ \end{array}\end{eqnarray} \tag{ A.18 }$$

$$\begin{eqnarray}&&=-{{r}^{M-1}}+{{p}_{M-2}}(r)+\frac{{{s}^{M}}}{s-r+{\rm i}\kappa },\end{eqnarray} \tag{ A.19 }$$

where ${{p}_{M-2}}(r)$ is a polynomial of the degree $M-2$, we establish the identity

$$\begin{eqnarray}&&\frac{{{\partial }^{M-1}}}{\partial {{r}^{M-1}}}{{\left[ \frac{{{(r-{\rm i}\kappa )}^{M}}}{s-r+{\rm i}\kappa } \right]}_{r=0}}=(M-1)!\left[ \frac{{{s}^{M}}}{{{(s+{\rm i}\kappa )}^{M}}}-1 \right],\end{eqnarray} \tag{ A.20 }$$

and hence

$$\begin{eqnarray}&&\begin{array}{l} \frac{{{\partial }^{M-1}}}{\partial {{r}^{M-1}}}{{\left[ \frac{{{(r-{\rm i}\kappa )}^{M}}}{(s-r+{\rm i}\kappa )(2\delta -s-r+{\rm i}\kappa )} \right]}_{r=0}} \\ =\frac{(M-1)!}{2(s-\delta )}\left[ \frac{{{(s-2\delta )}^{M}}}{{{(s-2\delta -{\rm i}\kappa )}^{M}}}-\frac{{{s}^{M}}}{{{(s+{\rm i}\kappa )}^{M}}} \right]. \\ \end{array}\end{eqnarray} \tag{ A.21 }$$

With its help we further transform (A.17) and obtain the contribution (18).

At $\delta =0$ equation (18) amounts to

$$\begin{eqnarray}\begin{array}{rcl} i\mathcal{T}_{M}^{(2)}(d,d^{\prime} ) & = & -\frac{\kappa {{{\rm e}}^{-\kappa D/2}}}{(M-1)!}\frac{{{\partial }^{M-1}}}{\partial {{s}^{M-1}}}{{\left\{ {{{\rm e}}^{-{\rm i}sD}}{{s}^{M-1}}\left[ 1-\frac{{{(s-{\rm i}\kappa )}^{M}}}{{{(s+{\rm i}\kappa )}^{M}}} \right] \right\}}_{s=0}} \\ {} & = & -2\kappa {{{\rm e}}^{-\kappa D/2}}\frac{1-{{(-1)}^{M}}}{2}, \\ \end{array}\end{eqnarray} \tag{ A.22 }$$

which explicitly proves the parity effect (23), (24).