Modelling spin-1 Rubidium-87 Bose Einstein Condensate to study ground states under inhomogenous magnetic fields

IN

DEGREE PROJECT ENGINEERING PHYSICS, SECOND CYCLE, 30 CREDITS

STOCKHOLM SWEDEN 2019,

Modelling spin-1 Rubidium-87 Bose Einstein Condensate to study ground states under

inhomogenous magnetic fields

SANDHYA GANESH

KTH ROYAL INSTITUTE OF TECHNOLOGY

Modelling spin-1 ⁸⁷ Rb Bose Einstein Condensate to study ground states under

inhomogenous magnetic fields

Sandhya Ganesh

Masters Thesis work submitted in partial fulfilment of the requirements for the award of the degree of

Master of Science in Electrophysics (TRITA-SCI-GRU 2019:401)

This work was carried out at the

Institute of Photonic Sciences - ICFO, Castelldefells - Spain Supervisor: Prof. Dr. Morgan. W. Mitchell,

Atomic Quantum Optics Group - ICFO Examiner: Prof. Dr. Valdas Pasiskevicius,

Department of Applied Physics - KTH

Abstract

Precision studies of ⁸⁷Rb spinor Bose Einstein Condensate quantum dynamics with non-demolition mea- surements provide opportunities for understanding fundamental physics of low temperature systems and explore its applications in high resolution magnetic field sensing. A cavity enhanced detection of the mag- netization in the cold atomic ensembles provide pathways for magnetometry with the spinor condensate, with the ensemble magnetization being the preserved observable for the quantum non-demolition measure- ments. This cavity enhancement results in the ’pancake’-shaped condensates in an optical lattice. Studying the ground states of such pancake condensates serve as the starting point for understanding variety of phenomena such as low-energy dynamics, spin textures, effects of magnetic dipole interactions. The aim of this thesis is to develop a computational framework using numerical methods to simulate the ground state solutions of a F=1 hyperfine manifold of⁸⁷Rb Bose Einstein Condensate by solving multi-component Gross-Pitaevskii equation. We study the mean field ground states under external bias fields and in the presence of the field generated due to ensemble magnetic moments. This thesis project was carried out in the Atomic Quantum Optics group at ICFO headed by Prof. Morgan. W. Mitchell.

Acknowledgements

I would like to take this opportunity to thank everyone who helped me during the course of this thesis and Masters. Pursuing higher education in Physics was a personal dream for over a decade, and I am grateful to everyone who taught me at all steps during this journey, and helped me learn better everytime.

Firstly, I would like to thank Prof. Morgan Mitchell, for providing me the opportunity to work with his team and offering me an exciting problem to work on. His scientific inputs and guidance are an invalu- able part of my learning curve in Physics. I thank Prof Valdas Pasiskevicius for offering me his kind support and confidence throughout the course of this thesis. I would like to thank Pau Gomez and Daniel Benedicto Orenes for offering me valuable guidance and perspectives during the project in terms of experimental and computational aspects. I would especially like to thank Chiara Mazzinghi for her immense support at all steps of this thesis, making me a better student in Physics. Discussions with her helped me develop a deeper understanding of this field, and develop better problem solving skills in terms of computations and experimental techniques. Learning experiences with her in the lab with the ’spinor BEC machine’, discussions during after-lunch coffee times and all the moments with white-board and marker, are things that I would always cherish. I would like to thank the entire AQO team for offering me a great time in the ’007’ office, during lunch, and during all the team events over the summer. Special thanks to all my ICFOnian friends in Castelldefells, who made me feel at home in Barcelona.

Finally I would like to thank my parents, relatives and friends in India for constantly supporting me with their love and care and letting me have the privilege of pursuing studies in fundamental science.

Contents

1 Introduction 2

2 Spinor Bose-Einstein Condensates 4

2.1 Spin-1 Operators . . . 4

2.2 Single-particle Hamiltonian . . . 4

2.3 Interaction Hamiltonian . . . 5

2.4 Zeeman Interaction . . . 6

2.5 Mean-field Description . . . 6

2.6 Gross-Pitaevskii Equation . . . 7

2.6.1 Single-Mode Approximation . . . 8

3 Numerical methods for simulating spin-1 BECs 9 3.1 Dimensionless GPE . . . 9

3.2 Stationary States . . . 10

3.3 Conjugate Normalised Gradient Flow . . . 10

3.3.1 Time and Space discretizations . . . 10

3.3.2 Initial Condition for CNGF and Iterative Solvers . . . 12

3.4 Extension to the three component case . . . 12

4 Ground State Computation for Static Magnetic Fields 14 4.1 Numerical ground state solutions at B=0 . . . 14

4.2 Numerical ground state solutions at B6=0 . . . 16

4.2.1 With non-zero linear and quadratic Zeeman terms . . . 16

4.2.2 With linear Zeeman term p = 0 and non-zero quadratic Zeeman term . . . 17

4.3 On the choice of computational parameters . . . 18

5 Ground State Computation for Inhomogenous Magnetic Fields 20 5.1 Magnetic field generated by a neighbouring pancake . . . 21

5.1.1 Magnetic field distribution . . . 21

5.2 Ground state solutions . . . 22

5.2.1 Ground State Solutions with magnetic field generated by one pancake . . . 22

5.2.2 Ground State Solutions with magnetic field generated by two pancakes . . . 24

5.3 Numerical estimate of inter-pancake dipole forces . . . 24

6 Conclusion and Outlook 26

Chapter 1

Introduction

Since the first experimental realization of Bose-Einstein condensation [1], ultra cold atoms have established themselves as a versatile platform to simulate the dynamics of many-body quantum systems, owing to their rich experimental control and isolation from the environment [2]. Further development in the technologies that followed, interests in understanding dilute and weakly interacting systems gained impetus, pushing the field to further improvement in techniques. This lead to the development of the optical dipole trap and all-optical evaporation strategy [3, 4]. The magnetic traps earlier employed in obtaining the condensate, pose limitations in studying a multi-component condensate as only weak-field seeking atomic states are confined; whereas in an optical trap based on the optical dipole force confines atoms in all hyperfine states.

This provided the the opportunity to explore a system of degenerate Bose gases in the presence of the new additional degree of freedom in the spin space, which is thus known as a spinor BEC. A spinor BEC is interesting and useful to explore for various reasons. Theoretical studies on them promise to shed insights on a range of topics such as the role of symmetry and topology in quantum-ordered materials, quantum phase transitions, nonequilibrium quantum dynamics, entanglement, spin squeezing, spin-wave formation and suppression of quantum phase diffusion [5, 6].

