Floquet-engineered nonlinearities and controllable pair-hopping processes: From optical Kerr cavities to correlated quantum matter

This work explores the possibility of creating and controlling unconventional nonlinearities by periodic driving, in a broad class of systems described by the nonlinear Schr\"odinger equation (NLSE). By means of a parent quantum many-body description, we demonstrate that such driven systems are well captured by an effective NLSE with emergent nonlinearities, which can be finely controlled by tuning the driving sequence. We first consider a general class of two-mode nonlinear systems - relevant to optical Kerr cavities, waveguides and Bose-Einstein condensates - where we find an emergent four-wave mixing nonlinearity, which originates from pair-hopping processes in the parent quantum picture. Tuning this drive-induced nonlinearity is shown to modify the phase-space topology, which can be detected through relative population and phase measurements. We then couple individual (two-mode) dimers in view of designing extended lattice models with unconventional nonlinearities and controllable pair-hopping processes. Following this general dimerization construction, we obtain an effective lattice model with drive-induced interactions, whose ground-state exhibits orbital order, chiral currents and emergent magnetic fluxes through the spontaneous breaking of time-reversal symmetry. We analyze these intriguing properties both in the weakly-interacting (mean-field) regime, captured by the effective NLSE, and in the strongly-correlated quantum regime. Our general approach opens a route for the engineering of unconventional optical nonlinearities in photonic devices and controllable drive-induced interactions in ultracold quantum matter.

