Two-photon paired solitons supported by medium polarization

We derive for the first time fundamental equations that describe soliton spatial profiles consisting of two-photon mode fields and macroscopic polarization of medium. Numerical solutions of this basic equation are presented to give single and multiple soliton chains taking an example of para-H$_2$ $v=0 \leftrightarrow 1$ vibrational system. Although effects of dissipative relaxation are included in the general form for two-level system, the existence of static soliton-condensate is established. Stability analysis of soliton solutions is given. Arguments given here support that soliton-condensates thus derived are final remnants of time evolution of macro-coherent two-photon process. These soliton-condensates are expected to be important in the proposed neutrino mass spectroscopy using atoms.

I Introduction Solitons of macroscopic size which are made of dynamical electromagnetic fields supported by medium polarization are of great interest both from points of fundamental physics and applications. Basic entity when medium atoms are excited by the usual electric dipole transition is known as polariton [1], [2]. The study of solitons made of the dipole field has a long history since the early works [3]. Polariton type of solitons are expected to play fundamental roles in a new type of lasing material and quantum information, and there have been remarkable experimental progress on polariton condensates [4], [5].
We investigate in the present work solitons of a different kind; the system of medium atoms that are coupled to two-photon mode. The transiton between two relevant atomic levels are E1 forbidden like J = 2 → 0, ∆P = + transitions, but E2 or two-photon allowed. Time evolving macro-coherent dynamical process |e ↔ |g + γγ has been studied in our previous works [6], [7] and existence of explosive two-photon events of characteristic features have been established. This dynamical process is termed paired superradiance (PSR). Highly entangled two-photon states that emerge in PSR may have an advantage in some aspects of applications over the usual type of polariton. Moreover, at the fundamental level PSR dynamics is closely related to the proposed macro-coherent radiative emission of neutrino pair (RENP) |e → |g +γ +νν [7], [8] when an applied static electric field induces admixture of different parity state appropriate for RENP. Once solitons related to PSR are formed, they may become an ideal target state of RENP, because these solitons are stable remnants and do not emit QED background two-photon pairs from target ends. Many advantages of atomic RENP for future neutrino physics have been reviewed in [7].
In the present work we derive for the first time correct form of soliton profile equations related to PSR. Effects of phase-coherence relaxation and population decay are both included in this derivation. By numerical calculations we show that both single soliton and multiple soliton chains are derived by this profile equation. Results are surprising in the sense that these stable soliton solutions exist despite of the most general form of dissipative terms in the two-level system. Furthermore, the stability of soliton profiles is analyzed. Significance of soliton formation in the neutrino mass spectroscopy is in enhancing the signal to the background raito and in providing a stable experimental means. For this experimental purpose a steady switching mechanism between PSR and RENP modes has to be devised.
The rest of the present work is organized as follows. In Section II, after a brief review of the modified master equation of PSR dynamics we derive the static profile equation for field components after the Bloch vector is solved in terms of field components. With the symmetry of PSR system fully taken into account we reduce the effective field degrees of freedom to two real fields, which may be taken as two counter-propagating field strengths. A classical mechanical analogue of particle in two space dimensions is emphasized. The acting force is not conservative and dissipation plays important roles in deriving multiple soliton-condensates. In Section III the soliton profile equation thus derived is numerically solved under the condition of infinitely long target. There exist both single and multiple soliton chain solutions. We present numerically obtained profiles of soliton-condensates for para-H 2 vibrational transition v = 0 ↔ 1, a system strictly E1 forbidden with the main de-excitation mode being two-photon process. In Section IV the finite size of targets is treated. In Section V we discuss stability of soliton solutions and its stability under a condition is shown. Under other conditions static solutions may be understood as related to explosive PSR events [6]. In Appendix we give some technical details on the master equation and input data of para-H 2 .
Throughout this work we use the natural unit of = c = 1.

