Document reducedspin

docs/Project.toml

@@ -1,7 +1,9 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
 
 [compat]
 Documenter = "~0.25"
 PyPlot = "2"
+QuantumOptics = "0.8"

docs/make.jl

@@ -8,6 +8,7 @@ pages = [
     "effective_interaction.md",
     "system.md",
     "collective_modes.md",
+    "reducedspin.md",
     "descriptions.md",
     "api.md"
 ]

docs/src/api.md

@@ -372,39 +372,39 @@ CollectiveSpins.mpc.squeezingparameter
 ### [Reduced Spin](@id API: Methods-reduced)
 
 ```@docs
-CollectiveSpins.reducedspin.reducedsigmap
+ReducedSpinBasis
 ```
 
 ```@docs
-CollectiveSpins.reducedspin.reducedspintransition
+reducedspintransition
 ```
 
 ```@docs
-CollectiveSpins.reducedspin.ReducedSpinBasis
+reducedsigmap
 ```
 
 ```@docs
-CollectiveSpins.reducedspin.reducedsigmam
+reducedsigmam
 ```
 
 ```@docs
-CollectiveSpins.reducedspin.reducedsigmaz
+reducedsigmax
 ```
 
 ```@docs
-CollectiveSpins.reducedspin.reducedspinstate
+reducedsigmay
 ```
 
 ```@docs
-CollectiveSpins.reducedspin.reducedsigmay
+reducedsigmaz
 ```
 
 ```@docs
-CollectiveSpins.reducedspin.reducedsigmax
+reducedsigmapsigmam
 ```
 
 ```@docs
-CollectiveSpins.reducedspin.sigmap_sigmam
+reducedspinstate
 ```
 
 ## [Collective Modes](@id API: Collective Modes)

docs/src/reducedspin.md

