Space-time quasicrystals in Bose-Einstein condensates

An autoresonant approach for exciting space-time quasicrystals in Bose-Einstein condensates is proposed by employing two-component chirped frequency parametric driving or modulation of the interaction strength within Gross-Pitaevskii equation. A weakly nonlinear theory of the process is developed using Whitham's averaged variational principle yielding reduction to a two-degrees-of-freedom dynamical system in action-angle variables. Additionally, the theory also delineates permissible driving parameters and establishes thresholds on the driving amplitudes required for autoresonant excitation.


I. INTRODUCTION
Quasicrystals in materials are structures ordered in space, but not exactly periodic [1][2][3].By analogy, time or space-time quasicrystals are ordered systems in spacetime, but not exactly periodic in space and/or time.By now, these systems have been studied in periodically driven magnon condensates [4,5] and ultracold atoms [6,7], where temporal symmetry was destroyed due to subharmonic response.In this study, we explore a different path to space-time quasicrystals.It is well known that a number of so-called integrable nonlinear partial differential equations (PDE's) have multiphase solutions [8] of the form f (θ 1 , θ 2 , ..., θ N ; λ 1 , λ 2 , ..., λ N ), where N phase variables θ i = k i x − ω i t have wave numbers k i which are multiples of some k 0 (the solution is spatially periodic), λ i are constants, while frequencies ω i are functions of k i and λ i .By choosing some set of λ i , one can make some or all of these frequencies incommensurate, forming an ideal space-time quasicrystalline structure, which is periodic in space and aperiodic in time, still having a complex long range time ordering.Examples of such PDEs are the Korteweg-de-Vries (KdV), sine-Gordon (SG), and nonlinear Schrodinger (NLS) equations, which find multiple application in physics [8].However, why these ideal space-time quasicrystals are not yet realized in experiments?The answer lies in complexity, as their analysis typically requires advanced mathematical methods, such as the inverse scattering transform (IST) [9], while experimental realization depends on forming complicated space-time dependent initial conditions, which is an unrealistic task.
In this work, we suggest a different approach to realizing space-time quasicrystals based on autoresonance [10].Autoresonance is an important nonlinear phenomenon, where a system phase-locks to a chirped frequency driving perturbation and remains phase-locked continuously for an extended period of time, despite variations of the driving frequency.As the driving frequency varies in time, so does the frequency of the excited solution, lead-ing to formation of a stable highly nonlinear state.We will focus on autoresonant formation and control of twophase excitations of the Gross-Pitaevskii (GP) equation describing Bose-Einstein condensates (BECs).The excitation proceeds from trivial initial conditions and uses a combination of two independent small amplitude wavelike driving perturbations.In contrast to existing applications using large amplitude pulsed drives or optical crystalline structures for containing BEC's, we can remove the driving perturbation after some time, remaining with a free, slightly perturbed, but stable nonlinear two-phase GP solution.These excited quasicrystalline structures are controlled by two independently chirped driving frequencies, thus exploring a continuous range of parameters λ i , i.e., a continuous set of space-time quasicrystals.The autoresonance approach has been used previously in excitation of multi-phase waves in different applications with the theory based on the IST method [11][12][13].Here we will apply a simpler analysis of the process of capture into a double autoresonance in the system using two chirped frequency drives, similar to recent studies on the formation of two-phase waves in plasmas [14,15].Our presentation will be as follows.In the next section, we illustrate the formation of space-time quasicrystals in BECs through numerical simulations.Section III presents the quasi-linear theory of formation of GP quasicrystals using Whitham's averaged variational approach [16].Section IV addresses the problem of the allowed parameter space for autoresonant excitations and with the associated threshold phenomenon on the driving amplitudes.Finally, Sec.V presents our conclusions.

II. QUASICRYSTALS IN A BEC VIA SIMULATIONS
The basic model for studying nonlinear dynamics of BECs is GP equation [17],written in dimensionless form Here, time is measured in units of the inverse transverse trapping frequency ω −1 ⊥ , while space and density are measure in units of l ⊥ = [ /(2mω ⊥ )] 1/2 and mω ⊥ /2π |a 0 |, respectively, where m represents the atomic mass.In Eq. ( 1), g = 2a(x, t)/|a 0 | is the normalized, space-time modulated interaction strength, where a represents the s-wave scattering length of interacting particles in the BEC.For condensates with repulsive interactions of particles a < 0 and for attractive interactions a > 0. U(x, t) denotes the longitudinal potential.We will assume that our system is perturbed by a combination of independent, small amplitude waves (2) where ψ i (t) = ω di (t)dt and ω di (t) = ω 0i −α i t are slowly chirped driving frequencies.We consider two driving options.The first is a parametric-type driving U(x, t) = f , |f | ≪ 1, while the interaction strength is not perturbed, i.e., g = 2σ and σ = ±1.The second driving scenario is a modulation of the interaction strength by external magnetic field,for example, near Feshbach resonance [18].In this case, we assume U = 0, and again |f | ≪ 1.For both driving options we can rewrite Eq. ( 1) as a weakly perturbed NLS equation where F is either −f or −2σ|ϕ| 2 f .In computer simulations below, we will use the parametric driving, assume periodic boundary conditions ϕ(x, t) = ϕ(x + l, t), thus k 1,2 are multiples of k 0 = 2π/l.Nevertheless, in Secs.III and IV, we will discuss both driving options.It should be mentioned that periodic boundary conditions are usually assumed in numerical simulations of an infinite domain for a spatially periodic driving as well as in 1D modeling along the torus-like BEC (ring-trap geometry) [19,20].
A special case when the driving (3) is a standing wave, i.e. ε 1 = ε 2 , ψ 1 = ψ 2 , and k 1 = −k 2 was studied recently [21], so the present investigation is a generalization to two independent driving components.
In the periodic case, the unperturbed ground state of a BEC is a spatially homogeneous solution of Eq. ( 4) with constant amplitude |ϕ| = U 0 .The frequency of a perturbation of the homogeneous state is [22] Condensates with repulsive interaction of particles, when σ = −1, are stable.In this case, frequency (6) is known as the Bogolubov frequency.Dark solitons are typical structures in these condensates.In the opposite case (σ = 1), bright solitons exists.In this case ω 0 can be imaginary, leading to modulational instability.This instability is well-known in plasma physics and nonlinear optics [23,24].If a condensate has a length l, then the wave number of the main mode is k 0 = 2π/l and the stability condition restricts the density of the condensate to U 2 0 < π 2 /l 2 .If the condensate has a cigar-like shape with the transverse dimension l ⊥ , then, in physical variables, the stability condition can be written as a restriction on the number of particles, n < (l ⊥ /l)(l ⊥ /|a 0 |).
We proceed to numerical simulations of autoresonant formation of a space-time GP quasicrystal by focusing on the case of σ = −1 and starting in the ground state (5) with U 0 = 1.The driving parameters are ε 1,2 = 0.01, k 1,2 = 1, −3, chirp rates α 1,2 = 0.0012, 0.0024, and ω 0i are given by Eq. ( 6) for each k i .The simulation begins at t 0 = −300 and both components of the drive pass the corresponding Bogolubov resonances at t = 0. Furthermore, we switch off the drives at t = t s = 250.The corresponding time dependence of the maximum (over x) |ϕ| 2 is shown in Fig. 1 by the upper blue line.One can observe that the excitation amplitude increases continuously until the drive is turned off, and the maximal amplitude remains nearly stationary afterwards.We also show the excited quasicrystalline structure in the time interval 380 < t < 400 (after the driving is switched off) using a colormap in the upper panel of Fig. 2. Furthermore, the ratio of the two driving frequencies at t = t s in this example is approximately 1 : 5.The second panel from the top in Fig. 2 shows |ϕ| 2 versus time in the same example at x = 0.One can see the short and long driving periods in the panel illustrating 1 : 5 quasi-periodicity and the two-phase locking with the drives.
To further illustrate the characteristics of autoresonant excitation in the system, we show, by the lower two lines in Fig. 1, the cases when only one of the two drives is present in the same example and illustrate the colormaps of the associated excitations in the lowest two panels in FIG. 2. Two-and one-phase autoresonant excitations in simulations.The upper panel shows the color map of a two-phase solution |ϕ| 2 in space-time for driving parameters ε1,2 = 0.01, k1,2 = 1, −3, chirp rates α1,2 = 0.0012, 0.0024.The second panel from the top shows |ϕ| 2 versus time at x = 0, illustrating 1:5 quasi-periodicity of the solution.The lowest two panels show color maps of autoresonant single-phase excitations |ϕ| 2 for the same parameters and initial conditions as in the upper panel, but when only one of the diving components is applied.Fig. 2. A single-phase parametric autoresonant excitation in this system was analyzed in Ref. [25].In this case, one forms a growing amplitude nonlinear wave traveling with the phase velocity ω di /k i of the corresponding drive.The directions of these propagation velocities are clearly seen in Fig. 2. We observe the same two characteristic directions in the upper panel in Fig. 2 corresponding to the two-phase autoresonant excitation, illustrating again the continuing phase locking with the two driving components.
One of the most important issues associated with the autoresonance is the threshold on the driving amplitudes for the continuing phase locking in the system.Figures 1 and 2 show that the two-phase autoresonant excitation (where both drives are present) is very different from that with a single drive.This means that the driving amplitudes must be sufficiently large to obtain a twophase quasicrystalline structure, leading to the problem of thresholds for the transition to autoresonance.We illustrate this sudden transition in Fig. 3, showing the maximal |ϕ| 2 versus time in three cases with the same parameters as in Fig. 1, but ε 2 = 0.006, 0.005, and 0.004, while keeping ǫ 1 = 0.01.One observes a sharp transition when passing from ǫ 2 = 0.005 to 0.004.The color maps of |ϕ| 2 in space-time in these three cases are shown in Fig. 4. One can see that the quasicrystalilne structure in the upper two panels in the Figure is similar to that shown in the upper panel in Fig. 2, but the struc- ture changes significantly as one passes to ǫ 2 = 0.004.The lowest panel in Fig. 4 has smaller amplitude and is closer to the lowest panel in Fig. 2, corresponding to the single phase excitation.We interpret this transition as the loss of the phase locking when ǫ 2 is below some threshold value between ǫ 2 = 0.004 and 0.005.
As in all single phase autoresonant interactions, the double autoresonant phase locking in the driven GP system starts in the initial excitation stage, as the two drives simultaneously pass through the linear (Bogolubov) frequencies in the problem.We also find that the autoresonant phase locking is a weakly nonlinear phenomenon.Therefore, the next section is devoted to the quasi-linear theory of two-phase GP autoresonance.

III. WEAKLY NONLINEAR THEORY
In developing the weakly nonlinear theory of autoresonant two-phase solutions of the GP equation, we will focus primarily on the parametric-type driving scenario.The case of driving by modulation of the interaction strength is obtained similarly and will be briefly discussed at the end of this Section.Therefore, we proceed from Eq. ( 4) for the ponderomotive case where θ di = k i x − ψ di (t) and seek solution of form ϕ = U exp(iV ) governed by the following set of real equations 9) The Lagrangian density for this problem is This Lagrangian representation suggests using Whitham's averaged variational approach [16] to analyze our problem.The first step in this direction is to assume constant frequency drives, ψ di = ω di t, and seek phase-locked solutions of the linearized problem of form [Note that our unperturbed solution is ϕ 0 = U 0 exp(2iσU 2 0 t)].Then, by linearization, and neglecting products of linear amplitudes U i and ε i , Eqs (8) and (9) become yielding solutions where the linear resonance frequencies Now, we proceed to chirped-driven problem, where ψ i = ω di (t)dt and extend Eqs.(11) and (12) to next nonlinear order Here U i and V i are small (viewed as first-order perturbations), while all other amplitudes are assumed to be of second-order in U i and V i .In these solutions θ i = k i x−ψ i and ψ i = ω i (t)dt is a new fast independent variable.At this stage, we do not assume phase-locking in the system, but view the difference Φ i (t) = ψ i − ψ di as slow function of time.Similarly, all the amplitudes are also assumed to be slowly varying functions of time.The reason for choosing the second-order ansatz of this form is consistent with the form of the Lagrangian density containing either different powers of U or products of derivatives of V and powers of U .The auxiliary phase ξ = γ(t)dt in Eq. ( 19) is necessary because V is the potential (it enters the Lagrangian density via derivatives only [16]).
The next step is to replace θ di = θ i + Φ i (t) in the driving part of the Lagrangian density (10), substitute the above ansatz into the Lagrangian density, and average it over θ i ∈ [0, 2π].This averaging is done via the Mathematica package in the Appendix.The resulting averaged Lagrangian density Λ ) is a function of all 13 slow first and second-order amplitudes and Φ 1,2 and ξ.The Lagrange's equations for all these 16 variables, form a system describing slow autoresonant evolution in the problem.Reducing this problem to a smaller set of evolution equations involves tedious algebra, which, nevertheless, can be performed using the Mathematica package.The details of this reduction are given in the Appendix, and here we present the final closed system of 4 equations for U 1 , U 2 , Φ 1 , and Φ 2 (see Eqs.((71),(72),(74)) in the Appendix): (ω All remaining dependent variables in the problem, i.e., V 1,2 , u 0 , u 11 , u 22 , u 12p , u 12m , v 11 , v 22 , v 12p , v 12m , Φ 1,2 , ξ are related to U 1 , U 2 (see Eqs. (59-68) in the Appendix).Before proceeding to the analysis of this system, we rewrite Eqs.(( 21),( 22)) explicitly as differential equations for Φ 1,2 .We approximate , which allows to write these equations as Equations ( 20), (23), and (24) comprise a complete set of differential equations for studying the dynamics in our problem.By solving this system, and calculating all second-order objects as described above, we obtain a full quasilinear two-phase solution (18,19) of the chirpeddriven GP equation.As an illustration, Fig. 5 shows the dynamics of U 1,2 and Φ 1,2 for two cases with the same parameters as in the two lower panels in Fig. 4. The figure illustrates the loss of phase-locking with one of the components of the drive as one passes from ǫ 2 = 0.005 to 0.004.The upper two panels in Fig. 5 show double phase locking of Φ 1,2 at π [mod(2π)] and a continuous autoresonant growth of U 1,2 , while in the lower two panels only U 1 continues to grow, while U 2 saturates as Φ 2 escapes.Thus, one requires both driving amplitudes to be above some minimal values for a persisting double autoresonance in the system.We discuss this threshold effect in the next section.
We conclude this section by discussing the case of driving via the modulation of the interaction strength (see the comments at the beginning of Sec.II).It is shown at the end of the Appendix that the reduced system of equations describing two-phase autoresonant formation of space time GP quasicrystals in this case is the same as for the ponderomotive drive (see Eqs. ( 20),( 23),( 24)), but U 0 in the driving terms must be replaced by 2σU 3 0 .We present an example of such a case in Fig. 6, where the parameters are the same as in the lower two panels FIG. 5. Threshold phenomenon.The upper two panels illustrate double phase-locking in the system for ǫ2 = 0.005, while in the lower two panels, where ǫ2 = 0.004, the phase locking with one of the driving components is lost.All other driving parameters in the two cases are the same as in Fig. 4, i.e., ǫ1 = 0.01, α1,2 = 0.0012, 0.0024.FIG. 6.Double phase-locking in the system driven by modulation of the interaction strength.The parameters are the same as in the lower two panels in Fig. 5, but U0 in the driving term is replaced by 2σU 3 0 .
in Fig. 5, but the driving is via the interaction strength in the GP equation, as discussed above.One can see two differences between between Fig. 5 and 6.One is that in the case of the interaction strength drive the double phase-locking is restored.This is because effectively the driving amplitude increased by a factor of two, and passed the threshold.The second difference is that the phase mismatches Φ 1,2 in Fig. 6 are locked near 0. The reason is that in this case, there is a new factor of σ in the driving perturbation, which for σ = −1 in this example changes the phase-locking location, as will be discussed in the next Section.

