Explicit expressions for dipolar integrals in the TDVP

A. Supplementary material for LSE construction 115

A.2. Supplementary material for the dipolar integrals

A.2.1. Explicit expressions for dipolar integrals in the TDVP

In the vector r of equation (3.15) the 7N² integrals hg^l|σV_d|g^ki with the abbre-viation σ = (1, x, z, x², y², z², xz) given in equation (3.50a) appear. The terms with σ 6= 1 can be obtained by derivation of hg^l|V_d|g^ki. The derivatives may be rewritten as

∂

∂p^kl_x =−

A^kl−1

∂

∂p¯^klijx −

A^kl−1

∂

∂p¯^klijz

∂

∂p^kl_z =−

A^kl−1

∂

∂p¯^klijz −

A^kl−1

∂

∂p¯^klijx

,etc. (A.4) We use the following abbreviations

−

A^kl−1

αβ = ¯G^kl_αβ ≡G¯_αβ , G¯_σσ ≡G¯_σ ,

∂

∂p¯^klijσ

≡∂_p_¯_σ , with α, β, σ ∈ {x, y, z}, α 6=β , (A.5)

A.2. Supplementary material for the dipolar integrals

whereα and β denote the elements of the inverse matrix of−A^kl. The derivatives then read

dp^kl_x = ¯G_x∂_p_¯_x + ¯G_xz∂_p_¯_z, (A.6a) d

dp^kl_y = ¯G_y∂_p_¯_y ≡0, (A.6b)

dp^kl_z = ¯G_z∂_p_¯_z + ¯G_xz∂_p_¯_x, (A.6c) d²

(dp^kl_x)² = ¯G²_x∂_p²_¯_x + 2 ¯G_xG¯_xz∂_p_¯_z∂_p_¯_x + ¯G²_xz∂_p²_¯_z, (A.6d) d²

dp^kl_y 2 = ¯G²_y∂_p²_¯_y, (A.6e)

d²

(dp^kl_z )² = ¯G²_z∂_p²_¯

z + 2 ¯G_zG¯_xz∂_p_¯_x∂_p_¯_z + ¯G²_xz∂_p²_¯

x, (A.6f)

d²

dp^kl_xdp^kl_z = ¯G_xG¯_xz∂_p²_¯_x + ¯G_zG¯_xz∂_p²_¯_z + ¯G_xG¯_z+G²_xz

∂_p_¯_x∂_p_¯_z . (A.6g)

The derivation ofJ_σ^klij leads to integrals of the form (3.65) which can be seen from

∂ⁿ

∂p¯^klijσ

nexp

−1

4k^TA¯^klijk+ 1

2i ¯p^klijT

= i

2k_σ n

exp

−1

4k^TA¯^klijk+1

2i ¯p^klijT

. (A.7)

In y-direction all integrals with odd powers of y vanish because p^klij_y = 0. Com-bining (A.6a), (3.65), and (A.7) and abbreviating ₂ⁱG¯^kl_σ ≡ G^kl_σ yields the dipolar integrals in the TDVP

g^lV_d,effg^k

i,j

I₀^klI₀^ijJ₁^klij, (A.8a)

g^lxV_d,effg^k

=−X

i,j

I₀^ij

∂_p^kl

xI₀^kl

J₁^klij+I₀^kl G^kl_xJ_x^klij +G^kl_xzJ_z^klij

, (A.8b) g^lzV_d,effg^k

=−X

i,j

I₀^ij

∂_pkl z I₀^kl

J₁^klij+I₀^kl G^kl_zJ_z^klij +G^kl_xzJ_x^klij

, (A.8c) g^lx²V_d,effg^k

i,j

I₀^ijh

∂_p²kl x I₀^kl

J₁^klij + 2 ∂_p^kl

z I₀^kl

G^kl_xJ_x^klij+G^kl_xzJ_z^klij +I₀^kl

G^kl_x2

J_xx^klij+ 2G^kl_xG^kl_xzJ_xz^klij + G^kl_xz2

J_zz^klij i , (A.8d) g^ly²V_d,effg^k

i,j

I₀^ij

∂_p²kl yI₀^kl

J₁^klij+I₀^kl G^kl_y2

J_yy^klij

, (A.8e)

g^lz²V_d,effg^k

i,j

I₀^ijh

∂_p²kl z I₀^kl

J₁^klij + 2 ∂_pkl x I₀^kl

G^kl_zJ_z^klij+G^kl_xzJ_x^klij +I₀^kl

G^kl_x2

J_xx^klij+ 2G^kl_xG^kl_xzJ_xz^klij + G^kl_xz2

J_zz^klij i

, (A.8f) g^lxzV_d,effg^k

i,j

I₀^ijh

∂_p^kl

x∂_p^kl

z I₀^kl

J₁^klij + ∂_p^kl

xI₀^kl

G^kl_zJ_z^klij +G^kl_xzJ_x^klij + ∂_p^kl_z I₀^kl

G^kl_xJ_x^klij +G^kl_xzJ_z^klij +I₀^kl

G^kl_xz G^kl_xJ_xx^klij +G^kl_zJ_zz^klij +

G^kl_xG^kl_z + G^kl_xz2

J_xz^kliji

. (A.8g)

A.2.2. Quadratures

The basic idea of numerical quadrature rules is the approximate computation of an integral Rb

a f(x)dx by interpolation of the integrand f(x). Quadratures rely on the evaluation of f(x) at a finite set of points that are called the abscissas or nodes and processing these values to produce an approximation of the integral (usually involving a weighted average). Widely-used quadratures are those relying on orthogonal polynomials (e.g. Gauß-Hermite and Gauß-Chebyshev).

A.2. Supplementary material for the dipolar integrals

Figure A.1.:Real part of the integrandK_J¹¹¹²

1 in equation(3.80)as a two-dimensional function of the integration coordinates ρand x. The variational parameters have been chosen such that there is a large distance between the Gaussians g¹ and g². This results in a highly oscillating integrand.

ρ-integration

As pointed out by P. J¨ackel in reference [122] it is not possible to use a Gauß-Laguerre quadrature for the ρ-integration (“Hold your horses right there!”[122]) by substituting ν =ρ²/2, because this introduces a square root function that has no polynomial representation. One should rather use a half-open Gauß-Hermite quadrature with the weight function ρe⁻¹²^ρ² for the integration of equation (3.80).

The corresponding roots and weights of this quadrature can be obtained by a procedure given e.g. in reference [79, 123] and are also given in reference [124].

In order to apply this quadrature, it is convenient to rewrite the ρ-integral in (3.80) into a form corresponding to [122]

Z∞ 0

f(ρ)e⁻¹²^ρ²ρdρ≈ Xm

j=1

ω_jf(ρj). (A.9)

Although this quadrature can be applied successfully in many cases, difficulties arise, particularly in regions of the variational parameter space when two GWPs have large spatial separation or possess large relative momenta. In figure A.1 the integrand of equation (3.80) is shown in such a case and we immediately see that here the number of nodes would have to be very high. Furthermore,

the integrand shown in figure A.1 is not restricted to the area plotted there and sometimes oscillations are present even for high values ofρ, with additional crucial dependency on x.

In regions of the parameter space, where the dipolar integral is not zero, but the distance of the Gaussians is still large, a multipole expansion of the DDI can be done. Yet, we found slow convergence with increasing order of the Taylor expansion. For these reasons it is more appropriate to apply a second method, the Taylor expansion of the integrand which we discussed in section3.3.1.

x-integration

The x-integration can be carried out with a Gauß-Chebyshev quadrature Z+1

−1

f(x)

√1−x²dx≈

nXquad

i=1

w_if(xi), (A.10)

as the weight function is (1−x²)^−1/2 and the roots x_i are given by

x_i = cos

2i−1 2nquad

→









 x_c =

r1 +x_i 2 = cos

2i−1 4nquad

x_s =

r1−x_i 2 = sin

2i−1 4nquad

(A.11)

wherei= 1, . . . , nquad and n_quad is the number of nodes. The weights only depend on the number of nodes

ω_i = π

n_quad . (A.12)

The combination of this integration with the Taylor expansion is discussed in section 3.3.1. It turns out that in most cases a number of n_quad ≈ 32 nodes is sufficient.

A.2.3. Recursion formula for the terms of the Taylor series

Our aim here is the computation of the terms B_σ^± in equation (3.86). We will now consider only the case σ = 1 → B_σ=1^± (see equation (3.50a)). All other cases can be handled in the same manner because the additional factor of ρ or ρ² in f_σ^± only leads to an index shift in the Taylor series. In the following we omit all dependencies of B, A, and p on x, the combination (k, l, i, j), the index σ on B,

A.2. Supplementary material for the dipolar integrals

and furthermore we omit the index ± at B^± and p^±. For the Taylor terms B(m) given in equation (3.86) we find

B(m) = Z∞ 0

dρ ρ^2+me⁻¹²^Aρ²^+pρ

= r2

m+3Z∞ 0

dy y^2+me⁻^y²⁻^2by², (A.13)

where b=−¹₂pp

2/A. By partial integration we find

B(m) = r2

m+3

