Spin-wave growth via Shapiro resonances in a spinor Bose-Einstein condensate

We theoretically study the resonant phenomenon in a spin-1 Bose-Einstein condensate periodically driven by a quadratic Zeeman coupling. This phenomenon is closely related to the Shapiro steps in superconducting Josephson junctions, and the previous experimental work [Evrard $et al.,$ Phys. Rev. A 100, 023604 (2019)] for a spin-1 bosonic system observed the resonant dynamics and then called it Shapiro resonance. In this work, using the spin-1 Gross-Pitaevskii equation, we study the Shapiro resonance beyond the single-mode approximation used in the previous work, which assumes that all components of the spinor wavefunction have the same spatial configuration. Considering resonant dynamics starting from a polar state, we analytically calculate the Floquet-Lyapunov exponents featuring an onset of the resonance under a linear analysis and find that spin waves with finite wavenumbers can be excited. This kind of non-uniform excitation cannot be described by the single-mode approximation. Furthermore, to study the long-time resonant dynamics beyond the linear analysis, we numerically solve the one-dimensional spin-1 Gross-Pitaevskii equation, finding that the nonresonant hydrodynamic variables also grow at wavelengths of even multiples of the resonant one due to the nonlinear effect.


I. INTRODUCTION
The engineering of quantum systems by periodic driving has drawn great attention over a decade, and ultracold atoms have become a promising platform for realizing such driven quantum systems due to their high experimental controllability [1][2][3]. Indeed, applying various periodic modulations to ultracold atoms, recent experiments have realized several topological models such as the Haldane model and the Hofstadter-Harper model [4][5][6], and have also observed exotic phases of matters such as a time crystal [7,8].
Such engineering by external driving is recently utilized to generate quantum entanglement in a spin-1 Bose-Einstein condensate (BEC) [9], which is comprised of spin-1 bosons characterized by the three magnetic sublevels m = 1, 0, and −1 [10,11]. Modulating a quadratic Zeeman (QZ) coupling by microwaves, Evrard et al. [9] induce resonance between the different magnetic sublevels and observe entangled spin states. This resonant phenomenon is essentially the same as the Shapiro steps originally investigated in Josephson junctions between superconductors [12] and thus is called Shapiro resonance. In the series of these works [9,13], they have investigated the Shapiro resonance mainly in a restricted situation where a single-mode approximation is valid. Under this approximation, all three magnetic sublevels have the same spatial configuration, and thus we cannot address spin-wave excitation with a finite wavelength as well as nontrivial magnetic pattern formations [14,15].
In this work, we theoretically study the Shapiro resonance in a spin-1 BEC beyond the single-mode approximation. Using the Gross-Pitaevskii equation (GPE) [16,17], we investigate the resonant dynamics starting from a polar state by periodically modulating the QZ coupling. First, we analytically derive the resonant condition by linearizing the spin hydrodynamic equations equivalent to the GPE [18] and subsequently applying the Floquet's theorem to the linearized equations. Here, the resonant condition is identified by a nonzero Floquet-Lyapunov (FL) exponent, which captures an onset of the resonance. Second, we numerically solve the one-dimensional (1D) GPE, demonstrating the validity of our linear analysis and investigating the nonlinear dynamics over a long time. Spin-wave excitations due to the nonlinear effects are explained from the constraints on the hydrodynamic variables. Finally, we discuss the experimental possibilities by using the parameters used in the previous experiments [19][20][21][22][23].
The rest of this article is organized as follows. In Sec. II, we introduce the GPE and the spin hydrodynamic equations for a spin-1 spinor BEC. In Sec. III, we derive the analytical expressions for the resonant conditions by a linear analysis of the spin hydrodynamic equations. In Sec. IV, we show the numerical results of the 1D GPE, showing that our analytical results work well and investigating the long-time resonant dynamics. In Sec. V, we discuss the experimental possibilities for observing nonuniform Shapiro resonances. The summary for this work is given in Sec. VI.

