A solvable model for symmetry-breaking phase transitions

Analytically solvable models are benchmarks in studies of phase transitions and pattern-forming bifurcations. Such models are known for phase transitions of the second kind in uniform media, but not for localized states (solitons), as integrable equations which produce solitons do not admit intrinsic transitions in them. We introduce a solvable model for symmetry-breaking phase transitions of both the first and second kinds (alias sub- and supercritical bifurcations) for solitons pinned to a combined linear-nonlinear double-well potential, represented by a symmetric pair of delta-functions. Both self-focusing and defocusing signs of the nonlinearity are considered. In the former case, exact solutions are produced for symmetric and asymmetric solitons. The solutions explicitly demonstrate a switch between the symmetry-breaking transitions of the first and second kinds (i.e., sub- and supercritical bifurcations, respectively). In the self-defocusing model, the solution demonstrates the transition of the second kind which breaks antisymmetry of the first excited state.


INTRODUCTION
The topic: spontaneous symmetry breaking in nonlinear systems Dynamics of collective excitations in physical systems is determined by the interplay of the underlying diffraction or dispersion, nonlinear self-interactions of the fields or wave functions, and potentials acting on the fields.In this context, it is commonly known that the ground state (GS) of linear systems reproduces the symmetry of the underlying potential, while excited states may realize other representations of the same symmetry [1].In particular, the wave function of a particle trapped in a symmetric double-well potential (DWP) is even, while the first excited state is odd.
While these basic properties are demonstrated by the linear Schrödinger equation, the dynamics of Bose-Einstein condensates (BECs) is governed, in mean-field approximation, by the Gross-Pitaevskii equation (GPE), which takes into regard interactions between particles, adding the cubic term to the Schrödinger equation for the single-particle wave function [2,3].The repulsive or attractive interactions are represented by the cubic term with the self-defocusing (SDF) or self-focusing (SF) sign.Essentially the same model is the celebrated nonlinear Schrödinger equation (NLSE), which governs the propagation of optical waves in nonlinear media [4] and finds other realizations, as the universal model to govern the interplay of the weak diffraction or dispersion and cubic SF nonlinearity [5].In optics, a counterpart of the trapping potential is the term in the NLSE which accounts for the waveguiding structure induced by a transverse profile of the refractive index.
The GS structure in models combining the DWP and SF nonlinearity follows the symmetry of the underlying potential only in the weakly nonlinear regime.A generic effect which occurs with the increase of the nonlinearity strength is the symmetry-breaking phase transition, which makes the GS asymmetric with respect to two wells of the DWP [6].This effect of the spontaneous symmetry breaking (SSB) implies, inter alia, that the commonly known principle of quantum mechanics, according to which GS cannot be degenerate [1], is no longer valid in the nonlinear models: obviously, the SSB gives rise to a degenerate pair of two mutually symmetric GSs, with the maximum of the wave function pinned to the left or right potential well of the underlying DWP.The same system admits a symmetric state coexisting with the asymmetric ones, but, above the SSB point, it does not represents the GS, being unstable against symmetry-breaking perturbations.
In systems with the SDF sign of the nonlinearity, the GS remains symmetric and stable, while the SSB transition breaks the antisymmetry of the first excited state (it is a spatially odd one, with precisely one zero of the wave function, located at the central point).The resulting state with the spontaneously broken antisymmetry keeps the zero point, which is shifted from the center to right or left.
The concept of the SSB in systems of the NLSE type with the SF nonlinearity was first proposed, in an abstract mathematical form, by Davies in 1979 [7].Another early realization of this concept was introduced in 1985 by Eilbeck, Lomdahl, and Scott, in the form of the "self-trapping model" [8].The latter work had actually initiated systematic studies of SSB phase transitions.
In optics, the SSB was observed experimentally in a photorefractive crystal with saturable SF nonlinearity and an effective DWP waveguiding structure [9].SSB was also predicted for photonic modes supported by a symmetrically designed plasmonic metamaterial [10].For the self-attractive BEC loaded into a DWP trap, the symmetry-breaking transition was elaborated in Refs.[11,12].In that context, tunnel-coupling oscillations between condensates trapped in two potential wells separated by a barrier represent the bosonic Josephson effect [13].Experimentally, the selftrapping of a stationary state with spontaneously broken antisymmetry in a self-repulsive condensate loaded in DWP, as well as Josephson oscillations in that setup, were reported in Ref. [14].
A ramification of the topic is SSB in dual-core systems, such as twin-core optical fibers, with the SF cubic nonlinearity acting in each core [15].In such systems, the interplay of the SF and linear coupling between the parallel cores gives rise to the SSB transition from the symmetric state to a spontaneously established one with unequal powers carried by the two cores.This type of the symmetry-breaking phase transition was studied in detail theoretically [16]- [21] and recently demonstrated experimentally in a twin-core fiber [22].In terms of the respective system of linearlycoupled NLSEs, the SSB transition is represented by the bifurcation which links symmetric and asymmetric solutions [23].Depending on the type of the intra-core nonlinearity and the wave form under the consideration (delocalized or self-trapped), the symmetry-breaking bifurcation may be of the supercritical (alias forward) or subcritical (backward) type.The corresponding bifurcations give rise to the destabilization of the symmetric state and creation of the pair of asymmetric ones at the SSB point, which go forward or backward as stable or unstable branches, respectively, following the variation of the SSB-driving nonlinearity strength.In the latter (subcritical) case, the unstable lower branches of the asymmetric states normally reverse into the stable forward-going upper ones at certain turning points (see, e.g., Fig. 7 below).As a result, stable asymmetric states emerge subcritically, at a value of the nonlinearity strength which is smaller than that at the SSB point.Accordingly, the system is bistable in the interval between the turning and SSB points, where the symmetric and upper asymmetric states coexist as stable ones.In terms of statistical physics, the super-and subcritical bifurcations are identified as symmetry-breaking phase transitions of the second and first kinds, respectively.In the latter case, the bistability corresponds to the hysteresis between the GS and overcooled or overheated phases.
Other varieties of optical SSB effects occur in dual-core laser setups combining the SF nonlinearity with gain and loss.The theoretical model of such setups is based on a pair of linearly coupled complex Ginzburg-Landau equations with the cubic-quintic nonlinearity [24].A spontaneously established asymmetric regime of the operation of a symmetric pair of coupled lasers was observed in Ref. [25].
The SSB phenomenology was also predicted in models with symmetric nonlinear potentials, induced by spatial modulation of the local SF or SDF coefficient [27].In optical waveguides, the modulation can be imposed by spatially inhomogeneous distributions of a resonant dopant, which gives rise to strong local nonlinearity [28].In experiments with BEC, the Feshbach resonance (FR) controlled by spatially nonuniform laser illumination of the condensate may be employed to build an effective nonlinearity landscape [29][30][31].Other techniques available to the experimental work with BEC make it possible to "paint" a necessary FR-induced nonlinear potential by a fast moving laser beam [32] or a spatial light modulator [33][34][35].

