Caustics in quantum many-body dynamics
CCaustics in quantum many-body dynamics
W. Kirkby, Y. Yee, K. Shi, and D. H. J. O’Dell
Department of Physics and Astronomy, McMaster University,1280 Main St. W., Hamilton, Ontario, Canada, L8S 4M1 (Dated: February 2, 2021)Caustics are a striking phenomena in natural optics and hydrodynamics: high-amplitude char-acteristic patterns that are singular in the short wavelength limit. We use exact numerical andapproximate semiclassical analytic methods to study quantum versions of caustics that form inFock space during far-from-equilibrium dynamics following a quench in the Bose Hubbard dimerand trimer models. Due to interparticle interactions, both models are nonlinear and the trimer isalso nonintegrable. The caustics exhibit a hierarchy of morphologies described by catastrophe the-ory that are stable against perturbations and hence occur generically. The dimer case gives rise todiscretized versions of folds and cusps which are the simplest two of Thom’s catastrophes. The extradimension available in the trimer case allows for higher catastrophes including the codimension-3hyperbolic umbilic and elliptic umbilic catastrophes which we argue are organized by, and projec-tions of, an 8-dimensional corank-2 catastrophe known as X . These results indicate a hithertounrecognized form of universality in quantum many-body dynamics organized by singularities. I. INTRODUCTION
The remarkable level of control that can be exertedover ultracold atomic gases has allowed, over the lasttwenty five years, the experimental realization of a num-ber of paradigmatic quantum many-body models [1].Long relaxation times, coupled with the ability to rapidlyvary external potentials and interparticle interactions,means that quantum gases are particularly suited tostudying nonequilibrium dynamics and this has reawak-ened interest in fundamental dynamical questions suchas how quantum systems reach thermal equilibrium [2, 3]and whether quantum dynamics can display universalityakin to that seen near equilibrium phase transitions [5–8].The Bose-Hubbard (BH) model is one of the key mod-els that has been realized. It describes interacting bosonsmoving on a lattice, and has a quantum phase transi-tion between superfluid and insulating phases at a crit-ical ratio of interaction strength to hopping amplitude[9, 10]. This transition was observed for the first timein a landmark experiment in 2002 [11] using ultracoldatoms in an optical lattice [12]. Subsequent theoreti-cal [13–27] and experimental [28–37] work on dynamicsin the BH model has addressed a range of topics in-cluding timescales for many-body quantum revivals andthe establishment of coherence, light-cone-like propaga-tion of correlations, effective Hamiltonians in periodicallydriven “Floquet” systems, and relaxation to equilibrium,to name just a few.The BH dimer and trimer models are particular casesthat consider two and three lattice sites (modes), respec-tively. The dimer describes bosonic Josephson junctions[38–45] that have been realized experimentally with Bose-Einstein condensates (BECs) trapped in double-well po-tentials [46–51] and also with spinor BECs exhibitingthe internal version of the Josephson effect [52]. Ad-ditionally, the same Hamiltonian describes trapped ionswith two internal states and long-range interactions [53–56]. The BH dimer displays a rich variety of behaviour including a transition from oscillatory population dy-namics to macroscopic quantum self-trapping (an excitedstate of the system with a self-maintained populationimbalance) when the ratio between interparticle inter-actions and intersite tunneling exceeds a threshold value[39, 46, 47, 52, 57, 58]. Similar “dynamical phase transi-tions” have been seen in other systems [59–63] and, liketheir equilibrium versions, they obey scaling laws in thecritical regime thereby providing examples of universalityin dynamics. The many-body caustics we will describein this paper are another example of universality in dy-namics characterized by a hierarchy of universal shapesin Fock space that have self-similar scaling [64–66].The BH trimer model describes spin-1 BECs [67–69]where the atoms share a common external state (as ina tight trap), as well as scalar BECs in triple-well po-tentials. Spin-1 BECs have been realized in experimentson Na [70, 71] and Rb [72–77] where the three in-ternal states are provided by the Zeeman sublevels ofthe F=1 hyperfine manifold (to explore the full trimermodel an integrability breaking RF field can be ap-plied [78]). These systems are also predicted to dis-play macroscopic quantum self-trapping [79–81], how-ever, unlike the dimer, the trimer is nonintegrable andits classical dynamics exhibits chaos [78, 79, 82–86] giv-ing behaviour qualitatively closer to the many-site model.Another feature of the trimer is that it allows next-to-nearest-neighbor interactions which are important whenthe atoms have dipole-dipole interactions [87–90].Caustics are best known as an optical phenomenonarising from natural focusing. Everyday examples in-clude rainbows, bright lines on the bottom of swim-ming pools, and twinkling starlight [91]. They dominatewave patterns because they are places where multiplewaves coalesce; in the classical limit caustics are the lo-cus of points where neighboring rays touch and appear asbright, sharp, lines on a screen or surfaces in space wherethe intensity diverges [92, 93]. As understood by Kelvin[94], caustics also occur in hydrodynamics as ship wakes a r X i v : . [ c ond - m a t . qu a n t - g a s ] J a n and tidal bores [95], and have been identified as one of thecauses of freak waves and extreme events [96–100]. Likethe BH model we study here, the hydrodynamic case is ingeneral nonlinear and this has inspired studies compar-ing freak waves in linear and nonlinear optical systems[101–109]. Another field where caustics play a role isastrophysics through gravitational lensing [92], and as adescription of the large scale structure of the universe[110, 111].Several key features of caustics are illustrated by Fig.1, which was made by shining laser beams through waterdroplets which act as natural (imperfect) lenses. Firstly,caustics take on certain characteristic shapes which arethe structurally stable morphologies described by catas-trophe theory (a theory describing singularities in gra-dient maps such as Fermat’s principle) [112]. Structuralstability means stability against perturbations and hencethey occur in the wild, i.e. in typical or generic situ-ations where all symmetries are broken. For example,an isolated point focus arising from a perfect lens is notstructurally stable in two or more dimensions and in-stead evolves into an extended caustic in the presence ofaberrations. Catastrophe theory predicts a hierarchy ofsuch characteristic shapes (catastrophes) which becomemore elaborate as the dimension of space increases, withthe higher catastrophes containing the lower ones [113–115]. Fig. 1 shows 2D slices through 3D caustics in spaceand in both panels (a) and (b) we can see that these 3Dcaustics contain one or more cusp shapes which are thestructurally stable catastrophe in 2D. Secondly, causticsdisplay different structure at different scales. At scaleslarge compared to the wavelength caustics appear as sin-gular lines or surfaces, but zooming down to the wave-length scale we see that caustics are not singular but areinstead softened by wave interference. In fact, each classof catastrophe is dressed by a particular interference pat-tern known as a “wave catastrophe” (in the case of thecusp this pattern is the Pearcey function [116]) and de-fined through a diffraction integral [117]. Furthermore,these patterns display self-similar scaling as the wave-length is changed so that each class of wave catastrophecomes with a set of scaling exponents that determine thescaling of the amplitude and fringe spacing in the dif-ferent possible directions. Zooming down even further,at the finest scales one discovers that the interferencepattern contains a network of nodes that are quantizedvortices that come in vortex-antivortex pairs [119].As a generic wave phenomenon, caustics also occur inquantum waves. They have long been studied in nu-clear scattering [120] and more recently have been ob-served in electron microscopy [121], atom optics [122–124], and in the experiment described in reference [123]a cusp caustic was recorded in the time-dependent atomdensity distribution of a dilute BEC moving in a 1D op-tical lattice. In all these examples the matter waves areadequately described by the single-particle Schr¨odingerwave equation, but caustics are not limited to this caseand also arise in solutions of the Gross-Pitaevskii equa- FIG. 1. Optical caustics made by shining a laser beamthrough a water droplet of radius ∼ Panel (a):
Hyperbolic umbilic,
Panel (b):
Elliptic umbilic,
Panel (c):
Elliptic umbilic near its most singular point. What we seein these photos are slices through three dimensional catas-trophes that also contain lower catastrophes: the hyperbolicumbilic contains a single cusp (which is the only structurallystable catastrophe in 2D and is dressed by the Pearcey func-tion wave catastrophe) and the elliptic umbilic contains threecusps. Taking a 1D slice across a cusp would give the simplestcatastrophe of all, the fold catastrophe (whose wave patternis the Airy function). These images and the methods usedto make them were inspired by the experiments reported inreference [118]. The different colors arise from using threedifferent color lasers. tion (nonlinear Schr¨odinger equation) which can equallydescribe a self-interacting condensate [125] or nonlinearoptics in a fibre [101]. Nevertheless, all these examplesare morally similar to the classical wave scenario found inoptics and hydrodynamics, whether linear or nonlinear.The case of quantum many-body systems is different be-cause in contrast to the mean-field theory described bythe Gross-Pitaevskii equation, the discrete nature of thefield quanta is taken into account.Inspired by ideas from quantum optics [126, 127], wehave initiated a theory of many-body (second quantized)caustics in a series of papers on spin chains [128] and theBH dimer model [64–66, 129]. In the former case, the cel-ebrated “light-cones” observed [32, 55, 56] to govern thepropagation of correlations following a quench are shownto in fact be caustics akin to ship wakes. Of greater rel-evance to the present paper is the latter case where thecaustics live in Fock space. This is a fundamentally dis-crete space and the natural arena for many-body physics.In particular, the dimer displays discretized versions ofthe fold and cusp catastrophes living in (1+1) dimen-sional Fock space plus time, with the role of Planck’s con-stant (or the inverse wavelength) played by 1 /N , where N is the total number of particles. In the semiclassicalregime N (cid:29) N ∼
150 particles, weare unable to obtain analytic mappings onto the canon-ical catastrophe wavefunctions for the trimer because itis non-integrable. In order to provide some analytic ex- amples, in Section VI we instead study δ -kicked dynam-ics where the interactions are flashed on and off suchthat the mapping can be achieved analytically. We notethat interactions can be engineered in cold atom exper-iments using Feshbach resonances [136]. In Section VIIwe go beyond the quantum phase model (rigid pendulummodel), which assumes all three modes are significantlyoccupied, and find corrections that break circular sym-metry in Fock space in favor of triangular symmetry andidentify the particular sub-family of the X catastropheat work. In order to do this we introduce a path inte-gral representation for the wavefunction. In Section VIIIwe compare repulsive and attractive interactions and dis-cuss the crucial role interactions play in the formation ofcaustics in BH dynamics. Finally, in Section IX we giveour conclusions. There are also three appendices whichcontain details of some of the calculations, including aderivation of the path integral. II. TWO- AND THREE-MODE BH MODELSA. Two-Mode Equations of Motion
We first review the BH dimer, which will form up to(1+1)-dimensional caustics in the dynamics. For a morecomplete analysis of the formation of caustics in two-mode models, see Refs. [64, 65]. The BH dimer has theHamiltonian [137, 138]ˆ H dimer = − J (cid:16) ˆ a † l ˆ a r + ˆ a † r ˆ a l (cid:17) + U (cid:16) ˆ a † l ˆ a l − ˆ a † r ˆ a r (cid:17) , (1)where ˆ a ( † ) l/r annihilates (creates) a particle in the left/rightwell, J is the hopping energy and U is the on-site interac-tion energy between atoms in each well. The annihilationand creation operators obey the usual bosonic commuta-tion relations [ˆ a i , ˆ a † j ] = δ ij , where i and j correspond toeither l or r .In this paper we study caustics that form in Fock space.For the dimer case the Fock states are the eigenstates ofthe half-number-difference operator ˆ n ≡ (ˆ a † l ˆ a l − ˆ a † r ˆ a r ) / | Ψ( t ) (cid:105) = (cid:88) n c n ( t ) | n (cid:105) . (2)Inserting Eq. (2) into the time-dependent Schr¨odingerequation, i (cid:126) ∂ t | Ψ( t ) (cid:105) = ˆ H | Ψ( t ) (cid:105) , we obtain a set of N + 1coupled differential equations for the Fock-space ampli-tudes c n ( t ) which we refer to as the generalized Raman-Nath (RN) equations because a similar set of differentialdifference equations were originally derived by Ramanand Nath in the context of dynamical diffraction [139–141],i (cid:126) ˙ c n ( t ) = 4 U n c n ( t ) − Jc n − ( t ) (cid:114) N N − n + n − Jc n +1 ( t ) (cid:114) N N − n − n (3)where the dot represents a time derivative. These arethe equations of motion which must be solved to find thequantum dynamics of the two-mode BH model.The role of geometric rays is played by the mean-fieldsolutions where both the number and phase difference aresimultaneously defined (in the quantum theory only oneof these observables can take a precise value at any mo-ment). The mean-field Hamiltonian can be determinedvia the Heisenberg substitution rules, replacing operatorswith complex amplitudes [83, 142]ˆ a r/l → (cid:113) N r/l e i θ r/l , (4)leading to H dimerMF = 4 U n − J (cid:112) N − n cos φ , (5)where φ = θ r − θ l is the phase difference between thetwo modes and is the conjugate variable to n . H dimerMF de-scribes a classical nonrigid pendulum where n plays therole of angular momentum and φ angular position. Thevariable length of the nonrigid pendulum is accounted forby the square root factor [39]. Solving Hamilton’s equa-tions of motion, better known in this context as Joseph-son’s equations for two coupled superfluids,˙ n = − (cid:126) ∂∂φ H dimerMF = − J (cid:126) (cid:112) N − n sin φ (6)˙ φ = 1 (cid:126) ∂∂n H dimerMF = 8 U (cid:126) n + J (cid:126) n √ N − n cos φ , (7)gives the mean-field trajectories (rays) [143].Attempts to semiclassically quantize the mean-fieldproblem are complicated by the appearance of both num-ber and phase variables in the potential energy term,meaning that the Hamiltonian is not separated into thesum of a ‘kinetic’ term proportional to n and a purely‘potential’ term V ( φ ), but can be pushed through withsome care [43, 144]. However, providing the populationdifference is always small in comparison to the total par-ticle number ( n (cid:28) N ) the square root term in H dimerMF can be set to unity reducing it to that of a standard rigidpendulum. This is the relevant Hamiltonian in super-conducting Josephson junctions where it is known as the quantum phase model (QPM) [145, 146], and also atomicBECs in optical lattices when there are many atoms persite [147–152] . In this paper we shall sometimes makeuse of the QPM for simplicity but will also consider cor-rections to it (we will see in Section VII that this canmake a difference to the caustics that occur). FIG. 2. A small region of Fock space for the BH trimer. Itcan be tiled by hexagonal cells in the coordinates { δn , n X } we use throughout this paper. Hopping terms shift δn byonly 1 unit, while n X ≡ n − n can change by 1 or 2 units. B. Three-mode Equations of Motion
The Hamiltonian for the BH trimer can be written asˆ H = − K L (ˆ a † ˆ a + ˆ a † ˆ a ) − K R (ˆ a † ˆ a + ˆ a † ˆ a ) (8) − K X (ˆ a † ˆ a + ˆ a † ˆ a ) + U (cid:88) i =1 ˆ n i (ˆ n i −
1) + (cid:88) i =1 (cid:15) i ˆ n i . The parameters K L , K R , K X correspond to the hoppingenergies between wells 1 and 2, 2 and 3, and 1 and 3,respectively. The (cid:15) i allow for different well depths.In its linear configuration ( K X = 0), the three-modeBH system has been used as a spatial model for rapid adi-abatic passage by controlling the well depths (cid:15) i ( t ) as func-tions of time [153–158]. Indeed, adding a magnetically-induced tilt to the lattice allows additional rich behaviourincluding control of correlations [167, 168]. The linearconfiguration has also been studied from the point ofview of an ultracold atom transistor-like device [159–162].In its fully-connected triangular configuration, the three-mode system provides a minimal model for superfluidcircuits and discrete vortices [163–166], and the triangu-lar geometry also allows for an extended BH model de-scribing a trapped long-range dipolar gas with nearest-neighbour interactions to be mapped instead to a BHmodel with effective on-site interactions [89]. Further-more, the three-site chain and triangle have been shownto exhibit the possibility of quantum steering [169–172].To the best of our knowledge the three-mode BH modelhas not yet been realized experimentally using, e.g.,triple-well BECs, although detailed experimental propos-als with tuneable hopping and interaction parameters viaFeshbach resonances certainly exist [173]. As mentionedin the Introduction, spin-1 BECs provide another physi-cal system where the three mode BH model might providethe natural theoretical description [78].When U = 0, the mean-field three-mode BH model ex-hibits regular dynamics, while for nonzero interactions itexhibits chaos, thereby indicating nonintegrability. Closeto the ground state these chaotic trajectories are mixedwith islands of regular dynamics [79, 142]. Energy levelstatistics in the quantum version of the trimer tell a sim-ilar story: in general they obey neither the Poisson norWigner distributions but are better described by a Berry-Robnic distribution [174] which is a signature of a clas-sical limit containing both regular and chaotic dynamics[175]. The significance of this for the present work is thatwe find caustics in a nonintegrable model even thoughcaustics are usually associated with integrable behav-ior [176]. We therefore conjecture that the caustics arestable against integrability breaking terms as describedby the famous Kolmogorov-Arnold-Moser theorem [177],although we will leave the detailed examination of thispoint to a future publication.The quantum many-body state can be expanded as avector in the Fock basis as | Ψ( t ) (cid:105) = (cid:88) n n n M n n n ( t ) | n , n , n (cid:105) (9)= (cid:88) n n X M δn n X ( t ) | δn , n X (cid:105) (10)where in the second line we have introduced new Fockspace coordinates δn ≡ n − N/ n X ≡ n − n ,and have assumed the total particle number is conserved,allowing one of the sums to be eliminated (the variables δn and n X are similar in spirit, but different in detail,to those used by Arwas et al [163]). The allowed Fockspace is then triangular in shape and the schematic inFig. 2 depicts a small region of it. Since each Fock stateis coupled to six others by all the possible hopping terms,Fock space can be tiled by a hexagonal pattern as shown.Inserting the general quantum state given in Eq. (10)into Schr¨odinger’s equation we obtain the following gen-eralized Raman-Nath equations for the Fock space am-plitudes M ij ( t ) i (cid:126) ˙ M δn ,n X ( t ) = − K L √ (cid:113)(cid:0) N + δn + 1 (cid:1) (cid:0) N − δn + n X (cid:1) M δn +1 ,n X − − K L √ (cid:113)(cid:0) N + δn (cid:1) (cid:0) N − n + n X + 2 (cid:1) M δn − ,n X +1 − K R √ (cid:113)(cid:0) N + δn (cid:1) (cid:0) N − δn − n X + 2 (cid:1) M δn − ,n X − − K R √ (cid:113)(cid:0) N + δn + 1 (cid:1) (cid:0) N − n − n X (cid:1) M δn +1 ,n X +1 − K X (cid:113)(cid:0) N − n + n X + 2 (cid:1) (cid:0) N − n − n X (cid:1) M δn ,n X +2 − K X (cid:113)(cid:0) N − n − n X + 2 (cid:1) (cid:0) N − n + n X (cid:1) M n ,n X − + U (cid:20) δn + n X (cid:21) M δn ,n X (11) where we have assumed all (cid:15) i = 0.The mean-field approximation is obtained as aboveby replacing operators with complex amplitudes, ˆ a i → √ n i e i φ i . The resulting Hamiltonian is, H MF = − K L √ n n cos( ϕ − ϕ ) − K R √ n n cos( ϕ − ϕ ) (12) − K X √ n n cos( ϕ − ϕ )+ U (cid:88) i =1 n i ( n i −
1) + (cid:88) i =1 (cid:15) i n i . Like in the quantum case, we can change our coordinatesto δn and n X and eliminate the third member due tonumber conservation. However, in the mean-field prob-lem we also require the phase variables conjugate to thenumber variables and these are φ X ≡ ( ϕ − ϕ ) and φ C ≡ (2 ϕ − ϕ − ϕ ), respectively. The third phasevariable Θ ≡ ϕ + ϕ + ϕ is irrelevant to the mean-fielddynamics studied here, and in the quantum wavefunc-tion becomes a global phase. The resulting mean-fieldequations of motion in these variables come out to be: ˙ n X = − K L (cid:113) (cid:0) N + δn (cid:1) (cid:0) N − δn + n X (cid:1) sin ( φ X − φ C ) − K R (cid:113) (cid:0) N + δn (cid:1) (cid:0) N − δn − n X (cid:1) sin ( φ X + φ C ) − K X (cid:113)(cid:0) N − n (cid:1) − n X cos (2 φ X ) (13)˙ δn = K L (cid:113) (cid:0) N + δn (cid:1) (cid:0) N − n + n X (cid:1) sin ( φ X − φ C ) − K R (cid:113) (cid:0) N + δn (cid:1) (cid:0) N − δn − n X (cid:1) sin ( φ X + φ C )(14)˙ φ X = U n X − K L (cid:0) N + δn (cid:1) cos( φ X − φ C ) (cid:113) (cid:0) N + δn (cid:1) (cid:0) N − δn + n X (cid:1) + K R (cid:0) N + δn (cid:1) cos( φ X + φ C ) (cid:113) (cid:0) N + δn (cid:1) (cid:0) N − δn − n X (cid:1) + K X n X cos(2 φ X ) (cid:113)(cid:0) N − δn (cid:1) − n X (15)˙ φ C = 3 U δn − K L (cid:0) N − δn + n X (cid:1) cos( φ X − φ C ) (cid:113) (cid:0) N + δn (cid:1) (cid:0) N − δn + n X (cid:1) − K R (cid:0) N − δn − n X (cid:1) cos( φ X + φ C ) (cid:113) (cid:0) N + δn (cid:1) (cid:0) N − δn − n X (cid:1) + K X (cid:0) N − δn (cid:1) cos(2 φ X ) (cid:113)(cid:0) N − δn (cid:1) − n X . (16) While the two-mode case gave mean-field equations de-scribing a non-rigid pendulum, the three-mode case canbe interpreted as describing three coupled anharmonic os-cillators [142]. Eqns. (13)-(16) must in general be solvednumerically, but a problem can potentially arise for tra-jectories that touch the boundaries of Fock space wherethe square root factors in the denominators vanish. Phys-ically, the boundaries correspond to situations where oneof the modes is empty. We find empirically that this be-comes less of a problem as N is increased and almostnever occurs in the semiclassical regime we consider in Catastrophe Symbol n Q Φ Q ( s ; C )Fold A s + Cs Cusp A s + C s + C s Swallowtail A s + C s + C s + C s Butterfly A s + C s + C s + C s + C s HyperbolicUmbilic D +4 s + s + C s s + C s + C s EllipticUmbilic D − s s − s + C ( s + s )+ C s + C s ParabolicUmbilic D s + s s + C s + C s + C s + C s TABLE I. Thom’s seven elementary catastrophes, their sym-bols, and generating functions Φ Q ( s ; C ), organized by corank n , and codimension Q [112]. this paper where N ∼
150 because trajectories spendmost of their time in the central region of Fock space.
III. WAVE CATASTROPHES
The seven elementary catastrophes originally discussedby Ren´e Thom [113] are listed in Table I, and some exam-ples of higher catastrophes are listed in Table II. Theyare organized by their corank (number of “state” vari-ables s = { s , s , . . . } that label paths) and codimension(number of control parameters C = { C , C , . . . } whichin our case is the dimension of Fock space plus time andany other parameters in the Hamiltonian). The key ob-jects are the generating functions Φ Q ( s ; C ) each of whichcorresponds to a particular catastrophe. The generatingfunctions act as potentials for gradient maps, and in ourphysical application they give the local action close tothe caustic. Stationary points ∂ s Φ Q = 0 specify classicalpaths or rays, which are the mean-field solutions. Caus-tics occur where the action is stationary to higher order,i.e. they are regions where classical paths either coalesceor are born (bifurcations). Remarkably, this can onlyhappen in certain ways if the bifurcation is to be struc-turally stable against perturbations which is what leadsto the distinct catastrophes. The art of applying catas-trophe theory lies in finding coordinate transformationsthat map the physical action onto one of the generatingfunctions.Each catastrophe has a ‘germ’, which is the part of thegenerating function that remains when it is evaluated atthe origin of control space C = 0. The germ characterizesthe order of the singularity. The other terms show howthe catastrophe ‘unfolds’ as one moves away from theorigin of control space. The corank 1 catastrophes, A Q +1 are called the cuspoids , and extend beyond the butterflyto the wigwam and star catastrophes (not listed). Theremaining catastrophe types, D Q +1 , E Q +1 , and above aretypically called umbilics , referring to the classification ofcubic forms near an umbilic point (a point on a surface Symbol n Φ Q ( s ; C ) A Q +1 s Q +2 + (cid:80) Q − i =1 C i s i D ± Q +1 s Q ± s s + C Q +1 s + (cid:80) Q − i =2 C i s i + C s E s + s + C s s + C s + C s s + C s + C s E s + s s + C s + C s + C s + C s s + C s + C s E s + s + C s s + C s s + C s + C s s + C s + C s + C s X ± s + Ks s ± s + C s s + C s s + C ( s + s ) + C ( s − s )+ C s s + C s + C s TABLE II. Catastrophe organization beyond Thom’s sevenelementary list [92], with general control space dimension Q .Many of the higher catastrophes do not have names, and arereferenced by their group-theoretic symbol, yet are often splitinto cuspoids A Q +1 , and umbilics { D Q +1 , E Q +1 , X } . with locally spherical curvature) which become the germsfor these catastrophes.Catastrophes obey projection identities: higher catas-trophes contain those of lower order. For example, theswallowtail catastrophe contains two cusps and three foldlines when projected into two dimensions. It is not, how-ever, guaranteed that catastrophes of high order (i.e.higher codimension and/or corank) contain all catastro-phes of lower order. In section VI, we will briefly discussdistinctions between families of the high-order catastro-phe X , which have different projection identities ( X seems to be the catastrophe that organizes all the struc-tures we see in the BH trimer dynamics). It is worthemphasizing that with special initial conditions and tun-ing of parameters, one could engineer a focusing eventwhich would result in any shape of caustic (not just theones covered by catastrophe theory). Catastrophe theoryinstead describes structurally stable caustics that resultfrom natural focusing, and so are more likely to appeargenerically.While geometric ray theory (mean-field theory) is suf-ficient to describe the shape of the caustic, it gives riseto a probability singularities in Fock space. To removethese one should include wave interference and we enterthe realm of wave catastrophes or diffraction integrals[92, 112]. Each ray catastrophe is dressed by a character-istic wave interference pattern described by a diffractionintegral, Ψ( C ) ∝ λ n/ (cid:90) ... (cid:90) d s e i λ Φ Q ( s ; C ) . (17)This wavefunction lives in the space of control parame-ters. It resembles a path integral where the generatingfunction plays the role of the action and we integrateover state variables s which label paths. The parame-ter λ acts as the inverse of Planck’s constant, and in theBH model is proportional to the total number of par-ticles N . The precise connection to path integrals willbe explained in Section VII and Appendix C. However,while wave theory removes geometric singularities it alsointroduces new ones, namely phase singularities. Theseare more commonly known as dislocations in optics andvortices in condensed matter systems. It turns out thatthe diffraction patterns given by Eq. (17) are organizedon the finest scale by a network of vortex lines (or pointvortices in two dimensions) where the phase is undefined(takes all values) [178].The term quantum catastrophe has been used by Leon-hardt [126] (see also Berry and Dennis [127]) to describethe blurring of optical phase singularities by quantumfield effects (vacuum fluctuations), although without a di-rect connection to catastrophe theory (the same term hasalso been applied in a different way to the Dicke model[179]). In our case, second quantization of the matterfield leads to discretized versions of wave catastrophes inFock space where the discretization also resolves phasesingularities by turning them into discrete vortices [66].The caustics we encounter in many-body dynamics aretherefore quantum catastrophes in a direct sense. In thelimit of N → ∞ we might expect to return to a continu-ous wave theory just as Maxwell’s equations describe theelectromagnetic field in the limit of large photon number.However, it is important to distinguish two limits in Fockspace:1. The wave-like theory where we preserve commu-tation relations between operators (interferencefringes preserved)2. The geometric-ray theory where we let the commu-tators vanish (interference fringes removed)In Fig. 3 we compare the mean-field (ray) and fullyquantum (second quantized) theories for a Fock-spacecaustic in the 3-mode BH model. Although present, thediscretization of Fock space is not visible and the quan-tum waves appear smooth. The mean-field theory weapply is the truncated Wigner approximation (TWA)in which an ensemble of classical trajectories is prop-agated in time using the classical equations of motionEqns. (13)–(16) with initial conditions sampled from aquantum quasiprobability distribution (the Wigner func-tion) [180–183]. The sum of these trajectories gives theproper mean-field approximation to the quantum dynam-ics. The initial state in Fig. 3 is a phase state, whichis a state with narrow relative phase distributions butwhich consequently has a flat probability distribution inFock space because number and phase are conjugate vari-ables. In the mean-field case shown in panel (a), wesee that the first time slice contains a representative setof points approximating an equal superposition of Fockstates. Panel (b) plots the absolute values of the quantumamplitudes of the Fock states found by solving the gener-alized Raman-Nath equations Eq. (11). In both cases thedynamics leads to focusing and clearly forms an ellipticumbilic caustic which can be compared to those shown inFig. 1. An equal superposition of Fock states was chosen FIG. 3. Time slices comparing mean-field and quantum dy-namics in the triangular Fock space found in the BH trimerunder the constraint of total number conservation.
Panel(a):
An ensemble of initial points are evolved in time usingthe classical equations of motion Eqns. (13)–(16) and forma space-time version of the elliptic umbilic catastrophe, withone particular trajectory picked out for illustration. In theTWA the initial conditions are drawn from a quantum prob-ability distribution and in this simple case the initial stateis a phase state meaning that the relative phases are sharplydefined (here taken to be φ X = φ C = 0) but with the conse-quence that their conjugate number differences take all possi-ble values with equal probability, thereby uniformly populat-ing the allowed triangular region of Fock space. Panel (b):
The quantum dynamics are obtained by solving Eq. (11). Theinitial phase state appears as a discretized plane wave in Fockspace (although here N = 150 and the discreteness of Fockspace is not visible). The correspondence between the quan-tum and classical dynamics is clear: the quantum wavefunc-tion is brightest where the density of trajectories is highest.In both images, K L = K R = K X ≡ J , U = − . J , and (cid:15) i = 0. FIG. 4. The wave catastrophes associated with the foldand the cusp.
Panel (a):
The fold catastrophe is decoratedby an Airy function, Eq. (19). The location of the classicalcaustic is indicated by a vertical dotted line. In the classical(ray) theory the intensity is divergent at this point and fallsoff as 1 / √− C on the bright side. By contrast, the Airy func-tion is finite at the caustic and two-wave interference givesrise to oscillations; these two waves coalesce at C = 0 andbecome evanescent on the dark side. Panel (b):
The cuspcatastrophe is decorated by the Pearcey function which is acomplex-valued function (here we plot the modulus) given byEq. (21), with the divergent classical cusp caustic shown as asolid black curve. because it is a plane wave-like state analogous to the ini-tial state often considered in optics, but here it is the BHdynamics that acts as an imperfect lens which focusesthe wave in Fock space. The key feature of this initialstate is that it is broad in Fock space, and structuralstability means that the caustics it generates will notbe qualitatively different from those generated by otherbroad states such as the gaussian-shaped ground state inthe case with strongly coupled sites such that the hop-ping dominates interactions (this case will be discussedin Section VI B).In the next two sections we investigate the hierarchy ofwave catastrophes that appear in BH dimer and trimerdynamics, building up to the high-order catastrophe X which ultimately organizes the lower catastrophes we see.It should be borne in mind that because the trimer is notintegrable our ability to analytically describe the appear-ance of each catastrophe is limited and for this reason wego to δ -kicked dynamics in Section VI. Furthermore, un-like free-space optics, classical rays in Fock space do nottravel in straight lines in the BH model even for the in-tegrable case ( U = 0). IV. CAUSTICS IN THE DIMER
Dynamical caustics in the BH dimer live in the (1+1)Dspace formed by Fock space and time. The only struc-turally stable catastrophes in two dimensions are foldlines which can meet at cusp points.
A. Fold
The simplest catastrophe is the fold. Folds are corank-1, codimension-1 objects with a cubic generating func-
FIG. 5. Cusp caustics in the two-mode BH model.
Panel(a):
Each curve is a mean-field configuration obtained fromEqns. (6)–(7) starting from n = 0 and with an initial phasesampled from the distribution φ = [0 ... π ) corresponding toa definite (equal) number of bosons in each well. Panel (b):
Modulus of the quantum wavefunction calculated using thetwo-mode RN equations Eq. (3) with N = 400 and wherethe initial state is the single Fock state | n = 0 (cid:105) . In bothcases we can identify a series of cusp caustics correspondingto partial revivals of the initial state. In the inset we see thatin the immediate vicinity of each cusp the wavefunction is aPearcey function [64, 66] (the discretization of Fock space isnot visible here so that wavefunction appears continuous). tion, Φ ( s ; C ) = s + Cs . (18)Folds arise where two families of rays coalesce and thiscan be at a point on a line, a line in a plane, a surface in3D etc. The corresponding wave catastrophe is2 π / Ai (cid:18) C / (cid:19) = (cid:90) ∞−∞ d s e i( s + Cs ) (19)and is plotted in Fig. 4(a). This function is the well-known Airy function introduced as the wave descrip-tion for light at rainbows in 1838 [184]. It not onlyremoves the singularity in the ray theory but its inter-ference fringes also explain the supernumerary arcs thatare sometimes visible inside the main bow in optical rain-bows. B. Cusp
In two dimensions, generic focusing events are fold linesthat meet at cusps. These have quartic generating func-tions, Φ ( s ; C , C ) = s + C s + C s . (20)Inside the cusp, three families of classical rays coexist,while only one family exists outside. This means thatalong each fold line, two sets of rays coalesce, while atthe highly singular cusp point, all three families coalesce.The wave catastrophe associated with the cusp is a 2Dwavefunction known as the Pearcey function [116, 117],Pe[ C , C ] = (cid:90) ∞−∞ d s e i( s + C s + C s ) , (21)and is plotted in Fig. 4(b). As can be seen, the Pearceyfunction consists of an interference pattern inside thecusp (due to three-wave interference), while outside itbecomes exponentially suppressed. Since the cusp catas-trophe consists of the meeting of two fold lines, a one-dimensional slice of the Pearcey function across one ofthese lines projects onto an Airy function.We show examples of these two wave catastrophes ap-pearing in Fock-space dynamics of the BH dimer follow-ing a quench in Fig. 5. The initial condition is a singleFock state | n = 0 (cid:105) (which is the opposite case to thephase state shown in Fig. 3). Physically, it describes thesituation where two independent BECs are suddenly cou-pled so that at t = 0 particles can begin to hop betweenthem, thereby building up coherence [185]. The struc-tural stability of catastrophes means that qualitativelysimilar behaviour is found for similar initial states suchas narrow gaussians in Fock space that can also be cen-tered away from n = 0. Fig. 5(a) shows the set of TWAtrajectories propagating from n = 0, each with a differ-ent relative phase φ . The complete certainty in n meanscomplete uncertainty in φ so the phases are drawn withequal probability from the interval φ = [0 ... π ). A seriesof cusps form at times Jt = mπ/ (2 (cid:112) N U/J ) for m = 1 , , , ... , due to imperfect focusing of trajectories.This is caused by the nonlinearity due to interactions:setting U = 0 leads to isolated focal points. In the quan-tum theory a set of repeating Pearcey patterns dress eachcusp. Moving away from the cusp tip the wavefunctionrapidly tends to back-to-back Airy functions describingthe two fold lines emanating from the cusp [64]. V. CAUSTICS IN THE TRIMER
Dynamics in the three-mode BH model can also yieldcusps if we restrict attention to just one of the Fock spacevariables, as illustrated in Fig. 6. The initial conditionwas this time chosen to be an equal superposition of allFock states. Like an individual Fock state, this does notcorrespond to an eigenstate of the Hamiltonian that weuse to propagate the system in time. Instead it is ahighly-excited state made up of a broad superpositionof eigenstates, and this allows the system to explore itsnonlinearity and produce clear caustics.When we examine the full (2+1)D space availablein the trimer then, as expected, we discover thecodimension-3 catastrophes, namely the elliptic and hy-perbolic umbilics. This is shown in Figs. 7 and 8 which
FIG. 6. Recurring cusps in the three-mode BH model, start-ing from an equal superposition of Fock states, the three-modeversion of a phase state.
Panel (a):
Classical trajectoriesalong the plane of fixed n X = 0. Panel (b):
Quantum dy-namics in the same plane as (a) where the classical cusps arenow dressed by interference fringes, locally approximated bythe Pearcey function. The bright streaks are finite-size effectsdue to reflections off the Fock-space boundary in the n X di-rection, but structural stability ensures the cusps survive. Inboth panels, K L = K R ≡ J = 100 U , K X = 0 (linear configu-ration), while for the quantum dynamics we used N = 180. compare numerical solutions of the generalized Raman-Nath equations given in Eq. (11) with the canonical wavecatastrophes. According to Table I, there is a thirdcodimension-3 catastrophe known as the swallowtail, butthis is a corank 1 catastrophe meaning that it only hasa single state variable. To realize the swallowtail catas-trophe in the trimer we could freeze one of the conjugatephases { φ X , φ C } , but we will not pursue this possibilityhere.It is important to emphasize that although we havechosen parameters where the match in Figs. 7 and 8 be-tween the numerical results and the canonical catastro-phes is quite good, we have not optimized the parame-ters and these catastrophes occur generically. If we couldsolve the trimer model analytically we could perform anexact mapping but due to its nonintegrability this is notpossible. We must therefore satisfy ourselves with quali-tative rather than quantitative matches which are in anycase fully within the spirit of catastrophe theory which isa topological theory [113–115]. In Section VI we considersimplified kick dynamics where we can precisely matchthe physical wavefunctions onto the canonical wave catas-trophes.0 BH Trimer Wave Catastrophe
FIG. 7. Comparison of BH trimer dynamics with the canon-ical elliptic umbilic catastrophe.
Panel (a):
Fock-space am-plitudes starting from an even superposition of Fock states[Eq. (28)] for K L = K R = K X ≡ J , U = 0 . J and N = 150at Jt/ (cid:126) = 0 . Panel (b):
Diffraction pattern amplitude inthe plane around the elliptic umbilic focus.
Panel (c):
Sameconditions as panel (a), but now at
Jt/ (cid:126) = 0 . Panel(d):
Diffraction pattern for a slice of the elliptic umbilicat C = 3 . Panel (e):
Caustic surface of the ellipticumbilic catastrophe. The triple-cusped intersection with a C = const. plane is highlighted in red. A. Elliptic Umbilic
The elliptic umbilic catastrophe,Φ D − ( C ) = 3 s s − s + C ( s + s ) + C s + C s , (22)is one of the two catastrophes of codimension 3 withcorank 2. As can be seen in Fig. 7, the dominant featureof the elliptic umbilic is its three-fold symmetry. Takingtwo-dimensional slices at fixed values of C (which forus is the time direction), the elliptic umbilic appears asthree curved fold lines that meet at three cusps. In itsfull three dimensional form we see that these are reallythree fold surfaces that meet at cusp shaped ribs. In thegeometric theory, there are four rays at all points insidethe caustic, while only two exist at any point outsidebecause two coalesce on the caustic.The corresponding diffraction catastrophe has beenstudied both theoretically and experimentally by Berry et al. in Ref. [118] by focusing light through a trian-gular water droplet lens (as repeated by us in Fig. 1). BH Trimer Wave Catastrophe
FIG. 8. Comparison of BH trimer dynamics with the canon-ical hyperbolic umbilic catastrophe.
Panel (a):
Wavefunc-tion amplitudes for hopping strengths K L = K R ≡ J = 4 U , K X = 0, and N = 150 at Jt/ (cid:126) = 0 .
24. The initial state was | δn , n X (cid:105) = | , (cid:105) . Panel (b):
Hyperbolic umbilic diffractionpattern on the plane C = 0. Panel (c):
Same conditionsas panel (a) now at
Jt/ (cid:126) = 0 .
34. The bottom umbilic is nowcompletely unfolded.
Panel (d):
Diffraction pattern for thehyperbolic umbilic in the plane C = 3. Panel (e):
Causticsurface of the elliptic umbilic catastrophe. A projection ontothe plane C = const. is highlighted in red. A Pearcey diffraction pattern locally dresses each cusp,as apparent by considering any particular corner of Fig.7(d). The three-fold cusp structure expands/contractsas the control parameter C is changed until it collapsesat C = 0. The most singular part of the caustic, theumbilic focus, is at the centre of control space ( C = ).In the wave theory the focus is dressed by a diffractionpattern composed of a three-fold symmetric fork with thebrightest patch at the center and surrounded by an Airyfringe pattern, as shown in Fig. 7(b). B. Hyperbolic Umbilic
The hyperbolic umbilic catastrophe,Φ D +4 = s + s + C s s + C s + C s , (23)is the remaining catastrophe of codimension 3 withcorank 2. The hyperbolic umbilic caustic surface cor-responds to an overlapping cusp and fold extended into1three-dimensional space, as shown in Fig. 8(e). Withinthe cusp, there are four classical rays at every point,two of which annihilate each other as the cusp surfaceis crossed. The remaining two rays annihilate as the foldis crossed. The resulting diffraction pattern projectedon a plane of constant C is a Pearcey-like function sur-rounded by an Airy fringe pattern as seen in Fig. 8(d).At the plane C = 0, the caustic is only partially un-folded, and now two sets of fold lines overlap to form aright-angle corner, dressed by a pattern described by aproduct of Airy functions in two dimensions, as shown inFig. 8(b).We can observe the hyperbolic umbilic in the dynam-ics of the trimer system by starting from the single Fockstate | δn , n X (cid:105) = | , (cid:105) and with K X = 0 (linear spa-tial configuration of the triple-well). The initial compactstate spreads non-uniformly into a two-fold symmetricpolygon with partially unfolded hyperbolic umbilic cor-ners as shown in Fig. 8(a). As the dynamics continue,the edges unfold completely, and the cusps separate fromthe fold lines (moving off the C = 0 plane), as is firstvisible in the bottom corner of Fig. 8(c). C. The X catastrophe In the absence of tilts ( (cid:15) i = 0), the BH trimercan be seen to have six independent control parame-ters: { U, K R , K X , δn , n X , N K L t } (or transformationsthereof), thus the codimension-3 catastrophes that wehave discussed are merely projections of a catastropheembedded in a higher dimensional space. We shall ar-gue in this and the following sections that the highercatastrophe that organizes the BH trimer dynamics isin fact the high order umbilic catastrophe known by itsgroup-theoretic symbol, X . This complicated object haspreviously been the subject of detailed theoretical anal-ysis by Borghi [186] and by Berry and Howls [187], andplays an important role in optical refraction through two-dimensional surfaces, such as water droplets [188], glassjunctions [189], gravitational lensing [92], and has alsobeen discussed in the context of stochastic resonance intwo dimensions [190]. X acts as an organizing centrefor a multitude of lower catastrophes and we refer thereader to Fig. 15 for a “bordering” diagram showing itsrelationship to these subcatastrophes.The X catastrophe is of corank 2 and its complete8-dimensional unfolding can be written as [92],Φ X ± = s + Ks s ± s + C s s + C s s + C ( s + s )+ C ( s − s ) + C s s + C s + C s . (24)Although X appears 8-dimensional, its control space hasonly seven parameters, while the eighth, K , is known asthe modulus . Catastrophes beyond codimension 5 (orabove corank 3) can contain moduli which are differentfrom regular control parameters in that they can only al-ter the caustic pattern geometrically [135] (rather than FIG. 9. 6-cusped diffraction pattern surrounded by a foldline, typical of systems with 3-fold symmetry also organizedby X . Panel (a):
Fock space amplitudes after a quenchstarting from all particles evenly distributed in each well, | δn , n X (cid:105) = | , (cid:105) . Here, N = 150, K L = K R = K X ≡ J , U/J = 0 .
04 at
Jt/ (cid:126) = 0 . Panel (b):
As panel (a), nowwith K R = J , K L = 1 . J , and K X = 0 . J , showing how anabsence of symmetry does not destroy the caustics. topologically by changing the number of critical points[188]), and cannot be removed via scaling arguments.This modulus has two excluded values: K (cid:54) = ±
2, wherethe singularities achieve infinite codimension, meaningthere are an infinite number of ways to unfold the sin-gularity, and perturbations can lead to any number ofcoalescing critical points. As we shall see in the follow-ing sections, the X germ Φ X ± [ C = arises naturallyas the base singularity in all the cases we study.The high dimension of X makes it hard to visualizeand we shall therefore concentrate on particular projec-tions. An example from BH trimer dynamics is shownin Fig. 9(a) where all the hopping amplitudes are equaland leads to two superimposed elliptic umbilics slightlyrotated from each other. More precisely, the character-istic features are a three-fold symmetric caustic featur-ing a six-cusped figure encapsulated by Airy-like fringes.This type of caustic is commonly encountered when light2passes through liquid drops, where a 2 n -cusped patternis formed by drops with n -fold symmetry, all surroundedby an ‘oval’ fold line [188]. Perturbing the symmetry willnot destroy this caustic structure, as shown in Fig. 9(b)where the symmetry of the hopping terms is broken. Itcan only be altered dramatically by changing the sym-metry in a fundamental way, such as setting a hoppingterm to 0.It may appear that the control space of X containsone too many parameters for the BH trimer. This isbecause the physical constraint of number conservationrestricts us from exploring the full mathematical controlspace. We shall see [Fig. 10] that loosening the restriction n + n + n = N allows us to access a complete sectionof the X catastrophe since a point in Fock space is thenspecified by three coordinates rather than two, bringingthe total number of control parameters to seven.The case K = +2 in Eq. (24) does not give a properstructurally stable catastrophe but is conceptually im-portant because it arises in models with circular symme-try, wherein a perfectly circular ‘spun cusp’ is realized,punctured by an unstable axial caustic line. This is infact the case we shall find within the QPM and for a per-fectly triangular trimer (all hopping amplitudes equal)to be discussed in Section VI B. The value K = − X into two distinct sub-families X and X which will come up in Section VII when we go be-yond the QPM. The catastrophe germ for K = 0, whichfor C = C = C = 0 is equal to Φ ( s ) + Φ ( s ), hasled to the X family being known as the ‘double cusp’,sometimes even for K (cid:54) = 0 [191–193]. VI. KICKED DYNAMICS
The patterns and shapes exhibited in Figs. 6–9 are ex-act numerical solutions for the BH trimer model which isin general analytically intractable. They clearly resem-ble the caustics that catastrophe theory predicts, andalso have the expected properties. However, it would bereassuring to have an analytical demonstration that insome tractable limit we really can map the dynamics tocatastrophes. This is what we now do using a simplifiedkicked Hamiltonian. This not only allows us to analyt-ically realize various versions of the X catastrophe butalso suggests an experimentally viable method for engi-neering precisely defined caustics.Following an optical analogy where caustics are formedafter light passes through ‘bad’ lenses which deform thewavefront (an ideal lens will produce a perfect hemispher-ical wavefront that results in a point focus), the role of δ -kicking is to produce a wavefront in Fock space that hasdistortions that upon further propagation will generatecaustics like those we have seen for the full Hamiltonian.Since catastrophes are stable to perturbations the caus-tics we find in the kicked case will survive under moregeneric conditions. In the examples that follow, the interaction term pro-portional to U will be flashed on instantaneously at t = 0 and the system will afterwards evolve purely un-der the hopping terms. In ultracold atom systems time-dependent manipulation of interactions can be achievedusing a Feshbach resonance [136]. In general, when pre-sented with Hamiltonians of the form,ˆ H δ ( t ) = ˆ H + ˆ H δ ( t ) , (25)time-evolution can be achieved via the Floquet operator,ˆ F ≡ e − i (cid:126) ˆ H t e − i (cid:126) ˆ H . (26)In previous work [65], we demonstrated the presenceof catastrophes of codimension 2 in (1+1)D kicked sys-tems. Here, we employ and generalize this frameworkto higher-dimensional catastrophes in the triple-well sys-tem, in particular to the different available unfoldings of X . The δ -kick is sometimes, but not always, necessaryto see caustics in each of the cases we study in this sec-tion, since Fock space trajectories can still be focused bythe hopping terms, a fact which can be demonstrated byignoring ˆ H . We include the kick as part of the calcula-tions since it is more general and it allows the ˆ H termto mimic a tuneable ‘lens’ as mentioned above.We begin by considering the mean-field Hamiltonian, H δ QPM = − JN φ X − φ C ) − JN φ X + φ C ) − K X N φ X ) + δ ( t ) ˜ U (cid:2) δn + n X (cid:3) , (27)which describes a kicked trimer within the quantumphase model and with a tuneable K X hopping term. Notethat the interaction strength in the kicked model has theunits of (cid:126) and is therefore given the symbol ˜ U in Eq. (27)in order to distinguish it from the original interaction en-ergy U in Eq. (8). We recall that the QPM is a valid ap-proximation to the BH model when the mode occupationnumbers are large, so that the square root factors can beneglected from H MF in Eq. (12). In fact, we shall seethat the QPM gives circularly symmetric caustics andwill consider small corrections to the QPM later on inSection VII in order to unfold the X catastrophe that isorchestrating the dynamics from the shadows. To quoteNye [188]: ‘...generic unfoldings may be best understoodas perturbations of symmetrical ones’.The integrability of the kicked model is evident fromthe fact that analytic solutions for the classical trajec-tories can be found and are given in Appendix A. Toobtain quantum dynamics under the Floquet operatorwith H δ QPM we re-promote the observables in Eq. (27)to operators such that they obey the Dirac number-phasecommutators [ ˆ φ, ˆ n ] = i.3 A. Linear trimer
We will first consider the dynamics of the linear BHtrimer which means we set K X = 0 in Eq. (27), andfurther specialize to starting from an equal superpositionof Fock states | Ψ(0) (cid:105) = (cid:88) δn (cid:48) ,n (cid:48) X | δn (cid:48) , n (cid:48) X (cid:105) . (28)This state is a plane wave in the number difference basisbut corresponds to a single phase state, | φ X = 0 , φ C = 0 (cid:105) .Being more concrete than we have in the earlier parts ofthis paper, we define the phase states as eigenstates ofthe phase operators and are the three-mode generaliza-tion of phase difference Bargmann states studied in Ref.[144]. These phase states are overcomplete and thus notstrictly orthogonal for finite N [142], but in what fol-lows we operate on the assumption that N (cid:29) | Ψ( t ) (cid:105) = (cid:88) δn (cid:48) ,n (cid:48) X e − i (cid:126) ˆΦ t e − i (cid:126) ˜ U (3 δn (cid:48) + n (cid:48) X ) | δn (cid:48) , n (cid:48) X (cid:105) , (29) where ˆΦ = Φ( ˆ φ X , ˆ φ C ) is the part of the Hamiltoniancontaining all the phase operators. To obtain the wave-function in Fock space we project | Ψ( t ) (cid:105) onto the Fockbasis by applying (cid:104) δn , δn X | and insert a resolution ofidentity = (cid:80) φX,φC | φ X ,φ C (cid:105)(cid:104) φ X ,φ C | between the exponentialterm and the ket | δn (cid:48) , n (cid:48) X (cid:105) . Making use of relations suchas (cid:104) φ X , φ C | δn (cid:48) , n (cid:48) X (cid:105) = exp[ − i( n (cid:48) X φ X + δn (cid:48) φ C )] and eval-uating Gaussian integrals over the variables n (cid:48) X and δn (cid:48) which appear at most quadratically, we arrive at Eq. (30).Next, since the phase variables { φ X , φ C } are localizedaround zero, at least for short times, we expand the co-sine terms to fourth order. Finally, under a change ofvariables, φ C → / φ C and φ X → / φ X , we obtainEq. (31) which is in a form recognizable as the diffrac-tion integral of X , although the symmetry of the QPMHamiltonian and the initial state we have chosen restrictsthe unfolding so that the terms C , C , and C do notyet appear (the remaining symmetries will be broken inSection VII). Note that both Eqns. (30) and (31) havebeen written in a slightly more general form than neces-sary for the linear trimer by including the K X term sothat they also apply to the triangular case discussed inSection VI B.Ψ( δn , n X , t ) = 4 π (cid:126) i √ U (cid:88) φ X ,φ C exp (cid:20) i 2 N Jt (cid:126) (cid:18) cos ( φ X − φ C ) + cos ( φ X + φ C ) + K X J cos (2 φ X ) (cid:19) + n X φ X + δn φ C + (cid:126) U ( φ C + 3 φ X ) (cid:21) (30) ψ ( δn , n X , t ) = A ( t ) (cid:90) (cid:90) d φ X d φ C exp (cid:34) i N Jt (cid:126) (cid:18) φ X + Kφ X φ C + φ C + αφ X + βφ C + ζφ X + ηφ C (cid:19)(cid:35) (31)For the linear trimer, the modulus in Eq. (31) comesout to be K = 6, and the control parameters are givenby α = 3 √ (cid:126) N Jt ˜ U − √ β = √ (cid:126) N Jt ˜ U − √ ζ = 2 / √ (cid:126) N Jt n X (34) η = 2 / √ (cid:126) N Jt δn (35)with, A ( t ) = 4 √ π (cid:126) U e − i π e i NJt (cid:126) . (36) The wave catastrophe described by Eq. (31) might ap-pear four dimensional, with coordinates { α, β, ζ, η } , butexamination of Eqns. (32) and (33) reveals that α and β are not independent: α = 3 β + 4 √
2, and thereforeit is really three dimensional. Let us consider the par-ticular case α = − β which occurs naturally at the time t = (cid:126) / ( JN ˜ U ). This gives a three-dimensional sectionof X that is hyperbolic and described by the generatingfunctionΦ HX = s + 6 s s + s + C ( s − s ) + C s + C s , (37)which has previously been studied by Berry and Howlsand is relevant to liquid-droplet lenses [187]. The cor-responding theoretical caustic surface is plotted in panel(a) of Fig. 10 where a two-dimensional section at α = C = √ t = (cid:126) / ( JN ˜ U ), then α (cid:54) = − β leading to the intro-duction of a C term.4 FIG. 10. Hyperbolic sections of X in the classical and quantum BH trimer. Panel (a):
Caustic surface for Φ HX , as givenin Eq. (37). The red highlights mark the intersection with the plane C = √ Panel (b):
Distribution of points arisingfrom classical trajectories of the δ QPM Hamiltonian Eq. (27) (see Appendix A for solutions of Hamilton’s equations) for thelinear trimer ( K X = 0) at t = (cid:126) / ( JN ˜ U ) and α = C = √
2, starting in an equal spread of number differences, and with phasedifferences φ C (0) = φ X (0) = 0. The resulting caustic is a partial section through the hyperbolic unfolding of X , restrictedby the triangular Fock space (black solid lines) of physical classical paths. Panel (c):
Same snapshot as in panel (b), nowincluding some unphysical paths (such as those which have started or ended with n X > N ), to demonstrate how the causticis restricted by the shape of Fock space (number conservation). Panel (d):
Wavefunction amplitude at the same moment aspanels (b) and (c), starting from an equal superposition of Fock states (28) and evolved using the RN equations correspondingto the δ QPM model (square root factors set to unity). Here, N = 150, and ˜ U/ (cid:126) = 0 . The actual caustic formed by the classical dynamics[with Hamiltonian Eq. (27)] is plotted in Fig. 10(b). Thisclearly resembles the analytical prediction but only halfof it is present. This is because of the restriction of phys-ical paths to always lie within the triangular Fock space − N ≤ δn ≤ N and − N ≤ n X ≤ N . If unphysi-cal paths are allowed (i.e. non-number conserving pathssuch as those with | n X | > N or n > N ), the full causticis captured, as shown in Fig. 10(c). Panel (d) shows theresulting quantum amplitudes in Fock space after timeevolution under the corresponding RN equations. Thequantum-classical correspondence is clear in this semi-classical regime and we can identify the hyperbolic X boundary, which is now dressed with characteristic inter-ference fringes across the fold lines. B. Triangular trimer
Let us now consider the δ QPM Hamiltonian Eq. (27)with K X = J , corresponding to a system of three sites ar-ranged in an equilateral triangle such that all three hop-ping amplitudes are equal. We follow the same procedureas the linear case, except that in order to bring the wave-function to the canonical form given in Eq. (31) the twophase variables must be scaled differently: φ C → / φ C and φ X → / φ X . We thereby obtain a modulus K = 2and find the following mapping between physical quanti-ties and abstract control parameters, α = β = √ (cid:126) ˜ U JN t − √ ζ = 2 / (cid:126) N Jt n X (39) η = 2 / √ (cid:126) N Jt δn (40) and, A ( t ) = 4 √ π (cid:126) ˜ U e − i π e i NJt (cid:126) . (41)Once again we find that α and β are not independentwhich this time leads to a spun cusp surrounding an axialcaustic line described by the generating functionΦ circ X = s + 2 s s + s + C ( s + s ) + C s + C s (42)and pictured in Fig. 11. The axial caustic line is notgeneric and is in fact unstable. It arises from the ro-tational symmetry we have assumed in the δ QPM with K X = J . Physically speaking, the combination of iso-metric hopping along with the simplified quantum phaseHamiltonian causes the system to not ‘feel’ the triangularsymmetry, and results in a perfectly circularly symmetricstructure in Fock space. From the parameters (38), wecan read off that the time at which the cusp point occursis t cusp = (cid:126) / (2 ˜ U JN ).Fig. 11 shows the dynamics of the wavefunction (31)with control parameters given in Eqs. (38)-(40) at t = (cid:126) / ( ˜ U JN ), twice the cusp time. Panel (a) shows the cir-cular caustic formed by the classical trajectories by mim-icking the initial state Eq. (28) with an even spread of ini-tial points in Fock space. Panel (b) shows the wavefunc-tion amplitude under the same conditions as (a), withclearly visible interference effects and a bright centralFock state amplitude | (cid:104) , | ψ ( t ) (cid:105) | corresponding to theaxial caustic. The canonical spun cusp caustic surface isshown in panel (c), with a two-dimensional circular caus-tic outlined in red. In Appendix B we give the controlparameters for the X wavefunction [Eq. (31)] for any in-termediate value of K X between the linear and triangulartrimer cases.5 FIG. 11. The spun cusp caustic generated by the triangularBH trimer ( K X = J ) within the δ QPM.
Panel (a):
Dis-tribution of points arising from classical trajectories of the δ QPM Hamiltonian Eq. (27) (see Appendix A for solutionsof Hamilton’s equations) at time t = (cid:126) / ( JN ˜ U ) (twice t cusp ),starting from an equal spread of number differences, and withphase differences φ C (0) = φ X (0) = 0. The circular causticfrom the K = 2 section of X is clearly visible. Panel (b):
Quantum wavefunction amplitude at the same time as panel(a), starting from an equal superposition of Fock states (28)and evolved using linearized RN equations. N = 150, and˜ U/ (cid:126) = 0 . Panel (c):
Caustic surface for the spun cusp,for positive C . The 2D circular caustic and the axial causticare highlighted in red. C. Gaussian Initial states
Using the same system parameters, let us revisit ourearlier choice of initial state as an equal superpositionof Fock states [Eq. (28)] which is the Fock space ana-logue of a plane wave. The primary goal of this choiceis to demonstrate the natural focusing effect of the BHdynamics. We claim that any sufficiently wide spread ofinitial Fock states will yield qualitatively similar resultswith only slightly altered coefficients. As an example,consider the ground state of the three-mode model (8)for U = 0, | ψ G (cid:105) = (cid:88) n + n + n = N (cid:114) N !3 N n ! n ! n ! | n , n , n (cid:105) (43)where the sum is over all Fock space occupation numberssuch that n + n + n = N . For more detail on the diag-onalization of this model, see, e.g., Ref. [194] and refer-ences therein. In the limit of large N , this coherent stateclosely resembles a Gaussian centered at δn = n X = 0. FIG. 12. A Gaussian initial state gives the same caus-tic as a ‘plane wave’.
Panel (a):
The amplitude of theground state of the noninteracting three-model model for K L = K X = K R ≡ J and N = 120 has a Gaussian form. Panel (b):
Time evolution of the ground state using thekicked Hamiltonian with ˜ U/ (cid:126) = 0 .
08, at time t = 2 (cid:126) / ( JN ˜ U ).Comparing with Fig. 11 we see a strong resemblance indicat-ing that initial states which differ significantly but have thesame general form (flat near the center of Fock space) willgive rise to qualitatively similar caustics. After changing variables, | ψ G (cid:105) ≈ √ πN (cid:88) δn , n X e − N (3 δn + n X ) | δn , n X (cid:105) . (44)Note that the coefficients of δn and n X are not the samesince n X ranges from − N to N and δn ranges from − N/ N/
3. The assumption of localized phase-statecontributions around φ = 0 still applies, since, (cid:104) φ X , φ C , Θ | ψ G (cid:105) ≈ √ πN ln 3 e iΘ N/ e − N ( φ C +3 φ X ) , (45)which for large N becomes narrowly peaked around φ X = φ C = 0, while widely spread in the number dif-ference values. Starting from this initial state, we usethe Floquet operator to propagate and the results areshown in Fig. 12 where we see that despite the fact thatthe Gaussian differs significantly from a plane wave westill obtain qualitatively the same caustic as in Fig. 11.Note that the result of applying the Floquet operator isto take U → U − i N , implying that for finite N thereis no longer a well-defined cusp point since the solution α = β = 0 cannot occur for a real value of t . The outercircular caustic surface does, however, remain and we seea diffraction pattern reminiscent of the K = 2 spun cusp.Furthermore, in the semiclassical regime N (cid:29) VII. BEYOND THE QUANTUM PHASEMODEL
In the QPM the effects of depletion of the modes [ac-counted for by the square root factors in the RN equa-tions Eq. (11)] are ignored. This is often a good approxi-mation in experimentally realizable superfluids, includingarrays of Josephson junctions [145, 146] and also some6 H δ QPM H δ MF H (cid:52) Quantum wavefunctionFIG. 13. Breaking circular symmetry by going beyond the QPM.
Panel (a):
The energy surface determining classicaltrajectories as provided by H δ QPM in Eq. (27) immediately following the δ -kick at t = 0. The initial values of the dynamicalvariables that generate this surface are φ X = φ C = 0 and an equal superposition of all number differences; the phase trajectoriesare bent by the δ -kick as explained in Appendix A. The overall scale is arbitrary and lighter colours represent higher energy.Note the circular symmetry as seen from the contour lines. Panel (b):
Same as panel (a) but now with the δ -kicked mean-fieldHamiltonian H δ MF given in Eq. (46) with all hopping amplitudes equal. The circular symmetry is broken by the square rootfactors and the resulting three-fold symmetric ‘valleys’ will act to pinch the caustic away from the corners of Fock space. Panel(c):
Same as panel (a) but now with H (cid:52) as given in Eq. (47). In this Hamiltonian the square root factors have been expandedto first order and this is enough break the circular symmetry and replace it with a three-fold symmetric energy surface. Panel(d):
Quantum wavefunction obtained using the full RN equations with δ -kicked interactions (quantum equivalent of H δ MF )and the same parameters as used in Fig. 11 (b). The chosen time is t = (cid:126) / ( JN ˜ U ) which is twice that of where the cusp pointappears. regimes of atomic BECs in optical lattices [147–152].However, in terms of caustics it can lead to some specialsituations such as the appearance of a K = 2 spun-cuspcaustic for the triangular trimer with K X = J , as dis-cussed in the previous section. This has circular symme-try which is of course non-generic and, indeed, the modu-lus K = 2 is excluded from the X catastrophe. The sym-metry is broken by any perturbation which is noncircular,such as choosing K X (cid:54) = J . Counter-intuitively, anotherway to break the symmetry is to go beyond the QPM byincluding the effects of the square root factors. For ourpurposes of illustrating generic many-body caustics, itsuffices to expand the square roots and keep just the firstorder corrections because this is enough to generate the X catastrophe (even when all the hopping amplitudesare equal), thereby providing an instructive example ofstructural stability, or lack thereof, for caustics that donot correspond to catastrophes. A. Triangular deformations: a path integralformulation
The mean-field Hamiltonian that retains the squareroot factors, has δ -kicked interactions and all-equal hop-ping is given by H δ MF = − J (cid:113) δn + N )( n X − δn + N ) cos( φ X − φ C ) − J (cid:113) δn + N )( − n X − δn + N ) cos( φ X + φ C ) − J (cid:113) ( δn − N ) − n X cos(2 φ X ) + δ ( t ) ˜ U (cid:2) δn + n X (cid:3) . (46) In Fig. 13(a) and (b) we compare the energy surfacesproduced by H δ QPM and H δ MF . It is evident from thecontour lines that H δ MF breaks the circular symmetryand replaces it with a triangular one. However, as ex-plained above, our interest is more in generic many-bodycaustics rather than specific models, so all we really needto do is perturb away from the QPM and hence we ex-pand the square roots in H δ MF and keep only the firstorder corrections: H (cid:52) = − J (cid:18) n X + δn + 2 N (cid:19) cos ( φ X − φ C ) − J (cid:18) − n X + δn + 2 N (cid:19) cos ( φ X + φ C ) − J (cid:18) N − δn (cid:19) cos (2 φ X ) + δ ( t ) ˜ U (cid:2) δn + n X (cid:3) ≡ N J Φ (cid:52) + δ ( t ) ˜ U (cid:2) δn + n X (cid:3) . (47)As can be seen from the energy surface in Fig. 13(c),this ‘triangular’ Hamiltonian retains the basic triangularsymmetry possessed by the full mean-field model. Thecone shape of the energy surface present in all three mod-els has an overall focusing effect on trajectories, but thetriangular shape of the latter two models channels thetrajectories away from the three corners of Fock spaceand effectively pinches the outer part of the caustic inthree spots. This can be seen in the quantum versionshown in panel (d) created using the RN equations underno approximations at t = 2 t cusp = (cid:126) / ( JN ˜ U ). The caus-tic evolves at roughly the same speed in both the classicaland quantum cases, and three bright spots appear bothbefore and after the cusp point and correspond to the7 FIG. 14. Stability of cusp formation in the n X = 0 planeagainst inclusion of beyond QPM effects. Panel (a):
Set ofclassical trajectories for the δ QPM when n X = 0 starting froma range of δn values and kicked with ˜ U/ (cid:126) = 0 .
02. The cusppoint here is
Jt/ (cid:126) ≈ . Panel (b):
Exact trajectories(without approximation to the square root factors) for thetriangular trimer with a kicked interaction term as given by H δ MF . The cusp structure remains intact and is only curvednear the edge of Fock space, showing that the δ QPM is a goodapproximation for this sub-catastrophe and that caustics havestructural stability against perturbations. three valleys on the triangular energy surfaces. Anothercrucial feature is that there is no bright patch at thecenter of Fock space, due to the destruction of the axialcaustic which is a feature peculiar to circular symmetry.In order to treat this problem analytically we needto reconsider our approach, since calculating an analyticform for the wavefunction as a single diffraction integralin the same manner as we did for the δ QPM is not possi-ble for the full Hamiltonian. More precisely, in the δ QPMthere is a separation between the phase and number dif-ference variables, like in a standard H = p / m + V ( x )type Hamiltonian, allowing us to integrate out the num-ber (momentum) variable using gaussian integrals leavinga diffraction integral purely in terms of the phase (posi-tion) variables. The presence of the square root factors effectively gives us a potential V ( x, p ) which depends onboth position and momentum variables. This difficultyis not reduced by expanding out the square root factorsand we instead resort to a phase space path integral-styleformulation for the wavefunction: ψ ( n X , δn , t ) = (cid:104) n X , δn | ˆ F| ψ (0) (cid:105) (48) ∝ (cid:90) D φ D n e i S [ n ( τ ) , φ ( τ )] (49)where we have introduced the bold vector notation φ =(Θ , φ X , φ C ) and n = ( N , n X , δn ). In the semiclassicalregime the dominant contributions to the path integralcome from the evaluation of the action S along classicalpaths which we denote by { ¯ n , ¯ φ } . However, the stan-dard WKB approximation blows up precisely at causticsbecause these are places where saddles of S coalesce. Fol-lowing Schulman [195] and Dangelmayr and Veit [196], aproper treatment of the problem shows that the leadingcontribution to the path integral close to caustics canbe factorized into a part involving the phase along theclassical path and part which is a diffraction integral ψ ( n , t ) ∝ e i S (¯ n , ¯ φ ,t ) (cid:90) (cid:90) d s d s e i λ Φ Q ( s ; C ) , (50)(in the above cited papers the diffraction integral is some-times referred to as a ‘generalized Airy function’). Forthe three-mode BH model, s is two-dimensional and de-pendent on φ , and the control parameters C depend onthe remaining parameters of the system, including n .The action that appears in Eq. 48 can be derived bybreaking up the time evolution operator into infinitesi-mal steps, i.e. applying the Trotter prescription to theoperator ˆ F with H (cid:52) . The details of this calculation arepresented in Appendix C where we find S [ n ( τ ) , φ ( τ )] = (cid:90) NJt (cid:126) d τ (cid:104) n · ˙ φ − Φ (cid:52) (cid:105) + (cid:126) U (cid:2) φ C (0) + 3 φ X (0) (cid:3) (51)= (cid:90) NJt/ (cid:126) d τ (cid:34) (cid:18) − N δn (cid:19) φ X + (cid:18)
118 + 124
N δn (cid:19) φ C + (cid:18)
13 + 14
N δn (cid:19) φ X φ C − N n X (cid:0) φ X φ C + φ X φ C (cid:1) − (cid:18) − N δn (cid:19) φ X − (cid:18)
23 + 12
N δn (cid:19) φ C + n X N φ X φ C + n · ˙ φ (cid:35) + (cid:126) U (cid:2) φ C (0) + 3 φ X (0) (cid:3) . Comparing to the standard relation L = p ˙ x − H , wesee that the generating function Φ plays the role of aLagrangian.The action given in Eq. (51) is not in a canonicalform for any unfolding of X , since the coefficients of thefourth order phase variables have not been appropriatelyscaled away yet, and the terms φ X φ C and φ X φ C are stillpresent. It is possible to remove these terms by an appro-priate change of variables, resulting in the introduction of the cubic unfolding terms φ X and φ X φ C , which leadto triangular symmetry and a stable caustic. Howeverwe shall not attempt this here since it involves the simul-taneous solution of five equations of quartic and cubicorder. We shall instead proceed by restricting ourselvesto the projection n X = 0. In Fig. 14 we plot the classicaltrajectories in the n X = 0 plane and compare betweenthe H δ QPM and the H δ MF cases. As expected, in (1+1)Dthe stable caustics are cusps and indeed the cusp point8 FIG. 15. Bordering (or abutment) diagram for the catastro-phe X , adapted from Nye [92] and obtained using theoremsof catastrophe projection. Each catastrophe of higher ordercontains (many, not necessarily all) catastrophes of lower or-der, shown by the direction of arrows. The subfamilies of X lead to two different sub-catastrophe sets. X containsthe singularities E and D , shown by the red dashed arrows. X does not contain E and will only contain E , shown bythe blue dotted arrow. occurs nearly simultaneously for H δ QPM and the H δ MF .This is because the square roots multiplying the cosinesdo not drastically affect the shape of the focusing surfacenear the centre of Fock space. Rather, the effects of thesquare roots only become significant near the edges ofFock space, a region in which the QPM approximationbecomes inaccurate. Indeed, we see from Fig. 14 that theeffect of the square root factors is to make fold lines ofthe full model, shown in panel (b), become curved nearthe boundaries. Thus, although the rotational symmetryof the caustic is removed by beyond QPM effects, otherfeatures of the caustic are robust against such changes, aresult that follows from structural stability.Nevertheless, the n X = 0 projection does allow us toaccount for some beyond-QPM effects and even obtainqualitative features such as the subfamily of X that re-sults from the breaking of the circular symmetry in H (cid:52) and H δ MF . In particular, after scaling phase variablesto put X in its canonical form, the modulus for the re-stricted problem for H (cid:52) becomes, K n X =0 ( δn ) = (8 N + 6 δn ) (cid:114) δn + 4 N (cid:114) N − δn . (52)This result warrants some explanation because the wave-function contains a range of the Fock space variables( δn , n X ), and hence there is no single well defined valueof the modulus. However, because all physical paths mustlie in the interval − N ≤ δn ≤ N , we have to first or-der, √ ≤ K ≤ K = 2 occurring only at δn = 0.This range of values for K indicates that the system sym-metry selects the X subfamily of X , and specificallyits X +9 variant, called compact by Callahan [192]. Aschematic plot, known as an ‘abutment’ or ‘bordering’diagram is given in Fig. 15 which summarizes the rela-tionships between the sub-catastrophes which can appearwithin X [92]. We see that the X subfamily does notcontain the umbilic catastrophes E or D given in TableII, both of which are instead members of the X sub- FIG. 16. Most singular sections of the diffraction patternsfor D − and D − showing how the foci of the two sub-familiesof X differ. Panel (a)
The focal plane of D − revealingan elliptic umbilic focus. Panel (b)
The focal plane of D − revealing another, but subtly different type of elliptic umbilic-like structure. In the case of D − the caustic in its canonicalform is threefold symmetric with the brightest ribs (tracedwith green dashed lines) meeting in straight lines at angles of2 π/
3. For D − the bright central ribs meet along curved linesobeying 27 C = 64 C . The presence of a D − focus can inprinciple be a diagnostic tool for determining the presence of X . family. Diagnosing the family of X experimentally viathe sub-catastrophes is in principle possible via a carefulanalysis of elliptic umbilic foci [Fig. 1(c), the last timeslice in Fig. 3(b), and Figs. 7(a) and (b) all show ellipticumbilic foci]. As illustrated in Fig. 16, elliptic umbilicfoci are three-fold symmetric about straight axes of sym-metry meeting at the origin, while for D − two of theselines become curved and the pattern is only two-fold sym-metric (see Nye [197] for more details).There is an important lesson to be learned from Eq.(52) and the fact there is no single value of the modulus(except locally, at the origin of Fock space). Catastro-phes have their origin in topology, and as such there isconsiderable flexibility as to their precise shape. In facteach catastrophe forms an equivalence class, where differ-ent specific realizations within each class are related bysmooth transformations (diffeomorphisms) of state vari-ables and control parameters (there is no smooth map-ping between different classes). The dynamics of a non-linear system such as the BH model results in a caustic ofthe X class which slowly varies in space-time but whichis not destroyed by nonlinearities. B. Unequal hopping amplitudes
So far our consideration of beyond-QPM dynamics hascentered on the effect of mode depletion accounted for bythe square root factors. However, one can also considerthe situation where we have already included the effectsof the mode depletion but then additionally break thesymmetry further by making the hopping amplitudes un-equal (a similar effect could be had by including a bias, (cid:15) i (cid:54) = 0). In Fig. 17 we see the effects of changing thehoppings so that they are slightly detuned from one an-other ( K R (cid:54) = K L (cid:54) = K X ). The elliptic umbilic caustic9 FIG. 17. Distorted elliptic umbilic due to unequal hop-ping strengths using the full three-mode Hamiltonian. K R =1 . K L , K X = 1 . K L , U/K L = 0 .
01 and N = 150, start-ing from an equal superposition of Fock states. Panel (a):
Distorted focus at K L t/ (cid:126) = 0 . Panel (b):
Unfolded at K L t/ (cid:126) = 0 . becomes distorted and asymmetric, but due to the struc-tural stability of the underlying catastrophe it remainsintact and recognizable. In particular, panel (a) showsthe effect on the focus, which can be shifted in space andtime but clearly retains its form. Likewise for panel (b)which shows a time slice somewhat after the focus. VIII. ROLE OF INTERACTIONS
In their original setting of natural optics (rainbowsetc.), caustics occur in a linear and hence integrable sys-tem. Similarly, the BH dimer is an integrable, albeit non-linear, system as long as energy is conserved. The trimer,by contrast, is not generally integrable. In our analyticcalculations in this paper we used δ -kicked interactionssuch that the time evolution is integrable. While thismade calculations possible, it does raise the question ofthe stability of caustics in the presence of constant inter-actions. It is therefore worth emphasizing again that allthe (numerical) examples of caustics shown in Section Vwere obtained with the interactions switched on through-out the time evolution. Moreover, interactions can some-times be necessary for catastrophes to fully manifest inboth the dimer and trimer. This is best understood withan example, as given in Fig. 18. In panel (a) U = 0, K L = K R = J , and K X = 0, for which the quantum re-vival time is Jt rev (cid:126) = π/ √ n X = 0plane, similar to Fig. 6(a) except that here the trajec-tories form isolated focal points. According to catastro-phe theory these are unstable in two dimensions and in-deed, in panel (b) where interactions are present the focalpoints are unfolded to cusps. For variety we have chosenattractive interactions U <
U >
0. The difference is that we findforward-opening cusps for
U >
U <
0. The effect is similar when startingfrom a highly focused state, such as a Fock state, wherea point focus will recur infinitely unless interactions are
FIG. 18. Instability of perfect focusing events to the in-troduction of interactions.
Panel (a):
Set of classical tra-jectories starting from an equal spread of Fock states, with K L = K R , U = K X = 0. Focusing events are isolated, mean-ing that all trajectories meet at a point. Panel (b):
Sameinitial state as panel (a), but now with weakly attractive in-teractions, K L = K R ≡ J and U/J = − .
03. A similarresult was found for repulsive interactions in panel (a) of Fig.6, except that the cusps open in the opposite direction forattractive interactions. introduced and structurally stable caustics form. Thenonlinearity introduced by interactions is therefore cru-cial to fully unfold the caustics (essentially by introducingdifferent periods for different amplitudes of excitation).Another example of dynamics with attractive interac-tions is shown in Fig. 19. The effect of negative U isto pull the caustic outwards towards the corners of Fockspace in comparison to the case of the repulsive U shownin Fig. 8. This means that the elliptic umbilic emergesmore clearly before the central focusing event, but alsoresults in fold ‘lips’ around the edges of each cusp point. FIG. 19. Formation of fold lines around cusps due to at-tractive interactions.
Panel (a):
Same as panel (c) of Fig.7, now with
U/J = − .
01. The elliptic umbilic diffractionpattern is now surrounded by bright fringes, which consist ofsmall fold lines.
Panels (b)-(c):
Same as panel (a) but at
Jt/ (cid:126) = 0 .
641 and
Jt/ (cid:126) = 0 . Panel (d):
Schematic of a section of the ellipticumbilic caustic, with ‘lips’ surrounding each cusp point. X with a negative modulus studied in[188] (see also Fig. 9). These lips continue to extend intolong fold lines as time progresses and ultimately intersectone another after the elliptic umbilic focus, reminiscentof the triple glass junction studied by Berry [189] andlater elaborated on by Nye [92].Finally, we note that our choice of parameters in thispaper has been guided empirically so as to make thecatastrophes as visible as possible. For example, the el-liptic and hyperbolic umbilic unfoldings shown in Sec-tion V were made using relatively weak interactions inthe range U/J = 0 . J to U/J = 0 . J . This regime iscalled the strong-tunneling (or Josephson [160]) regime(see Refs. [84, 157, 163, 165, 166, 174]), and exhibitsmoderate quantum revivals, allowing the catastrophesto be clearly identified because a single caustic stretchesacross a good fraction of Fock space while the nonzerointeractions prevent perfect focusing. Stronger interac-tions result in more powerful effective focusing poten-tials which, in the case of repulsive interactions, com-press the wavefunction into a small region in the centerof Fock space. The chaotic nature of the BH trimer inthis regime also will result in the eventual destructionof sharply visible diffraction patterns after longer times.In Fig. 20(a) we show that the canonical hyperbolic um-bilic catastrophe is still clearly visible at the strongercoupling U = J , although Fig. 20(b) indicates that atlonger times this wavefunction evolves into a highly in-tricate structure (which will eventually reveal the dis-creteness of Fock space as the wavefunction fringes reachsmall scales). The highly distorted wavefronts that arisein this situation would probably be best described usingthe statistical version of catastrophe theory developed inthe context of light passing through a turbulent atmo-sphere [91], where caustics manifest themselves as ex-treme amplitude events that occur more frequently thanexpected from random Gaussian fluctuations. This isthe freak/rogue wave paradigm recently explored in mi-crowave [97] and optical [101–109] experiments, althoughnow the discreteness of Fock space due to field quantiza-tion needs to be taken into account [66]. IX. DISCUSSION AND CONCLUSIONS
It is not clear a priori what a liquid drop lens and theBH trimer have in common, but the universality of wavecatastrophes means that similar striking caustic struc-tures can be seen in both. At a mathematical level thisis because, according to catastrophe theory, only certainshapes of singularity are structurally stable and henceoccur universally. More specifically, catastrophe theorydeals with singularities of gradient maps and in our phys-ical application the gradient map governing dynamics inboth systems is the minimization underlying Fermat’sprinciple (principle of least time/action). Singularities
FIG. 20. Catastrophe formation at high interaction strength.In this figure, we show exact (numerical) dynamics deep inthe chaotic regime with U = J in the linear BH trimer where K L = K R ≡ J , and K X = 0. Caustic formation followssimilarly to previously discussed cases with weaker interac-tions, however stronger interactions restrict the dynamics toa smaller region in Fock space. Panel (a):
Formation ofa fourfold diffraction pattern at
Jt/ (cid:126) = 0 .
13, similar to theone seen in Fig. 8(a), but its largest extent is smaller thanthe weaker interaction case.
Panel (b): At Jt/ (cid:126) = 0 . occur where multiple stationary points coalesce and thiscan only happen in a certain number of ways, dependingon the number of control parameters we have.The triangular nature of Fock space in the BH trimermeans that it behaves analogously to a liquid droplet lenswhose principal curvatures have been modified by a tri-angular perimeter. The optics of such lenses is knownto produce elliptic and hyperbolic umbilic catastropheswhich are three dimensional structures (see Fig. 1). Inthe BH trimer these catastrophes live in two dimensionalFock space plus time and could be generated by applyinga sudden quench to the interactions (or in fact the lat-tice) to generate a far-from-equilibrium state. The mean-field approximation corresponds to geometric rays thatare curved even in the absence of interactions. Many-1body states correspond to interfering superpositions ofmany rays, or more generally, a path integral as we sawin Section VII. Local to a caustic the action functionalin the path integral can be mapped to one of the catas-trophe generating functions Φ listed in Tables I and IIand in this way the wavefunction can be expressed interms of one of the canonical wave catastrophes. A caus-tic occurs where two or more rays coalesce (coalescenceof stationary points of the action). For example, theelliptic umbilic shown in Fig. 7 has a bright triangularborder inside of which four rays interfere to produce ahexagonal diffraction pattern. Two rays coalesce as wecross the boundaries (fold lines) except at the three cusppoints where three rays coalesce. If we choose a timeslice such that the elliptic umbilic is exactly in its focalplane ( C = 0) then all four rays coalesce, giving themost singular part of the caustic. The power of catastro-phe theory lies in the fact that not only can we predictthe geometric (classical) catastrophes that can occur ina system based on its dimension, but we can also pre-dict the qualitative form of the quantum wavefunctionbecause the wave patterns that dress catastrophes areuniversal.The closest analogue to the usual optical situation isan initial state which resembles a ‘plane wave’ in Fockspace. However, due to the the structural stability ofcatastrophes, the resulting caustic does not usually de-pend sensitively on the initial state. A moderately broadgaussian wavefunction (such as the ground state when thehopping amplitude is dominant, followed by a quench togenerate dynamics) with a flat region near the center ofFock space produces qualitatively similar results (SectionVI C). The caustics can be observed in an experiment bymeasuring the populations of each mode, e.g. by in situ absorption imaging of the atom density in the case of spa-tial wells [46]. Alternatively, caustics in the phase vari-ables can be measured by releasing the atoms from thewells and imaging after some time of flight to allow theatom clouds from each well to overlap and interfere [198].There is, however, an important difference between theclassical optical and quantum cases: in optics the com-plete diffraction pattern can be captured in a single shotlike in Fig. 1, whereas in the quantum case a single de-structive measurement will randomly collapse the wave-function. A single absorption image will yield just twonumbers: the two relative population differences betweenthe wells (assuming total number conservation). The ex-periment must be repeated many times under the sameconditions in order to build up a probability distributionfor the population differences thereby giving the modulusof each of the Fock space amplitudes. The whole schememust then be repeated for different evolution times in or-der to map out a three dimensional caustic. This wouldbe a laborious task, although less demanding than fullquantum tomography (see [66] for a brief discussion).Why have caustics not been discussed before in thecontext of many-body dynamics? The answer to this isprobably that many-body wavefunctions are high dimen- sional and hard to measure, as just mentioned. Althoughquantum gases can be probed in new and novel ways,information is traditionally obtained about many-bodysystems by measuring one- or two-point correlation func-tions. We note in this context that correlation functionsare in essence projections of wavefunctions and so maydisplay their own caustic structure; the possibility of asignature of caustics in correlations functions will be atopic for future investigation. Another possible futuredirection is to go over to the statistical version of catas-trophe theory mentioned in Section VIII, and that haspreviously been used to understand freak waves in hy-drodynamics and optics. This is particularly relevant inthe strong interaction regime where the wavefunction be-comes complicated and may appear random; caustics canbe detected through large amplitude fluctuations that oc-cur more frequently than expected from Gaussian statis-tics and whose moments obey certain characteristic scal-ing laws [91].We close by pointing out that the common origin ofuniversality in both equilibrium phase transitions andnonequilibrium caustics is the presence of underlying sin-gularities. It is a deep fact that close to singularities thebehaviour is highly organized and this leads to the re-markable level of caustic structure we see in the dynamicsof the BH dimer and trimer models. ACKNOWLEDGMENTS
We are grateful to the Natural Sciences and Engineer-ing Research Council of Canada (NSERC) for funding.
Appendix A: Classical Trajectories of the δ QPM
Using Hamilton’s equations of motion,˙ n = − (cid:126) ∂H∂φ C (A1)˙ n X = − (cid:126) ∂H∂φ X (A2)˙ φ C = 1 (cid:126) ∂H∂n (A3)˙ φ X = 1 (cid:126) ∂H∂n X (A4)The δ -kick allows classical trajectories to be calculateddirectly via integration, φ C ( t ) = φ C (0) + 3 U (cid:126) δn (0)Θ( t ) (A5) φ X ( t ) = φ X (0) + U (cid:126) n X (0)Θ( t ) , (A6)where we make use of the Heaviside function,Θ( t ) = (cid:40) t > t ≤ { δn (0) , n X (0) } , but all φ C (0) = φ X (0) = 0, then for t > δn ( t ) = δn (0) − JN t (cid:126) cos (cid:18) U (cid:126) n X (0) (cid:19) sin (cid:18) U (cid:126) δn (0) (cid:19) (A8) n X ( t ) = n X (0) − JN t (cid:126) (cid:34) cos (cid:18) U (cid:126) δn (0) (cid:19) sin (cid:18) U (cid:126) n X (0) (cid:19) (A9)+ 2 cos (cid:18) U (cid:126) n X (0) (cid:19) sin (cid:18) U (cid:126) n X (0) (cid:19)(cid:35) Note that the trajectories n X (0) = δn (0) = 0 have notime dependence and correspond exactly to the unstableaxial caustic in the K = 2 excluded family of X .In the case of the linear oligomer ( K X = 0), then δn ( t )remains the same, while, n X ( t ) = n X (0) − JN t (cid:126) cos (cid:18) U (cid:126) δn (0) (cid:19) sin (cid:18) U (cid:126) n X (0) (cid:19) . (A10) Appendix B: Dependence of δ QPM wavefunction on K X For general K X , the Hamiltonian (27) yields wavefunc-tions of the same form as in Eq. (31), but now withmodulus, K = 6 (cid:113) K X J + 1 (B1)and control parameters, α = √ (cid:0) (cid:126) − N JtU − N tU K X (cid:1) N JtU (cid:113) K X J + 1 (B2) β = √ (cid:0) (cid:126) − JN tU (cid:1)
JN tU (B3) ζ = 2 / √ (cid:126) N Jt (cid:0) K X J + 1 (cid:1) / n X (B4) η = 2 / √ (cid:126) N Jt δn (B5)and A ( t ) = (cid:18) K X /J + 1 (cid:19) / π (cid:126) √ U e − i π e i NJt (cid:126) . (B6) Appendix C: Derivation of the Path Integral
Starting with the Hamiltonian given in Eq. (47) whichgoes beyond the QPM by including the square root fac-tors to first order and hence includes effects due to the tri-angular boundaries of Fock space, we expand the cosineterms to fourth order in the phase difference coordinatesto give H (cid:52) ≈ − N J (cid:34)(cid:18) − N δn (cid:19) φ X + (cid:18)
118 + 124
N δn (cid:19) φ C + (cid:18)
13 + 14
N δn (cid:19) φ X φ C − N n X (cid:0) φ X φ C + φ X φ C (cid:1) (C1) − (cid:18) − N δn (cid:19) φ X − (cid:18)
23 + 12
N δn (cid:19) φ C + n X N φ X φ C (cid:35) + δ ( t ) ˜ U (cid:2) δn + n X (cid:3) ≡ N J Φ (cid:52) + δ ( t ) ˜ U (cid:2) δn + n X (cid:3) . (C2)Notice that relative to the expansion of H δ QPM , all of the circularly symmetric terms get a perturbation proportionalto δn and we also pick up some non-circularly symmetric terms proportional to n X . Applying the Floquet operatorincorporating this Hamiltonian to an initial state comprising of an equal superposition of all Fock states, and projecting3onto the Fock basis, we obtain the wavefunction ψ ( n X , δn , t ) = (cid:88) n (cid:48) X ,δn (cid:48) e − i U [ δn (cid:48) + n (cid:48) X ] (cid:104) n X , δn | e i NJt (cid:126) Φ (cid:52) (ˆ n X ,δ ˆ n , ˆ φ C , ˆ φ X ) (cid:122) (cid:125)(cid:124) (cid:123) ≈ (cid:80) Θ ,φX,φC | Θ ,φ X ,φ C (cid:105)(cid:104) Θ ,φ X ,φ C | | n (cid:48) X , δn (cid:48) (cid:105) (C3)= (cid:90) dΘd φ X d φ C (cid:104) n X , δn | e i NJt (cid:126) Φ (cid:52) (ˆ n X ,δ ˆ n , ˆ φ C , ˆ φ X ) t | Θ , φ X , φ C (cid:105) e − i Θ N (cid:88) n (cid:48) X ,δn (cid:48) e − i U (cid:126) [ δn (cid:48) + n (cid:48) X ]e − i [ n (cid:48) X φ X + δn (cid:48) φ C ] (cid:124) (cid:123)(cid:122) (cid:125) ≈ π (cid:126) i √ U exp [ i (cid:126) U ( φ C +3 φ X )] (C4)where in the first line we have indicated where a resolution of the identity in terms of phase states should be inserted,and in second line we have turned the resulting double sums over phase variables into integrals as well as indicatingthat the double sum over the primed number variables can be approximated by gaussian integrals. Let us now focuson the matrix elements of e i NJt (cid:126) ˆΦ (cid:52) , (cid:104) n X , δn | e i NJt (cid:126) Φ (cid:52) (ˆ n X ,δ ˆ n , ˆ φ C , ˆ φ X ) | Θ , φ X , φ C (cid:105) ≈ (cid:104) n X , δn | (cid:18) + i N JtM (cid:126) ˆΦ (cid:52) (cid:19) M | Θ , φ X , φ C (cid:105) (C5)= (cid:104) n X , δn | (cid:18) + i N JtM (cid:126) ˆΦ (cid:52) (cid:19) (cid:18) + i N JtM (cid:126) ˆΦ (cid:52) (cid:19) ... (cid:18) + i N JtM (cid:126) ˆΦ (cid:52) (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) M times | Θ , φ X , φ C (cid:105) , (C6)where M is an integer giving the number of infinitesimal time steps t/M into which the propagation is decomposed.We shall assume that M (cid:29) N (cid:29)
1. Switching to the bold vector notation φ = (Θ , φ X , φ C ) and n = ( N , n X , δn ) forbrevity, we insert resolutions of the identity = (cid:90) d φ | φ (cid:105) (cid:104) φ | (C7)between each set of parentheses (cid:104) n X , δn | e i NJt (cid:126) Φ (cid:52) (ˆ n X ,δ ˆ n , ˆ φ C , ˆ φ X ) | Θ , φ X , φ C (cid:105) ≈ (cid:90) d φ M − . . . d φ (cid:26) (cid:104) n X , δn | (cid:18) + i N JtM (cid:126) ˆΦ (cid:52) (cid:19) | φ M − (cid:105)× (cid:104) φ M − | (cid:18) + i N JtM (cid:126) ˆΦ (cid:52) (cid:19) | φ M − (cid:105) × . . . × (cid:104) φ | (cid:18) + i N JtM (cid:126) ˆΦ (cid:52) (cid:19) | φ (cid:105) (cid:27) , (C8)where we have used | φ (cid:105) to denote | Θ , φ X , φ C (cid:105) . In order to evaluate each matrix element in this product, we notethat ˆΦ (cid:52) has the formˆΦ (cid:52) = (cid:18) − N ˆ δn (cid:19) ˆ φ X + (cid:18)
118 + 124 N ˆ δn (cid:19) ˆ φ C + (cid:18)
13 + 14 N ˆ δn (cid:19) ˆ φ X ˆ φ C − N ˆ n X (cid:16) ˆ φ X ˆ φ C + ˆ φ X ˆ φ C (cid:17) + ... (C9)where the operator ordering is assumed to be the same as the classical expression given in Eq. (C1) such that allthe number operators lie to the left of the phase operators. We therefore insert resolutions of the identity over thenumber states = (cid:82) d n | n (cid:105) (cid:104) n | inside the matrix elements so that the number operators can act to the left and thephase operators to the right (cid:104) φ j +1 | (cid:18) + i N JtM (cid:126) ˆΦ (cid:52) (cid:19) | φ j (cid:105) −→ (cid:104) φ j +1 | (cid:90) d n j | n j (cid:105) (cid:104) n j | (cid:18) + i N JtM (cid:126) ˆΦ (cid:52) (cid:19) | φ j (cid:105) (C10) ≈ (cid:90) d n j exp [i n j · ( φ j +1 − φ j )] exp (cid:20) i N JtM (cid:126) Φ (cid:52) ( n j , φ j ) (cid:21) (C11)where to obtain the second line we have used the relation (cid:104) φ | n (cid:105) = exp[i n · φ ] twice and Φ (cid:52) ( n j , φ j ) is now a functionof ordinary variables rather than operators. Thus, all the matrix elements apart from the most lefthand one in Eq.(C8) contribute a phase factor exp (cid:2) i NJtM (cid:126) Φ (cid:52) ( n j , φ j ) + i n j · ( φ j +1 − φ j ) (cid:3) , and the wavefunction becomes ψ ( n X , δn , t ) = (cid:90) d φ M − (cid:90) M − (cid:89) j =0 d φ j d n j B ( n X , δn , φ M − ) (cid:124) (cid:123)(cid:122) (cid:125) boundary term exp i M − (cid:88) j =0 (cid:26) N JtM (cid:126) Φ (cid:52) ( n j , φ j ) + n j · ( φ j +1 − φ j ) (cid:27) × (cid:18) π (cid:126) i √ U exp (cid:20) i (cid:126) U (cid:0) { φ C } + 3 { φ X } (cid:1)(cid:21)(cid:19) . (C12)4The boundary term comes from the most lefthand matrix element in Eq. (C8) and is also a pure phase factor B ( n X , δn , φ M − ) = exp (cid:20) i N JtM (cid:126) Φ (cid:52) [ n X , δn , φ M − ] − n X { φ X } M − − δn { φ C } M − (cid:21) . (C13)in which the notation { φ X } M − and { φ C } M − is used for the individual components of φ M − and likewise { φ X } and { φ C } is used for the individual components of φ .A condition that the above expansion of the propagator into M terms is an accurate approximation is that N Jt/ ( (cid:126) M ) ≡ ε is a small quantity. Assuming this to be the case we can write B ( n X , δn , φ M − ) + ε M − (cid:88) j =0 (cid:26) Φ (cid:52) [ n j , φ j ] + n j · φ j +1 − φ j ε (cid:27) ∼ (cid:90) NJt/ (cid:126) d τ (cid:110) Φ (cid:52) [ n ( τ ) , φ ( τ )] + n ( τ ) · ˙ φ ( τ ) (cid:111) . (C14)This becomes exact when ε →
0, or equivalently M → ∞ so that we can express the wavefunction as a path integral ψ ( n X , δn , t ) = 4 π (cid:126) i √ U (cid:90) D φ D n exp (cid:34) i (cid:90) NJt/ (cid:126) d τ { Φ (cid:52) [ n ( τ ) , φ ( τ )] + n ( τ ) · ˙ φ ( τ ) } (cid:35) exp (cid:20) i (cid:126) U (cid:8) φ C (0) + 3 φ X (0) (cid:9)(cid:21) (C15)where D φ D n = lim M →∞ d φ M − M − (cid:89) j =0 d φ j d n j = lim M →∞ M − (cid:89) j =0 d φ j d n j δ [ n M − − ( n X , δn )] (C16) [1] I. Bloch, J. Dalibard, and W. Zwerger, Many-bodyphysics with ultracold gases, Rev. Mod. Phys. , 885(2008).[2] J. M. Deutsch, Quantum statistical mechanics in aclosed system, Phys. Rev. A.
888 (1994).[4] J. Eisert, M. Friesdorf, and C. Gogolin, Quantum many-body systems out of equilibrium, Nat. Phys. 11, 124(2015).[5] A. Polkovnikov, K. Sengupta, A. Silva, and M. Ven-galattore, Nonequilibrium dynamics of closed interact-ing quantum systems, Rev. Mod. Phys. , 863 (2011).[6] E. Nicklas, M. Karl, M. H¨ofer, A. Johnson, W. Mues-sel, H. Strobel, J. Tomkoviˇc, T. Gasenzer, and M. K.Oberthaler, Observation of Scaling in the Dynamics ofa Strongly Quenched Quantum Gas, Phys. Rev. Lett. , 245301 (2015).[7] S. Erne, R. B¨ucker, T. Gasenzer, J. Berges, J.Schmiedmayer, Universal dynamics in an isolated one-dimensional Bose gas far from equilibrium, Nature ,225 (2018).[8] M. Heyl, Dynamical quantum phase transitions: a re-view, Rep. Prog. Phys. , 054001 (2018).[9] M. P. A. Fisher, P. B. Weichman, G. Grinstein, D. S.Fisher, Boson localization and the superfluid-insulatortransition, Phys. Rev. B
546 (1989).[10] S. Sachev, Quantum Phase Transitions (CambridgeUniversity Press, New York, 2001).[11] M. Greiner, O. Mandel, T. Esslinger, T. W. H¨ansch, andI. Bloch, Quantum phase transition from a superfluidto a Mott insulator in a gas of ultracold atoms, Nature , 39 (2002).[12] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, andP. Zoller, Cold Bosonic Atoms in Optical Lattices, Phys.Rev. Lett. , 3108 (1998).[13] E. Altman and A. Auerbach, Oscillating Superfluidity ofBosons in Optical Lattices, Phys. Rev. Lett. , 250404(2002).[14] L. Isella and J. Ruostekoski, Nonadiabatic dynamics ofa Bose-Einstein condensate in an optical lattice, Phys.Rev. A , 011601(R) (2005).[15] R. Sch¨utzhold, M. Uhlmann, Y. Xu, and U. R. Fis-cher, Sweeping from the Superfluid to the Mott Phase inthe Bose-Hubbard Model, Phys. Rev. Lett. , 200601(2006).[16] C. Kollath, A. M. L¨auchli, and E. Altman, Quench Dy-namics and Nonequilibrium Phase Diagram of the Bose-Hubbard Model, Phys. Rev. Lett. 98, 180601 (2007)[17] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Re-laxation in a completely integrable many-body quantumsystem: An ab initio study of the dynamics of the highlyexcited states of 1d lattice hard-core bosons, Phys. Rev.Lett. , 050405 (2007).[18] M. Rigol, V. Dunjko, and M. Olshanii, Thermalizationand its mechanism for generic isolated quantum sys-tems, Nature , 854 (2008).[19] J. Dziarmaga, M. Tylutki, and W. H. Zurek, Quenchfrom Mott insulator to superfluid, Phys. Rev. B ,144521 (2012).[20] P. Barmettler, D. Poletti, M. Cheneau, and C. Kollath,Propagation front of correlations in an interacting Bosegas, Phys. Rev. A , 053625 (2012).[21] A. J. Daley, H. Pichler, J. Schachenmayer, and P. Zoller, Measuring Entanglement Growth in Quench Dynamicsof Bosons in an Optical Lattice, Phys. Rev. Lett. ,020505 (2012)[22] M. (cid:32)Lacki and J. Zakrzewski, Fast Dynamics for Atoms inOptical Lattices, Phys. Rev. Lett. , 065301 (2013).[23] G. Kordas, D. Witthaut, and S. Wimberger, Non-equilibrium dynamics in dissipative Bose-Hubbardchains, Ann. Phys. (Berlin) , 619 (2015).[24] F. Vicentini, F. Minganti, R. Rota, G. Orso, and C.Ciuti, Critical slowing down in driven-dissipative Bose-Hubbard lattices, Phys. Rev. A , 013853 (2018).[25] J. G. Cosme, Hierarchical relaxation dynamics in atilted two-band Bose-Hubbard model, Phys. Rev. A 97,043610 (2018).[26] M. R. C. Fitzpatrick and M. P. Kennett, Light-cone-like spreading of single-particle correlations in the Bose-Hubbard model after a quantum quench in the strong-coupling regime, Phys. Rev. A , 053618 (2018).[27] K. Nagao, M. Kunimi, Y. Takasu, Y. Takahashi, and I.Danshita, Semiclassical quench dynamics of Bose gasesin optical lattices, Phys. Rev. A , 023622, (2019).[28] L. Fallani, L. De Sarlo, J. E. Lye, M. Modugno, R. Saers,C. Fort, and M. Inguscio, Observation of dynamical in-stability for a Bose-Einstein condensate in a moving 1Doptical lattice, Phys. Rev. Lett. , 140406 (2004).[29] A. K. Tuchman, C. Orzel, A. Polkovnikov, and M. A.Kasevich, Nonequilibrium coherence dynamics of a softboson lattice, Phys. Rev. A , 051601(R) (2006).[30] S. Will, T. Best, U. Schneider, L. Hackerm¨uller, D.-S. L¨uhmann and I. Bloch, Time-resolved observationof coherent multi-body interactions in quantum phaserevivals, Nature , 197 (2010).[31] S. Trotzky, Y-A. Chen, A. Flesch, I. P. McCulloch, U.Schollw¨ock, J. Eisert and I. Bloch, Probing the relax-ation towards equilibrium in an isolated strongly corre-lated one-dimensional Bose gas, Nature Physics , 325(2012).[32] M. Cheneau, P. Barmettler, D. Poletti, M. Endres, P.Schauß, T. Fukuhara, C. Gross, I. Bloch, C. Kollath,and S. Kuhr, Light-cone-like spreading of correlations ina quantum many-body system, Nature , 484 (2012).[33] T. Fukuhara, P. Schauß, M. Endres, S. Hild, M. Che-neau, I. Bloch and C. Gross, Microscopic observation ofmagnon bound states and their dynamics, Nature ,76 (2013).[34] F. Meinert, M. J. Mark, E. Kirilov, K. Lauber, P. Wein-mann, A. J. Daley, and H.-C. N¨agerl, Quantum Quenchin an Atomic One-Dimensional Ising Chain, Phys. Rev.Lett. , 053003 (2013).[35] S. Braun, M. Friesdorf, S. S. Hodgman, M. Schreiber,J. P. Ronzheimer, A. Riera, M. del Rey, I. Bloch, J.Eisert, and U. Schneider, Emergence of coherence andthe dynamics of quantum phase transitions, PNAS ,3641 (2015).[36] J. G. Cosme, C. Georges, A. Hemmerich, and L.Mathey, Dynamical control of order in a cavity-BECsystem, Phys. Rev. Lett. , 153001 (2018).[37] T. Boulier, J. Maslek, M. Bukov, C. Bracamontes, E.Magnan, S. Lellouch, E. Demler, N. Goldman, andJ. V. Porto, Parametric Heating in a 2D PeriodicallyDriven Bosonic System: Beyond the Weakly Interact-ing Regime, Phys. Rev. X , 011047 (2019).[38] G. J. Milburn, J. Corney, E. M. Wright, and D. F. Walls,Quantum dynamics of an atomic Bose-Einstein conden- sate in a double-well potential. Phys. Rev. A , 4318(1997).[39] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R.Shenoy, Quantum coherent atomic tunneling betweentwo trapped Bose-Einstein condensates. Phys. Rev.Lett. , 4950 (1997).[40] A. Vardi and J. R. Anglin, Bose-Einstein Condensatesbeyond Mean Field Theory: Quantum Backreaction asDecoherence, Phys. Rev. Lett. , 568 (2001).[41] G.-S. Paraoanu, S. Kohler, F. Sols, and A. J. Leggett,The Josephson plasmon as a Bogoliubov quasiparticle,J. Phys. B: At. Mol. Opt. Phys. , 180402 (2001).[43] E. M. Graefe and H. J. Korsch, Semiclassical quantiza-tion of an N-particle Bose-Hubbard model, Phys. Rev.A , 032116 (2007).[44] M. Chuchem, K. Smith-Mannschott, M. Hiller, T. Kot-tos, A. Vardi, and D. Cohen, Quantum dynamics in thebosonic Josephson junction, Phys. Rev. A , 053617(2010).[45] H. Veksler and S. Fishman, Semiclassical analysis ofBose-Hubbard dynamics, New J. Phys. , 053030(2015).[46] M. Albiez, R. Gati, J. F¨olling, S. Hunsmann, M. Cris-tiani, and M. K. Oberthaler, Direct observation of tun-neling and nonlinear self-trapping in a single bosonicJosephson junction, Phys. Rev. Lett. , 010402 (2005).[47] S. Levy, E. Lahoud, I. Shomroni, and J. Steinhauer,The a.c. and d.c. Josephson effects in a Bose-Einsteincondensate, Nature , 579 (2007).[48] L. J. LeBlanc, A. Bardon, J. McKeever, M. Extavour,D. Jervis, J. Thywissen, F. Piazza, and A. Smerzi, Dy-namics of a Tunable Superfluid Junction, Phys. Rev.Lett. 106, 025302 (2011).[49] S. F¨olling, S. Trotzky, P. Cheinet, M. Feld, R. Saers, A.Widera, T. M¨uller, and I. Bloch, Direct observation ofsecond-order atom tunnelling, Nature , 1029 (2007).[50] C. Ryu, P. W. Blackburn, A. A. Blinova, and M. G.Boshier, Experimental realization of Josephson junc-tions for an atom SQUID, Phys. Rev. Lett. , 205301(2013).[51] A. Trenkwalder, G. Spagnolli, G. Semeghini, S. Coop,M. Landini, P. Castilho, L. Pezz`e, G. Modugno, M. In-guscio, A. Smerzi, and M. Fattori, Quantum phase tran-sitions with parity-symmetry breaking and hysteresis.Nat. Phys. , 826 (2016).[52] T. Zibold, E. Nicklas, C. Gross, and M. K. Oberthaler,Classicial bifurcation at the Transition from Rabi toJosephson dynamics. Phys. Rev. Lett. , 204101(2010).[53] A. Das, K. Sengupta, D. Sen, and B. K. Chakrabarti,Infinite-range Ising ferromagnet in a time-dependenttransverse magnetic field: quench and ac dynamics nearthe quantum critical point. Phys. Rev. B , 144423(2006).[54] J. W. Britton, B. C. Sawyer, A. C. Keith, C.-C. J.Wang, J. K. Freericks, H. Uys, M. J. Biercuk, and J.J. Bollinger, Engineered two-dimensional Ising interac-tions in a trapped-ion quantum simulator with hundredsof spins, Nature , 489 (2012).[55] P. Jurcevic, B. P. Lanyon, P. Hauke, C. Hempel, P.Zoller, R. Blatt, and C. F. Roos, Quasiparticle engineer- ing and entanglement propagation in a quantum many-body system, Nature , 202 (2014).[56] P. Richerme, Z.-X. Gong, A. Lee, C. Senko, J. Smith,M. Foss-Feig, S. Michalakis, A. V. Gorshkov, and C.Monroe, Non-local propagation of correlations in quan-tum systems with long-range interactions, Nature ,198 (2014).[57] P. Ribeiro, J. Vidal, and R. Mosseri, Exact spectrum ofthe Lipkin-Meshkov-Glick model in the thermodynamiclimit and finite-size corrections, Phys. Rev. E 78, 021106(2008).[58] L. F. Santos, M. T´avora, and F. P´erez-Bernal, Excited-state quantum phase transitions in many-body systemswith infinite-range interaction: Localization, dynamics,and bifurcation, Phys. Rev. A , 012113 (2016).[59] A. Reinhard, J.-F. Riou, L. A. Zundel, D. S. Weiss, S.Li, A. M. Rey, and R. Hipolito, Self-trapping in an arrayof coupled 1D Bose gases, Phys. Rev. Lett. , 033001(2013).[60] M. Abbarchi, A. Amo, V. G. Sala, D. D. Solnyshkov, H.Flayac, L. Ferrier, I. Sagnes, E. Galopin, A. Lemaitre,G. Malpuech, and J. Bloch, Macroscopic quantumself-trapping and Josephson oscillations of exciton-polaritons, Nat. Phys. , 275 (2013).[61] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P.Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C.Monroe, Observation of a many-body dynamical phasetransition with a 53-qubit quantum simulator, Nature , 601 (2017).[62] S. Smale, P. He, B. A. Olsen, K. G. Jackson, H. Sharum,S. Trotzky, J. Marino, A. M. Rey, and J. H. Thywissen,Observation of a transition between dynamical phases ina quantum degenerate Fermi gas, Sci Adv , eaax1568(2019).[63] J. A. Muniz, D. Barberena, R. J. Lewis-Swan, D. J.Young, J. R. K. Cline, A. M. Rey and J. K. Thompson,Exploring dynamical phase transitions with cold atomsin an optical cavity, Nature , 607 (2020).[64] D. H. J. O’Dell, Quantum catastrophes and ergodicityin the dynamics of bosonic Josephson junctions, Phys.Rev. Lett. , 150406 (2012).[65] J. Mumford, W. Kirkby, and D. H. J. O’Dell, Catastro-phes in non-equilibrium many-particle wave functions:universality and critical scaling, J. Phys. B: At. Mol.Opt. Phys. , 044005 (2017).[66] J. Mumford, E. Turner, D. W. L. Sprung, and D. H. J.O’Dell, Quantum Spin Dynamics in Fock Space Follow-ing Quenches: Caustics and Vortices, Phys. Rev. Lett. , 170402 (2019).[67] C. K. Law, H. Pu, and N. P. Bigelow, Quantum SpinsMixing in Spinor Bose-Einstein Condensates, Phys.Rev. Lett. , 5257 (1998).[68] W. Zhang, D. L. Zhou, M.-S. Chang, M. S. Chapman,and L. You, Coherent spin mixing dynamics in a spin-1atomic condensate, Phys. Rev. A , 013602 (2005).[69] D. M. Stamper-Kurn and M. Ueda, Spinor Bose gases:Symmetries, magnetism, and quantum dynamics, Rev.Mod. Phys. , 1191 (2013).[70] J. Stenger, S. Inouye, D. M. Stamper-Kurn, H.-J. Mies-ner, A. P. Chikkatur and W. Ketterle, Spin domainsin ground-state Bose-Einstein condensates, Nature ,345 (1998).[71] A. T. Black, E. Gomez, L. D. Turner, S. Jung, and P.D. Lett, Spinor dynamics in an antiferromagnetic spin-1 condensate, Phys. Rev. Lett. , 070403 (2007).[72] L. Sadler, J. Higbie, S. Leslie, M. Vengalattore, andD. Stamper-Kurn, Spontaneous symmetry breaking ina quenched ferromagnetic spinor Bose-Einstein conden-sate, Nature (London) , 312 (2006).[73] C. S. Gerving, T. M. Hoang, B. J. Land, M. Anquez,C. D. Hamley, and M. S. Chapman, Non-equilibriumdynamics of an unstable quantum pendulum exploredin a spin-1 Bose-Einstein condensate, Nature Commun. , 1169 (2012).[74] D. Linnemann, H. Strobel, W. Muessel, J. Schulz, R. J.Lewis-Swan, K. V. Kheruntsyan, and M. K. Oberthaler,Quantum-Enhanced Sensing Based on Time Reversalof Nonlinear Dynamics, Phys. Rev. Lett. , 013001(2016).[75] K. Lange, J. Peise, B. L¨ucke, I. Kruse, G. Vitagliano,I. Apellaniz, M. Kleinmann, G. T´oth, and C. Klempt,Entanglement between two spatially separated atomicmodes, Science , 416 (2018).[76] P. Kunkel, M. Pr¨ufer, H. Strobel, D. Linnemann,A. Fr¨olian, T. Gasenzer, M. G¨arttner, and M. K.Oberthaler, Spatially distributed multipartite entangle-ment enables EPR steering of atomic clouds, Science , 413 (2018).[77] P. Kunkel, M. Pr¨ufer, S. Lannig, R. Rosa-Medina, A.Bonnin, M. G¨arttner, H. Strobel, and M. K. Oberthaler,Simultaneous Readout of Noncommuting CollectiveSpin Observables beyond the Standard Quantum Limit,Phys. Rev. Lett. , 063603 (2019).[78] M. Rautenberg and Martin G¨arttner, Classical andquantum chaos in a three-mode bosonic system, Phys.Rev. A , 053604 (2020).[79] R. Franzosi and V. Penna, Chaotic behavior, collec-tive modes, and self-trapping in the dynamics of threecoupled Bose-Einstein condensates, Phys. Rev. E ,046227 (2003).[80] B. Liu, L.-B. Fu, S.-P. Yang, and J. Liu, Josephson oscil-lation and transition to self-trapping for Bose-Einsteincondensates in a triple-well trap, Phys. Rev. A ,033601 (2007).[81] P. Buonsante, V. Penna, and A. Vezzani, Quantum sig-natures of the self-trapping transition in attractive lat-tice bosons, Phys. Rev. A 82, 043615 (2010).[82] M. A. Garcia-March, S. van Frank, M. Bonneau, J.Schmiedmayer, M. Lewenstein, and L. F. Santos, Relax-ation, chaos, and thermalization in a three-mode modelof a Bose-Einstein condensate, New J. Phys. , 113039(2018).[83] Q. Thommen, J. C. Garreau, and V. Zehnl´e, Classicalchaos with Bose-Einstein condensates in tilted opticallattices, Phys. Rev. Lett. , 210405 (2003).[84] A. R. Kolovsky, Semiclassical quantization of the Bo-goliubov spectrum, Phys. Rev. Lett. , 020401 (2007).[85] M. Hiller, T. Kottos, and T. Geisel, Wave-packet dy-namics in energy space of a chaotic trimeric Bose-Hubbard system, Phys. Rev. A , 023621 (2009).[86] T. F. Viscondi and K. Furuya, Dynamics of a Bose-Einstein condensate in a symmetric triple-well trap, J.Phys. A: Math. Theor. , 053601 (2011).[88] D. Peter, K. Pawlowski, T. Pfau, and K. Rza˙zewski, Mean-field description of dipolar bosons in triple-wellpotentials, J. Phys. B: At. Mol. Opt. Phys. (2012).[89] L. Dell’Anna, G. Mazzarella, V. Penna, and L. Salas-nich, Entanglement entropy and macroscopic quantumstates with dipolar bosons in a triple-well potential,Phys. Rev. A , 053620 (2013).[90] S. Baier, M. J. Mark, D. Petter, K. Aikawa, L. Chomaz,Z. Cai, M. Baranov, P. Zoller, and F. Ferlaino, ExtendedBose-Hubbard models with ultracold magnetic atoms,Science , 201 (2016).[91] M. V. Berry, Focusing and twinkling: critical exponentsfrom catastrophes in non-Gaussian random short waves.J. Phys. A: Math. Gen. , 2061 (1977).[92] J. F. Nye, Natural Focusing and Fine Structure of Light (Institute of Physics, Philadelphia, 1999).[93] Y. A. Kravtsov and Y. I. Orlov,
Caustics, Catastro-phes and Wave Fields
Second Edition (Springer-Verlag,Berlin, 1993).[94] Lord Kelvin, Deep water ship-waves, Phil. Mag. 9, 733(1905).[95] M. V. Berry, Minimal analytical model for undular tidalbore profile; quantum and Hawking effect analogies.New J. Phys. , 053066 (2018).[96] B. S. White and B. Fornberg, On the chance of freakwaves at sea, J. Fluid Mechanics , 113 (1998).[97] R. H¨ohmann, U. Kuhl, H.-J. St¨ockmann, L. Kaplan,and E. J. Heller, Freak waves in the linear regime: amicrowave study. Phys. Rev. Lett. , 093901 (2010).[98] M. Onoratoab, S. Residoric, U. Bortolozzoc, A. Monti-nad, and F. T. Arecchi, Rogue waves and their gener-ating mechanisms in different physical contexts, Phys.Rep. , 47 (2013).[99] T. A. A. Adcock and P. H. Taylor, The physics ofanomalous (’rogue’) ocean waves, Rep. Prog. Phys. ,105901 (2014).[100] H. Degueldre, J. J. Metzger, T. Geisel and R. Fleis-chmann, Random focusing of tsunami waves, Nat. Phys. , 259 (2016).[101] D. R. Solli, C. Ropers, P. Koonath and B. Jalali, OpticalRogue Waves, Nature , 1054 (2007).[102] F. T. Arecchi, U. Bortolozzo, A. Montina, and S. Resi-dori, Granularity and Inhomogeneity Are the Joint Gen-erators of Optical Rogue Waves, Phys. Rev. Lett. ,153901 (2011).[103] N. Akhmediev, J. M. Dudley, D. R. Solli, and S. K.Turitsyn, Recent progress in investigating optical roguewaves, J. Opt. , 060201 (2013).[104] N. Marsal, V. Caullet, D. Wolfersberger, and M.Sciamanna, Spatial rogue waves in a photorefractivepattern-forming system, Opt. Lett. , 3690 (2014).[105] A. Mathis, L. Froehly, S. Toenger, F. Dias, G. Genty,and J. M. Dudley, Caustics and Rogue Waves in anOptical Sea, Sci. Rep. , 12822 (2015).[106] Spatial Rogue Waves in Photorefractive Ferroelectrics,D. Pierangeli, F. Di Mei, C. Conti, A.?J. Agranat, andE. DelRe, Phys. Rev. Lett. , 093901 (2015).[107] M. Mattheakis, J. Pitsios, G.P. Tsironis, S. Tzortzakis,Extreme events in complex linear and nonlinear pho-tonic media, Chaos, Solitons and Fractals , 73 (2016).[108] A. Safari, R. Fickler, M. J. Padgett, R. W. Boyd, Gen-eration of Caustics and Spatial Rogue Waves from Non-linear Instability, Phys. Rev. Lett. , 203901 (2017).[109] A. Zannotti, Caustic Light in Nonlinear Photonic Media (Springer, Cham, Switzerland, 2020). [110] V. I. Arnold, S. F. Shandarin, and Ya. B. Zeldovich, Thelarge scale structure of the Universe I. General proper-ties. One- and two-dimensional models, Geophys. As-trophys. Fluid Dyn. , 111 (1982).[111] J. Feldbrugge, R. van de Weygaert, J. Hidding, and J.Feldbrugge, Caustic skeleton & cosmic web, J. Cosmol.Astropart. Phys. 05 (2018) 027.[112] M. V. Berry and C. Upstill, Catastrophe optics: mor-phologies of caustics and their diffraction patterns,Prog. Opt. XVIII , 257 (1980).[113] R. Thom,
Structural Stability and Morphogenesis (Ben-jamin, Reading MA, 1975).[114] V. I. Arnol’d, Critical points of smooth functions andtheir normal forms
Russ. Math. Survs. , 1 (1975).[115] E. C. Zeeman, Catastrophe Theory: Selected Papers1972-1977 (Addison-Wesley, Reading, MA, 1977).[116] T. Pearcey, The structure of an electromagnetic fieldin the neighborhood of a caustic. Phil. Mag. , 311(1946).[117] NIST Handbook of Mathematical Functions , edited byOlver et al. (Cambridge University, New York, 2010),chapter 36. Available online at dlmf.nist.gov[118] M. V. Berry, J. F. Nye, and F. J. Wright, The ellipticumbilic diffraction catastrophe, Phil. Trans. R. Soc. A , 453 (1979).[119] M. Berry,
Singularities in Waves and Rays in LesHouches, Session XXXV, 1980
Physics of Defects , editedby R. Balian et al. (North-Holland Publishing, Amster-dam, 1981).[120] R. Da Silveira, Rainbow interference effects in heavy ionelastic scattering, Phys. Lett. B , 211 (1973).[121] T. C. Petersen, M. Weyland, D. M. Paganin, T. P. Sim-ula, S. A. Eastwood, and M. J. Morgan, Electron Vor-tex Production and Control Using Aberration InducedDiffraction Catastrophes, Phys. Rev. Lett. , 033901(2013).[122] W. Rooijakkers, S. Wu, P. Striehl, M. Vengalattore, andM. Prentiss, Observation of caustics in the trajectoriesof cold atoms in a linear magnetic potential, Phys. Rev.A , 063412 (2003).[123] J. H. Huckans, I. B. Spielman, B. L. Tolra, W. D.Phillips, and J. V. Porto, Quantum and classical dy-namics of a Bose-Einstein condensate in a large-periodoptical lattice. Phys. Rev. A , 043609 (2009).[124] S. Rosenblum, O. Bechler, I. Shomroni, R. Kaner, T.Arusi-Parpar, O. Raz, and B. Dayan, Demonstration ofFold and Cusp Catastrophes in an Atomic Cloud Re-flected from an Optical Barrier in the Presence of Grav-ity, Phys. Rev. Lett. , 120403 (2014).[125] R. Plestid, P. Mahon, and D. H. J. O’Dell, Violent re-laxation in quantum fluids with long-range interactions,Phys. Rev. E , 012112 (2018).[126] U. Leonhardt, A laboratory analogue of the event hori-zon using slow light in an atomic medium, Nature (Lon-don) , 406 (2002).[127] M. V. Berry and M. R. Dennis, Quantum cores of opticalphase singularities, J. Opt. A , S178 (2004).[128] W. Kirkby, J. Mumford, and D.H.J. O’Dell, Quantumcaustics and the hierarchy of light cones in quenchedspin chains, Phys. Rev. Research , 033135 (2019).[129] A. Z. Goldberg, A. Al-Qasimi, J. Mumford, and D.H. J. O’Dell, Emergence of singularities from decoher-ence: Quantum catastrophes, Phys. Rev. A , 063628(2019). [130] L. S. Schulman, Tricritical Points and Type-ThreePhase Transitions, Phys. Rev. B , 1960 (1973).[131] I. Stewart and T. Poston, Catastrophe Theory and ItsApplications (Dover, New York, 2012).[132] R. Gilmore, Structural Stability of the Phase Transitionin Dicke-like Models, J. Math. Phys. , 17 (1977).[133] R. Gilmore and L. M. Narducci, Relation between theEquilibrium and Nonequilibrium Critical Properties ofthe Dicke Model, Phys. Rev. A , 1747 (1978).[134] R. Gilmore and L. M. Narducci, Laser as Catastrophe,in: Coherence and Quantum Optics (L. Mandel and E.Wolf, Eds.), (Plenum, New York, 1977).[135] R. Gilmore, Catastrophe Theory for Scientists and En-gineers , (Wiley, New York, 1981).[136] F. A. van Abeelen and B. J. Verhaar, Time-dependentFeshbach resonance scattering and anomalous decay ofa Na Bose-Einstein condensate, Phys. Rev. Lett. ,1550 (1999).[137] A. J. Leggett, Bose-Einstein condensation in the alkaligases: Some fundamental concepts, Rev. Mod. Phys. , 307 (2001).[138] R. Gati and M. K. Oberthaler, A bosonic Josephsonjunction, J. Phys. B: At. Mol. Opt. Phys. R61 (2007).[139] C. V. Raman, N. S. N. Nath, The diffraction of lightby high frequency sound waves: part IV. Proc. IndianAcad. Sci. , 119 (1936).[140] M. V. Berry, The Diffraction of Light by Ultrasound (Academic, New York, 1966).[141] D. H. J. O’Dell, Dynamical diffraction in sinusoidal po-tentials: uniform approximations for Mathieu functions,J. Phys. A: Math. Gen. , 3897 (2001).[142] S. Mossmann and C. Jung, Semiclassical approach toBose-Einstein condensates in a triple well potential,Phys. Rev. A. , 033601 (2006).[143] S. Raghavan, A. Smerzi, S. Fantoni, and S. R. Shenoy,Coherent oscillations between two weakly coupled Bose-Einstein condensates: Josephson effects, π -oscillations,and macroscopic quantum self-trapping, Phys. Rev. A , 620 (1999).[144] J. R. Anglin, P. Drummond, and A. Smerzi, Exact quan-tum phase model for mesoscopic Josephson junctions,Phys. Rev. A , 063605 (2001).[145] A. van Otterlo, K. H. Wagenblast, R. Baltin, C. Bruder,R. Fazio, and G. Sch¨on, Quantum-Phase Transitionsof Interacting Bosons and the Supersolid Phase, Phys.Rev. B , 16176 (1995).[146] R. Fazio, H. van der Zant, Quantum phase transi-tions and vortex dynamics in superconducting networks,Phys. Rep. , 235 (2001).[147] F. S. Cataliotti, S. Burger, C. Fort, P. Maddaloni, F.Minardi, A. Trombettoni, A. Smerzi, and M. Inguscio,Josephson Junction Arrays with Bose-Einstein Conden-sates, Science , 843 (2001).[148] C. Orzel, A. K. Tuchman, M. L. Fenselau, M. Yasuda,and M. A. Kasevich, Squeezed States in a Bose-EinsteinCondensate, Science , 2386 (2001).[149] Z. Hadzibabic, S. Stock, B. Battelier, V. Bretin, andJ. Dalibard, Interference of an Array of IndependentBose-Einstein Condensates, Phys. Rev. Lett. , 180403(2004).[150] C. Schori, T. St¨oferle, H. Moritz, M. K¨ohl, and T.Esslinger, Excitations of a Superfluid in a Three-Dimensional Optical Lattice, Phys. Rev. Lett. 93,240402 (2004). [151] K. Xu, Y. Liu, D. E. Miller, J. K. Chin, W. Setiawan,and W. Ketterle, Observation of Strong Quantum De-pletion in a Gaseous Bose-Einstein Condensate, Phys.Rev. Lett. , 180405 (2006).[152] Schweikhard, V., S. Tung, and E. A. Cornell, Vor-tex proliferation in the Berezinskii-Kosterlitz-Thoulessregime on a two-dimensional lattice of Bose-Einsteincondensates, Phys. Rev. Lett. , 030401 (2007).[153] A. D. Greentree, J. H. Cole, A. R. Hamilton, and L. C.L. Hollenberg, Coherent electronic transfer in quantumdot systems using adiabatic passage, Phys. Rev. B ,235317 (2004).[154] J. H. Cole, A. D. Greentree, L. C. L. Hollenberg, andS. Das Sarma, Spatial adiabatic passage in a realistictriple well structure, Phys. Rev. B , 235418 (2008).[155] M. Rab, J. H. Cole, N. G. Parker, A. D. Greentree, L.C. L. Hollenberg, and A. M. Martin, Spatial coherenttransport of interacting dilute Bose gases, Phys. Rev. A , 061602(R) (2008).[156] T. Opatrn´y and K. K. Das, Conditions for vanishingcentral-well population in triple-well adiabatic trans-port, Phys Rev. A. , 012113 (2009).[157] C. J. Bradly, M. Rab, A. D. Greentree, and A. M. Mar-tin, Coherent tunneling via adiabatic passage in a three-well Bose-Hubbard system, Phys. Rev. A , 053609(2012).[158] M. K. Olsen, Quantum dynamics and entanglement incoherent transport of atomic population, J. Phys. B: At.Mol. Opt. Phys. (2014) 095301.[159] J. A. Stickney, D. Z. Anderson, and A. A. Zozulya, Tran-sistorlike behavior of a Bose-Einstein condensate in atriple-well potential, Phys. Rev. A , 013608 (2007).[160] K. W. Wilsmann, L. H. Ymai, A. P. Tonel, J. Links,and A. Foerster, Control of tunneling in an atomtronicswitching device, Commun. Phys. , 91 (2018).[161] Z. Zhang, V. Dunjko, and M. Olshanii, Atom transistorfrom the point of view of nonequilibrium dynamics, NewJ. Phys. , 025010(2016).[163] G. Arwas, A. Vardi, and D. Cohen, Triangular Bose-Hubbard trimer as a minimal model for a superfluidcircuit, Phys Rev A , 013601 (2014).[164] G. Arwas, A. Vardi, and D. Cohen, Superfluidity andChaos in low dimensional circuits, Sci. Rep. , 13433(2015).[165] A. Gallem´ı, M. Guilleumas, J. Martorell, R. Mayol, A.Polls, and B. Juli´a-D´ıaz, Fragmented condensation inBose-Hubbard trimers with tunable tunnelling, New J.Phys. , 073014 (2015).[166] C. Lee, T. J. Alexander, and Y. S. Kivshar, Melting ofdiscrete vortices via quantum fluctuations, Phys. Rev.Lett. , 180408 (2006)[167] Q. Guo, X. Chen, and B. Wu, Tunneling dynamics andband structures of three weakly coupled Bose-Einsteincondensates, Opt. Express , 19219 (2014).[168] S. Dutta, M. C. Tsatsos, S. Basu, and A. U. J. Lode,Management of the correlations of ultracold bosons intriple wells, New J. Phys. ,033627 (2015). [170] M. K. Olsen, Asymmetric steering in coherent transportof atomic population with a three-well Bose-Hubbardmodel, J. Opt. Soc. Am. B , A15 (2015).[171] J. K. Kalaga, A. Kowalewska-Kud(cid:32)laszyk, W. Leo´nski,and A. Barasi´nski, Quantum correlations and entangle-ment in a model comprised of a short chain of nonlinearoscillators, Phys. Rev. A , 032304 (2016).[172] J. K. Kalaga, W. Leo´nski, and R. Szcz¸e´sniak, Quantumsteering and entanglement in three-mode triangle Bose-Hubbard system, Quantum Inf. Process. , 265 (2017).[173] M. K. Olsen, T. W. Neely, and A. S. Bradley, Meso-scopic dynamical differences from quantum state prepa-ration in a Bose-Hubbard trimer, Phys. Rev. Lett. ,230406 (2018).[174] A. A. Bychek, P. S. Muraev, D. N. Maksimov, E. N.Bulgakov, and A. R. Kolovsky, Chaotic and regular dy-namics in the three-site Bose-Hubbard model, AIP Con-ference Proceedings , 020007 (2020).[175] M. V. Berry and M. Robnik, Semiclassical level spacingswhen regular and chaotic orbits coexist, J. Phys. A ,2413 (1984).[176] M. V. Berry, Semiclassical Mechanics of regular and ir-regular motion, p171 in Les Houches Lecture Series Ses-sion XXXVI, eds. G. Iooss, R. H. G. Helleman and R.Stora (North Holland, Amsterdam, 1983).[177] V. I. Arnold, Mathematical Methods of Classical Me-chanics, Second Edition (Springer, Berlin, 1989).[178] J. F. Nye, Wave dislocations in the diffraction patternof a higher-order optical catastrophe, J. Opt. (2010)015702.[179] C. Emary, N. Lambert, and T. Brandes, Entanglementin quantum catastrophes, Phys. Rev. A , 062302(2005).[180] A. Sinatra, C. Lobo, and Y. Castin, Classical-FieldMethod for Time Dependent Bose-Einstein CondensedGases, Phys. Rev. Lett. , 210404 (2001).[181] A. Sinatra, C. Lobo, and Y. Castin, The truncatedWigner method for Bose condensed gases: limits of va-lidity and applications, J. Phys. B , 3599 (2002).[182] A. Polkovnikov, Phase space representation of quantumdynamics, Ann. Phys. (Amsterdam) , 1790 (2010).[183] J. Javanainen and J. Ruostekoski, Emergent classicalityin continuous quantum measurements, New J. Phys. ,013005 (2013).[184] G. B. Airy, On the intensity of light in the neighbour- hood of a caustic, Trans. Cambridge Philos. Soc. 6, 379(1838).[185] I. Zapata, F. Sols and A. J. Leggett, Phase dynamicsafter connection of two separate Bose-Einstein conden-sates, Phys. Rev. A , 021603(R) (2003).[186] R. Borghi, Numerical computation of diffraction catas-trophes with codimension eight, Phys. Rev. E ,046704 (2012).[187] M. V. Berry and C. J. Howls, Axial and focal-planediffraction catastrophe integrals, J. Phys. A , 375206(2010).[188] J. F. Nye, The catastrophe optics of liquid drop lenses,Proc. R. Soc. Lond. A , 1 (1986).[189] M. V. Berry and J. F. Nye, Fine Structure in CausticJunctions, Nature , 34 (1977).[190] C. Nicolis, Stochastic resonance in multistable systems:The role of dimensionality, Phys. Rev. E. , 011133(2012).[191] J. Callahan, The double cusp has five minima, Proc.Cambridge Philos. Soc. , 537 (1978).[192] J. Callahan, Special bifurcations of the double cusp,Proc. London Math. Soc. , 227 (1982).[193] C. Upstill , F. J. Wright , J. V. Hajnal, and R. H. Tem-pler, The double-cusp unfolding of the X diffractioncatastrophe, Optica Acta , 1651 (1982).[194] A. R. Kolovsky, Bose-Hubbard Hamiltonian: QuantumChaos approach, A. R. Kolovsky, Int. J. Mod. Phys. B , 1630009 (2016).[195] L. S. Schulman, Techniques and Applications of PathIntegration , (Wiley, New York, 1981). Caustics and mul-tivaluedness: Two results of adding path amplitudes, in“Functional Integration and Its Applications” (A. M.Arthurs, Ed.), Oxford Univ. Press XXX(Clarendon),London, 1975.[196] G. Dangelmayr and W. Veit, Semiclassical Approxima-tion of Path Integrals on and near Caustics in Terms ofCatastrophes, Ann. Phys. , 108 (1979).[197] J. F. Nye, D. R. Hawst, and R. A. Smith, Use of diffrac-tion gratings with curved lines to study the opticalcatastrophes D +6 and D − . J. Mod. Opt. , 407 (1987).[198] M. R. Andrews, C. G. Townsend, H.-J. Miesner, D. S.Durfee, D. M. Kurn, and W. Ketterle, Observation ofInterference Between Two Bose Condensates, Science275