update docs and format source

docs/src/guide.md

@@ -81,13 +81,16 @@ N=301; # discretization size
 model=OneSpinModel(T,L,N,B)
 println(model)
 ```
+This provides us an overview of the model. It's a single spin shuttling problem with initial state `Ψ₀` and an Ornstein-Uhlenbeck noise. The total time of simulation is `T`, which is discretized into `N` steps.
 
-
-The fidelity of the spin state after shuttling can be calculated using numerical integration of the covariance matrix.  
-
-This provides us an overview of the model. It's a single spin shuttling problem with initial state `Ψ₀` and an Ornstein-Uhlenbeck noise. The total time of simulation is `T`, which is discretized into `N` steps.  
-
-The state fidelity after such a quantum process can be obtained by different numerical methods. 
+The effective noise of this spin qubit is completely characterized by its covariance matrix.  
+```@example quickstart
+heatmap(collect(sqrt.(model.R.Σ)), title="sqrt cov, 1-spin one-way shuttling", 
+size=(400,300), 
+xlabel="t1", ylabel="t2", dpi=300,
+right_margin=5Plots.mm)
+```
+The state fidelity after such a quantum process can be obtained using numerical integration of the covariance matrix.  
 ```@example quickstart
 f1=averagefidelity(model); # direct integration
 