In the field of quantum metrology applied to atomic sensors, a spinor BEC system offers a promising regime to detect magnetic fields with high spatial and temporal resolution. Several features of spinor Bose- Einstein condensates qualify them for use in such spatially resolving sensors. The Bose-Einstein condensate is a high-density, low-energy medium. The high density reduces the atomic shot noise[7] and low energy permits the use of shallow optical traps, which do not affect the atomic spin and from which spontaneous scattering is minimal. The rotational symmetry of the s-wave interactions is preserved, ensuring that Lar- mor precession in a fully magnetized gas under a linear Zeeman shift is unaffected by interactions [6].

Dispersive optical measurements, such as a Faraday rotation measurements on the cold atomic ensembles are used to convey information about their density as well as the spin polarization. Performing them non- destructively, serve as a means to spin-squeezing [8, 6], and this is particularly facilitated by the sensitivity of the optical cavities [9].

Motivated by these possibilities, an interesting probe to explore a spinor BEC system further, would be to employ a cavity aided detection and manipulation [9]. Such a cavity enhancement in the system, with a probe resonated in a linear optical cavity results in ‘pancake’ shaped condensates in the optical lattice, as illustrated in Fig. 1.1.

Figure 1.1: Illustration of ’pancake’-shaped condensates formed in optical lattice: a cavity designed to support two different wavelengths - one as a probe at 780 nm for detection and another for creating an optical lattice at 1560 nm. The pancakes are spaced at 0.78 µm as shown.

This work is mainly focused on developing a computational framework using numerical methods to model the ground state spin dynamics of such a⁸⁷Rb spin-1 BEC. With pancake shaped atomic ensembles in perspective, we will model the condensate in two dimensions. Using suitable numerical conditions, we will characterise the ground state solutions obtained at different external magnetic fields and test our model as we solve for mean-field ground state solutions, and look at the validity of single mode approximation.

The thesis is organized as follows:

• Chapter 2 contains the relevant theoretical framework. We will start with a review of a spin-1 system and develop the interaction Hamiltonian. We will review the fundamentals of mean-field description and derive the Gross-Pitaevskii equation(GPE) which is employed for numerical modelling.

• Chapter 3 contains the fundamentals of the numerical methods employed. We will discuss the suit- able time and space discretisations and iterative algorithms for solving GPE to obtain ground state solutions.

• In Chapter 4 we start testing for the validity of our model, and study the system under different applied external homogenous magnetic fields and observe the effects of Zeeman dynamics in the spin ensembles.

• In Chapter 5 we will study the ground states under the magnetization effect of the atomic ensembles.

We will look at the ground state solutions of the pancakes under the effect of a magnetic field produced by the magnetic moment in the neighbouring atomic ensembles.

Chapter 2

Spinor Bose-Einstein Condensates

In this chapter, we will review the fundamentals of the theoretical framework of a spin-1 Bose Einstein Condensate. We will start with reviewing relevant spin-1 operators and discuss the features of a spin-1 interacting system employing the mean-field description. We will then derive the Gross-Pitaevskii equation, which is solved numerically in the later sections to obtain the ground state solutions.

2.1 Spin-1 Operators

A spin-¹₂ system is described by a spin wave function given by the complex two component vector Ψ_1/2= (ψ_↑, ψ_↓)^T [10]. The set of possible physical transformations acting on these states are contained in the SU(2) group comprising of 2x2 complex unitary matrices - the Pauli matrices σ_j where j ∈ {x, y, z}. The spin operators given by these matrices are related by the commutation relation, [ ˆf_x, ˆf_y] = i ˆf_z(where ~ = 1).

When we get to the spin-1 system, as we have an additional degree of freedom, the possible physical transformations are described by the SU(3) group with 3x3 complex unitary matrices. The generators of the SU(3) group that span the space of the one-body observables of a spin-1 particle consists of nine operators - an identity operator, 3 spin vector operators( ˆf ) and 5 spin quadruple operators(ˆq). In the basis of ˆfz eigenstates, written as[11],

fˆx= 1

√2





0 1 0 1 0 1 0 1 0



 fˆy= i

√2





0 −1 0

1 0 −1

0 1 0



 fˆz=





1 0 0

0 0 0

0 0 −1



 (2.1)

ˆ qyz= i

√2





0 −1 0

1 0 1

0 −1 0



 qˆxz = 1

√2





0 1 0

1 0 −1

0 −1 0



 qˆxy= 1

√2





0 0 −i

1 0 0

i 0 0





ˆ q_x2−y² =





0 0 1 0 0 0 1 0 0



 qˆ_zz=1

3





−1 0 0

0 2 0

0 0 −1



 I =





1 0 0 0 1 0 0 0 1





2.2 Single-particle Hamiltonian

The fundamental Hamiltonian of a spinor BEC can be constructed quite generally based on the symmetry argument. We consider a system of identical bosons with mass M and spin F that are described by the field operators ˆψm(r), where m = −F, F + 1, ...F , denoting the magnetic quantum number. Here, F is the hyperfine spin of an atom, that is, the composition of the the electron orbital angular momentum L, electron spin S, and the nuclear spin I. At low magnetic fields, the hyperfine interaction between these spins dominates, and the ground-state subspace breaks into manifolds of states with definite total spin F= I ± 1/2. With a nuclear spin of I =3/2, ⁸⁷Rb has either F=1 or F=2 in the J = 1/2 electronic ground state.

In this work, we are focusing on modelling a spin-1 BEC, where in we are specifically considering the F=1 hyperfine manifold.

The non-interacting part of the Hamiltonian with the kinetic term and the trapping potential Utrap(r), without the Zeeman terms contributed by the external magnetic field is written as,

Hˆ0= Z

dr

1

X

m,m⁰=−1

ψˆ^†_m(r) −~²∇²

2M + Utrap(r)



ψˆm⁰(r) (2.2)

2.3 Interaction Hamiltonian