II. MODELS
In this work, we study the dynamics of a spin-1 spinor BEC within the mean-field approximation. Here, we introduce the mean field equation, namely the GPE, and subsequently explain the spin hydrodynamic form [18], which is quite convenient for analytically investigating the Shapiro resonance in Secs. III and IV.
A. Spin-1 Gross-Pitaevskii equation We consider a spin-1 BEC in a uniform system without a trapping potential. Under the mean-field approximation, the macroscopic wave functions ψ m of atoms in the magnetic sublevel m = 1, 0, −1 obey the following GPE: where M is the atomic mass, and c 0 and c 1 are the strength of the spin-independent and spin-dependent interaction, respectively. The sign of c 1 determines the magnetism of the system: The condensate is antiferromagnetic (AFM) for c 1 > 0 and ferromagnetic (FM) for c 1 < 0. The total particle number density ρ and the spin density F µ (µ = x, y, z) are defined by wheref µ (µ = x, y, z) is the µ-component of the spin-1 matrix given bŷ We assume that the QZ coupling strength q(t) consists of a static term and an oscillating term with frequency Ω as given by where q 0 and q osc are the strength of the static and oscillating QZ couplings. Here, we eliminate the linear Zeeman coupling without loss of generality because it can be removed when we move onto the rotating frame of reference in spin space with the Larmor frequency [10].

B. Spin hydrodynamic equations
We can rewrite GPE (1) without any approximation to the equations of motion for the particle number density ρ, the particle density current v µ , the spin density vector f µ , and the nematic density tensor n µν (µ, ν = x, y, z), which are defined by Eq. (2) and where ζ m ≡ ψ m / √ ρ is the normalized spinor wavefunction. By taking the time derivative of the above quantities and using GPE (1), we obtain the following spin hydrodynamic equations [18]: Here, we have defined the spin current v (s) µ and the ne- Note that the number of variables in the spin hydrodynamic equations is larger than in the three-component GPE (1). This is because some of the variables in the spin hydrodynamic equations are dependent on each other via the following constraints [18]: ν=x,y,z The detail is described in Ref. [18]. The previous works on Shapiro resonance [9,13,24] analyze the Madelung form of GPE (1), which is the equations of motion for the density and phase of each component (see Appendix A), under the single-mode approximation. In this paper, we use the spin hydrodynamic equations rather than the Madelung form since the linear analysis beyond the single-mode approximation becomes simpler for the former case as discussed in Sec. III.

C. Parameter setup
In this study, we prepare a polar BEC, where all atoms are spatially uniform and condensed in the m = 0 state, and investigate how the numbers of the atoms in the m = ±1 state increase via parametric resonance. We therefore choose parameters such that the polar state is stable when the oscillating frequency Ω is off-resonance. This condition is satisfied when q(t) moves in the polar phase region: whereρ is the mean particle number density of the condensate. In the following calculations, we choose q 0 , q osc > 0 and q 0 − q osc > max(0, −2c 1ρ ). We also assume Ω > 0 without loss of generality. The bulk chemical potential in the polar state is given by c 0ρ . The corresponding length and time are given by which we use as the characteristic scales of the system.

III. LINEAR ANALYSIS
In this section, we apply linear analysis to the spin hydrodynamic equations around the polar state. We employ the Floquet's theorem to obtain the resonant condi-tions and the FL exponents analytically.

