Dipolar potential

with

A˜^kl =A^kl+ Σg, (3.44)

˜

p^kl =p^kl−2Σgq_g, (3.45)

˜

γ^kl =γ^kl+q_g^TΣgq_g, (3.46) and where ˆv_g,Σ_g,q_g are the amplitude, width matrix and center coordinate of the Gaussian beam, respectively. Gaussian laser beams can be used for instance in the generation of multi-well potentials, where a superposition of several Gaussian beams is deployed. For the experimental setup see reference [76]. We use these Gaussian potentials in the description of the triple-well in chapter 5 and of the complex double-well in chapter 6, where ˆv_g is a complex number to describe gain and loss of particles.

The scattering potential in equation (2.14) is calculated analogously to the norm integral (3.38) and one obtains

g^lV_cg^k

= 8πaX

i,j

g^lg^jgⁱg^k

= 8πaI₀ A^klij,p^klij, γ^klij

, (3.47)

with

g^klij =g^j∗g^l∗g^kgⁱ, (3.48a) A^klij =A^l∗+A^k+A^j∗ +Aⁱ, (3.48b) p^klij =p^l∗+p^k+p^j∗+pⁱ, (3.48c)

γ^klij =γ^l∗ +γ^k. (3.48d)

The additional integrals

g^lxV_cg^k

, etc. in equation (3.22) can be treated the same way like the overlap integrals (3.41) above by

g^lg^jα^mβⁿgⁱg^k

=I_α^m_βⁿ A^klij,p^klij, γ^klij

. (3.49)

3.3. Dipolar potential

of elliptic integrals, which can be calculated very efficiently. We will see in the following that, with the ansatz (3.1), the integral (2.18) is neither analytically solvable nor reducible to elliptic integrals anymore. Furthermore, we will see in section 3.3.3 that a large number of these integrals has to be calculated in every time step. The numerical evaluation of the dipolar integrals has therefore to be performed as efficiently as possible. In fact, in our simulations the calculation of the dipole integral accounts for more than 80% of the run time, even with all optimizations included.

In this thesis systems are studied, where the dipoles are aligned either in the direction perpendicular or parallel to the direction of the strong confinement along the y-axis. In the following derivations we assume the dipoles to be aligned in a direction perpendicular to the y-direction of the strong confinement, which we choose without loss of generality as z-axis. In appendix A.4 the corresponding derivations are sketched, for the case that the dipoles are aligned parallel to the direction of strong confinement. In the following we use the abbreviations

σ=

1, x, z, x², y², z², xz , (3.50a)

∂_σ^kl= (

1,− ∂

∂p^kl_x ,− ∂

∂p^kl_z , ∂²

(∂p^kl_x)², ∂²

∂p^kl_y 2, ∂²

(∂p^kl_z )², ∂²

∂p^kl_x∂p^kl_z )

(3.50b) κ(σ) =

1, k_x, k_z, k_x², k_y², k²_z, k_xk_z , (3.50c) to accommodate the set of seven integrals required in the vector r of equation (3.22). By the use of the ansatz (3.1) for the condensate wave function Ψ in (2.18) we obtain for the elements of r

g^lσV_dg^k

=X

j,i

g^lσV_d^ijg^k

=X

j,i

ZZ

d³xd³x⁰ σg^kl(x)g^ij(x⁰)

1−3(z−z⁰)²

|x−x⁰|²

1

|x−x⁰|³ . (3.51) Note that in general (3.51) is a six-dimensional integral.

With the abbreviations (3.50) and by the use of the convolution theorem we get g^lσV_dg^k

=X

j,i

Z

d³x σg^kl(x)F⁻¹



F



Z

d³x⁰ g^ij(x⁰)1−3^(z−z0)2

|x−x⁰|²

|x−x⁰|³









=X

j,i

Z

d³x σg^kl(x)F⁻¹



(2π)^3/2F

g^ij(x⁰) F



1−3^(z−z0)2

|x−x⁰|²