Here, we again start with the case where there are no external applied fields. We consider the low density of the atomic gases where the interactions result almost exclusively from binary collision, and thus even in a many-body system the interactions can be treated through a two-body scattering problem. In this scenario, the interaction between two atoms conserves the total angular momentum, denoted by the quantum number, Tpair - this includes both the internal angular momenta of the two atoms and the orbital angular momen- tum of their relative motion. In case of low energy collisions between two atoms in short-range two-body potential, the resulting collisions will be in the s-wave channel. This results in the lowest order partial wave approximation, where the relative orbital angular momentum of the two interacting atoms Lpair = 0. [6, 12]

Thus the total angular momentum is just the sum of the two hyperfine spins, F_pair = F₁+ F₂ (F_pair is the spin quantum number). Collisions between identical bosons preserve exchange symmetry, implying that the total angular momentum of the colliding atoms take only even integers in the s-wave regime, ie., F_pair ∈ {0, 2}. Following the arguments in [11], we can write the collisional interaction in the form,

Vˆint= X

Fpair∈{0,2}

4π~²

M aF_pairPˆF_pairδ(r) (2.3)

where M is the mass of the single atom, r is the relative position between the atoms and aFpair is the scattering length. The projection operator ˆP_Fpair projects the pair of atoms into the total spin-F_pair state.

It is convenient to express (Eq. (2.3)) in terms of the spin operators. Considering spin ( ˆF) and identity ( ˆI) operator identities [6],

[ ˆI1⊗ ˆI2]S = ˆP0+ ˆP2 (2.4) [ ˆF₁· ˆF₂]_S = ˆP₂− 2 ˆP₀

where the subscripts {1,2} represent the two colliding atoms and the subscript S indicates that we are including only those states that symmetric under exchange. (Eq. (2.3)) then becomes,

Vˆint=1 2



c0[ ˆI1⊗ ˆI2]S+ c2[ ˆF1· ˆF2]S



(2.5) where the interaction potential is grouped as the spin independent part parametrized by c₀ and the spin dependent part parametrized by c₂, and they are given as,

c0= 4π~²

3M (a0+ 2a2), c2= 4π~²

3M (a2− a0) (2.6)

For the F=1 hyperfine manifold of ⁸⁷Rb, (a₀, a₂) are given as (101.8a_B, 100.4a_B) respectively [6], where a_B is the Bohr radius. c₂ > 0 correspond to anti-ferromagnetic interactions and c₂ < 0 correspond to ferromagnetic interactions, which is our case in F=1 hyperfine manifold of ⁸⁷Rb.

Considering the second quantization form, each identity operator becomes the number-density operator given as, ˆn(r) =P

mψˆ_m^†(r) ˆψm(r) and the spin vector is written as ˆFk =P1

m,m⁰=−1ψˆ_m^† (r)(f_(k)mm0) ˆψm⁰(r) (where k ∈ x, y, z). Thus ˆV_int can be written as,

Vˆint= Z

dr c0

2 : ˆn²: +c2

2 : ˆF1· ˆF2:



(2.7)

where for any operator ˆO, : ˆO: denotes normal ordering. Considering such interactions, we can re-write the general Hamiltonian as, ˆH0+ ˆVint

H =ˆ Z

dr

1

X

m,m⁰=−1

ψˆ_m^† (r) −~²∇²

2M + Utrap(r)



ψˆm⁰(r) + Z

dr c₀

2 : ˆn²: +c₂

2 : ˆF1· ˆF2:



(2.8)

Thus the two terms in the 1st integral are the single particle contributions - kinetic and potential energy, and the terms in the second integral account for the density and the spin-spin interactions.

2.4 Zeeman Interaction

We shall now consider the system in the presence of an applied external magnetic field. In this condition, we get the Hamiltonian contribution from the Zeeman terms as,

Hˆ_z≡ ˆH_LZ+ ˆH_QZ (2.9)

Hˆz= Z

dr

1

X

m,m⁰=−1

ψˆ^†_m(r)



p( ˆfz)mm⁰+ q( ˆf_z²)mm⁰

 ψˆm⁰(r)

The two terms on the right side are the linear Zeeman shift (LZS) and quadratic Zeeman shift (QZS) contributions respectively. fˆz is the z-component of the spin matrix, and p = gµBB and q = ^(gF_∆E^µBB)2

hf

are the linear and quadratic Zeeman terms. For the F=1 hyperfine manifold of ⁸⁷Rb the Land´e g-factor gF = −1/2 [13, 14], µBis the Bohr magneton, energy splitting between the hyperfine states ∆Ehf ' 6.8GHz.

We have taken the axis of quantization and the direction of applied magnetic field Bzz to be along the z-axis. The representation basis for the spin state is generally chosen for such quantization to be along z, and therefore | + 1i, |0i, | − 1i are the eigenstates of the spin operator ˆf_z. The LZS is a constraint on the longitudinal magnetization. A spinor gas can evolve by the spin-spin collisions that contribute to mixing of the hyperfine levels(spin mixing collisions). This process is insensitive to the LZS, but depends on the QZS, which shifts the average energy of the |m = ±1i states by an energy q with respect to that of the |m = 0i state. QZS comes about from the fact than an applied magnetic field mixes hyperfine levels, causing energy repulsion between states in, say, the upper and lower hyperfine states with the same value of m in the electronic ground state of alkali gases. It is good to note here that the Zeeman shifts are in general symmetry-breaking external influences. While they still preserve the symmetry of rotations about the longitudinal direction (z axis), more generally, the linear and quadratic Zeeman energies could, in principle, be applied along different axes, breaking rotational symmetry altogether [6]. Thus (Eq. (2.8)) becomes,

H =ˆ Z

dr

1

X

m,m⁰=−1

ψˆ_m^† (r) −~²∇²

2M + Utrap(r) + p( ˆfz)mm⁰+ q( ˆf_z²)mm⁰



ψˆm⁰(r) + Z

dr c₀

2 : ˆn²: +c₂

2 : ˆF1· ˆF2:



(2.10)

2.5 Mean-field Description

Following (Eq. (2.10)) we see that studying the spin dynamics in a gas with larger particle number, say 10⁴ as in our case, is generally far more complex and computationally expensive as we need to solve the many- body Schr¨odinger equation. In dilute gases this many-body problem can be simplified by considering the effect of individual particles as an average or mean action of the fluid, thus in this way ignoring quantum fluctuations and correlations between particles. This is obtained by replacing the field operators with expectation values,

ψˆm(r) →√

N ψm(r) (2.11)

ψˆ_m^†(r) →√

N ψ_m^∗(r)

Thus we get a condensate wavefunction defined to be a complex number with a definite phase. We get factor