The model
The use of the nonlinear potential suggests a possibility to design experimentally feasible solvable SSB settings, which admit exact analytical solutions for symmetric, antisymmetric, and asymmetric states.The key component of solvable models is the nonlinear term in the NLSE with coordinate x, which is concentrated at x = 0, being represented by the δ-function: This model is formulated in terms of optics, with the evolution along propagation distance z under the action of the real nonlinearity coefficient σ, scaled to be σ = +1 or −1, which corresponds, respectively, to the SF or SDF sign of the nonlinearity.In that case, the δ-function term represents a narrow layer of an optical material with strong cubic susceptibility (e.g., AlGaAs, whose susceptibility exceeds that of silica by a factor ≃ 700 [37]) embedded in the linear planar waveguide, provided that the width of the layer is small in comparison with that of self-trapped light beams propagating in the waveguide.This setting can be readily implemented in the experiment, as the typical width of spatial solitons is measured in tens of microns [36].In that case, the linear trapping potential, −εδ(x), present in Eq.
(1), is relevant too, as the linear refractive index of materials such as AlGaAs is much higher than the background value in the host material (silica).As concerns the sign of the nonlinearity, the consideration of the SDF layer is also interesting, as semiconductor materials may demonstrate negative nonlinear susceptibility.The same equation (1), with z replaced by time t, applies to BEC, with the δ-function potential induced by the FR-inducing laser beam tightly focused at x = 0.The same optical beam also induces the linear potential represented by coefficient ε.In a similar context, Eq. ( 1) with ε = 0 was first introduced, as a model of a nonlinear bosonic junction, in Ref. [38].Further, a model of the matter-wave soliton interferometer with a nonlinear soliton splitter corresponds to ε < 0 and σ = −1 in Eq. (1) [39].
Equation (1) gives rise to the exact solution for a family of self-trapped states (solitons) pinned to the deltafunctional potential: where k is an arbitrary propagation constant, and the shape function is The self-trapped modes are characterized by their integral power, which is a dynamical invariant of Eq. ( 1).The well-known Vakhitov-Kolokolov (VK) criterion, dP/dk > 0 [40][41][42][43], immediately implies that the family of solutions (3) in the case of the SF nonlinearity, σ = +1, and ε > 0 is stable in its entire existence region, k > ε 2 /2 (and completely unstable if the linear potential is repulsive, with ε < 0).For localized states supported by the SDF nonlinearity, with σ = −1, the VK stability criterion is replaced by the anti-VK one [44], dP/dk < 0. Accordingly, in this case the localized states (3) are also stable in their entire existence region, which is 0 < k < ε 2 /2.The definition of the power given by Eq. ( 4) demonstrates that the bound states pinned to δ-function potential with the SF sign of the nonlinearity exist in the interval of 0 < P 0 < 1, while the competition of the linear attractive potential and SDF nonlinear term gives rise to the bound states in the entire range of 0 < P 0 < ∞.
An exceptional case is the one corresponding to σ = +1 (SF) and ε = 0 (no linear potential), for which Eq. ( 4) demonstrates the degeneracy of the localized states, whose power takes the single value, P 0 = 1, which does not depend on k.This property implies that the corresponding family represents a specific example of Townes solitons (a commonly known family of Townes solitons is one produced by localized solutions of two-dimensional NLSE with the spatially uniform cubic SF nonlinearity [45]).Because Townes solitons, with their single value of the power, have dP/dk = 0, the VK criterion predicts that they correspond to a border between the stability and instability.It is known that, in fact, the Townes solitons are subject to the subexponentially commencing instability, which eventually leads to the onset of the critical collapse (emergence of a local singularity after a finite propagation distance) [42,43].
It is also worthy to mention the value of the Hamiltonian of the pinned state (3), [the Hamiltonian is another dynamical invariant of Eq. ( 1)].Note that the existence condition for solution (3), σ √ 2k − ε > 0, implies H 0 < 0 for ε > 0, hence the localized solution represents a true bound state with the negative energy.
The possibility to produce exact analytical solutions for localized states pinned to the δ-function nonlinear potential suggests a possibility to design a solvable DWP model based on a set of two δ-functions, separated by distance which may be set equal to 1 by means of rescaling: The equation for stationary states is produced by the substitution of expression (2) in Eq. ( 6): The Hamiltonian corresponding to Eq. ( 6) is cf. Eq. ( 5).The physical implementation of the model in optics and BEC is straightforward: in the former case, one can embed two parallel nonlinear layers in the linear waveguide, while in the former case the necessary configuration may be created by two tightly focused FR-inducing laser beams.
A particular case of Eq. ( 6) with ε = 0 was introduced, in the context of BEC, in Ref. [46].Exact solutions for symmetric, antisymmetric, and, which is most interesting, asymmetric stationary wave functions were produced in that work, demonstrating a very peculiar feature, namely, an SSB bifurcation of the extreme subcritical type, in which backward-going branches of unstable states never turn forward and, accordingly, never become stable.In other words, it is a unique example of the symmetry-breaking phase transition of the first kind which does not produce any stable phase past the transition point, and gives rise to a fully unstable overcooled phase, represented by the completely unstable asymmetric states.
Recently, another example of such an anomalous phase transition was found in Ref. [47] in the study of dual-core couplers with the SF nonlinearity and fractional diffraction, represented by operator −∂ 2 /∂x 2 α/2 , with Lévy index α [48], acting in each core.In that case, the extreme subcritical SSB bifurcation takes place at α = 1, which is the border between the normal symmetry-breaking phase transition of the first kind at 1 < α < 2, and full instability of the system, driven by the supercritical collapse, at α < 1.However, the fractional-coupler model cannot be solved analytically, on the contrary to Eq. (7).