A. Linearized spin hydrodynamic equations
We discuss the linear stability of the polar state, (ψ 1 , ψ 0 , ψ −1 ) = (0, √ρ , 0). In the absence of fluctuations, the particle current, the spin density vector, and the nematic density tensor for the polar state are given bȳ We introduce small fluctuations, δρ, δf µ , δn µν , and δv µ , and write the hydrodynamic variables as By substituting them into Eqs. (8)- (11) and expanding the equations up to the first order in the fluctuations, we obtain 13 linearized equations. Among them, the equations for δf x , δf y , δn xz , and δn yz include the driving QZ term and are divided into two sets of coupled equations: where (A, B) = (δf x , δn yz ) and (δf y , δn xz ). Here, the plus and minus signs on the right-hand side of Eq. (27) are for (A, B) = (δf x , δn yz ) and (δf y , δn xz ), respectively. These equations indicate that the transverse spin components, δf x and δf y , and the off-diagonal elements of the nematic density tensors, δn yz and δn xz , can be amplified by the driving QZ term. Thus, we call them the resonant variables. Note that Eq. (27) is valid only for the initial stages of the resonant dynamics since we have neglected the nonlinear terms.
B. Application of the Floquet's theorem to Eq. (27) Equation (27) can be transformed into an eigenvalue problem of an infinite-dimensional matrix by utilizing the Floquet's theorem, which enables us to derive the reso-nant conditions. First, we introduce the Fourier transform and rewrite Eq. (27) as where k = 2 k 2 /(2M ). Second, we perform the linear transformation such that the time-independent part of the matrix in Eq. (30) is diagonalized. The linear transformation is given by where ±iE k are the eigenvalues of the matrix in Eq. (30) with q osc = 0: Then, Eq. (30) reduces to Note that E k is identical to the spin-wave excitation energy obtained by the Bogoliubov analysis under the static QZ term. Equation (33) indeed reproduces the equation of motion for the spin wave when q osc = 0.
Because Eq. (33) has the periodicity of T = 2π/Ω, we can apply the Floquet's theorem: Without loss of generality, the functions S A and S B can be expanded as where λ F is the Floquet exponent. The summation over j ∈ Z comes from the Fourier expansion of T -periodic functions, and C j,A and C j,B are the Fourier components. By substituting Eq. (36) into Eq. (33), we obtain whereQ and |C are the infinite-dimensional matrix and the column vector, respectively, given bŷ Here, we omit the sign ± from Eq. (37) because both equations provide the same result about the real part of the Floquet exponents (see Appendix B). The FL exponent is defined as the real part of λ F : If there exists a positive λ FL , S A and S B exhibit exponential increase, that is, the Shapiro resonance. Therefore, we refer to the positive ones as the FL exponents in the following discussions unless otherwise noted. In the following sections, we solve the infinitedimensional eigenvalue problem of Eq. (37) by assuming b 1. This assumption is satisfied for large E k . Furthermore, with the parameters discussed in Sec. II C, the off-diagonal elements ofQ is proved to be smaller than E k /2, i.e., b(Ĝ) l,m < E k /2 with l, m = 1, 2, which also supports the validity of the perturbative expansion below. The proof is given in Appendix C.

C. Resonant conditions
We start from the simplest case of b = 0, at which the matrix (38) becomes diagonal and has the eigenvalues Note that all λ F 's are pure imaginary in this case, which means that the amplitudes of the spin density vectors and nematic density tensors do not grow in time.
A small finite b couples the above eigenmodes, and its effect becomes prominent when two eigenvalues are close to each other. Note thatQ in Eq. (37) is a pseudo-Hermitian matrix, i.e., there exists an Hermitian matrix η such thatQ † = η †Q η [25]; In the present case, η is given by η = Diag [· · · , 1, −1, 1, −1, · · · ]. In such a case, two close eigenvalues can coalesce and become a pair of complex conjugate values. Thus, λ F 's have nonzero real parts around the intersections of ω j,± (j = 0, ±1, ±2, · · · ) as a function of Ω [see Fig. 1(a)]. By solving ω j1,+ = ω j2,− with j 1 , j 2 ∈ Z, we obtain the intersection points at where J ≡ j 2 − j 1 . Since Ω is assumed to be positive, we restrict J to be a positive integer. Equation (44) is the approximated condition for the Shapiro resonance with spatial degrees of freedom.
To confirm the above argument, we numerically calculate λ FL from Eq. (33) by following the method used in Ref. [26]: For a fixed wavenumber, we numerically solve Eq. (33) from t = 0 to T with the initial conditions (S A , S B ) T = (1, 0) T and (0, 1) T , obtaining S 1 (T ) and is the numerically obtained one shown in Fig. 1(b). The parameters are the same as those in Fig. 1(b). For a better visibility, we color the curves according to the value of kξ. We can see that the FDMA works well, in particular, at the large wavenumbers.
S 2 (T ), respectively; The 2×2 matrixÛ = (S 1 (T ), S 2 (T )) is the one-period time-evolution operator and its eigenvalue λ U is related to the Floquet exponent as λ U = e λFT ; Thus, we obtain the FL exponent from the numerically obtained λ U as λ FL = Re [(ln λ U )/T ]. The color plot in Fig. 1(b) shows the numerical result of λ FL , where we obtain two FL exponents with opposite signs and plot only the positive one in Fig. 1(b). We also plot Eq. (44) with dotted curves in the same figure, which show an excellent agreement with the numerical result.
Note that although Eq. (37) has an infinite number of eigenvalues, we obtain only a pair of FL exponents in the numerical procedure, suggesting that the imaginary part of the eigenvalues ofQ can only take two values corresponding to the numerically obtained λ FL . Indeed, we will see below (and also in Appendix B) that the detailed analysis of the infinite-dimensional matrixQ results in a pair of FL exponents with opposite signs.