e^b² Z∞

(y−b)^m+2e⁻^y²dy (A.14)

= r2

m+3

e^b²

m+ 3(y−b)^m+3e⁻^y² _∞

− Z∞

−2y

m+ 3 (y−b)^m+3e^−y²dy

! .

For the first term the upper and lower limit vanishes, the second term yields a prefactor and we obtain

B(m) = r2

m+3 2

m+ 3e^b² Z∞

y(y−b)^m+3e^−y²dy . (A.15) The integral in equation (A.15) can be transformed to

Z∞ b

y(y−b)^m+3e⁻^y²dy= Z∞

((y−b) +b) (y−b)^m+3e⁻^y²dy

= Z∞

(y−b)^m+4e⁻^y²dy+ Z∞

b(y−b)^m+3e⁻^y²dy , (A.16)

and we obtain

B(m) = r2

m+3 2

m+ 3[Im+4+b I_m+3] , (A.17)

where the integralsI_m are given by the recursion relation I₀ =

√π

2 w(ib) , (A.18a)

I₁ = 1

2−b I₀, (A.18b)

I_m = m−1

2 I_m−2−b I_m−1. (A.18c)

Proof 1 Recursion formula for I_m

I₀ = e^b² Z∞

e^−y²dy =

√π

2 ω(ib) , (A.19)

I₁ = e^b² Z∞

(y−b)e^−y²dy= e^b² Z∞

ye^−y²dy−bI₀ = 1

2e^b²e^−b² −bI₀ = 1

2 −bI₀, (A.20)

I_m ≡e^b² Z∞

(y−b)^me⁻^y²dy= e^b² Z∞

(y−b)^me⁻^y²dy

