Multivariable Painleve'-II equation: connection formulas for asymptotic solutions

Multivariable Painle v ´ e-II equation: connection formulas for asymptotic solutions Nikolai A. Sinitsyn 1 , ∗ 1 Theor etical Division, Los Alamos National Laboratory , Los Alamos, New Me xico 87545, USA (Dated: October 2025) It is shown that a generalization of the Painlev ´ e-II equation (P-II) to a system of coupled equations with symmetry breaking terms is integrable. A Lax pair for this system is used to relate the asymptotic beha vior of the solutions at dif ferent infinities via an asymptotically exact WKB approach. The analysis relies on an exact solution of the quantum mechanical Demkov-Oshero v model (DOM). An application to the problem of unstable vacuum decay during a second order phase transition provides precise scaling of the number of excitations, including subdominant contributions. I. INTR ODUCTION Theoretical physics relies on a set of special functions, which are usually solutions of differential equations whose properties can be understood in great detail. Thus, through- out the 20th century , a physicist’ s standard toolbox included knowledge of solvable linear ordinary second order differen- tial equations. The commonly used property of such equations is the existence of connection formulas that relate the param- eters of the asymptotic solutions in different limits. T o wards the end of the 20th century , the connection formu- las for asymptotic solutions of certain second order nonlinear differential equations were also deriv ed. For e xample, the ho- mogeneous P-II u ′′ ( x ) = xu ( x ) − 2 u ( x ) 3 , (1) is one of the most frequently encountered nonlinear second or- der differential equations in theoretical physics. Connection formulas [ 1 , 2 ] relate the behavior of u ( x ) in the x → ±∞ limits. They were used to write integrals of P ainlev ´ e tran- scendents ov er x in terms of classical special functions [ 3 – 5 ]. Applications of such connection formulas hav e already been numerous, including in ultracold atoms [ 6 , 7 ], viscous flow [ 8 ], spectroscopy [ 9 ], quantum phase transitions [ 10 ], quan- tum information [ 11 ], plasma physics [ 12 , 13 ], liquid crystals [ 14 ], and stochastic processes [ 15 ]. In a seemingly different direction, there is a progress on solvable higher-or der linear systems, including the multistate Landau-Zener models [ 16 – 19 ] and solvable quantum quench models [ 20 – 26 ], which describe interacting quantum systems of potentially combinatorial complexity under e xplicitly time- dependent conditions. This raises the question of whether there e xist practically interesting higher -order systems of non- linear differential equations that are integrable and possess simple connection formulas for the asymptotic behavior of their solutions as well. The present Letter provides exact analytical connection for- mulas for a multiv ariable generalization of P-II. This nonlin- ear system is sho wn to be related to the Demk ov-Oshero v model (DOM) [ 10 , 27 ], which is one of the earliest known solvable generalizations of the quantum mechanical Landau- Zener model to an arbitrary number of interacting states with linearly time-dependent parameters. As the class of solv- able multistate Landau-Zener systems is kno wn to be large [ 28 , 29 ], its rele vance to solv ability of nonlinear systems sug- gests that the class of inte grable multi-v ariable generalizations of the Painlev ´ e equations with tractable asymptotic beha v- ior may be as large. Thus, the synergy between the two re- search directions in mathematical physics may enable analyt- ical treatment of more complex nonlinear problems. II. DEFINITION OF THE MODEL Let a system of n coupled nonlinear differential equations depend on n parameters e 1 < e 2 < . . . < e n : u ′′ 1 ( x ) = xu 1 ( x ) − 2 u 1 ( x ) n X k =1 u 2 n ( x ) + e 1 u 1 ( x ) , u ′′ 2 ( x ) = xu 2 ( x ) − 2 u 2 ( x ) n X k =1 u 2 n ( x ) + e 2 u 2 ( x ) , · · · (2) u ′′ n ( x ) = xu n ( x ) − 2 u n ( x ) n X k =1 u 2 n ( x ) + e n u n ( x ) . A shift, x → x − e 1 , sets e 1 = 0 , which will be assumed here. The system in Eq. ( 2 ) is integrable. The asymptotic behav- ior of the solution as x → ±∞ can be connected by analytic formulas. The first observation that enables this result is that the system in Eq. ( 2 ) is obtained as a consistency condition  ∂ H ∂ x − ∂ H 1 ∂ t  − i [ H , H 1 ] = 0 , (3) for two ( n + 1) -dimensional Hermitian matrices, H ( t, x ) and H 1 ( t, x ) , gi ven e xplicitly by 2 H = A − 4 tB + 2 C , H 1 =        t − u 1 · · · − u n − 1 − u n − u 1 − t 0 · · · 0 . . . . . . . . . . . . . . . − u n − 1 0 · · · − t 0 − u n 0 · · · 0 − t        , (4) where A =           4 t 2 + x − 2 n P k =1 u 2 k ( x ) 0 · · · · · · 0 0 − (4 t 2 + x + 2 e 1 ) + 2 u 2 1 2 u 1 u 2 · · · 2 u 1 u n . . . 2 u 2 u 1 − (4 t 2 + x + 2 e 2 ) + 2 u 2 2 . . . 2 u 2 u n . . . . . . . . . . . . . . . 0 2 u n u 1 ( x ) 2 u n u 2 · · · − (4 t 2 + x + 2 e n ) + 2 u 2 n           , B =      0 u 1 · · · u n u 1 0 · · · 0 . . . . . . . . . . . . u n 0 · · · 0      , C =      0 − iu ′ 1 · · · − iu ′ 3 ( x ) iu ′ 1 0 · · · 0 . . . . . . . . . . . . iu ′ n 0 · · · 0      , where u k ≡ u k ( x ) , u ′ k ≡ du k ( x ) dx . (5) Thus, if the functions u k ( x ) , k = 1 , . . . n , satisfy Eq. ( 2 ), then the Hamiltonians H and H 1 from Eq. ( 4 ) satisfy the consistency condition ( 3 ). Similar properties hav e been pre- viously used in studies of the P ainlev ´ e [ 2 ] and many other nonlinear dif ferential equations of mathematical physics [ 30 ]. Therefore, the fact that similar pairs of Hamiltonians can be found for higher-order nonlinear equations is not surprising. Howe ver , the consistenc y conditions do not guarantee that the WKB analysis that they enable leads to analytically tractable equations. The second observation that enables the following results is that such a WKB analysis for the system in Eq. ( 2 ) can be performed completely analytically . III. THE SYSTEM OF TWO V ARIABLES: MAIN RESUL TS Deriv ation of the connection formulas for the entire system in Eq. ( 2 ) is somewhat in volv ed. The goal of this brief Letter is, instead, to announce their e xistence and demonstrate the consequences of the integrability for the simplest nontri vial case of only two variables, n = 2 in Eq. ( 2 ), while leaving many details to future publications. For n = 2 , the system in Eq. ( 2 ) reduces to u ′′ 1 ( x ) = xu 1 ( x ) − 2 u 1 ( x )  u 2 1 ( x ) + u 2 2 ( x )  , (6) u ′′ 2 ( x ) = xu 2 ( x ) − 2 u 2 ( x )  u 2 1 ( x ) + u 2 2 ( x )  − eu 2 ( x ) , (7) where e > 0 . A. Connection formulas Standard perturbativ e analysis [ 31 ] fixes the asymptotic behavior of the solution for the system in Eqs. ( 6 ) and ( 7 ) as x → −∞ up to four parameters, α 1 , 2 and φ 1 , 2 , where α 1 , 2 > 0 . Namely , for x → −∞ : u 1 = α 1 ( − x ) 1 / 4 sin  2 3 ( − x ) 3 / 2 + (3 α 2 1 + 2 α 2 2 ) 4 ln( − x ) + φ 1  , u 2 = α 2 ( − x + e ) 1 / 4 sin[ 2 3 ( − x + e ) 3 / 2 + (3 α 2 2 + 2 α 2 1 ) 4 ln( − x ) + φ 2 ] . (8) The constant parameters α 1 , 2 and φ 1 , 2 will be referred to as initial amplitudes and phases, respectiv ely . As x → + ∞ , the trivial perturbation theory easily fix es only the leading x -dependent terms in the oscillation phases, while the sub-leading terms, starting with the terms behaving as ∝ ln x , are deri ved using the WKB approach described in Section IV . The solution is finally , i.e., as x → + ∞ , parametrized by two positiv e amplitudes, ρ and A , tw o phases, 3 FIG. 1. Numerical test of Eqs. ( 14 )-( 17 ). Solid curves correspond to theoretical predictions for the dependence on the equation parameter e ∈ (0 . 05 , 5) (a) of the action variables I 1 (brown) and I 2 (blue); (b) of the angle-related variables sin ϕ 1 (brown) and sin ϕ 2 (blue). Discrete points correspond to results of numerical simulations. Both for analytical predictions and numerical calculations the initial conditions were chosen to be α 1 = 0 . 9 , α 2 = 0 . 8 , φ 1 = π / 2 , φ 2 = π / 3 . Simulations were performed in the interv al x ∈ ( − 5000 , 5000) with a discretization step dx = 0 . 00001 . The algorithm is described in [ 32 ]. Both theory and numerical simulations agree on that σ = − 1 for all points here. ϕ 1 and ϕ 2 , and one sign parameter , σ = ± 1 : F or x → + ∞ : u 1 ( x ) = σ r x 2 + (9) σ ρ (2 x ) 1 / 4 cos " 2 √ 2 3 x 3 / 2 − 3 2 ρ 2 ln x + ϕ 1 # , u 2 ( x ) = σ A cos  √ ex − A 2 √ e 2 ln x + ϕ 2  . (10) Let me define the following combinations of the final and initial parameters: I 1 ≡ ρ 2 2 , I 2 ≡ A 2 √ e 2 , (11) p 1 = e − π α 2 1 , p 2 ≡ e − π α 2 2 , (12) including two phases: Φ 1 ≡ π 4 + arg Γ  i α 2 1 2  + φ 1 − 3 α 2 1 2 ln 2 + α 2 2 2 ln( e/ 4) , Φ 2 ≡ π 4 + arg Γ  i α 2 2 2  + φ 2 − 3 α 2 2 2 ln 2 + α 2 1 2 ln( e/ 4) . The main result of this article is the following connection for- mulas that relate the final parameters { σ, I 1 , I 2 , ϕ 1 , ϕ 2 } to the set of initial parameters { p 1 , p 2 , Φ 1 , Φ 2 } and the parameter of the equation, e : σ = sign [sin (Φ 1 )] , (13) I 1 = − 1 4 π ln 1 − p 2 1 p 2 2     1 + 1 − p 1 p 1 e 2 i Φ 1 + 1 − p 2 p 1 p 2 e 2 i Φ 2     2 ! , (14) ϕ 1 = − 3 π 4 − I 1 · 7 ln 2 + arg Γ (2 iI 1 ) − arg  1 + 1 − p 1 p 1 e 2 i Φ 1 + 1 − p 2 p 1 p 2 e 2 i Φ 2  , (15) I 2 = − 1 π ln  2 p p 2 p 1 (1 − p 1 ) | sin (Φ 1 ) |  − 2 I 1 , (16) ϕ 2 = 3 π 4 − 2 3 e 3 / 2 − I 2 ln(4 e 1 / 2 ) + arg [Γ( iI 2 )] − arg  e i Φ 2 + e − i Φ 2  p 1 + (1 − p 1 ) e 2 i Φ 1  . (17) B. T ests of connection formulas The first test for the validity of Eqs. ( 13 )-( 17 ) is the obser- vation that at α 2 = 0 , Eq. ( 7 ) has a tri vial solution u 2 ( x ) = 0 , which corresponds to I 2 = 0 and u 1 ( x ) satisfying Eq. ( 1 ). As expected, Eqs. ( 14 ) and ( 15 ) then reduce to the previously 4 FIG. 2. Numerical test of Eqs. ( 13 )-( 17 ). Solid curves correspond to theoretical predictions for the dependence on the angle φ 1 ∈ (0 , π ) (a,b) of the action variables I 1 (brown) and I 2 (blue), respecti vely; (c,d) of the angle-related v ariables sin ϕ 1 (brown) and sin ϕ 2 (blue), respecti vely . Discrete points correspond to results of numerical simulations. Both for analytical predictions and numerical calculations the initial conditions were chosen to be e = 1 , α 1 = 0 . 8 , α 2 = 0 . 6 , φ 2 = π / 2 . Simulations were performed in the interv al x ∈ ( − 5000 , 5000) with a discretization step dx = 0 . 00001 . Inset in (b) shows the test of Eq. ( 13 ) for the sign σ dependence on φ 1 (purple). known connection formulas for Eq. ( 1 ) listed in [ 2 ]. More generally , the connection formulas were tested nu- merically for x -ev olution with Eqs. ( 6 ) and ( 7 ). Simulations were performed using the approach described in [ 32 ], which treats x as a time variable. The e volution w as considered dur- ing x ∈ ( − 5000 , 5000) with initial parameters taken to be O (1) . After this interv al, a small piece of the trajectory was recorded and the final parameters were inferred by the best fit of this piece with Eqs. ( 9 ) and ( 10 ). The results are shown in Figs. 1 and 2 . Since the angles ϕ 1 , 2 are determined only up to an integer multiple of 2 π , to obtain smooth curves, sin ϕ 1 and sin ϕ 2 are plotted instead of the angles themselves. Apparently , despite the highly non- linear and sometimes singular dependence of the final set of parameters on the initial ones, the numerical results appear to be in perfect agreement with the analytical predictions. Many other similar tests (not sho wn) were performed for different values, e.g., of α 1 , 2 , with the same conclusion. 5 Finally , let me comment on corrections to the connection formulas. For large but finite x , they appear starting with the terms ∝ 1 / √ x , which decay relatively slo wly . They can be calculated using the same WKB approach, but the complete study of such terms appears to be dif ficult due to the large number of di verse relev ant contributions. Some of them are singular in the limit e → 0 . This is why the simulations were restricted to e > 0 . 05 – otherwise, a considerably larger in- tegration interval is needed in order to achieve a comparable agreement with theoretical predictions. Numerically , for the range of the tested parameters, with e > 0 . 05 , there were two specific corrections observed that dominated o ver other terms of the same order in x . These two corrections can be easily included into the connection formu- las, leading to noticeable accuracy improvements at large but finite x . First, the regular part of u 1 ( x ) in Eq. ( 9 ) is renormal- ized by interactions as σ r x 2 → σ r x − 2 u 2 ( x ) 2 2 , (18) where u 2 ( x ) is taken from Eq. ( 10 ). This correction vanishes as ∝ 1 / √ x , but it is only a factor ∼ 1 /x 1 / 4 smaller than the leading oscillatory term in Eq. ( 9 ). Therefore, its ef fects are well visible during simulations as a slowly decaying modula- tion of oscillations in u 1 ( x ) . Second, there is a correction to ϕ 1 in Eq. ( 15 ) that follows from the renormalization of the energy of the state | 0 ⟩ in the WKB approach applied to the e volution along the path P ∞ (see Section IV for details). It can be included in the connec- tion formulas by subtracting an additional term from the right hand side of Eq. ( 15 ): ϕ 1 → ϕ 1 − 2 π I 2 r e 2 x . (19) Although well suppressed, small effects of this correction were visible even for the longest simulations. This is at- tributed to its growth with e and the existence of singularities that make I 2 anomalously large at specific initial conditions. IV . DERIV A TION OF CONNECTION FORMULAS A. Strategy Consistency conditions ( 3 ) allo w one to define a state vec- tor , | Ψ( t, x ) ⟩ , as a solution of the tw o-time Schr ¨ odinger equa- tion [ 30 ] i d | Ψ ⟩ dt = H ( t, x ) | Ψ ⟩ , (20) i d | Ψ ⟩ dx = H 1 ( t, x ) | Ψ ⟩ , (21) where H and H 1 are 3 × 3 matrices gi ven by Eq. ( 4 ) for n = 3 . The ev olution operator U P along an arbitrary path P in the FIG. 3. An integration path P (dashed arro w) with x 0 → −∞ and t ∈ ( − t 0 , t 0 ) , where t 0 ≫ x 0 , is deformed into the path P ∞ , such that the horizontal segment of P ∞ lies at x 0 → + ∞ and t ∈ ( − t 0 , t 0 ) (dotted arrows). This deformation does not change the evolution operator in Eq. ( 22 ), since the initial and final points of the path remain unchanged. The vertical legs of P ∞ hav e t = ± t 0 → ±∞ , which makes the ev olution along them adiabatic. two-time space ( t, x ) can be written as a path-ordered expo- nent: U P = T P e − i R P { H dt + H 1 dx } , (22) where T P is the path ordering operator , such that factors cor- responding to earlier points along P appear further to the right inside the product of ev olution exponents e − i { H dt + H 1 dx } . Equation ( 22 ) can be written as U P = T P e R P A · d τ , (23) where τ ≡ ( t, x ) is the two-time point and A ≡ ( − iH , − iH 1 ) is a non-Abelian field. This field is flat (has zero curvature) [ 18 ], so the result of the evolution in Eq. ( 23 ) depends only on the endpoints, which are ( − t 0 , − x 0 ) and ( t 0 , − x 0 ) in Fig. 3 , but does not depend on the choice of the path connecting these points. Thus, the ev olution ov er the path P shown in Fig. 3 can be calculated in two ways. First, one can consider t changing in the interval t ∈ ( − t 0 , t 0 ) at fixed negati ve x = − x 0 . The ev olution operator then takes the form U P = T t e − i R t 0 − t 0 H ( t, − x 0 ) dt . (24) Alternativ ely , U P can be obtained as an evolution operator along the path P ∞ in Fig. 3 . The vertical legs of P ∞ corre- spond to e volution at fixed t . For example, ev olution along the left vertical le g is described by the operator U v l = T x e − i R x 0 − x 0 H 1 ( − t 0 ,x ) dx , where the corner points of the paths are tak en as t 0 , x 0 → ∞ , with the limit t 0 → ∞ taken first. 6 FIG. 4. (a) Eigen v alues of the Hamiltonian ( 26 ) as functions of t = p | x | τ at large negativ e x (blue curves). Here, x = − 16 ; α 1 = α 2 = 0 . 1 , φ 1 = π / 3 , φ 2 = 3 π / 4 . The values of u 1 ( x ) and u 2 ( x ) are approximated by asymptotic formulas in Eq. ( 8 ). The Hamiltonian ( 26 ) is written in the basis of diabatic states. The labels 0 , 1 , 2 define the conv ention for the diabatic state indices. The diabatic states coincide asymptotically with the eigenstates as t → ±∞ . The red and green arrows show two semiclassical trajectories that originate on diabatic lev el 0 as t → −∞ and terminate at le vel 1 as t → + ∞ . The dashed circle encloses the beha vior of the energy lev els near the pseudo-time point t = p | x | / 2 . (b) Near t = p | x | / 2 , the dynamics is described by the exactly solvable DOM for three levels that cross linearly at fixed couplings g + 1 and g + 2 . The diagram of diabatic lev els shows only the time dependence of diagonal elements of the Hamiltonian, which in the DOM are straight lines in the time-energy plot. (c) Eigen values of the Hamiltonian H ( t, x ) + x 3 / 2 (4 τ 2 + 1) ˆ 1 , as functions of t = √ xτ at large positiv e x (blue). Here, x = 10 , A = 0 . 2 , ρ = 0 . 12 , e = 1 . 2 . A term proportional to ˆ 1 was added to the Hamiltonian to expose the spectrum better . The values of u 1 ( x ) and u 2 ( x ) are approximated by asymptotic formulas in Eqs. ( 9 ) and ( 10 ), respectiv ely . The red arrows show the unique path that connects diabatic lev els 0 and 1 . At large positi ve x , it passes through two av oided crossing points near times t 0 and t + . Let | k ⟩ , where k = 0 , 1 , 2 , be the diabatic states, in which basis the Hamiltonians H and H 1 are written in Eq. ( 4 ). For t 0 → ∞ , these states coincide with the eigenstates of H 1 ( − t 0 , x ) . The evolution along the left and right legs of P ∞ is adiabatic, resulting only in phase factors e ± 2 it 0 x 0 . Thus, only the horizontal segment of P ∞ contributes to the inter- state transition amplitudes. Both the horizontal part of P ∞ and the path P , describe unitary quantum ev olution over a pseudo-time t ∈ ( − t 0 , t 0 ) at large positiv e and large negativ e x , respecti vely . Since the operator H ( t, x ) depends on u 1 , 2 ( x ) , the relations between the asymptotic values of u 1 , 2 ( x ) as x → ±∞ are obtained by equating U P to U P ∞ . Since the inter -state transition am- plitudes are determined by the horizontal paths, at either large positiv e or large negati ve x , both e volution operators can be ev aluated using the WKB approach, which becomes exact in the limits x → ±∞ . B. WKB analysis in the limit x → −∞ In this case, there are two resonant regions for t ∈ ( − t 0 , t 0 ) near t = ± p | x | / 2 , where all nonadiabatic transitions occur . For large | x | , these resonances are well separated. By rescal- ing time t = τ p | x | , (25) the Schr ¨ odinger equation ( 20 ) is transformed to i d dτ | Ψ ⟩ = H ( τ ) | Ψ ⟩ , where H ( τ ) =   | x | 3 / 2 (4 τ 2 − 1) − 2 √ x ( u 2 1 ( x ) + u 2 2 ( x )) − 4 | x | τ u 1 − 2 i p | x | u ′ 1 − 4 | x | τ u 2 − 2 i p | x | u ′ 2 − 4 | x | τ u 1 + 2 i p | x | u ′ 1 −| x | 3 / 2 (4 τ 2 − 1) + 2 p | x | u 1 ( x ) 2 2 p | x | u 1 ( x ) u 2 ( x ) − 4 | x | τ u 2 + 2 i p | x | u ′ 2 2 p | x | u 1 ( x ) u 2 ( x ) −| x | 3 / 2 (4 τ 2 − 1) + 2 p | x | ( e + u 2 ( x ) 2 )   . (26) The time-dependent spectrum of this Hamiltonian for lar ge negati ve x is sho wn in Fig. 4 (a). Near τ = ± 1 / 2 , additional simplifications follow . Let g − 1 , 2 ≡ 2( | x | u 1 , 2 − i p | x | u ′ 1 , 2 ) . After shifting the time v ariable to τ m = τ + 1 / 2 , the Hamilto- nian in the vicinity of τ m = 0 reduces to the three-state DOM Hamiltonian with one level crossing linearly two parallel lev- els: H − 1 / 2 ( τ m ) =   − 4 | x | 3 / 2 τ m g − 1 g − 2 ( g − 1 ) ∗ 4 | x | 3 / 2 τ m 0 ( g − 2 ) ∗ 0 4 | x | 3 / 2 τ m + 2 p | x | e   . 7 Similarly , consider the time interval near τ = 1 / 2 and keep only the relev ant x -dependent terms. Introduce the couplings g + 1 , 2 = − 2( | x | u 1 , 2 + i p | x | u ′ 1 , 2 ) . After a shift of time as τ p = τ − 1 / 2 , the Hamiltonian near τ p = 0 is again the DOM: H 1 / 2 ( τ p ) =   4 | x | 3 / 2 τ p g + 1 g + 2 ( g + 1 ) ∗ − 4 | x | 3 / 2 τ p 0 ( g + 2 ) ∗ 0 − 4 | x | 3 / 2 τ p + 2 p | x | e   . The diabatic le vel diagram for H 1 / 2 is shown in Fig. 4 (b). The two parallel le vels may be close to one another , so that all three states interact simultaneously near the a voided cross- ings. Ne vertheless, this model is solv able analytically . The DOM exact solution for the scattering amplitudes coin- cides with the prediction of the independent crossing approx- imation, in which each encountered diabatic level crossing is treated as an independent two-state Landau-Zener transi- tion between the crossing ener gy levels, followed by adiabatic ev olution between successiv e lev el crossings [ 10 ]. Therefore, the entire calculation of the elements of the ev olution matrix U P reduces to identifying semiclassical trajectories connect- ing the initial and final states and then finding the correspond- ing quantum amplitudes using the independent crossing ap- proximation. An e xample of two semiclassical trajectories connecting the initial state | 0 ⟩ to the final state | 1 ⟩ is shown by the red and green arro ws in Fig. 4 (a). The transition ampli- tudes between diabatic states are obtained using the scattering matrix for the Landau-Zener model (listed, e.g., in Chapter 5 of Ref. [ 10 ]) to describe transitions between pairs of states at any diabatic le vel intersection. Between such pairwise scatter - ings, the system undergoes adiabatic ev olution with energies that can be estimated within second-order quantum perturba- tion theory in the small parameter 1 / | x | . This WKB approach is equiv alent to the textbook indepen- dent crossing approximation, which is widely used in theoret- ical physics and described, e.g., in Section 5.4 of [ 10 ]. This version of WKB was used in the deriv ation of transition am- plitudes in multistate Landau-Zener models [ 28 ]. The advan- tage of this approach over the WKB analysis commonly used in studies of Painle v ´ e equations [ 2 ] is that the pseudo-time t is ne ver treated as a complex v ariable, so there is no need to analyze the full Stokes phenomenon to obtain the evolution matrix U P . In the present case, the full matrix U P is not required. It is sufficient to assume that initially the system is in state | 0 ⟩ , which has the highest energy as t → −∞ , thus leading to the deriv ation of transition amplitudes U P, 00 , U P, 10 , and U P, 20 . C. WKB analysis as x → + ∞ This limit is simpler because, as illustrated in Fig. 4 (c), the lev el crossing pattern leads only to three elementary av oided crossings between pairs of lev els. These crossings are well separated both in energy and in time t . For e xample, the times of diabatic lev el crossings are giv en by t 0 = 0 and t ± ≈ ± x/ (4 √ e ) . Therefore, each crossing is described by the scattering matrix for the two-state Landau-Zener model. Moreov er , for the initial conditions with only lev el 0 popu- lated, the av oided crossing at t − is irrelev ant, so only a single trajectory connects the initial state | 0 ⟩ to any final diabatic state. Therefore, path interference does not occur . Figure 4 (c) illustrates this property by sho wing the only semiclassical tra- jectory connecting lev els 0 and 1 . T o calculate the scattering amplitudes using the asymptot- ically exact WKB, the functions u 1 , 2 ( x ) as x → + ∞ were assumed to take the form given in Eqs. ( 9 ) and ( 10 ) with un- known parameters, including at logarithmic phase contribu- tions. Near each av oided crossing, the couplings were ex- pressed in terms of these unknown parameters, and the scat- tering amplitudes were estimated using the known scattering matrix for the Landau-Zener model. After obtaining the amplitudes of U P,k 0 , where k = 0 , 1 , 2 , separately for P and P ∞ , the WKB results were compared order by order in powers of x and t 0 , down to O (1) contribu- tions to the phase and amplitude. T erms that depended on t 0 and on high powers of x on both sides then canceled, leaving the connection formulas written in Eqs. ( 13 )-( 17 ) and fixing the logarithmic x -dependence in Eqs. ( 9 ) and ( 10 ). V . EXCIT A TIONS AFTER UNST ABLE V A CUUM DECA Y The system in Eqs. ( 6 ) and ( 7 ) w as previously introduced in [ 32 , 33 ], where it described a passage through a phase transi- tion in scalar field theory . It can be interpreted as Hamiltonian dynamics, in which x is associated with physical time, t (not to be confused with the pseudo-time introduced in the previ- ous section), and u 1 ( t ) and u 2 ( t ) are coordinates of a particle with two degrees of freedom. Let X = ( u 1 , u 2 ) , and the cor- responding time-dependent Hamiltonian is H ( t ) = P 2 2 − t X 2 2 + X 4 2 + eu 2 2 2 , (27) where P = ( P 1 , P 2 ) is the vector of momenta canonically conjugated to u 1 and u 2 , and where X 2 ≡ u 2 1 + u 2 2 . The system in Eqs. ( 6 ) and ( 7 ) follo ws from Ne wton’ s equations of motion with the Hamiltonian ( 27 ). In the limits t → ±∞ , the potential energy for this classi- cal motion is that of a harmonic oscillator, so the amplitudes α 1 , 2 are related to the asymptotically conserv ed, as t → −∞ , adiabatic in variants [ 10 , 12 ]. J 1 , 2 ≡ α 2 1 , 2 2 . Similarly , I 1 and I 2 from Eq. ( 11 ) are identified as asymp- totically conserved adiabatic inv ariants as t → + ∞ . Then, φ 1 , 2 and ϕ 1 , 2 are the angles canonically conjugated to, re- spectiv ely , J 1 , 2 and I 1 , 2 . 8 According to [ 32 , 33 ], this classical dynamics appears as a saddle-point solution for a quantum field theory describ- ing a time-dependent passage through a second-order quan- tum phase transition. The dif ference between J 1 , 2 and I 1 , 2 captures inevitable nonadiabatic effects as the classical phase space trajectory passes through a separatrix. Physically , this ev ent corresponds to crossing the critical point of a phase tran- sition in the corresponding field theory [ 33 ]. Semiclassically , the actions are related to the numbers of produced quasiparti- cles: I 1 / ℏ for the number of Higgs bosons and I 2 / ℏ for the number of Goldstone bosons. The latter acquire a small mass due to the symmetry breaking term ∝ e in Eq. ( 27 ). The anal- ysis in [ 32 ] did not provide explicit connection formulas but made sev eral observations that can no w be verified. Thus, a useful feature of Eqs. ( 14 ) and ( 16 ) is that I 1 and I 2 depend on e only through the angles Φ 1 and Φ 2 , which them- selves depend linearly on the initial phases φ 1 and φ 2 . There- fore, ⟨ I 1 ⟩ and ⟨ I 2 ⟩ are independent of e , where the brackets denote averaging ov er a uniform distribution of φ 1 and φ 2 . This in variance of the averaged actions was conjectured in [ 32 ] based on less rigorous analytical arguments and numeri- cal simulations, and is now confirmed. Next, the solution in Eqs. ( 13 )-( 17 ) e xhibits asymmetry am- plification, which was the main finding in [ 32 ]. No matter how small the symmetry breaking parameter e is, the dynam- ics for u 1 ( t ) and u 2 ( t ) strongly dif fer from one another as t → + ∞ . While u 1 ( t ) has a monotonically growing com- ponent, the asymptotic solution for u 2 ( t ) oscillates around zero. Note that if e changes sign, the asymptotic solutions in Eqs. ( 9 ) and ( 10 ) are obtained by switching indices 1 and 2 . The solution at e = 0 reduces to the standard P-II solution, which is unstable in the present context. The study in [ 32 ] also addressed the following question: Giv en that initially the system is nearly at the quantum v ac- uum state, which corresponds, in dimensionless variables here, to J 1 = J 2 ≡ J ≪ 1 , (28) what are the v alues of the adiabatic in variants as t → + ∞ ? A veraging over the initial phases φ 1 and φ 2 was required. Using that R 2 π 0 ln(2 | sin(Φ) | ) d Φ = 0 , av eraging Eq. ( 17 ) ov er φ 1 , 2 at conditions in Eq. ( 28 ) giv es that ⟨ I 2 ⟩ + 2 ⟨ I 1 ⟩ ≈ 1 2 π ln  1 2 π J  . (29) Similarly , Eq. ( 14 ) with Eq. ( 28 ) lead to ⟨ I 1 ⟩ ≈ 1 4 π  ln  1 2 π J  − c 1  , (30) where c 1 = 1 π 2 Z Z π 0 ln[4 − 2 cos(2Φ 1 ) − 2 cos(2Φ 2 )] d Φ 1 d Φ 2 ≈ 1 . 166 . (31) Substituting Eqs. ( 30 ) and ( 31 ) into Eq. ( 29 ) leads also to ⟨ I 2 ⟩ ≈ c 2 2 π , c 2 = c 1 ≈ 1 . 166 . (32) Expressions for ⟨ I 1 ⟩ and ⟨ I 2 ⟩ in Eqs. ( 30 ) and ( 32 ) were conjectured in [ 32 ] for J ≪ 1 without an analytical deriv a- tion of c 1 , 2 , while numerical simulations in [ 32 ] estimated that c 1 ≈ 1 . 14 and c 2 ≈ 1 . 19 , with an uncertainty in the last sig- nificant digit. This agrees well with the analytically deri ved values of these coefficients in Eq. ( 32 ). Curiously , the value of the double phase integral in Eq. ( 31 ), ≈ 1 . 166 , is the fa- mous square-lattice spanning tree constant [ 34 ]. VI. CONCLUSION When a system of nonlinear equations emerges from con- sistency conditions ( 3 ) for matrix operators, the nonlinear problem is ef fecti vely described by a linear one. The dif fi- culty for higher-order systems of nonlinear equations often lies in the lack of a fully analytical solution of the associated linear model, even when WKB methods are applicable. Ne v- ertheless, the class of known fully solvable multistate linear systems has been gro wing recently . Any of the newly solved linear systems may emerge in the WKB analysis enabled by the consistency conditions for a certain nonlinear system. The present article provides an example of such a situation, lead- ing to analytical connection formulas for asymptotic solutions of two coupled P-II equations with a minimally broken sym- metry . Author thanks Prof. Alexander Its for useful discus- sions. This work was carried out under the auspices of the U.S. DoE through the Los Alamos National Labora- tory , operated by Triad National Security , LLC (Contract No. 892333218NCA000001). ∗ nsinitsyn@lanl.gov [1] A. R. Its and A. A. Kapaev , Mathematics of the USSR-Izv estiya 31 , 193 (1988) . [2] A. S. Fokas, A. R. Its, A. A. Kapaev , and V . Y . Novok- shenov , P ainlev ´ e T ranscendents: The Riemann-Hilbert Ap- pr oach , Mathematical Surveys and Monographs (American Mathematical Society , 2006). [3] J. Baik, R. Buckingham, J. DiFranco, and A. Its, Nonlinearity 22 , 1021 (2009) . [4] P . K okocki, Studies in Applied Mathematics 144 , 504 (2020) . [5] T . J. Bothner and R. Buckingham, Communications in Mathe- matical Physics 359 , 223 (2017) . [6] A. P . Itin and P . T ¨ orm ¨ a, Phys. Re v . A 79 , 055602 (2009) . [7] V . G. Sadhasivam, F . Suzuki, B. Y an, and N. A. Sinitsyn, Nat. Commun. 15 , 10246 (2024) . [8] S.-Y . Lee, R. T eodorescu, and P . W iegmann, Physica D: Non- linear Phenomena 240 , 1080 (2011) , arXi v:1005.0369 [nlin.SI] . [9] S. Banerjee and N. A. Sinitsyn, Phys. Rev . E 108 , 034212 (2023) . 9 [10] N. A. Sinitsyn and V . L. Pokrovsky , Quasi-Adiabatic Ef fects: Intr oduction to Geometric Phases and Landau-Zener T ransi- tions , Oxford Graduate T exts (Oxford Univ ersity Press, Oxford, United Kingdom, 2026) hardcov er (Feb . 5 2026). [11] A. R. Its, F . Mezzadri, and M. Y . Mo, Communications in Math- ematical Physics 284 , 117 (2008) . [12] L. M. Zelenyi, A. I. Neishtadt, A. V . Artemyev , D. L. V ainchtein, and H. V . Malov a, Physics–Uspekhi 56 , 347 (2013) . [13] A. Neishtadt, A. Artemye v , and D. T uraev , Chaos 29 , 051104 (2019) . [14] M. G. Clerc, J. D. D ´ avila, M. Ko walczyk, P . Smyrnelis, and E. V idal-Henriquez, Calculus of V ariations and Partial Differ- ential Equations 56 , 93 (2017) . [15] T . Bothner and K. Liechty , Nonlinearity 26 , 1449 (2013) . [16] N. A. Sinitsyn and V . Y . Chernyak, J. Phys. A: Math. Theor . 50 , 255203 (2017) . [17] N. A. Sinitsyn and F . Li, Phys. Re v . A 93 , 063808 (2016) . [18] N. A. Sinitsyn, E. Y uzbashyan, V . Y . Chernyak, A. Patra, and C. Sun, Phys. Re v . Lett. 120 , 190402 (2018) . [19] V . Y . Chernyak, N. A. Sinitsyn, and C. Sun, Phys. Rev . B 100 , 224304 (2019) . [20] J. Lin and N. A. Sinitsyn, J. Phys. A: Math. Theor . 47 , 175301 (2014) . [21] F . Li, V . Y . Chernyak, and N. A. Sinitsyn, Phys. Rev . Lett. 121 , 190601 (2018) . [22] B. Y an and N. A. Sinitsyn, Nat. Commun. 13 , 2212 (2022) . [23] S. Barik, L. Bakker, V . Gritse v , J. Min ´ a ˇ r, and E. A. Y uzbashyan, arXiv preprint (2025), arXi v:2507.10856 . [24] P . R. Pasnoori and E. A. Y uzbashyan, Exact many-body wave- function of the Kondo model with time-dependent interaction strength (2025), arXiv:2509.05640 [cond-mat.str -el] . [25] P . R. Pasnoori, Phys. Re v . B 112 , L060409 (2025) . [26] P . R. Pasnoori, Phys. Re v . B 113 , 094514 (2026) . [27] Y . N. Demkov and V . I. Osherov , Zh. Eksp. T eor . Fiz. 53 , 1589 (1967), sov . Phys. JETP 26 , 916 (1968), English translation. [28] V . Y . Chernyak, N. A. Sinitsyn, and C. Sun, J. Phys. A: Math. Theor . 51 , 245201 (2018) . [29] V . Y . Chernyak, F . Li, C. Sun, and N. A. Sinitsyn, J. Phys. A: Math. Theor . 53 , 295201 (2020) . [30] L. D. Faddee v and L. A. T akhtajan, Hamiltonian Methods in the Theory of Solitons (Springer , New Y ork, 1987). [31] L. D. Landau and E. M. Lifshitz, Mechanics, Course of Theo- r etical Physics V ol. 1 , 3rd ed. (Perg amon, Oxford, 1980). [32] B. T yagi, F . Suzuki, V . Y . Chernyak, and N. A. Sinitsyn, Phys. Rev . A 111 , 032205 (2025) . [33] F . Suzuki, R. K. Malla, and N. A. Sinitsyn, Phys. Rev . A 112 , 013307 (2025) . [34] G. M. V iswanathan, Phys. Re v . E 95 , 062138 (2017) .