D. Finite-dimensional matrix approximation
We have obtained the resonant condition of Eq. (44), but our numerical results in Fig. 1(b) show that there is a width of the resonant frequency. In this section, we analytically calculate the FL exponent and the width of the resonance by approximating the infinite-dimensional matrix in Eq. (37) with a finite one [27,28].
Our procedure is as follows. Let |C j,± denote the normalized eigenmode at b = 0 with eigenvalue ω j,± . We consider the case when the oscillating frequency is close to the Jth resonance, Ω 2E k /J, at which |C j,+ and |C j+J,− are almost degenerate. In such a case, we can neglect the other modes, andQ is approximated by the 2×2 matrix a|Q |a with a, a = C j,+ and C j+J,− . This procedure is similar to what we do in calculating the energy bands of the nearly free electron model, where a band gap opens at the boundary of the Brillouin zone. However, the 2 × 2 matrix is not enough for J > 1 because the off-diagonal elements vanish for J > 1; In the perturbative expansion ofQ in powers of b, the coupling between |C j,+ and |C j+J,− first appears in the Jth order term. We therefore need to take into account the intermediate states that appear in the coupling between |C j,+ and |C j+J,− and approximateQ with the one projected onto the restricted Hilbert space. Below, we demonstrate the cases of J = 1 and 2.

J = 1
We consider the coupling between |C j,+ and |C j+1,− which is the first order in b. We approximate the matrix Q with the 2 × 2 matrix: from which we obtain the FL exponent The FL exponent takes a nonzero value when the frequency Ω is in the resonant region: Note that although j takes an arbitrary integer, the j dependence of the eigenvalue appears only in the real part of iλ F . That is, all the intersection points in Fig. 1(a) at Ω = 2E k give different Im[λ F ]'s but the same FL exponent. This is consistent with the fact that we have only two FL exponent in the numerical calculation.
To see the validity of the analytical result (46), we compare λ is the numerically obtained FL exponent shown in Fig. 1(b). One can see that the analytical result well reproduces the width and the amplitude of the resonance obtained by the numerical calculation.

J = 2
In the case of J = 2, we consider the coupling between |C j−1,+ and |C j+1,− . These states couple with each other via the intermediate state |C j,± in the second order of b. Note that in the same order of b, the couplings with |C j−2,± and |C j+2,± states shift the eigenvalues associated with |C j−1,+ and |C j+1,− , respectively. Thus, we need to solve the following 8 × 8 matrix eigenvalue problem: The analytical expression of the eigenvalue is so complicated that we show the detail in Appendix D. We compare the obtained analytical solution λ (ana) FL with the numerical one in Fig. 2(b). As in the case of J = 1, one can see that λ  Although we can calculate λ FL for higher J in a similar manner as in the cases of J = 1 and 2, it requires tedious calculations. We note that λ FL at Ω = 2E k /J can be obtained more easily by using the conventional perturbation theory for degenerate states. We regard b as a small perturbation parameter and rewrite Eq. (37) in the following form: whereÊ and bV are the matrices composed of the diagonal and off-diagonal elements ofQ, respectively. As we have discussed in the above, at Ω = 2E k /J, the states |C j,+ and |C j+J,− are degenerate at b = 0. These states are coupled in the Jth order perturbation, giving rise to a nonzero λ FL . Following the conventional perturbation theory for degenerate states, λ FL 's for J = 1, 2 and 3 are obtained as the imaginary parts of the eigenvalues of the 2 × 2 matrices, respectively, given by where l, l = (j, +) and (j + J, −), and we take the summation over l 1 , l 2 = (j, s) for all possible combinations of j ∈ Z and s = + and − other than (j, +) and (j + J, −).
The resulting FL exponents are given by The results for J = 1 and 2 coincide with Eqs. (46) and (D1), respectively. We can calculate the FL exponent for J ≥ 4 in the same manner, but we need to be careful that the diagonal terms of V (J ) ll with even J < J generally lifts the degeneracy of (j, +) and (j + J, −) states, which also contributes to Imλ F for J = 3. Figure 3 compares numerical results of λ FL (dotted lines) with Eqs. (53)-(55) (solid lines). As one can see, the approximation works quite well for larger kξ because the perturbation parameter b = q osc /2E k becomes smaller for larger kξ.