z }|0 {

e^b²

(y−b)^m+1 m+ 1 e⁻^y²

#_∞

−e^b² Z∞

−2y

m+ 1(y−b)^m+1e⁻^y²dy

= 2

m+ 1



e^b² Z∞

(y−b)^m+2e^−y²dy+be^b² Z∞

(y−b)^m+1e^−y²dy





= 2

m+ 1[I_m+2+bI_m+1] (A.21)

Therefore, we find

I_m+4 = m+ 3

2 I_m+2−b I_m+3 and ⇒I_m+4+b I_m+3 = m+ 3

2 I_m+2. (A.22)

We insert these relations in equation (A.17) and yield B(m) =

r2 A

m+3

I_m+2. (A.23)

A.2. Supplementary material for the dipolar integrals

A.2.4. Pad´ e approximation and -Wynn algorithm

For a given function f(x) the Pad´e approximant P_ij is defined by the quotient of two polynomials

P_ij = Pi s=0

p_ijsx^s Pj s=0

p⁰_ijsx^s

, i= 0,1, . . . , j= 0,1, . . . , (A.24)

that agrees with f(x) to the highest possible order of the Taylor expansion^∗ of f(x) atx= 0. The relation between the coefficients p_ijs andp⁰_ijs by this condition is given by the corresponding Pad´e table [125].

For given x the Pad´e approximants can be computed from the partial sums of the Taylor series by Wynn’s -algorithm [125] that makes use of relationships in the Pad´e table. The functionf(x) may also be a formal power series and the Pad´e approximation can therefore also be applied to the summation of divergent series.

A.2.5. Additional dipolar integrals reducible to elliptic integrals

Here the integralsJ_xx^klkl, J_yy^klkl, J_zz^klklandJ_xz^klkl, which can be expressed in terms of el-liptic integrals, will be calculated. At first we will discuss the integrals examplarily for the integral

J_exx^klkl = 1 2π²

d³kk²_xk²_z

|k|² exp

−1

4k^TA¯^klklk

= 1 2π²

d³k(c⁴₁+c⁴₀−4c²₁c²₀)k_x²k_z²+c²₀c²₁k_z⁴+c²₀c²₁k⁴_x

|k|² exp

−1

4k^TA˜¯^klklk

. (A.25) The terms with odd powers in eitherk_x ork_z have already been omitted since the integrals with an antisymmetric integrand vanish. The terms in equation (A.25) can be expressed by derivatives with respect to ˜¯A^klkl_σ

1 2π²

d³kk_α²k_β²

|k|² exp

−1

4k^TA˜¯^klklk

=−4 ∂

∂A˜¯^klkl_β

J_e^klkl(α). (A.26)

∗Actually, the Pad´e approximation is defined for any power series, yet our application deals with Taylor series.

The derivatives of equation (3.96) read

∂

∂A˜¯^klkl_σ

J_e^klkl(σ) = −2 3

vu ut 1

A˜¯^klkl_σ 7

3 ˜¯A^klkl_σ R_D

˜ κ^kl_α2

, ˜κ^kl_β2

,1 + 2 ˜¯A^klkl_α R_∂1

˜ κ^kl_α2

, ˜κ^kl_β2

,1 + 2 ˜¯A^klkl_β R_∂2

˜ κ^kl_α2

, ˜κ^kl_β2

,1# ,

(A.27a)

∂

∂A˜¯^klkl_α

J_e^klkl(σ) = 4 3

vu ut 1

π A˜¯^klkl_σ

5 R_∂1

˜ κ^kl_α2

, κ˜^kl_β2

, (A.27b)

with the derivatives of the elliptic integrals equation (A.25) then reads

J_xx^klkl=−4 c⁴₁+c⁴₀−4c²₁c²₀ ∂

∂A˜¯^klkl_x

J_e^klkl(z)

−4c²₀c²₁ ∂

∂A˜¯^klkl_z

J_e^klkl(z) + ∂

∂A˜¯^klkl_x

J_e^klkl(x)

. (A.28)

Thereby, the elliptic integral of second kind in Carlson form is given by

R_D(x, y, z) = 3 2

Z∞ 0

p dt

(x+t)(y+t)(z+t)³ , (A.29)

and the derivatives R_∂1 and R_∂2 with respect to the first and second argument, respectively, read [30]

R_∂1(x, y, z)≡ ∂

∂xR_D(x, y, z) = 1

2(z−x)[RD(x, y, z)−R_D(z, y, x)] , (A.30a) R_∂2(x, y, z)≡ ∂

∂yR_D(x, y, z) = 1

2(z−y)[RD(x, y, z)−R_D(x, z, y)] . (A.30b)

A.2. Supplementary material for the dipolar integrals

Analogous, the calculation for J_yy^klkl, J_zz^klkl and J_xz^klkl can be done

J_yy^klkl = 1 2π²

d³kk²_yk²_z

|k|² exp

−1

4k^TA¯^klklk

= 1 2π²

d³kc²₁k_x²k²_y+c²₀k²_zk_y²

|k|² exp

−1

4k^TA˜¯^klklk

, (A.31)

J_zz^klkl = 1 2π²

d³k k⁴_z

|k|² exp

−1

4k^TA¯^klklk

= 1 2π²

d³kc⁴₁k_x⁴+ 6c²₁c²₀k_x²k_z²+c⁴₀k_z⁴

|k|² exp

−1

4k^TA˜¯^klklk

, (A.32)

J_xz^klkl = 1 2π²

d³kk_xk³_z

|k|² exp

−1

4k^TA¯^klklk

= 1 2π²

d³kc³₀c₁k_z⁴−c³₁c₀k_x⁴+ 3 (c0c³₁−c³₀c₁)k²_xk²_z

|k|² exp

−1

4k^TA˜¯^klklk

. (A.33)

With the abbreviation ∂/

∂A˜¯^klkl_σ

≡∂¯_σ we can summarize

J₁^klkl =c²₁J_e^klkl(x) +c²₀J_e^klkl(z), (A.34a) J_xx^klkl =−4 c⁴₁+c⁴₀−4c²₁c²₀∂¯_xJ_e^klkl(z)−4c²₀c²₁ ∂¯_zJ_e^klkl(z) + ¯∂_xJ_e^klkl(x)

, (A.34b) J_yy^klkl =−4c¹₁∂¯_yJ_e^klkl(x)−4c²₀∂¯_yJ_e^klkl(z), (A.34c) J_zz^klkl =−4c⁴₁∂¯_xJ_e^klkl(x)−24c²₀c²₁∂¯_xJ_e^klkl(z)−4c⁴₀∂¯_zJ_e^klkl(z), (A.34d) J_xz^klkl = 4c0c³₁∂¯_xJ_e^klkl(x)−12 c₀c³₁−c³₀c₁∂¯_xJ_e^klkl(z)−4c³₀c₁∂¯_zJ_e^klkl(z). (A.34e)

Therefore only six elliptic integrals have to be calculated

R_D

A˜¯^klkl_x A˜¯^klkl_z ,

A˜¯^klkl_y A˜¯^klkl_z ,1

, R_D

A˜¯^klkl_x A˜¯^klkl_z ,1,

A˜¯^klkl_y A˜¯^klkl_z

, R_D 1, A˜¯^klkl_y A˜¯^klkl_z ,

A˜¯^klkl_x A˜¯^klkl_z

! ,

R_D A˜¯^klkl_z A˜¯^klkl_x ,

A˜¯^klkl_y A˜¯^klkl_x

, R_D A˜¯^klkl_z A˜¯^klkl_x

,1,A˜¯^klkl_y A˜¯^klkl_x

, R_D 1,A˜¯^klkl_y A˜¯^klkl_x

,A˜¯^klkl_z A˜¯^klkl_x

! . (A.35)

A.3. Error function and related functions

In the derivation of an appropriate form of the dipolar integral we find the complex error function

Φ(z) = 2

√π Zz

e^−t²dt , (A.36)

where the argumentzis a complex number and, therefore, Φ(z) is a complex-valued function. This function is connected to the complementary error function

erfc(z) = 1−Φ(z), (A.37)

and the Faddeeva function

w(z) = e⁻^z²erfc (−iz) , (A.38a)

⇔erfc(z) = e⁻^z²w(iz) , d

dzw(z) = 2i

√π −2zw(z) , (A.38b)

which is used for the numerical calculation of all forms of the complex error func-tion, as it is the most stable and efficient form [81, 126]. The Faddeeva function is visualized in figure A.2.

(a) (b)

Figure A.2.: (a) Plot of the Faddeeva function w(z), defined in equation (A.38a) in the complex plane. The real part is represented by the brightness and the phase is given by the hue, where in (b) z is plotted to show the color-scale.

A.4. Dipoles parallel to strong confinement

In section 3.3 we performed all calculations of the dipole integral for the dipole axis being perpendicular to the strong confinement. Here we briefly sketch the changes that need to be applied, when the dipoles are aligned parallel to the strong confinement. Therefore, we choose the y-axis as the dipole-axis. Then the dipolar integral reads

D g^l

σVˆ_d

g^kE

j,i

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

1− 3(y−y⁰)²

|x−x⁰|²

|x−x⁰|³ , (A.39) with g^kl defined in equation (3.39a) and σ given in equation (3.50a). By a calcu-lation analogous to the one in section 3.3 we find

Jˆ_σ^klij = 1 2π²

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

−1

4k^TA¯^klijk+1

2i ¯p^klijT

, (A.40)

with κ(σ) given in equation (3.50c). Both integrals J_s^klij and the prefactors in hg^l|σVˆ_d|g^ki do not change. The integral (A.40) can be rewritten as

Jˆ_σ^klij =− 4 2π²

∂

∂A¯^klij_y Z

d³k κ(σ) k² exp

−1

4k^TA¯^klijk+1

2i ¯p^klijT

, (A.41) and we obtain analogously to the calculation in section 3.3

J_σ^klij =−1 π

∂

∂A¯^klijy

Z∞ 0

dρ Z1

−1

dx X

±xs

±xc

fˆ_σ^±(ρ, x)

√1−x²w(iλρ) e⁻¹²^A(x)ρ²^+p^±^(x)ρ, (A.42)

with ˆf_σ^± = f_σ^±/f₁^± and the abbreviations (3.82). Note that the superscript ± indicates dependencies on ±x_c and ±x_s, again and that λ=λ( ¯A^klij_y ).

A.4.1. Taylor expansion of the integrand for polarization parallel to the strong confinement

For the derivation regarding the perpendicular orientation of the dipoles see section A.2.3. We start once again with the Taylor expansion of the Faddeeva function and apply the operator ∂/∂A¯^klij_y

∂

∂A¯^klijy

w iλ A¯^klij_y ρ

= X∞ m=0

ξ(m) (2ρ)^m m

2 A¯^klij_y ^m₂−1

= X∞ m=0

ξ(m)m

2 (2ρ)^m(2λ)^m⁻² , (A.43)

where we have used the abbreviation (3.84) and we obtain Jˆ_σ^klij =−1

π Z1

−1

dxX

±xs

±xc

Z∞ 0

dρ w(iλρ) e⁻¹²^A(x)ρ²^+p^±^(x)ρfˆ_σ^±(ρ, x)

√1−x²

=−1 π

−1

dxX

±xs

±xc

√ 1 1−x²

X∞ m=1

ξ(m)ˆ Z∞

dρfˆ_σ^±ρ^me⁻¹²^A(x)ρ²^+p^±^(x)ρ

| {z }

≡Bˆ^±(m,x)

. (A.44)

The terms of the expansion can be obtained analogously (omitting indices ± and dependencies of x):

B(m) =ˆ Z∞ 0

dρ ρ^me⁻¹²^Aρ²^+pρ

= r2

m+1Z∞ 0

dy y^me⁻^y²⁻^2by²

= r2

m+1 2

m+ 3e^b²



 Z∞

(y−b)^m+2e⁻^y²dy+ Z∞

b(y−b)^m+1e⁻^y²dy





= r2

m+1 2

m+ 1[Im+2+b I_m+1]

= r2

m+1

I_m, (A.45)

whereI_m is given by the recursion relation (A.18).

A.4.2. Elliptic integrals for polarization axis parallel to the strong confinement

The dipolar integrals given in section 3.3.2, which could be reduced to elliptic integrals, are given for the dipole polarization axis y (parallel to the strong con-finement) as follows

D g^l

Vˆ_d^kl

g^kE

= e⁻^2γ^kl

r π 8 detA^kl

4π 3

ˆκ^kl_xˆκ^kl_zR_D κˆ^kl_x,κˆ^kl_y ,1

−1

, (A.46)

A.4. Dipoles parallel to strong confinement

where the abbreviations now read ˆ

κ^kl_α2

= Aˆ¯^klkl_α Aˆ¯^klkl_β

, α=x, z , β =y . (A.47) As no diagonalization is necessary in this case, we obtain

Jˆ₁^klkl= ˆJ_e^klkl(y), (A.48) where

Jˆ_e^klkl(σ) = 4 3

s 1

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

, κˆ^kl_β2

, κˆ^kl_α,β2

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

(A.49) with α, β, σ =x, y, z , α6=β 6=σ ,

and

Jˆ_exx^klkl= 1 2π²

d³kk_x²k_y²

|k|² exp

−1

4k^TA¯^klklk

= 1

2π² Z

d³kc²₀k_x²k_y²+c²₁k_z²k_y²

|k|² exp

−1

4k^TA˜¯^klklk

, (A.50)

where in analogy to (A.25) the terms with odd powers of eitherk_x ork_z have been omitted. We follow (A.26) and obtain

Jˆ₁^klkl= ˆJ_e^klkl(y), (A.51)

Jˆ_xx^klkl=−4c²₀∂¯_xJˆ_e^klkl(y)−4c²₁∂¯_zJˆ_e^klkl(y), (A.52) Jˆ_yy^klkl=−4 ¯∂_yJˆ_e^klkl(y), (A.53) Jˆ_zz^klkl=−4c²₁∂¯_xJˆ_e^klkl(y)−4c²₀∂¯_zJˆ_e^klkl(y), (A.54) Jˆ_xz^klkl= 4c0c₁∂¯_xJˆ_e^klkl(y)−4c0c₁∂¯_zJˆ_e^klkl(y), (A.55) where the derivatives are given by (A.27). Thus, in the case of the polarization axis chosen to be the y-axis even less integrals, namely

R_D A˜¯^klkl_x A˜¯^klkl_y

,A˜¯^klkl_z A˜¯^klkl_y

, R_D A˜¯^klkl_x A˜¯^klkl_y

,1,A˜¯^klkl_z A˜¯^klkl_y

, R_D 1,A˜¯^klkl_z A˜¯^klkl_y

,A˜¯^klkl_x A˜¯^klkl_y

! , (A.56) have to be calculated.

A.5. Gaussian integrals over subspace

Here, integrals of Gaussian type over the half volume are calculated. We will use the following notation:

I_σ^klm+ ≡ Z∞

dσ Z∞

∞

dα Z∞

∞

dβ σ^mg^l^∗g^k, I_σ^klm− ≡ Z0

−∞

dσ Z∞

∞

dα Z∞

∞

dβ σ^mg^l^∗g^k, (A.57) where σ, α, β=x, y, z, σ 6=α 6=β 6=σ, and m = 0,1,2. We are only interested in integrals whereσ =x, z. Furthermore, the simplifications of the restricted ansatz (3.17) reduce the problem so that the y-integration

I_(y)^kl = r π

A^kl_y (A.58)

can always be pulled outside. The remaining two-dimensional integration yields for m= 0

I_σ^kl0± =I_(y)^kl × πe (^pkl^α)²

4Aklαα −γ^kl

2 q

A^kl_xxA^kl_zz −(A^kl_xz)² w





± i A^kl_ααp^kl_σ −A^kl_xzp^kl_α 2

r A^kl_αα

A^kl_xxA^kl_zz −(A^kl_xz)²





 , (A.59)

where σ, α = x, z, σ 6= α, and w(x) is the Faddeeva function defined in equa-tion (A.38a). The factor σ^m can be expressed by negative derivatives −∂/∂p^kl_σ of I_σ^klij0± once again. For the derivative of w(x) see equation (A.38b). If we use the abbreviations

x_± =± i A^kl_ααp^kl_σ −A^kl_xzp^kl_α 2

r A^kl_αα

A^kl_xxA^kl_zz −(A^kl_xz)², (A.60a)

∂x_±

∂p^kl_σ ≡x⁰_± =± iA^kl_αα 2

r A^kl_αα

A^kl_xxA^kl_zz −(A^kl_xz)², (A.60b)

A.6. Excerpt from the source code

for the argument of the Faddeeva function and its derivative with respect to p^kl_σ, respectively, we obtain

I_σ^kl1± = 2x⁰_±

x_±w(x_±)− i

√π

×I_(y)^kl πe (^pkl^α)²

4Aklαα −γ^kl

2 q

A^kl_xxA^kl_zz −(A^kl_xz)²

, (A.61)

I_σ^kl2± = 2 x⁰_±2

2 (x_±)²w(x_±)− 2ix_±

√π −w(x_±)

×I_(y)^kl πe (^pkl^α)²

4Aklαα −γ^kl

2 q

A^kl_xxA^kl_zz −(A^kl_xz)² , (A.62) for m= 1 and m= 2, respectively.

A.6. Excerpt from the source code

As discussed in section3.5, theFortran-code alone makes up for more than 10 000 lines of source code that is available from the author on request. In the following we present an excerpt from the source code for the calculation of the integralJ_σ^klij in equation (3.87) that represents one of the most crucial parts in the computation of the EOM (3.14). The following code has to be taken rather as an exemplary code snippet than as a complete part of a program, as the modulesglobaland Er-rorFunction are not given here. Furthermore, the -Wynn algorithm (subroutine EPSAL) is given in reference [125].

Subroutine for computation of J_σ^klij with dipoles in z-direction

modulemyKinds implicit none

integer,parameter ::sp=kind(1.0) !# single precision integer,parameter ::dp=kind(1.0d0) !# double precision end modulemyKinds

typeklij_combination_type

integer ::k,l,i,j !# index combination if combined type character::integral=’n’ !# kind of integral (n,e)

end typeklij_combination_type typevar

complex(dp),dimension(4)::A=1._dp !# width parameters of gaussian complex(dp),dimension(3)::p=1._dp !# complex p (translation and impuls) complex(dp)::g=1._dp !# amplitude and phase of gaussian type(klij_combination_type)::combination !# index combination if combined type end typevar

subroutineJ2_numerical_z(zklij,c,TOL,J2,convergence_order)

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

!%% Main Jsigma-calculation routine for dipoles in z-direction. Here the work is done. %%

!%% --- %%

!%% Variables %%

!%% INPUT: %%

!%% zklij : type(var), holds the variational parameters %%

!%% c : complex, entries of the C-matrix used for Jacobi-transformation %%

!%% TOL : tolerance for chebyshev-quadrature, if not constant order %%

!%% convergence_order : maximum oder of chebyshev integration %%

!%% OUTPUT: %%

!%% J2 : complex, dimension 7, hold the desired Jsigma-integrals on output %%

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

useglobal !# uses module myKinds

useErrorFunction !# contains the subroutine WOFZ (Faddeva function) implicit none

type(var),intent(in)::zklij !#(A^kl)⁻¹ + (A^ij)⁻¹ complex(dp),dimension(0:1),intent(in)::c !# entries of rotation matrix real(dp),intent(in)::TOL !# relative error

complex(dp),dimension(7),intent(out)::J2 !# result

integer,intent(inout)::convergence_order !# maximum order of chebyshev nodes used integer,dimension(7)::addPow=(/0,1,1,2,2,2,2/)!# add. power ofJ_σ int. (J_x:1,J_xx:2, ...) complex(dp),dimension(7)::fsigma !# terms dependent ofxc,xsfor all comb. of sign complex(dp)::AyS,Am,Am2S,pm !# abbreviations; if last letter=S: sqrt involved complex(dp),dimension(0:potential%dipole%taylorOrder+4)::AyS_pow !# powers of AyS

complex(dp)::b !# abbreviation for completing the square

complex(dp),dimension(0:potential%dipole%taylorOrder,7)::D,D2,J2pade !# temporary array

complex(dp),dimension(0:potential%dipole%taylorOrder+4)::Ik,Ik2 !# powers of b, Ik given by recursion complex(dp)::erf_arg !# error function argument, abreviations (see below) real(dp)::U,V !# real and imaginary part of error function result real(dp)::x,xc,xs !# cosroot,sinroot for all combinations of sign real(dp),dimension(2)::signSC=(/1._dp,-1._dp/) !# for sum over all combinations of sign integer::i,j,k,l,m,jj,tn,o,jjo !# loop integer;

logical::FLAG !# flag for Faddeeva function

integer,dimension(7)::PFLAG !# flag for pade approximation integer,dimension(3)::corder !# chebyshev sum

complex(dp)::J20yy !# intermediate results

complex(dp)::estlim !# result of pade approximation

complex(dp),dimension(0:potential%dipole%taylorOrder,7)::est !# result of pade approximation

complex(dp),dimension(potential%dipole%taylorOrder,7)::sumD !# sum of taylor terms multip. with Xi and weight real(dp),dimension(potential%dipole%taylorOrder,7)::diffD !# diff. betw. chebyshev results for diff. orders complex(dp),dimension(0:potential%dipole%taylorOrder)::Xi !# coefficients for taylor expansion in AyS*rho real(dp),dimension(0:30)::Xic=(/1._dp,-0.564189583547756287_dp,0.25_dp,-0.0940315972579593812_dp,0.03125_dp,&

&-0.00940315972579593812_dp,0.00260416666666666667_dp,-0.000671654266128281294_dp,0.000162760416666666667_dp,&

&-0.0000373141258960156274_dp,8.13802083333333333d-6,-1.69609663163707397d-6,3.39084201388888889d-7,&

&-6.52344858321951529d-8,1.21101500496031746d-8,-2.17448286107317176d-9,3.78442189050099206d-10,&

&-6.3955378266857993d-11,1.05122830291694224d-11,-1.6830362701804735d-12,2.6280707572923556d-13,&

&-4.0072292147154131d-14,5.97288808475535364d-15,-8.71136785807698499d-16,1.24435168432403201d-16,&

&-1.742273571615397d-17,2.39298400831544617d-18,-3.22643254002851296d-19,4.2731857291347253d-20,&

&-5.56281472418709131d-21,7.12197621522454217d-22 /)