|x−x⁰|³







 , (3.52)

where F and F⁻¹ denote the Fourier and inverse Fourier transform, respectively and g^kl defined in equation (3.39a). The terms in (3.52) can be calculated sep-arately. The first Fourier transform of g^ij can be computed analytically which yields

F

g^ij(x⁰) = 1

√2π³ Z∞

−∞

d³x⁰ e⁻(^{x0TAijx0+(pij)Tx0+γij})e^−ikTx0

= e^−γij

√2π³

r π³ detA^ij exp

1

4 p^ij + ikT

A^ij₋1

p^ij + ik

. (3.53)

This can be simplified and expressed by the overlap integral I₀^ij given in equation (3.40) and we obtain

F

g^ij(x⁰) = e^−γij

√2π³

r π³ detA^ij exp

1 4 p^ijT

A^ij₋1

p^ij

×exp

−1

4k^T A^ij−1

k+ 1

2i p^ijT

A^ij−1

k

= I₀^ij

√2π³ exp

−1

4k^T A^ij₋1

k+1

2i p^ijT

A^ij₋1

k

. (3.54)

The Fourier transform of the kernel yields [33]

(2π)^3/2F





1−3^(z−z0)2

|x−x⁰|²

|x−x⁰|³



= 4π 3

3k²_z k² −1

. (3.55)

Insertion of (3.54) and (3.55) in equation (3.52) yields g^lσV_dg^k

= 1 6π²

X

j,i

Z

d³x σg^kl Z

d³k 3k²_z

k² −1

e^ikTx

×I₀^ijexp

−1

4k^T A^ij₋1

k+ 1

2i p^ijT

A^ij₋1

k

. (3.56)

3.3. Dipolar potential

The factor σ in (3.56) can be expressed by derivatives of g^kl with respect to p^kl_σ as already done in the calculation of I_α^m_βⁿ in equation (3.36). Interchanging of the integrations with subsequent integration over xyields

g^lσV_dg^k

= ∂_σ^kl 6π²

X

j,i

Z

d³k I₀^klexp

−1

4k^T A^kl−1

k−1

2i p^klT

A^kl−1

k

×I₀^ijexp

−1

4k^T A^ij₋1

k+ 1

2i p^ijT

A^ij₋1

k

× 3k²_z

k² −1

, (3.57)

where we used (3.50b) and expressed the result in terms of I₀^kl, see (3.40), again.

The two exponentials can be combined by rewriting equation (3.57) as g^lσV_dg^k

= ∂_σ^kl 6π²

X

j,i

I₀^klI₀^ij Z

d³k exp

−1

4k^TA¯^klijk+ 1

2i ¯p^klijT

k

× 3k_z² k²

|{z}

(1)

−1

|{z}

(2)

!

, (3.58)

with the abbreviations

A¯^klij = (A^kl)⁻¹+ (A^ij)⁻¹, (3.59a)

¯

p^klij = (A^ij)⁻¹p^ij −(A^kl)⁻¹p^kl. (3.59b) Equation (3.58) can be split in two terms

g^lσV_dg^k

=X

j,i

∂_σ^klI₀^klI₀^ijh

J_σ=1^klij −J_s^kliji

, (3.60)

whereJ_σ=1^klij ∝term (1) andJ_s^klij ∝term (2) in equation (3.58) (see equation (3.65) below for the definition ofJ_σ^klij). The termJ_s^klij has the same form as the scattering term in equation (3.47) and we can therefore rewrite it as

J_s^klij =− 1 6π²

Z

d³k exp

−1

4k^TA¯^klijk+ i

2 p¯^klijT

k

=−4 3

r 1

πdet ¯A^klij exp

−1

4 p¯^klijT A¯^klij−1

¯ p^klij

=−8π·1 6

vu

ut 1 π³det

(A^kl)⁻¹+ (A^ij)⁻¹

×exp

−1

4 p¯^klijT A¯^klij−1

¯ p^klij

. (3.61)