IV. AUTORESONANCE CONDITIONS AND THRESHOLD PHENOMENON
In this section, we discuss conditions for the autoresonant evolution of space-time quasicrystals in BECs.We proceed by defining new (action) variables Then the system of Eqs.(20,23,24) describing two-phase BECs can be rewritten as where The action-angle Hamiltonian of this system is Note that as mentioned at the end of Sec.III, in the case of two-phase driving via the modulation of the interaction strength, the system of equations ( 26)-( 29) remains the same, but in the expression for η 1,2 , U 0 must be replaced by 2σU 3 0 .Finally, a similar system of equations was derived recently in application to Langmuir and ion acoustic waves in plasmas [14,15].
The autoresonant evolution corresponds to double phase-locking in the system as phase mismatches Φ 1,2 remain bounded continuously subject to small initial conditions on I 1,2 at large negative time t 0 .In the initial stage, the pairs of variables (I 1 , Φ 1 ) and (I 2 , Φ 2 ) decouple and are described by The phase-locking dΦi dt ≈ 0 in each of these decoupled systems is guaranteed at large negative times (see Ref. [26] for a detailed analysis) and yields and phase locking at either Φ 1,2 ≈ 0 or π if α 1,2 is negative or positive, respectively, and in both cases Next, assuming that the autoresonant phase-locking ( dΦi dt ≈ 0) continues as the system reaches large positive times, in double autoresonance, the actions I 1,2 are given by the solution of Since coefficients a, b, c are all positive, for having positive solutions for I 1,2 chirp rates α 1,2 must have the sign of −σ, i.e., can be written as α 1,2 = −σ|α 1,2 |.This yields autoresonant solutions varying linearly in time The examples of space-time quasicrystals for these two sets of parameters are shown in Fig. 9 where Then D = ac − b 2 must be positive.This is obviously the case for σ = −1, so we can always find some ratio |α 2 /α 1 | in this case for having positive I 1,2 , linearly increasing in time.The case σ = 1 is more complex.Since we still need D > 0 for having increasing I 1,2 at large positive times, we must satisfy where The last inequality yields the condition The region S in the (Y, X)-plane satisfying this inequality (the allowed region) is shown in Fig. 7.Note that in this region X, Y < 3/4, which guarantees the positivity of ω 2 r1,2 .However, U 2 0 can not be too large to satisfy X, Y < 3/4 and can be chosen as follows.Let k 2 2 > k 2 1 .Then, X < Y and we can choose some value Y 0 < 3/4 yielding . This guarantees that the point (X 0 , Y 0 ) is inside the allowed region if Y 0 < 1/2.However, if 1/2 < Y 0 < 3/4, for having (X 0 , Y 0 ) ∈ S we have a restriction The inequalities (41) and ( 42) above are based on the analysis at large positive times and comprise only necesary conditions for synchronized (autoresonant) evolution.We have already discussed the phase-locking at large negative times.However, for synchronized passage through the vicinity of t = 0, i.e., for having bounded Φ 1,2 at all times, in addition to the above, it requires η 1,2 be large enough (for both σ = +1 and −1).In dealing with this issue, we choose some value of' |α 1 |, η 1 and r = |α 2 /α 1 | (r must satisfy (38) as described above), which defines |α 2 |.We also fix ratio q = η 2 /η 1 , which defines η 2 , and we are left with the problem of finding the critical value of η 1th for autoresonant transition through the vicinity of t = 0 for this |α 1 |.Now we return to our original system ( 26)-( 29) and rewrite it as where . Therefore, for a given r and q, we are left with a single additional parameter µ, and there may exists some minimal value of µ th in the problem which still guarantees a continuous phase-locking in the system as it passes from large negative times through t = 0, to large positive times.This value can be found numerically, yielding the minimal (threshold) driving amplitude ε 1th for autoresonance in the system: the case of driving via modulation strength, the threshold formula remains the same, but U 0 is replaced by 2U 3 0 .We illustrate the characteristic 3/4 power scaling of ε 1th with |α 1 | in Fig. 8, for the parametric driving case, σ = −1, and two sets of parameters: Set A: k 1 = 2, k 2 = −3, U 0 = 0.775, r = 1.08, q = 0.7 (red line) and Set B: k 1 = 1, k 2 = −2, U 0 = 0.447, r = 0.805, q = 1.3 (blue line).The values of k 1,2 and U 0 in these two examples correspond to the two points in the allowed autoresonance region in Fig. 7.The circles and squares in Fig. 8 show the results of numerical simulations.Finally, Fig. 10 shows colormaps of |ϕ| 2 in space-time autoresonant space-time quasicrystals for the same two sets of parameters (the upper and lower panels correspond to sets A and B, respectively).In both cases, α 1 = 0.001, and the driving amplitudes are 10% above the thresholds in Fig. 8.