√

N as we normalize the single-particle wavefunction to unity. In this picture, we assume that the Bose- condensed bosons occupy single spatial mode. We also note here that a spinor condensate is a system where the atoms, although in the ground state, have the possibility to occupy different internal states and go from one state to another; in contrast to mixtures of several species where the populations in each component are strictly conserved. This implies that a spinor condensate cannot be described by a scalar field operator ψ(r), and instead calls for a vectorial field operator with 2F+1 components. For F=1 hyperfine groundstateˆ of ⁸⁷Rb , we have ˆψ(r) = ( ˆψ+1(r), ˆψ0(r), ˆψ₋₁(r))^T In mean-field description we write this as,

ψ(r) = h ˆψ(r)i = (ψ₊₁(r), ψ₀(r), ψ₋₁(r))^T (2.12) Thus, the interaction between the particles are described by an order parameter or a macroscopic wave- function given by the average inter-particle interaction.

2.6 Gross-Pitaevskii Equation

Using the mean-field description and ˆH developed so far, we can get the energy functional as, (in this section we mainly follow [13])

E(ψ_m) ≡ h ˆHi = Z

dr

" ₁ X

m=−1

ψ^∗_m −~²∇²

2M + U_trap(r) + pm + qm²



ψ_m+c₀ 2n²+c₂

2|F|²

#

(2.13)

where, ( ˆf_z)_mm⁰ = mδ_mm⁰ and ( ˆf_z²)_mm⁰ = m²δ_mm⁰. The local density,

n(r) = hˆn(r)i =

1

X

m=−1

|ψm(r)|² (2.14)

The spin density vector,

Fk(r) = h ˆFk.(r)i =

1

X

m,m⁰=−1

ψ^∗_m(r)(F(k)mm⁰)ψm⁰(r) (2.15)

where k ∈ {x, y, z}. The components of the spin vector F is given as, Fx= 1

√

2[ψ₁^∗ψ0+ ψ₀^∗(ψ1+ ψ−1) + ψ^∗₋₁ψ0] (2.16) F_y= i

√2[−ψ^∗₁ψ₀+ ψ^∗₀(ψ₁− ψ−1) + ψ₋₁^∗ ψ₀] (2.17)

F_z= |ψ²₁| − |ψ₋₁² | (2.18)

The time evolution of the mean field is given as,

i~∂ψm(r)

∂t = ∂δE

∂δψ_m^∗(r) (2.19)

Upon substituting the energy functional value(2.13) in (2.19), we get,

i~∂ψm(r)

∂t = −~²∇²

2M + U_trap(r) + pm + qm²



ψ_m+ c₀nψ_m+ c₂

1

X

m⁰=−1

F.f_mm⁰ψ_m⁰ (2.20)

(where m= (+1,0,-1)). (Eq. (2.20)) is the three component Gross-Pitaevskii equation of the spin-1 Bose Einstein Condensate. For the stationary state, we substitute,

ψ_m(r, t) = ψ_m(r)e^iµt/~ (2.21)

in (Eq. (2.20)), where µ is the chemical potential. We then get,

 −~²∇²

2m + Utrap(r) + pm + qm²



ψm+ c0nψm+ c2 1

X

m⁰=−1

F.fmm⁰ψm⁰ = µψm (2.22) Writing down the three components m = 1, 0, -1 explicitly, we obtain,

 −~²∇²

2m + Utrap(r) + p + q + c0n + c2Fz− µ



ψ1+ c2

√2F₋ψ0= 0 (2.23)

c₂

√2F₊ψ₁+ −~²∇²

2m + U_trap(r) + p + q + c₀n − µ



ψ₀+ c₂

√2F₋ψ₋₁= 0 (2.24)

c2

√2F+ψ0+ −~²∇²

2m + Utrap(r) − p + q + c0n + c2Fz− µ



ψ₋₁= 0 (2.25)

where F_±≡ Fx± iFy. Solving Eq. (2.23), Eq. (2.24) and Eq. (2.25) we investigate the properties of ground states and excited states in a spin-1 GPE.

2.6.1 Single-Mode Approximation

As discussed before in Section 2.5, one of the salient features that distinguishes a spinor BEC from a binary mixture of BECs is the spin-mixing dynamics. The initial population balance can change in a spinor BEC via the spin exchange collision, whereas it is fixed in a mixture. In a spin-1 BEC, two atoms in the magnetic sublevel m = 0 can coherently and reversibly scatter into a pair of atoms in the m = +1 and m = −1 states, and vice versa. To simplify this scenario, we consider the single mode approximation (SMA) - where we assume that all the spin components share the same spatial dependence, and that they vary only time.

So far we have considered groundstate in a uniform system, whereas to consider experiments, where there are confining optical traps, SMA is usually adopted to study spin dynamics. Hence the vectorial order parameter in SMA,

ψm(r, t) =√

N ζm(t)ψSMA(r)e^iµt/~ (2.26)

where ζ_m(t) is the space-independent normalized spinor. Under SMA, we substitute (Eq. (2.26)) in (Eq. (2.20)), and obtain the equation of motion of ζ_mas,

i~∂ζ_±

∂t = (±p + q)ζ_±+ ˜c2[(ρ_±+ ρ0− ρ_∓)ζ_±+ ζ₀²ζ_∓^∗] (2.27) i~∂ζ0

∂t = ˜c2[(ρ1+ ρ₋₁)ζ0+ 2ζ0ζ₋₁ζ₀^∗] (2.28) where ˜c₁ = c₁/V^eff, V^eff = (R dr|ψSMA|⁴)⁻¹, is the effective volume of the system and ρ_m = |ζ_m|². Thus under SMA, we keep track of two vital quantities, the number conservation N and the longitudinal magnetization M, (ie., M_(z)=R drFz(r)) while studying the spin mixing dynamics,

N = ρ1+ ρ0+ ρ₋₁ (2.29)

M = ρ₁− ρ₋₁ (2.30)

In the sections to follow, we will study the ground state solutions by minimizing the energy functional Eq. (2.13). An important aspect of that process is monitoring these constraints, and we shall now make a few comments on them. In a typical experiment, the last stage before condensation have atomic evaporations during which neither N nor M is conserved. In case of a scalar condensate, a ground state is obtained under only N conservation, which gives rise to the the chemical potential µ we encountered earlier, which is nothing but the Lagrange multiplier of this constrained minimization. Typically in case of the vector condensate we account for two Lagrange multiplier terms, µ and λ for N and M respectively [15]. While is N is always conserved, conservation of M is disturbed under the presence of externally applied fields(ref Section 4.2). Thus SMA is not exact even in ground states, when the linear and quadratic Zeeman effects are not compatible with the magnetism determined by the spin-dependent interactions [16].

