Efficient and accurate methods for solving the time-dependent spin-1 Gross-Pitaevskii equation

We develop a numerical method for solving the spin-1 Gross-Pitaevskii equation. The basis of our work is a two-way splitting of the spin-1 evolution equation that leads to two exactly solvable flows. We use this to implement a second-order and a fourth-order symplectic integration method. These are the first fully symplectic methods for evolving spin-1 condensates. We develop two nontrivial numerical tests to compare our methods against two other approaches.

Figure 1

The components of the spin density are shown for a continuous-wave solution with parameters: $k_{+}=5, k_{-}=3, {\theta }_{+}=0, {\theta }_{-}=0, A_{+}=3, A_{-}=1$ with $c_0=10, c_1=1$, and $q=0.5$. The remaining parameters are determined by the relations given in Ref. [18].

Figure 2

Relative errors in the conserved quantities (a) normalization, (b) magnetization, and (c) energy for the propagation of the continuous wave solution. The various algorithms, as labeled in (b), all use the same time step ($\tau =3.125\times {10}^{-4}$) and mesh ($M=256, L=2\pi$). The initial condition and other simulation parameters are given in Fig. 1. The dashed (dotted) line is a guide to the eye to indicate linear (square root) scaling with time $t$.

Figure 3

Global error in propagating the continuous wave solution. (a) Global error during a propagation using a time step of $\tau =6.25\times {10}^{-4}$. The dashed line is a guide to the eye showing linear scaling with time $t$. (b) Scaling of the global error after propagation to time $t_f=0.1$ as a function of the time step $\tau$. The line types for the various algorithms are indicated in (b), as is the estimated stable time step size ($t_{\mathrm{stab}}=3.8\times {10}^{-4}$). The dashed (dotted) lines in (b) indicate fourth-order (second-order) power-law scaling with respect to $\tau$. Other parameters for the solution are as indicated in Fig. 2.

Figure 4

Quasisoliton collision: [(a)–(d)] Snapshots of dynamics showing component densities during the collision dynamics at a range of times. We define the collision time as $t_{\mathrm{col}}=19$. These results were calculated using S4 with $\tau =0.01$ and match the results presented in Fig. 1 of Ref. [20]. Other parameters: $L=384, M=2048, x_0=1, c_1=0.314, q=0, \eta =3.091, \xi =1.54, \mu =2$.

Figure 5

Quasisoliton evolution: Relative errors in nonzero conserved quantities (a) $N$ and (c) $E$. Panel (b) shows the magnitude of $M_z$ that develops in the W2 scheme, noting that $|M_z|$ remains zero for the other methods. Dashed (dotted) lines are guides to the eye showing linear (square root) scaling with respect to time $t$. Simulations are for the scenario (initial condition, time, and mesh parameters) given in Fig. 4.

Figure 6

Global error in propagating the quasisolitons using reference S4 solution with $\tau =5\times {10}^{-4}$. (a) Global error during a propagation using a time step of $\tau =0.02$. The dashed line is a guide to the eye that scales linearly with time $t$, while vertical dash-dotted line shows the time of collision of the quasisolitons ($t_{\mathrm{col}}=19$). (b) Scaling of the global error with time step $\tau$ of the solitons propagated to $t=2$. The dashed (dotted) line indicates fourth-order (second-order) power-law scaling with respect to $\tau$, while the vertical dash-dotted line shows the estimated stable time step size ($t_{\mathrm{stab}}=0.022)$. Simulation parameters and the solution mesh are as discussed in Fig. 4.