Objectives of the work
Our purpose is to produce an analytical solution of the full model, with the combined linear-nonlinear δ-functional DWP in Eq. ( 6).The linear terms are represented by ε > 0 (the attractive potential), while both SF and SDF signs of the nonlinearity, σ = ±1, will be addressed.For σ = +1, the solution explicitly demonstrates gradual switch from the extreme subcritical bifurcation to the supercritical one via a regular subcritical bifurcation, in which the backward-going (lower) branches of unstable asymmetric states reverse into stable upper branches at turning points.For σ = −1 the results are more straightforward, corroborating the stability of the symmetric GS and the occurrence of the supercritical antisymmetry-breaking transition in the first excited state.
While the asset of the model is its analytical solvability, some results are produced in a numerical form, using Eqs.( 6) and ( 7) with a regularized δ-function, defined by a small width w (in most cases, w = 0.01 was used, which is 1/100 of the distance between the two δ-funtions).In this connection, it is relevant to mention that the realization of the present model as the optical waveguide implies that a characteristic value of the separation between the two narrow attractive layers may be ∼ 50 µm, hence w = 0.01 corresponds to the layer's thickness ∼ 0.5 µm.In view of the above-mentioned possibility to use a material with the nonlinear susceptibility exceeding that in the bulk waveguide by a factor ∼ 700, this thickness will be sufficient to provide the requires nonlinearity.In the case of the realization in BEC, a relevant size of the separation may be ∼ 10 µm.Then, the nearly-delta-functional potential may be induced by a laser beam focused on a spot of size ∼ 0.5 µm, which will correspond to w ≃ 0.05, in the scaled units.The comparison with the numerical solutions is relevant to check how well the solvable model represents a realistic one, with the finite width w of the potential wells, and also to test predictions for stability of the symmetric, antisymmetric, and asymmetric solitons pinned to the δ-functional DWP.The analytical and numerical results are summarized in the next section, and are discussed in the concluding one.