IV. NONLINEAR RESONANT DYNAMICS
We numerically solve the spin-1 spinor GPE (1) in a one-dimensioinal system and investigate dynamics of the Shapiro resonance beyond the linear analysis. We first confirm the validity of the linear analysis in Sec. III. We then study the nonlinear dynamics using the hydrodynamic variables and discuss how the long-time dynamics proceeds.
A. Numerical results using the GPE We solve Eq. (1) for a 1D system of a system size L x /ξ = 256 with a periodic boundary condition. The initial wavefunction is given by where R m,l (x) (l = 1, 2) is a uniform random number in [0,1]. We choose the parameters as (q 0 , q osc , Ω) = (0.5, 0.39, 2E kres ) with k res = 10 × 2π/L x . According to Eq. (46), the Ω is the resonance frequency of J = 1, and the spin waves with the wavelength L x /10 will grow. The interaction parameters are set to be c 1 /c 0 = 1/20 and −1/20 for AFM and FM BECs, respectively. In Fig. 4, we plot the time evolution of the spin density f x (x, t) which is expected to grow as e λFLt within the linear analysis. Figures 4(a) and (b) show the spatial and temporal distributions of f x (x, t) for the AFM and FM systems, respectively. As expected from the linear analysis, the spin waves with wavelength L x /10 grow in the early stage (see the enlarged views in Fig. 4). In the late stage (t/τ 1500), f x (x, t) in the AFM system and −1/20, respectively. The parameters for the QZ terms are (q0, qosc, Ω) = (0.5c0ρ, 0.39c0ρ, 2E kres ). The insets show the enlarged view of the early and late stages. The spin density in the AFM system tends to be finer structure than that in the FM one. The growing speed of fx in the AFM system is slightly slower than that in the FM one because λFL of the former is smaller than that of the latter.
exhibits finer structure than that in the FM system. This can be seen more clearly in the k-space. Figure 5 shows the dynamics in the k-space. Here, we define the Fourier components in the 1D system asÃ(k) ≡ Lx 0 A(x)e −ikx dx and plot |f x (k, t)| 2 averaged over 1000 samples of random initial states given by Eq. (56). Figures 5(a1) and (a2) show the time evolution of |f x (k res , t)| 2 for the AFM and FM systems, respectively, where the dashed lines are the analytical result e 2λFLt with λ FL given in Eq. (53). One can clearly see that, for both cases, |f x (k res , t)| 2 in the early stage agrees well with the linear analysis.
The linear analysis breaks down in the long-time dynamics, and |f x (k res , t)| 2 in Figs. 5 (a1) and (a2) deviates from the exponential growth. The breakdown of the linear analysis is also seen in the k dependence of |f x (k, t)| 2 shown in Fig. 5(b), where the results for AFM (FM) systems are shown with open (filled) circles. As one can see, in addition to the main peak at k = k res , small peaks appear around odd multiples of k res at t/τ = 400. They are consequences of the nonlinear effect as mentioned in the next subsection. As time evolves, these peaks become broader, and |f x (k, t)| 2 eventually distributes in the wide range of k. Figure 5(b) directly shows that the sign of c 1 significantly alters the late dynamics especially at t/τ ≥ 680: The spectra of the AFM system have the larger values in the high wavenumber region than those of the FM system. This is consistent with the real space configuration of f x (x, t) in Fig. 4.