Chapter 3

Numerical methods for simulating spin-1 BECs

Proceeding from the theoretical background we have obtained, we will now start developing the numerical framework to solve GPE, to study the ground state under varied conditions. We will employ the numerically developed MatLab opensource package GPELab [17] for our purposes. Following [17] In this section we will essentially build the numerical solver scheme for single component for simplicity, and then scale it to the multi-component system at the end which we employ for our studies.

3.1 Dimensionless GPE

As a general prescription in the numerical computations, we will non-dimensionalize the quantities involved, and we use the following time, length and energy scales.

t → t

ω_ho, r → raho, aho=

r ~

M ω_ho, E → ~ω^hoE, ψ → ψa^3/2_ho/√

N (3.1)

where ω_ho and a_ho are the harmonic oscillator frequency and length scale respectively, M is the atomic mass, and N is the total number of atoms in the system. We then get the dimensionless energy functional,

End(ψ) = Z

R²

h1

2|∇ψ|²+ V |ψ|²+β 2|ψ|⁴i

dr (3.2)

and the dimensionless GPE,

i∂ψ

∂t =



−1

2∇ + V + β|ψ|²



ψ (3.3)

β is the interaction strength, which in case of our 3 component system will be given by the c₀and c₂terms.

Since we solve for pancake-shaped condensates, we are solving for two-dimensional GPE - where we assume that the condensate is confined in the z direction., ie., excitations along the x and y axes require less energy than in the z direction [18]. Hence the potential,

V (r) =1

2(γxx²+ γyy²) (3.4)

where γx= ωx/min(ωx, ωy), γy= ωx/min(ωx, ωy) Under the two-dimensional system and symmetric har- monic frequency, γx≡ 1 and γy ≡ 1. Thus we have decomposed the three dimensional wavefunction Ψ, in x, y and the third z dimensions as,

Ψ(r, t) = ψ(x, y, t)ϕ(z) (3.5)

where,

ϕ(z) =

 Z

R²

|Ψ0(x, y, z)|²dxdy

^1/2

(3.6)

with Ψ0 being the 3D stationary state; as Ψ0 is normalized we have Z

R

|ϕ(z)|²dz = 1 (3.7)

3.2 Stationary States

The stationary states of a quantum system are the eigenfunctions of the Hamiltonian describing the system and the eigenvalues are the quantified energies associated with the Hamiltonian. Thus we essentially find the solution to,

ψ(r, t) = e^−iµtφ(r) (3.8)

where µ is the chemical potential of the condensate and φ is the time independent function, such that, it satisfies the normalization constraint, the L²-Norm

||φ||²₀= Z

R²

|φ(r)|²dr = 1 (3.9)

Under such constraint, the solution is given as the solution to the equation, µφ(r) = −1

2∇φ(r) + V (r)φ(r) + βφ|(r)|²φ(r) (3.10) This is a nonlinear eigenvalue problem and is solved by computing the chemical potential µ The energy functional (3.2) becomes,

End(φ) = Z

R²

h1

2|∇φ|²+ V |φ|²+β 2|φ|⁴i

dr (3.11)

Thus the eigenfunctions are the critical points of the energy functionl End(φ) [19]. The global minima of End(φ) under the normalization constraint gives us the ground state solution.

3.3 Conjugate Normalised Gradient Flow

There is a vast class of minimization methods developed in the literature for solving such nonlinear eigen- value problems [20, 21, 22]. However, as developed in [23] we use the approach of Conjugate Normalized Gradient Flow (CNGF) or commonly known as the imaginary time evolution [24, 23], which consists in building a minimizing sequence of the energy functional End(φ). The algorithm to compute the sequence of iterates (φ(x, tn))_n∈N in CNGF is given by,























∂tφ(r, t) = 1

2∇φ(r, t) − V (x)φ(x, t) − β|φ|²φ(r, t), (tn< t < tn+1) φ(r, tn+1) = φ(r, t⁺_n+1) = φ(r, t⁺_n+1)

||φ(r, t⁺_n+1)||₀, φ(r, 0) = φ₀(r), r ∈ R², (with||φ||₀= 1)

(3.12)

Here we have time discretised the system, with (t_n)_n∈N, t₀ = 0, and with the defined local time step δt_n = t_n+1− t_n, ∀n ∈ N. As we can observe, this method involves a good choice of initial data φ0 as

’Ansatz’ of the underlying GPE, in an iterative solver, to get the convergence (See Section 3.3.2).

3.3.1 Time and Space discretizations

Semi-Implicit Backward Euler Scheme in Time

Different discretisation schemes have been explored over time for the ground state computations. Following the work of Bao and Du in [23] we see that the semi-implicit Euler scheme is a good discretisation scheme

for CNGF, as it leads to a decaying modified energy End(φ) at each step of the projected steepest descent algorithm. Accordingly (Eq. (3.12)) modified as,























∂_tφ(r, t) = 1

2∇φ(r, t) − V (r)φ(r, t) − β|φ|²φ(r, t), (t_n< t < t_n+1) φ(r, t_n+1) = φ(r, t⁺_n+1) = φ(r, t⁺_n+1)

||φ(r, t⁺_n+1)||0

,

φ(r, 0) = φ0(r), r ∈ R², (with||φ||0= 1)

(3.13)

Remembering the fact that the L²-Norm is strictly preserving and that there is a minimizing sequence of energy functional,we get that,

t→∞lim φ(r, t) = φg(r) (3.14)

where φg(r) is a stationary state. The long time computation here is fixed according to a stopping criterion, given as,

||φⁿ⁺¹− φⁿ||_∞< δt (3.15)

where ||.||∞ is the uniform norm. The choice of  should be small enough to obtain a good accuracy of stationary state, particularly when we consider a highly accurate solutions based on pseudo-spectral approximation techniques.

Backward Euler psesudoSPectral Scheme in Space(BESP)