Our analysis is based on a parent quantum many-body Hamiltonian, which describes two species of interacting bosons subjected to a periodic pulse sequence; see the sketch in Fig. 1.This theoretical framework captures the physics of driven nonlinear optical settings in the classical (mean-field) limit [65,77,78].Specifically, we first derive an effective quantum Hamiltonian that well describes the stroboscopic dynamics of the driven parent quantum system, in the highfrequency regime of the pulse sequence [13,[79][80][81][82]. From this, we then derive the effective classical equations of motion, hence revealing the effective optical nonlinearities generated by the driving sequence [Section III].We explore the validity of both approximations (i.e. the high-frequency approximation related to the drive and the mean-field approximation associated with the classical limit) through numerical simulations of the quantum and classical dynamics, comparing the full time dynamics generated by the pulse sequence to the effective descriptions [Section IV].As a by-product, this analysis illustrates how the effective nonlinearities can be detected through the dynamics of simple observables: the relative intensity and relative phase of the two optical modes.We then discuss how the strength of the effective optical nonlinearities can be tuned by simply adjusting the pulse sequence [Section V].This control over drive-induced optical nonlinearities is directly reflected in the topology of the system's phase space, which can be detected through light intensity and phase measurements. .Schematic of the approach.We consider a general class of two-mode nonlinear systems, driven by a periodic pulse sequence and described by the nonlinear Schrödinger equation (NLSE).To analyse these settings, we introduce a parent quantum many-body Hamiltonian, Ĥ0 + V (t), which describes two species of interacting bosons driven by a periodic pulse sequence.From this, we derive an effective quantum Hamiltonian Ĥeff in the high-frequency limit of the drive (ω → ∞).We then derive the effective nonlinear Schrödinger equation upon taking the classical limit, N → ∞, where N is the number of bosons, hence revealing the effective nonlinearities generated by the driving pulse sequence.

II. THE DRIVEN NONLINEAR SYSTEM AND ITS EFFECTIVE DESCRIPTION
We consider a class of two-mode nonlinear systems, described by the (possibly coupled) nonlinear Schrödinger equations Here, ψ 1,2 (x, t) denote the complex amplitude of the fields corresponding to the modes σ = 1, 2; they depend on the evolution "time" t and the "spatial" coordinate x.The focus of this work is set on the "internal" dynamics associated with the two modes, such that the "spatial" coordinate x [and the related kinetic-energy term ∼ γ in Eq. ( 1)] does not play any role in the following.For the sake of generality, the equations of motion (1) contain two types of nonlinearities, which are generically present in optical cavities [59][60][61]: the so-called self-phase modulation and the cross-phase modulation, whose respective strengths are set by the parameter β; we have also included a static linear coupling of strength Ω 0 /2.We point out that the nonlinear equations (1) are decoupled in the limit Ω 0 = β = 0, i.e. in the absence of linear coupling and crossphase modulation.While Eq. ( 1) naturally describes the two polarization modes of a light field propagating in a lossless cavity [59][60][61], or light propagating in a pair of adjacent waveguides [37,46], it should be noted that Eq. ( 1) also captures the physics of quantum gases trapped in a double well potential and twocomponent Bose-Einstein condensates [63,64].
In order to create effective nonlinearities in Eq. ( 1), we now include a time-periodic pulse sequence of period T , which mixes the two modes in a stroboscopic manner: , where n ∈ N, the two components undergo the mixing operation • Pulse : at times t n = T ×n, the system undergoes the reverse operation In a two-mode optical cavity [59][60][61], these pulsed operations would correspond to a coupling between the two polarization eigenmodes of the cavity, as directly realized by means of quarter-wave plates [83,84]; see the sketch in Fig. 2(a).In this case, the "time" coordinate t should be interpreted as the propagation distance along the cavity [85].
More generally, when the mixing processes in Eqs. ( 2)-(3) cannot be directly performed by a device, they can also be realized by activating a linear coupling between the two modes, during a short pulse duration τ T , such that the equations of motion of the driven system can be written in the form Here, the function Ω(t) = Ω 0 − f pulse (t) includes the pulse sequence defined by the function [Fig.3(a)] To verify that the drive in Eqs. ( 4)-( 5) indeed realizes the mixing operations in Eqs. ( 2)-(3), we restrict ourselves to the (linear) driving terms in the coupled Schrödinger equations (4) and we obtain the time-evolution operators corresponding to the first and second pulses, respectively: where σx is the standard Pauli matrix.The operators Ûmix and Û † mix in Eq. ( 6) indeed realize the mixing operations in Eqs. ( 2)-(3), respectively.We note that these mixing operations are reminiscent of π/2 pulses in quantum optics [86][87][88].We also point out that the choice of the pulse function in Eq. ( 5) is not unique, e.g. the amplitude of the second pulse ( ) can be set to the value (−π/2 + 2πp)/τ , with p ∈ Z, without affecting the dynamics.The advantage of the pulse function proposed in Eq. ( 5) is that the linear coupling does not change sign over time, which can be convenient for certain physical realizations.
In optical-waveguides settings [37], the two modes (1 and 2) would describe light propagating in two adjacent waveguides.In this case, the pulsed linear couplings in Eqs. ( 4)-( 5) can be realized by abruptly changing the spatial separation between the two waveguides; see Fig. 2(b) for a sketch and Refs.[38,39,41,42] for experimental realizations using ultrafast-laser-inscribed waveguides.Such opticalwaveguides settings could benefit from the state-recycling technique of Refs.[40,89], where light is re-injected into the waveguides (and possibly modified) at every roundtrip; see also Refs.[90][91][92] regarding setups based on recirculating fiber loops.
In the context of quantum gases trapped in a double well potential, the linear coupling between the two neighboring sites (orbitals) can be activated in a pulsed manner through a dynamical variation of the potential barrier [93]; see also Ref. [88], where fast π/2 pulses were implemented in a twocomponent Bose-Einstein condensate. ) with modulated inter-waveguide separation, realizing a "time-periodic" linear coupling between the two optical modes [Eq.( 4)].In both cases, the "time" coordinate corresponds to the propagation direction [37,85].
In the limit of a fast pulse sequence, namely, when the period of the drive T T eff is much smaller than the effective "time" scale of the system (to be discussed below), the stroboscopic time-evolution of the nonlinear system is found to be well described by the effective equations of motion where we introduced the quantity χ = 1 − β, and where the system is assumed to be measured stroboscopically (t = T ×n, with n ∈ N).Comparing Eq. ( 7) with the original Eq. ( 1), we find that the repeated mixing processes in Eqs. ( 2)-( 3) effectively produce a new form of nonlinearity, commonly known in optics as four-wave mixing [94,95].The drive also renormalizes the self-phase-modulation and it effectively annihilates the cross-phase modulation.We point out that the effective four-wave mixing is induced even in the limit of two initially decoupled modes (β = Ω 0 = 0).It is the aim of the following Sections III-IV to demonstrate the effective description displayed in Eq. ( 7) and to explore its regimes of validity.We then introduce an "imbalanced" pulse sequence in Section V, which allows one to tune the relative strengths of effective nonlinearities and induce topological changes in phase space.While this work sets the focus on the Floquet engineering of classical nonlinear systems, the results obtained below also apply to the quantum dynamics of driven interacting bosonic systems.
As a technical note, we point out that the mixing processes in Eqs. ( 2)-( 3) do not modify the kinetic-energy terms in Eq. ( 1).For the sake of presentation, we henceforth set γ = 0 (except otherwise stated), but we do keep in mind that these terms can be readily added in the description without affecting the results [96].4)-( 5).(b) The pulse sequence associated with the time-evolution operator in Eq. ( 18), which involves stroboscopic mixing operations Û ( †) mix separated by "free" time evolution ( Ĥ0).

