Home

GPELab A Matlab toolbox for computing stationary solutions

1. n w x Fi W x Fin W x F W x 21 W x 22 W x 2n W x Fy 1 W x Fon x sd Fy n W x setting Fe W x Fem Di k X5 4 1 lt fm lt Ne To be consistent with the one k EO 1 x component case we consider MRo71 2 BF Wo x x The operator AR Myny C is defined by ARengy ARR Y ARGY I 1 Re n EN Re n 1 2 Ap t v z vme 2 y 5 28 5 318 e e ela Y 5 2 RELAXATION PSEUDO SPECTRAL SCHEME RESP 91 The finite dimensional operator Aa Myn C is explicitly given through the matrices I0 0 Va Vi lt Van O I 0 Va Vas PS Van In ee me Myn R V S Myn R 0 0 aai I YN 1 Von SS Vn Ne where the diagonal matrices are defined by Ilem 9 5 1 00 1 Mm R and Ven Vem X5k 51000 1 Mm R The block diagonal matrix AR in 5 28 is implicitly given by the discrete differentiation operators via the FFTs iFFTs ANE Adel ea CL and 0 7 levee E C yl yellow CUP 5 29 For k 1 2 we also define Gh Gh suis Gin GF Sa Gx sis Gon EMO Sha Gky 5 Chew setting Gi Gm e Geos Mmu C The right hand side operator BRen MN CMNe is defined by pBRen y BREN Row Ren In 1 Re n 1 2 m E eRe VAMO w 5 30 1 1 By v 3 Ata G 0 Pla Y For solving the second equati
2. e 77 4 7 Ground state of a 3d Gross Pitaevskii equation with quadratic potential 5 1 1 The Lie ADI TSSP scheme for the rotating GPE 84 5 1 2 The Strang ADI TSSP scheme for the rotating GPE 86 5 1 3 Extension to the multi components case 2 e 87 5 2 Relaxation pseudo SPectral scheme ReSP 5 2 1 A Relaxation pseudo SPectral scheme ReSP for the one component case 5 2 2 Extension of the ReSP scheme to the multi components case 90 5 3 Integration of a stochastic potential 2 0 20 0 200000002 92 DW ue ee Oe ee di Se A a 92 Be ei hae ty se or Ear amp arc e 93 bonas as Gea a aca 93 5 5 GPELab functions for the dynamics e 94 ia An Rowe A ened ee a A ee 94 oe Oe ee E Se A S 95 eS tee heed oe ee ee ee 97 linearity Te ag cog ene peeve Seley gee de reek Gl oy eds a dia Bee es 98 5 6 2 Dynamic of a dark soliton in a Bose Einstein condensate in 2D 100 5 6 3 Dynamic of a rotating Bose Einstein condensate in 2D 2 103 CONTENTS 5 A Copyrights amp credits 107 109 CONTENTS List of Figures 2 1 Two possible gaussian initial data for initializing the iterative scheme with y MEL My 00 eis ae oe a eee ae A Ge me oe A 23 23 2 3 Examples Thomas Fermi approximation for potentials 2 13 left and 2 19 right 2 4 Examples Thomas Fermi approximation for potentials 2 20 left
3. 0 a 80 4 36 Choosing the initial data e 80 4 37 Initializing the outputs setting the informations to print and launching the compu he See ee A ee ee Se 81 9 1 The Method_Var2d function ooa e 94 5 2 An example of initialization and use of the Method_Var2d function 96 Ware Bag ww El be eh ee Ae ee 96 5 4 An example to use the TimePotential_Var2d function 97 fe a ween he E E E ee a 97 5 6 An example to use the StochasticPotential_Var2d function 98 EN N 99 5 8 Setting the Physics1D structure and adding the nonlinear operator 99 al A eee deba o ee hk ee Se 99 agen dun st used Ad 100 shes 100 ed oe ERE eee ean ee 2 eS 100 AAA te ee ee es Ge a ec ee ee eee re ae ee se 101 A 101 ai 102 5 16 Launching the computation of the ground state o 102 rasa losas 103 E ps ple 103 5 19 Building the Method and Geometry2D structures for the simulation 104 mr A 105 a di a oe es ee AS 105 A atau deg E 105 12 LIST OF TABLES Chapter 1 Introduction amp informations 1 1 Overview GPELab Gross Pitaevskii Equation LABoratory is a Matlab toolbox devoted to the numerical computation of stationary and dynamical solutions of the 1d 2d 8d GPE for general nonlinearities including rotation terms and possibly multi components problems This user guide explains which problems GPELab can solve and which algorit
4. 3 7 Launching the simulation setting the outputs and informations on the computation 3 7 1 The OutputsINI_Var2d function Outputs OutputsINI_Var2d Method Evo_outputs save userdef_outputs userdef_outputs_names globaluserdef_outputs globaluserdef_outputs_names Table 3 27 The OutputsINI_Var2d function The OutputsINI_Var2d function initializes the outputs of a simulation by building the Outputs structure Outputs are scalar values computed using each component of the wave function during the simulation In GPELab the predefined outputs are the modulus of the wave function at the center of the domain the root mean square in each direction the energy the chemical potential and the angular momentum More outputs can be computed by using user defined functions The outputs are computed and displayed in the command window at each iteration incremented by the value of the Evo_outputs variable They are also stored after the simulation in the Outputs 54 CHAPTER 3 HOW TO USE GPELAB STATIONARY SOLUTIONS structure see the GPELab2d function Section page 60 The Method structure is a required argument of this function Concerning the optional arguments we have e Evo_outputs N 5 is a variable corresponding to the number of iterations between each computation of the outputs It must be smaller or equal to Evo from the Print_Var2d see Section page e save N 0 is a variable corresponding to the choice of wherever
5. Table 3 17 An example to use the Physics2D_Var2d function 3 6 2 The Potential_Var2d function Physics2D Potential_Var2d Method Physics2D Potential G Table 3 18 The Potential_Var2d function The Potential_Var2d function allows to define the time independent potential operator i e V t x V x in the problem by modifying the Physics2D structure It must be provided with the Method and Physics2D structures The optional arguments are the following e Potential If a function Potential in F M y 1 R M y 1 C is provided the physical potential is defined as follow for each j k 1 Ne Potential z y if j k If Potential is a cell array of functions in CN Ne F Mwy 2 RJ Mn 1 C then the potential is defined by Vir x y Potential j kx y for j k 1 Ne The default argument is quadratic_potential2d which corresponds to 3 1 y if j k Note that in the case of a stationary state computation the potential operator should be time independent e G My n C ones N_c is a complex variable that multiplies the potential element by element leading to the following potential Vj k x y G j k Potential j kj x y for j k 1 Ne 46 CHAPTER 3 HOW TO USE GPELAB STATIONARY SOLUTIONS For example we want to set a quadratic potential for the computation of a ground state for a multi components BEC with internal atomic Josephson junction as in where
6. e the variable coefficients matrices in front of the gradient Cue GG Cg e Ghi x Gh Gy x Gio Gale Chay 00 the diagonal gradient On U t x Oz Y t X 1 1 No the nonlinearity matrix Fui U t x x Fia V t x x Fin W t x x F W t x x _ Ea ERE F V t x x o Fon U tx Pal Dd PEO lt lt EAU 2 6 EXTENSION TO THE MULTI COMPONENTS CASE 31 Moreover we have the following mass normalization constraint Ne Ne Ne N Y Y NW Y f W t x 2dx Y f ly 0 x dx U 1 j l aie T We also introduce the energy E W zoi Pit a 22 a J f R vex where Fenergy is an operator related to the nonlinearity F by the differentiation relation U Fenergy v where 6 designates the G teaux derivative For example in the case of a decoupled cubic nonlin earity Fenergy is already defined in GPELab and is given by 2 50 d V x Gn BFenergy V t x x k vex dx Lar 0 1 0 a t x 0 Fonergy W t x xX 2 za 0 0 b t x For dipolar gazes when a nonlocal integral form of the nonlinearity must be considered then the user must define himself the corresponding function Fenergy in GPELab 2 6 2 Stationary states and the CNGF As in the single component case we consider the problem of finding stationary states for the system 5 14 More specifically we are looking for a solution Y such that U t x e G x whe
7. A V x Boala Bial P Ye K 10 2 Va where V x gt a y is the quadratic potential j are the interactions constants and is the intensity of the Rashba coupling For defining the operator effect we have to set the derivation operators such that the gradient operators are Gi a y a pal G x y C rue In GPELab we thus have to create a cell array of functions and use the Gradientx_Var2d and Gradienty_Var2d functions like in Table to add the Rashba coupling to a system of two Gross Pitaevskii equations function Gradx Example_gradientx Kappa Gradx cel1 2 Gradx 1 1 x y 0 Gradx 1 2 x y 1i Kappa Gradx 2 1 x y 1i Kappa Gradx 2 2 x y 0 end function Grady Example_gradienty Kappa Grady cell 2 Grady 1 1 x y 0 Grady 1 2 x y Kappa Grady 2 1 x y Kappa Grady 2 2 x y 0 end Kappa 1 Physics2D Gradientx_Var2d Method Physics2D Example_gradientx Kappa Physics2D Gradienty_Var2d Method Physics2D Example_gradienty Kappa Table 3 24 An example to use the Gradientx_Var2d m and Gradienty_Var2d functions 3 7 LAUNCHING THE SIMULATION 53 3 6 6 The InitialData_Var2d function Phi_O InitialData_Var2d Method Geometry2D Physics2D InitialData_Choice XO YO gamma_x gamma_y Table 3 25 The InitialData_Var2d function The InitialData_Var2d function builds an initial wave
8. CUP that is assumed to approximate X k CM for a function 14 computes a vector Y CAP such that Y AP Arp Y ARDY aero GL vg ater 2 39 ARE ga A N y 28 CHAPTER 2 COMPUTATION OF STATIONARY STATES FOR THE GPE The evaluation of the two above operators is made as follows For Pd the application is direct since it is realized pointwize in the physical space by setting lia j k V Vix en leo ye 107x510 2 40 for j k Pz The symbol 0 denotes the Dirac delta symbol which is equal to 1 if and only if j k and O otherwise Let us note that the discrete operator is represented by a diagonal matrix after reshaping The label TF refers to the fact that this operator is related to the discretization of the Thomas Fermi approximation On another hand using 2 36 and 2 37 the differential operators in the x or y direction are discretized as V j k Pik J 2 1 1 ane Hell P 7 gt iLpdp YR te Mp aj p J 2 1 K 2 1 ss q K 2 Therefore we obtain the following pseudo spectral approximation of the operator L on the spatial grid VG k PJK MENE i 2AN ye lArl B in 2 41 Another differentiation leads to the second order differential operators in the x or y direction V j k PJK J 2 1 2 1 eee Moz 5 S12 dp yu then est p J 2 1 K 2 1 a AZA y M ap 1 Urta q K 2 which gives the pseudo spectral approxi
9. Table 5 10 Creating a function to locate the center of the soliton We finally set the outputs and the printing informations then we launch the simulation We choose to print informations in the command window every 15 iterations and to draw a figure of the square of the module of the solution Moreover we add the previous function as an output to compute with the name Center of the soliton and we want it to be computed each 10 iterations We do not save the solution during the simulation This is done in Table 5 11 Solution_save 0 Outputs_iterations 10 Output_function 1 Phi X Soliton_center Phi X Output_name 1 Center of the soliton Outputs OutputsINI_Varid Method Outputs_iterations Solution_save Output_function Output _name Printing 1 Evo 15 Draw 1 Print Print_Varid Printing Evo Draw Phi Outputs GPELab1d Phi_0 Method Geometry1D Physics1D Outputs Print Table 5 11 Setting the outputs and the Print structure then launching the simulation At the end of the simulation we can retrieve the informations about the soliton using the Outputs structure For instance we have computed the center of the soliton at each 10 iterations We can print it using the plot function from Matlab This is coded in Table Time 0 0 01 1 plot Time Outputs User_defined_local 1 xlabel Time ylabel Center of the soliton Table 5 12 Plottin
10. 020 92 t x 5 6 pe tn x PO tat x A second step which corresponds in fact to consider 5 3 is to solve the following PDE from time t tn tot tn41 ip t x V t x b Ex BlY t x PW E x 5 7 PO tnx Y tn 1 X l We know that this equation can be solved exactly and the solution is given by YOE x YO tnp 1 YM ataata Si V sds 5 8 which provides an approximation of y l x VO tai x We remark that the ADI technique implies a loss of symmetry of the scheme when solving the partial differential operators of Eq 5 2 Indeed we first integrate in the x direction in 5 5 and next in the y direction according to 5 6 To avoid this problem and symmetrize the scheme we alternate the ordering of the derivative directions any two time steps Most specifically from t to tn 1 we solve 5 5 and next equation 5 6 followed by 5 7 For the next step i e from t 1 to tn 2 we first solve 5 6 and then equation 5 5 and finally again Eq 5 7 Space discretization in 2D and implementation We consider now the notations related to the direct and inverse Fourier series transforms 2 36 2 37 page Fourier transforming Eq 5 5 with respect to the x variable and by using the orthogonality of the Fourier functions we obtain E 1 4 10 t y Gu ump WD t y 1 J 2 lt p lt 3 2 tn St lt tny 86CHAPTER 5 COMPUTATION OF THE DYNAMICS OF GPE METHODS FUNCTIONS EXAMPLES
11. 10 We know that the default nonlinearity is the cubic nonlinearity Therefore we only have to build the Physics1D structure with the desired coefficients and to add the default non linear operator to the physics of the problem We report on Table 5 8 how it is coded Delta 1 Beta 10 Physics1D Physics1iD_Varid Method Delta Beta PhysicsiD Nonlinearity_Varid Method Physics1D Table 5 8 Setting the Physics1D structure and adding the nonlinear operator We then set the initial function to use for the computation We want to simulate a soliton and we would like the initial data to be like 2 Wo x 4 aap ech var exp i ibo To do so we have to use the mesh grid contained in the geometry structure and define the initial data We choose the following parameters a 1 c 5 and 07 0 This is done in Table 5 9 a 1 c 5 theta0 0 X Geometry1D X Phi_011 sqrt 2 a 100 sech sqrt a X exp 1i c X 2 1i theta0 Table 5 9 Building the initial data We want to save the position of the center of the soliton corresponding to the position of its highest module as an output Therefore we define the function that locates the center of the soliton which done in Table 1OOCHAPTER 5 COMPUTATION OF THE DYNAMICS OF GPE METHODS FUNCTIONS EXAMPLES function X_center Soliton_center Phi X Max_Phi I_center max abs Phi X_center X I_center
12. ifft2 1ixffty fft2 phi Grad_x_norm sqrt Geometry2D dx Geometry2D dy sum sum abs Grad_x 2 Grad_y_norm sqrt Geometry2D dx Geometry2D dy sum sum abs Grad_y 2 Grad_norm Grad_x_norm Grad_y_norm end Outputs OutputsINI_Var2d Method 1 phi x y fftx ffty Example_outputs Geometry2D phi x y fftx ffty Table 3 28 An example to use the OutputsINI_Var2d function for a user defined function with a single component However if one wants to compute the root mean square of the sum of two components nus w WPa tox Yalt x Pedy one has to proceed differently because a function computing the root mean square of the sum of two components takes the cell vector of the two wave functions as argument Therefore we have to use globaluserdef_outputs We show how to do this in table function RMS Example_outputs Geometry2D Phi x y fftx ffty RMS_local x 72 y 72 abs Phi 1 Phi 2 2 RMS sqrt Geometry2D dx Geometry2D dy sum sum abs RMS_local 72 end Outputs OutputsINI_Var2d Method 1 Phi x y fftx ffty Example_outputs Geometry2D Phi x y fftx ffty Table 3 29 An example to use the OutputsINI_Var2d function for a user defined function with multi components 3 7 2 The Continuation_Var2d function Continuation Continuation_Var2d Coefficient_name Coefficient Continuation Table 3 30 The Continuation_Var2d function The Continuation
13. At y Yer 7yy 2 2 15 x e x e Phol VT Phol VT With the above initial data ground states for rotating gazes can be obtained for Q lt yy while this is not e g the case with 2 12 when the rotating term is too large In GPELab all these possibilities can be found in the GaussianInitialData2d function where the following extended version is coded y 1 2 bno x xe NG x xe x l 2 16 0 2 TA rol x NB x Xe on with 1 4 1 4 Pho X tg ent Pio X a a iqyy e A 2 17 Two examples of initial data using 2 16 2 17 are presented on Figure 2 1 for the discrete L norm defined by 2 58 and function L2_norm2d Iphi ey Iphi ey 01 008 0 08 0 04 0 02 5 4 3 2 1 a 1 2 3 4 5 o 5 a Q 0 b Q 0 99 Figure 2 1 Two possible gaussian initial data for initializing the iterative scheme with yz yy 1 and Xy 0 0 22 CHAPTER 2 COMPUTATION OF STATIONARY STATES FOR THE GPE In the case of a strong linearity one may also consider the Thomas Fermi TF approximation of the ground state as initial data For the 2d case and a quadratic potential the TF approximate function is such that qm A if TF gt V x eS d 0 otherwise The eigenvalue approximation wpe is given by 1 4872 yy m 2 2 The corresponding func tion is Thomas_Fermi2d An example is given on figure More details about these functions as well as InitialData_Var2d can be found in
14. Opp t 2 y Ap t 2 y a2 Yylyl a yell alu b t x y BIp t x y YCE x y iQ yd 10y Y t y with the parameters a 1 2 k 0 3 Y Yy 1 B 1000 and Q 3 5 We know that the default nonlinearity is the cubic nonlinearity and the default gradients define the rotational operator Therefore we only have to define the quadratic plus quartic potential by using a Matlab script Table 4 18 and then construct the Physics2D structure with the desired coefficients function Potential quadratic_plus_quartic_potential2d gamma_x gamma_y alpha kappa X Y Potential 1 alpha gamma_x X 72 gamma_y Y 72 2 kappa gamma_x X 72 gamma_y Y 72 2 72 Table 4 18 Defining the quadratic plus quartic potential We can now construct the Physics2D structure by using the Physics2D_Var2d function and then add the quadratic plus quartic potential that we have defined the default nonlinear operators and the default gradients to the Physics2D structure thus including them in the physics of the problem see Table 4 19 4 4 GROUND STATE OF 2D GPE EQUATION WITH A FAST ROTATION 73 Delta 0 5 Beta 1000 Omega 3 5 Physics2D Physics2D_Var2d Method Delta Beta Omega Alpha 1 2 Kappa 0 3 Gamma_x 1 Gamma_y 1 Physics2D Potential_Var2d Method Physics2D X Y quadratic_plus_quartic_potential2d Gamma_x Gamma_y Alpha Kappa X Y Physics2D Nonlinearity_Var2d Me
15. Physics1D Nonlinearity_Varid Method Physics1D Table 4 2 Creating the Physics1D structures We then set the initial function to use for the computation We choose a centered gaussian by setting the InitialData_choice to 1 in the InitialData_Varid function as in Table 4 3 InitialData_choice 1 Phi_O InitialData_Varid Method Geometry1D Physics1D InitialData_choice Table 4 3 Building the initial data We finally set the outputs and the printing informations and then launch the simulation We choose to print informations in the command window every 15 iterations and to draw a figure of the square of the modulus of the solution This is done in Table 4 4 Outputs OutputsINI_Varid Method Printing 1 Evo 15 Draw 1 Print Print_Varid Printing Evo Draw Phi Outputs GPELab1d Phi_0 Method Geometry1D Physics1D Outputs Print Table 4 4 Setting the outputs the Print structure and launching the computation 4 1 GROUND STATE OF A 1D GROSS PITAEVSKI EQUATION WITH AN OPTICAL POTENTIAL AND At the end of the simulation we obtain the following informations on the command window of Matlab Iteration 249 on 1000000 Qutputs of component 1 Square at the origin 0 15144001489396 x radius mean square 3 36090790724700 Energy 26 08386211010322 Chemical potential 38 06922564389800 Energy evolution 0 00000000000002 The solution as also been
16. This ODE can be integrated in time exactly and we obtain Vt tn tn 1 and Vp 1 J 2 lt p lt J 2 BY e y clans tone em e y Similarly for Eq 5 6 one gets Vt tn tn 1 and Yq 1 K 2 lt q lt K 2 pe La AA GH ta Thus the first step of the Lie scheme where we solve equation 5 5 and then equation 5 6 on the time interval tn tn 1 and the domain Oj x will be implemented in the following way J 2 1 1 VO tati Tj Yu a gt 1 ain yp etait La p J 2 y ERA ya tni Tj Yk K y A O54 Ent tm GH AE ae alurtLa q K 2 Of course these operations are based on the fftQ and ifft Matlab functions The frequency de pendent quantities 15 2 QyK Mp and A2 2 Oar Aq are obtained through the Delta_grad_Fourier to compute 4p and q and the meshgrid functions Then the exponential matrix is computed by the usual exponential Matlab function To discretize 5 8 we use the standard Simpson s quadrature rule in time Q tn 1 1 1 V s 25 Yn ds 6 V tn j Yk 6V tn41 2 Zj Yk V tn 1 j Yk tn 1 tn tn Vales Y ni ta with tn 1 2 tn tn 1 2 leading to PO tr Ej yk BO tna y yale AO a e This corresponds to only change the phase of the solution by a suitable shift Let us also remark that all we have done above extend directly to the case of a general nonlinearity f x This scheme is globally first order in time and spectr
17. pa A V n ERL g ati Fete In 16 the authors prove that the scheme is energy diminishing if and only if a CFL type condition is satisfied which imposes a restriction on the time step Indeed At must be sufficiently small in this case This is no such restriction for the Backward Euler scheme for the CNGF Indeed it can be proved that BESP is unconditionally energy diminishing this is typically drastically different from what is usually met in a time dependent problem with a real time and not an imaginary time Concerning 2 44 2 45 the fully discrete system is ANG peN 2 46 blo ort where ACN and bEN2 are such that B AN ATP ANG agora 1 Suvi gelee y 2 47 ao FUAN 392 y and pom a HAN IV 3806 Olza y 2 48 The corresponding numerical method is coded in CNSP_CNGF2d for the two dimensional case with a pseudospectral spatial discretization and CNFD_CNGF2d for the finite difference scheme The output physical quantities are again the same as these provided by the BEFD scheme see subsection 2 3 2 page f 2 5 One and three dimensional cases 2 5 1 One dimensional case Al the functions that are found in the two dimensional case have been developed for the one dimensional case In this situation there is clearly no rotation term All functions can be found in the directory Code1D and have the same corresponding names as for the two dimensional case but wi
18. 1 Ne For example we want to set a quadratic potential with an intensity which varies in time We would like to set V t x 3 cos t x for each component To this end we first need to create a square cell array of functions such that the diagonal part is the desired potential and 0 otherwise This is done in table 5 4 where we set the cell array of functions and then modify the Physics2D by using the TimePotential_Var2d function 5 5 GPELAB FUNCTIONS FOR THE DYNAMICS 97 function P Example_timepotential Method Ncomp Method Ncomponents P cell Ncomp for j 1 Ncomp for k 1 Ncomp if j k P j k t x y 3 2 cos t r x 72 y 72 else P j k t x y O end end end Physics2D TimePotential_Var2d Method Physics2D Example_timepotential Method Table 5 4 An example to use the TimePotential_Var2d function Physics2D StochasticPotential_Var2d Method Physics2D StochasticPotential G StochasticProcess Table 5 5 The StochasticPotential_Var2d function 5 5 3 The StochasticPotential_Var2d function In the case where one wants to compute the dynamics of GPE GPELab offers the possibility of handling a stochastic potential That is a potential that is defined using the derivative of a stochastic processes denoted here by w t The StochasticPotential_Var2d function allows to define the stochastic potential operator i e V w t x in the problem by modify
19. 1 1 Phi X Y 1 2 Beta_coupled 1 1 abs Phi 1 72 1 2 Beta_coupled 1 2 abs Phi 2 2 CoupledCubicEnergy 2 2 Phi X Y 1 2 Beta_coupled 2 2 abs Phi 2 72 1 2 Beta_coupled 2 1 abs Phi 1 2 CoupledCubicEnergy 1 2 Phi X Y 0 CoupledCubicEnergy 2 1 Phi X Y 0 Table 4 25 Defining the energy associated with the coupled nonlinearity We can now build the Physics2D structure accordingly to our problem We set the coefficients for the Laplacian the nonlinearity and the rotational operator by using the Physics2D_Var2d function Then we add the default potential the default gradients and the coupled cubic nonlinearity to the physics of the problem thanks to the functions associated to each operator see Table 4 26 76 CHAPTER 4 EXAMPLES OF SIMULATIONS FOR STATIONARY SOLUTIONS Delta 0 5 Beta 200 Beta_coupled 2 1 1 1 Omega 0 8 Physics2D Physics2D_Var2d Method Delta Omega Physics2D Potential_Var2d Method Physics2D Physics2D Nonlinearity_Var2d Method Physics2D Coupled_Cubic2d Method Beta_coupled Coupled_Cubic_energy2d Method Beta_coupled Physics2D Gradientx_Var2d Method Physics2D Physics2D Gradienty_Var2d Method Physics2D Table 4 26 Building the Physics2D structure We then set the initial function as the Thomas Fermi approximation see Table 4 27 InitialData_choice 2 Phi_O InitialData_Var
20. 1 5 Bug reports In case of problems or bug please send an email to the following address gpelab univ lorraine fr We will do our best to provide an answer rapidly 1 6 Copying conditions GPELab is a free software this means that everyone is free to use it and to redistribute it on a free basis GPELab is not in the public domain it is copyrighted and there are restrictions on its distribution but these restrictions are designed to permit everything that a good cooperating citizen would want to do What is not allowed is to try to prevent others from further sharing any version of GPELab that they might get from you 1 6 COPYING CONDITIONS 17 Specifically we want to make sure that you have the right to give away copies of GPELab that you receive source code or else can get it if you want it that you can change GPELab or use pieces of GPELab in new free programs and that you know you can do these things To make sure that everyone has such rights we have to forbid you to deprive anyone else of these rights For example if you distribute copies of GPELab you must give the recipients all the rights that you have You must make sure that they too receive or can get the source code And you must tell them their rights Also for our own protection we must make certain that everyone finds out that there is no warranty for GPELab If GPELab is modified by someone else and passed on we want their recipients to know that what they h
21. 107 108 APPENDIX A COPYRIGHTS amp CREDITS Appendix B License GNU GENERAL PUBLIC LICENSE Version 3 29 June 2007 Copyright 2007 Free Software Foundation Inc http fsf org Everyone is permitted to copy and distribute verbatim copies of this license document but changing it is not allowed Preamble The GNU General Public License is a free copyleft license for software and other kinds of works The licenses for most software and other practical works are designed to take away your freedom to share and change the works By contrast the GNU General Public License is intended to guarantee your freedom to share and change all versions of a program to make sure it remains free software for all its users We the Free Software Foundation use the GNU General Public License for most of our software it applies also to any other work released this way by its authors You can apply it to your programs too When we speak of free software we are referring to freedom not price Our General Public Licenses are designed to make sure that you have the freedom to distribute copies of free software and charge for them if you wish that you receive source code or can get it if you want it that you can change the software or use pieces of it in new free programs and that you know you can do these things To protect your rights we need to prevent others from denying you these rights or asking you to surrender the rights Therefore y
22. GED r Finally the norm lo of a discrete vector of Ne components is defined by Ne vB e C Ne Bllo lleill 2 58 l 1 where the discrete norm for each component is given by 2 58 34 CHAPTER 2 COMPUTATION OF STATIONARY STATES FOR THE GPE For solving 5 28 we again use the preconditioned BiCGStab Concerning the TF precondi tioner since we have coupling between gazes through V and F 6 Anen is a nondiagonal matrix We propose here to extract the diagonal part of Ann to build the preconditioner which means that we only consider potential and nonlinear effects in each single component Concretely we build the following diagonal TF preconditioner P s given by PRE ag Eal 4 vain llas where Vaiag Vi r 1 ve and Faiag Fiu i 1 ne We can also directly inverse the matrix App First we remark that this matrix is composed of diagonal block matrices The following formula can be used to directly compute the inverse of a block matrix A BY _ AT ATBSICATL A 1BS71 C D E SCA se where 4 B C D are diagonal matrices and S D C A7 B is the Schur complement of A assuming that A can be inverted We note that the result of this computation is still a matrix composed of diagonal block matrices By using recursively this formula for the smaller blocks of the matrix i 1 1 1 1 1 1 ahe eB Eee ABS P Cp Dp p CpA Sp where Sp Dp Cp Ap 1Bp and
23. RISK AS TO THE QUALITY AND PERFOR MANCE OF THE PROGRAM IS WITH YOU SHOULD THE PROGRAM PROVE DEFECTIVE YOU ASSUME THE COST OF ALL NECESSARY SERVICING REPAIR OR CORRECTION 16 Limitation of Liability IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRIT ING WILL ANY COPYRIGHT HOLDER OR ANY OTHER PARTY WHO MODIFIES AND OR CONVEYS THE PROGRAM AS PERMITTED ABOVE BE LIABLE TO YOU FOR DAMAGES INCLUDING ANY GENERAL SPECIAL INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPER ATE WITH ANY OTHER PROGRAMS EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES 17 Interpretation of Sections 15 and 16 If the disclaimer of warranty and limitation of liability provided above cannot be given local legal effect according to their terms reviewing courts shall apply local law that most closely approx imates an absolute waiver of all civil liability in connection with the Program unless a warranty or assumption of liability accompanies a copy of the Program in return for a fee END OF TERMS AND CONDITIONS 119 How to Apply These Terms to Your New Programs If you develop a new program and you want it to be of the greatest possible use to the public the best way to achieve this is to make it free softwa
24. Tmax is not known a priori but rather fixed by a stopping criterion to check the convergence of the iterative scheme towards the ground state solution see Eq for example Until now GPELab only includes uniform time stepping in time for a fixed time step At For the numerical purpose the scheme still requires to be space discretized To this end we use a finite difference discretization here Since the domain is R we have to set suitable boundary conditions Here we impose the homogeneous boundary condition d x 0 for x on the boundary of a large enough computational box O Lz Lz x Ly Ly assuming that the physics takes place inside this box Moreover we introduce the indices of the spatial grid points 25 Yn for j k Dyk setting Dix U k EN 1 lt j lt J land1 lt k lt K 1 with J K gt 3 and uniform discretization steps hy and hy in the x and y directions respectively Therefore for 1 lt j lt J 1 he 23 By 2E J 2 3 GRADIENT FLOW FORMULATION AND DISCRETIZATION 25 and forl lt k lt K 1 hy Yk Yk 1 2Ey E The rotating term L is discretized by a two point second order scheme We associate a matrix L to this discrete operator and denote by q OTG k GED I the unknown vector where we assume that the global numbering is made by a local to global reordering procedure based on I j k j J 1 k 1 which corresponds to using the reshape Matlab function when c
25. _ fJ Gradientx x y ifj k for each j k 1 Ne If Gradientx is a cell array of functions in Cn Ne UF Mn N RJ Mwy C then the variable coefficients are G z y Gardientx j k x y for j k 1 Ne The default argument is the part of the rotational operator corresponding to f iQy if j k for the Gradientx_Var2d function where Q is to the rotational speed defined in the Physics2D structure i e Omega R In the case of the Gradienty_Var2d function we have Note that in the 1d case the default argument is In the 3d situation the default operator is the following rotational operator 1 Q3y Q22 if j k Gj Y 2 8 0 if j e k 2 i Q z E O32 if j k Gin 1 Y 2 0 if j E k 3 _ f i Qex Qiy ifj k G p z Y 2 0 if j E k where 2 corresponds to the rotation vector defined in the Physics3D structure i e Omega M 3 R 52 CHAPTER 3 HOW TO USE GPELAB STATIONARY SOLUTIONS e G My n C ones N_c is a variable that multiplies the gradient operator element by element Gj p z y G j k Gradientx j k x y for j k 1 Ne An interesting case of coupling between Gross Pitaevskii equations is the Rashba coupling For example in the case of a system of two Gross Pitaevskii equations with a quadratic potential a coupled cubic nonlinearity and a Rashba coupling we obtain the following equations ily E OOS Balti Boll n G2 Ove 2 i a
26. and 2 21 right 3 1 Ground state computed with GPELab using the parameters from Section 3 2 38 3 2 The FFTNonlinearity_Var2d function 2 0 000000 eee eee 48 ye ees CEN eee es eee ee 65 1 2 Evolution of the energy and the chemical potential during the computation 65 Ty Ser Cree ee 69 Se a ace cies a oars E tox oe 69 a Gabe cg ee hee heehee ene 70 ao hdd 6 Line be ee Be ee sed 71 tee Y pee Oe SER ee ee ee eee ee ee 74 Eis bee oe ee od 76 4 10 107 isovalues of modulus of the ground state o o 78 4 11 107 isovalues of the modulus of the ground state o o 81 5 1 Some outputs computed during the simulation 0 00004 101 5 2 Ground state computed with GPELab using the parameters from Section 5 6 2 102 5 3 Evolution of a dark soliton in a Bose Einstein condensate o 104 5 4 Evolution of a fast rotating Bose Einstein condensate when changing the intensity of the potential s nie fe ake Po aS a e A a Be ee Pete we gies 106 5 5 Evolution of the root mean square in the x and y direction 106 LIST OF FIGURES List of Tables du 36 foe a ee ee ot ae ee ee ee aaa A eo 37 cacas fennel a ath 37 eer 37 Meee eee oboe ot ese see 37 ret oie eet eae pe ee be oot 38 du 3 7 An example of element wise multiplication for matrices 39 3 8 An example of cell ee 39 3 9 Ac
27. arrangement you convey or propa gate by procuring conveyance of a covered work and grant a patent license to some of the parties receiving the covered work authorizing them to use propagate modify or convey a specific copy of the covered work then the patent license you grant is automatically extended to all recipients of the covered work and works based on it A patent license is discriminatory if it does not include within the scope of its coverage prohibits the exercise of or is conditioned on the non exercise of one or more of the rights that are specifically granted under this License You may not convey a covered work if you are a party to an arrangement with a third party that is in the business of distributing software under which you make payment to the third party based on the extent of your activity of conveying the work and under which the third party grants to any of the parties who would receive the covered work from you a discriminatory patent license a in connection with copies of the covered work conveyed by you or copies made from those copies or b primarily for and in connection with specific products or compilations that contain the covered work unless you entered into that arrangement or that patent license was granted prior to 28 March 2007 Nothing in this License shall be construed as excluding or limiting any implied license or other defenses to infringement that may otherwise be available to you unde
28. define the dipolar interaction which can be computed by using a FFT via Qin a ef ale x Pa F FEGE F 921 29 2 Using the previous formula we are able to efficiently compute the dipole dipole interaction and we use the FFTNonlinearity_Var3d function to add it to the Physics3d structure In Table 4 34 we can see how to define this type of nonlinearity in Matlab We can now proceed to build the Physics3D structure with the desired coefficients add the default nonlinear operator the quadratic plus quartic potential and the defined dipolar interaction to the Physics3D structure We can see in Table 1 35 how we add each operator to the Physics3D structure by using the Potential_Var3d Nonlinearity_Var3d and FFTNonlinearity_Var3d func tions We choose the Thomas Fermi approximation as the initial function see Table 4 36 80 CHAPTER 4 EXAMPLES OF SIMULATIONS FOR STATIONARY SOLUTIONS function Dipolar_interaction_nonlinearity Dipolar_interaction3d Phi FFTX FFTY FFTZ Dipolar_direction d Cross_norm sqrt FFTY Dipolar_direction 3 FFTZ Dipolar_direction 2 72 FFTZ Dipolar_direction 1 FFTX Dipolar_direction 3 72 FFTX Dipolar_direction 2 FFTY Dipolar_direction 1 72 Scalar_prod FFTX Dipolar_direction 1 FFTY Dipolar_direction 2 FFTZ Dipolar_direction 3 Angle atan2 Cross_norm Scalar_prod NLFFT fftn abs Phi 2 V d 2 4 3 pix 3 cos Angle 72 1 Dipolar_interaction_n
29. dle gr Bl n l _ n n 1 n n 1 n 1 n n LA ja A y CTT 5 24 n 1 2 Yet pyr se EA with the uniform time step At tn 1 tn and where 2 Pltn 1 2 X Y W tn x and V V tn x The initial conditions are 4 x y x and 9 x Blivo x The extension to a general nonlinearity is direct We now have to discretize the operator A QL To this end we use the pseudospectral approximation of the spatial derivatives based on the Fourier series transforms 2 36 2 37 page Keeping the same notations we have the following discretization prt Ren ARE n a pesan 5 25 where APRS pRen and cR amp are such that Ren 1 1 lanti Lipgn 1 2 1 aten iz GUAN vr he Son Ren _ I 1 1 n 1 n 1 2 ol n 5 26 phen ig GIAN 5v RAN OL 9 cRen 28 1p 1 pr 90CHAPTER 5 COMPUTATION OF THE DYNAMICS OF GPE METHODS FUNCTIONS EXAMPLES Like for the stationary case the evaluation of the differential operators is made through the FFT iFFTs while the diagonal matrices in the physical space are applied directly The linear system in 5 27 is solved at each iteration with the BiCGStab A preconditioner in the spirit of the Laplace Thomas Fermi preconditioners see is used to accelerate the convergence of the iterative Krylov subspace solver The resulting scheme is called Relaxation pseudo SPectral scheme ReSP The scheme is second order
30. function i e Wo x for the simu lations Already defined initial data corresponding to the Thomas Fermi approximation or the centered Gaussian are existing in GPELab Note that the user can also create its own initial wave function without using this function The Method Geometry2D and Physics2D structures are needed arguments for the function Optional arguments are e InitialData_Choice N 1 is a variable that must be either 1 if one uses centered gaussians or 2 for Thomas Fermi approximations as initial data The option 3 allows to use the imaginary time method with the BESP scheme to compute ground states for each component where the operators are restricted to their diagonal parts i e the components are decoupled e X0 YO M n R 0 are variables corresponding to the coordinates of the center of the gaus sian or Thomas Fermi approximation as initial data We note that in the 1d case we only have to define XO and in the 3d case ZO is required e gamma_x gamma_y R 1 are variables corresponding to the parameters of the centered gaus sian We note that in the 1d case we only have to define gamma_x and in the 3d case we have to add gamma_z For example if we want to compute a Thomas Fermi approximation for initial data we proceed as in table InitialData_Choice 2 Phi_O InitialData_Var2d Method Geometry2D Physics2D InitialData_Choice Table 3 26 An example to use the InitialData_Var2d function
31. of the kernel K Another example is the derivation operator y which has symbol 7 A full example is given in Section page We remark that this function is effective only in the case of spectral schemes The FFTNonlinearity_Var2d function allows to define the nonlinear operator i e F W t x x in our equation in the problem To this end it modifies the Physics2D structure The Method and Physics2D structures are required arguments and the optional arguments are the following 3 6 SETTING THE PHYSICAL PROBLEM 49 e FFTNonlinearity If we have a function FFTNonlinearity in F Cr nAMN Na C Mn Na R Mn Ne C Mn Na C the physical nonlocal nonlinearity is such that for each j k 1 Ne FFTNonlinearity W t x x Y Ex ifj k Fj U t x y 0 if j 4 y y E Ey J where and are the discrete Fourier frequencies in the x and y directions respectively If FFTNonlinearity is given by a cell array of functions in Cn No F CN N MNN C Ma RP Maye C My C the nonlocal nonlinearity is Fj V t x x y FFTNonlinearity j k W t x Y En Ey for j k 1 Ne The default argument is Cubic2d which corresponds to Df e G Mn n C ones N_c is a variable that multiplies the nonlocal nonlinearity element by element leading to Fy x W t x 2 y G j k FFTNonlinearity j k t x Y Ex Ey for j k 1 Ne e FFTNonlinearity_e
32. outputs computed during the simulation Computation Ground Ncomponents 1 Type BESP Deltat ie 1 Stop_time Stop_crit 1e 8 Method Method_Var2d Computation Ncomponents Type Deltat Stop_time Stop_crit xmin 10 xmax 10 ymin 10 ymax 10 Nx 279 1 Ny 279 1 Geometry2D Geometry2D_Var2d xmin xmax ymin ymax Nx Ny Table 5 13 Building the Method and Geometry2D structures for the computation of a stationary state following Gross Piteavskii equation 1 1 idt x y t AW 2 y t 5 Ixl y U x y t LE z y t U x y t where we choose 8 10000 Therefore we have to add the potential and nonlinear operators using the Potential_Var2d and Nonlinearity_Var2d The resulting code in available in Table Delta 0 5 Beta 10000 Physics2D Physics2D_Var2d Method Delta Beta Physics2D Potential_Var2d Method Physics2D x y 1 2 x 72 y 72 Physics2D Nonlinearity_Var2d Method Physics2D phi x y abs phi 72 Table 5 14 Setting the Physics2D structure to compute the stationary state We wish to use the Thomas Fermi approximation as an initial wave function for the computation 102CHAPTER 5 COMPUTATION OF THE DYNAMICS OF GPE METHODS FUNCTIONS EXAMPLES Therefore we set the initial data using the InitialData_Var2d as in Table InitialData_choice 2 Phi_O InitialData_Var2d Method Geometry2D Physi
33. p is the level of recursity we are able to compute the inverse of the block matrix AREN Therefore we can build the full TF preconditioner paip aal TF full At ad V BIFE Chapter 3 How to use GPELab stationary solutions GPELab is a flexible Matlab toolbox which is able to compute the stationary solution and dynamics of Gross Pitaevskii equations The physical description of the equations that are defined and the mathematical algorithms are those described in the previous Chapters for the stationary situation The code offers the possibility to compute solutions to multi components BECs and to get physical important quantities GPELab works through a main program that calls solvers corresponding to efficient and robust accurate numerical methods developed in the present document From the point of view of the standard user only this main program has to be modified Furthermore for more complex problems the advanced user can defined his own physical inputs potential nonlinearity number of components and outputs new physical interesting quantities The user can also manipulate the computed quantities and draw figures or create movies in relations to its calculations thanks to Matlab functions and already defined and well adapted GPELab visualization functions Finally GPELab also provides the possibility to include stochastic effects into the computations for example in the potential The aim of this chapte
34. printed out during the simulation At the end we obtain the solution given on Figure 4 1 Iphi x of component 1 Figure 4 1 Modulus of the ground state Moreover we can draw the evolution of the energy and the chemical potential during the com putation by plotting Outputs Energy 1 and Outputs Chemical_potential 1 see Figure 4 2 Chemical potential Figure 4 2 Evolution of the energy and the chemical potential during the computation 66 CHAPTER 4 EXAMPLES OF SIMULATIONS FOR STATIONARY SOLUTIONS 4 2 Ground state of a system of 1d Gross Pitaevskii equations with a quadratic potential and a Josephson junction We would like to produce the numerical results in where a system of two Gross Pitaevskii equations coupled with a Josephson junction is considered The system of Gross Pitaevskii equations is the following Op 4A V x 6 Buly BralWal di Ay 4 1 1002 5A V x Ba2IWa1 Bra W11 v2 Ay where is the detuning constant for the Raman transition pjg are the interactions constants and A is the effective Rabi frequency We can already identify the operators that we will need to define the potential operator with the quadratic potential the detuning constant the effective Rabi frequency and the nonlinear operator with the cubic nonlinearities We already have defined these operators by using the script in tables and for the 2d case We create two
35. similar scripts for our case in tables and function P quadratic_potential_Josephsonid Detuning_constant Rabi_frequency P cell 2 P 1 1 x 1 2 x 72 Detuning_constant P 1 2 x Rabi_frequency P 2 1 x Rabi_frequency P 2 2 0 x 1 2 x 72 end Table 4 5 The potential function used for the Josephson junction function NL Josephson_Nonlinearity Beta_11 Beta_22 Beta_12 NL cel1 2 NL 1 1 Phi x Beta_11 abs Phi 1 72 Beta_12 abs Phi 2 72 NL 2 2 Phi x Beta_22 abs Phi 2 72 Beta_12 abs Phi 1 72 NL 1 2 Phi x 0 NL 2 1 Phi x 0 end Table 4 6 The nonlinearity function used for the Josephson junction Moreover we define in Table the energy 2 50 see page which is associated to the nonlinearity function NLE Josephson_Nonlinearity_Energy Beta_11 Beta_22 Beta_12 NLE cel1 2 NLE 1 1 Phi x 1 2 Beta_11 abs Phit1 72 1 2 Beta_12 abs Phi 2 2 NLE 2 2 Phi x 1 2 Beta_22 abs Phi 2 2 1 2 Beta_12 abs Phi i 2 NLE 1 2 Phi x 0 NLE 2 1 Phi x 0 end Table 4 7 The function computing the energy 2 50 page associated to the nonlinearity used for the Josephson junction We now show how to compute the ground state of this system of coupled Gross Pitaevskii equations by using similar parameters as in 17 First we have to build the Method and Geometry1D 4 2 GROUND STATE OF
36. state solution at the end of the computa tions We can observe the creation of an annulus with uniformly spaced vortices inside Iphicy I of component 1 oors ots A 0o14 or oor 0 008 3 0 006 E 0 004 4 0 002 o a E 2 o 2 o a 10 2 x Figure 4 8 Modulus of the ground state 4 5 Ground state of a system of 2d Gross Pitaevskii equations with quadratic potentials rotational operators and coupled cubic nonlinearites The aim of this section is to consider the computation of the ground state of a 2d system composed of two Gross Pitaevskii equations with quadratic potentials rotational operators and coupled cubic nonlinearities The two first structures that need to be defined are the Method and Geometry2D structures In our case we have to set GPELab to simulate two components Moreover we use the BESP scheme for a time step equal to 107 and a spatial grid of 28 1 points in each direction of the domain 15 15 x 15 15 We set the stopping criterion to 107 see Table 4 23 Computation Ground Ncomponents 2 Type BESP Deltat 1e 2 Stop_time Stop_crit 1e 6 Method Method_Var2d Computation Ncomponents Type Deltat Stop_time Stop_crit xmin 10 xmax 10 ymin 10 ymax 10 Nx 278 1 Ny 278 1 Geometry2D Geometry2D_Var2d xmin xmax ymin ymax Nx Ny Table 4 23 Using the Method_Var2d and Geometry2D_Var2d functions to build the Meth
37. terms permissive or non permissive may be stated in the form of a separately written license or stated as exceptions the above requirements apply either way 8 Termination You may not propagate or modify a covered work except as expressly provided under this License Any attempt otherwise to propagate or modify it is void and will automatically terminate your rights under this License including any patent licenses granted under the third paragraph of section 11 However if you cease all violation of this License then your license from a particular copyright holder is reinstated a provisionally unless and until the copyright holder explicitly and finally terminates your license and b permanently if the copyright holder fails to notify you of the violation by some reasonable means prior to 60 days after the cessation Moreover your license from a particular copyright holder is reinstated permanently if the copy right holder notifies you of the violation by some reasonable means this is the first time you have received notice of violation of this License for any work from that copyright holder and you cure the violation prior to 30 days after your receipt of the notice Termination of your rights under this section does not terminate the licenses of parties who have received copies or rights from you under this License If your rights have been terminated and not permanently reinstated you do not qualify to receive new license
38. the Omega parameter Thus we first have to create the continuation structure for a parameter and then update it by inserting a second parameter This is done in Table where we first define the continuation on the Omega parameter and then we add the G matrix of the nonlinearity Continuation Continuation_Var2d Omega 1 2 3 3 5 4 GNL1 1 0 0 1 GNL2 0 1 1 0 for i 1 5 GNL i GNL1 i 5 GNL2 end Continuation Continuation_Var2d GNonlinearity GNL Continuation Table 3 33 An example to use the Continuation_Var2d function for 2 parameters 3 7 3 The Print_Var2d function Print Print_Var2d Printing Evo Draw Table 3 34 The Print_Var2d function The Print_Var2d function builds the Print structure The aim is to provide to the program the printing informations displayed during the computation The following optional arguments are e Printing N 1 is a variable equals to 1 for printing informations during the computation and 0 otherwise e Evo N 5 is a variable corresponding to the number of iterations between each displayed information including drawing some figures It must be bigger or equal to Evo_outputs from the OutputsINI_Var2d function see Section 3 7 1 page 53 e Draw N 1 is a variable equal to 1 if the modulus and the phase of the wave functions are drawn during the simulation and 0 if not For example if one wants to print informations every 10 itera
39. the system of two components Bose Einstein condensate is modeled by the following system of equations 10101 za V x 8 Bly bral pi Aya 101 V2 za V x B22 W2 Bl Ya y where is the detuning constant for the Raman transition P are the interactions constants and A is the effective Rabi frequency Thus we have to build a potential operator where the diagonal terms are quadratic potentials plus the detuning constant for the first component and the extradiagonal terms are the effective Rabi frequency A To this end we have to create a cell array of functions and then we modify the Physics2D structure to define the potential operator as it is done in table function P Example_potential Detuning_constant Rabi_frequency P cell 2 P 1 1 x y 1 2 x 72 y 72 Detuning_constant P 1 2 x y Rabi_frequency P 2 1 x y Rabi_frequency P 2 2 x y 1 2 x 72 y 72 end Detuning_constant 1 Rabi_frequency 5 Physics2D Potential_Var2d Method Physics2D Example_potential Detuning_ constant Rabi_frequency Table 3 19 An example to use the Potential_Var2d function 3 6 3 The Nonlinearity_Var2d function Physics2D Nonlinearity_Var2d Method Physics2D Nonlinearity G Nonlinearity_energy Table 3 20 The Nonlinearity_Var2d function The Nonlinearity_Var2d function allows to define the nonlinear term i e F W t x x in the problem b
40. tt Pt 2 22 p x tn 1 x tit I x tt llo x 0 0 x x R with o 1 In the above equations we set x t lim yit x t Hence iterations in times correspond to iterations in the projected gradient It is proved in that the CNGF is normalization conserving and energy diminishing if 0 and the potential is positive Another interpretation is that the 24 CHAPTER 2 COMPUTATION OF STATIONARY STATES FOR THE GPE CNGF is a first order time splitting scheme with discontinuous coefficients When t tends towards infinity gives an approximation of the steady state solution which is a critical point of the energy when the assumption on V is fulfilled V gt 0 The initial guess dp is chosen according to the possible choices provided in section Finally let us remark that in our notations we write x t and not t x like for the dynamical case to insist on the fact that t is not a real time but rather a continuation parameter 2 3 1 Time and space discretizations the Backward Euler scheme Different schemes can be considered for computing ground states In 16 the authors show that the Time Splitting sine Spectral TSSP and the Backward Euler Finite Difference schemes BEFD are well adapted when no rotation is included TSSP is supposed to be fast since this is an explicit scheme with FFT based spatial discretization but it requires very small time steps when it is used for ground states co
41. w t x 1 2 x wWw1 t y i2 t for each component where w and wz are brownian motions To this end we first need to create a square cell array of functions such that the diagonal part is the desired potential and 0 otherwise This is done in table where we set the cell array of functions then create a cell vector of two brownian motions and then modify the Physics2D by using the StochasticPotential_Var2d function function P Example_stochasticpotential Method Ncomp Method Ncomponents P cell Ncomp for j 1 Ncomp for k 1 Ncomp if j k P j k W x y 1 2 x 2 W 1 y 2 W 2 else P j k W x y O end end end functionW Example_stochasticprocess Method W cell 1 2 W 1 Brownian_Process2d Method W 2 Brownian_Process2d Method end Physics2D TimePotential_Var2d Method Physics2D Example_stochasticpotential Method Example_stochasticprocess Method Table 5 6 An example to use the StochasticPotential_Var2d function The initial data can be set using the InitialData_Var2d function see section page just like in the stationary case The outputs of a computation are set using the OutputsINI_Var2d function see section 3 7 1 page 53 We remark that the outputs are computed like in the stationary case that is we have to set the variable Evo_outputs which corresponds to the number of iterations between each computation of the outputs Here the number of iteratio
42. where we multiply the ground state Phi_1 by the previous function We finally set the outputs and the printing informations then we launch the simulation We choose to print informations in the command window every 15 iterations and to draw a figure of 5 6 EXAMPLES OF COMPUTATIONS 103 Computation Dynamic Ncomponents 1 Type Relaxation Deltat ie 3 Stop_time 1 5 Stop_crit Method Method_Varid Computation Ncomponents Type Deltat Stop_time Stop_crit Table 5 17 Building the Method structure for a dynamical problem X_0 5 Delta_theta_0 pi 3 s 0 2 Phi_1 1 Phi_1 1 exp 1i Delta_theta_0 2 1 tanh Geometry2D X X_0 s Table 5 18 Phase imprinting the ground state with a dark soliton save 1 Evo_save 100 Outputs OutputsINI_Varid Method save Evo_save Printing 1 Evo 15 Draw 1 Print Print_Varid Printing Evo Draw Figure Figure_Varid hot Phi Outputs GPELab1d Phi_1 Method Geometry1D Physics1D Outputs Print Figure the square of the module of the solution Moreover we choose to draw the figure using the hot colormap which is done using the function Figure_Var2d and to save the solution At the end of the simulation we obtain the saved solution in the outputs which can be used to print the modulus of the solution and show the evolution of the dark soliton in the Bose
43. work from the predecessor in interest if the predecessor has it or can get it with reasonable efforts You may not impose any further restrictions on the exercise of the rights granted or affirmed under this License For example you may not impose a license fee royalty or other charge for exercise of rights granted under this License and you may not initiate litigation including a cross claim or counterclaim in a lawsuit alleging that any patent claim is infringed by making using selling offering for sale or importing the Program or any portion of it 11 Patents A contributor is a copyright holder who authorizes use under this License of the Program or a work on which the Program is based The work thus licensed is called the contributor s contributor version A contributor s essential patent claims are all patent claims owned or controlled by the con tributor whether already acquired or hereafter acquired that would be infringed by some manner permitted by this License of making using or selling its contributor version but do not include claims that would be infringed only as a consequence of further modification of the contributor version For purposes of this definition control includes the right to grant patent sublicenses in a manner consistent with the requirements of this License Each contributor grants you a non exclusive worldwide royalty free patent license under the contributor s essen
44. 0 004 6 0 005 0 002 3 o 10 o o 8 6 4 2 o 2 4 6 8 10 o 8 6 4 2 o 2 4 6 8 10 x x a Ye Yy 1 Xo 1 0 d 1 wo 10 b Ye Ww l a 1 2 k 0 3 10 R 6 4 2 gt 0 2 4 6 8 10 1 Figure 2 3 Examples Thomas Fermi approximation for potentials 2 13 left and 2 19 right Thomas Fermi Approximation Thomas Fermi Approximation x10 10 R 6 4 2 gt 0 10 12 0 08 0 07 10 0 06 8 0 05 gt 6 0 04 0 08 4 d 0 02 4 2 4 0 01 sl 10 o o 8 6 4 2 o 2 4 6 8 10 o 8 6 4 2 o 2 4 6 8 10 x x a Ye Yy 1 a1 a2 25 d da 4 b Ye Yy 1 2 Vo 40 d 0 5 2 Figure 2 4 Examples Thomas Fermi approximation for potentials 2 20 left and 2 21 right 2 3 Gradient flow formulation and discretization One classical solution for computing the solution to 2 9 2 10 is through the projected gradient method which is also called imaginary time method in the Physics community This is the basic method that is coded in GPELab for computing minimal solutions to 2 10 The method consists in i computing one step of a gradient method and then ii project the solution onto the unit sphere S Let us denote by t lt lt tn lt the discrete times and by Atn tn 1 tn the local time step The Continuous Normalized Gradient Flow CNGF is given by p V ox Ep alh se Vo Blol 9 QL tn lt t lt trys gt b x
45. 1 At x Then we set Wt x Ya2 n 1 At x http techmath uibk ac at mecht research SpringSchool manuscript_thalhammer pdf 5 1 ALTERNATE DIRECTION IMPLICIT ADI TIME SPLITTING PSEUDO SPECTRAL ADLTSSP SC Time discretization and the ADI technique Throughout this section we consider the bounded domain O Lz Le x Ly Ly used in the BESP scheme and the time discretization to lt ty lt lt tn lt with At At tn 1 tn Let us fix n N and Y L O that we assume to be numerically compactly supported inside O The first step of the Lie splitting consists in solving between times t t and t try 10401 t x Ab t x _ OD t x 5 4 p tn x y x in O For a non rotating BEC Q 0 this can be done efficiently by means of the FFT since the Laplacian can be inverted in the Fourier space However for Q 0 and since the coefficients of L are not constant this can no longer be done Here we apply the solution proposed in Bao et al based on a ADI scheme The idea is to decouple the global 2d equation as two one dimensional equations where the coefficients are now constant with respect to the directional derivatives Therefore this leads to splitting the equation 5 4 as La first for t tn to b nti 1 iep t x 0000 t x 100 90 t x 65 VO tr x Y x 1 b and then for t tn to t tn 1 f 1 f ibp t x 3 yp t x
46. 1 Beta_22 Beta_12 Table 4 9 Building the Physics1D structure We then set the initial function to use for the computation that we choose as a centered gaussian see Table 4 10 68 CHAPTER 4 EXAMPLES OF SIMULATIONS FOR STATIONARY SOLUTIONS InitialData_choice 1 Phi_O InitialData_Varid Method Geometry1D Physics1D InitialData_choice Table 4 10 Gaussian initial data We finally set the outputs and the printing informations and launch the simulation We display the informations in the command window every 15 iterations We furthermore draw a figure of the square of the modulus of the solution The corresponding GPELab code is detailed in Table Outputs OutputsINI_Varid Method Printing 1 Evo 15 Draw 1 Print Print_Varid Printing Evo Draw Phi Outputs GPELab1d Phi_0 Method Geometry1D Physics1D Outputs Print Table 4 11 Main GPELab code At the end of the simulation we obtain the informations concerning the ground state in Table and the modulus of each component in figure Moreover we can plot the two modulii on the same figure and obtain the same result as in see Figure 4 4 Iteration 164 on 1000000 Qutputs of component 1 Square at the origin 0 03512761887793 x radius mean square 2 67309309073190 Energy 9 97806793424214 Chemical potential 32 84198217411518 Energy evolution 0 00000000000000 Dutputs of component 2 Squ
47. 2 32 The corresponding function is alpha_rms e Es olg The energy Eg o 0y defined at the continuous level by Bald f ZIO VOP Bli ROLA amp 235 All the discretizations use the previous schemes at the interior nodes based on the trapezoidal rule Our function is called energy_beta_psi_Diff2d e go The chemical potential ugo is defined by usa Esa Slocltax 2 34 The GPELab function is chemical_potential and the output is the chemical potential for the converged function This function is coded in such a way that we do not have to precise the discretization scheme e The angular momentum L is defined by Llo f i e vd 0 2 1 dedy 2 35 The GPELab function is Angularl_momentum_Diff2d and the output is the angular momen tum for the function at the converged state e The total cputime until convergence for computing the ground state solution by the scheme 2 3 GRADIENT FLOW FORMULATION AND DISCRETIZATION 27 2 3 3 BESP Spatial discretization pseudospectral scheme based on FFTs Rather than a finite difference scheme pseudospectral approximation of the spatial derivatives can be used to get high order accuracy We consider in GPELab an approach based on Fourier series through FFTs We now impose a periodic boundary condition on the boundary of a large enough computational box O Lz Le x Ly Ly assuming that the physics takes place inside this box and the solution is
48. 2d Method Geometry2D Physics2D InitialData_choice Table 4 27 Constructing the initial data Table provides the informations for the outputs and printing options Outputs OutputsINI_Var2d Method Printing 1 Evo 15 Draw 1 Print Print_Var2d Printing Evo Draw Phi Outputs GPELab2d Phi_0 Method Geometry2D Physics2D Outputs Print Table 4 28 Printing options and launching the computation We report on figure 4 9 the moduli of the stationary state solutions at the end of the computa tions Iphi x y of component x10 ig Iphi x y I of component 2 10 7 0 025 8 3 6 6 6 0 02 4 E 4 2 2 a 0 015 gt 0 gt o 2 2 0 01 4 4 2 e 0 005 1 8 8 10 o 10 o 10 8 6 4 2 o 2 4 6 3 10 10 8 6 4 2 o 2 4 6 8 10 x x a Component 1 of the ground state b Component 2 of the ground state Figure 4 9 Ground state obtained at the end of the simulation 4 6 GROUND STATE OF 3D GPE EQUATION WITH ROTATION TT 4 6 Ground state of a 3d Gross Pitaevskii equation with a quadratic potential a cubic nonlinearity and a rotational operator In this section we consider the problem of computing the ground state of a Gross Pitaevskii equation with quadratic potential cubic nonlinearity and a rotational operator in 3d We first build the Method and Geometry3D structures We use a BESP scheme and a step time of 107 with a spatial grid of 2 1 points in the z y and z dire
49. A SYSTEM OF 1D GROSS PITAEVSKII EQUATIONS WITH A QUADRATIC P structures We use the BESP scheme and a time step equal to 107 with a spatial grid of 21 1 points in 16 16 Moreover the stopping criterion is fixed to 107 see Table 4 8 Computation Ground Ncomponents 2 Type BESP Deltat ie 1 Stop_time Stop_crit le 6 Method Method_Varid Computation Ncomponents Type Deltat Stop_time Stop_crit xmin 16 xmax 16 Nx 2710 1 Geometry1D Geometry1D_Varid xmin xmax Nx Table 4 8 Building the Method and Geometry1D structures We now define the physics of the problem We build the Physics1D structure according to equation 4 1 Moreover one wishes to reproduce the numerical results in with the set of physical parameters S24 del B 500 bu 8B Piz 0 948 B22 0 978 We proceed as in table 4 9 where we use the two operators defined in Tables and and the function defining the energy associated to the nonlinearity in Table 4 7 Delta 0 5 Rabi_frequency 1 Detuning_constant 0 Beta 500 Beta_11 1 Beta_12 0 94 Beta_22 0 97 Physics1D Physics1iD_Varid Method Delta Beta Physics1D Potential_Varid Method Physics1D Josephson_Potential Detuning constant Rabi_frequency Physics1D Nonlinearity_Varid Method Physics1D Josephson_Nonlinearity Beta_11 Beta_22 Beta_12 Josephson_Nonlinearity_Energy Beta_1
50. ATIONS 105 Delta 0 5 Beta 1000 Omega 3 5 Physics2D Physics2D_Var2d Method Delta Beta Omega Alpha 1 2 Kappa 0 7 Gamma_x 1 Gamma_y 1 Physics2D Potential_Var2d Method Physics2D X Y quadratic_plus_quartic_potential2d Gamma_x Gamma_y Alpha Kappa X Y Physics2D Nonlinearity_Var2d Method Physics2D Physics2D Gradientx_Var2d Method Physics2D Physics2D Gradienty_Var2d Method Physics2D Table 5 20 Building and defining the Physics2D structure default Print structure using the Print_Var2d function Then we launch the computation using the GPELab2d function This is done in Table Outputs_iterations 100 Outputs_save 1 Outputs OutputsINI_Var2d Method Outputs_iterations Outputs_save Printing 1 Evo 15 Draw 1 Print Print_Var2d Printing Evo Draw Phi_2 Qutputs GPELab2d Phi_1 Method Geometry2D Physics2D Outputs Print Table 5 21 Launching the simulation At the end of the simulation we obtain the saved solution in the outputs which can be used to show the evolution of the rotating Bose Einstein condensate This evolution can be seen on Figure We can also plot the evolution of the x_rms and y_rms This is done in Table 5 22 The resulting figures can be seen on Figure 5 5 Time 0 0 01 1 plot Time Outputs x_rms 1 xlabel Time ylabel Root Mean Square in the x direction Table 5 22 Pl
51. At 5 34 Similarly for a time step At 2 one gets Ager o9p l At P 0 At 5 35 The Richardson s extrapolation method consists in combining the approximations 5 34 and 5 35 to eliminate the error term C At It is easy to see that At 2 2P 1 2PUA1 2 UAt VAt 2 9p 1 u Co 1 o At t This equation provides an approximation va 2 which is at least of order p 1 In fact depending on the scheme used it is possible to even gain a higher order of accuracy This process can be used iteratively to get high order schemes This results in the following triangular table 94CHAPTER 5 COMPUTATION OF THE DYNAMICS OF GPE METHODS FUNCTIONS EXAMPLES Approximations Extrapolations Order p Order p 1 Order p 2 Order p 3 UAt 2Purr 2 UAt UAt 2 VAt 2 a i e Pak UA OP vAarja Vat 2 UAt 4 VAt 4 o WA 4 api _ 2Puat s Uat 4 Pt varjg vat 4 O waryg Wat a UAt 8 VAt 8 2P 1 WaAt 8 ptI TAt 8 a In GPELab the Richardson s extrapolation method is used for the TSSP and ReSP schemes to get an increased accuracy in time 5 5 GPELab functions for the dynamics As for the computation of ground states GPELab contains functions that are used to define variables or functions in order to set the problem in the first place then to compute interesting quantities A typical script to compute a dynamical problem is very similar to a script
52. BESP to use the Backward Euler SPectral discretization scheme see section or CNSP to use the Crank Nicolson SPectral discretization scheme see section P 4 e Deltat I R 1e 3 is a variable corresponding to the time step of the method The time discretization is always uniform e Stop_time R 1 is a variable corresponding to the final time of computation in the case of a dynamic problem see section 5 5 1 e Stop_crit R 1e 8 is a variable corresponding to the stopping criterion 2 31 e Max_iter N 1e6 is a variable corresponding to the maximum number of iterations for a stationary state computation 42 CHAPTER 3 HOW TO USE GPELAB STATIONARY SOLUTIONS Preconditioner S ThomasFermi is a variable that must be either None for a calculation without preconditioner Laplace for the Laplace preconditioner and ThomasFermi for the Thomas Fermi preconditioner Output N 1 is a variable that must either be 1 if one computes outputs during the compu tations or 0 if not Splitting S Strang is a variable corresponding to the type of splitting in the case of a dynamic computation see section 5 5 1 BESP N 0 is a variable that must be either 1 if one uses the Jacobi method or 0 for the Krylov method for the BESP scheme Solver_FD N 0 is a variable that must be either 1 if one uses the direct Gauss solver from Matlab i e backslash 1 or 0 for the Krylov method Iterative_
53. Bose Einstein Condensates BECs based on the Gross Pitaevskii Equation GPE We do not want here to describe the complex physics behind the BECs and GPE but only to state a few well known facts about GPE and explain how to rewrite the physical GPE into a dimensionless GPE which is the model used in GPELab It is also developed in such a way that the user can define its own equations and compute its proper physical outputs of interest Let us assume that the temperature T is much smaller than the critical temperature TT Then we can describe a BEC under a rotation effect through a macroscopic wave function y which depends on the spatial variable x x y z R and the time t gt 0 This function has a dynamic which is governed by a specific nonlinear Schr dinger equation the so called Gross Pitaevskii Equation given by 00 _ SE aE Sy 2 0 Hg A V x N 1UoP OL a The atomic mass is m h is the Planck constant and H is the hamiltonian operator of the system The notation 55 designates the functional derivatives of the energy E of the system y being the complex conjugate of y The number of atoms in the condensate is N and Q is the angular velocity The potential function V is an external trap which depends on x but may also depend on t according to the physical situation The typical example of potential V is the confining harmonic or quadratic potential trap Y 4 wz 1 2 where wz wy and wz a
54. ELab2d function http www mathworks fr fr help matlab ref daspect html 3 7 LAUNCHING THE SIMULATION 61 Phi Outputs GPELab2d Phi_0 Method Geometry2D Physics2D 0utputs Table 3 41 An example of call to the GPELab2d function 3 7 7 The MakeVideo2d function MakeVideo2d Method Geometry2D 0utputs Function VideoName Figure Table 3 42 The MakeVideo2d function The MakeVideo2d function creates a movie from the saved wave functions during the computa tion It assembles snapshots from each saved component s wave function and creates a video Due to the fact that Matlab takes snapshots of the figures the user must be extremely careful when building the movie by not covering the figure during the construction process This function must be provided with the Method structure the Geometry2D structure and the Outputs structure We have the following optional arguments e Function If a function Function in F My n2 C Mp R Mny wz R is provided then the snapshots are considered by applying Function to each component Wilt x x y gt Function y t x y for j f1 Nc If Function is a cell array of functions in Cin F Mn n C My 1 R Mn 11 R then the snapshot are taken from the following matrix Y t x T y gt Function j q t x T y for j 1 Ne The default value is phi X Y abs phi 2 e VideoName S or C1 y S MyVid
55. Einstein condensate This can be seen on Figure We can also plot the evolution of the energy 5 6 3 Dynamic of a rotating Bose Einstein condensate in 2D We now show how to compute the dynamic of a fast rotating Bose Einstein condensate when changing the intensity of the potential The initial data for this simulation is the ground state computed in section Therefore we assume that Phi_1 is the ground state that is already computed First we build the Method and Geometry2D structures We remark that the geometry must be the same as in section We choose to use the fourth order splitting scheme to compute the dynamic of the solution This is done by directly changing the variable Splitting in the Method structure We fix the time step At 107 and the stopping time T 1 The geometry is set on a computational domain O 10 10 x 10 10 and the number of grid points J K 2 This is coded in Table Now we define the physical problem In this case we want to keep the same physical operators as the one used for the computation of the stationary state but we change the parameters of the 104CHAPTER 5 COMPUTATION OF THE DYNAMICS OF GPE METHODS FUNCTIONS EXAMPLES Iti of com phic of component 1 ponent 1 x10 Ip of component 1 x10 x10 20 20 5 15 5 15 5 10 10 4 4 a 5 5 3 a 3 gt 0 gt 0 2 h 2 Ed 2 10 10 1 1 1 15 15 lo 29 lo 29 lo 5 10 15 20 20 15 10 5 0 5 10 150 20 5 o 5 101580 20 15 10
56. GPELab A Matlab toolbox for computing stationary solutions and dynamics of Gross Pitaevskii Equations GPE Copyright 2013 Xavier Antoine Romain Duboscq Universit de Lorraine Institut Elie Cartan de Lorraine UMR CNRS 7502 F 54506 Vandoeuvre l s Nancy Cedex FRANCE Permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and this permission notice are preserved on all copies Contents 1 Introduction amp informations 13 gee ee ee ee So ee Se ea ae eae ee ee ne ts es ee 13 1 2 The Gross Pitaevskii Equation GPE ee eee ee 14 aeee i Bae bok A amp oe A 14 12 2 The dimensionless GPE ee 15 1 3 Which problems can GPELab solve 2 0 0 0 2 0020000 eee 15 1 4 How to read this manyall se s s 2 2 0 000000 a E E a ee 16 To BUE reports a ak ba Pelee SE Re Gg a id eee te bed 16 1 6 Copying conditions s s t e sos p oea e e aE a E E e e a o e n 16 19 pit e anys 19 beh ta ra a be eS 19 pa a A A aa oe gO 20 2 2 Approximate solutions as initial guess and potentials 20 LA bee eee Rd ee bale OR 23 3 discretizations the Backward 24 2 3 2 Backward Euler Finite Difference BEFD 24 oe 27 2 4 Crank Nicolson schemes 2 0 0 ee 29 2 5 One and three dimensional cases o a a a a 29 2 5 1 One dimensional case ooa aa 29 E as eee Widen le a a 30 a A e eae ee ae ao ae ek
57. NS FOR STATIONARY SOLUTIONS Chapter 5 Computation of the dynamics of GPE methods functions examples After having computed the stationary states of a GPE one can be interested in the dynamics of the stationary solutions for a modified GPE for example by changing the potential More generally the dynamics of an initial data is very important to analyze the dynamical properties of a BEC For this reason we first introduce the main schemes that are included in GPELab The first class of methods is Time Splitting SPectral TSSP schemes We briefly explain how these schemes work in Section and how an increased accuracy can be obtained by using suitable weighted approximations of the exponential operators appearing in the splitting formulae Another scheme that is discussed is the relaxation scheme in time coupled with pseudospectral approximation methods based on FFTs This approach is described in Section In Section we explain the Richardson s extrapolation method that is coded in GPELab This simple technique allows to increase easily the order of the time scheme by linear combinations between computations for fractional uniform time steps We next provide in Section 5 5 the different functions that are met in GPELab to compute the dynamics of a given system of GP equation s Finally in Section 5 6 we end the Chapter by a few computational examples to show how to use the GPELab functions for the dynamics and provide the solutions that
58. ONARY SOLUTIONS e 1076 Concerning the geometry the computational domain is O 10 10 x 10 10 and the number of grid points including the boundary points is set to Ny 28 1 and N 28 1 We can see on Figure 3 1 how it is coded in GPELab Computation Ground Ncomponents 1 Type BESP Deltat ie 2 Stop_time Stop_crit 1e 6 Method Method_Var2d Computation Ncomponents Type Deltat Stop_time Stop_crit xmin 10 xmax 10 ymin 10 ymax 10 Nx 278 1 Ny 278 1 Geometry2D Geometry2D_Var2d xmin xmax ymin ymax Nx Ny Table 3 1 An example of Method and Geometry2D in GPELab for computing a ground state The next step is to define the physical problem In our case we want to compute the ground state of the Gross Pitaevskii equation 1 0 W z y t SAW z y t 5 lal ly Yey 1 BIW 2 y t El y t 0 yO a 20y P x y t with 6 0 5 6 300 and Q 0 7 GPELab is designed in such a way that the user may define and add operators of the following types a potential operator a nonlinear operator and gradient operators The potential operator and the nonlinear operator are functions of the space variables and the wave function for the nonlinear operator that are multiplied by the wave function The gradient operators are defined by functions that are multiplied by the partial derivative of the wave function in the space directio
59. OS 59 3 38 The Figure_Var3d function oaa a a 59 TOPESTE ica an Ed 60 E AA E E ds de 60 3 41 An example of call to the GPELab2d function rts tyra 61 bk a a A 61 3 43 An example to use the MakeVideo2d function o o e 61 E E e 63 Lhe eh di ee Pale ERE e AR e BEM eS we 64 E SE A Ge OG ue Re eS 64 ee 64 e E 66 4 6 The nonlinearity function used for the Josephson junction 66 EEE 66 A eee ees 67 Fux ge ae Pe dee e ee we ee 67 4 10 Gaussian initial datali s sosse s i a ane a i a e a a N e a aaa a a a 68 4 11 Main GPELab code 0 aa a a a a 68 A Bioware atte 68 eee eee 69 PRI Cee eee eee ee 70 ee ga Sue Teen tea Aue us STR tue ee eee ee ey a ees 70 aie oe 71 AW r e de 72 Eee Eos Sao yates 72 ERENS EDER ande 73 e oe 6 eo oh a ee Se E 73 4 21 Defining the Outputs and Print structures and then launching the computation using GPELAbD2d J e oe gs ake ete eee Re SR A ye 73 4 22 Printed outputs for the 2d GPE with cubic nonlinearity rotation and quadratic plus CMabe pe gw ete tae faa te gees teas Gee Aa 73 05 Se eh we Ra a a Meee a ae 74 backs hE SA eB e eh he 75 eee 75 a E A Bele ee Se 76 oi A be ERE e do a be e 76 dee e da 76 o hte peo a S 77 78 soe 78 4 32 Creating the Outputs and Print structure and launching the computation 78 ihe e PR Roe REE Re ek Oe ee GS ee 79 O 80 LIST OF TABLES 11 4 35 Constructing the physics of the problem
60. Omega 0 0 0 9 Physics3D Physics3D_Var3d Method Delta Beta Omega Physics3D Potential_Var3d Method Physics3D Physics3D Gradientx_Var3d Method Physics3D Physics3D Gradienty_Var3d Method Physics3D Physics3D Gradientz_Var3d Method Physics3D Physics3D Nonlinearity_Var3d Method Physics2D Table 4 30 Setting the coefficients and adding the default operators to the Physics3D structure We then set the initial function by using the InitialData_Var3d We choose the Thomas Fermi approximation see Table 4 31 InitialData_choice 2 Phi_O InitialData_Var3d Method Geometry3D Physics3D InitialData_choice Table 4 31 Building the initial data as the Thomas Fermi approximation We finally set the outputs and the printing informations then we launch the simulation following Table Outputs OutputsINI_Var3d Method Printing 1 Evo 15 Draw 1 Print Print_Var3d Printing Evo Draw Phi Outputs GPELab3d Phi_0 Method Geometry3D Physics3D Outputs Print Table 4 32 Creating the Outputs and Print structure and launching the computation At the end of the computation we obtain on Figure the isovalues of the modulus of the stationary state solution Ipbitxy 2 of component 1 Figure 4 10 107 isovalues of modulus of the ground state 4 7 GROUND STATE OF 3D GPE EQUATION WITH DIPOLE DIPOLE INTERACTION 79 4 7 Ground state of a 3d Gros
61. Section 3 6 6 page 53 Thomas Fermi Approximation x107 16 14 12 10 8 6 4 R 2 o 8 6 4 2 o 2 4 6 8 10 id x 10 f a 8 1000 Figure 2 2 Thomas Fermi initial data for initializing the iterative scheme for potential 2 11 with Ya Vy 1 As already said many kinds of potentials may be used This is for example the case of a harmonic trap 2 11 quadratic_potential2d function with a possible added exponential term like in 2 13 quadratic_plus_exp_potential2d Examples of these two potentials are given below on Figure 2 4 Other possibilities include e Quadratic plus quartic potential quadratic_plus_quartic_potential2d function 1 K V x 1 0 5 Yer y qa ny Y 2 19 e Quadratic plus sin optical potential quadratic_plus_sin_potential2d function 1 99 22 ar Try an ny Vix see Wy gt sin a ml 2 20 e Double well trapping potential double_well_trapping_potential2d function ER V x 5 va gy oe 2 21 2 Examples of these potentials are given on Figures 2 4 Clearly any new initial data or potential can be added by just following the way the functions are written More details about the potential functions in GPELab are given in Subsection page 2 3 GRADIENT FLOW FORMULATION AND DISCRETIZATION 23 Thomas Fermi Approximation Thomas Fermi Approximation 10 0 016 R 0 025 0 014 8 4 0 012 9 02 0 01 a 0 015 0 008 2 0 006 0 01 4
62. Vj x tbj s ds Y gt Vix w t w tn V w t w tn x tn j l vin j 1 This leads to the following integration for equation 5 32 Wo t x pi n DAt x e ADAL Pt wlt wln 5 4 RICHARDSON S EXTRAPOLATION 93 5 3 2 The relaxation scheme In section 5 2 page 89 we have introduced the relaxation scheme The main problem in the case of our stochastic equation lies in the time discretization of the noise Considering a time discretization to lt ti lt lt tn lt with At At tn 1 tn and an initial wave function w x we propose the following time discretization Nstoch f H Vals xls x ds X gt Vi E A tn tn N stoch tn f tn 3 y V x ll 1 al yl x j l wj tn 1 W3 tn Therefore for our stochastic Gross Pitaevskii equation we obtain the following relaxation scheme grti 2 J gr 1 2 bly 2 py agp l rn y yn n 1 2 yntt yr 5 33 where p 1 2 Pltn 1 2 X Y P tn x and V V cee x The initial conditions are P x Yo x and 1 x Blvo x 1 5 4 Richardson s extrapolation The Richardson s extrapolation is a simple method that linearly combines two solutions given by a numerical scheme on two different time steps to obtain a new solution with an increased accuracy Let us assume that we can numerically compute an approximation ua of an exact quantity u by using a scheme of order p N un u Cp At Cpp AHP o
63. _Var2d function allows to define what continuation means in GPELab when you want to compute a stationary state By continuation we mean that the computation is initialized by an initial state and then doing loops on a chosen parameter p E Cy com C or Ci neont M N Na C where n is the total number of iterations At each iteration k the new initial state is taken to be equal to the stationary state computed at the previous iteration k 1 The continuation continues until the end of the loop and gives a stationary state solution at the final iteration n This way of proceeding is expected to provide at the end of the calculations a ground state solution through small increments of the parameter p at the iteration k In addition the Continuation_Var2d function allows to define a multi parameter continuation path More precisely let us assume that the Continuation_Var2d function is called successively nPerem times with nP different continuation parameters pj Cy neont C or Cincom M N N C j E nParam where n must be the same for all pj Then the multi parameter continuation path is defined by cont path RP 3 7 LAUNCHING THE SIMULATION 57 with path pi k pnparm k and 1 lt k lt net Concerning the Matlab function itself the following arguments are needed e Coefficient_name S identifies the parameter on which the continuation applies This pa rameter can be one of the followi
64. able portion of the object code whose source code is excluded from the Corresponding Source as a System Library need not be included in conveying the object code work A User Product is either 1 a consumer product which means any tangible personal property which is normally used for personal family or household purposes or 2 anything designed or sold for incorporation into a dwelling In determining whether a product is a consumer product doubtful cases shall be resolved in favor of coverage For a particular product received by a particular user normally used refers to a typical or common use of that class of product regardless of the status of the particular user or of the way in which the particular user actually uses or expects or is 114 APPENDIX B LICENSE expected to use the product A product is a consumer product regardless of whether the product has substantial commercial industrial or non consumer uses unless such uses represent the only significant mode of use of the product Installation Information for a User Product means any methods procedures authorization keys or other information required to install and execute modified versions of a covered work in that User Product from a modified version of its Corresponding Source The information must suffice to ensure that the continued functioning of the modified object code is in no case prevented or interfered with solely because modification has been ma
65. al as Nstoch vaw D gt Veo j 1 and ww t gt 0 is the time derivative of a multidimensional noise w t w1 t Wi t Here we will focus ourselves on the two previous schemes that we have introduced to solve such equations the time splitting scheme and the relaxation scheme 5 3 1 The time splitting scheme In the time splitting scheme see section page 83 we essentially split the equation in two equations that we solve separately Here we follow the same procedure In our case given a time discretization to lt ti lt lt tn lt with At At tn 1 tn and an initial wave function x we obtain the following scheme 1 let At gt 0 n N 0 and nAt lt t lt n 1 At we solve ide t x 5 Ad t x ALA t x nAt lt t lt n 1 At 5 31 wi nAt x Y x 2 and next 04 2 t x V w t x yo t x T Blivale x palt x nAt lt t lt n 1 At 5 32 wo nAt x 11 n 1 At x Equation 5 31 is solved using the ADI technique and fast Fourier transforms see section 5 1 1 page 85 Concerning equation 5 32 we have seen in section that we can exactly integrate the nonlinearity and the time dependent potential in time Therefore we obtain for all n t lt t lt n 1 At a t x da n 1 At x e BP DAL tn fen Vl 30 9 The time integration of the stochastic potential gives in fact t Nstoch pt Nstoch f V w s x ds Y gt
66. al in space 5 1 2 The Strang ADI TSSP scheme for the rotating GPE Let us briefly explain the Strang splitting scheme the implementation being quite similar to the Lie scheme for a rotating BEC The approximation of the solution on a time step At will either be e4 4t 2 BAt A At 2 or eB At 2 cAAteB At 2 setting A SA iQL and B iV t x iB w t x Using the ADI method and the first exponential splitting the TSSP Strang scheme is the following 1 solve from t tn to t tn41 2 l idp t x 00000 t x yO t x 69 b tn X Yn x 2 Solve from t tn to t tn41 2 1 iH t x FOyyh t x ad t x 5 10 ya tn x e po tn41 25 x 5 1 ALTERNATE DIRECTION IMPLICIT ADI TIME SPLITTING PSEUDO SPECTRAL ADLTSSP SC 3 Solve from t tn to t tn4i id pp t x VIO Bl Y t 20 0 x a ye En x ya tn 1 2 x by the change of phase 4 Solve from t ty to t th41 2 1 iO t x Oy UD t x Qrp t x 5 12 pO tr x y tny1 x 5 Solve from t ty to t tn41 2 1 pl tn x y En 1 2 x This last step finally gives yp x 40 tn41 2 Xx Each partial differential equation with respect to x or y is solve through FFTs and iFFTs The Lie ADI TSSP scheme is first order in time and spectral in space while the Strang ADI TSSP scheme is second order in time and spectrally accurate in space Their computational costs are O M log M The extension to the three d
67. andards body or in the case of interfaces specified for a particular programming language one that is widely used among developers working in that language The System Libraries of an executable work include anything other than the work as a whole that a is included in the normal form of packaging a Major Component but which is not part of that Major Component and b serves only to enable use of the work with that Major Component or to implement a Standard Interface for which an implementation is available to the public in source code form A Major Component in this context means a major essential component kernel window system and so on of the specific operating system if any on which the executable work runs or a compiler used to produce the work or an object code interpreter used to run it The Corresponding Source for a work in object code form means all the source code needed to generate install and for an executable work run the object code and to modify the work including scripts to control those activities However it does not include the work s System Libraries or general purpose tools or generally available free programs which are used unmodified in performing those activities but which are not part of the work For example Corresponding Source includes interface definition files associated with source files for the work and the source code for shared libraries and dynamically linked subprograms that t
68. arate and independent works which are not by their nature extensions of the covered work and which are not combined with it such as to form a larger program in or on a volume of a storage or distribution medium is called an aggregate if the compilation and its resulting copyright are not used to limit the access or legal rights of the compilation s users beyond what the individual works permit Inclusion of a covered work in an aggregate does not cause this License to apply to the other parts of the aggregate 6 Conveying Non Source Forms You may convey a covered work in object code form under the terms of sections 4 and 5 provided that you also convey the machine readable Corresponding Source under the terms of this License in one of these ways a Convey the object code in or embodied in a physical product including a physical distribution medium accompanied by the Corresponding Source fixed on a durable physical medium custom arily used for software interchange b Convey the object code in or embodied in a physical product including a physical distribution medium accompanied by a written offer valid for at least three years and valid for as long as you offer spare parts or customer support for that product model to give anyone who possesses the object code either 1 a copy of the Corresponding Source for all the software in the product that is covered by this License on a durable physical medium customarily used for
69. are at the origin 0 04853615539806 x radius mean square 2 99189318517412 Energy 13 16187704973455 Chemical potential 38 16552799804442 Energy evolution 0 00000000000000 CPU time 8 28 gt Table 4 12 Printed outputs for the 1d Josephson problem 4 3 Ground state of a 2d Gross Pitaevskii equation with an optical potential and a cubic nonlinearity We now consider the computation of the ground state of a Gross Pitaevskii equation with an optical potential and a cubic nonlinearity in the two dimensional case We use the parameters that are chosen in 21 First we build the Method and Geometry2D structures We use the BESP scheme The time step is At 107 and the spatial grid involves 28 1 points both in the x and y directions The computational domain is 16 16 x 16 16 Moreover the stopping criterion is 1076 We refer to Table 4 13 for the details 4 3 GROUND STATE OF 2D GPE EQUATION WITH AN OPTICAL POTENTIAL 69 phi x of component 2 r T r phi x of component 1 r T T T T T T T T 0 045 0 035 4 0 04 0 03 4 o L 1 1 1 L Y o 1 Y L L L Y 20 15 10 5 o 5 10 15 20 20 15 10 5 o 5 10 15 20 x x a Component 1 of the ground state b Component 2 of the ground state Figure 4 3 Ground state obtained at the end of the simulation phi x of all components Phi CyT Phi x y Figure 4 4 Both components of the ground state Com
70. are expected 5 1 Alternate Direction Implicit ADI Time Splitting pseudo SPec tral ADI TSSP schemes for the rotating GPE Let us introduce A and B two self adjoint operators such that D A c L D B c L and A B is a self adjoint operator on D A ND B We designate by D A and D B the domains of the operators A and B respectively We consider the following time dependent Partial Differential Equation PDE p t x Aw t x By t x v0 x wo x and we denote by y t x e 4 to x its solution t gt 0 and x RY The Time Splitting TS schemes consist in approximating the solution of this PDE by the exponential operators of A and B which means that we approximately split A and B into the exponential representation In fact a general approximation of the solution by a time splitting scheme can be written in the following form W tn At x prtt x em AAt on BAt a2 AAt baBA pap AAt by BAt yn x where ak bx 1 lt k lt p C R are some weights that must be computed in such a way that the approxi mation of e 4 4 is of a given order for a time step At which is supposed to be relatively small 83 84CHAPTER 5 COMPUTATION OF THE DYNAMICS OF GPE METHODS FUNCTIONS EXAMPLES Two of the most commonly used time splitting schemes are the Lie a b 1 cf subsection 5 1 1 and the Strang a a2 1 2 and b 1 b2 0 cf subsection 5 1 2 schemes which are first and second order time schemes res
71. at material or requiring that modified versions of such material be marked in reasonable ways as different from the original version or d Limiting the use for publicity purposes of names of licensors or authors of the material or e Declining to grant rights under trademark law for use of some trade names trademarks or service marks or f Requiring indemnification of licensors and authors of that material by anyone who conveys the material or modified versions of it with contractual assumptions of liability to the recipient for any liability that these contractual assumptions directly impose on those licensors and authors All other non permissive additional terms are considered further restrictions within the meaning of section 10 If the Program as you received it or any part of it contains a notice stating that it is governed by this License along with a term that is a further restriction you may remove that term If a license document contains a further restriction but permits relicensing or conveying under this License you may add to a covered work material governed by the terms of that license document provided that the further restriction does not survive such relicensing or conveying If you add terms to a covered work in accord with this section you must place in the relevant source files a statement of the additional terms that apply to those files or a notice indicating where to find the applicable terms Additional
72. ate numerical methods for computing ground and first excited states in Bose Einstein condensates J Comp Phys 219 pp 836 854 2006 An efficient and spectrally accurate numerical method for computing dynamics of rotating Bose Einstein condensates J Comp Phys 217 pp 612 626 2006 22 L Wen H Xiong and B Wu Hidden vortices in a Bose Einstein condensate in a rotating double well potential Phys Rev A 82 053627 2010 23 X Antoine and R Duboscq In preparation 2012 24 Denschlag J and Simsarian JE and Feder DL and Clark Charles W and Collins LA and Cubizolles J and Deng L and Hagley EW and Helmerson K and Reinhardt WP and others Generating solitons by phase engineering of a Bose Einstein condensate Science 287 pp 97 101 2000
73. ave is not what we distributed so that any problems introduced by others will not reflect on our reputation The precise conditions of the license for GPELab are found in the General Public License that accompanies the source code see Appendix B License page 109 Further information about this license is available from the GNU Project webpage Detailed copyright information can be found in Appendix A Copyright amp credits page 107 If you want to integrate parts of GPELab into a closed source software or want to sell a modified closed source version of GPELab you will need to obtain a different license Please contact us directly for more information GNU License 18 CHAPTER 1 INTRODUCTION amp INFORMATIONS Chapter 2 Computation of stationary states for the GPE 2 1 The Gross Pitaevskii Equation and its properties The aim of this section is to present the physical model the mathematical notations and numerical methods that are developed in GPELab for computing stationary states solution like ground states and excited states Furthermore we explain which Matlab functions correspond to each method and other physical quantities already existing in GPELab Let us note that this part is mathematically a little bit technical If you directly want to make computations Chapter 3 provides the full details on how to use the code Many examples are provided to help the beginner to directly launch some computations in Cha
74. ax x_min y_max y_min We tk 2 gt 2 where x_max x_min y_max and y_min have been defined by the Geometry2D_Var2d function see subsection page 43 x_rms y_rms C1 Na M1 Nou R is a cell array of vectors containing the values of the root mean square of each wave function Ye with respect to the x and y directions They are computed by x_rms k f a hulin 9 Pedy and os ym 000 otro Pray Energy Ci v M1 n R is a cell array of vectors such that Energy k Ep W tx see Equation 3 2 page 44 Chemical_potential Ci v M1 N R is a cell array of vectors such that Chemical_potential k uy W t see Equation 3 3 page 44 User_defined_local Cj you M1 Nou R are the user defined functions userdef_outputs User_defined_global Cy neou M1 Nous R are the user defined functions globaluserdef_outputs For example if one wants to compute the L norm of the gradient of each component of a Bose Einstein condensate on the computational domain Grad_norn Vy t x dedy O one has to first define a function that computes the L norm of the gradient by using a FFT and then create the Outputs structure by using the OutputsINI_Var2d function with the function as argument This is done in table 56 CHAPTER 3 HOW TO USE GPELAB STATIONARY SOLUTIONS function Grad_norm Example_outputs Geometry2D phi x y fftx ffty Grad_x ifft2 li fftx fft2 phi Grad_y
75. bit them from making any copies of your copyrighted material outside their relationship with you Conveying under any other circumstances is permitted solely under the conditions stated below Sublicensing is not allowed section 10 makes it unnecessary 3 Protecting Users Legal Rights From Anti Circumvention Law No covered work shall be deemed part of an effective technological measure under any applicable law fulfilling obligations under article 11 of the WIPO copyright treaty adopted on 20 December 1996 or similar laws prohibiting or restricting circumvention of such measures When you convey a covered work you waive any legal power to forbid circumvention of techno logical measures to the extent such circumvention is effected by exercising rights under this License with respect to the covered work and you disclaim any intention to limit operation or modification of the work as a means of enforcing against the work s users your or third parties legal rights to forbid circumvention of technological measures 4 Conveying Verbatim Copies You may convey verbatim copies of the Program s source code as you receive it in any medium provided that you conspicuously and appropriately publish on each copy an appropriate copyright notice keep intact all notices stating that this License and any non permissive terms added in accord with section 7 apply to the code keep intact all notices of the absence of any warranty and give all
76. cessing to an element of a cell 0 0 oo 0 02 e 39 E e Shoe E a E A dd bee 39 3 11 An example of structure o ee 40 3 12 The Method_Var2d function 2 a Al 3 13 An example of initialization and use of the Method_Var2d function 42 3 14 The Geometry2D_Var2d function 2 2 2 e 43 3 15 An example of how to use the Geometry2D_Var2d function 43 3 16 The Physics2D_Var2d function oaa e 44 patrones iii 45 3 18 The Potential_Var2d function o o a e 45 MEETER 46 3 20 The Nonlinearity_Var2d function oaoa e e 46 3 21 An example to use the Nonlinearity_Var2d function 48 3 22 An example of application of the FFTNonlinearity_Var2d function with the dipolar A AAA 50 3 23 The Gradientx_Var2d functi0M e o 51 Lo 52 La di A A A be oe 53 3 26 An example to use the InitialData_Var2d function o 53 a a ee y eos oe ee ee 53 AAA 5 ae eae aaa chess 56 Je O O 3 30 The Continuation_Var2d function 02 000 ee eee 56 3 31 An example for using the Continuation_Var2d function when the parameter is a a A a EO ae re ee E eS ee ee 57 AR E NE e wa eee he PUES OER 58 3 33 An example to use the Continuation_Var2d function for 2 parameters 58 9 10 LIST OF TABLES wee eee abe bed ee be eee a ARS OR 58 aah eee ee een 59 3 36 The Figure_Var2d function ooa ea 59
77. confined in the box Moreover we introduce the indices of the spatial grid points j Yk for j k Pjg setting Pak j k EN 0 lt j lt J 1and0 lt k lt K 1 with J K gt 2 and uniform discretization steps hy and hy in the x and y directions respectively Therefore for 0 lt j lt J 1 hy j yzi 2Lz J and fr 0 lt k lt K 1 hy yk Yk 1 2Ly K The partial Fourier pseudospectral discretizations in the x and y directions are respectively given by J 2 1 Q T Yk t gt dp Yk tje Mpeg rt J 2 Eo _ 2 36 aj Yk t ape del Tj t yer urta a K 2 where dbp and bq are respectively the Fourier coefficients in the x and y directions pa ea Pp Yk t S D t Yk t e ee 0 E aN 2 37 A Cainer k 0 with uy i and Ay T For the backward Euler scheme this implies that we have the following spatial approximation ABEng pBEn E 2 38 ms where Pr GW EP K is the discrete unknown array in M m C and the right hand side is bP o At with xjk G k Pz Mm C Here Mm C designates the set of 2D respectively 1D and 3D arrays with complex coefficients with Mp JK respectively Mp J and Mp JKL in 2D respectively 1D and 3D For conciseness let us remark that we do not make any distinction between an array in M m C and the corresponding reshaped vector in C P The operator APE is given by the map which for any vector y
78. cs2D InitialData_choice Table 5 15 Initialization by the Thomas Fermi approximation We finally launch the simulation We set the defaults Outputs using the OutputsINI_Var2d function and the default Print structure using the Print_Var2d function Then we launch the computation using the GPELab2d function and store the ground state under the variable Phi_1 This is done in Table 5 16 Outputs OutputsINI_Var2d Method Printing 1 Evo 15 Draw 1 Print Print_Var2d Printing Evo Draw Phi_1 Qutputs GPELab2d Phi_0 Method Geometry2D Physics2D Outputs Print Table 5 16 Launching the computation of the ground state At the end of the computation we obtain the ground state whose modulus is depicted in Figure Iphi x y I of component 1 x10 5 4 3 o 2 1 20 15 10 5 o 5 10 15 20 x Figure 5 2 Ground state computed with GPELab using the parameters from Section 5 6 2 It is now possible to phase imprint a dark soliton in the Bose Einstein condensate and simulate lts time evolution First we have to rebuild the Method structure for this dynamical problem We want to use a Relaxation scheme and a step time of 107 Moreover we want to set the stopping time to 1 5 We coded this in GPELab as in Table We now have to phase imprint the dark soliton 24 This is done by multiplying the initial data by O Oe where A o 7 3 zo 5 and s 0 2 This is done in Table 5 18
79. ctions on the domain 10 10 x 10 10 x 10 10 Moreover we set the stopping criterion to 107 We can see in Table how to proceed to set these parameters in the Method and Geometry3D structures by using the Method_Var3d and Geometry3D_Var3d functions Computation Ground Ncomponents 1 Type BESP Deltat ie 2 Stop_time Stop_crit le 4 Method Method_Var3d Computation Ncomponents Type Deltat Stop_time Stop_crit xmin 10 xmax 10 ymin 10 ymax 10 zmin 10 zmax 10 Nx 276 1 Ny 276 1 Nz 276 1 Geometry3D Geometry3D_Var3d xmin xmax ymin ymax zmin zmax Nx Ny Nz Table 4 29 Building the Method and Geometry3D structures The ground state that we want to compute is associated to the following Gross Pitaevskii equa tion 1 1 itt eyz AW t 2 y 2 gt Aelel vylyl 7 121 E t 2 y 2 6 V t T Y z T t T Y z EE Q x x V V t T Y z with Ys Yy Y 1 8 500 and Q 0 0 0 9 We keep in mind that the default nonlinearity is the cubic nonlinearity the default potential is the quadratic potential and the default gradients operators are the rotational operators This implies that we only have to build the Physics3D structure with the desired coefficients and then add each operator see Table for more details 78 CHAPTER 4 EXAMPLES OF SIMULATIONS FOR STATIONARY SOLUTIONS Delta 0 5 Beta 500
80. d of splitting strategy For the Lie scheme this leads to solving the equation in two steps Let At gt 0 and Y the approximation of the solution at time tn Then the splitting yields 88CHAPTER 5 COMPUTATION OF THE DYNAMICS OF GPE METHODS FUNCTIONS EXAMPLES 1 solve from tn to tn 1 d 00D t x AO x X GI x rj On VO x a 5 15 UO tn x U x 2 Solve from tn to tn 1 1d Y t x V x V t x BFW DIA ea VO tp x VO t 41 x l We furthermore have W x Yo x We remark now that 5 16 leads to Vt gt tn VW t x tn x Indeed we have Ne Y G U tx 2d R HP x UO e x m 1 2 ES WP t x Vinol Fino Y t x VP E x Using the fact that Vino x Vom x and Eno LO t x x Fom t x x we obtain that Ne gt atp t x m 1 2 y A Ne gt o gt m gt 1 2 Y mts Frm Y 39 39 VO 21 N 2m gt 1 4 Y cimas N gt o gt m gt 1 5 17 As a conclusion we have that VO t x V tn x for t gt tn Let n N and W x L O compactly supported in O For the sake of simplicity we explicit the method in the two dimensional case e g d 2 A first step is to solve from t tp to t tn4i 2 0 00 t x SAW x Y GAO VO x 5 18 UO tn x U x As in the one component case an ADI method must be used to decouple the effects of GJ Therefore we split Eq 5 15 in two equations and solv
81. de If you convey an object code work under this section in or with or specifically for use in a User Product and the conveying occurs as part of a transaction in which the right of possession and use of the User Product is transferred to the recipient in perpetuity or for a fixed term regardless of how the transaction is characterized the Corresponding Source conveyed under this section must be accompanied by the Installation Information But this requirement does not apply if neither you nor any third party retains the ability to install modified object code on the User Product for example the work has been installed in ROM The requirement to provide Installation Information does not include a requirement to continue to provide support service warranty or updates for a work that has been modified or installed by the recipient or for the User Product in which it has been modified or installed Access to a network may be denied when the modification itself materially and adversely affects the operation of the network or violates the rules and protocols for communication across the network Corresponding Source conveyed and Installation Information provided in accord with this section must be in a format that is publicly documented and with an implementation available to the public in source code form and must require no special password or key for unpacking reading or copying 7 Additional Terms Additional permissions a
82. e indexation must be done using curly brackets However due to the fact that cells gather elements such as strings basic operations are of course not implemented in Matlab A 1 string 3 1 1 2 1 A 1 2 ans string gt Table 3 9 Accessing to an element of a cell e Function Functions are scripts written in a simple way They are defined as variables making them more direct to create to use and to manipulate than a script written in a m function file To define a function one has first to give the input arguments then write the script f x y a b x strcmp a b y 2 3 world phone ans 3 gt Table 3 10 An example of simple function e Structure type Structures are variables that simply gather variables Contrary to cells there is no indexation and the access is done by directly naming the variables 3 4 Notations and preliminary remarks First let us introduce some general notations to understand the types of the input and output arguments in the GPELab functions Let us define 40 CHAPTER 3 HOW TO USE GPELAB STATIONARY SOLUTIONS S temperature 32 S computer office S temperature 32 computer office S computer ans office Table 3 11 An example of structure e Nz Ny Nz these parameters are equal to the number of degrees of freedom dof of the numerical method that is considered in the z y and z dir
83. e by one components They are evaluated through nt Matlab functions that must be of the form W t x 2 Y Ex Ey gt globaluserdef_outputs j U t x Y Ex Ey for j 1 n vt where and are the discrete Fourier frequencies in the x and y directions We remark that globaluserdef_outputs must have Ye t Xx Y Ex Ey as arguments only in the case where a spectral scheme is used Otherwise the arguments are we t x 2 y By default there is no predefined output quantity in GPELab which means that the user must define its own functions e globaluserdef_outputs_names C cou S User defined function is a variable that has the same role as userdef_outputs_names but for globaluserdef_outputs Let us assume that we launch a simulation that ends after Niter iterations Therefore the outputs are computed N Nout Int no 1 Evo_outputs times at ty k Evo_outputs At 1 lt k lt Nou and ty 41 NiterAt In the above equation Int r designates the integer part of a real valued number r The resulting Outputs structure has the following variables 3 7 LAUNCHING THE SIMULATION 55 Solution C1 Nou C1 N1 MN N C contain the computed solutions for times t if save 1 phi_abs_0 C1 N M1 N a R is a cell array of vectors that contains the values of the square modulus of each wave function Y at the center of the domain for times tk phi_abs_0 4 k i x_m
84. e e Go ae ee 30 2 6 1 The multi components GPE 0 0 0 00 0000000 ee ee ee 30 2 6 2 Stationary states and the CNGF 0 o e 20040 31 A A II 32 3 How to use GPELab stationary solutions 35 3 1 How to get and install GPELab 0 20 02 0 2 020000 0000 2 eee 35 a baw CARB o be 4b AAA 35 3 3 Variables types and various notions required for Matlab 38 bags 2G 6 ak ee a See Shee eid e ae 39 3 5 Setting the numerical scheme and the geometry 2 2004 41 3 5 1 The Method_Var2d function 2 2 2 2 20 0000 02 ee eee 41 aa a ale a a eps 43 ad e a ds ap A a 43 pl a a arb hohe ace Eke 44 Lal eb ak a wa ee Re Ee oO 45 3 6 3 The Nonlinearity_Var2d function 20 02000 46 3 4 CONTENTS 3 6 4 The FFTNonlinearity_Var2d function 00 48 3 6 5 The gradient functions o aaa ee 51 a le A 53 A Hae A A Se ee 53 3 7 1 The OutputsINI_Var2d function e e 53 DEGREE 56 a os Sk oe rr A a 58 3 7 4 The Figure_Var2d function 2 0 59 3 7 5 The Figure_Var3d function 2 e 59 3 1 6 The GPELab2d function 60 3 7 7 The MakeVideo2d function 0 0 00000 ee 61 4 Examples of simulations for stationary solutions 63 4 1 Ground state of a 1d Gross Pitaevskii equation with an optical potential and a cubic Sit iwy Bead a e ee ee ee Yaa eee Ee ee ee ae Ss 63 Wa nonlinearity and a rotational operator
85. e from t tn to t tn 1 idt UD t x O VOD t x G y 0 t x 5 19 TOD ta xX x and next from t tn to t tn41 i O LOD t x 5 py T t x G7 a U0 t x 5 20 WO ta x WO t41 X 5 2 RELAXATION PSEUDO SPECTRAL SCHEME RESP 89 A second step is to solve Eq 5 16 for a well chosen initial data from t tn to t tn 1 e g 10 2 t x V x WO t x BFW t x x W t 20 VD ta Xx TE ta41 x l This equation that can be solved exactly has for solution VO t x BEA tn41 1 tn iV x ttn yy 1 2 tn41 xX 5 22 This finally gives Y x WO t 41 x Let us remark that we have to compute the exponential of a matrix to effectively calculate 5 22 For the full approximation we adapt to each component the spectral approximation based on fft and ifft in space For both the Lie and Strang schemes the spectral approximation is written under a symmetrical form as for the one component case 5 2 Relaxation pseudo SPectral scheme ReSP Introduced by Besse in the relaxation method is a scheme that looks like the Crank Nicolson scheme but avoiding the resolution of the nonlinear term through a fixed point or a Newton Raphson method 5 2 1 A Relaxation pseudo SPectral scheme ReSP for the one component case If we consider the 2d cubic GPE te 5Ap t x V x tb t x Blw t 0l24 t x OL xER2 t gt 0 5 23 then the relaxation scheme is gnti 2
86. ections respectively We emphasize here on the fact that these are not equal to Ny Ny and N which designate the total number of grid points including the boundary points For the FD scheme the number of dof is Ny Nz 2 in the x direction and Nz Nx 1 for the SP scheme In example page Ns 28 1 but the number of dof is N 28 which optimizes FFTs computations e Ne is the number of components in GPE Furthermore let us consider the different sets of variables below that must be used when con sidering the corresponding Matlab variables in GPELab e N denotes the positive integers e R designates the real numbers e R R 0 is the set of strictly positive real numbers e C denotes the set of complex numbers We also need the set of strings of characters that we designate by S and the set of Matlab structures denoted by S We now introduce K xK jand L x A 0 where K and L are two sets of variables like the ones defined above In the sequel we use the following notations e F K L is the set of Matlab functions f from K gt L of the form f x1 x2 xN gt 0 x1 x2 xN x1 x2 xMN where 21 12 1N K x K x x Ky More generally we also sometimes use the notation F AK L F K K L if K is repeated p times e Mym K designates a N x M Matlab matrix with values in K for N and M eN e Cym K isa N x M Matlab cell array with values in K N M N Let u
87. eme e BESP N 0 is a variable corresponding to the type of method used to invert the linear system in BESP see section 3 5 1 page 41 e Solver_FD N 0 is a variable corresponding to the type of method used to invert the linear system in BEFD see section 3 5 1 page 41 e Iterative_tol IRF 1e 9 is a variable corresponding to the stopping criterion related to the difference between two successive iterates in the Krylov solver e Iterative_maxit N 1e3 is a variable corresponding to the stopping criterion related to the maximum number of iterations in the Krylov solver For example we want to compute the dynamic solution for a single component BEC by using a splitting scheme We choose a time step At 1073 a stopping time T 1 and the Laplace preconditioner Moreover we choose to compute outputs during the simulation This gives the code in table We now want to set the physical problem This is done exactly like in section 3 6 page 43 Here we list some of the physical operators that are used only in a dynamic problem 5 5 2 The TimePotential_Var2d function In the case where one wants to compute the dynamics of GPE GPELab offers the possibility of handling a time dependent potential The TimePotential_Var2d function allows to define the time dependent potential operator i e V t x in the problem by modifying the Physics2D structure It must be provided with the Method and Physics2D structures and takes the following
88. eo is a variable corresponding to the name of the video of each component If a single character string is provided the names of the videos are for each j 1 Ne VideoName followed by j If a cell vector of character strings is given the names of the videos are for each j 1 Ne VideoName j e Figure S Figure_Var2d is the figure structure If one wants to make a video of the angle of a wave function after a simulation then one should proceed as in table Function phi X Y angle phi MakeVideo2d Method Geometry2D Outputs Function Table 3 43 An example to use the MakeVideo2d function 62 CHAPTER 3 HOW TO USE GPELAB STATIONARY SOLUTIONS Chapter 4 Examples of simulations for stationary solutions This Section is devoted to different physical examples treated by GPELab for computing station ary states The user has a direct access to the source codes in the subdirectories 1d examples 2d examples and 3d examples of GPELab We will be happy to add some of your examples that are different from the ones given here 4 1 Ground state of a 1d Gross Pitaevskii equation with an optical potential and a cubic nonlinearity We now show how to compute the ground state of a Gross Pitaevskii equation with an optical potential and a cubic nonlinearity in 1d This example comes from 21 First we have to build the Method and Geometry1D structures We use the BESP scheme and a time step equal t
89. esponds to the Omega parameter in the Physics2D structure see 3 6 1 page 44 Therefore a continuation on 2 where one wants to increase it from 1 to 5 for example could be defined like in Table where we have set the evolution of the parameter as 2 1 2 3 3 5 4 4 5 5 Continuation Continuation_Var2d Omega 1 2 3 3 5 4 4 5 5 Table 3 31 An example for using the Continuation_Var2d function when the parameter is a scalar Another interesting case is when there is a coupling between multiple components In this scenario a continuation method could be interesting to smoothly couple the two components For example if the two components have coupled nonlinearities one can use a continuation method on the G matrix of the nonlinearity to increment the extradiagonal terms This is done in Table where we have defined a cell vector of G matrices with increasing extradiagonal terms from 0 to 1 in 5 steps 58 CHAPTER 3 HOW TO USE GPELAB STATIONARY SOLUTIONS GNL1 1 0 0 1 GNL2 0 1 1 0 for i 1 5 GNL i GNL1 i 5 GNL2 end Continuation Continuation_Var2d GNonlinearity GNL Table 3 32 An example for using the Continuation_Var2d function when the parameter is a matrix Finally let us consider for example the case of a coupled fast rotating two components Bose Einstein condensate We can proceed by using a continuation simultaneously on the G matrix for the coupling and on
90. fied Program or a work based on the Program To propagate a work means to do anything with it that without permission would make you directly or secondarily liable for infringement under applicable copyright law except executing it on a computer or modifying a private copy Propagation includes copying distribution with or without modification making available to the public and in some countries other activities as well To convey a work means any kind of propagation that enables other parties to make or receive copies Mere interaction with a user through a computer network with no transfer of a copy is not conveying An interactive user interface displays Appropriate Legal Notices to the extent that it includes a convenient and prominently visible feature that 1 displays an appropriate copyright notice and 111 2 tells the user that there is no warranty for the work except to the extent that warranties are provided that licensees may convey the work under this License and how to view a copy of this License If the interface presents a list of user commands or options such as a menu a prominent item in the list meets this criterion 1 Source Code The source code for a work means the preferred form of the work for making modifications to it Object code means any non source form of a work A Standard Interface means an interface that either is an official standard defined by a recog nized st
91. g the evolution of the solition center We obtain Figure 5 1 a We can also plot the evolution of the energy and the chemical potential which are conserved quantities and can be seen on Figure 5 1 b 5 6 2 Dynamic of a dark soliton in a Bose Einstein condensate in 2D We now show how to compute the dynamic of a dark soliton for Gross Pitaevskii equation with cubic nonlinearity in 2D with a single component We first compute a ground state for the Gross Pitaevskii equation with quadratic potential and cubic nonlinearity which corresponds to a stationary Bose Einstein condensate We have to build the Method and Geometry2D structures Here we choose to use the BESP scheme to compute a stationary state We fix the time step At 107 and the stopping criterion e 1078 The geometry is set on a computational domain O 10 10 x 10 10 and the number of grid points J K 2 This is coded in Table 5 13 Now we have to define the physical problem In this case we want to compute a stationary state corresponding to a Bose Einstein condensate in a quadratic trap This corresponds to the 5 6 EXAMPLES OF COMPUTATIONS 101 ES Energy Chemical Potential Center of the soliton N 0 0 1 0 2 0 3 0 4 0 6 0 7 0 8 0 9 A 0 0 1 0 2 0 3 0 4 0 6 0 7 0 8 0 9 1 0 5 Time 0 5 Time a Evolution of the center of the soliton b Evolution of the energy of the soliton Figure 5 1 Some
92. he spatial domain is always rectangular with a uniform mesh grid The output is the Geometry2D structure As for the Method_Var2d function this function includes default values for the input arguments The optional arguments are the following e xmin R 10 is a variable corresponding to the left endpoint of the computational domain in the x direction e xmax R 10 is a variable corresponding to the right endpoint of the computational domain in the x direction e ymin R 10 is a variable corresponding to the lower endpoint of the computational domain in the y direction e ymax IR 10 is a variable corresponding to the upper endpoint of the computational domain in the y direction e Nx N 277 1 is a variable corresponding to the number of points in the x direction e Ny N 277 1 is a variable corresponding to the number of points in the y direction In the case of a 1d simulation one has to discard ymin ymax and Ny Moreover in the case of a 3d simulation one must add zmin and zmax after ymax and Nz after Ny If one wants to take a larger computational box 15 15 x 15 15 and a large number of grid points Nz Ny 2 1 to use a spectral scheme then one builds the Geometry2D structure as in table xmin 15 xmax 15 ymin 15 ymax 15 Nx 279 1 Ny 279 1 Geometry2D Geometry2D_Var2d xmin xmax ymin ymax Nx Ny Table 3 15 An example of how to use the Geometry2D_Var2d func
93. he work is specifically designed to require such as by intimate data communication or control flow between those subprograms and other parts of the work The Corresponding Source need not include anything that users can regenerate automatically from other parts of the Corresponding Source The Corresponding Source for a work in source code form is that same work 2 Basic Permissions All rights granted under this License are granted for the term of copyright on the Program and are irrevocable provided the stated conditions are met This License explicitly affirms your unlimited permission to run the unmodified Program The output from running a covered work is covered by this License only if the output given its content constitutes a covered work This License acknowledges your rights of fair use or other equivalent as provided by copyright law You may make run and propagate covered works that you do not convey without conditions so long as your license otherwise remains in force You may convey covered works to others for the sole purpose of having them make modifications exclusively for you or provide you with facilities for running those works provided that you comply with the terms of this License in conveying all material for which you do not control copyright Those thus making or running the covered works for you must do so exclusively on your behalf under your direction and control on terms that 112 APPENDIX B LICENSE prohi
94. hms are used to effectively solve them In addition we discuss and detail all the functions available in the code Finally complete examples are given in such a way that any user can see how easy it is to manipulate it GPELab is developed to be installed on any computer that has a basic version of Matlab The methods are based on pseudospectral approximation techniques and therefore provide highly accurate solutions The structure of the document is the following Section of Chapter 1 shows which kinds of equations are solved and which dimensionless form of the GPE is used We then discuss rapidly the classes of problems that GPELab can solve A few remarks end Chapter Chapter 2Jexplains the numerical techniques that are coded in GPELab for stationary solutions Mainly the method is based on an imaginary time formulation of the problem The discretization of the resulting equation uses a semi implicit Backward Euler BE time approximation and a SPectral scheme in space based on Fast Fourier Transforms FFTs The method called BESP is extended to the case of multi components systems and high rotations values can be considered The nonlinearities can be classical like for the cubic case but can also be more complicate when considering dipole dipole problems that handle a convolution kernel Other methods included in GPELab are based on the Crank Nicolson approximation scheme in time and the standard 3 points and 5 points second order finite diffe
95. i 2 2 NL 2 2 Phi x y Beta_22 abs Phi 2 2 Beta_12 abs Phit1 72 NL 1 2 Phi x y 0 NL 2 1 Phi x y 0 end Beta_11 2 Beta_12 1 Beta_22 2 Physics2D Nonlinearity_Var2d Method Physics2D Example_nonlinearity Beta_11 Beta_22 Beta_12 Table 3 21 An example to use the Nonlinearity_Var2d function 3 6 4 The FFTNonlinearity_Var2d function Physics2D FFTNonlinearity_Var2d Method Physics2D FFTNonlinearity G FFTNonlinearity_energy Figure 3 2 The FFTNonlinearity_Var2d function In the case where one wants to consider a nonlocal integral type nonlinearity in the GPE this can be done in GPELab by using the Fourier transform The aim is to be able to define general nonlinear interactions through the representation formula Fin W t x 0 y ifft Kj n x E fft Hj W t x for j k 1 Ne For each j k the symbol Kj x is a Matlab function in F My 1 RJ My 1 C My C Each function H is an element of F Cy n Mn Nx C My C Here above fft and ifft designate the FFT and inverse FFT in Matlab respectively The symbol x is the point wise multiplication of arrays in Matlab For example if one wants to compute a convolution form like in a one component dipolar gaz E t x 2 y K x Y t x then you can use the relation F u t x 2 y ifft K t d t x where K is the Fourier transform
96. imensional case is direct Other properties are related to the fact that these schemes are time reversible mass conserving if it is true at the continuous level time transverse invariant and the dispersion relation holds However there is no energy conservation Indeed it can be proved that only the energy of a closely related system is conserved These explicit schemes are unconditionally stable More details can be found in 5 1 3 Extension to the multi components case We consider the notations of Section 2 6 page 30 The TSSP schemes can also be extended to the multi components case i e a system of Ne coupled GPEs 1 Lae iW t x 7 AY t x V x W t x 5 G x aj Or t x j l 6F U t x x U t x t gt 0 x Rf 5 14 where we set W t x a Um t x 1 and x x 21 E j 1 Tj 1 a The operators are defined in Section We furthermore assume that V and F are symmetric operators to get the mass conservation property Vem Vm and Fem Fme 1 lt Lm lt Ne Let us remark that for the TSSP scheme we assume that F only depends on the modulus of Y e g F W t x x F W t x x The main difference is related to the form of the variable coefficients matrices in front of the gradients GA Chey Gly x a5 Gi x 2 Ga Aj Gy x zj Gay AZ FLA A F 6 Cy Az Chatin Gh Ay The initial data is Y t 0 x Vo x To derive the TSSP scheme for this system we choose the same kin
97. in time and spectrally accurate in space It computational cost is O M log M The extension to the three dimensional case and other nonlinearities is direct Other properties are related to the fact that it is time reversible mass and energy for a cubic nonlinearity conserving when the property holds at the continuous level It is not time transverse invariant and the dispersion relation is not preserved The scheme is unconditionally stable We refer to for more details 5 2 2 Extension of the ReSP scheme to the multi components case The relaxation scheme can be extended to the multi components situation in a similar way as the CNGF 2 6 We consider the same notations as before see also page 87 system 5 14 We have the following time discretization of system 5 14 based on the relaxation scheme for a general nonlinearity gn 1 2 gr 1 2 2 prti yn 1 d yn 1 y a A AA NG x 3 O71 y R j 1 BF U x 1 lt n lt M x Rf for n gt 0 Concerning the spatial discretization we use again a pseudo spectral method based on the FFTs iFFTs Let us detail the two dimensional case e g d 2 For the relaxation scheme we have mRen 1 2 2BF W x MRen 1 2 ARengrtl Ren gyn 5 27 where BW 1i Xj k 000 1 gt WN Xie G h 07 x is the discrete unknown array in qxe with M Jk The nonlinear operator MRen 1 2 My y C is updated through the matrix
98. ing to use a splitting scheme see section 5 1 page 83 or Relaxation to use the relaxation scheme see section 5 2 page e Deltat I R 1e 3 is a variable corresponding to the time step of the method The time discretization is always uniform 5 5 GPELAB FUNCTIONS FOR THE DYNAMICS 95 e Stop_time Rt 1 is a variable corresponding to the final time of computation in the case of a dynamic problem This corresponds to the following stopping criterion Iter 0 while Deltat Iter lt Stop_time Compute the solution of the dynamic problem at iteration Iter for Deltat Iter Iter 1 end e Stop_crit R 1e 8 is a variable corresponding to the stopping criterion in the case of a ground state computation see section 3 5 1 page 41 e Max_iter N 1e6 is a variable corresponding to the maximum number of iterations for a stationary state computation e Preconditioner S ThomasFermi is a variable that must be either None for a calculation without preconditioner Laplace for the Laplace preconditioner and ThomasFermi for the Thomas Fermi preconditioner e Output N 1 is a variable that must either be 1 if one computes outputs during the compu tations or 0 if not e Splitting S Strang is a variable that must be either Lie to use the Lie type of splitting Strang to use the Strang type of splitting or Fourth to use the Fourth order type of splitting for the splitting sch
99. ing the Physics2D structure It must be provided with the Method and Physics2D structures and takes the following optional arguments e StochasticPotential Ifa function StochasticPotential in F R M n n R Mw n C is provided the physical time dependent potential will be defined as follows for each j k ls Ney StochasticPotential w t 1 y if j k Wapeni aa w t x y if j If StochasticPotential is a cell array of functions in Cn N F R My Na Ry MN Na C then the potential will be defined by Vir w t x y StochasticPotential j k w t x y for j k 1 Ne The default argument is quadratic_potential2d which corresponds to 1 2 DE er si HO Lig Vir w t y y 0 if j k e G Mn n C ones N_c is a complex valued variable that is used to multiply each compo nent of the potential Vi w t jx y G j k StochasticPotential j k ww t x y for j k 1 Ne 98CHAPTER 5 COMPUTATION OF THE DYNAMICS OF GPE METHODS FUNCTIONS EXAMPLES e StochasticProcess is a function or a cell vector of functions that is used as the stochas tic process w t when computing the StochasticPotential If a function is defined the StochasticPotential will be computed using a scalar value w t If a cell vector of functions is defined the StochasticPotential will be computed using a cell vector of scalar value Wj t For example we want to set a stochastic quadratic potential We would like to set V
100. l matrix with diagonal terms 7 j 4 Pr For the sake of conciseness we denote by p the discrete 2 norm of a vector In the finite difference context the norm of a vector is simply defined by lello h hP D loj 2 29 H ED IK The corresponding function is named L2_norm2d Now we come to the main function BEFD_CNGF2d which solves the ground states computation problem by the BEFD scheme Essentially the description of the input parameters of the function follows the physics and the discretization Furthermore one can add drawing option as well as 26 CHAPTER 2 COMPUTATION OF STATIONARY STATES FOR THE GPE solver choice This last parameter proposes to solve the linear system by a direct gaussian solver or by a preconditioned restarted GMRES Krylov subspace solver The preconditioner is based on the TF approximation related to the diagonal matrix 1 MEF Ai 1 V ello 2 30 The outputs are e the ground state solution at the final time according to the stopping criterion lo loo lt Eat 2 31 wthere the tolerance is e 1078 for example and where is the infinity norm of the function defined at the discrete level by J oo SUPG ED x lj e Time is the discrete vector of time steps 0 the value of the square of the ground state at the origin used in physics papers e Tms ANd Yrms the radius mean square in the z and y directions according to Orms lag llo
101. mation of the Laplacian operator A MANG 119219 1197119 2 42 jk The discrete Laplace operator A is diagonal in the Fourier space but not L Finally the discrete o norm is given by ve E CMP lollo PRP SY lo 2 43 EP 5 K In practice the linear system is efficiently solved by a Krylov solver BiCGStab precon ditioned by the TF operator AREN see In GPELab the function for computing the ground state by BESP is BESP_CNGF2d for the two dimensional case The operator A applied to a vector is given in Delta_Fourier2d the gradient operator is Grad_Fourier2d Other Matlab files can be found in the directory Code2D FFT but are nevertheless not useful when you use GPELab They are transparent for the user and do not really need to be explained They generally contain helpful functions for optimizing the BESP method for instance preconditioners BESP is the default method used to solve the CNGF in GPELab see Section 3 5 1 page 41 The output physical quantities are the same as these provided by the BEFD scheme see subsection 2 3 2 page pa 2 4 CRANK NICOLSON SCHEMES 29 2 4 Crank Nicolson schemes Another possibility is to use a semi implicit Crank Nicolson scheme 16 This results in the dis cretization Sek an Mn Mn ANG Nn nti 2 44 e iek ACNm BON are such that where and I 1 1 1 1 ANT A4 n2 OL 5 Gar VU shel 5 i 2 45 I 1 1 1 1
102. mputations For this reason we will not use this scheme for the stationary states but it will be used for the dynamics Here we rather consider the BEFD scheme with rotating term The scheme is implicit and therefore it requires at each step the solution to a linear system It however can be solved efficiently by using a direct solver or a preconditioned Krylov subspace methods e g a Bi Conjugate Gradient Stabilized BiCGStab The interesting property is that the scheme is energy diminishing for the non rotating case and larger time steps can be used The BEFD scheme is however only second order accurate in space which is a limitation for computing fast rotating condensates Higher order schemes could be used This is the goal of subsection 2 3 3 where we present a new scheme called BESP based on BE in time but on a SPectral FFT scheme in space to capture accurately the creation of vortices for fast rotating condensates The BESP is the scheme that you should prefer to use in GPELab when you consider fast rotating gazes 2 3 2 Backward Euler Finite Difference BEFD Concerning the time discretization of 2 22 the application of the Backward Euler scheme leads to the semi discrete semi implicit linear scheme BE scheme 228 lag v g peP ALg 1 lt n lt N xeR A 2 2 23 n 1 _ o 2 l XERS lello setting M At Tnax where Thrax is the maximal time of computation and N is the number of time steps Let us remark here that
103. ne the previous nonlocal nonlinearity which is then added to the Physics2D structure by using the FFTNonlinearity_Var2d m function function Dipolar2d Example_fftnonlinearity phi x y fftx ffty square_xi fftx 2 ffty 2 K 2 3 sqrt pi sqrt square_xi exp square_xi erfc sqrt square_xi Dipolar2d ifft2 K fft2 abs phi 72 end Physics2D FFTNonlinearity_Var2d Method Physics2D phi x y fftx ffty Example_nonlinearity phi x y fftx ffty Table 3 22 An example of application of the FFTNonlinearity_Var2d function with the dipolar operator http www mathworks fr fr help matlab ref erfc html 3 6 SETTING THE PHYSICAL PROBLEM 51 3 6 5 The gradient functions Physics2D Gradientx_Var2d Method Physics2D Gradientx G Table 3 23 The Gradientx_Var2d function The gradient functions allow to define the derivation operators a Gi x 0 in the problem by modifying the Physics2D structure Here we take for example the function Gradientx_Var2d as the other gradient functions work similarly We remark that we can only define Gradientx in 1d Gradientx and Gradienty in 2d and Gradientx Gradienty and Gradientz in 3d The Method and Physics2D structures are required arguments It is possible to include the following optional arguments e Gradientx Let us provide a function Gradientx in F Mn n RJ Mn 1 C then the variable coefficients in front of the gradient are defined as 1
104. nergy is a nonlinear operator used to compute the energy associated to the physical nonlocal nonlinearity i e it corresponds to Fenergy t x x in the energy 3 2 Note that it must be the same type of variable as FFTNonlinearity Again if G is defined it is multiplied element by element by FFTNonlinearity_energy If a function FFTNonlinearity_energy in F Cn Ne M Ny No C MN No Ry MN Na cr MN No C is given the nonlocal nonlinear energy operator is for each j k 1 Ne FFTNonlinearity_energy W t x Y Ex ifj k aS a If FFTNonlinearity_energy is a cell array of function in Cro NF Cn N MNN C Mn Na RJ Mn 12 C Mn 1 C then the nonlocal nonlinear energy operator is defined by Fenergy j k U t x x y FFTNonlinearity_energy j k U t x y Ex Ey for j k 1 Ne The default argument is Cubic_energy2d 1 Des Fenergy j k Y t x 2 Y 0 if j k 50 CHAPTER 3 HOW TO USE GPELAB STATIONARY SOLUTIONS A standard example of nonlocal nonlinearity comes from the dipolar interaction This interaction is usually described in 3d however a reduction to the 2d case is possible The dipolar interaction in the 2d case with a single component is given by see 18 F Y x x F 2 3Vml lelerfe lel F w t xP where erfc is the complementary error functior This nonlinearity can be easily implemented by using a Matlab script This is done in table where we defi
105. ng arguments Delta Beta or Omega see Section 3 6 1 page 44 GPotential where G is the matrix associated to the potential see Section 3 6 2 page T GNonlinearity with G is the matrix related to the nonlinearity see Section page 16 GFFTNonlinearity G is the matrix for the FFT based nonlinearity see Section 3 6 4 page HE or GGradientx or GGradienty if one wants to change the G matrix for the gradients see Section 3 6 5 page 51 e Coefficient Cy ncont C or Ci ner My Ne C is the p variable containing the values of the parameter that will be changed during the continuation method For a scalar parameter Delta Beta or Omega it must be a cell array of scalars In the case of a ma trix parameters GPotential GNonlinearity GFFTNonlinearity GGradientx or GGradienty it must be a cell array of matrices We have the following optional argument e Continuation S the continuation structure must be taken as an argument if one wants to add parameters in the continuation method for multi parameter paths A simple example for the continuation method is to increment a single parameter We know that if one wants to compute the ground state of a fast rotating Bose Einstein condensate one should do this by a continuation method on the rotation speed In GPELab the rotation speed is denoted by 2 which corr
106. ns In our case we identify e the potential function V x y 3 z yl e the nonlinear function F Y x y W t z y e the gradient function in the x direction G x y iQy e the gradient function in the y direction G x y Mz In this particular case the operators are predefined in GPELab but for clarity we define them again in our script To define the physical problem we first need to build the Physics2D structure and set the values of the parameters 6 6 and Q The Physics2D structure contains all the informations related to the physical problem that is among other the functions related to the operators There fore we also have to add the operators by using functions of GPELab to the Physics2D structure We remark that as we said the operators are already predefined and set as default in the functions Potential_Var2d Nonlinearity_Var2d Gradientx_Var2d and Gradienty_Var2d The resulting code in available in Table We now have to set the initial data Initial data in GPELab are defined as a cell array each cell containing a complex matrix which is the initial wave function of a component GPELab users can 3 2 A SIMPLE BUT COMPLETE EXAMPLE 37 Delta 0 5 Beta 300 Omega 0 7 Physics2D Physics2D_Var2d Method Delta Beta Physics2D Potential_Var2d Method Physics2D x y 1 2 x 72 y 72 Physics2D Nonlinearity_Var2d Method Physics2D phi x y abs phi 72 Physics2D Gradie
107. ns corresponds to the number of time steps We also need the Print and Figure structures which are built using the Print_Var2d function see section page 58 and Figure_Var2d function see sections 3 7 4 page and page 59 respectively Finally we use the GPELab2d function to launch the simulation see section 3 7 6 page 60 5 6 Examples of computations 5 6 1 Dynamic of a bright soliton for the Gross Pitaevskii equation with cubic nonlinearity in 1D We now show how to compute the dynamic of a soliton for Gross Pitaevskii equation with cubic nonlinearity in 1D with a single component First we have to build the Method and Geometry1D 5 6 EXAMPLES OF COMPUTATIONS 99 structures We want to use a Splitting scheme and a step time of 1078 with a spatial grid of 21 4 1 points on the interval 15 15 Moreover we want to set the stopping time to 1 All of this is coded in Table Computation Dynamic Ncomponents 1 Type Splitting Deltat ie 3 Stop_time 1 Stop_crit Method Method_Varid Computation Ncomponents Type Deltat Stop_time Stop_crit xmin 15 xmax 15 Nx 2710 1 Geometry1D Geometry1D_Varid xmin xmax Nx Table 5 7 Building the Method and Geometry1D structures We now have to define the physics of the problem We want to compute the soliton for the following Gross Pitaevskii equation O D t 0 AV t 2 BIU t z U E 2 with 8
108. nstant of the Laplacian operator i e 6 in equation 3 1 e Beta R 0 is a variable corresponding to the constant in front of the nonlinearity 6 in equation 3 1 e Omega IR 0 or M1 R 0 is a variable that defines the rotation speed if the default gradient operators 10 y r in 2d or x x V in 3d are present in the equation otherwise it has no effect It is a real valued parameter Q R in the 2d case or a vector Q Mj 3 R in the 3d case We note that this variable does not exist in the 1d situation If the default gradient operators are set then we have the following rotation operators Q 10 y r 0 a 0 Ss 0 Q 10 yOr 0 De G x Oz l i 2 z 0 0 OQ xOy y r in the 2d physical problem and Q x x V 0 vee 0 oe 0 N xxV 0 S Gi x l i E 0 0 sini ERY for the 3d case 3 6 SETTING THE PHYSICAL PROBLEM 45 We remind that in the case where one puts a rotational operator with a rotation speed ranging from 0 to 1 excluded then a quadratic potential V x x is enough to compensate the centrifugal force However for a rotation speed larger than 1 one must use a stronger potential for example a quartic potential V x x In table 3 17 we show the code where we create the Physics2D structure with 6 5 8 1000 and Q 0 7 Delta 0 5 Beta 1000 Omega 0 7 Physics2D Physics2D_Var2d Method Delta Beta Omega
109. ntx_Var2d Method Physics2D x y 1i xOmega y Physics2D Gradienty_Var2d Method Physics2D x y 1i Omega x Table 3 2 An example of how to define the Physics in GPELab through the Physics2D structure set the initial data themselves However the function InitialData_Var2d is helpful if one wants to use standard initial data like the centered gaussian or the Thomas Fermi approximation see Section page 20 We would like in our case to use the Thomas Fermi approximation as an initial wave function for the computation This is done by setting the InitialData_choice variable to 2 as we can see this in Table InitialData_choice 2 Phi_O InitialData_Var2d Method Geometry2D Physics2D InitialData_choice Table 3 3 Initialization by the Thomas Fermi approximation Finally we want to launch the simulation We would like to get informations about the wave function during the computations for example to be sure that the energy is diminishing Thus we need to build the Outputs structure that contains all the outputs computed during the simulation Some outputs like the energy or the mean square radius are already defined and stored when computing Moreover to print these informations during the calculations we have to build the Print structure that explains how to print the outputs For example in Table 3 4 we ask to print the informations every 15 iterations and to draw the solution Outputs Ou
110. o 10 15 20 20 15 10 a5 10S o a Soliton at time t 0 25 b Soliton at time t 0 5 c Soliton at time t 0 75 d Soliton at time t 1 e Soliton at time t 1 25 Soliton at time t 1 5 Figure 5 3 Evolution of a dark soliton in a Bose Einstein condensate Computation Dynamic Ncomponents 1 Type Splitting Deltat 1e 3 Stop_time 1 Method Method_Var2d Computation Ncomponents Type Deltat Stop_time Method Splitting Fourth xmin 10 xmax 10 ymin 10 ymax 10 Nx 278 1 Ny 278 1 Geometry2D Geometry2D_Var2d xmin xmax ymin ymax Nx Ny Table 5 19 Building the Method and Geometry2D structures for the simulation potential We consider the following Gross Piteavskii equation 1 l a K 2 iOpw t x y Avtt z y 9 4 12 21 a ylyl p t x y B w t x y Yl e y iQ yO 20 W t e y with the parameters a 1 2 k 0 7 Yz Yy 1 2 1000 and Q 3 5 Here we have changed the parameter from 0 3 for the stationary to 0 7 for the dynamical problem We define the potential operator and add the defaults nonlinear and gardients operators The resulting code can be seen in Table 5 20 We finally launch the simulation We set the defaults Outputs that we compute each 100 iterations and store the solution using the OutputsINI_Var2d function Moreover we build the 21 yylyl 5 6 EXAMPLES OF COMPUT
111. o 5 107 with a spatial grid of 21 1 points in the interval 16 16 Moreover we set the stopping criterion to 1078 We refer to Table 4 1 below for the corresponding GPELab code Computation Ground Ncomponents 1 Type BESP Deltat 5e 2 Stop_time Stop_crit le 8 Method Method_Varid Computation Ncomponents Type Deltat Stop_time Stop_crit xmin 16 xmax 16 Nx 2710 1 Geometry1D Geometry1D_Varid xmin xmax Nx Table 4 1 Defining the Method and Geometry1D structures We now define the physics of the problem to compute the ground state of the following Gross Pitaevskii equation Eli ilt a ZAVE a EL 25sm TE y a Blott 10 63 64 CHAPTER 4 EXAMPLES OF SIMULATIONS FOR STATIONARY SOLUTIONS with 8 250 Let us remark that since we consider the cubic case this is the default nonlinearity in GPELab and we therefore do not have to define it We have e to build the Physics1D structure with the desired coefficients e to define and add the optical potential e and finally to add the default nonlinear operator to the physics of the problem We follow this construction in Table 4 2 where we directly define the optical potential in the argu ments of the Potential_Varid function Delta 0 5 Beta 250 Physics1D Physics1iD_Varid Method Delta Beta Physics1D Potential_Varid Method Physics1D x x 72 2 25 sin pi x 4 2
112. o Bala y OL 1 12 a Te where b2 85 r and 1 Vale y 122 14 1 13 As a conclusion one may write the GPE in d dimensions through o iE LA Va Bale oI QL 1 14 for x Rf t gt 0 63 8 and V3 a y z V a y z 1 3 Which problems can GPELab solve GPELab Gross Pitaevskii Equation Laboratory is a versatile and robust Matlab code that provides highly accurate mainly spectral numerical methods for computing the solutions of Gross Pitevskii equations often used to simulate the physics of Bose Enstein Condensates for superfluids The main features contained in GPELab are 16 CHAPTER 1 INTRODUCTION amp INFORMATIONS e computation of stationary states and dynamics of solutions for GPE e one two and three dimensional problems e general potentials general local and nonlocal dipolar nonlinearities fast rotating gazes e multi components GPE e inclusion of stochastic terms at different places Furthermore Matlab offers different visualization tools to observe and to represent the solutions of your calculations GPELab also provides physical quantities of interest for the user For beginners GPELab is easy to use in many standard situations For more advanced users the flexibility of the code allows you to have your own defined output quantities that are not already contained in the code In terms of performance GPELab is developed by using an optimized implementation in Matlab Fur
113. o redis tribute it under certain conditions type show c for details The hypothetical commands show w and show c should show the appropriate parts of the General Public License Of course your program s commands might be different for a GUI interface you would use an about box You should also get your employer if you work as a programmer or school if any to sign a copyright disclaimer for the program if necessary For more information on this and how to apply and follow the GNU GPL see lt http www gnu org licenses gt The GNU General Public License does not permit incorporating your program into proprietary programs If your program is a subroutine library you may consider it more useful to permit linking proprietary applications with the library If this is what you want to do use the GNU Lesser General Public License instead of this License But first please read http www gnu org philosophy why not lepl html 120 APPENDIX B LICENSE Bibliography 1 F Kh Abdullaev B B Baizakov and V V Konotop Dynamics of a Bose Einstein condensate in optical trap Nonlinearity and Disorder Theory and Applications edited by F Kh Abdullaev O Bang and M P Soerensen NATO Science Series vol 45 Kluwer Dodrecht 2001 2 F Kh Abdullaev J C Bronski and R M Galimzyanov Dynamics of a trapped 2D Bose Einstein condensate with periodically and randomly varying atomic scattering length Physica D
114. od and Geometry2D structures 4 5 GROUND STATE OF 2D GPE SYSTEM WITH COUPLED NONLINEARITIES 75 We now define the physics of the problem Let us consider the problem of computing the ground state of the following system of Gross Pitaevskii equations Onn t 2 y svat ay 5 Izl lvl Bila t z y Bualdatt 2 0 eae 10 yO xy Y t Y Onbalt x y 5 Ayalt 2 y G jal yl Bolbe t x y bralda t 2 9 palt x y 10 yO xy alt x y with 61 2 400 B12 200 and Q 0 8 We recall that since the default potential is the quadratic potential and the default gradients define the rotational operator we only have to define the coupled cubic nonlinearity This is done in Table function CoupledCubicNonlinearity Coupled_Cubic2d Beta_coupled CoupledCubicNonlinearity cel1 2 CoupledCubicNonlinearity 1 1 Phi X Y Beta_coupled 1 1 abs Phit1 72 Beta_coupled 1 2 abs Phi 2 2 CoupledCubicNonlinearity 2 2 Phi X Y Beta_coupled 2 2 abs Phit2 72 Beta_coupled 2 1 abs Phi 1 2 CoupledCubicNonlinearity 1 2 CoupledCubicNonlinearity 2 1 Phi X Y 0 Phi X Y 0 Table 4 24 Defining the coupled nonlinearity We also have to impose the energy associated to the coupled cubic nonlinearity Table 4 25 function CoupledCubicEnergy Coupled_Cubic_energy2d Method Beta_coupled CoupledCubicEnergy cel1 2 CoupledCubicEnergy
115. oding Each discrete x and y derivative uses the two points scheme adapted to the homogeneous Dirichlet boundary condition Ft k Pik Piety Pj k 1 2hy Pr Nig by Pik 2 25 Concerning the GPELab toolbox the discrete rotational operator can be found in the Two_Points_Rotation2d function The derivative operators are the Two_Points_Gradientx2d and Two_Points_Gradienty2d functions These three operators provide square matrices in M m C with respect to the discretization with Mp J 1 K 1 The Laplacian is discretized thanks to the five points scheme with homogeneous Dirichlet bound ary conditions The interior scheme is based on n ES n n eg Pitta Oik TO k a xr jk h2 y y FG k fe Joy 5 eT OR j k Dix Phar 20 Gr Oik 2 h2 2 26 The corresponding function for A is named FivePoints_Laplacian2d It provides the matrix A Myp C which must be applied to a global vector q ir NED of size Mp Finally the potential is only considered at the discretization points and leads to a diagonal matrix V Mm C with diagonal elements V y 1 As a consequence the spatial discretization of leads to the Mp x Mp linear system with normalization step Alo b grt e 2 27 IIPllo with i i A gl lAl V elg 0 L 2 28 b ae At Hereabove we set I Myp C as the identity matrix and Mmp C as the diagonal potentia
116. on of 5 27 we again use the preconditioned BiCGStab Concern ing the TF like preconditioner since we have some coupling effects between the gazes through V and MRen 1 2 the matrix Azer is non diagonal We propose here to only keep the diagonal part of Aa for preconditioning that is to include the potential and nonlinear self interactions in each gas Concretely we build the following diagonal TF preconditioner Pre vee given by si Re n In 1 1 Re n 1 2 PTF diag At 5 V diag 5 Maiag Re n 1 2 Re n 1 2 Mee e 1 Ne where Vaiag Vee e 1 N and Maiag We can also build the full TF preconditioner pren by including the full potential and nonlinear operators and inverting the matrix This preconditioner is therefore given by pRen Me ly gRent1 2 TF full At 2 2 Finally we have the Laplace preconditioner pre which is built by inverting the laplacian which is a diagonal operator in the frequencies space by using a fast Fourier transform It is given by pRe iz man 92CHAPTER 5 COMPUTATION OF THE DYNAMICS OF GPE METHODS FUNCTIONS EXAMPLES 5 3 Integration of a stochastic potential In this section we consider a Gross Pitaevskii equation with a stochastic potential We are interested on how to discretize in time equations of the form i0 w t x 5 A t x V w t x v t x Bl t x v t x OL x R t gt 0 where we define the stochastic potenti
117. onlinearity ifftn V NLFFT Table 4 34 Defining the dipolar interaction via a FFT Delta 0 5 Beta 1000 Dipolar_direction 0 0 1 d 0 5 Physics3D Physics3D_Var3d Method Delta Beta Physics3D Potential_Var3d Method Physics3D Physics3D Nonlinearity_Var3d Method Physics2D Physics3D FFTNonlinearity_Var3d Method Physics3D Phi X Y Z FFTX FFTY FFTZ Dipolar_interaction3d Phi FFTX FFTY FFTZ Dipolar_direction d Table 4 35 Constructing the physics of the problem InitialData_choice 2 Phi_O InitialData_Var3d Method Geometry3D Physics3D InitialData_choice Table 4 36 Choosing the initial data 4 7 GROUND STATE OF 3D GPE EQUATION WITH DIPOLE DIPOLE INTERACTION 81 We finally set the outputs and the printing informations then we launch the simulation see Table E37 Outputs DutputsINI_Var3d Method Printing 1 Evo 15 Draw 1 Print Print_Var3d Printing Evo Draw Phi Outputs GPELab3d Phi_0 Method Geometry3D Physics3D Outputs Print Table 4 37 Initializing the outputs setting the informations to print and launching the computation At the end of the computation we obtain see Figure 4 11 the isovalues of the modulus of the stationary state solution Iphi x y 2 12 of component 1 Figure 4 11 107 isovalues of the modulus of the ground state 82 CHAPTER 4 EXAMPLES OF SIMULATIO
118. optional arguments 96 CHAPTER 5 COMPUTATION OF THE DYNAMICS OF GPE METHODS FUNCTIONS EXAMPLES Computation Dynamic Ncomponents 1 Type Splitting Deltat ie 3 Stop_time 1 Stop_crit Max_iter Precond_type Laplace Output 1 Method Method_Var2d Computation Ncomponents Type Deltat Stop_time Stop_crit Max_iter Precond_type Output Table 5 2 An example of initialization and use of the Method_Var2d function Physics2D TimePotential_Var2d Method Physics2D TimePotential G Table 5 3 The TimePotential_Var2d function e TimePotential If a function TimePotential in F R My w R My C is provided the physical time dependent potential will be defined as follows for each j k 1 Ne TimePotential t x y if j k Vatan f Dijk If TimePotential is a cell array of functions in Cn N F R My na R Mn 01 C then the potential will be defined by Vi k t 2 y TimePotential j k t x y for j k 1 Ne The default argument is quadratic_potential2d which corresponds to 1 2 Divs ee 7 5 a7 ifj k Viet 2 y a re a Note that in the case of a stationary state computation the potential operator should be time independent e G Mn n C ones N_c is a complex valued variable that is used to multiply each compo nent of the potential Virlt 1 y G j k TimePotential j k t x y for j k
119. or computing dynamics of rotating Bose Einstein condensates J of Comp Physics Vol 217 pp 612 626 2006 13 J Lee and B Fornberg Some unconditionally stable time stepping methods for the 3 D Maxwell s equations J of Comp and Applied Mathematics Vol 166 Issue 2 pp 497 523 2004 14 E Faou Analysis of splitting methods for reaction diffusion problems using stochastic calculus Math Comp Vol 78 pp 1467 1483 2009 121 122 BIBLIOGRAPHY 15 W Bao H Wang and P Markowich Ground symmetric and central vortex states in rotating Bose Einstein condensates Comm Math Sci 3 1 pp 57 88 2005 16 W Bao and Q Du Computing the ground state solution of Bose Einstein condensates by a normalized gradient flow SIAM J Sci Comput 25 5 pp 1674 1697 2004 17 W Bao and Y Cai Ground States of Two component Bose Einstein Condensates with an Internal Atomic Josephson Junction East Asia Journal on Applied Mathematics Vol 1 pp 49 81 2010 18 P Pedri and L Santos Phys Rev Lett 95 200404 2005 19 R Zeng and Y Zhang Efficiently computing vortex lattices in rapid rotating Bose Einstein condensates Comput Phys Comm 180 pp 854 860 2009 20 W Bao and H Wang An efficient and spectrally accurate numerical method for computing dynamics of rotating Bose Einstein condensates J Comp Phys 217 pp 612 626 2006 21 W Bao I L Chern and F Y Lim Efficient and spectrally accur
120. or not to save the computed wave functions in the output structure every Evo It must be either 1 if one saves the wave functions or 0 otherwise e userdef_outputs is a cell array of functions in Cy tout F Mn Nz MN No R M Ny Na oe R that allows the user to define itself n relevant physical output quantities These quantities are computed through n t Matlab functions that the user must write himself under the form welt x ZY En Ey userdef_outputs j e t x ZY Ez Ey for j 1 n where and y are the discrete Fourier frequencies in the z and y directions We remark that userdef_outputs must have Y t xX Y En Ey as arguments only in the case where a spectral scheme is used Otherwise the arguments are q t x y By default there is no other output computed than the predefined ones e userdef_outputs_names Cj nout S User defined function is a cell array of character strings where the 7 th component corresponds to the name displayed in the command window of the j th physical quantity appearing in userdef_outputs e globaluserdef_outputs is a cell array of functions in C1 nco F C1 N LM No CT MN Na Ry MN No Cy R that defines n relevant physical output quantities We remark that compared with the previous variable userdef_outputs these physical quantities can be defined through expres sions involving the full wave function Y and not only its on
121. orm of 4 Another quantity which must be conserved is the energy Eg that is given by for a cubic nonlinearity here and a time independent potential Bats 5IVul Veola Slit ayt ax 2 4 2 1 2 Stationary states One important question in the numerical solution of GPE is the computation of stationary states The problem consists in finding a solution Y t x e G x 2 5 where y is called the chemical potential of the condensate and is a time independent function This solution is given as the solution to the nonlinear elliptic equation 1 ud x 5AG x V x x ble Ax QL246 x 2 6 under the normalization constraint lol f loborax 1 2 7 R2 This nonlinear eigenvalue problem can be solved by computing the chemical potential p Hg alp Esol ale o x ax 2 8 with i Eaa d VA Vio Blof 06 L dx 2 9 This also means that the eigenfunctions are the critical points of the energy functional Eg n over the unit sphere S o 1 Furthermore can be seen as the Euler Lagrange equations associated with the constraint minimization problem 2 7 Computing the global minimal solutions y to the energy functional under the normalization constraint bg argmin Eg olo 2 10 pes corresponds to obtain a ground state solution while local minima are excited metastable states 2 2 Approximate solutions as initial guess and potentials When one wants to compute numerically sol
122. otting the evolution of x y rms 106CHAPTER 5 COMPUTATION OF THE DYNAMICS OF GPE METHODS FUNCTIONS EXAMPLES Iphi x y I of component 1 x10 Iphi x y I of component 1 Iphi x y l of component 1 2 12 2 x lo 02 a BEC at time t 0 b BEC at time t 0 12 c BEC at time t 0 23 Iphicey I of component 1 Iphitc I of component 1 0 028 k 0 035 a 0 08 0 02 8 4 0 025 0 015 2 0 02 ao e 2 1 018 4 0 01 0 005 6 0 005 a lo o 10 8 6 4 2 0 2 4 6 8 w Iphi ey 1 of component 1 10 6 2 10 0 022 a a 0 02 6 0018 4 4 ois ora 2 ore gt 0 0 01 A 0 008 4 10 008 6 0 004 10 002 e o 2 4 6 a o o 0 S 6 4 2 0 2 4 6 8 10 10 8 6 4 2 0 2 4 6 s w d BEC at time t 0 34 e BEC at time t 0 45 f BEC at time t 0 56 Figure 5 4 Evolution of a fast rotating Bose Einstein condensate when changing the intensity of the potential 48 r r T r r r T r r 4 8 x rms y rms 01 02 03 04 05 06 07 08 09 1 0 oi 02 03 04 05 06 07 08 09 Time Time a Evolution of the mean root square in the x b Evolution of the mean root square in the y direction direction Figure 5 5 Evolution of the root mean square in the x and y direction Appendix A Copyrights amp credits GPELab is copyright C 2013 Xavier Antoine lt xavier antoine at univ lorraine fr gt and Romain Duboscq lt romain duboscq at univ lorraine fr gt
123. ou have certain responsibilities if you distribute copies of the software or if you modify it responsibilities to respect the freedom of others For example if you distribute copies of such a program whether gratis or for a fee you must pass on to the recipients the same freedoms that you received You must make sure that they too receive or can get the source code And you must show them these terms so they know their rights Developers that use the GNU GPL protect your rights with two steps 1 assert copyright on the software and 2 offer you this License giving you legal permission to copy distribute and or modify it 109 110 APPENDIX B LICENSE For the developers and authors protection the GPL clearly explains that there is no warranty for this free software For both users and authors sake the GPL requires that modified versions be marked as changed so that their problems will not be attributed erroneously to authors of previous versions Some devices are designed to deny users access to install or run modified versions of the software inside them although the manufacturer can do so This is fundamentally incompatible with the aim of protecting users freedom to change the software The systematic pattern of such abuse occurs in the area of products for individuals to use which is precisely where it is most unacceptable Therefore we have designed this version of the GPL to prohibit the practice for those prod
124. pectively Higher order schemes t can be constructed for appropriately chosen weights az bk i lt k lt p The main idea behind the splitting of operators is that each equation associated with each operator A or B can be solved easily 5 1 1 The Lie ADI TSSP scheme for the rotating GPE The Lie TSSP scheme Let us first begin by the Lie splitting scheme Let us recall the rotating cubic GPE equation in 2d 1 i t x 534V x V t x v t x blYlt x Y t x OL y t gt 0 x R 5 1 In terms of splitting a natural decomposition of the above equation consists in solving 1 the partial differential equation related to A A iQL 2 and next the ODE with respect to the potential and nonlinear terms B iV t x a Blp t x Indeed as we will see later the equation associated to the partial differential operator can be solved with the help of spectral methods Furthermore the equation related to the nonlinearity and the potential parts can be integrated exactly Indeed for a time step At it can be proved that an approximation of the solution is given by pat x AGAHOLJA V tx 6 P tx Atyn x Hence for solving 5 1 we can proceed in two successive steps 1 let At gt 0 n N 0 and nAt lt t lt n 1 At we solve l iOnpr t x Ada t x QL t x NAL lt t lt n 1At ES vi nAt x y x 2 and next alex Vitali Altea MAIS DAS Ga wba nAt x di n
125. pter 2 1 1 The time dependent GPE Our aim is to solve the Gross Pitaevskii Equation GPE with rotating terms We essentially consider the dimensionless equation even if the equation with the usual physical parameters could be used In the two dimensional framework the equation writes down ip t x Ami x V t x y t x f y t x y t x OL p t x t x RFXR 2 1 where 7 is the condensate wave function and R 0 00 The Laplace operator is defined as A V Function V is the external usually trapping potential for example harmonic Parameter B is the nonlinearity strenght which describes the interaction between atoms in the condensate It is related to the s scattering length a It is positive for a repulsive interaction and negative for attractive interactions Function f describes the nonlinearity arising in the problem which is usually cubic f W but other cases will be considered later For vortices creation a rotating term is added The operator L is defined as the z component of the angular momentum L Pzr Py pz xA P with the momentum operator P iV In the two dimensional case and after some manipulations its expression is L i x y y r 2 2 Two invariants must be satisfied by w after normalization The mass must be conserved NW f be tx Pax f Oxa Iyl 1 gt 0 2 3 19 20 CHAPTER 2 COMPUTATION OF STATIONARY STATES FOR THE GPE where 7 o is the Z R n
126. putation Ground Ncomponents 1 Type BESP Deltat 1le 2 Stop_time Stop_crit 1e 6 Method Method_Var2d Computation Ncomponents Type Deltat Stop_time Stop_crit xmin 15 xmax 15 ymin 15 ymax 15 Nx 278 1 Ny 278 1 Geometry2D Geometry2D_Var2d xmin xmax ymin ymax Nx Ny Table 4 13 Creating the Method and Geometry1D structures 70 CHAPTER 4 EXAMPLES OF SIMULATIONS FOR STATIONARY SOLUTIONS We now define the physics of the problem Let us consider the computation of the ground state of the following Gross Pitaevskii equation 2 BlwW t 2 y Y t e y 1Op t z y sav peie i r sin sin w t x y with x 100 and 6 1000 Let us recall that since the default nonlinearity is the cubic nonlinear ity we do not have to redefine it Therefore we only have to impose the optical potential and then to build the Physics2D structure We end by adding the potential operator to Physics2D which is defined as the optical potential and the default nonlinear operator see Table 4 14 Delta 0 5 Beta 1000 Kappa 100 Optical_Potential x y 1 2 x 72 y 72 Kappa sin pi x 4 2 sin pix y 4 72 Physics2D Physics2D_Var2d Method Delta Beta Physics2D Potential_Var2d Method Physics2D 0ptical_Potential Physics2D Nonlinearity_Var2d Method Physics2D Table 4 14 Construction of the Physics2D st
127. r applicable patent law 12 No Surrender of Others Freedom If conditions are imposed on you whether by court order agreement or otherwise that con tradict the conditions of this License they do not excuse you from the conditions of this License If you cannot convey a covered work so as to satisfy simultaneously your obligations under this License and any other pertinent obligations then as a consequence you may not convey it at all For example if you agree to terms that obligate you to collect a royalty for further conveying from those to whom you convey the Program the only way you could satisfy both those terms and this License would be to refrain entirely from conveying the Program 13 Use with the GNU Affero General Public License Notwithstanding any other provision of this License you have permission to link or combine any covered work with a work licensed under version 3 of the GNU Affero General Public License into a single combined work and to convey the resulting work The terms of this License will continue to apply to the part which is the covered work but the special requirements of the GNU Affero General Public License section 13 concerning interaction through a network will apply to the combination as such 14 Revised Versions of this License The Free Software Foundation may publish revised and or new versions of the GNU General Public License from time to time Such new versions will be similar in spirit to
128. r is to introduce the way GPELab works when you want to compute stationary solutions and to detail some functions related to this class of problems 3 1 How to get and install GPELab You can freely download the GPELab solver at the following address a mettre cor rectement The installation process is simple You get GPELab from the website save the archive geplab zip on your laptop and unzip the file Then add the GPELab directory and its subdirecto ries in the Matlab setpath menu You can launch Matlab which now knows the correct paths for using the GPELab functions Once its done GPELab can be directly used You can open any of the examples files and test it to check that everything works well 3 2 A simple but complete example We now present how to compute the ground state of a one component Gross Pitaevskii equation with quadratic potential cubic nonlinearity and rotational operator in 2D The following program is an example of how you write a script in Matlab that will launch the computation of a ground state for such a physical configuration The first part of the script consists in building two structures named Method and Geometry2D that will contain all the informations related to the method and the geometry respectively In this example we choose the BESP scheme to compute a ground state Moreover we fix the time step At such that At 107 and the stopping criterion e in 2 31 to 35 36 CHAPTER 3 HOW TO USE GPELAB STATI
129. rd Euler discretization that is BEFD and BESP Using the operators that we have introduced for the one component case the extension is direct even from the point of view of the Krylov solver solution We have the following time discretization of system based on the semi implicit Backward Euler scheme d 6 09 1 5 s P a AD V x 2 x 0 D BF x 1 lt n lt M xeR grt 2 XK R Plo setting MAt Tmax where Tmax is the time of computation to get the solution satisfying the con vergence criterion Integer M is the number of related time steps For the spatial discretization we extend what is done in the one component case to the multi component case by considering the discrete Laplacian and gradients The functions are evaluated pointwise on a rectangular uniform discretization grid according to the space dimension For the Finite Difference respectively SPec tral scheme the resulting method is again called BEFD respectively BESP The semi implicit Crank Nicolson scheme is also implemented resulting in the CNFD and CNSP computational meth ods Let us consider for example the two dimensional problem and let us denote x by x x y We assume that the support of the evolving field of each component is inside a box O Lo Lalx Ly Lyl We consider the uniform grid Osk X k Weis je HLS Ki J and K being two even positive integers Furthermore we introduce the indices of the spatial g
130. re 1 N is a time independent function which is a solution of the following problem d 1 iu x A9 x V x 0 x Y G x 0 x BF P x x x 2 51 j l under the total mass constraint N 1 and where 0 is the chemical potential given by the formula l l x dx x 10 925 fa Ives ya Jonas Like in the one component case we propose to use the CNGF for the multi components problem which is a direct extension Let us denote by t lt lt tn lt the discrete times and by Atn tn 1 tn the local time step The Continuous Normalized Gradient Flow is given by d V x X G x s BF x 2 dx k 1 d 1 A AB Va E 5Ad V x 0 Y G 00 0 BF D x D tn lt t lt tnt j 1 ier 2 52 a EG 1 Mo x 0 Po x 32 CHAPTER 2 COMPUTATION OF STATIONARY STATES FOR THE GPE In the above equations we set x t lim x t Hence iterations in times correspond tott to iterations in the projected gradient When t tends towards infinity gives an approximation of the steady state which is solution to 2 51 that is supposed to exist here The ground state is again computed as a solution of the minimization problem of the energy functional E under the normalization constraint argmin E 2 53 Sllo 1 2 6 3 Time and space discretizations For the same reason as in the single component case we focus on schemes based on the Backwa
131. re terms that supplement the terms of this License by making excep tions from one or more of its conditions Additional permissions that are applicable to the entire Program shall be treated as though they were included in this License to the extent that they are valid under applicable law If additional permissions apply only to part of the Program that part may be used separately under those permissions but the entire Program remains governed by this License without regard to the additional permissions When you convey a copy of a covered work you may at your option remove any additional permissions from that copy or from any part of it Additional permissions may be written to require their own removal in certain cases when you modify the work You may place additional permissions on material added by you to a covered work for which you have or can give appropriate copyright permission Notwithstanding any other provision of this License for material you add to a covered work you may if authorized by the copyright holders of that material supplement the terms of this License with terms a Disclaiming warranty or limiting liability differently from the terms of sections 15 and 16 of this License or 115 b Requiring preservation of specified reasonable legal notices or author attributions in that material or in the Appropriate Legal Notices displayed by works containing it or c Prohibiting misrepresentation of the origin of th
132. re the trap frequencies in the directions x y and z respectively The quantity Uo defined by Ath as gas 1 3 m describes the interaction between atoms of the condensate as being the scattering length which is positive for a repulsive interaction and negative for an attractive interaction The operator L is such that Lz py YDx th xOy yOr 1 4 This is the z component of the angular momentum L x x P where the momentum operator P hV Pr Py pz The energy of the functional is defined by h 2 2 NO 14 Eb 5V4 Vly IY Qy L p dx 1 5 R3 2m 2 The wave function is generally normalized f wexlx 1 1 6 R3 1 3 WHICH PROBLEMS CAN GPELAB SOLVE 15 1 2 2 The dimensionless GPE Let us introduce the following changes of variables t i t gt Wm Min wr Wy wz Wm PR X gt X0Q0 ag A mw Y ne 1 7 yr 372 a Q gt Qum One obtains the dimensionless GPE OU ORB a spp oa de eo 2 where U N ra N a 0 Nas B hip ore 1 9 and L i x0y yO The potential is now 1 V x 0 ay 22 1 10 setting Yz y 2 Wa y 2 Wm The dimensionless energy functional Eg o is defined by Baot VOR Viv Slot otp dx 1 11 In the special case of a disk shaped condensation wg wy and wz gt wg This means that we have Ye 1 yy 1 and y gt 1 with w wy Then the 3d GPE reduces to a 2d GPE given by o 1 OP N V
133. re which everyone can redistribute and change under these terms To do so attach the following notices to the program It is safest to attach them to the start of each source file to most effectively state the exclusion of warranty and each file should have at least the copyright line and a pointer to where the full notice is found lt one line to give the program s name and a brief idea of what it does gt Copyright C lt year gt lt name of author gt This program is free software you can redistribute it and or modify it under the terms of the GNU General Public License as published by the Free Software Foundation either version 3 of the License or at your option any later version This program is distributed in the hope that it will be useful but WITHOUT ANY WAR RANTY without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PAR TICULAR PURPOSE See the GNU General Public License for more details You should have received a copy of the GNU General Public License along with this program If not see http www gnu org licenses Also add information on how to contact you by electronic and paper mail If the program does terminal interaction make it output a short notice like this when it starts in an interactive mode lt program gt Copyright C lt year gt lt name of author gt This program comes with ABSOLUTELY NO WARRANTY for details type show w This is free software and you are welcome t
134. recipients a copy of this License along with the Program You may charge any price or no price for each copy that you convey and you may offer support or warranty protection for a fee 5 Conveying Modified Source Versions You may convey a work based on the Program or the modifications to produce it from the Program in the form of source code under the terms of section 4 provided that you also meet all of these conditions a The work must carry prominent notices stating that you modified it and giving a relevant date b The work must carry prominent notices stating that it is released under this License and any conditions added under section 7 This requirement modifies the requirement in section 4 to keep intact all notices c You must license the entire work as a whole under this License to anyone who comes into pos session of a copy This License will therefore apply along with any applicable section 7 additional terms to the whole of the work and all its parts regardless of how they are packaged This License gives no permission to license the work in any other way but it does not invalidate such permission if you have separately received it d If the work has interactive user interfaces each must display Appropriate Legal Notices however if the Program has interactive interfaces that do not display Appropriate Legal Notices your work 113 need not make them do so A compilation of a covered work with other sep
135. rences scheme in space for the 1d and 2d cases but not for the 3d case We also discuss the possible classical choices of potentials that are met in the physics of Bose Einstein condensates but any other potential can be defined by the user itself and how to choose the initial data for evolving the imaginary time method In Chapter 3 a simple example of code is given for one complete problem It explains the general philosophy of GPELab and shows step by step the model code that a user has to define Next we give the full definition of all existing functions inside GPELab and describe them in details In particular the way the arguments must be defined is carefully explained In Chapter 4 we describe how different 1d 2d 3d stationary problems can be solved by GPELab This provides some generic codes for considering other problems that any user can meet In par ticular we give a few 1d 2d and 3d examples including different potentials nonlinearities coupled systems This offers the possibility of slightly modifying the code if your problem is not far from one of these examples In the case of a different problem we also explain what must be done in conjunction with Chapter 3 that gives more insight into the function use 13 14 CHAPTER 1 INTRODUCTION amp INFORMATIONS 1 2 The Gross Pitaevskii Equation GPE 1 2 1 The GPE equation coming from physics The aim of GPELab is to compute both stationary solutions and the dynamics of
136. rid points j Yk for j k Dyk setting Dik j k EN 31 lt 57 lt J landl lt k lt K 1 with J K gt 3 and uniform discretization steps hz and hy in the x and y directions respectively given by Ax Uj41 Tj hy 2La J Ayk Yk 1 Yk hy 2Ly E Since we assume that all the components are compactly supported in O then each 1 Ne satisfies a periodic boundary condition which can in fact be put to zero on OO and we can use discrete Fourier transforms For BESP the following approximation holds BEng _ HBEn APEn pBEn antx 2 55 2 54 IIS 2 6 EXTENSION TO THE MULTI COMPONENTS CASE 33 where Pl mero Ne j k 5 4 EDy x is the discrete unknown array in CAN and the right hand side is ben B At with B DT x k g kEDy r 9 OH k ED E CAN The operator ABE CMNe yw e CMN is defined by ABENG Arr P ARE ape Ex 4 vy pram a 2 56 age FUAN CNAN NCPN 8 The finite dimensional operator Ane is explicitly given through the matrices pi 0 i n Mall Ma e Vil geje o e gaf E Ded oo Mra 0 0 H Vna Var lt wen and Faa F E Fin penp EA en Fna 8 Fav gt FNN E In the above equations we set Fim Fim Eie i men q and 0 18 lrg w 1104118 9 6111 2 57 We also define on Kee Gi G G a AS et yyy Gra G n GK setting Ghnl Gin j k
137. ructure We use a centered gaussian as initial function see Table 4 5 InitialData_choice 2 Phi_O InitialData_Var2d Method Geometry2D Physics2D InitialData_choice Figure 4 5 Setting the initial data The outputs and the printing informations are then defined see Table 4 15 Outputs OutputsINI_Var2d Method Printing 1 Evo 15 Draw 1 Print Print_Var2d Printing Evo Draw Phi Outputs GPELab2d Phi_0 Method Geometry2D Physics2D Outputs Print Table 4 15 Setting the outputs and the printing of informations and then launching the simulation The results of the simulation are reported in Table and the solution is given in Figure 4 6 The evolution of the energy and the chemical potential are plotted in Figure 4 4 Ground state of a 2d Gross Pitaevskii equation with a quadratic plus quartic potential a cubic nonlinearity and a rotational op erator We consider the computation of the ground state of a 2d Gross Pitaevskii equation with a quadratic plus quartic potential a cubic nonlinearity and a rotational operator This is a typical example of 4 4 GROUND STATE OF 2D GPE EQUATION WITH A FAST ROTATION 71 Iteration 157 on 1000000 Dutputs of component 1 Square at the origin 0 06528267039451 x radius mean square 3 96548348047485 y radius mean square 3 96548346927846 Energy 51 22028604002639 Chemical potential 66 24901258432952 Angular moment
138. s Pitaevskii equation with quadratic potential cubic nonlinearity and dipole dipole interaction We now show how to compute the ground state of a Gross Pitaevskii equation with quadratic potential cubic nonlinearity and dipole dipole interaction in 3d We first build the Method and Geometry3D structures In Table we set GPELab to use the BESP scheme for a time step At equal to 107 with a spatial grid of 22 1 points in the z y and z direction on the domain 15 15 x 15 15 x 15 15 The stopping criterion is fixed to 107 Computation Ground Ncomponents 1 Type BESP Deltat 1le 2 Stop_time Stop_crit 1e 6 Method Method_Var3d Computation Ncomponents Type Deltat Stop_time Stop_crit xmin 10 xmax 10 ymin 10 ymax 10 zmin 10 zmax 10 Nx 276 1 Ny 276 1 Nz 276 1 Geometry3D Geometry3D_Var3d xmin xmax ymin ymax zmin zmax Nx Ny Nz Table 4 33 Setting the Method and Geometry3D We consider here the following Gross Pitaevskii equation with a dipole dipole interaction 1 1 iO V t x y 2 sAW t wy z 5 vela wlyl 1 1211 t y 2 1 3 cos a X ee ee HeT wy 2 W t a y 2 a f ts Paz w x y 2 re o y 2 X 9 with Ys Yy Ye 1 8 2000 and a 0 0 1 In GPELab the default nonlinearity is the cubic nonlinearity and the default potential is the quadratic potential We only have to
139. s consider any input variable xj in a set K of a function f In GPELab all inputs of f have already default values xj ult that can be modified For clarity in the sequel we designate this by the notation xj X xjdebult We essentially detail the Matlab functions for the two dimensional case Unless precise the extension from the 2d to the 1d and 3d cases is done by changing the functions names For example the Method_Var2d function corresponds to the Method_Varid function in 1d and to the Method_Var3d function in 3d If changing the dimension implies any modification of the number or 3 5 SETTING THE NUMERICAL SCHEME AND THE GEOMETRY 41 the nature of the input or output arguments of the function then we will precise it Concerning the form of the variables x y and z we use the standard meshgrid ordering of variables More precisely this means that x M1 y K in the 1d case z y My n K in the 2d case and zx y z are in Mn NN K for the 3d case Here K R The same situation occurs when computing the set of frequencies for example to compute nonlocal nonlinear interactions like for dipolar gazes but K C 3 5 Setting the numerical scheme and the geometry First the user has to define the geometry and the numerical method There exists two variables that need to be defined Method and Geometry Those are created by using the two following functions Method_Var2d and Geometry2D_Var2d 3 5 1 The Method_Var2d func
140. s for the same material under section 10 116 APPENDIX B LICENSE 9 Acceptance Not Required for Having Copies You are not required to accept this License in order to receive or run a copy of the Program Ancillary propagation of a covered work occurring solely as a consequence of using peer to peer transmission to receive a copy likewise does not require acceptance However nothing other than this License grants you permission to propagate or modify any covered work These actions infringe copyright if you do not accept this License Therefore by modifying or propagating a covered work you indicate your acceptance of this License to do so 10 Automatic Licensing of Downstream Recipients Each time you convey a covered work the recipient automatically receives a license from the original licensors to run modify and propagate that work subject to this License You are not responsible for enforcing compliance by third parties with this License An entity transaction is a transaction transferring control of an organization or substantially all assets of one or subdividing an organization or merging organizations If propagation of a covered work results from an entity transaction each party to that transaction who receives a copy of the work also receives whatever licenses to the work the party s predecessor in interest had or could give under the previous paragraph plus a right to possession of the Corresponding Source of the
141. s informations needed to display the wave functions in 3d The square modulus of the wave function is drawn in 3d by using an isovalues surface and the phase of the wave function which is displayed by using some slices along the x y and z directions The function has the following optional arguments e view Nor M 3 R 3 is a variable corresponding to the viewing angle For more informations see the view Matlab function e iso R 0 001 is the isovalue chosen for drawing the modulus of the wave functions e alpha R 0 6 is the transparency of the isovalue surface 60 CHAPTER 3 HOW TO USE GPELAB STATIONARY SOLUTIONS e aspect Mj 3 R 1 1 1 is a variable corresponding to the data aspect ratio For more informations see the daspect Matlab function e Sx R 0 respectively Sy Sz is a coordinate on the z axis where an orthogonal slice to the x axis respectively y axis z axis will be drawn e map S jet is a variable corresponding to the colormap of the figures It must be either jet hsv hot cool spring summer autumn winter gray bone gt copper pink or lines For example suppose that one wants to compute a rotating Bose Einstein condensate where the rotation is along the z direction and then to see the condensate along the z axis to get a view from above Thus the view angle as set by Matlab should be 0 0 1 Moreover i
142. software interchange for a price no more than your reasonable cost of physically performing this conveying of source or 2 access to copy the Corresponding Source from a network server at no charge c Convey individual copies of the object code with a copy of the written offer to provide the Corresponding Source This alternative is allowed only occasionally and noncommercially and only if you received the object code with such an offer in accord with subsection 6b d Convey the object code by offering access from a designated place gratis or for a charge and offer equivalent access to the Corresponding Source in the same way through the same place at no further charge You need not require recipients to copy the Corresponding Source along with the object code If the place to copy the object code is a network server the Corresponding Source may be on a different server operated by you or a third party that supports equivalent copying facilities provided you maintain clear directions next to the object code saying where to find the Corresponding Source Regardless of what server hosts the Corresponding Source you remain obligated to ensure that it is available for as long as needed to satisfy these requirements e Convey the object code using peer to peer transmission provided you inform other peers where the object code and Corresponding Source of the work are being offered to the general public at no charge under subsection 6d A separ
143. t would be convenient to set a slightly larger isovalue than the default one for the surface of the modulus of the wave function say 0 01 to better distinguish the vortices Then the Figure structure should be defined the same way as in table view 0 0 1 iso 0 01 Figure Figure_Var3d view iso Table 3 39 An example to use the Figure_Var3d function 3 7 6 The GPELab2d function Phi Outputs GPELab2d Phi_0 Method Geometry2D Physics2D 0utputs Continuation Print Figure Table 3 40 The GPELab2d function The GPELab2d function is the main function to launch a full simulation with respect to the given configuration The output arguments are e Phi Ci y Mn n C the wave functions computed at the final time with respect to a stopping criterion for a stationary state or a fixed time for a dynamical computation e and Outputs S which is a structure that contains all the outputs computed during the simulation The initial data Phi_0 C1 v Mvw n C and the Method Geometry2D Physics2D and Outputs structures are required arguments The optional arguments that can be considered are the following e Continuation S is the continuation structure if one wants to use a continuation method e Print S Print_Var2d is the printing structure e Figure S Figure_Var2d is the structure that fixes the parameters for drawing the figures We report in Table a model example of use of the GP
144. th the suffix 1d instead of 2d BESP BEFD CNSP and CNFD methods are coded From the user point of view the example file see subsection 4 1 page 63 and subsection 1 3 page 68 shows that considering a one or a two dimensional problem does not need a lot of modifications in the main GPELab file that is launched for the simulations 30 CHAPTER 2 COMPUTATION OF STATIONARY STATES FOR THE GPE 2 5 2 Three dimensional case For the three dimensional case only BESP and CNSP methods are coded In the same spirit as for the one and two dimensional functions the suffix is 3d An example is given in subsection 4 7 page 79 2 6 Extension to the multi components case 2 6 1 The multi components GPE The continuous normalized gradient flow can also be extended to the multi components case i e a system of coupled Gross Pitaevskii equations For the sake of conciseness the spatial variable x is defined by x 21 24 R We denote by Y 11 Yp with Ne N N 0 a vector of Ne wave functions and consider the following generic system of Gross Pitaevskii equations d 0 W t x AE x V x U t x Y GI x Ip W t x j 1 6F W t x x U t x ERT x R 2 49 with initial condition Y t 0 x Vo x and where the operators are defined by e the diagonal Laplacian AW t x Av t X j 1 Nes e the potential matrix Vii x V x Vin x ma oe 7 a A Vn i x Vn 2 x Vin n x
145. the present version but may differ in detail to address new problems or concerns 118 APPENDIX B LICENSE Each version is given a distinguishing version number If the Program specifies that a certain numbered version of the GNU General Public License or any later version applies to it you have the option of following the terms and conditions either of that numbered version or of any later version published by the Free Software Foundation If the Program does not specify a version number of the GNU General Public License you may choose any version ever published by the Free Software Foundation If the Program specifies that a proxy can decide which future versions of the GNU General Public License can be used that proxy s public statement of acceptance of a version permanently authorizes you to choose that version for the Program Later license versions may give you additional or different permissions However no additional obligations are imposed on any author or copyright holder as a result of your choosing to follow a later version 15 Disclaimer of Warranty THERE IS NO WARRANTY FOR THE PROGRAM TO THE EXTENT PERMITTED BY APPLICABLE LAW EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPY RIGHT HOLDERS AND OR OTHER PARTIES PROVIDE THE PROGRAM OAS ISO WITH OUT WARRANTY OF ANY KIND EITHER EXPRESSED OR IMPLIED INCLUDING BUT NOT LIMITED TO THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE THE ENTIRE
146. thermore the algorithms use recent developments in numerical methods GPELab can also be used on parallel platforms so that it uses the power of parallel architecture The time of computation of one given problem depends strongly of your physical configuration for example strong nonlinear interactions or fast rotating gazes and the computer resources that you can access 1 4 How to read this manual One way to read the manual is the following First Chapter 1 provides the general notations the equations and the physical quantities that are involved into the equations Then you can read Sections and to understand what we wish to compute in the case of stationary solutions Next you can go to Chapter 4 to see examples of computations for different physical situations and the associated GPELab codes If you are interested in the numerical methods that are coded for the stationary states computation you can read Sections 2 3 and 2 4 for the two dimensional case and next Section 2 5 for the 1d and 3d cases Furthermore if you want to understand the numerical schemes that are used as well as the notations for a multi component case then you have to take a look at Section Finally Section 3 explains in details the different functions that are included in GPELab Reading again the examples of Section 4 will give you more understanding into the codes Another way of reading the user manual is the standard linear one chapter after chapter
147. thod Physics2D Physics2D Gradientx_Var2d Method Physics2D Physics2D Gradienty_Var2d Method Physics2D Table 4 19 Building and defining the Physics2D structure We choose the Thomas Fermi approximation as initial data see Table 4 20 InitialData_choice 2 Phi_O InitialData_Var2d Method Geometry2D Physics2D InitialData_choice Table 4 20 Building the initial data The outputs and printing options are defined in Table Outputs OutputsINI_Var2d Method Printing 1 Evo 15 Draw 1 Print Print_Var2d Printing Evo Draw Phi Outputs GPELab2d Phi_0 Method Geometry2D Physics2D Outputs Print Table 4 21 Defining the Outputs and Print structures and then launching the computation using GPELab2d We obtain in Table the informations on the ground state at the end of the computation Iteration 46766 on 1000000 Qutputs of component 1 Square at the origin 0 00000000000000 x radius mean square 4 57951169686043 y radius mean square 4 57951071463754 Energy 115 52164061561449 Chemical potential 122 58168418655728 Angular momentum 146 32747911959200 Energy evolution 0 00000000141087 gt Table 4 22 Printed outputs for the 2d GPE with cubic nonlinearity rotation and quadratic plus quartic potential 74 CHAPTER 4 EXAMPLES OF SIMULATIONS FOR STATIONARY SOLUTIONS We report on figure 4 8 the modulus of the stationary
148. tial patent claims to make use sell offer for sale import and otherwise run modify and propagate the contents of its contributor version In the following three paragraphs a patent license is any express agreement or commitment however denominated not to enforce a patent such as an express permission to practice a patent or covenant not to sue for patent infringement To grant such a patent license to a party means to make such an agreement or commitment not to enforce a patent against the party If you convey a covered work knowingly relying on a patent license and the Corresponding Source of the work is not available for anyone to copy free of charge and under the terms of this 117 License through a publicly available network server or other readily accessible means then you must either 1 cause the Corresponding Source to be so available or 2 arrange to deprive yourself of the benefit of the patent license for this particular work or 3 arrange in a manner consistent with the requirements of this License to extend the patent license to downstream recipients Knowingly relying means you have actual knowledge that but for the patent license your conveying the covered work in a country or your recipient s use of the covered work in a country would infringe one or more identifiable patents in that country that you have reason to believe are valid If pursuant to or in connection with a single transaction or
149. tion Method Method_Var2d Computation Ncomponents Type Deltat Stop_time Stop_crit Max_iter Precond_type Output Splitting BESP Solver_FD Iterative_tol Iterative_maxit Table 3 12 The Method_Var2d function The Method_Var2d function creates the Method structure that contains all the parameters rel ative to the method By method we mean the solver which is used to compute a solution This includes the kind of computation dynamics or ground state the number of components the type of scheme BESP BEFD CNSP CNFD for the ground state and Relaxation Splitting for the dynamics the semi discretization parameters and other inputs that we explain below The only output is the structure Method As seen above the input variable of Method_Var2d have already default values that may be modified The optional arguments are the following e Computation S Ground is a variable that must be Ground to compute a ground state by using the Continuous Normalized Gradient Flow imaginary time method e Ncomponents N 1 is a variable corresponding to the number of components that describe the condensate e Type S BESP is a variable corresponding to the scheme used in the computation In the case of a ground state computation it must be either BEFD to use the Backward Euler Finite Difference scheme see section 2 3 2 gt CNFD to use the Crank Nicolson Finite Difference scheme see section P 4 gt
150. tion 3 6 Setting the physical problem We now explain how to set the physical problem We consider the following general GPE with Ne components each one being defined in the d dimensional space by d 0 V t x AV t x V x t x Y G x dn V t x j 1 6F WU t x x U t x Gs E RF x R 3 1 44 CHAPTER 3 HOW TO USE GPELAB STATIONARY SOLUTIONS This system corresponds to the one developed in Section 2 6 page Here 6 8 are two real valued constants in R The energy for each component is set as E U ft Vy t x 2ax d R v t x v S G x dn Fenergy U t x x vex dx Rd k 1 j 3 2 for each j 1 Ne and the chemical potential for each component is set as u af Vy t x dx 3 3 d 2 w t x veo Y G x de F U t x x k 1 v x ax j for each j 1 Ne 3 6 1 The Physics2D_Var2d function Physics2D Physics2D_Var2d Method Delta Beta Omega Table 3 16 The Physics2D_Var2d function The Physics2D_Var2d function builds the Physics2D structure and enables the user to define the basic physical constants 6 and the rotation speed Q if the gradients operators are set as default see section 8 6 5 The Physics2D structure also contains the physical operators as explained below The Method structure is a required argument and the optional arguments are the following e Delta R 1 2 is a variable corresponding to the co
151. tion Nonlinearity_energy in is given the nonlinear energy operator is defined as follows for each j k 1 N Nonlinearity_energy W t x z y if j k Fenergy jk Ute 0 if j E k If Nonlinearity_energy is a cell array of function in Ne Ne F Cro N A Mn 04 C Ma Ne R My na C then the nonlinear energy operator is Fenergy j k U t x x y Nonlinearity_energy j k W t x 2 y for j k 1 N The default argument is Cubic_energy2d which corresponds to 1 DA Fenergy j E b xX x y zW t x if j k Dif j 4k This way of proceeding allows us to continue the example from Section page where we set a potential in the case of an internal atomic Josephson junction We also need to define the coupled nonlinearities if we want to effectively take into account all the effects in the system of equations 17 In the case of a two components Gross Pitaevskii equation with a Josephson junction we have Foo W t x x bl Bilal and F o W t x x Fa t x x 0 This is done in table 3 21 where we create a cell array of functions corresponding to the previous nonlinearities and then define the nonlinear operator by using the Nonlinearity_Var2d function Fy W t x x buly Br2 1a1 48 CHAPTER 3 HOW TO USE GPELAB STATIONARY SOLUTIONS function NL Example_nonlinearity Beta_11 Beta_22 Beta_12 NL cel1 2 NL 1 1 Phi x y Beta_11 abs Phit1 72 Beta_12 abs Ph
152. tions but does not want to slow the program by drawing the wave function s modulus and angle then one can define the Print structure by using the Print_Var2d function like in table 3 7 LAUNCHING THE SIMULATION 59 Printing 1 Evo 10 Draw 0 Print Print_Var2d Printing Evo Draw Table 3 35 An example of use for the Print_Var2d function 3 7 4 The Figure_Var2d function Figure Figure_Var2d map Table 3 36 The Figure_Var2d function The Figure_Var2d function builds the Figure structure which contains informations needed to draw figures in 2d We have the following optional argument e map S jet is a variable corresponding to the colormap of the figures It must be either jet hsv hot cool spring summer autumn winter gray bone copper pink or lines see the Matlab documentation for further informations about colormay If one wants to draw figures using the hot colormap for example then it can be done by defining the Figure structure as in table map hot Figure Figure_Var2d map Table 3 37 An example of how to use the Figure_Var2d function 3 7 5 The Figure_Var3d function Figure Figure_Var3d view iso alpha aspect Sx Sy Sz map Table 3 38 The Figure_Var3d function The Figure_Var3d function builds the Figure structure which contain
153. to compute a ground state problem However some changes have to be made In this section we ll state the functions that need to be used and how they differ from the one used in the stationary case First we have to define the Method and Geometry variables This is done by using the Method_Vard2d and Geometry2D_Var2d functions The Geometry2D_Var2d is handled in the same way as in the stationary state we refer to section page on how to use it The Method_Var2d function is the same as the one in the stationary state see section page 41 however we recall the arguments that are important when computing a dynamical problem 5 5 1 The Method_Var2d function Method Method_Var2d Computation Ncomponents Type Deltat Stop_time Stop_crit Max_iter Precond_type Output Splitting BESP Solver_FD Iterative_tol Iterative_maxit Table 5 1 The Method_Var2d function The Method_Var2d function creates the Method structure that contains all the parameters rela tive to the method We refer to section page for further details The optional arguments are the following e Computation S Ground is a variable that must be Dynamic to compute the dynamic of the GPE e Ncomponents N 1 is a variable corresponding to the number of components that describe the condensate e Type S BESP is a variable corresponding to the scheme used in the computation In the case of a dynamical computation it must either be Splitt
154. tol R 1e 9 is a variable corresponding to the stopping criterion related to the difference between two successive iterates in the Krylov solver Iterative_maxit N 1e3 is a variable corresponding to the stopping criterion related to the maximum number of iterations in the Krylov solver For example we want to compute a stationary solution for a single component BEC by using the BESP scheme We choose a time step At 107 and a stopping criterion for e 1078 We set the maximal number of iterations to 10 and we choose to compute outputs during the simulation This gives the code in table Computation Ground Ncomponents 1 Type BESP Deltat ie 2 Stop_time Stop_crit le 8 Max_iter 10e6 Precond_type None Output 1 Method Method_Var2d Computation Ncomponents Type Deltat Stop_time Stop_crit Max_iter Precond_type Output Table 3 13 An example of initialization and use of the Method_Var2d function 3 6 SETTING THE PHYSICAL PROBLEM 43 3 5 2 The Geometry2D_Var2d m function The call to the function Geometry2D_Var2d m is Geometry2D Geometry2D_Var2d xmin xmax ymin ymax Nx Ny Table 3 14 The Geometry2D_Var2d function The aim of the Geometry2D_Var2d m function is to create the Geometry2D structure which contains the size of the computational box and the number of points in each spatial direction including the boundaries Note that t
155. tputsINI_Var2d Method Printing 1 Evo 15 Draw 1 Print Print_Var2d Printing Evo Draw Table 3 4 Printing drawing informations during the computations We are now ready to launch the simulation This is done by using the GPELab2d function which gather all the previous structure and thus the informations about the simulation The command is given in Table 3 5 Phi Outputs GPELab2d Phi_0 Method Geometry2D Physics2D Outputs Print Table 3 5 Launching the computation of the solution 38 CHAPTER 3 HOW TO USE GPELAB STATIONARY SOLUTIONS At the end of the simulation we obtain the following figures angle phi x y of component 1 a 2 1 o 4 2 lt o 2 6 4 2 o 2 4 6 8 10 x a Modulus of the ground state b Phase of the ground state Iphi x y of component 1 10 10 8 6 4 2 0 025 8 6 0 02 4 2 0 015 gt 0 2 0 01 4 4 0 005 s 8 o 10 o 8 6 4 2 o 2 4 6 3 10 1 x gt o 2 6 8 19 4 Figure 3 1 Ground state computed with GPELab using the parameters from Section 3 2 As we have seen with this example using GPELab is quite easy and direct We report more complex examples in Chapter 4 for different physical situations involving ground states dynamics 1d 2d 3d problems multi components gazes 3 3 Variables types and various notions required for Matlab To help the user who is not familiar with Matlab we give quickl
156. ucts If such problems arise substantially in other domains we stand ready to extend this provision to those domains in future versions of the GPL as needed to protect the freedom of users Finally every program is threatened constantly by software patents States should not allow patents to restrict development and use of software on general purpose computers but in those that do we wish to avoid the special danger that patents applied to a free program could make it effectively proprietary To prevent this the GPL assures that patents cannot be used to render the program non free The precise terms and conditions for copying distribution and modification follow TERMS AND CONDITIONS O Definitions This Licens refers to version 3 of the GNU General Public License opyright also means copyright like laws that apply to other kinds of works such as semicon 2C ight al ight like 1 that ly to other kinds of k h i ductor masks The Program refers to any copyrightable work licensed under this License Each licensee is addressed as you Licensees and recipients may be individuals or organizations To modify a work means to copy from or adapt all or part of the work in a fashion requiring copyright permission other than the making of an exact copy The resulting work is called a modified version of the earlier work or a work based on the earlier work A covered work means either the unmodi
157. um 0 00000000000000 Energy evolution 0 00000000003103 CPU time 172 47 gt Table 4 16 Printed outputs for the 2d GPE with cubic nonlinearity and optical lattice My ot component 1 20 0 08 15 i 0 05 sL ol 003 sL 0 02 10 oor 18 1 L 1 fi 1 fi 1 j as a0 o g 10 zo 20 ho y Figure 4 6 Modulus of the ground state 110 Energy Chemical potential Figure 4 7 Evolution of the energy and the chemical potential during the computation computation of a ground state where a fast rotation can be considered We begin by constructing the Method and Geometry2D structures We use the BESP scheme with a time step equal to 1073 to solve the problem on a spatial grid involving 2 1 points in the z and y directions for the 72 CHAPTER 4 EXAMPLES OF SIMULATIONS FOR STATIONARY SOLUTIONS domain 10 10 x 10 10 The stopping criterion is fixed to 1075 see Table 4 17 Computation Ground Ncomponents 1 Type BESP Deltat 1e 3 Stop_time Stop_crit le 5 Method Method_Var2d Computation Ncomponents Type Deltat Stop_time Stop_crit xmin 10 xmax 10 ymin 10 ymax 10 Nx 278 1 Ny 278 1 Geometry2D Geometry2D_Var2d xmin xmax ymin ymax Nx Ny Table 4 17 Setting the Method and Geometry2D structures In this section the Gross Pitaevskii equation that we consider is 1 l a K 2
158. utions to the minimization problem 2 9 then an iterative procedure is of course needed This means that an initial guess has to be given to the method to initialize it and then the minimization process compute or try to compute a minimal solution through iterations The aim of this section is to give informations concerning the choice of this initial guess which is often built analytically as an approximate minimal state for simplified problems For a non rotating BEC it can be proved that the global minimal solution is unique and gives a ground state g gt 0 for a positive initial data fp Therefore one usually choose the solution to the linear Schr dinger equation with harmonic potential when we are under the critical frequency Q lt Yay with yey min yz yy for a harmonic trap V x a yy 2 11 2 2 APPROXIMATE SOLUTIONS AS INITIAL GUESS AND POTENTIALS 21 The initial data is then given by x 1 A e7 It yyy 2 2 12 T This choice can also be considered for the non rotating harmonic potential and a potential of a stirrer corresponding to a far blue detuned Gaussian laser beam toroidal trap 1 Vita e sy tae he 2 13 When a rotation is taken into account the choice of the initial data is less clear In 15 Bao et al propose to choose for yz yy 1 1 2 Gnolx QGho X II Q bno x 2467 Ilo o x 2 14 with ee eee 2 Yat iYyy in z2 2 eet yy 2 v
159. vol 184 pp 319 332 2003 3 J Garnier F Kh Abdullaev and B B Baizakov Collapse of a Bose Einstein condensate induced by fluctuations of laser intensity Phy Rev A vol 69 053607 pp 369 386 2004 4 R Marty On a splitting scheme for the nonlinear Schr dinger equation in a random media Commun Math Sci 4 679 705 2006 5 Engel Klaus Jochen Nagel Rainer One parameter semigroups for linear evolution equations Springer 2000 6 C Besse B Bid garay and S Descombes Order estimates in time of splitting methods for the nonlinear schr dinger equation SIAM J Numer Anal Vol 40 No 1 pp 26 40 2002 7 C Besse A relaxation scheme for the non linear schr dinger equation SIAM J Numer Anal Vol 42 No 3 pp 934 952 2004 8 G Da Prato and J Zabczyk Stochastic equations in infinite dimensions Encyclopedia of math ematics and its applications Cambridge University Press Cambridge England 1992 9 A De Bouard and R Fukuizumi Stochastic fluctuations in the Gross Pitaevskii equation Non linearity 20 pp 2821 2844 2007 10 A De Bouard and A Debussche The stochastic nonlinear Schr dinger equation in H1 Stochas tic Analysis and Appl vol 21 p 97 126 2003 11 A De Bouard and R Fukuizumi Representation formula for stochastic Schrodinger evolution equations and applications preprint 12 W Bao and H Wang An efficient and spectrally accurate numerical method f
160. y a few basic notions For more informations we refer for example to the Matlab user guide and help These examples are just to fix what is a matrix a cell a function or a structure and do not intend to give more informations In practice GPELab uses the following variables types e Matrix type Matrices are the basic variables of Matlab A matrix is a multi dimensional array of integers real or complex numbers To access a certain element of a matrix the indexation must be done using parentheses y wre Pe Pe 1 3 2 3 6 4 3 w 2 4 3 O Table 3 6 An example of 2 x 2 real valued matrix Basic operations such as additions or multiplications for matrices are implemented in Matlab Moreover the element wise operations are done by adding a dot before the symbol of the operation To initialize a matrix one can choose the zeros function that will initialize a matrix of zeros e Cell type Cells are variables that gather elements of Matlab such as matrices functions or strings http www mathworks fr fr help matlab 3 4 NOTATIONS AND PRELIMINARY REMARKS 39 A 1 2 3 4 B 1 2 2 1 A B ans 14 6 4 gt Table 3 7 An example of element wise multiplication for matrices A 1 string 3 1 1 2 1 1 string 3 2x2 double Table 3 8 An example of cell They use the same indexation as matrices and to access a certain element of a matrix th
161. y modifying the Physics2D structure Note that in GPELab the solution of the system is defined as a cell array of matrices Cy v Mwny n2 C This function must be provided with the Method and Physics2D structures and it has the following optional arguments e Nonlinearity If a function Nonlinearity in F Cn Ne MNy No C MN No Ry MN No C is given the physical nonlinearity will be defined as follows for each j k 1 N Nonlinearity V t x 1 y if j k Peal Savas y W t x 1 y if j 3 6 SETTING THE PHYSICAL PROBLEM 47 If Nonlinearity is a cell array of functions in Cn No F Cno NAM na CJ Muy Ne R M no 1 C then the nonlinear operator will be defined by Fav x 2 y Nonlinearity j k P t x x y for j k 1 N The default argument is Cubic2d which corresponds to Dir gt e G My n C ones N_c is a complex valued variable that multiplies the nonlinearity element by element leading to the following nonlinearity definition Es W t x 2 y GU k Nonlinearity j k P t x x y for j k 1 N e Nonlinearity_energy is a nonlinear operator used to compute the energy associated to the physical nonlinearity It corresponds to Fenergy Y t x x in the energy definition 3 2 Note that it must be the same type of variable as the variable Nonlinearity If the variable G is defined it will also be multiplied element by element by Nonlinearity_energy If a func

GPELab A Matlab toolbox for computing stationary solutions

Contents

Download Pdf Manuals

Related Search

Related Contents