This leads after some transformations including the replacements of ¯A,p¯by A,p and the use of equations (3.59) to

∂_σ^klI₀^klI₀^ijJ_s^klij =−8π 6

r π³

detA^klij exp 1

4 p^klijT

A^klij−1

p^klij−γ^klij

×P h

A^klij−1

,p^klij i

=−8π

6 I_α^m_βⁿ A^klij,p^klij, γ^klij

, (3.62)

where P[(A^klij)⁻¹,p^klij] is a polynomial of the matrix and vector elements of (A^klij)⁻¹ and p^klij (identical with P in equation (3.37)). In equation (3.62) α^m andβⁿhave to be chosen according toσ in equation (3.60). This term can thus be included by a change of the scattering lengtha →a_eff =a− ¹₆ in equation (3.47) for the computation of the matrix elements of V_c, i.e.

g^lσ(Vc+V_d)g^k

=

g^l σ

a− 1

6 V_c

a +V_d,eff g^k

. (3.63)

Then, the calculation of the dipole integral (3.51) reduces to the calculation of an effective dipole potential

g^lσV_d,effg^k

=X

j,i

∂_σ^klI₀^klI₀^ijJ_σ=1^klij , (3.64)

with

J_σ^klij = 1 2π²

Z

d³k κ(σ)k_z² k² exp

−1

4k^TA¯^klijk+ i

2 p¯^klijT

k

. (3.65)

Here, the seven integralsJ_σ^klij occur after the application of ∂_σ^kl in equation (3.64).

The combination of these integrals in order to obtain (3.64) is described in appendix A.2.1, where the application of the operators ∂_σ^kl is discussed, as well.

Up to here the calculation was completely general. Henceforth, we consider the strong confinement in y-direction and restrict the problem to displacements only in x and z direction (p^klij_y = 0) and to rotations about the y-axis, which leads to matrices of the form

A¯^klij =





A¯^klij_x 0 A¯^klij_xz 0 A¯^klij_y 0 A¯^klij_xz 0 A¯^klij_z



 . (3.66)

3.3. Dipolar potential

Now, a transformation of the 2×2 subsystem inxandz into principle components is possible. This can be achieved by a complex rotation of the complex symmetric matrices ¯A^klij and we obtain

J_σ^klij = 1 2π²

Z

d³k⁰ κ(σ) (k_z⁰ cosα−k_x⁰ sinα)² k⁰² exp

−1

4k^0TA˜¯^klijk⁰+ i

2 p˜¯^klijT

k⁰

= 1 2π²

Z

d³k⁰ κ(σ) (c₀k_z⁰ −c₁k⁰_x)² k⁰² exp

−1

4k^0TA˜¯^klijk⁰+ i

2 p˜¯^klijT

k⁰

, (3.67) where ˜¯A^klij is the diagonalized complex matrix and

A˜¯^klijT

=C^TA¯^klijC , (3.68a)

˜¯

p^klijT

= ¯p^klijT

C , (3.68b)

k=Ck⁰, (3.68c)

with the complex orthonormal transformation matrix C =



 cosα 0 sinα

0 1 0

−sinα 0 cosα



= (c1,c_y,c₂) =





c₀ 0 c₁ 0 1 0

−c₁ 0 c₀



. (3.69)

The complex rotation angle α can be obtained by evaluation of one of the eigen-vectors c₁ or c₂, or equivalently, by a Jacobi transformation [79]. The tilde on A˜¯^klij and the prime on k⁰ will be omitted in the further calculations. After the principle component analysis equation (3.67) takes the form

J_σ^klij = 1 2π²

Z

d³k κ(σ) (c₀k_z−c₁k_x)² k² exp

−1

4 A¯^klij_x k_x²+ ¯A^klij_y k²_y+ ¯A^klij_z k_z²

×exp i

2 p¯^klij_x k_x+ ¯p^klij_z k_z

. (3.70)

Carrying out the k_y-integration yields (c.f. [80] 3.466 1.) forσ 6=y² J_σ^klij = 1

