Update to QOBase v0.4 and switch to ClausenFunctions for collective modes

Project.toml

@@ -3,24 +3,24 @@ uuid = "abf2e61e-3022-5fcf-a868-69d409dee72b"
 version = "0.1.5"
 
 [deps]
+ClausenFunctions = "1cc96f36-ef02-40ef-a648-5e2c13df3076"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
 QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
-Polylogarithms = "43cc4a2f-1a2b-4fbb-87e3-95c794af13de"
 
 [compat]
+ClausenFunctions = "1"
 Combinatorics = "1"
 DiffEqCallbacks = "2"
 Optim = "1"
 OrdinaryDiffEq = "6"
 QuantumOptics = "1"
-QuantumOpticsBase = "0.3"
-Polylogarithms = "0.1"
-julia = "1.4"
+QuantumOpticsBase = "0.4"
+julia = "1.6"
 
 [extras]
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

docs/Project.toml

@@ -4,6 +4,6 @@ PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
 
 [compat]
-Documenter = "~0.25"
+Documenter = "1"
 PyPlot = "2"
 QuantumOptics = "1"

docs/src/api.md

@@ -54,6 +54,21 @@ geometry.box
 geometry.cube
 ```
 
+```@docs
+geometry.ring
+```
+
+```@docs
+geometry.lhc1
+```
+
+```@docs
+geometry.lhc1_p
+```
+
+```@docs
+geometry.excitation_phases
+```
 
 ## [Dipole-Dipole Interaction](@id API: Dipole-Dipole Interaction)
 

docs/src/geometry.md

@@ -9,6 +9,7 @@ In order to simplify creation of various particle distributions, a few helper fu
 * [`geometry.hexagonal`](@ref)
 * [`geometry.box`](@ref)
 * [`geometry.cube`](@ref)
+* [`geometry.ring`](@ref)
 
 They can be used directly to create a [`SpinCollection`](@ref):
 

src/collective_modes.jl

@@ -6,10 +6,9 @@
 # atomic chains, see Asenjo-Garcia et al 10.1103/PhysRevX.7.031024.
 
 import LinearAlgebra
-import Polylogarithms
+using ClausenFunctions: cl2, cl3
 
 const la = LinearAlgebra
-const pl = Polylogarithms
 const k0 = 2pi
 const Rot90 = [0 -1; 1 0]
 
@@ -58,11 +57,9 @@
 """
 function Omega_k_chain(k::Real, a::Real, polarization::Array{<:Number, 1}) 
  polarization, x_comp, yz_comp = polarization_renorm(polarization)
- exp_sum = exp(im*(k0+k)*a)
- exp_diff = exp(im*(k0-k)*a)
- diag = pl.polylog(3, exp_sum) + pl.polylog(3, exp_diff) - im*k0*a*(pl.polylog(2, exp_sum) + pl.polylog(2, exp_diff))
- perp = (k0*a)^2*(log(1-exp_sum) + log(1-exp_diff))
- return 3/(4*(k0*a)^3) * (-2*x_comp*real(diag) + yz_comp*real(diag+perp))
+ diag = cl3((k0 + k)*a) + cl3((k0 - k)*a) + k0*a*(cl2((k0 + k)*a) + cl2((k0 - k)*a))
+ perp = (k0*a)^2*(log(1 - exp(im*(k0+k)*a)) + log(1 - exp(im*(k0-k)*a)))
+ return 3/(4*(k0*a)^3) * (-2*x_comp*diag + yz_comp*(diag+real(perp)))
 end
 
 
@@ -224,63 +221,63 @@
 * `k_vec`: xy-axis quasimomentum of collective mode in first BZ
 * `a_vec1, a_vec2`: 2D Bravais lattice vectors (real-space).
 """
 function Green_Tensor_k_2D(k_vec::Array{<:Number, 1}, a_vec1::Array{<:Number, 1}, a_vec2::Array{<:Number, 1}; Lambda=nothing, N1=nothing, N2=nothing)
  V = abs(a_vec1[1]*a_vec2[2] - a_vec1[2]*a_vec2[1]) #BZ volume
 
  #reciprocal lattice vectors
  b1 = 2π * Rot90*a_vec2 / la.dot(a_vec1, Rot90*a_vec2)
  b2 = 2π * Rot90*a_vec1 / la.dot(a_vec2, Rot90*a_vec1)
 
  #cutoff parameters
  Lambda = 10/sqrt(V)
  N1 = Int(ceil(5*Lambda/la.norm(b1))) #number of terms in momentum space sum
  N2 = Int(ceil(5*Lambda/la.norm(b2))) #number of terms in momentum space sum
 
  G_tensor = zeros(ComplexF64,3,3)
 
  G_sum_xx = 0
  G_sum_yy = 0
  G_sum_zz = 0
 
  G_sum_xy = 0
  G_sum_yx = 0
  for n1 = -N1:N1
  for n2 =-N2:N2
  k_eff = k_vec + n1*b1 + n2*b2
  kz = sqrt(Complex(k0^2 - la.norm(k_eff)^2))
  if abs(kz) > 0
  cutoff = 1im*exp(-(la.norm(k_eff)/Lambda)^2) / (2*kz)
  else
  cutoff = 0
  end
  G_sum_xx += (k0^2 - k_eff[1]^2)/k0^2 * cutoff
  G_sum_yy += (k0^2 - k_eff[2]^2)/k0^2 * cutoff
  G_sum_zz += la.norm(k_eff)^2/k0^2 * cutoff
 
  G_sum_xy += -k_eff[1]*k_eff[2]/k0^2 * cutoff
  G_sum_yx += -k_eff[1]*k_eff[2]/k0^2 * cutoff
  end
  end
 
  Re_G0_xx = -(Lambda/(32*sqrt(pi)*k0^2)) * exp(-(k0/Lambda)^2) * (Lambda^2 - 2*k0^2)
  Re_G0_yy = Re_G0_xx
  Re_G0_zz = -(Lambda/(32*sqrt(pi)*k0^2)) * exp(-(k0/Lambda)^2) * (-2*Lambda^2 - 4*k0^2)
 
  Gxx = 1/V*G_sum_xx - Re_G0_xx
  Gyy = 1/V*G_sum_yy - Re_G0_yy
  Gzz = 1/V*G_sum_zz - Re_G0_zz
 
  G_xy = 1/V*G_sum_xy
  G_yx = G_xy
 
  G_tensor[1,1] = Gxx
  G_tensor[2,2] = Gyy
  G_tensor[3,3] = Gzz
 
  G_tensor[1,2] = G_xy
  G_tensor[2,1] = G_yx
 
  return G_tensor
 end
 
 end # module