diff --git a/docs/make.jl b/docs/make.jl
index 6992cd9..926e655 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -7,6 +7,7 @@ pages = [
     "interaction.md",
     "effective_interaction.md",
     "system.md",
+    "collective_modes.md",
     "descriptions.md",
     "api.md"
 ]
diff --git a/docs/src/api.md b/docs/src/api.md
index addf502..b4713bb 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -409,12 +409,12 @@ CollectiveSpins.reducedspin.reducedsigmay
 CollectiveSpins.reducedspin.reducedsigmax
 ```
 
-### [Collective modes](@id API: Collective modes)
+## [Collective Modes](@id API: Collective Modes)
 
 ```@docs
-CollectiveSpins.collective_modes.Omega_k_chain(k, a, polarization)
+CollectiveSpins.Omega_k_chain
 ```
 
 ```@docs
-CollectiveSpins.collective_modes.Gamma_k_chain(k, a, polarization)
+CollectiveSpins.Gamma_k_chain
 ```
diff --git a/docs/src/collective_modes.md b/docs/src/collective_modes.md
new file mode 100644
index 0000000..91da1ea
--- /dev/null
+++ b/docs/src/collective_modes.md
@@ -0,0 +1,10 @@
+# Collective Modes
+
+A module that calculates the band structure and decay rates of different collective modes of quasimomentum k with functions:
+
+* [`CollectiveSpins.Omega_k_chain`](@ref)
+* [`CollectiveSpins.Gamma_k_chain`](@ref)
+
+```julia
+CollectiveSpins.Omega_k_chain(0.2*pi, 0.2, [1,0,0])
+```
diff --git a/docs/src/collective_modes_demo.jl b/docs/src/collective_modes_demo.jl
new file mode 100644
index 0000000..a7226c5
--- /dev/null
+++ b/docs/src/collective_modes_demo.jl
@@ -0,0 +1,38 @@
+import CollectiveSpins
+
+# Define chain parameters
+a = 0.2   # spin-spin distance
+polarization_par = [1, 0, 0]
+polarization_perp = [0, 1, 1]
+
+k_max = 10000
+k_list = [iii*pi/(a*k_max) for iii=-k_max:k_max]
+println()
+
+band_structure_par = []
+band_structure_perp = []
+decay_par = []
+decay_perp = []
+for k in k_list
+    append!(band_structure_par, CollectiveSpins.Omega_k_chain(k, a, polarization_par))
+    append!(band_structure_perp, CollectiveSpins.Omega_k_chain(k, a, polarization_perp))
+    append!(decay_par, CollectiveSpins.Gamma_k_chain(k, a, polarization_par))
+    append!(decay_perp, CollectiveSpins.Gamma_k_chain(k, a, polarization_perp))
+end
+
+using PyPlot
+
+subplot(211)
+plot(k_list*a/pi, band_structure_par)
+plot(k_list*a/pi, band_structure_perp, color="red")
+xlabel("k d/pi")
+ylabel("Omega(k)/gamma0")
+
+subplot(212)
+plot(k_list*a/pi, decay_par)
+plot(k_list*a/pi, decay_perp, color="red")
+xlabel("k d/pi")
+ylabel("Gamma(k)/gamma0")
+
+tight_layout()
+savefig("band_structure_decay.svg")
diff --git a/docs/src/index.md b/docs/src/index.md
index a4f2dd9..f4fa2ad 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -2,4 +2,4 @@
 
 **CollectiveSpins.jl** is a numerical framework that can be used to simulate quantum systems consisting of spatially distributed spins interacting via [Dipole-Dipole interaction](@ref), optionally coupled to a cavity.
 
-The [Geometry](@ref) module allows the rapid creation of arbitrarily placed spins in order to build up very general systems as explained in the [System](@ref) documentation. These can then be investigated using either a complete quantum description or cumulant expansions up to second order. The details are presented in [Theoretical Descriptions](@ref).
+The [Geometry](@ref) module allows the rapid creation of arbitrarily placed spins in order to build up very general systems as explained in the [System](@ref) documentation. These can then be investigated using either a complete quantum description or cumulant expansions up to second order. The details are presented in [Theoretical Descriptions](@ref). The [Collective Modes](@ref) provides functions for the calculation of collective spin band structures and decay rates.
diff --git a/src/CollectiveSpins.jl b/src/CollectiveSpins.jl
index 86dbd8d..e22dc33 100644
--- a/src/CollectiveSpins.jl
+++ b/src/CollectiveSpins.jl
@@ -1,7 +1,7 @@
 module CollectiveSpins
 
 export Spin, SpinCollection, CavityMode, CavitySpinCollection,
-        geometry, GreenTensor, G_ij, Gamma_ij, Omega_ij
+        geometry, GreenTensor, G_ij, Gamma_ij, Omega_ij, Omega_k_chain, Delta_k_chain
 
 
 include("system.jl")
diff --git a/src/collective_modes.jl b/src/collective_modes.jl
index 01945a9..d0c6392 100644
--- a/src/collective_modes.jl
+++ b/src/collective_modes.jl
@@ -1,25 +1,27 @@
 module collective_modes
 
+export Omega_k_chain, Gamma_k_chain
+
 # For a particularly elegant derivation of collective frequency shifts and decay rates in 1D
 # atomic chains, see Asenjo-Garcia et al 10.1103/PhysRevX.7.031024.
 
 import LinearAlgebra
 import Polylogarithms
 
-pl = Polylogarithms
+const pl = Polylogarithms
 
 """
     Omega_k_chain(k, a, polarization)
 
 Collective frequency shift Omega_k of mode k for an infinite chain of atoms along x-axis.
 
-WLOG, this calculation scales natural atomic linewidth lambda0=1 and decay rate gamma0=1.
+WLOG, this calculation scales natural atomic frequency wavelength lambda0=1 and decay rate gamma0=1.
 
 # Arguments
 * `k`: x-axis quasimomentum of collective mode in first BZ such that |k|<= pi/a.
 * `a`: Atomic lattice spacing.
 * `polarization`: 3D, complex vector of atomic polarization.
-"""  
+"""
 
 function Omega_k_chain(k::Real, a::Real, polarization::Array{<:Number, 1})  
     polarization = LinearAlgebra.normalize(polarization)
@@ -38,7 +40,7 @@ end
 
 Collective decay rate Gamma_k of mode k for an infinite chain of atoms along x-axis.
 
-WLOG, this calculation scales natural atomic linewidth lambda0=1 and decay rate gamma0=1.
+WLOG, this calculation scales natural atomic frequency wavelength lambda0=1 and decay rate gamma0=1.
 
 # Arguments
 * `k`: x-axis quasimomentum of collective mode in first BZ such that |k|<= pi/a.
@@ -47,7 +49,7 @@ WLOG, this calculation scales natural atomic linewidth lambda0=1 and decay rate
 """
 
 function Gamma_k_chain(k::Real, a::Real, polarization::Array{<:Number, 1})
-    if k > 2pi
+    if abs(k) > 2pi
         return 0
     end
     polarization = LinearAlgebra.normalize(polarization)
@@ -55,12 +57,12 @@ function Gamma_k_chain(k::Real, a::Real, polarization::Array{<:Number, 1})
     k0 = 2pi
     g = 2pi/a
     N = floor((k0 - k)/g)
-    sum = 0
+    Gamma_sum = 0
     for n = -N:N
-        sum += (k + n*g)^2
+        Gamma_sum += (k + n*g)^2
     end
-    sum = sum/k0^2
-    return 3pi/(4*k0*a) * (2*x_par*(2*N+1-sum) + (1-x_par)*(2*N+1+sum))
+    Gamma_sum = Gamma_sum/k0^2
+    return 3pi/(4*k0*a) * (2*x_par*(2*N+1-Gamma_sum) + (1-x_par)*(2*N+1+Gamma_sum))
 end
 
 end # module