B. Excitations beyond linear analysis
We discuss nonlinear effects not captured under the linear analysis by taking nonlinear terms on the hydrodynamic variables. Here, we first focus on f z and n zz which are never resonated under the linear analysis.
from which we find that δf z and δñ zz have peaks at k = 0 and ±2k res because the resonant variables grow around k = ±k res . Note thatf z is almost zero at k = 0 because of the conservation law of the total longitudinal magnetization f z (r, t)d d r which is a small value coming from the initial noise.
The spectra of other nonresonant variablesñ xx ,ñ yy , andñ xy can be understood from directly approximating Eq. (10) in the same manner as the above: ∂ ∂t δn xy c 1ρ (δf y δn xz − δf x δn yz ), which implies that these variables grow around k = 0 and ±2k res .
We systematically investigate the peak structures of all the spin and nematic variables. Our numerical calculations in Fig. 6 demonstrate that the resonant variables (f x ,f y ,ñ xz , andñ yz ) have peaks around odd multiples of k res and the others (f z ,ñ xx ,ñ yy ,ñ zz , andñ xy ) have peaks around zero and even multiples of k res in the early stage of the nonlinear dynamics. This is consistent with the results predicted by Eqs. (59)-(65). These spectra of the numerical results are similar to the typical pumped spectra seen in classical fluid [29]. By taking the fluctuation terms to higher orders, we expect that all the spectra in the early stage in the figure are explained.
Using these results, we also understand the time evolution of ρ ±1 , which is related to n zz and f z via Thus,ρ ±1 grows around zero and even multiples of k res . This means that the increase in the particle number of the m = ±1 component, which characterize the onset of the Shapiro resonances in the previous experimental works [13,24], is attributed to the nonlinear effect.   [19][20][21][22][23] for the J = 1 Shapiro resonance. In column c1ρ(0), the values marked with TF are calculated using the TF approximation, and the others are the values specified in the literature. In the evaluation of Ω, ∆Ω, and τFL, we use the effective spin-dependent interaction energy c 1ρ instead of c1ρ, where c 1ρ is related to c1ρ(0) via Eqs. (E2) and (E4).

V. DISCUSSION
We discuss experimental possibility for observing spinwave excitations due to the Shapiro resonance on the basis of our linear analysis of Sec. III. As shown in Sec. III D, the J = 1 mode has a larger FL exponent than those with higher J. This means that the nonuniform Shapiro resonance with J = 1 is experimentally more accessible. We, therefore, address the resonance with J = 1 by using realistic experimental parameters in what follows.
To consider the experimental possibility, we introduce two important quantities: One is an inverse of the maximum FL exponent τ FL ≡ (λ FL | max ) −1 , and the other is the width ∆Ω for the resonant frequency. Here, τ FL is the characteristic time scale for the growth of spin waves, thus should be sufficiently longer than the time resolu-tion of the experiments and shorter than a lifetime of a BEC. On the other hand, ∆Ω should be narrow enough to specify the wavenumber of the resonant spin waves, but not too narrow so that it is easy to adjust in experiments. Here, we show two constraints for τ FL and ∆Ω derived from Eqs. (46) and (47): The inequality (67) indicates that the lower bound of τ FL is the characteristic time scale of spin dynamics, which is a few to a few tens of milliseconds in typical experiments. The constraint (68) is a trade-off relation between ∆Ω and τ FL . We have used these constraints to choose the experimental parameters in the following discussion.
Using the obtained analytical results of Eqs. (46) and (47), we evaluate τ FL at Ω = 2E kres and ∆Ω for the parameters in the experiments [19][20][21][22][23]. We fix q 0 to be the QZ energy at the magnetic field of 700 mG and choose q osc such that the polar state becomes stable for off-resonant Ω. The values of the hyperfine splitting energy and the interaction strengths c 0 and c 1 are given in Refs. [10,23,30,31]. As for the resonant wavenumber k res , which determines Ω = 2E kres , we choose k res = 10π/R x for a quasi-1D BEC and k res = (10π/R x , 10π/R y ) for a quasi-two-dimensional (quasi-2D) BEC, where R x and R y are the largest and the second-largest Thomas-Fermi (TF) radii, respectively. Note that, in low-dimensional BECs, the spin interaction energy c 1ρ in Eqs. (46) and (47) and inequalities (67) and (68) is replaced with the effective one, c 1ρ , as derived in Appendix E. Here, c 1ρ is related to the peak density ρ(0) of the 3D distribution via Eqs. (E2) and (E4) for quasi-1D and quasi-2D BECs, respectively.
We summarize the estimated values of τ FL and ∆Ω in Table I, together with the values of c 1 ρ(0), q 0 , q osc , k res , and Ω. The obtained sets of values can be experimentally accessible.