V. CONCLUSIONS
We have shown that a combination of two independent, small amplitude and chirped frequency parametric-type drivings or modulations of the interaction strength in the GP equation allows controlled nonlinear two-phase excitation of BECs via the process of autoresonance.The amplitude of these excitations grows continuously as the phases of the excited solution follow those of the driving perturbations.The phase-locked excitation process starts from trivial ground state at large negative times and continues as the system passes linear resonances with both driving components at t = 0 and moves to large positive times.If both driving components are switched off at some large positive time, an ideal, stable space-time quasicrystal is formed, which is periodic in space and aperiodic in time, but preserves long time ordering (see examples of numerical simulations in Figs.2,4, and 9).
We have developed a weakly nonlinear theory for these two-phase periodic excitations using Whitham's averaged variational principle, yielding a two degrees of dynamical problem in action-angle variables (see Eqs. ( 26)-( 29)).The analysis of this system at large positive times limits the parameter space for autoresonant excitations (see inequalities (41) and (42).However, these are only necessary conditions for autoresonance in the system.A continuing phase-locking by passage through the linear resonance near t = 0 requires the driving amplitude to surpass a certain threshold, yielding the sufficient condition for autoresonant excitation.These thresholds scale with the driving frequency chirp rate α as ∼ |α| 3/4  (see Eq.( 47)), a relationship corroborated by numerical simulations (see Fig. 8).
Following concepts explored in this study, it seems interesting to investigate autoresonant excitation involving more complex multi-phase quasi-crystalline structures by simultaneously traversing linear resonances with three and more parametric drivings.Furthermore, given that numerous integrable nonlinear PDEs (such as KDV, Sine-Gordon, and more) describing various physical systems allow multiphase solutions, investigating autoresonant formation of space-time quasicrystals in these additional systems through two or more drivings using a similar approach seems interesting.Finally, previous investigations [27,28] have demonstrated that a sufficiently small dissipation in other applications did not destroy single-phase autoresonant synchronization, but modified the threshold for transition to autoresonance.Investigating the effects of dissipation on autoresonant two-phase BECs is also an important goal for future research.