Exact analytical solutions for the symmetric, asymmetric, and antisymmetric states
The self-focusing nonlinearity The fact that Eq. ( 7) is linear at x = ±1/2 makes it possible to construct obvious solutions in these areas, as exp − √ 2k |x + 1/2| and exp − √ 2k |x − 1/2| at x < −1/2 and x > +1/2, respectively, and a combination of these terms at |x| < 1/2.At points x = ±1/2, the solutions are matched by the continuity condition for U (x) and the jump condition for the derivative dU/dx, The generic solution satisfying these conditions can be looked for as where amplitudes U 1 (k) and U 2 (k) should be found from the substitution of ansatz ( 11)- (13) in Eq. (10).For this stationary solution, the value of Hamiltonian ( 8) is where P is the integral power, defined as per Eq. ( 4).
First, it is easy to find the exact solutions for symmetric states in the SF model (σ = +1), with equal amplitudes where A typical example of a symmetric bound state (soliton), for ε = 2 and k = 2.1, is displayed in Fig. 1(c).This plots is produced by the numerical solution of Eq. ( 7), being virtually indistinguishable from its counterpart given by the analytical solution, as provided by Eqs. ( 11)-( 13) and (15).Because S, as defined by Eq. ( 16), is always positive, the solution given by Eq. ( 15) for σ = +1 and −1 exists for E(ε, k) > 1 and E(ε, k) < 1, respectively.As it follows from Eq. ( 17), this condition implies that, in the case of SF nonlinearity, the symmetric state with given propagation constant k exists if the strength of the linear δ-function potential does not exceed a maximum value, In other words, for given ε the symmetric state exists for k exceeding a value (k min ) symm determined by Eq. ( 18) with < replaced by =, i.e., beneath the red curve in Fig. 2(a).In particular, In the SDF case, the existence area for the symmetric states is opposite, ε > (ε max ) symm .The existence boundary (18) is shown by the red curve in Fig. 2(a).