@@ -119,7 +122,6 @@ We can check that the fidelity between the initial and final state is consistent
 f=(Ψ'*ρt*Ψ)
 ```
 
-
 ## Dephasing of entangled spin pairs during shuttling. 
 Following the approach above, we can further explore the multi-spin system. 
 The general abstraction on such a problem is given by the data type `ShuttlingModel`.  
@@ -148,7 +150,7 @@ plot!(model.R.P[N+1:2N,1], label="x2(t)")
 
 
 ```@example quickstart
-heatmap(collect(model.R.Σ)*1e3, title="covariance matrix, two spin EPR", 
+heatmap(collect(model.R.Σ)*1e3, title="covariance, 2-spin sequential shuttling", 
 size=(400,300), 
 xlabel="t1", ylabel="t2", dpi=300,
 right_margin=5Plots.mm)

docs/src/manual.md

@@ -31,7 +31,6 @@ TwoSpinParallelModel
 dephasingmatrix
 ```
 
-
 ```@docs
 fidelity
 ```

src/SpinShuttling.jl

@@ -11,9 +11,9 @@ include("integration.jl")
 include("analytics.jl")
 include("stochastics.jl")
 
-export ShuttlingModel, OneSpinModel, TwoSpinModel, 
-OneSpinForthBackModel, TwoSpinParallelModel, RandomFunction, CompositeRandomFunction,
-OrnsteinUhlenbeckField, PinkBrownianField
+export ShuttlingModel, OneSpinModel, TwoSpinModel,
+    OneSpinForthBackModel, TwoSpinParallelModel, RandomFunction, CompositeRandomFunction,
+    OrnsteinUhlenbeckField, PinkBrownianField
 export averagefidelity, fidelity, sampling, characteristicfunction, characteristicvalue
 export dephasingmatrix, covariance, covariancematrix
 export W
@@ -49,9 +49,9 @@ function Base.show(io::IO, model::ShuttlingModel)
     println(io, "Time Discretization: N=$(model.N)")
     println(io, "Process Time: T=$(model.T)")
     println(io, "Shuttling Paths:")
-    t=range(0, model.T, model.N)
-    fig=lineplot(t, model.X[1].(t); width=30, height=9,
-    name="x1(t)")
+    t = range(0, model.T, model.N)
+    fig = lineplot(t, model.X[1].(t); width=30, height=9,
+        name="x1(t)")
     for i in 2:model.n
         lineplot!(fig, t, model.X[i].(t), name="x$i(t)")
     end
@@ -101,12 +101,12 @@ with total time `T` in `μs` and length `L` in `μm`.
 - `B::GaussianRandomField`: Noise field
 - `v::Real`: Velocity of the shuttling
 """
-function OneSpinForthBackModel(t::Real, T::Real, L::Real, N::Int, B::GaussianRandomField) 
-    x(t::Real, v::Real, L::Real)::Real = (t=t%(2L/v); v*t < L ? v*t : 2L-v*t)
-    return OneSpinModel(1 / √2 * [1, 1], t, N, B, τ -> x(τ, 2L/T, L))
+function OneSpinForthBackModel(t::Real, T::Real, L::Real, N::Int, B::GaussianRandomField)
+    x(t::Real, v::Real, L::Real)::Real = (t = t % (2L / v); v * t < L ? v * t : 2L - v * t)
+    return OneSpinModel(1 / √2 * [1, 1], t, N, B, τ -> x(τ, 2L / T, L))
 end
 
-function OneSpinForthBackModel(T::Real, L::Real, N::Int, B::GaussianRandomField) 
+function OneSpinForthBackModel(T::Real, L::Real, N::Int, B::GaussianRandomField)
     return OneSpinForthBackModel(T, T, L, N, B)
 end
 
@@ -174,7 +174,7 @@ The total shuttling time is `T` and the length of the path is `L` in `μm`.
 """
 function TwoSpinParallelModel(T::Real, D::Real, L::Real, N::Int,
     B::GaussianRandomField)
-    @assert length(B.θ)>=3 
+    @assert length(B.θ) >= 3
     x₁(t::Real)::Tuple{Real,Real} = (L / T * t, 0)
     x₂(t::Real)::Tuple{Real,Real} = (L / T * t, D)
     Ψ = 1 / √2 .* [0, 1, -1, 0]
@@ -233,8 +233,8 @@ function sampling(samplingfunction::Function, M::Int)::Union{Tuple{Real,Real},Tu
     Q = copy(A)
     for k in 1:M
         x = samplingfunction(k)::Union{Real,Vector{<:Real}}
-        Q = Q +(k-1)/k*(x-A).^ 2
-        A = A + (x-A)/k
+        Q = Q + (k - 1) / k * (x - A) .^ 2
+        A = A + (x - A) / k
     end
     return A, Q / (M - 1)
 end
@@ -248,15 +248,15 @@ function parallelsampling(samplingfunction::Function, M::Int)::Union{Tuple{Real,
         end
         A = mean(cache, dims=2)
         Q = var(cache, dims=2)
-        return A, Q 
+        return A, Q
     else
         cache = zeros(M)
         @threads for i in 1:M
             cache[i] = samplingfunction(i)
         end
         A = mean(cache)
         Q = var(cache)
-        return A, Q 
+        return A, Q
     end
 end
 
@@ -282,34 +282,34 @@ end
 - `p::Int`: index of spin state, range from (1,2^n)
 - `n::Int`: number of spins
 """
-function m(i::Int,p::Int,n::Int)
-    1/2-digits(p, base=2, pad=n)[i]
+function m(i::Int, p::Int, n::Int)
+    1 / 2 - digits(p, base=2, pad=n)[i]
 end
 
 """
 Calculate the dephasing matrix of a given spin shuttling model.
 """
 function dephasingmatrix(model::ShuttlingModel)::Symmetric{<:Real}
-    n=model.n
-    W=zeros(2^n,2^n)
+    n = model.n
+    W = zeros(2^n, 2^n)
     for j in 1:2^n
-        W[j,j]=1
+        W[j, j] = 1
         for k in 1:j-1
-            c=[trunc(Int,m(i,j-1,n)-m(i,k-1,n)) for i in 1:n]
+            c = [trunc(Int, m(i, j - 1, n) - m(i, k - 1, n)) for i in 1:n]
             R = CompositeRandomFunction(model.R, c)
-            W[j,k] = characteristicvalue(R)
-            W[k,j] = W[j,k]
+            W[j, k] = characteristicvalue(R)
+            W[k, j] = W[j, k]
         end
     end
     return Symmetric(W)
 end
 
 function dephasingcoeffs(n::Int)::Array{Real,3}
-    M=zeros(2^n,2^n, n)
+    M = zeros(2^n, 2^n, n)
     for j in 1:2^n
         for k in 1:2^n
-            c=[m(i,j-1,n)-m(i,k-1,n) for i in 1:n]
-            M[j,k, :] = c
+            c = [m(i, j - 1, n) - m(i, k - 1, n) for i in 1:n]
+            M[j, k, :] = c
         end
     end
     return M
@@ -330,11 +330,11 @@ function fidelity(model::ShuttlingModel, randseq::Vector{<:Real}; vector::Bool=f
         Z = A
     elseif model.n == 2
         # only valid for two-spin EPR pair, ψ=1/√2(|↑↓⟩-|↓↑⟩)
-        Z = A[1:N] - A[N+1:end] 
+        Z = A[1:N] - A[N+1:end]
     else
         Z = missing
     end
-    ϕ = vector ? cumsum(Z)* dt : sum(Z) * dt
+    ϕ = vector ? cumsum(Z) * dt : sum(Z) * dt
     return (1 .+ cos.(ϕ)) / 2
 end
 
@@ -349,41 +349,41 @@ Analytical dephasing factor of a one-spin shuttling model.
 - `B<:GaussianRandomField`: Noise field, Ornstein-Uhlenbeck or Pink-Brownian
 - `path::Symbol`: Path of the shuttling model, `:straight` or `:forthback`
 """
-function W(T::Real,L::Real,B::OrnsteinUhlenbeckField; path=:straight)::Real
-    κₜ=B.θ[1]
-    κₓ=B.θ[2]
-    σ =B.σ
-    β = κₜ*T
-    γ = κₓ*L
+function W(T::Real, L::Real, B::OrnsteinUhlenbeckField; path=:straight)::Real
+    κₜ = B.θ[1]
+    κₓ = B.θ[2]
+    σ = B.σ
+    β = κₜ * T
+    γ = κₓ * L
     if path == :straight
-        return exp(- σ^2/(4*κₜ*κₓ)/κₜ^2*P1(β, γ)/2)
+        return exp(-σ^2 / (4 * κₜ * κₓ) / κₜ^2 * P1(β, γ) / 2)
     elseif path == :forthback
-        β/=2
-        return exp(- σ^2/(4*κₜ*κₓ)/κₜ^2*(P1(β, γ)+P4(β,γ)))
+        β /= 2
+        return exp(-σ^2 / (4 * κₜ * κₓ) / κₜ^2 * (P1(β, γ) + P4(β, γ)))
     else
         error("Path not recognized. Use :straight or :forthback for one-spin shuttling model.")
     end
 end
 
 function W(T::Real, L::Real, B::PinkBrownianField)::Real
-    β= T.*B.γ
-    γ= L*B.θ[1]
-    return exp(-B.σ^2*T^2*F3(β,γ))
+    β = T .* B.γ
+    γ = L * B.θ[1]
+    return exp(-B.σ^2 * T^2 * F3(β, γ))
 end
 
 
 """
 Analytical dephasing factor of a sequenced two-spin EPR pair shuttling model.
 """
-function W(T0::Real,T1::Real,L::Real,B::OrnsteinUhlenbeckField; path=:sequenced)::Real
-    κₜ=B.θ[1]
-    κₓ=B.θ[2]
-    σ =B.σ
-    τ = κₜ*T0
-    β = κₜ*T1
-    γ = κₓ*L
+function W(T0::Real, T1::Real, L::Real, B::OrnsteinUhlenbeckField; path=:sequenced)::Real
+    κₜ = B.θ[1]
+    κₓ = B.θ[2]
+    σ = B.σ
+    τ = κₜ * T0
+    β = κₜ * T1
+    γ = κₓ * L
     if path == :sequenced
-        return exp(-σ^2/(4*κₜ*κₓ)/κₜ^2*(F1(β, γ, τ)-F2(β, γ, τ)))
+        return exp(-σ^2 / (4 * κₜ * κₓ) / κₜ^2 * (F1(β, γ, τ) - F2(β, γ, τ)))
     elseif path == :parallel
         missing("Parallel path not implemented yet.")
     else

src/analytics.jl

@@ -1,73 +1,72 @@
-function C1(β::Real,γ::Real,τ::Real)::Real
-    if β<=τ
-        return ℯ^(-β-γ-τ)*(1+ℯ^(2β)-2ℯ^τ+2ℯ^(β+τ)*(-β+τ))
+function C1(β::Real, γ::Real, τ::Real)::Real
+    if β <= τ
+        return ℯ^(-β - γ - τ) * (1 + ℯ^(2β) - 2ℯ^τ + 2ℯ^(β + τ) * (-β + τ))
     else
-        return ℯ^(-β-γ-τ)*(-1+ℯ^τ)^2
+        return ℯ^(-β - γ - τ) * (-1 + ℯ^τ)^2
     end
 end
 
-function C2(β::Real,γ::Real,τ::Real)::Real
-    if β<=τ
-        return ℯ^(-γ)*β*((ℯ^(-τ)*(-ℯ^β+ℯ^γ))/(β-γ)+(ℯ^(-β)*(γ-2ℯ^β*(β+γ)+ℯ^(β+γ)*(2β+γ)))/(γ*(β+γ)))
+function C2(β::Real, γ::Real, τ::Real)::Real
+    if β <= τ
+        return ℯ^(-γ) * β * ((ℯ^(-τ) * (-ℯ^β + ℯ^γ)) / (β - γ) + (ℯ^(-β) * (γ - 2ℯ^β * (β + γ) + ℯ^(β + γ) * (2β + γ))) / (γ * (β + γ)))
     else
-        return (ℯ^(-(((β+γ)*(β+3τ))/β))*β*(2ℯ^(β+γ+3τ+(2*γ*τ)/β)*(-1+ℯ^((γ*τ)/β))*β^2-ℯ^((2+(3*γ)/β)*τ)*(-1+ℯ^τ)*γ*(ℯ^τ*(β-γ)+ℯ^(β+γ)*(β+γ))))/((β-γ)*γ*(β+γ))
+        return (ℯ^(-(((β + γ) * (β + 3τ)) / β)) * β * (2ℯ^(β + γ + 3τ + (2 * γ * τ) / β) * (-1 + ℯ^((γ * τ) / β)) * β^2 - ℯ^((2 + (3 * γ) / β) * τ) * (-1 + ℯ^τ) * γ * (ℯ^τ * (β - γ) + ℯ^(β + γ) * (β + γ)))) / ((β - γ) * γ * (β + γ))
     end
 end
 
-function C3(β::Real,γ::Real,τ::Real)::Real
-    if β<=τ
-        return (ℯ^(-β-γ-τ)*β^2*(2(-1+ℯ^(2β))*β*γ+(1+ℯ^(2β)+2ℯ^(β+γ)*(-1+γ))*γ^2+β^2*(1+ℯ^(2β)-2ℯ^(β+γ)*(1+γ))))/(β^2-γ^2)^2
+function C3(β::Real, γ::Real, τ::Real)::Real
+    if β <= τ
+        return (ℯ^(-β - γ - τ) * β^2 * (2(-1 + ℯ^(2β)) * β * γ + (1 + ℯ^(2β) + 2ℯ^(β + γ) * (-1 + γ)) * γ^2 + β^2 * (1 + ℯ^(2β) - 2ℯ^(β + γ) * (1 + γ)))) / (β^2 - γ^2)^2
     else
-        return (1/((β^2-γ^2)^2))ℯ^(-(((β+γ)*(β+τ))/β))*β^2*(ℯ^((γ*τ)/β)*(β-γ)^2+ℯ^((2+γ/β)*τ)*(β-γ)^2-2ℯ^(β+γ+(γ*τ)/β)*(-((-1+γ)*γ^2)+β^2*(1+γ))+2ℯ^(β+γ+τ)*(β^3-β*(-2+γ)*γ-β^2*τ+γ^2*τ))
+        return (1 / ((β^2 - γ^2)^2))ℯ^(-(((β + γ) * (β + τ)) / β)) * β^2 * (ℯ^((γ * τ) / β) * (β - γ)^2 + ℯ^((2 + γ / β) * τ) * (β - γ)^2 - 2ℯ^(β + γ + (γ * τ) / β) * (-((-1 + γ) * γ^2) + β^2 * (1 + γ)) + 2ℯ^(β + γ + τ) * (β^3 - β * (-2 + γ) * γ - β^2 * τ + γ^2 * τ))
     end
 end
 
-function C4(β::Real,γ::Real,τ::Real)::Real
-    if β<=τ
-        return ℯ^(-γ)*β*((ℯ^-τ*(-ℯ^β+ℯ^γ))/(β-γ)+(ℯ^(-β)*(γ-2ℯ^β*(β+γ)+ℯ^(β+γ)*(2β+γ)))/(γ*(β+γ)))
+function C4(β::Real, γ::Real, τ::Real)::Real
+    if β <= τ
+        return ℯ^(-γ) * β * ((ℯ^-τ * (-ℯ^β + ℯ^γ)) / (β - γ) + (ℯ^(-β) * (γ - 2ℯ^β * (β + γ) + ℯ^(β + γ) * (2β + γ))) / (γ * (β + γ)))
     else
-        return (ℯ^(-(((β+γ)*(β+τ))/β))*β*(2ℯ^(β+γ+τ)*(-1+ℯ^((γ*τ)/β))*β^2-ℯ^((γ*τ)/β)*(-1+ℯ^τ)*(ℯ^(β+γ)+ℯ^τ)*β*γ+ℯ^((γ*τ)/β)*(-1+ℯ^τ)*(-ℯ^(β+γ)+ℯ^τ)*γ^2))/(β^2*γ-γ^3)
+        return (ℯ^(-(((β + γ) * (β + τ)) / β)) * β * (2ℯ^(β + γ + τ) * (-1 + ℯ^((γ * τ) / β)) * β^2 - ℯ^((γ * τ) / β) * (-1 + ℯ^τ) * (ℯ^(β + γ) + ℯ^τ) * β * γ + ℯ^((γ * τ) / β) * (-1 + ℯ^τ) * (-ℯ^(β + γ) + ℯ^τ) * γ^2)) / (β^2 * γ - γ^3)
     end
 end
 
-function P1(β::Real,γ::Real)::Real
-    return -((2β^2 *(1-ℯ^(-β-γ)-β-γ))/(β+γ)^2)
+function P1(β::Real, γ::Real)::Real
+    return -((2β^2 * (1 - ℯ^(-β - γ) - β - γ)) / (β + γ)^2)
 end
 
-function P2(β::Real,γ::Real,τ::Real)::Real
-    return 2*(-1+ℯ^-τ+τ)
+function P2(β::Real, γ::Real, τ::Real)::Real
+    return 2 * (-1 + ℯ^-τ + τ)
 end
 
-function P3(β::Real,γ::Real,τ::Real)::Real
-    return (ℯ^(-β-γ-τ)*(-1+ℯ^(β+γ))*(-1+ℯ^τ)*β)/(β+γ)
+function P3(β::Real, γ::Real, τ::Real)::Real
+    return (ℯ^(-β - γ - τ) * (-1 + ℯ^(β + γ)) * (-1 + ℯ^τ) * β) / (β + γ)
 end
 
-function P4(β::Real,γ::Real)::Real
-    return ((ℯ^(-2β)-1)*(γ/β)-2*ℯ^(-β-γ)+ℯ^(-2β)+1)/(1-γ^2/β^2)
+function P4(β::Real, γ::Real)::Real
+    return ((ℯ^(-2β) - 1) * (γ / β) - 2 * ℯ^(-β - γ) + ℯ^(-2β) + 1) / (1 - γ^2 / β^2)
 end
 
 """
 Ancillary function for the dephasing of the sequential shuttling pf 
 Bell state under OU sheets of noise
 """
-function F1(β::Real,γ::Real,τ::Real)::Real
-    return P1(β,γ)+P2(β,γ,τ)+2*P3(β,γ,τ)
+function F1(β::Real, γ::Real, τ::Real)::Real
+    return P1(β, γ) + P2(β, γ, τ) + 2 * P3(β, γ, τ)
 end
 
-function F2(β::Real,γ::Real,τ::Real)::Real
-    return C1(β,γ,τ)+C2(β,γ,τ)+C3(β,γ,τ)+C4(β,γ,τ)
+function F2(β::Real, γ::Real, τ::Real)::Real
+    return C1(β, γ, τ) + C2(β, γ, τ) + C3(β, γ, τ) + C4(β, γ, τ)
 end
 
 """
 Ancillary function for the dephasing of the Pink-Brownian noise.
 """
-function F3(β::Tuple{Real,Real},γ::Real)::Real
-    F(β::Real)=1/2*(expinti(-β)+(1-exp(-β))/β^2+(exp(-β)-2)/β)
-    F(β::Real,γ::Real)::Real=1/γ^2*(exp(-γ)*expinti(-β)+(γ-1)*(expinti(-β-γ)-log((β+γ)/β))-γ*((1-exp(-β-γ))/(β+γ)))
-    if γ==0 # pure 1/f noise
-        return (F(β[2])-F(β[1]))/log(β[2]/β[1])
+function F3(β::Tuple{Real,Real}, γ::Real)::Real
+    F(β::Real) = 1 / 2 * (expinti(-β) + (1 - exp(-β)) / β^2 + (exp(-β) - 2) / β)
+    F(β::Real, γ::Real)::Real = 1 / γ^2 * (exp(-γ) * expinti(-β) + (γ - 1) * (expinti(-β - γ) - log((β + γ) / β)) - γ * ((1 - exp(-β - γ)) / (β + γ)))
+    if γ == 0 # pure 1/f noise
+        return (F(β[2]) - F(β[1])) / log(β[2] / β[1])
     else
-        return (F(β[2],γ)-F(β[1],γ))/log(β[2]/β[1])
+        return (F(β[2], γ) - F(β[1], γ)) / log(β[2] / β[1])
     end
 end
-