From c3b3f571c98af98302f824144e6d7c75557d338f Mon Sep 17 00:00:00 2001 From: Matias Bundgaard-Nielsen <56378038+mabuni1998@users.noreply.github.com> Date: Wed, 5 Jun 2024 11:36:10 +0200 Subject: [PATCH 1/2] Normalization update Waveguide states are no longer by default initialized as normalized states. This is to ensure compatibility with creating superpositions / beamsplitter operations. In relation to this, the twophoton initializer was adjusted to so that twophotons in the same waveguide now requires a prefactor of 1/sqrt(2) in order to be normalized. Documentation was also updated. Precompilation updated to be faster. --- Manifest.toml | 46 ++++++++------ Project.toml | 3 +- docs/Manifest.toml | 4 +- docs/src/beamsplitter.md | 53 +++++++++++++--- docs/src/theoreticalbackground.md | 3 +- docs/src/tutorial.md | 12 ++-- src/WaveguideQED.jl | 1 + src/basis.jl | 100 +++++++++++++++--------------- src/precompile.jl | 9 ++- src/solver.jl | 8 ++- 10 files changed, 147 insertions(+), 92 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index f2ca72d..e36ea85 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,8 +1,8 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.9.3" +julia_version = "1.10.0" manifest_format = "2.0" -project_hash = "1c494bef7d7f4b16e46d8f7d4eae31049f5ce1bf" +project_hash = "436ff8abc8c4ebed1fdacfb2505a5c6812e54643" [[deps.ADTypes]] git-tree-sha1 = "e58c18d2312749847a74f5be80bb0fa53da102bd" @@ -139,7 +139,7 @@ weakdeps = ["Dates", "LinearAlgebra"] [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.0.5+0" +version = "1.0.5+1" [[deps.ConstructionBase]] deps = ["LinearAlgebra"] @@ -505,21 +505,26 @@ version = "1.0.0" [[deps.LibCURL]] deps = ["LibCURL_jll", "MozillaCACerts_jll"] uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" -version = "0.6.3" +version = "0.6.4" [[deps.LibCURL_jll]] deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" -version = "7.84.0+0" +version = "8.4.0+0" [[deps.LibGit2]] -deps = ["Base64", "NetworkOptions", "Printf", "SHA"] +deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" +[[deps.LibGit2_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] +uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" +version = "1.6.4+0" + [[deps.LibSSH2_jll]] deps = ["Artifacts", "Libdl", "MbedTLS_jll"] uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" -version = "1.10.2+0" +version = "1.11.0+1" [[deps.Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" @@ -620,7 +625,7 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" [[deps.MbedTLS_jll]] deps = ["Artifacts", "Libdl"] uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.28.2+0" +version = "2.28.2+1" [[deps.Missings]] deps = ["DataAPI"] @@ -633,7 +638,7 @@ uuid = "a63ad114-7e13-5084-954f-fe012c677804" [[deps.MozillaCACerts_jll]] uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2022.10.11" +version = "2023.1.10" [[deps.MuladdMacro]] git-tree-sha1 = "cac9cc5499c25554cba55cd3c30543cff5ca4fab" @@ -677,12 +682,12 @@ version = "1.12.9" [[deps.OpenBLAS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.21+4" +version = "0.3.23+2" [[deps.OpenLibm_jll]] deps = ["Artifacts", "Libdl"] uuid = "05823500-19ac-5b8b-9628-191a04bc5112" -version = "0.8.1+0" +version = "0.8.1+2" [[deps.OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] @@ -722,7 +727,7 @@ version = "0.12.3" [[deps.Pkg]] deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.9.2" +version = "1.10.0" [[deps.PoissonRandom]] deps = ["Random"] @@ -811,7 +816,7 @@ deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" [[deps.Random]] -deps = ["SHA", "Serialization"] +deps = ["SHA"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [[deps.Random123]] @@ -981,6 +986,7 @@ version = "1.1.1" [[deps.SparseArrays]] deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +version = "1.10.0" [[deps.SparseDiffTools]] deps = ["ADTypes", "Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "Graphs", "LinearAlgebra", "Reexport", "Requires", "SciMLOperators", "Setfield", "SparseArrays", "StaticArrayInterface", "StaticArrays", "Tricks", "VertexSafeGraphs"] @@ -1041,7 +1047,7 @@ version = "1.4.0" [[deps.Statistics]] deps = ["LinearAlgebra", "SparseArrays"] uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -version = "1.9.0" +version = "1.10.0" [[deps.StatsAPI]] deps = ["LinearAlgebra"] @@ -1098,9 +1104,9 @@ deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" [[deps.SuiteSparse_jll]] -deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"] +deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" -version = "5.10.1+6" +version = "7.2.1+1" [[deps.SymbolicIndexingInterface]] deps = ["DocStringExtensions"] @@ -1206,7 +1212,7 @@ version = "2.0.0" [[deps.Zlib_jll]] deps = ["Libdl"] uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.13+0" +version = "1.2.13+1" [[deps.ZygoteRules]] deps = ["ChainRulesCore", "MacroTools"] @@ -1217,14 +1223,14 @@ version = "0.2.3" [[deps.libblastrampoline_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" -version = "5.8.0+0" +version = "5.8.0+1" [[deps.nghttp2_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" -version = "1.48.0+0" +version = "1.52.0+1" [[deps.p7zip_jll]] deps = ["Artifacts", "Libdl"] uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" -version = "17.4.0+0" +version = "17.4.0+2" diff --git a/Project.toml b/Project.toml index 47b052d..e34e5df 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "WaveguideQED" uuid = "c4555495-0e8d-488d-8e6a-965ecefe9055" authors = ["Matias Bundgaard-Nielsen , Mikkel Heuck , and Stefan Krastanov "] -version = "0.2.4" +version = "0.2.5" [deps] FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" @@ -12,6 +12,7 @@ QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c" QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6" [compat] diff --git a/docs/Manifest.toml b/docs/Manifest.toml index 663368c..56027be 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -1425,10 +1425,10 @@ uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f" version = "0.2.0" [[deps.WaveguideQED]] -deps = ["FFTW", "LinearAlgebra", "Parameters", "PrecompileTools", "QuantumOptics", "QuantumOpticsBase", "SparseArrays", "Strided", "UnsafeArrays"] +deps = ["FFTW", "LinearAlgebra", "Parameters", "PrecompileTools", "QuantumOptics", "QuantumOpticsBase", "SparseArrays", "Strided", "Test", "UnsafeArrays"] path = ".." uuid = "c4555495-0e8d-488d-8e6a-965ecefe9055" -version = "0.2.4" +version = "0.2.5" [[deps.WignerSymbols]] deps = ["HalfIntegers", "LRUCache", "Primes", "RationalRoots"] diff --git a/docs/src/beamsplitter.md b/docs/src/beamsplitter.md index 03e782d..ac029d6 100644 --- a/docs/src/beamsplitter.md +++ b/docs/src/beamsplitter.md @@ -1,6 +1,6 @@ # [Beamsplitters](@id Beamsplitter) -Having introduced multiple waveguides in [Multiple Waveguides](@id multiple), it is natural to implement the beamsplitter operation and consider some of the common measurements one would make on states having undergone a beamsplitter transformation. +Having introduced multiple waveguides in [`Multiple Waveguides`](@ref multiple), it is natural to implement the beamsplitter operation and consider some of the common measurements one would make on states having undergone a beamsplitter transformation. We start by introducing how a beamsplitter can be implemented using two waveguides that interact. For simplicity, we start by considering only a single photon in one waveguide. First we create the basis, operators of the multiple waveguides and an initial single photon state with a gaussian wavefunction residing in waveguide 1: @@ -23,7 +23,7 @@ psi = onephoton(bw,1, ξ) nothing #hide ``` -We now want to consider the following beamsplitter transformation: $w_{k,1} \rightarrow \cos(V) w_{k,1} - i \sin(V) w_{k,2}$ and equivalently for waveguide 2: $w_{k,2} \rightarrow - i \sin(V) w_{k,2} + \cos(V) w_{k,1}$. Having access to our initial Gaussian wavefunction we could just create the transformed state as: +We now want to consider the following beamsplitter transformation: $w_{k,1} \rightarrow \cos(V) w_{k,1} - i \sin(V) w_{k,2}$ and equivalently for waveguide 2: $w_{k,2} \rightarrow - i \sin(V) w_{k,1} + \cos(V) w_{k,2}$. Having access to our initial Gaussian wavefunction we could just create the transformed state as: ```@example bs V = pi/4 @@ -34,6 +34,11 @@ nothing #hide A more automatic and equivalent method is, however, instead to let the waveguide state undergo evolution under the Hamiltonian: $H = V( w_{k,1}^\dagger w_{k,2} + w_{k,2}^\dagger w_{k,1})$, which performs the same transformation. See [Section 4.3](https://github.com/qojulia/WaveguideQED.jl/blob/main/Thesis/Master_s_thesis__Modeling_Tools_For_Quantum_Networks%20(9).pdf) for details of the derivation. We can confirm this by: ```@example bs +using PyPlot +rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams") #hide +rcParams["font.size"] = 20 #hide +rcParams["font.family"] = "serif" #hide +rcParams["mathtext.fontset"] ="cm" #hide Vs = 0:pi/32:pi reflection = zeros(length(Vs)) transmission = zeros(length(Vs)) @@ -48,7 +53,6 @@ for (i,V) in enumerate(Vs) transmission[i] = expect_waveguide(n1,psi_trans) reflection[i] = expect_waveguide(n2,psi_trans) end -using PyPlot #hide fig,ax = subplots(1,1,figsize=(9,4.5)) ax.plot(Vs/pi,reflection,"b-",label="Waveguide a") ax.plot(Vs/pi,transmission,"r-",label="Waveguide b") @@ -64,13 +68,13 @@ nothing #hide Here we see this population of waveguide 1 and 2 after the transformation vary as cosines and sines as we change the interaction parameter V. Thus, we confirm that we are applying the desired transformation. For an even beamsplitter, we thus choose $V=\pi/4$ ## [Hong Ou Mandel with twophotons](@id hom) -As a more advanced example, we now consider a Hong Ou Mandel setup, where we have one photon in each waveguide impinging on a beamsplitter. If the two photons equivalent, we will see the Hong Ou Mandel effect and thus expect no photons in both waveguide simultanouesly after the transformation. As a measure of this, we calcuate the chance of having a coincedence count where one photon is in waveguide 1 while the other is in waveguide 2. This calculated using the two projection operators: +As a more advanced example, we now consider a Hong Ou Mandel setup, where we have one photon in each waveguide impinging on a beamsplitter. If the two photons equivalent, we will see the Hong Ou Mandel effect and thus expect no photons in both waveguide simultanouesly after the transformation. As a measure of this, we calcuate the chance of having a coincidences count where one photon is in waveguide 1 while the other is in waveguide 2. This calculated using the two projection operators: $$P_1 = \int_0^T dt w_1^\dagger(t) |0\rangle\langle0| w_1(t) \qquad P_2 = \int_0^T dt w_2^\dagger(t) |0\rangle\langle0| w_2(t)$$ where $w_1(t)$ and $w_2(t)$ are lowering operators for two waveguides. The chance of coincidence count is computed by $\langle\psi|P_1 P_2 |\psi\rangle$. -To compute the councidence count expectation we create our own custom expectation value function: +To compute the coincidence count expectation we create our own custom expectation value function: ```@example bs n1 = wd1*w1 @@ -96,21 +100,21 @@ nothing #hide Here we evaluate $w_1^\dagger(t) w_1(t)$ at one timeidex `i`, while we evaluate $w_2^\dagger(t) w_2(t)$ at another timeidex `j`. Together this gives us the total coincedence count chance. In the following, we use the Hamiltonian from the previous section with $V=\pi/4$ and consider two Gaussian photons in each their waveguide with different centers of time $t_0$. By changing the difference in $t_0$, we can see the transition from a perfect overlap meaning no coincedence count, to no overlap meaning that the two photons never interact. In this case, the two photons will split up randomly and $1/4$ of the time they will end up in waveguide 1, similarly $1/4$ of the time they will end up in waveguide 2, and the remaining $1/2$ time they will end up in each of their waveguides. Thus, we expect a coincedence count of $1/2$ when the two pulses are fully seperated. Note that in the above function, we can just use the waveguide operators as projectors as we never have twophotons in both waveguides. ```@example bs -taus = 0:0.2:4 +taus = -3:0.1:3 ξ_twophoton(t1, t2, t01, t02) = ξ(t1, t01) * ξ(t2, t02) V = pi/4 H = V/dt*(wd1*w2 + wd2*w1) -coincedences = zeros(length(taus)) +coincidences = zeros(length(taus)) t01 = 5 for (i, τ) in enumerate(taus) t02 = 5 + τ psi_pre = twophoton(bw, [1,2], ξ_twophoton, t01,t02) ψ = waveguide_evolution(times,psi_pre,H) - coincedences[i] = expect_waveguide2(expval_op, ψ,times) + coincidences[i] = expect_waveguide2(expval_op, ψ,times) end -fig, ax = subplots(1,1, figsize=(6,4)) -ax.plot(taus, coincedences) +fig, ax = subplots(1,1, figsize=(9,4.5)) +ax.plot(taus, coincidences,"r-") ax.set_xlabel(L"Delay between pules $\tau$") ax.set_ylabel("Coincedence chance") tight_layout() @@ -118,3 +122,32 @@ savefig("hom.svg") #hide nothing #hide ``` ![hom_plot](hom.svg) + +We could also have created this plot by performing the the beamsplitter operation by hand and instead initializing the state directly in this state. The initial state we consider is a single photon in each waveguide: $$|\psi \rangle = \int \int dt_1 dt_2 \xi_1(t_1) \xi_2(t_2) w_1^\dagger(t_1) w_2^\dagger(t_2) \ket{\emptyset}$$, where $$\xi_1(t_1)$$ and $$\xi_2(t_2)$$ denote the wavefunction of the photon in waveguide 1 and 2, respectively. Notice that there is not factor of $$1/\sqrt(2)$$ in front of the initial state as the two photons occupy each their waveguide. If they initially occupied the same waveguide, we would need a factor of $$1/\sqrt(2)$$ for the state to be normalized. Performing the beamsplitter operation $w_1(t) \rightarrow 1/\sqrt(2) ( w_1(t) - i w_2(t))$ and $w_2(t) \rightarrow 1/\sqrt(2) ( - i w_1(t) + w_2(t))$, we arrive at the transformed state: + +$$\begin{equation*} +|\psi \rangle_{BS} = \frac{1}{2}\left ( i \int \int dt_1 dt_2 \xi_1(t_1) \xi_2(t_2) w_1^\dagger(t_1) w_1^\dagger(t_2) \ket{\emptyset} - i \int \int dt_1 dt_2 \xi_1(t_1) \xi_2(t_2) w_2^\dagger(t_1) w_2^\dagger(t_2) \ket{\emptyset} + \int \int dt_1 dt_2 \xi_1(t_1) \xi_2(t_2) w_1^\dagger(t_1) w_2^\dagger(t_2)\ket{\emptyset} \right - \int \int dt_1 dt_2 \xi_1(t_1) \xi_2(t_2) w_2^\dagger(t_1) w_1^\dagger(t_2)\ket{\emptyset} \right )\end{equation*}$$ + + +```@example bs +taus = -3:0.1:3 +coincidences_manual = zeros(length(taus)) +t01 = 5 +for (i, τ) in enumerate(taus) + t02 = 5 + τ + psi_trans = 1/2*( im*twophoton(bw, 1, ξ_twophoton, t01,t02) - im * twophoton(bw, 2, ξ_twophoton, t01,t02) + + twophoton(bw, [1,2], ξ_twophoton, t01,t02) - twophoton(bw, [2,1], ξ_twophoton, t01,t02)) + coincidences_manual[i] = expect_waveguide2(expval_op, psi_trans,times) +end + +fig, ax = subplots(1,1, figsize=(9,4.5)) +ax.plot(taus, coincidences,"r-",label="Hamiltonian transformation") +ax.plot(taus, coincidences_manual,"b--",label="Manual transformation") +ax.set_xlabel(L"Delay between pules $\tau$") +ax.set_ylabel("Coincedence chance") +ax.legend(fontsize=10) +tight_layout() +savefig("hom2.svg") #hide +nothing #hide +``` +![hom_plot2](hom2.svg) diff --git a/docs/src/theoreticalbackground.md b/docs/src/theoreticalbackground.md index 105cbc2..3314120 100644 --- a/docs/src/theoreticalbackground.md +++ b/docs/src/theoreticalbackground.md @@ -128,8 +128,7 @@ $$\begin{align*} \frac{1}{\sqrt{2}}\left[W^\dagger(\xi)\right]^2|0\rangle &= \frac{1}{\sqrt{2}} \int_{t_0}^{t_{end}} d t^{\prime} \int_{t_0}^{t_{end}} d t \ \xi(t) \xi\left(t^{\prime}\right) w^\dagger(t) w^\dagger\left(t^{\prime}\right)|0\rangle \\ & \rightarrow \frac{1}{\sqrt{2}} \sum_{i=1}^N \sum_{k=1}^N \xi\left(t_i\right) \xi\left(t_k\right) w^\dagger\left(t_i\right) w^{\dagger}\left(t_k\right)|0\rangle \\ & =\frac{1}{\sqrt{2}} \sum_{i=1}^N \sum_{k \neq i}^N \xi\left(t_i\right) \xi\left(t_k\right) w^{\dagger}\left(t_i\right) w^{\dagger}\left(t_k\right)|0\rangle+\sum_{i=1}^N \xi\left(t_i\right) \xi\left(t_i\right)\left|2 t_i\right\rangle \\ -& =\frac{2}{\sqrt{2}} \sum_{i=1}^N \sum_{k>i}^N \xi\left(t_i\right) \xi\left(t_k\right)\left|1_{t_i} 1_{t_k}\right\rangle+\sum_{i=1}^N \xi\left(t_i\right) \xi\left(t_i\right)\left|2 t_i\right\rangle \\ -& =\sqrt{2} \sum_{i=1}^N \sum_{k > i}^N \xi\left(t_i\right) \xi\left(t_k\right) \mid 1_{t_i} 1_{t_k}\rangle + \sum_{i=1}^N \xi\left(t_i\right) \xi\left(t_i\right)\left|2 t_i\right\rangle +& =\frac{1}{\sqrt{2}} \sum_{i=1}^N \sum_{k>i}^N (\xi\left(t_i\right) \xi\left(t_k\right) + \xi\left(t_k\right) \xi\left(t_i\right)) \left|1_{t_i} 1_{t_k}\right\rangle+\sum_{i=1}^N \xi\left(t_i\right) \xi\left(t_i\right)\left|2 t_i\right\rangle \end{align*}$$ diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 74060cb..b746e00 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -40,7 +40,9 @@ With this, we can now simulate the scattering of a single photon with a Gaussian nothing #hide ``` -Assuming the cavity is empty, the combined initial state is then: +Note that the wavefunction $$\xi(t)$$ is assumed to be normalized: $$\int_0^\infty \xi(t) dt = 1$$ when creating waveguide states. If one wants to use a non-normalized function, the keyword `norm=true` can be passed to [`onephoton`](@ref) to ensure normalization. + +Considering an initially empty cavity, the combined initial state is then: ```@example tutorial ψ_in = fockstate(bc,0) ⊗ ψ_waveguide @@ -148,17 +150,19 @@ H_twophoton = im*sqrt(γ/dt)*( ad ⊗ w_twophoton - a ⊗ wd_twophoton ) nothing #hide ``` -If we want an initial two-photon state, we instead use the function [`twophoton`](@ref) to create a two-photon state $\frac{1}{\sqrt{2}}\left[W^\dagger(\xi)\right]^2|0\rangle = \frac{1}{\sqrt{2}} \int_{t_0}^{t_{end}} d t^{\prime} \int_{t_0}^{t_{end}} d t \ \xi^{(2)}(t,t') w^\dagger(t) w^\dagger\left(t^{\prime}\right)|0\rangle $ (see [Theoretical Background](@ref theory) for details). In the following, we define the two-photon wavefunction $\xi^{(2)}(t,t') = \xi^{(1)}(t)\xi^{(1)}(t')$ which is thus a product state of two single-photons. +If we want an initial two-photon state, we instead use the function [`twophoton`](@ref) to create a two-photon state $\frac{1}{\sqrt{2}}\left[W^\dagger(\xi)\right]^2|0\rangle = \frac{1}{\sqrt{2}} \int_{t_0}^{t_{end}} d t^{\prime} \int_{t_0}^{t_{end}} d t \ \xi^{(2)}(t,t') w^\dagger(t) w^\dagger\left(t^{\prime}\right)|0\rangle $ (see [`Theoretical Background`](@ref theory) for details). In the following, we define the two-photon wavefunction $\xi^{(2)}(t,t') = \xi^{(1)}(t)\xi^{(1)}(t')$ which is thus a product state of two single-photons. ```@example tutorial ξ2(t1,t2,σ,t0) = ξ(t1,σ,t0)*ξ(t2,σ,t0) σ,t0 = 1,5 -ψ_twophoton = twophoton(bw,ξ2,σ,t0) +ψ_twophoton = 1/sqrt(2)*twophoton(bw,ξ2,σ,t0) ψ_in_twophoton = fockstate(bc,0) ⊗ ψ_twophoton nothing #hide ``` -Notice the structure of `ξ2(t1,t2,σ,t0)`, it now has two time-arguments and the remaining arguments are parameters. If we wanted to allow for two different widths of the single-photon states in the product state, we could have also defined: `ξ2(t1,t2,σ1,σ2,t0) = ξ(t1,σ1,t0)*ξ(t2,σ2,t0)`. In the following, we consider the more simple case of equivalent photons. We solve the two-photon scattering in the following. +Notice the structure of `ξ2(t1,t2,σ,t0)`, it now has two time-arguments and the remaining arguments are parameters. If we wanted to allow for two different widths of the single-photon states in the product state, we could have also defined: `ξ2(t1,t2,σ1,σ2,t0) = ξ(t1,σ1,t0)*ξ(t2,σ2,t0)`. Another important detail is the normalization. [`twophoton`](@ref) only creates $$\int_{t_0}^{t_{end}} d t^{\prime} \int_{t_0}^{t_{end}} d t \ \xi^{(2)}(t,t') w^\dagger(t) w^\dagger\left(t^{\prime}\right)|0\rangle $$ and we thus need the factor of $$1/\sqrt(2)$$ for the state to be normalized. + +In the following, we consider the more simple case of equivalent photons. We solve the two-photon scattering in the following. ```@example tutorial ψ_out = waveguide_evolution(times,ψ_in_twophoton,H_twophoton) diff --git a/src/WaveguideQED.jl b/src/WaveguideQED.jl index ccb9631..47f2f6c 100644 --- a/src/WaveguideQED.jl +++ b/src/WaveguideQED.jl @@ -4,6 +4,7 @@ using QuantumOptics using Strided using UnsafeArrays using FFTW +using Test import LinearAlgebra: axpy!, dot, mul!, rmul!,I import QuantumOpticsBase: create, dagger, destroy, expect, identityoperator, tensor,set_time! diff --git a/src/basis.jl b/src/basis.jl index 9317dc2..6bda58f 100644 --- a/src/basis.jl +++ b/src/basis.jl @@ -210,14 +210,14 @@ OnePhotonView(ψ,2) == ones(length(times)) function onephoton end """ - onephoton(b::WaveguideBasis{T,1},ξ::Function,args...,norm=True) where {T} + onephoton(b::WaveguideBasis{T,1},ξ::Function,args...,norm=false) where {T} * Since `b` only contains a single waveguide, the index of the waveguide is not needed. * `ξ` should be broadcastable as ξ.(times,args...), where `times = 0:b.dt:(b.nsteps-1)*b.dt` (the times used to define the waveguidebasis) * `args...`: additional arguments to be passed to ξ if it is a function. -* `norm::Bool=true`: If true normalize the resulting wavepacket. +* `norm::Bool=false`: If true normalize the resulting wavepacket. """ -function onephoton(b::WaveguideBasis{T,1},ξ::Function,args...;norm=true) where {T} +function onephoton(b::WaveguideBasis{T,1},ξ::Function,args...;norm=false) where {T} state = Ket(b) view = OnePhotonView(state) view .= ξ.(0:b.dt:(b.nsteps-1)*b.dt,args...)*sqrt(b.dt) @@ -228,13 +228,13 @@ function onephoton(b::WaveguideBasis{T,1},ξ::Function,args...;norm=true) where end """ - onephoton(b::WaveguideBasis{T,1},ξ::AbstractArray;norm=true) where {T} + onephoton(b::WaveguideBasis{T,1},ξ::AbstractArray;norm=false) where {T} * Since `b` only contains a single waveguide, the index of the waveguide is not needed. * `ξ` should be `AbstractArray` with length(ξ) == b.nsteps -* `norm::Bool=true`: If true normalize the resulting wavepacket. +* `norm::Bool=false`: If true normalize the resulting wavepacket. """ -function onephoton(b::WaveguideBasis{T,1},ξ::AbstractArray;norm=true) where {T} +function onephoton(b::WaveguideBasis{T,1},ξ::AbstractArray;norm=false) where {T} state = Ket(b) view = OnePhotonView(state) view .= ξ*sqrt(b.dt) @@ -245,15 +245,15 @@ function onephoton(b::WaveguideBasis{T,1},ξ::AbstractArray;norm=true) where {T} end """ - onephoton(b::WaveguideBasis{T,Nw},i::Int,ξ::Function,args...; norm=true) where {T,Nw} + onephoton(b::WaveguideBasis{T,Nw},i::Int,ξ::Function,args...; norm=false) where {T,Nw} * Since `b` contains `Nw` waveguides, the index of the waveguide needed. * `i` is the index of the waveguide at which the onephoton wavepacket is created in * `ξ` should be broadcastable as ξ.(times,args...), where `times = 0:b.dt:(b.nsteps-1)*b.dt` (the times used to define the waveguidebasis) * `args...`: additional arguments to be passed to ξ if it is a function. -* `norm::Bool=true`: If true normalize the resulting wavepacket. +* `norm::Bool=false`: If true normalize the resulting wavepacket. """ -function onephoton(b::WaveguideBasis{T,Nw},i::Int,ξ::Function,args...; norm=true) where {T,Nw} +function onephoton(b::WaveguideBasis{T,Nw},i::Int,ξ::Function,args...; norm=false) where {T,Nw} state = Ket(b) view = OnePhotonView(state,i) view .= ξ.(0:b.dt:(b.nsteps-1)*b.dt,args...)*sqrt(b.dt) @@ -264,14 +264,14 @@ function onephoton(b::WaveguideBasis{T,Nw},i::Int,ξ::Function,args...; norm=tru end """ - onephoton(b::WaveguideBasis{T,Nw},i,ξvec;norm=true) where {T,Nw} + onephoton(b::WaveguideBasis{T,Nw},i,ξvec;norm=false) where {T,Nw} * Since `b` contains a `Nw` waveguides, the index of the waveguide needed. * `i` is the index of the waveguide at which the onephoton wavepacket is created in * `ξ` should be `AbstractArray` with length(ξ) == b.nsteps -* `norm::Bool=true`: If true normalize the resulting wavepacket. +* `norm::Bool=false`: If true normalize the resulting wavepacket. """ -function onephoton(b::WaveguideBasis{T,Nw},i::Int,ξ::AbstractArray;norm=true) where {T,Nw} +function onephoton(b::WaveguideBasis{T,Nw},i::Int,ξ::AbstractArray;norm=false) where {T,Nw} state = Ket(b) view = OnePhotonView(state,i) view .= ξ*sqrt(b.dt) @@ -280,8 +280,8 @@ function onephoton(b::WaveguideBasis{T,Nw},i::Int,ξ::AbstractArray;norm=true) w end return state end -onephoton(b::WaveguideBasis{T,Nw},ξ::Function,args...;norm=true) where {T,Nw} = throw(ArgumentError("WaveguideBasis contains multiple waveguides. Please provide the index of the waveguide in which the excitation should be created.")) -onephoton(b::WaveguideBasis{T,Nw},ξ::AbstractArray;norm=true) where {T,Nw} = throw(ArgumentError("WaveguideBasis contains multiple waveguides. Please provide the index of the waveguide in which the excitation should be created.")) +onephoton(b::WaveguideBasis{T,Nw},ξ::Function,args...;norm=false) where {T,Nw} = throw(ArgumentError("WaveguideBasis contains multiple waveguides. Please provide the index of the waveguide in which the excitation should be created.")) +onephoton(b::WaveguideBasis{T,Nw},ξ::AbstractArray;norm=false) where {T,Nw} = throw(ArgumentError("WaveguideBasis contains multiple waveguides. Please provide the index of the waveguide in which the excitation should be created.")) @@ -569,16 +569,16 @@ bw = WaveguideBasis(2,3,times) function twophoton end """ - twophoton(b::WaveguideBasis{T,1},ξ::Function,args...,norm=True) where {T} + twophoton(b::WaveguideBasis{T,1},ξ::Function,args...,norm=false) where {T} # Arguments - `b::WaveguideBasis{T,Nw}` contains only one waveguide and no index needed. - ξ given as a function should follow `ξ(t1,t2,args...)`, where `t1` and `t2` are the input times. - `args...`: additional arguments to be passed to ξ if it is a function. -- `norm::Bool=true`: normalize the resulting wavepacket. +- `norm::Bool=false`: Toggles whether to normalize the resulting wavepacket. """ -function twophoton(b::WaveguideBasis{T,1},ξ::Function,args...;norm=true) where {T} +function twophoton(b::WaveguideBasis{T,1},ξ::Function,args...;norm=false) where {T} state = Ket(b) nsteps = get_nsteps(b) viewed_data = TwoPhotonView(state.data,nsteps,nsteps+1) @@ -586,11 +586,11 @@ function twophoton(b::WaveguideBasis{T,1},ξ::Function,args...;norm=true) where for l in 1:nsteps for m in 1:nsteps if l Date: Wed, 5 Jun 2024 12:01:55 +0200 Subject: [PATCH 2/2] Test fixed after normalization update Detection documentation is now no longer included as this is deprecated functions. --- test/helper_functions.jl | 8 ++++---- test/runtests.jl | 2 -- test/test_operators.jl | 8 ++++---- 3 files changed, 8 insertions(+), 10 deletions(-) diff --git a/test/helper_functions.jl b/test/helper_functions.jl index d2a8a8f..78c19f6 100644 --- a/test/helper_functions.jl +++ b/test/helper_functions.jl @@ -138,7 +138,7 @@ function test_multiplewaveguides(b,Nw,idx,order,wda1,adw1) nsteps = length(times) dt = times[2] - times[1] ξfun(t1) = 1 - psi = tensor([i == bw_idx ? onephoton(b.bases[i],idx,ξfun) : i==bc_idx ? fockstate(b.bases[i],1) : fockstate(b.bases[i],0) for i in 1:3]...) + psi = tensor([i == bw_idx ? onephoton(b.bases[i],idx,ξfun,norm=true) : i==bc_idx ? fockstate(b.bases[i],1) : fockstate(b.bases[i],0) for i in 1:3]...) tmp = copy(psi) tmp2 = copy(psi) @@ -172,7 +172,7 @@ function test_multiplewaveguides(b,Nw,idx,order,wda1,adw1) psi_view = view_waveguide(tmp,ground_idx) for k in filter(x -> x != idx, 1:Nw) - psi = tensor([i == bw_idx ? onephoton(b.bases[i],k,ξfun) : i==bc_idx ? fockstate(b.bases[i],1) : fockstate(b.bases[i],0) for i in 1:3]...) + psi = tensor([i == bw_idx ? onephoton(b.bases[i],k,ξfun,norm=true) : i==bc_idx ? fockstate(b.bases[i],1) : fockstate(b.bases[i],0) for i in 1:3]...) mul!(tmp,wda1,psi,1,0) i,j = min(k,idx),max(k,idx) index = (i-1)*Nw + j - (i*(i+1))÷2 @@ -199,7 +199,7 @@ function test_interaction(b,Nw,idx1,idx2,order,wd2w1) nsteps = length(times) dt = times[2] - times[1] ξfun(t1) = 1 - psi = tensor([i == bw_idx ? onephoton(b.bases[i],idx1,ξfun) : fockstate(b.bases[i],0) for i in 1:3]...) + psi = tensor([i == bw_idx ? onephoton(b.bases[i],idx1,ξfun,norm=true) : fockstate(b.bases[i],0) for i in 1:3]...) tmp = copy(psi) testvec2 = ones(nsteps) .* sqrt(2)/get_nsteps(psi.basis) @@ -212,7 +212,7 @@ function test_interaction(b,Nw,idx1,idx2,order,wd2w1) end ξfun(t1,t2) = 1 - psi = tensor([i == bw_idx ? twophoton(b.bases[i],idx1,ξfun) : fockstate(b.bases[i],0) for i in 1:3]...) + psi = tensor([i == bw_idx ? twophoton(b.bases[i],idx1,ξfun,norm=true) : fockstate(b.bases[i],0) for i in 1:3]...) tmp = copy(psi) mul!(tmp,wd2w1,psi,1,0) diff --git a/test/runtests.jl b/test/runtests.jl index 07f2b99..e8df8d5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -26,8 +26,6 @@ println("Starting tests with $(Threads.nthreads()) threads out of `Sys.CPU_THREA @doset "operators" @doset "singlecavity" -@doset "detection_single" -@doset "detection_double" # TODO doctests need fixing #VERSION == v"1.8" && @doset "doctests" get(ENV,"QUANTUMOPTICS_JET_TEST","")=="true" && @doset "jet" diff --git a/test/test_operators.jl b/test/test_operators.jl index db7feaa..b5af45e 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -18,7 +18,7 @@ include("helper_functions.jl") adw = ad ⊗ w ξfun(t1,t2) = 1 - input = twophoton(bw,ξfun) + input = twophoton(bw,ξfun,norm=true) psi = fockstate(bc,0) ⊗ input psi_out = copy(psi) @@ -42,7 +42,7 @@ end adw = ad ⊗ w ξfun(t1) = 1 - input = onephoton(bw,ξfun) + input = onephoton(bw,ξfun,norm=true) psi = fockstate(bc,1) ⊗ input psi_out = copy(psi) mul!(psi_out,wda,psi,1,0.0) @@ -67,7 +67,7 @@ end adw = ad ⊗ w ξfun(t1) = 1 - input = onephoton(bw,ξfun) + input = onephoton(bw,ξfun,norm=true) psi = fockstate(bc,0) ⊗ input wda_out_2 = copy(psi) adw_out_2 = copy(psi) @@ -84,7 +84,7 @@ end wda = a ⊗ wd adw = ad ⊗ w - input = onephoton(bw,ξfun) + input = onephoton(bw,ξfun,norm=true) psi = fockstate(bc,0) ⊗ input adw_out_1 = copy(psi) wda_out_1 = copy(psi)