2π² ZZ

dkxdkz 1−Φ r1

4A¯^klij_y (k_x²+k_z²)

!!

×πκ(σ) (c0k_z−c₁k_x)²e^14A¯klijy (^k2x+k_z²) pk²_x+k_z²

×exp

−1

4 A¯^klij_x k_x²+ ¯A^klij_z k_z² + i

2 p¯^klij_x k_x+ ¯p^klij_z k_z

, (3.71)

where Φ(x) denotes the complex error function (see appendix A.3). For σ = y² the k_y-integration yields (c.f. [80] 3.466 2.)

J_y^klij2 = 1 2π²

Z

d³k k_z²k_y² k² exp

−1

4 A¯^klij_x k²_x+ ¯A^klij_y k_y²+ ¯A^klij_z k_z²

×exp i

2 p¯^klij_x k_x+ ¯p^klij_z k_z

= s π³

A¯^klij_y Z

dkxdkz(−c₁k_x+c₀k_z)²e^{−14A¯klijx kx2−14A¯klijz k2z+2ip¯klijx kx+2ip¯klijz kz}

− 1 2π

Z

dkxdkz(−c₁k_x+c₀k_z)²p

k_x²+k_z²erfc 1

2

qA¯^klijy (k_x²+k²_z)

×e^{14A¯klijy (k2x+kz2)}exp

−1

4A¯^klij_x k_x²− 1

4A¯^klij_z k_z²+ i

2p¯^klij_x k_x+ i

2p¯^klij_z kz

, (3.72)

where erfc(x) is the complementary error function (see appendix A.3). The first summand in equation (3.72), in the following denoted by J_y^klij2,0, is a standard Gaussian-type integral which yields

J_y^klij2,0 = −4 ¯A^klij_x p¯_zc₀−A¯^klij_z p¯_xc₁2

+ 8 ¯A^klij_x A¯^klij_z A¯^klij_z c²₁+ ¯A^klij_x c²₀ rA¯^klijx

5

A¯^klijy

A¯^klijz

5

π

×exp (

−1 4

¯ p^klij_x 2

A¯^klijx

+ p¯^klij_z 2

A¯^klijz

!)

, (3.73)

and the second summand in equation (3.72), named with J_y^klij2,1 has the same form as J_σ^klij in equation (3.65) (we therefore set ˜J_y^klij2 ≡ J_y^klij2,1 and omit the tilde from here, however, note that whenJ_y^klij2,1 has been calculated it still has to be added the first summand according toJ_y^klij2 =J_y^klij2,0+J_y^klij2,1). In “polar coordinates”

k_x =ρcosφ , (3.74)

k_z =ρsinφ , (3.75)

3.3. Dipolar potential

we can write equation (3.71) as

J_σ^klij = 1 2π

ZZ

dρdφ erfc

qA¯^klijy

ρ 2

ρ²f_σ(ρ, φ)

×exp

−1

4 −A¯^klij_y + ¯A^klij_x cos²φ+ ¯A^klij_z sin²φ ρ²

×exp

+i

2 p¯^klij_x cosφ+ ¯p^klij_z sinφ ρ

, (3.76)

where f_σ is a function of both integration variables, and

J₁^klij : f₁ = (c0sinφ−c₁cosφ)² , (3.77a) J_x^klij : f_x =ρf₁(c0cosφ+c₁sinφ) , (3.77b) J_z^klij : f_z =ρf₁(c0sinφ−c₁cosφ) , (3.77c)

J_x^klij2 : f_x²=f_x², (3.77d)

J_y^klij2,1 : f_y²=ρ²f₁, (3.77e) J_z^klij2 : f_z²=ρ²f₁², (3.77f)

J_xz^klij : f_xz=f_xf_z. (3.77g)

We rewrite the complementary error function erfc(z) in terms of the Faddeeva func-tion ω(z) = exp(−z²) erfc(−iz) (the exponentially-scaled complementary complex error function, see appendix A.3). The numerical evaluation of w(z) is one of the numerically most expensive computations in the evaluation of the dipole integral.