logical,dimension(7)::converged !# set if converged

integer,dimension(7)::last_good_estlim,estlim_to_take !# test, if last value is bad

real(dp),dimension(0:potential%dipole%taylorOrder,7)::estlimDiff !# difference of estlim(m) to estlim(m-2) jjo=convergence_order; converged=.false.; convergence_order= 0

PFLAG= 0; J2pade= 0._dp; diffD= 0._dp; sumD = 0._dp; estlimDiff= 0._dp D= 0._dp; D2= 0._dp

AyS =sqrt(zklij%A(2)) !# abbreviation of sqrt(ay)

AyS_pow= 1._dp !# AyS**0

Xi(0)= 1._dp !# zero order

taylor_coefficients:dom=1,potential%dipole%taylorOrder

AyS_pow(m)=AyS*AyS_pow(m-1) !# powers of AyS

Xi(m)=Xic(m)*AyS_pow(m) !# coefficients for taylor expansion in AyS*rho end dotaylor_coefficients

raise_chebyshev_order:doo=1,10 !# 1:1; 2:2; 3:4, 4:8, 5:16; 6:32; 7:64; 8:128 if(o== 1)then

corder(1)= 0 !# starting point of chebyshev nodes

corder(2)= 1 !# end point of chebyshev nodes

corder(3)= 1 !# increment for chebyshev nodes

else

if(o== 2)D=D/ 2._dp if(o== 2)D2=D2/ 2._dp

corder(1)= 1 !# starting point of chebyshev nodes

corder(2)=corder(2)*2 !# end point of chebyshev nodes

corder(3)= 2 !# increment for chebyshev nodes

end if

Chebyshev:dok=corder(1),corder(2),corder(3)

x=cos(k*Pi/corder(2)) !# chebyshev node

Am =(zklij%A(1)+zklij%A(3)+(zklij%A(1)-zklij%A(3))*x)/4._dp !# quadratic term ofexp[−A_mr²/2 +p_mr]

Am2S=sqrt(2._dp/Am) !# another abbreviation (for expansion)

cos_sum:doi= 1,2 !# sum over±cos

sin_sum:doj= 1,2 !# sum over±sin

xc=signSC(i)*sqrt((1._dp+x)/2) !# xc with sign∓ xs=signSC(j)*sqrt((1._dp-x)/2) !# xs with sign∓

pm=imaginary*(zklij%p(1)*xc+zklij%p(3)*xs)/2._dp !# linear term ofexp[−A_mr²/2 +pmr]

b= -sqrt(2._dp/Am)*pm/2._dp !# abbreviation for completing the square

erf_arg=imaginary*b !# Faddeeva function argument

CALLWOFZ(dreal(erf_arg),aimag(erf_arg),U,V,FLAG) !# Faddeeva function

A.6. Excerpt from the source code