II Derivation of soliton profile equation
Since we assume that the target is irradiated by excitation and trigger lasers in one and its counterpropagating direction, the development of medium polarization occurs only in a direction taken here as x axis. We further assume that target atoms are distributed uniformly and take the continuous limit of target atom distribution. A further simplification we adopt for the moment is to take an infinitely large target.
The main part of two-photon interaction with two level atoms is described by an effective hamiltonian in 1+1 dimensional field theory: where c i (x) are wave functions for atoms in the level |i and α are 2×2 matrix given by product of two dipole transition elements d pe , d pg to an upper level |p coupled to two lower levels |e , |g by E1 dipole moments [6]. The hamiltonian is given in an interaction picture, hence the off-diagonal elements contain oscillating functions e ±iǫeg which may be eliminated by slowly varying envelope approximation (SVEA) done later. This hamiltonian is derived by taking the Markov approximation and eliminating |p . Numerical parameters of α matrix elements are illustrated for the example of para-H 2 vibrational system of v = 1 ↔ 0 in Appendix. In a typical situation the field E(x) = E R (x) + E L (x) is decomposed into a sum of right and left moving waves of different frequencies whose sum coincides with the energy difference, ω 1 + ω 2 = ǫ eg . From this hamiltonian one may derive the Bloch equation for bilinear forms of wave functions, namely the density matrix elements, R i = ψ † σ i ψ (with the normalization of the number density n for ψ † ψ). We include the most general dissipation term given by a rigorous analysis [9] for the two-level system, which agrees with description in terms of T 1 , T 2 due to Bloch. The medium polarization P is calculated from this density matrix element and inserted in the right hand side (RHS) of Maxwell equation, ( The resulting system of non-linear partial differential equations for dynamical variables, R i , E i , is called Maxwell-Bloch equation and constitutes the master equation for PSR dynamics after matched phase components are projected out. The master equation after reduction to the first order differential equations with slowly varying envelope approximation (SVEA) is given in [7].
We seek static solutions of the Maxwell-Bloch equation. An important departure taken here from [7] is that we rely on second order differential equations of space-time derivatives. This is because experience in known soliton solutions such as the real field Goldstone model [10] shows that the second order formalism is essential for derivation of solitons. With the spirit of SVEA we replace the derivative operation in RHS of Maxwell equation by −∂ 2 t → ω 2 i . The resulting new master equation for PSR is given in Appendix. We seek static remnant solutions of this master equation.
An important change in the second order formalism is in basic scale units of time, length and field strength. They are given by These differ from t * , E 2 * , the units defined in [7], in the important n (target number density) dependence. For the para-H 2 transition, ct 0 ∼ 0.3mm with n = 10 21 cm −3 .
The static version of Bloch equation is solved in the presence of inhomogeneous term (which reflects that the system ultimately approaches the ground state), resulting in the Bloch vector components written in terms of fields. This is given by an 3 × 3 matrix inversion for r = R/n, where A is an anti-symmetric 3 × 3 matrix, where e i , i = R, L are dimensionless fields in the unit of E 0 of right and left moving field envelopes and γ ± = (α ee ± α gg )/2α ge . In the para-H 2 example, γ − ∼ 0.64 , ǫ eg ∼ 0.52eV, and ǫ eg n ∼ 180GWmm −2 for n = 10 21 cm −3 . The resulting Bloch vector is inserted into RHS of Maxwell equation and the soliton profile equation is formulated, using real fields, e T = (ℜe R , ℑe R , ℜe L , ℑe L ) . The profile equation is given by in the 2 × 2 block-diagonal form. 4 × 4 matrix M is real and symmetric. The middle point of frequency It is useful to think of classical mechanical analogy by taking e as 4-vector of particle position, ξ as a fictitious time and RHS of (5) as a force. This force is not conservative, since Eigenvalues of the matrix M play important roles and it turns out that these are functions of R = |e R | 2 and L = |e L | 2 alone. We may thus expect that classical mechanics of this system is governed by these two variables and that there exist two-fold symmetry. This symmetry is associated with two conserved quantities, with X = ℜ(e R e L ), Y = ℑ(e R e L ). ′ indicates ξ derivative. Imposing these conservation makes it possible to write the soliton profile equation as ordinary second order differential equations in terms of right and left moving fluxes R, L: with two conserved quantities h, l given by We have neglected a small term in RHS of these equations that vanishes in the T 1 ≫ T 2 limit, which is valid in most targets.
III Single soliton and multiple soliton solutions Single soliton solutions are derived by solving a simpler version of the profile equation. We present results of soliton solutions using the para-H 2 vibrational transition.
The force-free region as defined by R ′′ = 0, L ′′ = 0 gives the velocity field R ′ , L ′ : Solutions of these simplified equations are illustrated in Fig(1). Profiles of right and left moving pulse energy densities start from vanishing derivative at a large distance from the center and end with vanishing derivative at another end. There are four types of solutions, depending on sign combinations in two spatially distant regions of ±∞. These single solitons have topological quantum numbers [10] and are expected to be stable [6]. We shall prove the stability below. Soliton sizes are of order several times ct 0 , hence are of macroscopic size, of several mm's in the example of Fig(1). In the approach here there is no limitation on the field strength, but it is actually limited by the stored energy density in the upper atomic level |e . We now turn to the full soliton profile equation, eqs(10), (11). At the point of vanishing field derivative, namely the end point of single soliton, the second derivative of field vanishes according to the full equation and beyond this point the second derivative of field can change its sign. This suggests that single soliton solutions of different type may be matched by solving the full equation.
Results for para-H 2 solutions are illustrated in Fig(2) and Fig(3). The quantity that gives directly the potentiality of both PSR and RENP is This quantity is also shown in presented figures. An example given in Fig(3) shows a large activity factor (17) of soliton chain. The largest value of this quantity is found to be 1 later in Section V. Corresponding to this largest value RENP rate (proportional to the activity factor) becomes maximal. For instance, in the Cs P 1/2 (1.39 eV) example studied in relation to the nuclear monopole pair emission [8] the rate becomes of order 10 4 Hz for the target number density n = 10 21 cm −3 (rate proportional to n 3 ) and for the target volume V = 10 2 cm 3 (proportional to V ).
IV Solitons for targets of finite length So far we assumed an infinitely large target. Finite size targets may be treated by devising appropriate boundary condition.  Let us first discuss the boundary condition for a single soliton solution in a finite size target. The appropriate boundary condition is to impose that both field and field derivative vanish along with the Bloch vector, r 1 = r 2 = 0, r 3 = −1 at the target end. These are satisfied by R = 0, R ′ = 0 or L = 0, L ′ = 0 in two field profile equation, (10) and (11). A practical way to satisfy this for targets of finite size would be to solve approximately the profile equation in the region near target ends. For simplicity we shall deal with the case of l = 0. For the right target end the approximate profile equation is which is solved as with K some integration constant and the target region of length L, −L/2 < ξ < L/2. Thus, the appropriate behavior near the right target end is R ′ /R ∼ 2/(ξ − L 2 ). A similar boundary condition is set up at the left target end.
Numerically, it is easier to impose this condition (for the end point behavior) at the middle point of targets simultaneously for right and left moving modes, because a periodic structure of condensate solitons are anticipated for multiple soliton-condensate. An example of finite size soliton thus obtained is illustrated in Fig(4). There is a tendency that the peak field energy density is larger for finite size soliton-condensates than for the infinitely large one, but it is difficult to accurately calculate the target length dependence of the peak power of soliton-condensates by means of numerical computations. The PSR activity factor, η−integrand (17), is of order, 10 −5 ∼ 10 −6 , but this quantity is roughly proportional to 1/τ 2 1 and may increase with smaller τ 1 . Thus, the estimate of this quantity depends much on target conditions.

V Stability analysis
The following general relation is important to stability analysis which is derived from the master equation in Appendix: This relation holds irrespective of any values of real fields e i , i = 1 ∼ 4. From its construction solutions of soliton profile equation are static and the right hand side (RHS) quantity of this equation vanishes for soliton-condensates. Solitons are thus characterized by This equation defines an elliptic manifold of two parameter (r 3 , |r T |) space characterized solely by two relaxation times T i . See Fig(5) for this manifold. Stability analysis may conveniently be divided depending on the sign of the RHS quantity. If the Bloch vector components (|r T |, r 3 ) fall within the elliptic manifold, the components r 3 , |r T | continue to vary with increasing magnitude of |r T | 2 + r 2 3 . If the Bloch components are initially outside the manifold, they continue to vary with the decreasing magnitude. Incidentally, we note a bound on the activity factor dη/dξ of two-photon process. The largest |r T | value T 2 /T 1 /2 is limited by the condition for relatxation times T 2 ≤ 2T 1 [9], which gives the largest allowed value 1/ √ 2 for |r T |. This may be combined with the largest stored energy in medium ǫ eg n, hence |e i | ≤ 1, to give the largest activity factor |r T | 2 (|e R | 2 + |e L | 2 ) = 1, as mentioned previously.
To examine more detailed behavior of perturbation, it is convenient to parametrize r 3 , r T in a way which makes it easy to identify the elliptic manifold, δ = 1 corresponds to the elliptic manifold of eq.(21), hence δ − 1 is a distance measure to the manifild, while δθ is a component parallel to the manifild.
The Bloch equation may be recast as differential equations for quantities, δ, θ, ϕ. Stability analysis can be made using the Bloch equation in terms of these variables near the elliptic manifold. A general equation that perturbation satisfies exist: to the leading T 2 /T 1 order, it is At the points of cos(2θ) = 0 solutions have a singular behavior, This means that time dependent perturbations maintain nodal points at cos(2θ) = 0. Otherwise, they follow bounded behavior as illustrated in Fig(6). We conjecture existence of breather solutions around soliton solutions. Perturbed soliton Figure 6: Example of perturbed soliton. The variable r 3 + 1/2 in solid black is a response to to a given perturbation of δ away from soliton manifold δ = 1 in dashed red.
We shall make a few further comments on the linear stability analysis. Equations linear in deviations from static profiles δ r(ξ, τ ), δ e(ξ, τ ) do not contain inhomogeneous terms (the unique inhomogeneous term ∝ −1/τ 1 does not contribute in perturbation equation). We shall assume that field perturbation does not exist and take fields that are given by a relation for soliton-condensates, e R e L = r T /(4τ 2 r 3 ) where r T is real and can now be either positive or negative. The coupled equation of perturbations δr T , δr 3 obeys We may assume the exponential behavior (δr T , δr 3 ) ∝ e λτ in order to locate the stability region in (r T , r 3 ). This gives a secular equation for eigenvalues λ: There are two possibilities for stability: two eigenvalues of λ are both real and negative, or are complex having negative real part. These two stability regions are plotted in Fig(7). It turns out that around the manifold curve of solitoncondensates only the exponential stability may be realized. This indicates that the exponential stability requires a range of parameters of r 3 > −0.70. Inside the stability region near the condensate boundary (21) a typical relaxation time towards soliton-condensate is of order,T 2 2 /T 1 , assuming all r i of order unity. This decay time is usually much shorter than the phase de-coherence time T 2 .
Outside the stability region solutions contain unstable objects. This may be understood as related to existence of explosive PSR events found in [6]. A momentarily static object may be created, followed by its fast decay emerging as large output fluxes of explosive events. The following long time behavior of perturbation makes clear relation to the self-induced transparency (SIT) or electromagnetically induced transparency (EIT). Stability analysis that contain interacting field perturbation and Bloch vector is beyond the scope of this work.

Long term stability
At t ≫ T i one may ignore the Bloch equation and can take δ r = 0 in the Maxwell equation. The homogeneous equation for δ e is where M 0 is the matrix when a static r of soliton field equation is inserted in the formula (5). To be precise, we note that this inserted static r is allowed to differ from the one obtained by solving the static profile equation in actual situations, since the precise time evolution of linear perturbation at τ ≤ τ 1 is unknown. The main feature of the long term stability is then governed by propagation effects of fields under the created medium Bloch vector r. Right after formation of soliton-condensates the medium polarization matrix given by M 0 has a simple structure. This matrix has eigenvalues of where r is the Bloch vector given by soliton solutions. Eigen-modes consist of combinations of fields given by an orthogonal mixture of δe i of 4 × 4 matrix, From the soliton profile equation one can determine the combination r 2 1 + r 2 2 , but not the ratio r 2 /r 1 that appears in this transformation. This ratio, hence the parameter ϕ in mode functions, is related to the conserved quantity, l = θ ′ RL/(R + L), and satisfies where θ = arctan L/R. When l = 0, the mixing parameter ϕ is independent of target position of solitons. Refractive indexes n may be defined by (∂ 2 t − ∂ 2 x )δ e = n 2 δ e for eigen-modes δ e. Hence in this case of space independent ϕ refractive indexes may be identified for each eigen-mode, giving two indexes, In the large relaxation limit τ 1 ≫ 1 which holds in most cases, C term here may be neglected compared to A term, and one has two refractive indexes, 1 and 1 + √ 2A. In the first case perturbation propagates freely, and even the second case describes a simple propagation effect. The first case is the analogue of self induced transparency in super-radiance [3].
These Bloch vector components are illustrated in Fig(8) in which the smaller refractive indexe is also shown. Two refractive indexes for right and left pulses can differ in principle, but in all our examples they are identical. In the example of Fig(2) the refractive index is nearly 1, with its deviation less than −0.003.
With this respect it would be useful to decompose perturbation δ e into eigen-modes of M 0 whose eigenvectors are described by the mixing parameter ϕ of eq.(30). We computed the ϕ spatial variation in the example of soliton target for the case given in Fig(10) of non-vanishing l. The mixing phase ϕ changes monotonically, δϕ ∼ −2.9ξ. The monotonic change implies that eigen-modes have periodic field change in soliton targets.
The larger refractive index is illustrated in Fig(11). Sample numerical computations show that perturbation do not grow, and follow faithfully input pulses.
In summary, we demonstrated for the first time that despite of existence of relaxation effects stable soliton-condensates exist in medium atomic system coupled to two-photon mode. Stability analysis indicates that these condensates are stable objects under some condition and furthermore that they are generically final remnants of time evolution in macro-coherent medium after trigger irradiation. The largest PSR activity factor of condensate-solitons given by the product of polarization and field energy density has been identified. Dynamical process of formation of soliton-condensates ought to be worked out by solving full space-time evolution equation.      The basic time and length unit t 0 and ct 0 are identified from the Maxwell equation (∂ 2 t − ∂ 2 x )E = −∂ 2 t P . RHS of the Maxwell equation has α ge as the coupling parameter and linear both in Bloch vector components R i = nr i and E. We thus identify the basic units and dimensionless time τ and length ξ as where t * is the unit defined in PTEP paper [7]. The field strength unit is identified from the Bloch equation, which is first order in time derivative and in RHS contains α ge E 2 R j . Hence the unit field strength E 2 0 and dimensionless quantity e i are defined by These modifications of basic units give scale dimensions when solitons are formed. In particular, length and time units are considerably smaller than those of [7]. The RENP rate unit ∝ E 2 0 L (L = the target length) is unchanged and coincides with previous results, because the scale change from [7] remains this combination invariant. The dynamical factor, in particular, after solitons are formed, is calculated using the soliton profile obtained by solving the profile equation written in the new units.