An efficient algorithm is given e.g. in reference [81]. With the Faddeeva function inserted we obtain

J_σ^klij = 1 2π

Z∞ 0

w

i

qA¯^klij_y ρ 2

ρ²e^{−14A¯klijz ρ2}

× Z2π

0

f_σ(ρ, φ) e⁻¹⁴(^{A¯klijx −A¯klijz} )^ρ2cos2φe²ⁱ(^{p¯klijx ρcosφ+¯pklijz ρsinφ})dφdρ . (3.78)

Equation (3.78) can be transformed to J_σ^klij = 1

2π Z∞

0

dρ w

i

qA¯^klij_y ρ 2

ρ²e^{−14A¯klijz ρ2} X

±sinφ

±cosφ

Zπ/2 0

dφ f_σ(ρ,±sinφ,±cosφ)

×e⁻¹⁴(^A¯klijx −A¯^klijz )^ρ2cos2φe²ⁱ(±p¯^klijx ρcosφ±p¯^klijz ρsinφ), (3.79) where it is summed over all 4 combinations of signs of cosφ and sinφ (note the dependency of f_σ(ρ,±sinφ,±cosφ)). With the substitution x = cos(2φ) = 2 cos²φ−1 we obtain after some short transformations

J_σ^klij = 1 4π

Z∞ 0

dρ Z1

−1

dx X

±xs

±xc

f_σ(ρ,±x_s,±x_c)

√1−x² w(iλρ)ρ²e^{−12A(x)ρ2+p(±xs,±xc)ρ}, (3.80)

with

x_c =x_c(x) =

r1 +x

2 , x_s =x_s(x) =

r1−x

2 , (3.81)

and the abbreviations

λ≡

qA¯^klij_y

2 , (3.82a)

A(x)≡ 1

4 A¯^klij_x + ¯A^klij_z + ¯A^klij_x −A¯^klij_z x

, (3.82b)

p^±(x) = p(±x_s,±x_c)≡ 1

2i ±p¯^klij_x x_c±p¯^klij_z x_s

. (3.82c)

In equation (3.80) it has to be summed over all four combinations of sign of the quantities x_c and x_s. For diagonal matrices A^klij the integral has a simpler form as the factor c₁ = 0. Note that in equation (3.80) all A^klij_σ and all p^klij_σ are the versions after the principle component analysis (3.68).

Up to here all derivations have been performed analytically. However, in general the two-dimensional integral (3.80) cannot be reduced further. In order to solve it we have to employ numerical methods, such as quadratures and expansions. The objective here is to minimize the amount of computational effort as much as possi-ble. In our studies we found a combination of a Taylor expansion of the Faddeeva

3.3. Dipolar potential

function with a Gauß-Chebyshev quadrature to be the appropriate method. The – quite subtle – details of this method will be discussed in section 3.3.1, whereas in appendix A.2.2 a review on numerical quadratures is given and their application to the integral (3.80) is discussed.

Before we address the general case and thereby the numerical solution in section 3.3.1 we notice that some special cases exist, where further simplifications are applicable that allow for avoiding the numerically expensive methods. In those cases that we deal with in section3.3.2, a reduction of the two-dimensional integral (3.80) to elliptic integrals is possible.

3.3.1. Numerical evaluation of J

_σ^klij

In general, the integrals (3.80) have to be calculated numerically. As the r.h.s. vec-torr in equation (3.22) has to be constructed several times for every time step and it contains a large number of such integrals, a very efficient computational method is required. Concerning the x-integration, the application of Gauß-Chebyshev quadrature is self-evident due to the factor (1−x²)^−1/2 that is equivalent to the weight function of this quadrature rule (see appendixA.2.2). We show in appendix A.2.2 that although for the ρ-direction a quadrature is, in principle, possible nu-merical difficulties arise and we find a Taylor expansion of the Faddeeva function to be the fastest and most stable method.