fsigma(1)=(c(0)*xs-c(1)*xc)**2 !#f₀ : factor fromk_z² fsigma(2)=fsigma(1)*(c(0)*xc+c(1)*xs) !#f_x : factor fromk²_z k_x fsigma(3)=fsigma(1)*(c(0)*xs-c(1)*xc) !#f_z : factor fromk³_z fsigma(4)=fsigma(2)*(c(0)*xc+c(1)*xs) !#f_xx : factor fromk_z²k²_x fsigma(5)=fsigma(1) !#f_{f yy}: factor fromk_z²k²_y fsigma(6)=fsigma(1)**2 !#f_zz : factor fromk_z⁴ fsigma(7)=fsigma(2)*(c(0)*xs-c(1)*xc) !#f_xz : factor fromk³_z k_x Ik(0)=sqrt(Pi)/2._dp*cmplx(U,V,dp) !# starting point for recursion (even) Ik(1)= 0.5_dp-b*Ik(0) !# starting point for recursion (odd) dom=2,potential%dipole%taylorOrder+4

Ik(m)=(m-1)/2._dp*Ik(m-2)-b*Ik(m-1) !# Recursion forR∞

b [x^kexp[−x²]]

end do

taylor_terms:dom=0,potential%dipole%taylorOrder

D(m,1:7)=D(m,1:7)+Ik(m+2+addPow(1:7))*Am2S**(m+3+addPow(1:7))*fsigma(1:7) end dotaylor_terms

end dosin_sum end docos_sum end doChebyshev

!###### convergence check #######################################################################

if(o>3)then

converge_seven:dojj=1,7

if(converged(jj))cycleconverge_seven

sumD(o-3,jj)=sum(D(0:potential%dipole%taylorOrder,jj)*Xi(0:potential%dipole%taylorOrder))*Pi/corder(2) if(o>4)then

diffD(o-4,jj)=abs(sumD(o-3,jj)-sumD(o-4,jj)) if(diffD(o-4,jj)<TOL*max(1._dp,abs(sumD(o-3,jj))))then

converged(jj)=.true.

if(corder(2)>convergence_order)convergence_order=corder(2) end if

end if

end doconverge_seven end if

if(constantCebyshevOrder)then

if(o==potential%dipole%chebyshevOrder)exitraise_chebyshev_order else

if(all(converged))exitraise_chebyshev_order end if

end doraise_chebyshev_order est= 0._dp

pade:dom=0,potential%dipole%taylorOrder !# Pade approximation dojj=1,7

if(m>0)then

D(m,jj)=Xi(m)*D(m,jj)

D(m,jj)=D(m,jj)+D(m-1,jj) !# partial sums

end if

CALLEPSAL2(D(m,jj),m,J2pade(:,jj),m,estlim,PFLAG(jj)) !# -Wynn extrapolation if(mod(m,2)== 0)est(m,jj)=estlim

if(printexp)then

if(jj==jjo)print*,D(m,jjo),J2pade(m,jjo),estlim,PFLAG(jjo) end if

end do end dopade

dojj=1,7 !# quick analysis for failors of

last_good_estlim(jj)= 0 !# the pade-approximation at high orders analyze_up:dom=2,potential%dipole%taylorOrder,2

estlimDiff(m,jj)=abs(est(m,jj)-est(m-2,jj)) if(estlimDiff(m,jj)< 100*abs(D(m,jj)))then

last_good_estlim(jj)=m else

exitanalyze_up end if

end doanalyze_up

if(last_good_estlim(jj)== 0)then last_good_estlim(jj)= 2

J2(jj)=D(potential%dipole%taylorOrder,jj) else

estlim_to_take(jj)=last_good_estlim(jj) analyze_down:dom=last_good_estlim(jj),2,-2

if(estlimDiff(m,jj)>estlimDiff(m-2,jj) )then estlim_to_take(jj)=m-2

else

exitanalyze_down end if

end doanalyze_down

J2(jj)=est(estlim_to_take(jj),jj) end if

end do

J2=J2/(corder(2)*4._dp) !#J₂ · weight/( 4 π)

!# special case ofJ_y2-Integral :

J20yy=exp(-(zklij%p(1)**2/zklij%A(1)+zklij%p(3)**2/zklij%A(3))/4._dp)*&

( 8._dp*zklij%A(1)**2*zklij%A(3)*c(0)**2&

&+ 8._dp*zklij%A(3)**2*zklij%A(1)*c(1)**2&

&- 4._dp*(-zklij%A(3)*zklij%p(1)*c(1)+zklij%A(1)*zklij%p(3)*c(0))**2 ) &

&/(sqrt(zklij%A(1))*zklij%A(1)**2*AyS*sqrt(zklij%A(3))*zklij%A(3)**2*sqrt(Pi) ) J2(5)=J20yy-J2(5)

end subroutineJ2_numerical_z

Bibliography

[1] P. J. Mohr, B. N. Taylor, and D. B. Newell. CODATA recommended values of the fundamental physical constants: 2010. Rev. Mod. Phys. 84, 1527 (2012).

[2] S. Bose. Plancks Gesetz und Lichtquantenhypothese. Zeitschrift f¨ur Physik 26, 178 (1924).

[3] A. Einstein. Quantentheorie des einatomigen idealen Gases.Sitzungsberichte der Preußischen Akademie der Wissenschaften 1924, 261 (1924).

[4] A. Einstein. Quantentheorie des einatomigen idealen Gases II. Sitzungs-berichte der preußischen Akademie der Wissenschaften 1925, 3 (1925).

[5] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A.

Cornell. Observation of Bose-Einstein Condensation in a Dilute Atomic Va-por. Science 269, 198 (1995).

[6] C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet. Evidence of Bose-Einstein Condensation in an Atomic Gas with Attractive Interactions.

Phys. Rev. Lett. 75, 1687 (1995).

[7] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle. Bose-Einstein Condensation in a Gas of Sodium Atoms. Phys. Rev. Lett. 75, 3969 (1995).

[8] B. DeMarco and D. S. Jin. Onset of Fermi Degeneracy in a Trapped Atomic Gas. Science 285, 1703 (1999).

[9] A. G. Truscott, K. E. Strecker, W. I. McAlexander, G. B. Partridge, and R. G. Hulet. Observation of Fermi Pressure in a Gas of Trapped Atoms.

Science 291, 2570 (2001).

[10] F. Schreck, L. Khaykovich, K. L. Corwin, G. Ferrari, T. Bourdel, J. Cubi-zolles, and C. Salomon. Quasipure Bose-Einstein Condensate Immersed in a Fermi Sea. Phys. Rev. Lett. 87, 080403 (2001).

[11] Z. Hadzibabic, C. A. Stan, K. Dieckmann, S. Gupta, M. W. Zwierlein, A. G¨orlitz, and W. Ketterle. Two-Species Mixture of Quantum Degener-ate Bose and Fermi Gases. Phys. Rev. Lett. 88, 160401 (2002).

[12] A. Fetter and J. Walecka. Quantum Theory of Many-particle Systems. Dover Books on Physics. Dover Publications (2003).

[13] J. F. Annett. Superconductivity, Superfluids, and Condensates. Oxford Mas-ter Series in Condensed MatMas-ter Physics. Oxford University Press (2004).

[14] M. Ueda. Fundamentals and new frontiers of Bose-Einstein condensation.

World Scientific Publishing Co. Pte. Ltd. 1st ed. (2010).

[15] S. Inouye, M. R. Andrews, J. Stenger, H.-J. Miesner, D. M. Stamper-Kurn, and W. Ketterle. Observation of Feshbach resonances in a Bose-Einstein condensate. Nature 392, 151 (1998).

[16] S. L. Cornish, N. R. Claussen, J. L. Roberts, E. A. Cornell, and C. E.

Wieman. Stable ⁸⁵Rb Bose-Einstein Condensates with Widely Tunable In-teractions. Phys. Rev. Lett. 85, 1795 (2000).

[17] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. Pfau. Bose-Einstein Condensation of Chromium. Phys. Rev. Lett. 94, 160401 (2005).

[18] M. Lu, S. H. Youn, and B. L. Lev. Trapping Ultracold Dysprosium: A Highly Magnetic Gas for Dipolar Physics. Phys. Rev. Lett. 104, 063001 (2010).

[19] K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R. Grimm, and F. Fer-laino. Bose-Einstein Condensation of Erbium. Phys. Rev. Lett.108, 210401 (2012).

[20] K.-K. Ni, S. Ospelkaus, M. H. G. de Miranda, A. Pe’er, B. Neyenhuis, J. J.

Zirbel, S. Kotochigova, P. S. Julienne, D. S. Jin, and J. Ye. A High Phase-Space-Density Gas of Polar Molecules. Science 332, 231 (2008).

[21] R. E. Rosensweig. Ferrohydrodynamics. Dover Publications Inc., New York/US (2014).

[22] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau. The physics of dipolar bosonic quantum gases.Reports on Progress in Physics72, 126401 (2009).

[23] W. Krauth. Quantum Monte Carlo Calculations for a Large Number of Bosons in a Harmonic Trap. Phys. Rev. Lett.77, 3695 (1996).

Bibliography