The CNGF scheme requires spatial descritizations and two widely used schemes are Backward Euler Finite Difference(BEFD) and Backward Euler psesudoSPectral(BESP) schemes. Here we employ the psesudoSPec- tral approximation for the spatial derivatives for their higher order accuracy in space[18]. This scheme is based on FFTs. We consider a large enough computational box : O :=] − a^x; ax[×] − ay; ay[. We define the spatial grid indices as (xj, yk) with (j, k) ∈ I^J,K,

IJ,K = (j, k) ∈ N²; 0 ≤ j ≤ J − 1 and 0 ≤ k ≤ K − 1

(3.16) with J, K ≥ 2 and uniform discretisation steps h_xand h_yin x, y directions respectively. Given that this is a Fourier Series based method with FFTs, we have imposed a periodic boundary condition, on the fictitious boundary of O. The partial Fourier psesudoSPectral discretizations in x, y directions are given as,

φ(x˜ i, yj, t) = 1 J

J/2−1

X

p=−J/2

c˜

φp(yk, t)e^icp(xj+ax), φ(x˜ i, yj, t) = 1 K

K/2−1

X

q=−K/2

c˜

φq(xj, t)e^icq(yk+ay) (3.17)

where cφ˜p, cφ˜q are the respective Fourier coefficients in x and y directions,

c˜

φp(yk, t) =

J −1

X

j=0

φ(x˜ j, yk, t)e^{−icp(xj+ax)}, cφ˜q(yk, t) =

K−1

X

k=0

φ(x˜ j, yk, t)e^{−icq(yk+ay)} (3.18)

where cp= ^πp_a

x, cq =^πq_a

y. Then the first and second order spatial derivatives along x and y is given as, ([∂_x])_j,k= 1

J

J/2−1

X

p=−J/2

ic_pφc˜_p(y_k, t)e^icp(xj+ax), ([∂_x])_j,k= 1 K

K/2−1

X

q=−K/2

ic_qcφ˜_q(x_j, t)e^icq(yk+ay)

([∂_x²])_j,k= 1 J

J/2−1

X

p=−J/2

−c²_pcφ˜_p(y_k, t)e^icp(xj+ax), ([∂_x²])_j,k= 1 K

K/2−1

X

q=−K/2

−c²_qcφ˜_q(x_j, t)e^icq(yk+ay)

(3.19)

Accordingly Laplacian, potential and nonlinear operators become,

([∆])_j,k= ([∂_x²] + [∂_y²])_j,k, ([V ])_j,k= V (x_j,k), ([f (|φⁿ|²)])_j,k= f (|φⁿ_j,k|²) (3.20) Thus the CNGF scheme (Eq. (3.12)) the produces a sequence of vectors solution to







Aⁿφ = b˜ ⁿ φⁿ⁺¹= ^φ˜

|| ˜φ||₀

φ⁰:= φ0

(3.21)

where,

Aⁿ= −1

2[[∇]] + [[V ]] + [[β|φⁿ|²]] (3.22) bⁿ:= φⁿ

δt (3.23)

||.||0 here is the L²-Norm in the grid O for each vector φ,

||φ||0= h^1/2_x h^1/2_y X

(j,k)∈IJ,K)

|φj,k|²
1/2

(3.24)

Finally the discrete energy functional is given as, E(φ) = h^1/2x h^1/2_y X

(j,k)∈IJ,K)

R

 φ^∗_j,k



−1

2[∇]φ + [V ]φ + [f (|φ|²)]φ



j,k



(3.25)

3.3.2 Initial Condition for CNGF and Iterative Solvers

As discussed earlier, to obtain solution to the minimization problem (Eq. (3.11)) an initial guess has to be given to the iterative method to initialize it and then the minimization process tries to compute a minimal solution through iterations. Considering the scheme (Eq. (3.21)) we observe that at each iteration n, the minimization method requires a solution of the linear system Aⁿφ = b˜ ⁿ. We use a pseudoSpectral approx- imation method, where operator Aⁿ is given implicitly through a FFT. Thus a direct matrix solver cannot be used in this context. And also the matrix solvers are computationally expensive, given the process of computing inverse of the matrix that is involved. An alternative method is using a matrix-free iterative solvers method. A common iterative solver employed is the Krylov subspace iterative method, specifically the type BiConjugate Gradient Stabilized (BiCGStab) method in our case [25].The introduction of iter- ative solvers accelerate our computation by simple operator-based preconditioners, and can be efficiently implemented with the help of bicgstab() routine available in MatLab. Since the numericals of iterative solvers are beyond the scope of this thesis, the reader is guided to refer to [26] for more details.

For any numerical computation based on iterative schemes, it is quite natural to have an initial condi- tion that is closer to the solution of the equation that we look for, GPE in this context. More on the prescription of the initial conditions are available in the physics literature [19, 27]. We are employing the Gaussian function as the initial data [17] for our solver, given as,

φ0(x) = (γxγy)^(1/4)

√π e^{−(γxx2+γyy2)/2} (3.26)

(where γ_x, γ_y equals to 1 in our case)

3.4 Extension to the three component case

The same scheme of numerical computations apply to the three components system too, where we will define the Laplacian, potential and nonlinearity as 3×3 matrices as follows: Diagonal Laplacian,

∇ψ(r, t) → (∇ψn(r, t))n=1,2,3 (3.27)

Diagonal potential matrix (diagonal as we consider a symmetric trap),

V (r) → (Vnn(r))n=1,2,3 (3.28)

Non-linearity matrix,

f (|ψ|²)ψ →





f₁₁ f₁₂ f₁₃ f₂₁ f₂₂ f₂₃ f₃₁ f₃₂ f₃₃



ψ (3.29)

Looking at coupled GPE equations, (Eq. (2.23), Eq. (2.24), Eq. (2.25)), we get the nonlinearity matrix as,

f (|ψ|²)ψ →







(−p + q + C0n + C2Fz) ^√C2

2F₋ 0

C₂

√2F₊ C₀n ^C√2

2F₋

0 ^√C2

2F+ (p + q + C0n − C2Fz)









 ψ₊₁

ψ₀ ψ₋₁



 (3.30)

We can obtain (Eq. (3.29)) by scaling dimensions in corresponding values of (Eq. (3.30)). The non-linearity enters here via n from Eq. (2.14). The linear and quadratic Zeeman terms are scaled as

p → p

~ωho

, q → q

~ωho

(3.31) The scaled interaction constants are given in two-dimensions as [28],

c0→ N c0

a³_ho~ωho

=4πN (a0+ 2a2) 3aho

√1