In this sense, we make use of the Taylor expansion of the Faddeeva function at z = 0

w(iz) = X∞ m=0

ξ(m) (2z)^m = 1− 2

√πz+z² − 4 3√

πz³+ 1

2z⁴+. . . , (3.83) with [82]

ξ(m) = −¹₂m

Γ ^m₂ + 1π^{mmod 22 −1−(−1)}

m

4 , m≥0, (3.84)

where Γ(x) is the gamma function [83], and find forw(iλρ) in equation (3.80)

w(iλρ) = X∞ m=0

ξ(m) (2λρ)^m . (3.85)

In terms of the abbreviations given in equations (3.82) we can express equation (3.80) by

J_σ^klij = 1 4π

Z1

−1

dxX

±xs

±xc

Z∞ 0

dρ w(iλρ)ρ²e^{−12A(x)ρ2+p±(x)ρ}f_σ^±(ρ, x)

√1−x²

= 1 4π

Z1

−1

dxX

±xs

±xc

√ 1 1−x²

X∞ m=0

ξ(m) Z∞

0

dρ f_σ^±(ρ, x)ρ^2+me^{−12A(x)ρ2+p±(x)ρ}

| {z }

≡Bσ^±(m,x)

,

(3.86) where we replaced the dependency on (±x_s,±x_c) by the dependency on x and indicate the±by a superscript. The Taylor terms B_σ^±(m, x) can be calculated an-alytically, and we find a recursion formula for the particular terms, see appendix A.2.3. Thereby the starting point of the recursion formula is given in terms of the Faddeeva function again and it is a remarkable result that thus only one (compu-tational expensive) evaluation of w(z) is necessary to gain all Taylor terms. With this result the ρ-integration can be considered given and only the x-integration remains. A Gauß-Chebyshev quadrature is suited best for the numerical compu-tation of this integral, see appendix A.2.2. The actual computation of the Taylor sum does, however, require more effort as described here. We have to consider the application of the Gauß-Chebyshev quadrature and we have to ensure the conver-gence of the Taylor series. This non-trivial task will be handled in the following.

In the practical use of the Taylor expansion (3.86) we sum over a finite number of terms P_∞

m=0 → Pmmax

m=0. To provide for the convergence of the series in all possible configurations of the variational input parameters A^klij and p^klij and to lower the number of m_max two further steps have to be taken. The first one is the application of a nonlinear summation method, the Pad´e approximation, the

‘best’ approximation of a function by a rational function of given order. Details about the Pad´e approximation and the numerical computation with the -Wynn algorithm can be found in appendix A.2.4. The second step that enhances the convergence properties even in cases, where the sum Pmmax

m=0 contains terms with alternating sign, is the interchanging of the Chebyshev quadrature with the Taylor summation

J_σ^klij = 1 4π

X∞ m=0

Z1

−1

dx 1

√1−x²ξ(m)X

±xs

±xc

B^±_σ(m, x). (3.87)

3.3. Dipolar potential

Although two finite sums are interchanged, the step is non-trivial as we apply the Pad´e approximation to the Taylor series, what implies the interchanging of this nonlinear summation method with the Chebyshev sum. To confirm the validity of this step, we applied extensive numerical tests for a wide range of external parameters ¯A^klij and ¯p^klij used in equation (3.80). With the described procedure we obtained best results regarding both the stability of the algorithm and speed of convergence with a usual Taylor order of m_max≈12 and n_quad ≈32 Chebyshev nodes.

3.3.2. Cases with analytic solutions for J

_σ^klij

Equation (3.59b) shows that ¯p^klij = 0 for all integrals with (i, j) = (k, l). In these cases the numerically integrated J_σ^klij integrals reduce to elliptic integrals^†. The following derivations are addressing the integrals with σ = 1 in equation (3.50a), whereas the integrals with σ 6= 1 are discussed at the end of this section. From reference [30] we know the result of the dipolar integral forp^k = 0 andA^klassumed to be diagonal, which reads for (i, j) = (k, l)