[24] S. Giorgini, J. Boronat, and J. Casulleras. Ground state of a homogeneous Bose gas: A diffusion Monte Carlo calculation.Phys. Rev. A60, 5129 (1999).

[25] I. Tikhonenkov, B. A. Malomed, and A. Vardi. Anisotropic Solitons in Dipo-lar Bose-Einstein Condensates. Phys. Rev. Lett. 100, 090406 (2008).

[26] R. Eichler, J. Main, and G. Wunner. Variational calculations for anisotropic solitons in dipolar Bose-Einstein condensates. Phys. Rev. A 83, 053604 (2011).

[27] S. Rau, J. Main, and G. Wunner. Variational methods with coupled Gaus-sian functions for Bose-Einstein condensates with long-range interactions. I.

General concept. Phys. Rev. A 82, 023610 (2010).

[28] S. Rau, J. Main, H. Cartarius, P. K¨oberle, and G. Wunner. Variational methods with coupled Gaussian functions for Bose-Einstein condensates with long-range interactions. II. Applications. Phys. Rev. A 82, 023611 (2010).

[29] S. Rau. Variational methods with coupled Gaussian functions for Bose-Einstein condensates with long range interaction. Diplomarbeit, Universit¨at Stuttgart (2010).

[30] R. Eichler. Variationsrechnungen zu anisotropen Solitonen in dipolaren Bose-Einstein-Kondensaten. Diplomarbeit, Universit¨at Stuttgart (2010).

[31] P. Pedri and L. Santos. Two-Dimensional Bright Solitons in Dipolar Bose-Einstein Condensates. Phys. Rev. Lett. 95, 200404 (2005).

[32] P. K¨oberle, D. Zajec, G. Wunner, and B. A. Malomed. Creating two-dimensional bright solitons in dipolar Bose-Einstein condensates. Physical Review A 85, 023630 (2012).

[33] P. K¨oberle. Ground-state structures and dynamics of dipolar Bose-Einstein condensates in single and multi- layered traps. Dissertation, Universit¨at Stuttgart (2011).

[34] R. Eichler, D. Zajec, P. K¨oberle, J. Main, and G. Wunner. Collisions of anisotropic two-dimensional bright solitons in dipolar Bose-Einstein conden-sates. Phys. Rev. A 86, 053611 (2012).

[35] E. J. Heller. Time-dependent approach to semiclassical dynamics. Chem.

Phys. 62 (1975).

[36] E. J. Heller. Classical S-matrix limit of wave packet dynamics. The Journal of Chemical Physics 65, 4979 (1976).

[37] E. J. Heller. Frozen Gaussians: A very simple semiclassical approximation.

The Journal of Chemical Physics 75, 2923 (1981).

[38] T. Fabˇciˇc. Wave packet dynamics in atomic systems and Bose-Einstein con-densates. Dissertation, Universit¨at Stuttgart (2008).

[39] S. Rau, J. Main, P. K¨oberle, and G. Wunner. Pitchfork bifurcations in blood-cell-shaped dipolar Bose-Einstein condensates. Phys. Rev. A 81, 031605(R) (2010).

[40] C. Menotti, C. Trefzger, and M. Lewenstein. Metastable States of a Gas of Dipolar Bosons in a 2D Optical Lattice. Phys. Rev. Lett. 98, 235301 (2007).

[41] M. Baranov, L. Dobrek, K. G´oral, L. Santos, and M. Lewenstein. Ultracold Dipolar Gases – a Challenge for Experiments and Theory. Physica Scripta 2002, 74 (2002).

[42] G. J. Milburn, J. Corney, E. M. Wright, and D. F. Walls. Quantum dynamics of an atomic Bose-Einstein condensate in a double-well potential.Phys. Rev.

A55, 4318 (1997).

[43] G. Theocharis, P. G. Kevrekidis, D. J. Frantzeskakis, and P. Schmelcher.

Symmetry breaking in symmetric and asymmetric double-well potentials.

Phys. Rev. E 74, 056608 (2006).

[44] M. Asad-uz Zaman and D. Blume. Aligned dipolar Bose-Einstein condensate in a double-well potential: From cigar shaped to pancake shaped.Phys. Rev.

A80, 053622 (2009).

[45] B. Xiong, J. Gong, H. Pu, W. Bao, and B. Li. Symmetry breaking and self-trapping of a dipolar Bose-Einstein condensate in a double-well potential.

Phys. Rev. A79, 013626 (2009).

[46] T. Lahaye, T. Pfau, and L. Santos. Mesoscopic Ensembles of Polar Bosons in Triple-Well Potentials. Phys. Rev. Lett. 104, 170404 (2010).

[47] D. Peter, K. Paw lowski, T. Pfau, and K. Rza˙zewski. Mean-field description of dipolar bosons in triple-well potentials. J. Phys. B: At. Mol. Opt. Phys.

45, 225302 (2012).

[48] R. Fortanier, D. Zajec, J. Main, and G. Wunner. Dipolar Bose-Einstein con-densates in triple-well potentials. Journal of Physics B: Atomic, Molecular and Optical Physics 46, 235301 (2013).

Bibliography

[49] C. M. Bender and S. Boettcher. Real Spectra in Non-Hermitian Hamiltonians Having PT Symmetry. Phys. Rev. Lett. 80, 5243 (1998).

[50] C. M. Bender, S. Boettcher, and P. N. Meisinger. PT-symmetric quantum mechanics. Journal of Mathematical Physics 40, 2201 (1999).

[51] A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G. A. Siviloglou, and D. N. Christodoulides. Observation ofPT -Symmetry Breaking in Complex Optical Potentials. Phys. Rev. Lett. 103, 093902 (2009).

[52] C. E. R¨uter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, and D. Kip. Observation of parity-time symmetry in optics. Nat. Phys. 6, 192 (2010).

[53] A. Regensburger, C. Bersch, M.-A. Miri, G. Onishchukov, D. N.

Christodoulides, and U. Peschel. Parity-time synthetic photonic lattices.

Nature 488, 167 (2012).

[54] E. M. Graefe, H. J. Korsch, and A. E. Niederle. Mean-Field Dynamics of a Non-Hermitian Bose-Hubbard Dimer. Phys. Rev. Lett. 101, 150408 (2008).

[55] E. M. Graefe, U. G¨unther, H. J. Korsch, and A. E. Niederle. A non-Hermitian PT–symmetric Bose-Hubbard model: eigenvalue rings from unfolding higher-order exceptional points.Journal of Physics A: Mathematical and Theoretical 41, 255206 (2008).

[56] E.-M. Graefe, H. J. Korsch, and A. E. Niederle. Quantum-classical cor-respondence for a non-Hermitian Bose-Hubbard dimer. Phys. Rev. A 82, 013629 (2010).

[57] H. Cartarius, D. Haag, D. Dast, and G. Wunner. Nonlinear Schr¨odinger equation for aPT-symmetric delta-function double well. Journal of Physics A: Mathematical and Theoretical 45, 444008 (2012).

[58] D. Dast, D. Haag, H. Cartarius, G. Wunner, R. Eichler, and J. Main. A Bose-Einstein condensate in a PT symmetric double well. Fortschritte der Physik 61, 124 (2013).

[59] D. Dast, D. Haag, H. Cartarius, G. Wunner, R. Eichler, and J. Main. De-scription of Bose-Einstein condensates in PT-symmetric double wells. In Proceedings of the International Symposium on Self-Organization in Com-plex Systems: The Past, Present, and Future of Synergetics (2013). In press.

[60] H. Cartarius, D. Dast, D. Haag, G. Wunner, R. Eichler, and J. Main. Sta-tionary and dynamical solutions of the Gross-Pitaevskii equation for a Bose-Einstein condensate in a PT-symmetric double well. Acta Polytechnica 53, 259 (2013).

[61] R. Fortanier, J. Main, and G. Wunner. Dipolar Bose-Einstein condensates inPT-symmetric double-well potential. To be published.

[62] A. D. McLachlan. A variational solution of the time-dependent Schr¨odinger equation. Molecular Physics 8, 39 (1964).

[63] C. A. Regal, M. Greiner, and D. S. Jin. Observation of Resonance Conden-sation of Fermionic Atom Pairs. Phys. Rev. Lett.92, 040403 (2004).

[64] T. Fließbach. Statistische Physik. Lehrbuch zur Theoretischen Physik IV.

Spektrum Akademischer Verlag Heidelberg 5th ed. (2010).

[65] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga. Feshbach Resonances in Ultracold gases. Rev. Mod. Phys. 82, 1225 (2010).

[66] E. Gross. Structure of a quantized vortex in boson systems.Il Nuovo Cimento Series 10 20, 454 (1961).

[67] L. P. Pitaevskii. Vortex lines in an imperfect Bose gas. Soviet Physics JETP-USSR 13 (1961).

[68] R. Gut¨ohrlein, J. Main, H. Cartarius, and G. Wunner. Bifurcations and exceptional points in dipolar Bose-Einstein condensates. Journal of Physics A: Mathematical and Theoretical 46, 305001 (2013).

[69] H. Cartarius, J. Main, and G. Wunner. Discovery of exceptional points in the Bose-Einstein condensation of gases with attractive 1/r interaction. Phys.