VI. CONCLUSION
In this work, considering the polar state in the uniform spin-1 BECs, we theoretically investigated the Shapiro resonance driven by a periodic forcing of the QZ term. Unlike the previous works which discuss the Shapiro resonance in a strongly confined BEC without spatial degrees of freedom [13,24], we took into account the spatial dependence of the condensate and investigated the growth of spin-wave excitations.
First, we studied the Shapiro resonance within the linear analysis by employing the spin hydrodynamic equations equivalent to the GPE. Applying the Floquet's theorem to the equations, we analytically obtained the real parts of the Floquet exponents, i.e., the FL exponents, featuring the growth rate of the resonant spin and ne-matic variables by using the two kinds of approximations: the finite-dimensional matrix approximation and the degenerate perturbation approximation. From these results, we identified the resonant conditions and found the spin-wave excitation with finite wave-numbers, which cannot be described by the single-mode approximation.
Second, we numerically solved the spin-1 GPE and study the validity of the linear analysis and the nonlinear dynamics in the late stage. In the k-space time evolution of the hydrodynamic variables, we confirmed that the resonant variables are excited in the early stage of the dynamics as expected in the linear analysis. However, as time goes by, we found the emergence of spin-wave excitations at wavenumbers of integer multiples of the resonant one, which cannot be predicted by our linear analysis. We explained this nonlinear effect by using the constraints on the hydrodynamic variables and the spin hydrodynamic equations. Continuing to drive the QZ term, we investigated the long-time dynamics and observed that, in the AFM spinor BEC, the finer spin distributions emerged compared with the FM spinor BEC.
In the final section, we discussed the experimental possibilities for observing spin-wave excitations due to the Shapiro resonance on the basis of our linear analysis. Using the parameters used in the previous experiments [19][20][21][22][23], we showed that a driving QZ field at a frequency in the order of 100 Hz under a bias field of 700 mG can induce the growth of spin waves at the experimentally accessible length and time scales. Our results give an experimental procedure to excite spin waves of a specific wavenumber selectively, which would be useful for future studies of spin dynamics. of the GPE: where ∆θ = θ 1 + θ −1 − 2θ 0 and These equations of motion elucidate that the phase difference ∆θ induces changes in the number fraction in each magnetic sublevel. This mechanism is similar to the Josephson effect. Since the left-hand side of Eq. (A1) includes the number densities of all components, the particle flow between spin components vanishes at the place where at least one of the densities ρ 0,±1 becomes zero. One can also see from the Madelung form that a spatially uniform linear Zeeman term p does not affect the dynamics of ρ m .
The equations of motion (A1) and (A2) certify the following identity: where N m ≡ ρ m d d r. This is due to the conservation of the longitudinal magnetization, F z d d r. Therefore, no resonance occurs when we start from a fully polarized state along +z or −z direction.
That is, ε * is also an eigenvalue ofQ, and the corresponding eigenstate is given by η|C . Since (iλ F ) is an eigenvalue ofQ [see Eq. (37)], (iλ F ) * = i(−λ * F ) is also an eigenvalue ofQ. It follows that if there is a nonzero FL exponent λ FL = Re[λ F ], there is always another nonzero FL exponent −λ FL = Re[−λ * F ].
Appendix C: Validity of the finite-dimensional matrix approximation and the perturbation approximation The off-diagonal elements of the matrixQ in Eq. (38) is proved to be less than the half of the diagonal ones under the polar regime given by Eq. (17): b(Ĝ) l,m E k < q osc ( k + q 0 + c 1ρ ) 2E 2 k = q osc ( k + q 0 + c 1ρ ) 2( k + q 0 )( k + q 0 + 2c 1ρ ) <      q osc ( k + q osc + |c 1 |ρ) 2( k + q osc )( k + q osc + 2|c 1 |ρ) (c 1 < 0) q osc 2( k + q 0 ) (c 1 > 0) Here, we use the condition (17) from the second to the third line. This confirms the validity of the finitedimensional matrix approximation and the perturbation approximation in this paper.
Here, we define α = k + q 0 , β = q osc , and γ = c 1ρ . The above results were obtained using Mathematica. In Fig. 2(b), the curves labeled with λ (ana) FL represents the positive imaginary part of δλ F in Eq. (D1), showing good agreement with the numerically obtained one. Note that although we have calculated for (j 1 , j 2 ) = (j − 1, j + 1), the imaginary part of Eq. (D1) does not depend on j, suggesting that the FL exponent, i.e., Im(λ F ), is the same for all combinations of (j 1 , j 2 ) satisfying j 2 − j 1 = 2. We can also confirm that the FL exponent at the resonance point Ω = E k coincides with Eq. (54) by substituting Ω = E k into Eq. (D1) and expanding it with respect to b up to the second order.