g^lV_d^klg^k

= e^−2γkl

r π³ 8 detA^kl

4π 3

hκ˜^kl_x˜κ^kl_yR_D

˜ κ^kl_α2

, κ˜^kl_β2

,1

−1i

, (3.88) where the abbreviations

˜ κ^kl_α2

= A˜¯^klkl_α A˜¯^klkl_β

, α=x, y , β =z , (3.89) are used andR_Ddenotes the elliptic integral of second kind in Carlson form [84,85], see appendixA.2.5

R_D(x, y, z) = 3 2

Z∞ 0

p dt

(x+t)(y+t)(z+t)³ . (3.90) The tilde denotes the diagonal version of ¯A^klkl. The second summand in (3.88) has been dealt with already by introducing the effective dipolar integral in equation (3.64). The difference to the result of reference [30], given in equation (3.88), is made by the presence of the off-diagonal elements of ¯A^klij. It is therefore reasonable to express the quantities J_σ^klkl by elliptic integrals. If we compare the first term in equation (3.88) with the result of equation (3.64)

g^lV_d,eff^kl g^k

= I₀^kl2

J₁^klkl, (3.91)

†This would also be the case, if all GWPs are located at the origin and possess zero momentum, however, in the scenarios dealt with in this thesis this is very unlikely to happen.

we immediately obtain J˜₁^klkl= 4

3

rdetA^kl

8π ˜κ^kl_xκ˜^kl_yR_D

˜ κ^kl_α2

, κ˜^kl_β2

,1

, (3.92)

which is, however, only true for diagonal matrices ¯A^klkl. We may transform the square root in the prefactor as follows

rdetA^kl

8π =

vu utdet

2 ¯A^klkl₋1

8π =

r 1

πdet ¯A^klkl =

s 1

πdet ˜¯A^klkl. (3.93) The integral (3.92) can be split in two terms after the diagonalization of ¯A^klij

J₁^klkl=c²₁J_e^klkl(x) +c²₀J_e^klkl(z), (3.94) with

J_e^klkl(σ) = 4 3

s 1

πdet ˜¯A^klkl˜κ^kl_ακ˜^kl_βR_D

˜ κ^kl_α2

, κ˜^kl_β2

,1

, ˜κ^kl_α,β2

= A˜¯^klkl_α,β A˜¯^klkl_σ , (3.95) where{α, β, σ}is an arbitrary permutation of {x, y, z}. Equation (3.95) simplifies to

J_e^klkl(σ) = 4 3

vu ut 1

π

A˜¯^klkl_σ 3R_D

˜ κ^kl_α2

, ˜κ^kl_β2

,1

. (3.96)

The integrals J_xx^klkl, J_yy^klkl, J_zz^klkl and J_xz^klkl can also be expressed in terms of elliptic integrals. The integrals J_x^klkl and J_z^klij vanish as can be seen easily when taking into account that the antisymmetric integrals yield zero. The calculation of these integrals can be found in appendix A.2.5.

The reduction to elliptic integrals can be done for (i, j) = (k, l) i.e. for the total number of N((N −1)/2 + 1) integrals.

3.3.3. Optimizations

The numerical implementation of all integrals has to be efficient for the reasons given above. We therefore want to avoid unnecessary evaluations of any integral and it is reasonable to apply the following procedure for every computation of the EOM. First, the integrated parameters have to be rewritten in the form of equation (3.31). Then, all integrals I₀^kl have to be calculated and with these the remaining

3.3. Dipolar potential

basic integrals given in section 3.2.1. With these integrals saved, it is now easy to calculate the external, scattering, and dipolar potentials. If all integrals are pre-calculated the matrix K and the vector r in equation (3.15) can be constructed and the EOM can be solved. One further advantage of this procedure is that all quantities of interest, e.g. the mean-field energy, can be calculated easily with this data.