ACKNOWLEDGEMENT
This work was supported by the US-Israel Binational Science Foundation Grant No. 2020233.

APPENDIX: REDUCTION TO TWO DEGREES OF FREEDOM
We proceed from the Lagrangian density (10) for the ponderomotive drive case The case of driving via the modulation of the interaction strength will be discussed at the end of the Appendix.We average the Lagrangian over θ 1 and θ 2 between 0 and 2π using the weakly nonlinear ansatz (18) and (19).The resulting averaged Lagrangian density Λ consists of the following five terms +k 2 1 (4u This change affects only the driving term Λ 5 (see Eq. ( 53) in the averaged Laplacian, transforming it to Formally, this is a replacement of U 0 in the coefficient from the ponderomotive drive case by 2σU 3 0 .The same replacement should be done in the case of driving via interaction strength in all driving terms in the reduced system of equations ( 71),(72),(74).

FIG. 4 .
FIG. 4. The passage through threshold.The colormaps of |ϕ| 2 over x in space-time in three examples shown in Fig. 3.

FIG. 7 .
FIG. 7. Allowed double phase-locking region for BECs in the case σ = −1 is below the blue line in the figure.The circle and the square points correspond to parameters k1 = 2, k2 = −3, U 0 = 0.775 and k1 = 1, k2 = −2, U0 = 0.447, respectively.The examples of space-time quasicrystals for these two sets of parameters are shown in Fig. 9