III. A QUANTUM MANY-BODY APPROACH
Our approach consists in three successive steps [Fig.1]: • We introduce a parent quantum many-body Hamiltonian, whose semiclassical dynamics reproduces the time evolution of the driven nonlinear system in Eq. ( 4); • Within this quantum framework, we derive the effective (Floquet) Hamiltonian that well captures the long time dynamics in the high-frequency limit (2π/T → ∞); • We then obtain the effective classical equations of motion from the effective quantum Hamiltonian.
The validity of this approach will then be verified in Section IV, through numerical studies of both quantum and classical dynamics.
A. The parent quantum many-body system Our starting point is the quantum many-body Hamiltonian where â † σ (resp.âσ ) creates (resp.annihilates) a boson in the mode σ = 1, 2. These operators satisfy the bosonic commutation relations, [â σ , â † σ ] = δ σ,σ .The first line in Eq. ( 8) describes intra-mode (Hubbard) interactions, while the second line describes inter-mode (cross) interactions of strength β; the Hamiltonian also includes single-particle hopping processes of amplitude Ω 0 /2; see Fig. 4(a)-(c) for a sketch of the processes and Refs.[97,98].Henceforth, the Hubbard interaction strength U = 1 sets our unit of energy, as well as our unit of time t unit = /U .First of all, we note that the classical equations of motion in Eq. ( 1) are readily obtained from Heisenberg's equations, dâ σ /dt = i[ Ĥ0 , âσ ], upon taking the classical limit â1,2 → ψ 1,2 ; see Refs.[66,77,78].Specifically, the selfphase modulation in Eq. ( 1) stems from the intra-mode (Hubbard) interaction terms in Eq. ( 8), while the cross-phase modulation stems from the inter-mode (cross) interaction term.Hence, this justifies the choice of Eq. ( 8) as a proper parent quantum Hamiltonian for our initial (non-driven) system.Note that we set = 1 throughout this work.
In fact, for the sake of later convenience, it is instructive to derive the classical equations of motion in Eq. ( 1) using a different approach.Indeed, this will allow us to introduce central notions and quantities, which will be used throughout this work.Let us introduce a set of angular momentum (Schwinger) operators [99], defined as These operators satisfy the spin commutation relations [ Ĵµ , Ĵν ] = iε µνλ Ĵλ , and the operator N counts the total number of bosons in the system (assumed to be constant); note that Ĵµ = σµ /2 for a single boson (N = 1), where σx,y,z denote the Pauli matrices.Using the operators in Eq. ( 9), the parent Hamiltonian in Eq. ( 8) simply reads and we henceforth neglect the constant terms (proportional to N and N 2 ); see Appendix A. We note that the Hamiltonian in Eq. ( 10) has been extensively studied in the context of the bosonic Josephson effect [64,100] and nuclear physics [101].
The equations of motion associated with Eq. ( 10) are readily obtained from Heisenberg's equations In order to connect Eq. ( 11) to the classical nonlinear Schrödinger equation in Eq. ( 1), we take the classical limit and introduce the Bloch-Poincaré sphere representation (θ, ϕ) through the mapping Injecting this into Eq.( 11), one obtains the classical equations of motion for the two canonical conjugate variables z(t) and ϕ(t) [63,64].We point out that Eq. ( 13) is equivalent to the nonlinear Schrödinger equation in Eq. ( 1) upon representing the complex amplitudes ψ 1,2 on the Bloch-Poincaré sphere [78,102] where we introduced the relative phase ϕ between the two modes, the relative population (or relative light intensity) and the total population (or total light intensity) We emphasize that the dynamics in phase space, i.e. the trajectories (z(t), ϕ(t)), can be simply monitored in an optical setting by measuring the light intensity and the relative phase of the two modes.
For the sake of completeness, we note that the equations of motion in Eq. ( 13) can be derived from Hamilton's equation, using the classical Hamiltonian [63,64,103] The classical dynamics of the non-driven system hence relies on a competition between the "mean-field" interaction parameter g = χN and the linear coupling Ω 0 .This competition is at the core of bifurcations and symmetry breaking in the bosonic Josephson effect [59,64,78].