FIG. 2. (a)
In the model with the SF nonlinearity, σ = +1, the symmetric bound states with amplitudes (15) exist beneath the boundary in the plane of (k, ε) displayed by the red curve, which is produced by Eq. ( 18).The asymmetric states, with the amplitudes given by Eqs. ( 22) and ( 23), exist beneath the green boundary, which is produced by Eq. ( 24).(b) In the model with the SDF nonlinearity, σ = −1, the antisymmetric bound states with amplitudes (15) exist above the brown boundary, which is defined by Eq. ( 18).The states with broken antisymmetry and amplitudes given by Eqs. ( 39) and ( 40) exist above the blue boundary, which is defined by Eq. ( 41).
The bound states (solitons) are characterized by the total power defined as per Eq. ( 4).For the symmetric states in the model with the SF nonlinearity, it is As k varies from the minimum value (k min ) symm [see Eq. ( 19)] towards k → ∞, the power (20) grows from 0 to 2, so that Examples of this dependence for ε = 1 and 2 are displayed in Fig. 3.Note that it satisfies the above-mentioned VK criterion, dP/dk > 0.
An essential fact is that the substitution of ansatz ( 11)-( 13) in Eq. ( 10) produces, as well, an exact solution for asymmetric bound states in the model with the SF nonlinearity, with the following values of amplitudes U 1 and U 2 : FIG. 3. The dependence of the integral power of the symmetric bound states on the propagation constant, in the model with the SF sign of the nonlinearity, as given by Eq. ( 20), for ε = 1 and 2. As shown by Eq. ( 21), with the increase of k the power is slowly approaching the limit value, Psymm(k = ∞) = 2. (or . Typical examples of stable asymmetric states are presented in Fig. 1(b).They are produced as numerical solutions of Eq. ( 7), being indistinguishable from their analytically found counterparts.
Obviously, the solution given by Eqs. ( 22) and (23) bifurcates from the symmetric one (15) (with σ = +1) at E = 2, and exists at E > 2. For a given propagation constant, the asymmetric solution exists if ε does not exceed a respective maximum value, cf. Eq. ( 18).The boundary (24) is shown by the green curve in Fig. 2(a).For fixed ε, the asymmetric solution exists in the region beneath this boundary, and only the symmetric state exists in the stripe between the red and green curves in Fig. 2(a).In particular, (ε max ) symm (k → 0) = 0, i.e., at ε = 0 the symmetric states exist in the entire region of 0 < k < ∞, while it follows from Eq. ( 24) that, in the same limit of ε → 0, the asymmetric state exists, in agreement with Ref. [46], at In accordance with generic properties of the SSB bifurcation [23], the symmetric states are stable solely in the stripe between the red and green curves in Fig. 2(a), being destabilized by the SSB bifurcation beneath the green one.These expectations are corroborated below by direct simulations of the perturbed evolution of the symmetric modes displayed in Fig. 12.
The asymmetry degree of stationary states is defined, in terms of the respective integral power, as Full analytical expressions for the integral power of the asymmetric states, P asy (k), and the respective value of Θ are very cumbersome.Nevertheless, it is easy to find that, while k grows from the minimum value (k min ) asy (ε) at the SSB bifurcation point, which is determined by the left inequality in Eq. ( 24) replaced by the equality [see, in particular, Eq. ( 25) for ε = 0], towards k → ∞, P asy (k) varies from the bifurcation-point value, [with P symm (k) given by Eq. ( 20)] to Actually, Eq. ( 29) is the same value as given above by Eq. ( 4) with σ = +1 and k → ∞.It follows from the above expressions that, as ε increases from zero towards infinity, P bif monotonously decreases from to P bif (ε → ∞) = 0.In particular, P bif (ε) is exponentially small for large ε: Comparison of limit values ( 28) and ( 29) of the integral power for the asymmetric states makes it possible to identify a threshold value ε thr for the switch of the SSB phase transition between the first and second kinds (i.e., the switch between the sub-and supercritical SSB bifurcation): the phase transition may only be of the first kind for P bif > P asy (k → ∞) ≡ 1, and it becomes the second-order transition for P bif < 1.The corresponding equation, P bif = 1, combined with Eq. ( 24), in which ε < (ε max ) asy is replaced, as said above, by ε = (ε max ) asy , amounts to where k thr ≡ (k min ) asy (ε = ε thr ).Numerical solution of Eq. ( 32) produces the single root, k thr ≈ 0.298, with the respective threshold value of ε produced by Eq. ( 24): This result is corroborated by comparison with numerically generated SSB diagrams, in the form of Θ(P ) dependences, which are displayed in Fig. 4. In a detailed form, the numerical data demonstrate that the threshold value belongs to interval 0.07 < ε thr < 0.08, while it is difficult to extract ε thr from the data with higher accuracy.Note that narrow intervals of the variation of P for branches of the asymmetric states in panels (a-f) of Fig. 4 correspond to the analytical results presented here [see, e.g., the limits of the variation of P given by Eqs. ( 29) and (30)]; in panels (g-i), the Θ(P ) curves are partly cut.The range of the variation of P for the branch of the symmetric states, with Θ ≡ 0, is chiefly determined by the limit value (21).For ε = 0, the detailed analysis, reported for this case in Ref. [46], demonstrates, in agreement with Fig. 4(a), that the largest power of the symmetric solitons is (P symm ) max ≈ 2.08, which is attained at k ≈ 1.40.
The asymmetric solitons are completely stable in the area ε thr < ε < (ε max ) asy , as illustrated below by Fig. 13.At ε < ε thr , the asymmetric solutions belonging to the lower branches in Figs.4(a-d), with dΘ/dP < 0, are unstable, while the upper branches, with dΘ/dP > 0, are stable.Actually, the instability intervals for the asymmetric solitons are very narrow.
In addition to the symmetric and asymmetric stationary states, Eqs. ( 6) and ( 6) with the SF sign of the nonlinearity, σ = +1, give rise to antisymmetric ones, with U (−x) = −U (x), see an example in Fig. 1(a).However, as well as in the case of ε = 0 [46], the antisymmetric states are completely unstable because, for the same value of integral power P , they correspond to higher values of Hamiltonian (8) than the symmetric bound states.The instability of the antisymmetric states is illustrated below by Fig. 14.

The self-defocusing nonlinearity
Typical examples of antisymmetric, broken-antisymmetry, and symmetric states produced by Eq. ( 7) with the SDF nonlinearity, i.e., σ = −1 in (7), are displayed in Figs.1(d-f), respectively.In this case, the symmetric state, given by solution (15) with σ = −1, which exists, as said above, at ε > (ε max ) symm [see Eq. ( 18)], is always stable, realizing FIG. 4. The asymmetry parameter (27) for the numerically produced solutions of Eq. ( 7) with the SF nonlinearity (σ = +1) vs. the integral power (26) at different values of strength ε of the linear δ-functional potential, which are indicated in panels.The switch between the symmetry-breaking phase transitions of the first and second kinds, alias sub-and supercritical SSB bifurcations, takes place between ε = 0.07 and 0.08, in agreement with analytical result (33).the model's GS.Accordingly, it is not subject to SSB.More interesting is the first excited state above the GS, i.e., the antisymmetric one, given by the Eqs.( 11)-( 13) (with σ = −1) where S and E is defined, as above, as per Eqs.( 16) and (17).Because S is always positive, this solution exists under condition E < −1.The substitution of Eq. ( 17) demonstrates that this condition amounts to cf. Eq. ( 18).The antisymmetric state exists, at ε > 1, in the area of the (µ, ε) plane above the brown boundary shown in Fig. 2(b).Because Eq. ( 35) yields ε ≥ 1 in the limit of k → 0, there are no antisymmetric states at ε < 1.
The integral power of the antisymmetric state is This dependence is displayed in Fig. 5 for ε = 2.Note that expression (36) with all values of ε satisfies the abovementioned anti-VK criterion, dP/dk < 0, which is necessary for the stability of bound states supported by the SDF nonlinearity [44].
FIG. 5.The dependence of the integral power of the antisymmetric bound states on the propagation constant k, in the case of the SDF sign of the nonlinearity (σ = −1), as given by Eq. ( 36), for ε = 4.At k → 0 the power diverges according to Eq. ( 37).The power vanishes at k ≈ 7.686, which is determined by Eq. ( 35) with ε = 4.
With the variation of k from the largest value, k max , which is determined by Eq. ( 35), towards k → 0, the power (36) monotonously increases from P anti (k max ) = 0 to values diverging as at k → 0. The divergence is explained by the fact that, in the limit of k → 0, there is an antisymmetric delocalized state with divergent power, The bound state with broken antisymmetry is given by Eqs. ( 11)-( 13) with amplitudes This solution exists under condition E < −2.The substitution of expression (17) in the latter condition leads to the following existence area for the solutions with broken antisymmetry: cf. Eq. ( 35).This area is located above the blue boundary in Fig. 2(b).Because Eq. ( 41) yields ε ≥ 3/2 in the limit of k → 0, there are no states with broken antisymmetry at ε < 3/2.In agreement with the existence of the delocalized antisymmetric state (38), there is also a delocalized state with broken antisymmetry, viz., where The mirror image of this solution is also a delocalized state with broken antisymmetry.Note that the delocalized antisymmetric state and the one with the broken antisymmetry exist, according to Eqs. ( 38) and ( 43), at ε > 1 and ε > 3/2, respectively, in accordance with what is said above for the generic solutions of the same types.

NUMERICAL SOLUTIONS
Numerical investigation of Eqs. ( 6) and ( 7) with the δ-functions approximated as per Eq. ( 9) is relevant for checking the analytical results reported above.Because the width of the linear and nonlinear potential wells in real systems is finite, the numerical results are also relevant for the verification of the relevance of the analytical predictions, which are obtained above with the use of the ideal δ-functions.

The self-focusing nonlinearity
Numerically found examples of bound states of the symmetric and antisymmetric types, as well as ones with broken symmetry and antisymmetry, in the cases of the SF and SDF signs of the nonlinearity, are displayed above in Fig. 1.In a systematic way, the evolution of the antisymmetric, asymmetric, and symmetric solitons produced by Eq. ( 7) with σ = +1, following the increase of propagation constant k, is summarized in Fig. 6.
The most essential results in the form of the SSB diagrams for the SF model, which corroborate the basic analytically predicted property of the model, viz., the switch of the character of the symmetry-breaking phase transition from the first to second kind (in other words, the switch from the subcritical SSB bifurcation to the supercritical one) at the threshold point (33), are also demonstrated above in Fig. 4. In addition to that, it is relevant to plot the bifurcation diagrams in the plane of k and asymmetry parameter Θ.These are presented in Fig. 7, for the same set of values of ε as in Fig. 4. The branch of the symmetric states commences at k = (k min ) symm [see Eq. ( 19)], while the value of k(ε) at the SSB bifurcation point is determined by Eq. (24).
Further, the families of symmetric and asymmetric solitons are characterized, as physical states, by the respective dependences P (k) and H(P ), where H is the Hamiltonian defined by Eq. ( 8).These dependences are displayed, respectively, in Figs. 8 and 9.In the former figure, the branches of the symmetric states commence at k = (k min ) symm , see Eq. (19).In panels (a-f) of Fig. 8, P varies between limit values 0 and 2 along the symmetric branches, and between P bif [see Eq. ( 28)] and P = 1 along the the asymmetric ones [in panels (g-i), the variation range of P is truncated; it is also partly cut in Fig. 9(f)].In Fig. 9(i), dependences H(P ) for the symmetric and asymmetric states FIG. 7. The asymmetry parameter (27) for the numerically produced solutions of Eq. ( 7) with σ = +1 vs. propagation constant k at the same values of ε which are presented in Fig. 4.
are indistinguishable.Note also that, in latter case, the value P bif at the SSB point is extremely small, in agreement with Eq. (31).On the other hand, coordinates of the SSB points in Figs.8(a) and 9(a) are correctly predicted by Eqs. ( 25) and (30).
In the range where the asymmetric states exist, they realizes the minimum of H, i.e., the system's GS.A specific feature of the system is that it has no true GS at larger values of P , where only the unstable symmetric states exist, and there are no stationary states whatsoever at P > 2.
The self-defocusing nonlinearity Dependences P (k), H(P ), Θ(k), and θ(P ) for the families of antisymmetric solitons and those with broken antisymmetry, as produced by the numerical solution of Eq. ( 7) with σ = −1, are collected, severally, in panels (a-c), (d-f), (g-i), and (j-l) of Fig. 10, for three different values of the strength of the linear δ-function potential, viz., ε = 2, 3, and 5.These sets of plots are counterparts of those for the model with σ = +1 which are displayed above in Figs. 8, 9, 7, and 4, respectively.In particular, the P (k) curves and SSB points on all the curves plotted in Fig. 10 are correctly predicted by Eqs.(36) and (35), respectively.
An obvious difference from the case of the SF nonlinearity is that the bifurcation of the spontaneous breaking of antisymmetry in the SDF case is always supercritical, as seen in Figs.10(j-l).In other words, the model with the SDF nonlinearity always gives rise to the antisymmetry-breaking phase transition of the second kind.It is also worthy to note that the soliton branches with both unbroken and broken antisymmetry always satisfy the above-mentioned anti-VK criterion, dP/dk < 0, which is the necessary (but not sufficient) condition for their stability.
Finally, the evolution of the antisymmetric, broken-antisymmetry, and symmetric bound states produced by Eq. ( 7) with σ = −1 and ε = 2, following the increase of propagation constant k, is summarized in Fig. 11.Note that, in agreement with the analytical solutions, the evolution is opposite to that in the model with the SF nonlinearity (σ = +1), which is displayed above in Fig. 6.Namely, the amplitude and integral power of the solitons increase/decrease with the growth of k in the SF/SDF system.

The evolution of unstable bound states
It is relevant to test the expected (in)stability of symmetric and antisymmetric bound states, as well as ones with broken symmetry and antisymmetry, in direct simulations of Eq. ( 6) with the ideal δ-function replaced by its regularized version (9), for both the SF and SDF signs of the nonlinearity, i.e., σ = +1 and −1.
First, Fig. 12 collects typical examples which demonstrate the perturbed evolution of stable [panels (a,d,f,h,i)] and unstable [panels (b,c,e,g)] symmetric bound states in the model with the SF nonlinearity.These results are compatible with the prediction of the stability area for the symmetric states, in the form of the stripe between the lower and upper boundaries in Fig. 2(a).It is observed that, naturally, all the unstable symmetric states demonstrate manifestations of the SSB instability, leading to spontaneous formation of asymmetric states.In some cases, such as the one displayed in panel 2(f), the unstable symmetric state, which is located close to the instability boundary, features conspicuous persistent oscillations, while in other cases the stronger instability creates nearly stationary modes with strong asymmetry.
Another expected result corroborated by the direct simulations of the perturbed evolution is that nearly all the asymmetric solitons are stable in the case of the SF nonlinearity, as shown in Fig. 13 for strongly asymmetric solutions.Unstable are asymmetric solitons belonging to the backward-going (lower) branch in Figs.4(a-d).In fact, they exist only in a very narrow parameter region, and the development of the instability pulls them towards a stable counterpart existing at the same value of P (not shown here in detail, as this is a known feature of the subcritical SSB bifurcation).
In addition to the above results, direct simulations displayed in Fig. 14  antisymmetric bound states in the case of the SF nonlinearity.In the cases shown in panels (e) and (f) of the figure, the instability is barely visible as the interaction between two power peaks of the antisymmetric modes is very weak.
Lastly, characteristic examples of the perturbed evolution of the bound states of the symmetric, antisymmetric, and broken-antisymmetry types in the model with the SDF nonlinearity are collected in Fig. 15.In particular, in agreement with the above predictions all symmetric states are stable in this case, as shown in Figs.15(a-c).Further, panels (d,e) and (f) of Fig. 15 present, respectively, examples of moderately and weakly unstable antisymmetric states, in agreement with the boundaries plotted in Fig. 2(b).It is seen that the instability leads to spontaneous replacement of the corresponding states by oscillatory ones with broken antisymmetry.Finally, panels (g-i) demonstrate stability of the stationary states with weakly or strongly broken antisymmetry, also in agreement with the boundary plotted in Fig. 2(b).

DISCUSSION
Starting from the two-dimensional Ising lattice [49,50], exactly solvable models serve as benchmarks for studies of phase transitions in diverse physical settings [51]- [55].Transitions from a paramagnetic phase to a ferromagnetic one in spin systems, and similar transitions in many other media are intrinsically related to spontaneous breaking of the symmetry of the underlying setting.It is well known that phase transitions are classified as ones of the first and second kinds.Hysteresis and bistability, in the form of the coexistence of the GS (ground state) with a metastable overcooled or overheated phase, are possible in the former case.
Similar phenomenology is exhibited by nonlinear dynamical systems, in the form of bifurcations, i.e., transitions between different stable states of the system, caused by variation of the system's control parameter(s) [23].Counterparts of the phase transitions of the first and second kinds are identified as bifurcations of the subcritical and supercritical types in dynamical systems.The subcritical bifurcation creates stable states prior to the destabilization of the symmetric one, thus the bifurcation of this type admits the bistability and hysteresis, like phase transitions of the first kind.
In most cases, phase transitions in statistical physics, as well as bifurcations in dynamical systems, are studied between spatially uniform states.On the other hand, transitions between spatially localized (self-trapped) modes, such as solitons, are possible too.The analysis of the latter topic may benefit from the consideration of models admitting exact solutions for symmetry-breaking transitions in self-trapped states.However, finding solvable models is a challenging task, because basic integrable models that give rise to solitons, such as the one-dimensional NLSE, do not admit intrinsic transitions in the solitons.
The objective of the present work is to introduce a solvable nonlinear model with the DWP (double-well potential) which makes it possible to produce exact solutions for localized states with full and broken symmetries, that are linked by symmetry-breaking transitions of both first and second kinds.In other words, the states with unbroken and broken symmetries may be linked by bifurcations of the sub-and supercritical types.The integrability of the present model is possible due to the fact that the nonlinearity is represented by the symmetric set of two δ-functions.A prototype of this model was introduced previously in Ref. [46], but it had produced a very limited result, viz., the SSB (spontaneous-symmetry-breaking) bifurcation of the extreme subcritical form.That bifurcation gave rise to completely unstable asymmetric states, represented by backward-going solution branches which never turned forward.In the present work, we have introduced the solvable DWP model including both nonlinear and linear potentials, which are based on the symmetric pair of δ-functions.The respective nonlinear potential is considered with both the SF (self-focusing) and SDF (self-defocusing) signs.
The analytical solutions, confirmed by their numerically found counterparts (which were produced replacing the   functions.Starting from the above-mentioned extreme subcritical bifurcation at ε = 0, the switch is found analytically to occur at the point given by Eqs. ( 32) and (33), which is corroborated by the numerical findings.To the best of our knowledge, no previously studied model made it possible to predict the change of a symmetry-breaking phase transition between the first and second kinds (or the change of the sub/supercritical character of the SSB bifurcation) in an analytical form.The analytical solution is also reported here for the model with the SDF nonlinearity, where the situation is simpler: the GS is always represented by the completely stable symmetric localized state, while the antisymmetry-breaking phase transition of the second kind (i.e., the supercritical bifurcation) destabilizes the lowest excited state (a spatially odd stationary one) at the critical point given by Eq. ( 41).These analytical results for the SDF model are confirmed by the numerical solution too.
The solvable models elaborated in the present work suggest possibilities for analytical studies of SSB phase transitions and bifurcations in more complex settings.In particular, it may be interesting to address a two-component system with the combined linear-nonlinear DWP potential.A degenerate form of the two-component system, with the nonlinear-only SF potential, based on the symmetric pair of δ-functions, was introduced in Ref. [56].In that model, the SF nonlinearity includes self-interaction in each component and cross-interaction between the components.Note that the two-component SF model, unlike the single-component one, admits an antisymmetry-breaking phase transition in spatially odd localized states, and it also opens the way to the consideration of the SSB transition in a state which combines spatially symmetric and antisymmetric components.
Another new possibility is offered by a model with three equidistant δ-functions set on a circle, unlike the infinite onedimensional domain considered in the present work (the circle with the purely nonlinear SDF potential, represented by a symmetric pair of δ-functions set at diametrically opposite points, was addressed in Ref. [57], in which case it did not give rise to SSB transitions, the respective GS being always symmetric).Various setups with a triangle of potential wells embedded in a nonlinear medium were studied for BEC [58,59] and multicore optical fibers [60]- [64].In the circular setting with three δ-function wells, one can construct exact solutions carrying vorticity and address feasible SSB transitions in them.
On the other hand, as a step towards the consideration of two-dimensional models, where full solvability is not plausible, one can introduce a set of two parallel linearly-coupled one-dimensional lines, each bearing a DWP represented by the symmetric pair of the δ-functions.In all these extensions, analytical solutions will take an essentially more cumbersome form than the one addressed in the present work, but the analysis may still be possible.

FIG. 1 .
FIG. 1. (a-c) Typical examples of antisymmetric, asymmetric, and symmeric bound states (solitons) produced by the numerical solution of Eq. (7) with the δ-functions approximated by expression (9), the SF nonlinearity (σ = +1), and parameters ε = 2, k = 2.1.In panel (b) two asymmetric states are plotted, which are mirror images of each other.The antisymmetric and symmetric states are unstable, while the asymmetric one is stable.(d-f) Typical examples of antisymmetric, brokenantisymmetry, and symmetric bound states for the SDF nonlinearity (σ = −1) and ε = 2, k = 1.In panel (e), two states with broken antisymmetry are mirror images of each other.The antisymmetric state is unstable, while the ones with broken antisymmetry and unbroken symmetry are stable.

FIG. 8 .
FIG. 8.The integral power of the symmetric and asymmetric bound states in the case of the SF nonlinearity (σ = +1) vs. the propagation constant for the same values of ε as in Figs. 4 and 7.
FIG. 9.The Hamiltonian of the symmetric and asymmetric bound states, calculated as per Eq.(8) with σ = +1, vs. the integral power, P , for the same values of ε as in Figs. 4, 7, and 8.