Rev. A77, 013618 (2008).

[70] J. J. Rasmussen and K. Rypdal. Blow-up in Nonlinear Schroedinger Equations-I A General Review. Physica Scripta 33, 481 (1986).

[71] M. D. Feit, J. A. Fleck, Jr., and A. Steiger. Solution of the Schr¨odinger Equation by a Spectral Method. J. Comp. Phys. 47, 412 (1982).

[72] P. K¨oberle, D. Zajec, G. Wunner, and B. A. Malomed. Creating two-dimensional bright solitons in dipolar Bose-Einstein condensates. Phys. Rev.

A85, 023630 (2012).

Bibliography

[73] D. Zajec. Dipolar Bose-Einstein condensates. Dissertation, Universit¨at Stuttgart (in preperation).

[74] M. Kreibich, J. Main, and G. Wunner. Relation between the eigenfrequencies of Bogoliubov excitations of Bose-Einstein condensates and the eigenvalues of the Jacobian in a time-dependent variational approach. Phys. Rev. A 86, 013608 (2012).

[75] E. Faou and C. Lubich. A Poisson integrator for Gaussian wavepacket dy-namics. Computing and Visualization in Science 9, 45 (2006).

[76] P. Weinmann. Setting up a high-resolution optical system for trapping and imaging ultracold chromium atoms. Diplomarbeit, Universit¨at Stuttgart (2011).

[77] L. Santos, G. V. Shlyapnikov, P. Zoller, and M. Lewenstein. Bose-Einstein Condensation in Trapped Dipolar Gases. Phys. Rev. Lett. 85, 1791 (2000).

[78] P. Wagner. Bose-Einstein condensates with electromagnetically induced 1/r-potential. Diplomarbeit, Universit¨at Stuttgart (2007).

[79] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numer-ical Recipes in Fortran. Cambridge University Press (1997).

[80] I. Gradshteyn and I. W. Ryzhik. Tables of Integrals, Series and Products.

Academic Press New York (1965).

[81] G. P. M. Poppe and C. M. J. Wijers. Algorithm 680: Evaluation of the Complex Error Function. ACM Transactions on Mathematical Software 16, 47 (1990).

[82] W. Lang and E. W. Weisstein. The Online Encyclopedia of Integer Se-quences. http://oeis.org/ (2013).

[83] M. Abramowitz and I. A. Stegun. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover New York (1964).

[84] B. C. Carlson and J. L. Gustafson. Asymptotic approximations for symmetric elliptic integrals. SIAM J. Math. Anal. 25, 288 (1994).

[85] B. C. Carlson. Numerical computation of real or complex elliptic integrals.

Numerical Algorithms 10, 13 (1995).

[86] R. Brankin, I. Gladwell, and L. Shampine. RKSUITE. Numerical Al-gorithms Group Ltd. Wilkinson House, Oxford and Department of Math-ematics Southern Methodist University, Dallas; 1st ed. (1991). http:

//www.netlib.org/ode/rksuite/.

[87] C. C. Bradley, C. A. Sackett, and R. G. Hulet. Analysis of in situ images of Bose-Einstein condensates of lithium. Phys. Rev. A55 (1997).

[88] B. G. Burton, E. H. Kenneth, J. Jorge, and more. Powell hybrid method from minpack. University of Chicago (1999).http://www.netlib.org/minpack/.

[89] The GNOME Project and PyGTK Team. PyGTK: GTK+ for Python.

http://www.pygtk.org/(2004-2010).

[90] J. D. Hunter. Matplotlib: A 2D graphics environment. Computing In Science

& Engineering 9, 90 (2007).

[91] P. Ramachandran and G. Varoquaux. Mayavi: Making 3D data visualization reusable. In G. Varoquaux, T. Vaught, and J. Millman, editors, Proceed-ings of the 7th Python in Science Conference pp. 51–56 Pasadena, CA USA (2008).

[92] W. Schroeder, K. Martin, and B. Lorensen. The Visualization Toolkit, Third Edition. Kitware Inc. (2007).

[93] V. E. Zakharov and A. B. Shabat. Interaction between solitons in a stable medium. JETP 37, 1627 (1972).

[94] K. E. Strecker, G. B. Partridge, A. G. Truscott, and R. G. Hulet. Formation and propagation of matter-wave soliton trains. Nature 417, 150 (2002).

[95] L. Khaykovich, F. Schreck, G. Ferrari, T. Bourdel, J. Cubizolles, L. D. Carr, Y. Castin, and C. Salomon. Formation of a Matter-Wave Bright Soliton.

Science 296, 1290 (2002).

[96] S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Sengstock, A. Sanpera, G. V. Shlyapnikov, and M. Lewenstein. Dark Solitons in Bose-Einstein Con-densates. Phys. Rev. Lett. 83, 5198 (1999).

[97] J. Denschlag, J. E. Simsarian, D. L. Feder, C. W. Clark, L. A. Collins, J. Cubizolles, L. Deng, E. W. Hagley, K. Helmerson, W. P. Reinhardt, S. L.

Rolston, B. I. Schneider, and W. D. Phillips. Generating Solitons by Phase Engineering of a Bose-Einstein Condensate. Science 287, 97 (2000).

Bibliography

[98] C. Becker, S. Stellmer, P. Soltan-Panahi, S. Dorscher, M. Baumert, E.-M.

Richter, J. Kronjager, K. Bongs, and K. Sengstock. Oscillations and inter-actions of dark and dark-bright solitons in Bose-Einstein condensates. Nat.

Phys. 4, 496 (2008).

[99] B. Eiermann, T. Anker, M. Albiez, M. Taglieber, P. Treutlein, K.-P. Marzlin, and M. K. Oberthaler. Bright Bose-Einstein Gap Solitons of Atoms with Repulsive Interaction. Phys. Rev. Lett. 92, 230401 (2004).

[100] S. K. Adhikari and P. Muruganandam. Two-dimensional dipolar Bose-Einstein condensate bright and vortex solitons on a one-dimensional optical lattice . J. Phys. B: At. Mol. Opt. Phys. 45, 045301 (2012).

[101] S. Giovanazzi, A. G¨orlitz, and T. Pfau. Tuning the Dipolar Interaction in Quantum Gases. Phys. Rev. Lett. 89, 130401 (2002).

[102] L. E. Young-S, P. Muruganandam, and S. K. Adhikari. Dynamics of quasi-one-dimensional bright and vortex solitons of a dipolar Bose-Einstein con-densate with repulsive atomic interaction. J. Phys. B: At. Mol. Opt. Phys.

44 (2011).

[103] N. F. Mott and H. S. W. Massey. The Theory of Atomic Collisions. Oxford University Press 3rd ed. (1965).

[104] S. Raghavan, A. Smerzi, S. Fantoni, and S. R. Shenoy. Coherent oscillations between two weakly coupled Bose-Einstein condensates: Josephson effects, πoscillations, and macroscopic quantum self-trapping. Phys. Rev. A59, 620 (1999).

[105] A.-X. Zhang and J.-K. Xue. Dipolar-induced interplay between inter-level physics and macroscopic phase transitions in triple-well potentials. J. Phys.

B: At. Mol. Opt. Phys.45, 145305 (2012).

[106] T. Lahaye, T. Koch, B. Fr¨ohlich, M. Fattori, J. Metz, A. Griesmaier, S. Gio-vanazzi, and T. Pfau. Strong dipolar effects in a quantum ferrofluid. Nature 448, 672 (2007).

[107] N. Moiseyev. Non-Hermitian quantum mechanics. Cambridge University Press (2011).

[108] N. Moiseyev and L. S. Cederbaum. Resonance solutions of the nonlinear Schr¨odinger equation: Tunneling lifetime and fragmentation of trapped con-densates. Phys. Rev. A 72, 033605 (2005).

[109] P. Schlagheck and T. Paul. Complex-scaling approach to the decay of Bose-Einstein condensates. Phys. Rev. A73, 023619 (2006).

[110] T. Paul, M. Hartung, K. Richter, and P. Schlagheck. Nonlinear transport of Bose-Einstein condensates through mesoscopic waveguides. Phys. Rev. A 76, 063605 (2007).

[111] F. K. Abdullaev, V. V. Konotop, M. Salerno, and A. V. Yulin. Dissipative periodic waves, solitons, and breathers of the nonlinear Schr¨odinger equation with complex potentials. Phys. Rev. E 82, 056606 (2010).

[112] Y. V. Bludov and V. V. Konotop. Nonlinear patterns in Bose-Einstein con-densates in dissipative optical lattices. Phys. Rev. A 81, 013625 (2010).

[113] H. Cartarius and G. Wunner. Model of a PT–symmetric Bose-Einstein con-densate in a δ-function double-well potential. Phys. Rev. A 86, 013612 (2012).

[114] P. W¨urtz, T. Langen, T. Gericke, A. Koglbauer, and H. Ott. Experimental Demonstration of Single-Site Addressability in a Two-Dimensional Optical Lattice. Phys. Rev. Lett. 103, 080404 (2009).