@@ -0,0 +1,82 @@
+# Reduced Spin
+
+A common approach to treating many spins quantum mechanically is to truncate the Hilbert space at a specific number of excitations. This widely reduces the number of possible states. **CollectiveSpins.jl** offers generic functionality to truncate a system of spins at an arbitrary number of excitations. The following functions can be used, and are largely equivalent to the `SpinBasis` implemented in **QuantumOptics.jl**.
+
+
+* [`ReducedSpinBasis`](@ref)
+* [`reducedspintransition`](@ref)
+* [`reducedsigmap`](@ref)
+* [`reducedsigmam`](@ref)
+* [`reducedsigmax`](@ref)
+* [`reducedsigmay`](@ref)
+* [`reducedsigmaz`](@ref)
+* [`reducedsigmapsigmam`](@ref)
+* [`reducedspinstate`](@ref)
+
+
+## Example
+
+Let us illustrate how to work with the `reducedspin` submodule. In the following example, we will compute the intensity radiation pattern of a regular chain of atoms in a subradiant state from A. Asenjo-Garcia *et al*, [10.1103/PhysRevX.7.031024](https://doi.org/10.1103/PhysRevX.7.031024) (Fig. 3b).
+
+```@example
+using CollectiveSpins
+using QuantumOptics
+using PyPlot
+
+# Parameters
+N = 50
+M = 1 # Number of excitations
+d = 0.33
+pos = geometry.chain(d,N)
+μ = [[1.0,0,0] for i=1:N]
+S = SpinCollection(pos,μ)
+
+# Collective effects
+Ωmat = OmegaMatrix(S)
+Γmat = GammaMatrix(S)
+
+# Hilbert space
+b = ReducedSpinBasis(N,M,M) # Basis from M excitations up to M excitations
+
+# Effective Hamiltonian
+spsm = [reducedsigmapsigmam(b, i, j) for i=1:N, j=1:N]
+H_eff = dense(sum((Ωmat[i,j] - 0.5im*Γmat[i,j])*spsm[i, j] for i=1:N, j=1:N))
+
+# Find the most subradiant eigenstate
+λ, states = eigenstates(H_eff; warning=false)
+γ = -2 .* imag.(λ)
+s_ind = findmin(γ)[2]
+ψ = states[s_ind]
+
+# Compute the radiation pattern
+function G(r,i,j) # Green's Tensor overlap
+    G_i = GreenTensor(r-pos[i])
+    G_j = GreenTensor(r-pos[j])
+    return μ[i]' * (G_i'*G_j) * μ[j]
+end
+function intensity(r) # The intensity ⟨E⁻(r)⋅E⁺(r)⟩
+    real(sum(expect(spsm[i,j], ψ)*G(r,i,j) for i=1:N, j=1:N))
+end
+
+y = -50d:2d:50d
+z = 5d
+x = y .+ 0.5d*(N-1)
+I = zeros(length(x), length(y))
+for i=1:length(x), j=1:length(y)
+    I[i,j] = intensity([x[i],y[j],z])
+end
+
+# Plot
+figure(figsize=(9,4))
+contourf(x./d,y./d,I',30)
+for p in pos
+    plot(p[1]./d,p[2],"o",color="w",ms=2)
+end
+xlabel("x/d")
+ylabel("y/d")
+colorbar(label="Intensity",ticks=[])
+tight_layout()
+savefig("radiation-pattern.png") # hide
+```
+
+![](radiation-pattern.png)

src/CollectiveSpins.jl

@@ -5,7 +5,7 @@ export System, Spin, SpinCollection, CavityMode, CavitySpinCollection,
         interaction, GreenTensor, OmegaMatrix, GammaMatrix,
         reducedspin, ReducedSpinBasis, reducedspintransition, reducedspinstate,
                 reducedsigmap, reducedsigmam, reducedsigmax, reducedsigmay,
-                reducedsigmaz,
+                reducedsigmaz, reducedsigmapsigmam,
         collective_modes, Omega_k_chain, Gamma_k_chain, Omega_k_2D, Gamma_k_2D
 
 

src/reducedspin.jl

@@ -3,7 +3,8 @@ module reducedspin
 using QuantumOptics, Base.Cartesian
 using Combinatorics: combinations
 
-export ReducedSpinBasis, reducedspintransition, reducedspinstate, reducedsigmap, reducedsigmam, reducedsigmax, reducedsigmay, reducedsigmaz
+export ReducedSpinBasis, reducedspintransition, reducedspinstate, reducedsigmap,
+    reducedsigmam, reducedsigmax, reducedsigmay, reducedsigmaz, reducedsigmapsigmam
 
 import Base: ==
 
@@ -225,13 +226,13 @@ reducedspinstate(b::ReducedSpinBasis, n) = reducedspinstate(b, convert(Vector{In
 reducedspinstate(b::ReducedSpinBasis, n::Int) = reducedspinstate(b, [n])
 
 """
-    sigmap_sigmam(b::ReducedSpinBasis, i, j)
+    reducedsigmapsigmam(b::ReducedSpinBasis, i, j)
 
-Create the product operator product σᵢ⁺σⱼ⁻ directly on a [`ReducedSpinBasis`](@ref).
+Create the operator product σᵢ⁺σⱼ⁻ directly on a [`ReducedSpinBasis`](@ref).
 Useful especially for a basis where only one excitation is included, since then the
 single operators are zero (do not conserve excitation number), but the product is not.
 """
-function sigmap_sigmam(b::ReducedSpinBasis, i, j)
+function reducedsigmapsigmam(b::ReducedSpinBasis, i, j)
     op = SparseOperator(b)
     for idx1 in b.indexMapper
         if j in idx1[1]
@@ -260,7 +261,7 @@ function Hamiltonian(S::SpinCollection, M::Int=1)
     H = SparseOperator(b)
     for i=1:N, j=1:N
         if i != j
-            H += OmegaM[i, j]*sigmap_sigmam(b,i,j)
+            H += OmegaM[i, j]*reducedsigmapsigmam(b,i,j)
         end
     end