Reference

SuperBosonMPS

A finite-size matrix-product state type that represents mixed states in the superboson formalism [4].

References

• [4] Schmutz, M. Z Physik B 30, 97–106 (1978)

nmodes(v::SuperBosonMPS)

Return the number of "true" bosonic modes contained in the SuperBosonMPS v, i.e. half its actual length.

enlargelocaldim(v::AbstractMPS, newdim)

Return a new MPS with the local dimension increased to newdim.

Note that the new MPS will be defined on a new set of site indices, so it will be incompatible with the original one.

sb_siteinds(; nmodes, maxnumber)

Return a list of ITensor site indices suited to define a superbosonic state, with alternating "physical" and "ancillary" sites.

Keyword arguments (mandatory)

• nmodes: the number of "real-world" bosonic modes of the system (the actual MPS will have 2nmodes sites).
• maxnumber: maximum number of bosons allowed on each site. The local Hilbert spaces will be truncated to a dimension of maxnumber+1.

Example

julia> sb_siteinds(; nmodes=2, maxnumber=4)
4-element Vector{Index{Int64}}:
 (dim=5|id=962|"Boson,Site,n=1,phy")
 (dim=5|id=41|"Boson,Site,anc,n=1")
 (dim=5|id=794|"Boson,Site,n=2,phy")
 (dim=5|id=198|"Boson,Site,anc,n=2")

sb_outer(ψ::MPS)

Compute the projection operator $|ψ⟩⟨ψ| / ‖ψ‖²$, from the MPS ψ representing a pure state, returning a SuperBosonMPS object (of double the size).

expect(v::SuperBosonMPS, op::AbstractString...; kwargs...)
expect(v::SuperBosonMPS, op::Matrix{<:Number}...; kwargs...)
expect(v::SuperBosonMPS, ops; kwargs...)

Given a superbosonic MPS v and a single operator name or matrix, returns a vector of the expected value of the operator on each site of the SuperBosonMPS.

If multiple operator names are provided, returns a tuple of expectation value vectors.

If a container of operator names is provided, returns the same type of container with names replaced by vectors of expectation values.

Optional keyword arguments

• sites = 1:nmodes(v): compute expected values only for modes in the given range

Examples

N = 10

s = sb_siteinds("Boson", N)
v = sb_outer(random_mps(s; linkdims=4))
expect(v, "N")  # compute for all sites
expect(v, "N"; sites=2:4)  # compute for sites 2, 3 and 4
expect(v, "N"; sites=3)  # compute for site 3 only (output will be a scalar)
expect(v, ["A*Adag", "N"])  # compute A*Adag and N for all sites
expect(v, [0 0; 0 1])  # same as expect(v, "N") if maxnumber == 1

correlation_matrix(
    v::SuperBosonMPS, A::AbstractString, B::AbstractString; kwargs...
)

correlation_matrix(
    v::SuperBosonMPS, A::Matrix{<:Number}, B::Matrix{<:Number}; kwargs...
)

Given a SuperBosonMPS v representing a state $ρ$ and two strings or matrices A and B denoting operators (as recognized by the op function), computes the two-point correlation function matrix $C_{ij} = \textrm{tr}(A_i B_j ρ)$ using efficient MPS techniques. Returns the matrix C.

Optional keyword arguments

• sites = 1:nmodes(v): compute correlations only for sites in the given range
• ishermitian = false : if false, force independent calculations of the matrix elements above and below the diagonal, while if true assume they are complex conjugates.

Examples

julia> s = siteinds("Boson", 3; dim=4);

julia> vac = sb_outer(MPS(ComplexF64, s, "0"));

julia> correlation_matrix(vac, "x", "x")
3×3 Matrix{ComplexF64}:
 0.5+0.0im  0.0+0.0im  0.0+0.0im
 0.0-0.0im  0.5+0.0im  0.0+0.0im
 0.0-0.0im  0.0-0.0im  0.5+0.0im

julia> correlation_matrix(vac, "p", "x")
3×3 Matrix{ComplexF64}:
 0.0-0.5im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0-0.5im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0-0.5im

sample(m::SuperBosonMPS)

Given a "superbosonic" MPS m with unit trace, compute a Vector{Int} of nmodes(m) elements corresponding to one sample of the probability distribution defined by the components of the density matrix that the MPS represents.

firstmoments(v::Union{MPS,SuperBosonMPS}; warn_atol=1e-14)

Compute the first moments of the state represented by v (either a pure state, if v is an MPS, or a mixed one, if v is a SuperBosonMPS), displayed in the xpxp order, i.e. the vector $(⟨x_1⟩, ⟨p_1⟩, ⟨x_2⟩, ⟨p_2⟩, ..., ⟨x_N⟩, ⟨p_N⟩)$.

The warn_atol keyword argument can be used to adjust the threshold used by the function to warn when the computed moments are not real. The imaginary part of the result is always truncated.

covariancematrix(v::Union{MPS,SuperBosonMPS}; warn_atol=1e-14)

Compute the covariance matrix of the state represented by v (either a pure state, if v is an MPS, or a mixed one, if v is a SuperBosonMPS), displayed in the xpxp order.

The warn_atol keyword argument can be used to adjust the threshold used by the function to warn when the computed moments are not real. The imaginary part of the result is always truncated.

attenuate(v::SuperBosonMPS, attenuation, n)

Apply the attenuator channel $ρ ↦ \sum_{k=0}^{+∞} B_k ρ \adj{B_k}$ to the n-th mode of the state $ρ$ represented by the SuperBosonMPS v, where