2π (3.32)

c2→ N c2

a³_ho~ωho

= 4πN (a2− a0) 3aho

√1

2π (3.33)

Chapter 4

Ground State Computation for Static Magnetic Fields

In this chapter, we will start testing the numerical framework. All the results presented here are for 10⁴ number of atoms, with a trapping frequency of 2π×50 Hz. We will start with testing for different external bias field conditions and look at the effects of linear and quadratic Zeeman shifts on the ground state solutions, and discuss the choice of computational parameters for the simulations.

4.1 Numerical ground state solutions at B=0

The simple case to start with is without the Zeeman term contributions- ie., no applied external fields.

We consider a symmetric trapping frequency of 2π×50 Hz. We start with the Gaussian function as the initial ansatz as described earlier, (Fig. 4.2a) for initiating the numerical iterator. We find the ground state solutions by letting the GPE propagate in the imaginary time evolution, sufficiently long enough, terminating with a suitable stopping criterion. We obtain a Thomas-Fermi density profile for the condensate as expected (Fig. 4.2b). In the case where the spin dependent collisional constant c2= 0, SMA is strictly valid in the ground states. The computations are performed in two dimensions Fig. 4.1.

(a) m=+1 (b) m=0 (c) m=-1

Figure 4.1: Ground states density distributions of individual spin states at the end of imaginary time evolution in two dimensions - with Zeeman terms p=0,q=0, and c2=0. Parameters of spatial extent in x and y axes directions, and the colour scale indicating the density distribution are scaled as given in Eq. (3.1).

Spatial extent here indicates the choice of computational box size. We have a condensate of Thomas Fermi radius ≈ 10µm. Without the c₂ term the states are equally populated as observed from the density scale.

-10 -8 -6 -4 -2 0 2 4 6 8 10 x

0 0.02 0.04 0.06 0.08 0.1 0.12

|(x,y=0)|2

m= +1 m= 0 m= -1

(a)

-10 -8 -6 -4 -2 0 2 4 6 8 10

x 0

0.005 0.01 0.015

|(x,y=0)|2

m= +1 m= 0 m= -1

(b)

Figure 4.2: Density profiles plotted against spatial dependence - 1D plot of Fig. 4.1 to illustrate the relative density profiles of individual states - (Fig. 4.2a): Initial Condition φ₀chosen as ansatz for initiating computation (Fig. 4.2b): Ground state distribution at end of imaginary time evolution for a system with Zeeman terms- p=0, q=0 and constant of spin dependence c₂=0

In the case when we consider the c2 term too- ie., the contribution from the asymmetric part of the Hamiltonian, we observe that the validity of SMA is disturbed, but the symmetry between the |m = +1i and |m = −1i states are preserved. We observe that the collisional interactions result in a state with the populations favouring the |m = 0i state as in Fig. 4.3, Fig. 4.4. We will observe that this spin exchange phenomena, will become more pronounced as we include the effect of the external magnetic field.

(a) m=+1 (b) m=0 (c) m=-1

Figure 4.3: Ground states distributions of individual spin states in two dimensions - with Zeeman terms p=0,q=0, and spin dependent constant c₂<0 for the ferromagnetic condensate of F=1 hyperfine manifold.

We note from the density scales that the |m = 0i state is favoured.

-10 -8 -6 -4 -2 0 2 4 6 8 10

x 0

0.005 0.01 0.015 0.02 0.025

|(x,y=0)|2

m= +1 m= 0 m= -1

Figure 4.4: Ground state distributions in Fig. 4.3 plotted in one dimension. As opposed to Fig. 4.2b, there is difference in amplitude between the |m = 0i and |m = +1i, |m = −1i states.

4.2 Numerical ground state solutions at B6=0

4.2.1 With non-zero linear and quadratic Zeeman terms

We will now consider the case where we have an externally applied homogeneous magnetic field along the z direction. In an uniform bias field, we get a net magnetization along the field direction, resulting in a shift in population to the |m = +1i level, ie.,in a (1,0,0)^T state. This is characterised as a ferromagnetic phase, resulting in a net longitudinal magnetization, Mz =R drFz(r), of +1, as we observe in (Fig. 4.8), under a static B field of 0.1 Gauss along the z direction.

(a) m=+1 (b) m=0 (c) m=-1

Figure 4.5: Ground state solutions of individual spin states in two dimensions - computed with linear(p) and quadratic(q) Zeeman terms at Bz=0.1 Gauss

-10 -8 -6 -4 -2 0 2 4 6 8 10

x 0

0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045

|(x,y=0)|2

m= +1 m= 0 m= -1

Figure 4.6: Ground state distributions in Fig. 4.5 plotted in one dimension - note that the application of the magnetic field overwhelms the atomic mean-field interactions and the ground state corresponds to all the atoms condensing into the |m = +1i Zeeman sublevel. (the red and yellow lines corresponding to |m = 0i and |m = −1i respectively overlap each other at the bottom)

The presence of B-field enters the system through the Zeeman Hamiltonian with linear and quadratic Zeeman terms as,

p = gµBB · ~~ F = gµB





Bz 0 0

0 0 0

0 0 Bz



 (4.1)

q = (gµB)²( ~B · ~F )²

∆Ehf

= (gµB)²

∆Ehf





B_z² 0 0

0 0 0

0 0 B_z²



 (4.2)

which is obtained by considering the ˆf_zspin-1 matrix , as the external homogeneous field is along z.

We will discuss a bit more on the role of these two terms. The Zeeman dynamics of the atoms in the magnetic field consists of the Larmor precession given by the linear zeeman term. This Larmor oscillation amplitudes are modulated by the qudratic Zeeman shift term, q. Fig. 4.7 shows the relation between the linear and quadratic Zeeman energies for different magnetic field values, illustrating that the linear term is over four orders of magnitude higher than the qudratic part in lower order magnetic fields. Thus the presence of the linear Zeeman term, results in the system being longitudinally magnetized, ie., we get a spinor condensate in the mean-field ground state of (1, 0, 0)^T with the energy of the system minimized when all the atoms are aligned along the magnetic field. While the addition of external fields contribute to symmetry breaking mechanisms, the spin rotational symmetry is still preserved along the z-axis, the axis of quantization.

10^-3 10^-2 10^-1 10⁰ 10¹ 10²

B field (Gauss) 10^-4

10^-2 10⁰ 10² 10⁴ 10⁶ 10⁸ 10¹⁰

