diff --git a/notebooks/atomic_dipole_arrays.ipynb b/notebooks/atomic_dipole_arrays.ipynb new file mode 100644 index 0000000..136abb4 --- /dev/null +++ b/notebooks/atomic_dipole_arrays.ipynb @@ -0,0 +1,265 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Excitation transfer between dipole-coupled nanorings" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example we look at the radiation pattern and transfer of collective states in rings of two-level emitters with transition frequency $\\omega_0$ following closely the work of M. Moreno-Cardoner et al, 10.1103/PhysRevA.100.023806.\n", + "\n", + "By considering only a single excitation for the system the Hamiltonian can be cast in non-Hermitian form and the dynamics follow the Schroedinger equation $\\dot{|\\Psi\\rangle} = \\frac{1}{i\\hbar} H_\\mathrm{nh} {|\\Psi\\rangle}$ with\n", + "\n", + "$$ H_\\mathrm{nh} = \\sum_{i\\neq j} (\\Omega_{ij}-\\frac{i}{2}\\Gamma_{ij}) \\sigma^-_i \\sigma^+_j$$\n", + "\n", + "For a symmetric ring, where the dipole orientations preserve the rotational symmetry (e.g. if the dipoles are\n", + "oriented perpendicularly to the plane of the ring, or tangentially to the ring), the collective modes in the\n", + "single-excitation manifold are perfect spin waves given by $|\\Psi_m\\rangle = S^+_m |g\\rangle$, with\n", + "\n", + "$$S^+_m = \\frac{1}{\\sqrt{N}} \\sum_{j=1}^N e^{i m \\phi_j} \\sigma_j^+.$$\n", + "\n", + "Here $\\phi_j = 2\\pi (j-1)/N$ is the angle associated with position $j \\ (j=1,...,N)$, and $m = 0,\\pm 1,\\pm 2,...,\\lceil \\pm (N-1)/2 \\rceil$ corresponds to the angular momentum of the mode.\n", + "\n", + "First we load the necessary libraries." + ] + }, + { + "cell_type": "code", + "execution_count": 98, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "true" + ] + }, + "execution_count": 98, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using QuantumOptics, PyPlot, CollectiveSpins;\n", + "pygui(true)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we define the parameters, operators, basis and the geometry for two rings with 10 emitters each. Here we choose a tangential polarization for the dipole transitions." + ] + }, + { + "cell_type": "code", + "execution_count": 99, + "metadata": {}, + "outputs": [], + "source": [ + "N1 = 10 # number of emitters in first ring\n", + "Γ0 = 1.0 # decay rate\n", + "d = 0.1 # dipole distance\n", + "k0 = 2pi # corresponds to λ=1\n", + "\n", + "# Ring Geometry\n", + "ϕ(j) = 2π/N1*(j-1) # Angle separating atoms\n", + "Rot(φ) = [cos(φ) -sin(φ) 0;sin(φ) cos(φ) 0; 0 0 1] # Tangential polarization\n", + "R = 0.5d/sin((ϕ(2)/2))\n", + "pos = [[R*cos(ϕ(j)), R*sin(ϕ(j)), 0.0] for j=1:N1];\n", + "dips = [normalize(Rot(pi/2)*pos[i]) for i=1:N1];\n", + "\n", + "# add second ring of same size\n", + "for i = 1:N1\n", + " push!(pos,pos[i]+[0.0,2R+2d,0.0])\n", + " push!(dips,normalize(Rot(-pi/2)*pos[i]))\n", + "end\n", + "\n", + "# create N-Level basis for the whole system\n", + "N2 = length(pos)\n", + "b = NLevelBasis(N2)\n", + "sm(i) = transition(b, N2, i)\n", + "sp(i) = transition(b, i, N2)\n", + "σ = sm.([1:N2;])\n", + "σp = sp.([1:N2;]);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The collective dipole and decay matrix elements are related to the Green's tensor $G(\\vec{r}-\\vec{r}_i,\\omega_0)$ in free space and can be generated via the CollectiveSpins.jl package. In this example we start with the most subradiant state of the first ring corresponding to $m = \\lceil \\pm (N-1)/2 \\rceil$ and let the state evolve in time.\n", + "\n", + "We also calculate the output field $E^+(\\vec{r}) = \\frac{|\\vec{\\mu}|k_0^2}{\\epsilon_0}\\sum_i G(\\vec{r}-\\vec{r}_i,\\omega_0)\\cdot \\vec{\\mu}_i\\sigma^-_i$ of the system and then the intensity $I(\\vec{r}) = \\langle E^+ E^- \\rangle$ at a point $\\vec{r}$, which is generated by all the dipoles." + ] + }, + { + "cell_type": "code", + "execution_count": 100, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "I_exp (generic function with 1 method)" + ] + }, + "execution_count": 100, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "G(i,j) = transpose(dips[i]) * GreenTensor(pos[i] - pos[j], k0) * dips[j]\n", + "\n", + "function Ω(i,j)\n", + " if i==j\n", + " return 0.0\n", + " else\n", + " return -3pi/k0*0.5*real(G(i,j))\n", + " end\n", + "end\n", + "\n", + "function Γ(i,j)\n", + " if i==j\n", + " return 1.0\n", + " else\n", + " return 6pi/k0*0.5*imag(G(i,j))\n", + " end\n", + "end\n", + "\n", + "# Eigenstates for the 1. ring\n", + "H_nh = sum((Ω(i,j)-1im*Γ(i,j)*0.5)*sp(i)*sm(j) for i=1:N1,j=1:N1) \n", + "λ = eigenstates(dense(H_nh);warning = false)\n", + "psi0 = λ[2][end] # initial state, most subradiant\n", + "\n", + "\n", + "# Schroedinger time evolution for both rings\n", + "tspan = [0.0:0.2:100;]\n", + "Heff = sum((Ω(i,j)-1im*Γ(i,j)*0.5)*sp(i)*sm(j) for i=1:N2,j=1:N2)\n", + "tout,psit = timeevolution.schroedinger(tspan,psi0,H_nh)\n", + "\n", + "function Ep(r,N) # electric field\n", + " out = [0*transition(NLevelBasis(N), N, N) for i=1:3]\n", + " for i=1:N\n", + " G_ = GreenTensor(r - pos[i], k0)\n", + " Gpi = (G_ * dips[i])\n", + " @inbounds for j=1:3\n", + " out[j].data .+= sm(i).data .* Gpi[j]\n", + " end\n", + " end\n", + " out\n", + "end\n", + "\n", + "function Intensity(r,N)\n", + " E = Ep(r,N)\n", + " sum(dagger(e)*e for e=E)\n", + "end\n", + "\n", + "I_exp(r, ψ) = real(expect(Intensity(r,N2), ψ))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally a grid in real space for the radiation field is defined and the instances of time are chosen such that the excitation is either almost completely in ring one or two." + ] + }, + { + "cell_type": "code", + "execution_count": 101, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Define grid\n", + "xmin = -0.75+R ; xmax = 0.75-R;\n", + "x = collect(range(-0.5, stop=0.5, length=61))\n", + "y = collect(range(-0.5, stop=2R+0.5, length=61))\n", + "z = 1.5R\n", + "grid = [[x1, y1, z1] for x1=x, y1=y, z1=z];\n", + "\n", + "figure(figsize=(9, 3))\n", + "kspan = [1,81,181]\n", + "for k = 1:3\n", + " t = kspan[k];t2 = tout[t]\n", + " I_profile = [I_exp(gr, psit[t]) for gr=grid];\n", + " subplot(1,3,k)\n", + " title(\"t/Γ₀ = $t2\")\n", + " contourf(y, x, I_profile, 100,cmap=\"viridis\")\n", + " for i = 1:N2\n", + " plot(pos[i][2],pos[i][1],\"wo\",markersize=6,alpha=0.5)\n", + " end\n", + " xlabel(\"x\")\n", + " ylabel(\"y\")\n", + " tight_layout()\n", + "end\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To look more closely at the transport behaviour of the two rings we define a projector which counts only excitations in the second ring and look at its time evolution for the 5 most subradiant states starting in the first ring." + ] + }, + { + "cell_type": "code", + "execution_count": 102, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "function projector(psi)\n", + " acc = 0.0\n", + " for i = N1+1:N2\n", + " acc += expect(sp(i)*sm(i),psi)\n", + " end\n", + " return acc\n", + "end\n", + "\n", + "figure(figsize=(7, 4))\n", + "for i = 1:4\n", + " psi0 = λ[2][end+1-i]\n", + " tout,psit = timeevolution.schroedinger(tspan,psi0,Heff)\n", + " subplot(2,2,i)\n", + " plot(tout,projector.(psit))\n", + " xlabel(\"Γ₀t\")\n", + " ylabel(\"2.Ring population\")\n", + " tight_layout()\n", + " ylim(0,1)\n", + "end\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.6.1", + "language": "julia", + "name": "julia-1.6" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}