B. The pulse sequence and the effective Floquet Hamiltonian
We now introduce the quantum-many-body analogue of the pulse sequence introduced in Eqs. ( 2)- (5).We write the timeevolution operator over one period T in the form [Fig. 3(b)] where the mixing operator is defined as We note that this indeed corresponds to the π/2-pulse operator in Eq. ( 6) for a single boson (N = 1), which is consistent with the fact that the mixing operation is a single-particle process.
We also point out that we explicitly took the limit τ → 0, where τ is the pulse duration; see Eq. ( 5).
The state of the quantum many-body system at time t n = T ×n is then obtained as where |ψ(0) denotes the initial state of the system.We now derive the effective (Floquet) Hamiltonian [13,79,80], which captures the stroboscopic dynamics of the driven system, and hence, its time evolution over long time scales t n T .The effective Hamiltonian is defined through the time-evolution operator over one period [13,104] and it can be evaluated explicitly through a 1/ω-expansion, where ω = 2π/T denotes the drive frequency; see Refs.[13,[79][80][81][82].In order to reach convergence of this infinite series expansion, we partially resum the series [13] by splitting the time-evolution operator in Eq. ( 18) into two parts where we introduced the operator Ĥ1 defined as Then, assuming that T ω eff 1, where ω eff is the characteristic frequency associated with the processes included in the Hamiltonians Ĥ0 and Ĥ1 , we apply the Trotter approximation to Eq. ( 22), from which we directly obtain the effective Hamiltonian [Eq.( 21)] Our problem of finding the effective Hamiltonian thus reduces to the calculation of Ĥ1 defined in Eq. ( 23).This step can be performed exactly, by noting that where we used the definition of Ĥ0 in Eq. (10).Using the Baker-Campbell-Hausdorff formula, one obtains [71] The effective Hamiltonian in Eq. ( 25) finally reads From this result, we find that the Trotter approximation [Eq.(24)] is valid for a sufficiently short driving period satisfying T 1/χ and T 1/Ω 0 .It is instructive to rewrite the effective Hamiltonian in Eq. ( 29) using the original bosonic operators [Appendix A], A comparison with the initial Hamiltonian Ĥ0 in Eq. ( 8) indicates that the driving pulse sequence has effectively generated novel interaction terms; see the second line of Eq. ( 30).These "pair-hopping" terms [98,105] describe processes by which two particles in mode σ collide and end up in the other mode σ = σ; see Fig. 4(d).As we now discuss below, these pairhopping terms are at the origin of the four-wave mixing nonlinearity announced in Eq. (7).We also note that the effective interaction strength is given by U eff = χ/8 = (1−β)/8, where β sets the strength of the inter-mode (cross) interaction in the initial Hamiltonian Ĥ0 in Eq. ( 8).

C. Effective classical equations of motion
First of all, we find that the effective nonlinear Schrödinger equation in Eq. ( 7) is directly obtained from the effective Hamiltonian Ĥeff in Eq. ( 30), using Heisenberg's equations dâ σ /dt = i[ Ĥeff , âσ ], and upon taking the classical limit â1,2 → ψ 1,2 .In particular, the effective four-wave mixing in Eq. ( 7) originates from the effective pair-hopping terms in Eq. (30).
In analogy with Eqs. ( 11)-( 13), we explicitly derive the classical equations of motion for the two canonical conjugate variables z(t) and ϕ(t), describing the relative population and phase of the two modes.Using the effective Hamiltonian in Eq. ( 29) and Heisenberg's equations, we find Finally, applying the Bloch-Poincaré-sphere mapping in Eq. ( 12), we obtain the classical equations of motion We find that the equations of motion (31) can be derived from Hamilton's equation, using the classical Hamiltonian We stress that the classical equations of motion in Eq. ( 31) are physically equivalent to the effective nonlinear Schrödinger equation announced in Eq. ( 7), through the mapping provided by Eq. ( 14).

IV. NUMERICAL ANALYSIS
This Section aims at exploring the validity of the effective-Hamiltonian analysis developed in Section III B and its classical limit presented in Section III C.

A. Validating the effective quantum Hamiltonian
First, we demonstrate that the dynamics associated with the effective Hamiltonian in Eq. ( 30) reproduces the stroboscopic dynamics of the driven system described by Eqs. ( 18)- (20).To this end, we choose a coherent spin state as an initial state [87,100] which corresponds to a macroscopic occupation of the singleparticle state, defined on the Bloch sphere.Here we introduced the singleparticle states |1 = â † 1 |∅ and |2 = â † 2 |∅ , associated with the two modes, as well as the creation operator â † θ,ϕ |∅ = |θ, ϕ .We note that the chosen initial state in Eq. ( 33) behaves classically in the limit N → ∞ [100], which will be convenient for later purposes (i.e. when comparing quantum and classical dynamics).
We analyze the quantum dynamics through the evaluation of the population imbalance where the time-evolved state |ψ(t n ) is obtained from: (i) the full time dynamics of the driven system [Eqs.( 18)-( 20)], and (ii) the effective Hamiltonian [Eq.(30)].Figure 5 compares these two results for both N = 10 and N = 50 bosons, and the same interaction parameter g = χN = 5.In both cases, one obtains that the effective description well captures the stroboscopic dynamics when the driving period is sufficiently small, T 0.1 in the current units [see Eq. ( 8)].This analysis validates the effective Hamiltonian in Eq. ( 30) in the highfrequency regime.For each case, the full time dynamics of the driven system is generated using the sequence in Eq. ( 18) with a period T = 0.2, T = 0.1 and T = 0.05.Here the interaction parameter is set to g = χN = 5 and the static linear coupling is set to Ω0 = 0; the initial coherent spin state |N, θ, ϕ corresponds to z = cos θ = 0.4 and ϕ = 2.25.In all plots, the time-evolved state is evaluated at stroboscopic times tn = T ×n.

B. The effective semiclassical dynamics
As a next step, we now show that the effective Hamiltonian Ĥeff in Eq. ( 30) well captures the classical dynamics generated by the equations of motion in Eq. ( 31).We remind that the latter classical description is associated with the Hamiltonian function H eff (z, ϕ) displayed in Eq. ( 32), where z and ϕ describe the relative population and phase of the two modes; see Eqs. ( 14)- (15).The agreement between the quantum and classical descriptions is expected to be reached in the limit N → ∞, where quantum fluctuations become negligible [63,64,77,100,106,107].We also remind the reader that the classical equations of motion in Eq. ( 31), which are analyzed in this Section, are equivalent to the effective nonlinear Schrödinger equation in Eq. ( 7), through the mapping defined in Eq. ( 14).
First of all, let us analyze the dynamics generated by the effective classical equations of motion in Eq. (31).In order to highlight the role of nonlinearities, we hereby set the static linear coupling to Ω 0 = 0.In Fig. 6, we display a few representative trajectories over the energy landscape H eff (z, ϕ) defined in Eq. ( 32).These trajectories reflect the presence of two stable fixed points at (z = 0, ϕ = 0) and (z = 0, ϕ = π).We stress that this configuration of fixed points radically differs from that associated with the non-driven system [see H 0 (z, ϕ) in Eq. ( 17)] for the same choice of Ω 0 = 0. We now compare these classical predictions to the quantum dynamics associated with the effective Hamiltonian Ĥeff in Eq. ( 30), using a coherent spin state |N, θ, ϕ as an initial condition; see Eq. ( 33). Figure 7 shows the trajectories z(t) for N = 5, 10, 80, 170 bosons, while keeping the "mean-field" interaction parameter χN = 5 constant.From these results, we confirm that a good agreement between the effective classical and quantum descriptions is indeed obtained in the large N limit.
In order to further appreciate the residual deviations between the quantum and classical dynamics in the small N regime, we depict the time-evolving Husimi function Q(z, ϕ; t) in Fig. 8 for the case N = 80.The Husimi function [64,100,[108][109][110][111][112] is obtained by evaluating the squared overlap of the time-evolving state |ψ(t) with the coherent spin states defined over the Bloch sphere (with same particle number N ), Here the state |ψ(t) is evolved according to the effective Hamiltonian Ĥeff in Eq. ( 30), so that the evolution of the Husimi function in Fig. 8 is to be compared with the quantum dynamics displayed in Fig. 7(c) for N = 80 bosons.The timeevolution of the Husimi function Q(z, ϕ; t) shown in Fig. 8 indicates that the initial coherent spin state |ψ(0) = |N, θ, ϕ becomes substantially squeezed [87] around time t ≈ 3, which also corresponds to the time around which the classical trajectory starts deviating from the effective-Hamiltonian quantum dynamics in Fig. 7(c).At later times, t ≈ 12, the state becomes oversqueezed and it exhibits Majorana stars in the Husimi distribution [109,111,112].We find that these non-classical features are postponed to later evolution times upon increasing the number of bosons N while keeping the interaction parameter g = χN fixed.Despite these non-classical features, the center of mass of the Husimi function is found to approximately follow a classical orbit around the stable fixed point (z = 0, ϕ = π), as depicted in Fig. 6.

C. The driven nonlinear Schrödinger equation and its effective description
In this Section, we analyze the agreement between the classical dynamics associated with the driven nonlinear Schrödinger equation [Eqs.(4)-( 5)] and the dynamics generated by the effective classical equations of motion [Eq.(31)], which derive from the Hamiltonian H eff (z, ϕ) in Eq. (32).We remind that these effective equations of motion are equivalent to the effective nonlinear Schrödinger equation announced in Eq. (7).30).Here, the number of bosons is N = 80, and the other parameters are the same as in Fig. 7(c).The initial coherent spin state |ψ(0) = |N, θ, ϕ , at z = cos θ = 0 and ϕ = 2.7, becomes substantially squeezed around time t ≈ 3, hence signaling the breakdown of its classical description.An oversqueezed state, exhibiting Majorana stars, appears around t ≈ 12.The trajectory predicted by the effective classical equations of motion [Eq.( 31)] is depicted in white.
In practice, we numerically solve the classical equations of motion [Eq.( 13)] where the pulse function f pulse (t) is defined in Eq. ( 5).These equations of motion are equivalent to the driven nonlinear Schrödinger equation in Eqs. ( 4)-( 5) through the mapping provided by Eq. ( 14).
The resulting dynamics are displayed in Fig. 9, together with the dynamics generated from the effective classical Hamiltonian H eff (z, ϕ) in Eq. (32).The results in Fig. 9 confirm that the effective classical description very well captures the dynamics of the driven nonlinear system at stroboscopic times t = t n , while a finite micromotion is observed at intermediate times t = t n .
Altogether, the numerical studies presented in this Section IV validate the effective description announced in Eq. (7) [see also Sections III B and III C], and hence, confirm the creation of effective interactions and nonlinearities through the repeated pulse sequence.

V. TUNING INTERACTION PROCESSES AND CLASSICAL NONLINEARITIES
A. The imbalanced pulse sequence The pulse sequence introduced in Eq. ( 18) corresponds to a balanced four-step sequence, with a free-evolution duration  set to T /2 during the first and third steps of the sequence [Fig.3(b)].However, it is instructive to consider the "imbalanced" sequence where the parameter α quantifies the imbalance; note that α = 1/2 for the balanced sequence in Eq. (18).Following the approach of Section III B, the effective Hamiltonian in Eqs. ( 25) and ( 29) is then generalized to At this stage, it is important to consider two limiting cases: when α = 1, one finds Ĥeff = Ĥ0 , which reflects the triviality of the sequence in Eq. ( 37) in this case.When α = 0, one finds the effective Hamiltonian which is thus strictly equivalent to the non-driven Hamiltonian Ĥ0 up to a unitary transformation [see Eq. ( 27)]: the Hamiltonians Ĥ0 and Ĥ(0) eff share the same spectrum.In fact, in the "pathological" case α = 0, the driving sequence simply generates an initial and final kick [13], as can be deduced by explicitly writing the time-evolving state at some arbitrary stroboscopic time t = t n [Eq.(20)] The long-time dynamics in Eq. ( 40) is indeed dictated by the static Hamiltonian Ĥ0 , but it is also affected by the initial kick e i π 2 Ĵx and the final kick e −i π 2 Ĵx .Altogether, one finds that non-trivial interaction (or nonlinearity) effects are generated by driving sequences corresponding to Ûα (T ; 0) in Eq. ( 37) with α = 0 and α = 1.One indeed verifies that Ĥ(α) eff in Eq. ( 38) and Ĥ0 do not have the same spectrum in this case.
B. The classical analysis: phase-space transitions and spontaneous symmetry breaking Following Section III C, we obtain the generalized classical equations of motion for the relative population and phase, These equations of motion are found to derive from the classical Hamiltonian We represent the corresponding trajectories in Fig. 10, for various values of the imbalance parameter α; in order to highlight the role of effective nonlinearities, we again set the static linear coupling to Ω 0 = 0. Interestingly, the system undergoes a succession of transitions as the imbalance parameter α is varied, which are characterized by changes in the topology of phase space [64]: When α = 0, the system is characterized by two stable classical fixed points at (z = 0, ϕ = π/2) and (z = 0, ϕ = 3π/2) [Fig.10(a)]; increasing α then generates two new stable fixed points at (z = 0, ϕ = 0) and (z = 0, ϕ = π) [Fig.10(b)]; the two initial fixed points at (z = 0, ϕ = π/2) and (z = 0, ϕ = 3π/2) then become unstable at α = 0.5 [Fig.10(c)], giving rise to two new stable fixed points located at the poles of the Bloch-Poincaré sphere z = ±1 [Fig.10(d)].We note that the emerging fixed points at z = ±1 are associated with the notion of spontaneous symmetry breaking, and were previously investigated in the context of ultracold gases [64,100] and in optical microcavities [59,78].In the present context, the symmetry breaking occurs as soon as the Ĵ2 z interaction term dominates over the Ĵ2 y interaction term; see the effective Hamiltonian Ĥ(α) eff in Eq. (38).We finally point out that the fixed points at (z = 0, ϕ = 0) and (z = 0, ϕ = π) become unstable when α = 1 (not shown in Fig. 10).

C. Tunable interactions and nonlinearities
We now rewrite the generalized effective Hamiltonian in Eq. ( 38) in terms of the original bosonic operators [Eq.( 9)], where the interaction strengths are given by From this, we find that the imbalanced pulse sequence offers an efficient method to control the strength and sign of the different interaction processes.This picture also offers an insightful view on the transition to spontaneous symmetry breaking discussed in Section V B: this transition, which takes place at α = 1/2, results from a competition between the intra-mode (Hubbard) interaction strength U 1 and the pair-hopping strength U 3 in Eq. (42), in the absence of linear coupling (Ω 0 = 0).This is different from the transition discussed in Refs.[59,64], which involves a competition between the Hubbard interaction strength U 1 and the linear coupling Ω 0 .Finally, in the classical limit (N → ∞), the effective nonlinear Schrödinger equation associated with the imbalanced pulse sequence reads where the three types of nonlinearities are controlled by the parameters U 1,2,3 displayed in Eq. ( 43).Consequently, the three types of effective nonlinearities can be tuned by adjusting the imbalanced pulse sequence.
In an optics setting, the modification of (effective) optical nonlinearities will directly manifest in the topology of phase space [Fig.10], which can be explored by extracting the trajectories (z(t), ϕ(t)) through light intensity and phase measurements.

VI. CONCLUDING REMARKS
This work proposed a method to engineer and tune nonlinearities in two-mode optical devices, using a designed pulse sequence that couples the optical modes in a fast and periodic manner.These repeated mixing operations simply correspond to the pulsed activation of a linear coupling between the two modes, and they can thus be implemented in a broad range of two-mode nonlinear systems, ranging from microresonators [59,78] and two-waveguide couplers [37,46] to circuit-QED platforms [11].While we considered a generic setting that includes both self-phase and cross-phase modulations in the absence of the periodic drive [Eq.( 1)], we found that effective nonlinearities emerge even when a single type of nonlinearity is present.Importantly, we demonstrated that the strength (and sign) of effective nonlinearities can be tuned by simply adjusting the pulse sequence.
To detect the emergence of drive-induced nonlinearities, we proposed to study changes in the phase space's topology [64], which can be explored by monitoring the dynamics of the relative intensity z(t) and phase ϕ(t) of the two optical modes.According to our numerical studies, these properties could already be revealed over "time" scales of the order of 5 − 10T , where T denotes the period of the driving sequence.This is particularly appealing for waveguide settings [37], where the "evolution time" associated with the propagation distance -and hence the number of driving periods -is limited.In this context, it would be interesting to combine such driving schemes with a state-recycling protocol [40].
While we considered a simple pulse sequence, characterized by the alternance of linear mixing operations and "free" evolution, we note that more complicated protocols and configurations could be envisaged.For instance, different types of mixing processes could be activated within each period of the drive, including nonlinear processes.Moreover, we note that similar driving schemes could be designed for Nmode devices, such as realized in arrays of ultrafast-laserinscribed waveguides [37].In the simplest scenario, one would couple the N modes by pairs, i.e. apply the pulse sequence in Eq. ( 18) for the neighboring modes 1 ↔ 2, 3 ↔ 4, . . ., (N − 1) ↔ N , and then for the complementary pairs 2 ↔ 3, 4 ↔ 5, . . ., (N − 2) ↔ (N − 1), and repeat this whole sequence in a fast and time periodic manner.In this context, it would be exciting to study the interplay of driveinduced nonlinearities and topological band structures; see for instance Ref. [113], where edge solitons were studied in the presence of four-wave mixing.
It would also be intriguing to explore the applicability of our scheme in the context of superconducting microwave cavities [114], where optical nonlinearities originate from the coupling to transmon ancillas.Indeed, it was recently shown that such optical nonlinearities can be modified by applying an offresonant drive on the transmon ancillas [58].Moreover, in circuit-QED platforms, the linear coupling between neighboring qubits can be modulated in a time-periodic manner [11]; applying our pulse protocol to such settings could be used to modify the nonlinearity of the qubits, and hence, the interaction between microwave photons.In general, we anticipate that drive-induced nonlinearities, such as the effective fourwave mixing studied in this work, could be useful for nonlinear optics applications [94,95].
Finally, we remark that the present work relies on a non-dissipative theoretical framework, where the number of bosons is conserved.Our scheme could nevertheless be applied to driven-dissipative optical devices [77], such as fiber ring cavities or microresonators described by the Lugiato-Lefever equation [61,[115][116][117][118], upon treating dissipation within the Floquet analysis [75,119,120].

Figure 1
Figure1.Schematic of the approach.We consider a general class of two-mode nonlinear systems, driven by a periodic pulse sequence and described by the nonlinear Schrödinger equation (NLSE).To analyse these settings, we introduce a parent quantum many-body Hamiltonian, Ĥ0 + V (t), which describes two species of interacting bosons driven by a periodic pulse sequence.From this, we derive an effective quantum Hamiltonian Ĥeff in the high-frequency limit of the drive (ω → ∞).We then derive the effective nonlinear Schrödinger equation upon taking the classical limit, N → ∞, where N is the number of bosons, hence revealing the effective nonlinearities generated by the driving pulse sequence.

Figure 2 .
Figure 2. Two possible realizations in optics: (a) Two modes in an optical ring cavity (1 and 2), repeatedly undergoing mixing operations (⊕ and ) along the ring.These operations correspond to a coupling between the two polarization eigenmodes of the cavity, as realized by means of quarter-wave plates; see Eqs. (2)-(3).(b) Two optical waveguides (1 and 2) with modulated inter-waveguide separation, realizing a "time-periodic" linear coupling between the two optical modes [Eq.(4)].In both cases, the "time" coordinate corresponds to the propagation direction[37,85].

Figure 4 .
Figure 4. Processes in the Hamiltonian in Eq. (8): (a) Intra-mode (Hubbard) interactions; (b) inter-mode (cross) interactions; and (c) single-particle hopping processes.(d) The effective Hamiltonian in Eq. (30) includes pair-hopping processes, by which two interacting particles in the same mode simultaneously change mode.In this illustration, the two modes 1 and 2 correspond to the low-energy orbitals of a double-well potential, and the bosons are represented by green atoms.

Figure 5 .
Figure 5. Population imbalance z as a function of time, as obtained from the quantum dynamics of the driven system (blue curve) and the effective-Hamiltonian quantum dynamics (red curve) for: (a) N = 10 bosons and (b) N = 50 bosons.For each case, the full time dynamics of the driven system is generated using the sequence in Eq. (18) with a period T = 0.2, T = 0.1 and T = 0.05.Here the interaction parameter is set to g = χN = 5 and the static linear coupling is set to Ω0 = 0; the initial coherent spin state |N, θ, ϕ corresponds to z = cos θ = 0.4 and ϕ = 2.25.In all plots, the time-evolved state is evaluated at stroboscopic times tn = T ×n.

Figure 6 .
Figure 6.Energy landscape associated with the classical Hamiltonian Heff(z, ϕ) displayed in Eq. (32), for Ω0 = 0.A few trajectories are indicated as thin blue curves (equipotential lines of the energy landscape).

Figure 7 .
Figure 7. Population imbalance z as a function of time, as obtained from the effective-Hamiltonian quantum dynamics (red curve) and the effective classical equations of motion (blue curve).The number of bosons is: (a) N = 5; (b) N = 10; (c) N = 80; (d) N = 170.Here the interaction parameter is set to g = χN = 5, while the static linear coupling is set to Ω0 = 0; the initial coherent spin state |N, θ, ϕ corresponds to z = cos θ = 0 and ϕ = 2.7; the same initial condition is chosen for the effective classical dynamics.

Figure 8 .
Figure 8. Time-evolving Husimi function Q(z, ϕ; t) for a state |ψ(t) that evolves according to the effective Hamiltonian Ĥeff in Eq.(30).Here, the number of bosons is N = 80, and the other parameters are the same as in Fig.7(c).The initial coherent spin state |ψ(0) = |N, θ, ϕ , at z = cos θ = 0 and ϕ = 2.7, becomes substantially squeezed around time t ≈ 3, hence signaling the breakdown of its classical description.An oversqueezed state, exhibiting Majorana stars, appears around t ≈ 12.The trajectory predicted by the effective classical equations of motion [Eq.(31)] is depicted in white.

Figure 9 .
Figure 9.The driven nonlinear Schrödinger equation versus its effective description: (a) Population imbalance z(t) as a function of time, as obtained from the driven nonlinear Schrödinger equation in Eq. (36) (blue curve) and from the effective classical equations of motion in Eq. (31) (red curve).(b) Zoom in the panel (a): the blue dots highlight the stroboscopic dynamics at times tn; note the micromotion at arbitrary times t = tn.(c) Stroboscopic dynamics z(tn) obtained from the driven nonlinear Schrödinger equation (blue curve and dots), compared with the effective classical description (red curve).(d) Same as in panel (c) but for the other canonical variable ϕ.In all panels, the period of the drive is set to T = 0.1 and the pulse duration to τ = T /20; the interaction parameter is set to g = χN = 5, while the static linear coupling is set to Ω0 = 0; the initial condition corresponds to z = cos θ = 0 and ϕ = 2.7 as in Fig. 7.

Figure 10 .
Figure 10.Energy landscape associated with the effective classical Hamiltonian Heff(z, ϕ; α) displayed in Eq. (41), for four values of the imbalance parameter: (a) α = 0, (b) α = 0.25, (c) α = 0.5 and (d) α = 0.5.A few trajectories are indicated as thin blue curves (equipotential lines of the energy landscape) for each case.The static linear coupling is set to Ω0 = 0. Note the emergence and disappearance of stable fixed points on the Bloch-Poincaré sphere, as the imbalance parameter α is varied.