The calculation of the integrals obtained by the Taylor expansion and the Gauß-Chebyshev quadrature – together with the elliptic integrals – are responsible for most of the time cost in the numerical integration. However, the integrals are not independent. Hence, it is appropriate to reduce the number of computed integrals to the minimal number.

In this section we will use the abbreviation “k” for the Gaussian g^k and “l” for g^l∗ etc. If one takes a look at equation (3.51), interchangingk withl and simulta-neously i with j leads to the complex conjugate of the whole integral hg^l|V_d^ij|g^ki. This also holds for any other of the six integrals

g^lσV_d^jig^k

=

g^lσV_d^ijg^k∗

. (3.97)

For theJ_σ^klij-integrals we have to consider that with equation (3.59a) and equation (3.59b)

A¯^lkji = A¯^klij∗

, p¯^lkji= ¯p^klij∗

. (3.98)

Therefore equation (A.7) leads to a different prefactor

∂ⁿ

∂p¯^klij_σ nexp

−1

4k^T A¯^klij∗

k+1

2i ¯p^klij†

k

= ∂ⁿ

∂p¯^klij_σ nexp

−1

4k^T A¯^klij∗

k+ 1

2 p¯^klijT

k

= 1

2k_σ n

exp

−1

4k^TA¯^klijk+1

2i ¯p^klijT

k _∗

. (3.99) Now the J_σ^lkji-integrals can be expressed by J_σ^klij

J₁^lkji= J₁^klij

_∗

, J_x^lkji =− J_x^klij∗

, J_z^lkji =− J_z^klij∗

,

J_xx^lkji = J_xx^klij∗

, J_yy^lkji = J_yy^klij∗

,

J_zz^lkji = J_zz^klij∗

, J_xz^lkji = J_xz^klij∗

. (3.100) Changing the index pairs (k, l) and (i, j) yields a minus sign for ¯p^klij in equation (3.59b). This leads to a minus sign for J_x^ijkl and J_z^ijkl.

J₁^ijkl=J₁^klij, J_x^ijkl =−J_x^klij, J_z^ijkl =−J_z^klij,

J_xx^ijkl =J_xx^klij, J_yy^ijkl =J_yy^klij,

J_zz^ijkl =J_zz^klij, J_xz^ijkl =J_xz^klij.

(3.101)

Applying both index swaps yields J₁^jilk =

J₁^klij∗

, J_x^jilk = J_x^klij∗

, J_z^jilk = J_z^klij∗

,

J_xx^jilk = J_xx^klij∗

, J_yy^jilk = J_yy^klij∗

,

J_zz^jilk= J_zz^klij∗

, J_xz^jilk= J_xz^klij∗

. (3.102) We notice that for a given set (k, l, i, j) the integrals for the combinations (l, k, j, i), (i, j, k, l) and (j, i, l, k) are given. Therefore a total number ofC_numerical = (N⁴+ N²−2N)/4 integrals have to be calculated numerically and C_elliptic =N(N+ 1)/2 can expressed in terms of elliptic integrals. The remaining N⁴ integrals can be obtained from these by the relations given above.

If the initial wave is an eigenfunction of the parity operator, the number of Gaussians can be reduced. In this case only half of the number of integrals have to be calculated. The vector r = (r₀,r_x,r_z,r_xx,r_yy,r_zz,r_xz)^T in equation (3.22) is then given by

r₀^l = XN k=1

g^lV(x)g^k

= XN/2 k=1

g^lV(x)g^k +

g^lV(x)Πg^k

, (3.103)

r₀^l+N/2 = XN k=1

Πg^lV(x)g^k

= XN/2 k=1

Πg^lV(x)g^k +

Πg^lV(x)Πg^k

=r^l₀ , (3.104) and as can be shown easily

r^l+N/2_σ =s·r_σ^l , (3.105)

wheres =−1 for σ= (x, z) ands = 1 otherwise.

Im Dokument Variational approaches to dipolar Bose-Einstein condensates (Seite 42-56)