\[B_k = \sum_{m=0}^{+∞} \binom{m+k}{k}^{\frac12} (1-η^2)^{\frac{k}{2}} η^m |m⟩⟨m+k|\]

and $η$ is the attenuation coefficient, such that $\sum_{k=0}^{+\infty} B_k ρ \adj{B_k}$ is equal to $|0⟩⟨0|$ when $η = 0$ and to $ρ$ when $η = 1$.

displace(v::Union{MPS,SuperBosonMPS}, α, n; kwargs...)

Apply the displacement operator

\[D(α) = \exp(α \adj{a} - \conj{α} a)\]

with $α ∈ ℂ$, to the n-th mode of the state represented by v.

displace(v::Union{MPS,SuperBosonMPS}, α; kwargs...)

Apply the displacement operator

\[⨂_{i=1}^{n} D(α_i), \quad D(α) = \exp(α \adj{a} - \conj{α} a)\]

with $α_i ∈ ℂ$, to all modes of the state represented by v.

squeeze(v::Union{MPS,SuperBosonMPS}, z, n; kwargs...)

Apply the squeezing operator

\[S(z) = \exp\bigl(\tfrac12 z (\adj{a})^2 - \tfrac12 \conj{z} a^2\bigr)\]

with $z ∈ ℂ$, to the n-th mode of the state represented by v.

squeeze(v::Union{MPS,SuperBosonMPS}, z; kwargs...)

Apply the squeezing operator

\[⨂_{i=1}^{n} S(z_i), \quad S(z) = \exp\bigl(\tfrac12 z (\adj{a})^2 - \tfrac12 \conj{z} a^2\bigr),\]

with $z_i ∈ ℂ$, to all modes of the state represented by v.

squeeze2(v::Union{MPS,SuperBosonMPS}, z, n1, n2; kwargs...)

Apply the two-mode squeezing operator

\[S₂(z) = \exp(z \adj{a} ⊗ \adj{a} - \conj{z} a ⊗ a)\]

with $z ∈ ℂ$, on modes n1 and n2 of the state represented by v.

beamsplitter(v::Union{MPS,SuperBosonMPS}, z, n1, n2; kwargs...)

Apply the beam-splitter operator

\[B(z) = \exp(z \adj{a} ⊗ a - \conj{z} a ⊗ \adj{a})\]

with $z ∈ ℂ$, to modes n1 and n2 of the state represented by v.

normal_mode_decomposition(g::GaussianState, N, maxnumber; kwargs...)

Compute the normal-mode decomposition [6, Eq. (3.60) at page 47] of the Gaussian state g up to the N largest eigenvalues, calculated by considering only the k-particle sectors for k ≤ maxnumber on each mode. Return vals, nums, S where vals contains the eigenvalues in decreasing order, nums contains the occupation numbers of the Fock basis vector associated to the corresponding value in vals, and S is the symplectic matrix from the Williamson decomposition of g.

The output satisfies the following identities:

julia> D, S = williamson(Symmetric(g.covariance_matrix));

julia> d = diag(D)[1:2:end];

julia> normalmode_evs, nums, S' = normal_mode_decomposition(g, N, maxnumber);

julia> all(
    val == prod(normal_mode_eigenvalue(d[j], m[j]) for j in eachindex(m)) for
    (val, m) in zip(normalmode_evs, nums)
)
true

julia> S == S'
true

References

• [6] Serafini, A. (2023)

franckcondon(m, α, Wₗ, S, Wᵣ, n)
franckcondon(m, Wₗ, S, Wᵣ, n)

Compute the matrix element $⟨m| D(α) U(W_l) U(S) U(W_r) |n⟩$ according to the algorithm presented in Ref. [2]. Wₗ, S and Wᵣ are symplectic matrices, in particular Wₗ and Wᵣ are also orthogonal, and α is a complex vector (which defaults to the zero vector). If $N$ is the number of modes of the system, then

• m and n are tuples of $N$ natural numbers,
• α is a vector of $N$ complex numbers,
• Wₗ, S and Wᵣ are $2N × 2N$ matrices.

References

• [2] Quesada, N. J Chem Phys 150.16 (2019)

MPS(g::GaussianState; maxdim, maxnumber, kwargs...)

Build an MPS representation of the Gaussian state g with bond dimension up to maxdim, truncating the Fock space of each mode at the maxnumber-particle sector.

optimise(g::GaussianState; verbose=false, scs_eps=nothing)

Return gₚ, W where gₚ is a new Gaussian state and W is a positive semi-definite matrix such that W + gₚ.covariance_matrix == g.covariance_matrix and gₚ contains a smaller number of photons.

Set verbose = true to print the information provided by SCS about the optimisation.

sample_displaced(ψ::MPS, W; nsamples, nsamples_per_displacement, eval_atol=0, kwargs...)

Apply random displacements sampled from the positive semi-definite matrix W to the pure state ψ and sample nsamples elements from the resulting state.

Return s, αs such that:

• s is a matrix of UInt8 elements such that each column is a sample drawn from the final state, with a new displacement vector for each nsamples_per_displacement draws. (The matrix will actually have a number of columns equal to nsamples_per_displacement * floor(nsamples / nsamples_per_displacement).)
• αs is a matrix whose columns are the displacement vectors drawn from W, such that αs[:, k] is used for the samples from (k-1)N to kN, with N = floor(nsamples / nsamples_per_displacement).

The eval_atol keyword argument is used as threshold to decide whether an eigenvalue of W must be considered zero (usually it should be of the same order of the eps tolerances of the SCS optimiser).