[115] Y. Shin, G.-B. Jo, M. Saba, T. A. Pasquini, W. Ketterle, and D. E. Pritchard.

Optical Weak Link between Two Spatially Separated Bose-Einstein Conden-sates. Phys. Rev. Lett. 95, 170402 (2005).

[116] R. Gati, M. Albiez, J. F¨olling, B. Hemmerling, and M. Oberthaler. Realiza-tion of a single Josephson juncRealiza-tion for Bose-Einstein condensates. Applied Physics B 82, 207 (2006).

[117] D. Dast, D. Haag, H. Cartarius, J. Main, and G. Wunner. Eigenvalue struc-ture of a Bose-Einstein condensate in aPT-symmetric double well. J. Phys.

A46, 375301 (2013).

[118] E.-M. Graefe. Stationary states of a PT symmetric two-mode Bose-Einstein condensate. Journal of Physics A: Mathematical and Theoretical 45, 444015 (2012).

[119] D. O’Dell, S. Giovanazzi, G. Kurizki, and V. M. Akulin. Bose-Einstein Condensates with 1/r Interatomic Attraction: Electromagnetically Induced

“Gravity”. Phys. Rev. Lett. 84, 5687 (2000).

[120] M. Abad, M. Guilleumas, R. Mayol, M. Pi, and D. M. Jezek. Dipolar con-densates confined in a toroidal trap: Ground state and vortices. Phys. Rev.

A81, 043619 (2010).

Bibliography

[121] W. Johnson. The Curious History of Fa`a di Bruno’s Formula. The American Mathematical Monthly 109, 217 (2002).

[122] P. J¨ackel. A note on multivariate Gauss-Hermite quadrature. http://

www.jaeckel.org/ANoteOnMultivariateGaussHermiteQuadrature.pdf (2005). Website status: Aug. 15, 2013.

[123] G. H. Golub and J. Kautsky. Calculation of Gauss Quadratures with Multiple Free and Fixed Knots. Numer. Math.41, 147 (1983).

[124] P. J¨ackel. Roots and weights for radial

Gauss-Hermite quadrature. http://www.jaeckel.org/

RootsAndWeightsForRadialGaussHermiteQuadrature.cpp (2005). Web-site status: Aug. 15, 2013.

[125] P. Wynn. Upon systems of recursions which obtain among the quotients of the Pad´e table. Numerische Mathematik 8, 264 (1966).

[126] W. Gautschi. Efficient computation of the complex error function. SIAM Journal on Numerical Analysis 7, 187 (1970).

Zusammenfassung in deutscher Sprache

In dieser Arbeit werden dipolare Bose-Einstein-Kondensate (BECs) auf Grundla-ge der Gross-Pitaevskii-Gleichung (GPE) mit Hilfe des zeitabh¨angiGrundla-gen Variations-prinzips (TDVP) untersucht. BECs sind seit den theoretischen Vorhersagen von Bose und Einstein [2–4] und ihrer experimentellen Realisierung [5–7] Gegenstand umfangreicher theoretischer wie experimenteller Forschung. In Experimenten mit

52Cr-Atomen [17], welche ein großes magnetisches Dipolmoment aufweisen, gelang es, BECs mit dipolarer Wechselwirkung herzustellen. Das Feld dipolarer BECs entwickelte sich seitdem rasant und zahlreiche physikalische Effekte konnten in der Theorie beschrieben und im Experiment realisiert werden. Die Dipol-Dipol-Wechselwirkung (DDI) zeigt ein vollst¨andig anderes Verhalten als beispielsweise die Streuwechselwirkung, die in der Beschreibung nicht-dipolarer BECs verwendet werden kann. Faszinierende Effekte treten bereits in klassischen Systemen auf, wie zum Beispiel die Rosensweig-Instabilit¨at [21] in Ferrofluiden. Die neue, interessante und manchmal ¨uberraschende Physik wird durch die besonderen Eigenschaften der DDI hervorgerufen. Es handelt sich um eine langreichweitige, lokale, nicht-lineare und anisotrope Wechselwirkung. Diese speziellen Besonderheiten sind Ur-sprung einer Vielzahl interessanter Effekte wie z.B. der Verformung der atomaren Wolke, strukturierter Grundzust¨ande, des Roton-Maxon Spektrums, unterschied-licher Kollapsdynamik, von Solitonen, der Strukturbildung, neuer Quantenphasen und vielem mehr.

Die theoretische Beschreibung in dieser Arbeit basiert auf der Mean-Field N¨aher-ung, welche auf die dipolare Gross-Pitaevskii-Gleichung f¨uhrt. Dabei handelt es sich um eine nichtlineare Schr¨odingergleichung. Eine etablierte Methode zur L¨o-sung der GPE ist die Repr¨asentation der Wellenfunktion auf einem numerischen Gitter und eine anschließende L¨osung in Real- oder Imagin¨arzeit mit der Split-Operator-Methode. Dieses mathematisch einfache Konzept ist f¨ur die dipolare GPE mit mehreren Schwierigkeiten verbunden. So muss nicht nur die Wellenfunk-tion, sondern auch das langreichweitige Dipolpotential auf das Gitter abgebildet werden. Dies macht die Methode sehr rechenaufw¨andig und erfordert zum Teil spezielle Hardware zum Erreichen angemessener Rechenzeiten. Des weiteren sind station¨are Zust¨ande jenseits des Grundzustandes nicht ohne weiteres verf¨ugbar.

Ein zentraler Bestandteil dieser Arbeit ist daher die Entwicklung einer alternati-ven Methode zur L¨osung der zeitabh¨angigen dipolaren GPE. Diese Methode wird auf verschiedene physikalische Systeme und Fragestellungen angewandt werden.

Insbesondere die Kollision zweier quasi zweidimensionaler Solitonen, die Unter-suchung dipolarer Kondensate im Dreimuldenpotential, sowie die Analyse eines BECs im PT-symmetrischen Doppelmuldenpotential stellen Anwendungen dar.

In Kapitel 2 wird in K¨urze eine Einf¨uhrung in die theoretische Beschreibung von Bose-Einstein-Kondensaten gegeben. Dazu wird inAbschnitt 2.1dieses Ph¨a-nomen zun¨achst durch thermodynamische Argumentation am idealen Bose-Gas beschrieben. Wir zeigen hier, dass es sich bei der Bose-Einstein Kondensation um einen thermodynamischen Phasen¨ubergang handelt. In Abschnitt 2.2 f¨uh-ren wir in die Beschreibung ultra-kalter atomarer Gase ein, erl¨autern das Konzept der “off-diagonal long-range order” (ODLRO) und der makroskopischen Wellen-funktion. Dabei stellen wir die N¨aherungen dar, welche im folgenden Abschnitt zur Gross-Pitaevskii-Gleichung f¨uhren. Nach der Rechtfertigung der Mean-Field-N¨aherung wenden wir uns in Abschnitt 2.3 der zentralen Gleichung in dieser Arbeit, der dipolaren Gross-Pitaevskii-Gleichung zu. Dabei werden in den Unter-abschnitten 2.3.1 und 2.3.2 die kurzreichweitige Kontaktwechselwirkung und die langreichweitige Dipol-Dipol-Wechselwirkung (DDI) diskutiert. Die erweiter-te zeitabh¨angige Gross-Pitaevskii-Gleichung, wie sie der weierweiter-teren Arbeit zugrun-de liegt, wird in Unterabschnitt 2.3.3 angegeben und mit ihren Eigenschaften und Konsequenzen, die sich beispielweise aus ihrer Nichtlinearit¨at ergeben, im fol-genden Unterabschnitt 2.3.4 dargestellt. In Abschnitt 2.4 werden etablierte Methoden zur L¨osung der GPE vorgestellt. Diese beruhen meist auf der Diskreti-sierung der Wellenfunktion auf einem, im allgemeinen dreidimensionalen, Gitter.

F¨ur die Zeitentwicklung eines Wellenpakets in Real- oder Imagin¨arzeit (siehe Un-terabschnitt 2.4.1) wird dann h¨aufig die Split-Operator-Methode, welche auf der Aufteilung des Zeitentwicklungsoperators f¨ur kleine Zeitschritte und jeweils einer Multiplikation im Orts- und Impulsraum beruht, verwendet. Die Transformation der Wellenfunktion von einem in den anderen Raum erfolgt mittels Fast-Fourier-Transformationen (FFTs) und ist f¨ur ein dreidimensionales Gitter mit entsprechen-der Aufl¨osung numerisch sehr aufw¨andig. Wir diskutieren daher die Implementie-rung f¨ur Grafikkarten inUnterabschnitt 2.4.2, die einen hohen Grad an Paralle-lisierbarkeit erm¨oglichen. Zur Untersuchung der Stabilit¨at und zur Berechnung des Anregungsspektrums werden h¨aufig die Bogoliubov-de Gennes-Gleichungen (Bd-GE) verwendet, welche sich aus einer Linearisierung der GPE ergeben. Wir stellen diese kurz in Unterabschnitt 2.4.3vor.

Im Dokument Variational approaches to dipolar Bose-Einstein condensates (Seite 116-158)