Zeeman Shifts(Hz)

p q

Figure 4.7: Linear and quadratic Zeeman effects for different magnetic fields B for ⁸⁷Rb . We observe that the linear Zeeman shift is higher than quadratic Zeeman shift by over four orders of magnitude in magnetic fields lower than 1 Gauss

We note that in this case, we allow for the system to be magnetized, accounting for a non M conserving process observed incase of free propagation into ground state solution. Incase of studying processes with a given magnetization, Mz=R drFz(r), where magnetization is conserved, we account for the conservation in computations by introducing the Lagrange multiplier λ in Eq. (2.19) as E − λMz, as discussed earlier. Thus we modify µ in Eq. (2.23) and Eq. (2.24) as (µ − λ) and (µ + λ) respectively. This constraint is essentially applied in studying systems where we consider an M conserving evaporation process, and studying the ground state solutions when a condensate is subjected to external manipulations conserving both N and M [29].

4.2.2 With linear Zeeman term p = 0 and non-zero quadratic Zeeman term

We shall now consider the case, where the linear Zeeman term is set 0 and there exists a quadratic Zeeman term contribution. In case of F=1 ⁸⁷Rb , we have c2 < 0 and q > 0. Under this case, we observe the atoms favouring |m = 0i state, with the longitudinal magnetization remaining zero, owing to the symme- try between the +1 and -1 states - ie., the quadratic Zeeman effect causes minimization of total energy with two atoms in |m = 0i state created for every collision of |m = +1i atom with that of a |m = −1i atom.

When p=0, the atomic ensemble is magnetized perpendicular to the direction of the applied magnetic field, and its amplitude monotonically varies as a function of q, until |m = 0i is maximally populated. This phase is characterized as the polar phase [13]. Fig. 4.8 shows the changing population between between the |m = 0i and |m = +1i,|m = −1i states with changing values of the B-field along z. Fig. 4.9 shows the different ground state energies obtained for varying B-field values, with the trend agreeing with the quadratic shift contribution.

Setting p = 0, is essentially equivalent to considering the evolution of a spinor Bose gas in a frame rotating with Larmor precession frequency ωL = −γBz (where γ is the gyromagnetic ratio). This is equivalent to a typical experimental condition, where the atoms prepared in (1, 0, 0)^T state, defined by the quantization

axis along a set bias field is rotated with a π/2 radio-frequency(RF) pulse to a state (1/√

2, 1/2, 1/√ 2)^T. This way, we assume a fixed magnetization and we study the spin dynamics in the transverse plane.

-10 -8 -6 -4 -2 0 2 4 6 8 10

x 0

0.005 0.01 0.015 0.02 0.025 0.03

|(x,y=0)|2

m= +1 m= 0 m= -1

(a) Bz=0.6 Gauss

-10 -8 -6 -4 -2 0 2 4 6 8 10

x 0

0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04

|(x,y=0)|2

m= +1 m= 0 m= -1

(b) Bz=1 Gauss

-10 -8 -6 -4 -2 0 2 4 6 8 10

x 0

0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045

|(x,y=0)|2

m= +1 m= 0 m= -1

(c) Bz=1.2 Gauss

Figure 4.8: Ground states distributions without linear Zeeman shift term p, with non-zero q value simulated at different external magnetic fields ( one-dimensional illustration of computations done in two-dimensions).

We observe increase in population in the |m = 0i state with increase in magnetic fields

0 0.2 0.4 0.6 0.8 1 1.2

B(Gauss) 6.8

6.9 7 7.1 7.2 7.3 7.4 7.5 7.6 7.7

Ground State Energies

Figure 4.9: Variation in Ground state energies obtained for different values of external magnetic fields - the energy values are scaled as described in Eq. (3.1). We observe the quadratic trend owing two the q term as expected

4.3 On the choice of computational parameters

In this section we will make a comment on the choice of computational parameters such as time and space discretization parameters and stopping criterion for simulations. Fig. 4.10 shows the ground state energies for the different choice of time and space discretizations and stopping criteria.

The choice of such parameters depend on the type of physics we are simulating to obtain ground state solutions at. A very small grid size, as in cases shown in Fig. 4.10c, might result in over sampling of a system, where the Fourier components at each step in an iterative subspace results an overly converged system, costing more computational time [18]. On the other hand a good spatial resolution is required to simulate situations when parameters involved are spatially dependent (say in case of gradient or inhomo- geneous magnetic field which we will explore in the following chapter). A loose stopping criterion constant

, means that the we have not let the system globally converge enough to its stationary state(Fig. 4.10b).

Larger time steps might not resolve the physics under study, ie., for example, in a case where we compute the ground state with both linear and quadratic Zeeman terms, a sufficiently smaller time step is required to resolve the contribution from p that accounts for rapid Larmor precessions and q resolves the mean-field dynamics; whereas in a case with only q, we do not resolve Larmor precessions, but study the mean-field dy- namics directly. With these in perspective, we have performed the computations in this section for a 2⁷× 2⁷ number of grid points, in a computational box O := [−10, 10] × [−10, 10] in x and y (ie.,dx=dy=0.1563,

obtained by dividing the spatial extent by the number of grid points), dt = 10⁻³ for computations without p term, and dt = 10⁻⁴ for computations with p term, and a stopping criterion ( dt) with  = 10⁻³.

10^-4 10^-3 10^-2 10^-1 10⁰

dt 6.867

6.8675 6.868 6.8685 6.869 6.8695

Ground State Energies

(a) dt vs E

(at  = 10⁻³,dx=dy=0.1563)

10^-4 10^-3 10^-2 10^-1 10⁰ 6 7 8 9 10 11 12 13 14

Ground State Energies

(b)  vs E

(at dt=10⁻³,dx=dy=0.1563)

0 0.2 0.4 0.6 0.8 1 1.2 1.4

dx,dy 4.5

5 5.5 6 6.5 7 7.5 8

Ground State Energies

(c) dx,dy vs E (at dt=10⁻³, = 10⁻³)

Figure 4.10: Ground state energies plotted for different choices of computational parameters - (Fig. 4.10a)time step, (Fig. 4.10b)stopping criterion constant and (Fig. 4.10c)spatial step respectively.

Plots indicate the trend the ground state energies obtained when each parameter is weighed for differed values, keeping the other two parameters constant. Choice of optimal computational parameters is made looking at the convergence obtained.